Add log space fallback for PPF func (#2)

* Add log space fallback for PPF func

* PR feedback

* small fixes

* add global fallback flag to avoid scipy instantiation problems

* poetry lock

* pull function definition out of loop

* randomize test order to catch any global fallback flag issues

* Bump version=0.1.7

* don't unbox to python type for a single scalar argument

* formatting

Author: hbmartin

.github/workflows/pytest-poetry.yml

@@ -25,9 +25,9 @@ jobs:
           cache-dependency-path: poetry.lock
       - name: Install dependencies
         run: poetry install --no-interaction
-      - name: Run tests
-        run: poetry run pytest 
       - name: run black
         run: poetry run black . --check
       - name: Run ruff
-        run: poetry run ruff check .
+        run: poetry run ruff check .
+      - name: Run tests
+        run: poetry run pytest --random-order

betapert/__init__.py

@@ -9,6 +9,8 @@
 
 from betapert import funcs
 
+FALLBACK = None
+
 
 class PERT(scipy.stats.rv_continuous):
     """The `PERT distribution <https://en.wikipedia.org/wiki/PERT_distribution>`_ is defined by the
@@ -60,7 +62,7 @@ def _stats(self, mini, mode, maxi):
         return funcs.stats(mini, mode, maxi)
 
     def _ppf(self, q, mini, mode, maxi):
-        return funcs.ppf(q, mini, mode, maxi)
+        return funcs.ppf(q, mini, mode, maxi, fallback=FALLBACK)
 
     def _rvs(self, mini, mode, maxi, size=None, random_state=None):
         return funcs.rvs(mini, mode, maxi, size=size, random_state=random_state)
@@ -116,7 +118,7 @@ def _stats(self, mini, mode, maxi, lambd):
         return funcs.stats(mini, mode, maxi, lambd)
 
     def _ppf(self, q, mini, mode, maxi, lambd):
-        return funcs.ppf(q, mini, mode, maxi, lambd)
+        return funcs.ppf(q, mini, mode, maxi, lambd, fallback=FALLBACK)
 
     def _rvs(self, mini, mode, maxi, lambd, size=None, random_state=None):
         return funcs.rvs(mini, mode, maxi, lambd, size=size, random_state=random_state)

betapert/funcs.py

@@ -2,12 +2,63 @@
 
 This module contains the core mathematical functions used by the PERT and modified PERT distribution
 classes. Each function takes the distribution parameters (minimum, mode, maximum, and optionally
-lambda) and implementsa specific statistical operation like pdf, cdf, etc.
+lambda) and implements a specific statistical operation like pdf, cdf, etc.
 """
 
 import numpy as np
+import scipy.optimize
 import scipy.stats
 
+# Avoid log(0) or log(1) which would cause -inf or 0
+_CLIP_EPSILON = 1e-15
+_BRENTQ_BOUND = 1e-10
+
+
+def _ppf_fallback_log_space(q, mini, mode, maxi, lambd):
+    """Use log-space to avoid numerical issues with extreme probabilities"""
+    alpha, beta = _calc_alpha_beta(mini, mode, maxi, lambd)
+
+    # Handle scalar and array inputs consistently
+    _q = np.atleast_1d(q)
+    results = np.zeros_like(_q, dtype=float)
+
+    # Define the equation to solve: log(CDF(x)) - log(q) = 0
+    def make_log_cdf_eq(qi_val):
+        log_qi = np.log(np.clip(qi_val, _CLIP_EPSILON, 1 - _CLIP_EPSILON))
+
+        def log_cdf_eq(x_normalized):
+            # Ensure x_normalized stays in [0,1]
+            x_clamped = np.clip(x_normalized, _CLIP_EPSILON, 1 - _CLIP_EPSILON)
+            return scipy.stats.beta.logcdf(x_clamped, alpha, beta) - log_qi
+
+        return log_cdf_eq
+
+    for i, qi in enumerate(_q.flat):
+        try:
+            # Use brentq instead of fsolve, guaranteed convergence within bounds
+            x_normalized = scipy.optimize.brentq(
+                make_log_cdf_eq(qi),
+                _BRENTQ_BOUND,
+                1 - _BRENTQ_BOUND,
+            )
+            results.flat[i] = mini + (maxi - mini) * x_normalized
+
+        except (ValueError, RuntimeError):
+            # ValueError: Invalid function values, convergence issues, or invalid bounds
+            # RuntimeError: Maximum iterations exceeded, numerical problems
+            # Fallback to clamped ppf if log-space fails
+            qi_safe = np.clip(qi, _CLIP_EPSILON, 1 - _CLIP_EPSILON)
+            x_normalized = scipy.stats.beta.ppf(qi_safe, alpha, beta)
+            results[i] = mini + (maxi - mini) * x_normalized
+
+    # Returns scalar for scalar input, array for array input
+    return results[0] if np.isscalar(q) else results
+
+
+_ppf_fallbacks = {
+    "log": _ppf_fallback_log_space,
+}
+
 
 def _calc_alpha_beta(mini, mode, maxi, lambd):
     """Calculate alpha and beta parameters for the underlying beta distribution.
@@ -42,9 +93,13 @@ def sf(x, mini, mode, maxi, lambd=4):
     return scipy.stats.beta.sf((x - mini) / (maxi - mini), alpha, beta)
 
 
-def ppf(q, mini, mode, maxi, lambd=4):
+def ppf(q, mini, mode, maxi, lambd=4, *, fallback=None):
     alpha, beta = _calc_alpha_beta(mini, mode, maxi, lambd)
-    return mini + (maxi - mini) * scipy.stats.beta.ppf(q, alpha, beta)
+    _beta_ppf = mini + (maxi - mini) * scipy.stats.beta.ppf(q, alpha, beta)
+    # Use fallback if any values are NaN
+    if fallback is not None and np.any(np.atleast_1d(np.isnan(_beta_ppf))):
+        return _ppf_fallbacks[fallback](q, mini, mode, maxi, lambd)
+    return _beta_ppf
 
 
 def isf(q, mini, mode, maxi, lambd=4):

poetry.lock

@@ -1,7 +1,7 @@
 [tool]
 [tool.poetry]
 name = "beta-pert-dist-scipy"
-version = "0.1.6"
+version = "0.1.7"
 homepage = "https://github.com/hbmartin/betapert"
 description = "Top-level package for beta-PERT distribution."
 authors = ["Tom Adamczewski <[email protected]>", "Harold Martin [email protected]"]
@@ -19,14 +19,15 @@ packages = [
 
 [tool.poetry.dependencies]
 python = ">=3.11"
-scipy = "^1.14.0"
+scipy = ">=1.14.1"
 
 [tool.poetry.group.dev.dependencies]
 coverage = "*"
 pytest = ">=7.2.0"
 black = {extras = ["d"], version = "*"}
 matplotlib = "^3.9.0"
 ruff = "^0.12.3"
+pytest-random-order = "^1.2.0"
 
 
 
@@ -54,7 +55,7 @@ select = ["ALL"]
 ignore = ["ANN001", "ANN201", "ANN202", "ARG001", "D203", "D205", "D213", "D400", "D415", "PLR0913"]
 
 [tool.ruff.lint.per-file-ignores]
-"tests/*.py" = ["ARG002", "D100", "D101", "D102", "D103", "D104", "D200", "D212", "D401", "D404", "E501", "E731", "NPY002", "PLR2004", "RET503", "S101"]
+"tests/*.py" = ["ANN002", "ANN003", "ARG002", "ANN206", "D100", "D101", "D102", "D103", "D104", "D200", "D212", "D401", "D404", "E501", "E731", "EM101", "NPY002", "PLR2004", "PT011", "RET503", "S101", "SLF001", "TRY003"]
 "betapert/funcs.py" = ["D103", "RUF002", "RUF003"]
 
 [tool.ruff.format]

pyproject.toml

@@ -1,2 +1,3 @@
 [pytest]
-addopts = --doctest-modules
+addopts = --doctest-modules
+filterwarnings = error::pytest.PytestCollectionWarning