From f6555d902a422f913a397d83a6e102a96844bd12 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 4 Aug 2026 09:00:51 +0000 Subject: [PATCH 01/13] Say why a Turnbull fit does not equal a Kaplan-Meier fit MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Turnbull.fit defaults to turnbull_estimator="Fleming-Harrington" while KaplanMeier.fit is KM. The EM recovers the same r and d either way; the three options differ only in how those become a survival curve. So comparing a default Turnbull fit against KaplanMeier and reading the gap as a defect is an easy mistake, and the docstring gave no hint of it. It is the mistake #260 was filed on, and the mistake made again while checking whether #260 was still open — twice is enough to write it down. On x=[2,3,3,4,5,6], tl=[0,0,1,1,2,2] the survival at 2 is 0.750 under KM, 0.765 under FH, 0.779 under NA. Matched, Turnbull agrees with KaplanMeier to ~1e-9 on sf and cb across right-censored and left-truncated data. Only KM is the NPMLE. Maximising the truncated likelihood directly over the mass vector gives 0.750; FH's 0.765 scores worse on that same likelihood, which is what an exp(-H) construction should do. FH is the default for tail and zero-inflation behaviour (v0.8.0), not because it maximises anything. The new test pins the three figures and the NPMLE identity against a brute-force Nelder-Mead maximisation, so the docstring cannot quietly stop being true. No behaviour change. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts --- docs/changelog.rst | 27 ++++++++ .../univariate/nonparametric/test_turnbull.py | 66 +++++++++++++++++++ .../nonparametric/nonparametric_fitter.py | 20 ++++++ 3 files changed, 113 insertions(+) diff --git a/docs/changelog.rst b/docs/changelog.rst index 3b009543..692f77b5 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -1,6 +1,33 @@ Changelog ========= +v0.19.1 (unreleased) +-------------------- + +- **Documented why a Turnbull fit does not equal a Kaplan-Meier fit.** + ``Turnbull.fit`` defaults to ``turnbull_estimator="Fleming-Harrington"`` + while ``KaplanMeier.fit`` is, unsurprisingly, KM. The Turnbull EM + recovers the same ``r`` and ``d`` either way; the three estimator + options then differ in how those become a survival curve. Comparing the + default against ``KaplanMeier`` and reading the gap as a defect is an + easy mistake — it is the one #260 was filed on, and the one made again + while checking whether #260 was still open. + + On ``x=[2,3,3,4,5,6], tl=[0,0,1,1,2,2]`` the survival at 2 is 0.750 + under KM, 0.765 under FH and 0.779 under NA. With the estimator matched, + Turnbull agrees with ``KaplanMeier`` to around 1e-9 on both ``sf`` and + ``cb``, across right-censored and left-truncated data. + + Only the KM option is the non-parametric MLE. Maximising the truncated + likelihood directly over the mass vector gives 0.750; FH's 0.765 scores + worse on that same likelihood, as an ``exp(-H)`` construction should. + FH is the default because it behaves better in the far tails and on + zero-inflated data (v0.8.0), not because it maximises anything. The + docstring now says all of this, and a test pins the three figures and + the NPMLE identity against a brute-force maximisation. + + No behaviour change. + v0.19.0 (4 August 2026) ----------------------- diff --git a/surpyval/tests/univariate/nonparametric/test_turnbull.py b/surpyval/tests/univariate/nonparametric/test_turnbull.py index e72cdedf..261e42c7 100644 --- a/surpyval/tests/univariate/nonparametric/test_turnbull.py +++ b/surpyval/tests/univariate/nonparametric/test_turnbull.py @@ -367,3 +367,69 @@ def test_max_iter_zero_raises(): surpyval.Turnbull.fit( np.array([1.0, 2.0, 3.0]), c=np.array([0, 1, 0]), max_iter=0 ) + + +def test_only_the_km_option_is_the_npmle(): + # The three turnbull_estimator options are not three ways of computing + # one number: the EM recovers the same r and d, and they then differ in + # how those become a survival curve. Only Kaplan-Meier is the + # non-parametric MLE; Nelson-Aalen and Fleming-Harrington are exp(-H) + # constructions that are not maximising anything. + # + # This pins the figures quoted in the ``turnbull_estimator`` docstring, + # which exist because comparing a default (FH) Turnbull fit against + # KaplanMeier and reading the gap as a defect is an easy mistake -- it + # is the mistake #260 was originally filed on. + from scipy.optimize import minimize + + x = np.array([2.0, 3.0, 3.0, 4.0, 5.0, 6.0]) + tl = np.array([0.0, 0.0, 1.0, 1.0, 2.0, 2.0]) + support = np.array([2.0, 3.0, 4.0, 5.0, 6.0]) + + def neg_ll(u): + # masses on the support, softmax-parameterised so they stay simplex + p = np.exp(u - u.max()) + p = p / p.sum() + total = 0.0 + for xi, ti in zip(x, tl): + mass = p[support == xi].sum() + at_risk = p[support > ti].sum() # (entry, exit] + if mass <= 0 or at_risk <= 0: + return 1e6 + total += np.log(mass) - np.log(at_risk) + return -total + + best = None + for seed in range(6): + res = minimize( + neg_ll, + np.random.default_rng(seed).normal(size=support.size), + method="Nelder-Mead", + options={"maxiter": 40000, "fatol": 1e-14, "xatol": 1e-12}, + ) + if best is None or res.fun < best.fun: + best = res + + p = np.exp(best.x - best.x.max()) + p = p / p.sum() + npmle_sf2 = p[support > 2.0].sum() + assert npmle_sf2 == pytest.approx(0.75, abs=1e-5) + + got = { + est: float( + np.ravel( + surpyval.Turnbull.fit(x, tl=tl, turnbull_estimator=est).sf(2.0) + )[0] + ) + for est in ("Kaplan-Meier", "Fleming-Harrington", "Nelson-Aalen") + } + + # KM reaches the NPMLE; the other two are elsewhere, and are ordered. + assert got["Kaplan-Meier"] == pytest.approx(npmle_sf2, abs=1e-5) + assert got["Fleming-Harrington"] == pytest.approx(0.7652, abs=1e-3) + assert got["Nelson-Aalen"] == pytest.approx(0.7788, abs=1e-3) + assert ( + got["Kaplan-Meier"] + < got["Fleming-Harrington"] + < got["Nelson-Aalen"] + ) diff --git a/surpyval/univariate/nonparametric/nonparametric_fitter.py b/surpyval/univariate/nonparametric/nonparametric_fitter.py index ecec6d77..e9011f8f 100755 --- a/surpyval/univariate/nonparametric/nonparametric_fitter.py +++ b/surpyval/univariate/nonparametric/nonparametric_fitter.py @@ -117,6 +117,26 @@ def fit( KM, NA, or FH estimator with the Turnbull estimates of r, and d. Defaults to FH. + **This default is why a Turnbull fit does not equal a + KaplanMeier fit on data both can handle.** The Turnbull EM + recovers the same ``r`` and ``d``; the three options then differ + in how they turn those into a survival curve, so the difference + is the estimator, not the data or the EM. On + ``x=[2,3,3,4,5,6], tl=[0,0,1,1,2,2]`` the survival at 2 is + 0.750 under KM, 0.765 under FH and 0.779 under NA. Pass + ``turnbull_estimator='Kaplan-Meier'`` to compare like with like + -- it then agrees with :code:`KaplanMeier` to around 1e-9 on + both ``sf`` and ``cb``, on right-censored and left-truncated + data alike. + + Only the KM option is the non-parametric MLE. NA and FH are + ``exp(-H)`` constructions and are not trying to maximise the + likelihood: on the data above the direct NPMLE of the truncated + likelihood is 0.750, and FH's 0.765 scores worse on it by + design. FH is the default because it behaves better than KM in + the far tails and on zero-inflated data (see the v0.8.0 notes), + not because it is the maximum likelihood answer. + tol : float, optional Turnbull only. The EM stops once the largest change in any interval's probability mass falls below this. Defaults to 1e-10. From 5798d94cd2f2e973be6f0fc3dc710cfd875dfb7d Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 4 Aug 2026 09:11:57 +0000 Subject: [PATCH 02/13] Reformat the new assertion to black's layout The local check was piped through tail, so the pipeline reported tail's exit status and the failure was invisible. Exactly the pipefail hazard described a few commits ago; running black with its own exit status surfaces it. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts --- surpyval/tests/univariate/nonparametric/test_turnbull.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/surpyval/tests/univariate/nonparametric/test_turnbull.py b/surpyval/tests/univariate/nonparametric/test_turnbull.py index 261e42c7..dfe487ef 100644 --- a/surpyval/tests/univariate/nonparametric/test_turnbull.py +++ b/surpyval/tests/univariate/nonparametric/test_turnbull.py @@ -429,7 +429,5 @@ def neg_ll(u): assert got["Fleming-Harrington"] == pytest.approx(0.7652, abs=1e-3) assert got["Nelson-Aalen"] == pytest.approx(0.7788, abs=1e-3) assert ( - got["Kaplan-Meier"] - < got["Fleming-Harrington"] - < got["Nelson-Aalen"] + got["Kaplan-Meier"] < got["Fleming-Harrington"] < got["Nelson-Aalen"] ) From ac43467f545a90f0cf197e0495b2385818c91b00 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 4 Aug 2026 11:40:17 +0000 Subject: [PATCH 03/13] Correct nine wrong examples in the distribution docstrings MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit pytest --doctest-modules over distributions/ gives 39 failures. Thirty are the numpy-2 scalar repr and are cosmetic. Nine were not. Uniform.ff's example called Uniform.sf, and ExpoWeibull.cs's called ExpoWeibull.sf. In both the printed values were correct for the function being documented and wrong for the one being called, so each example read as though the two functions agreed. LogLogistic.sf carried values from some other parameterisation (0.622 where the answer is 0.988), LogLogistic.mean(3, 4) claimed 3 against 3.3322 — the closed form is alpha (pi/beta) / sin(pi/beta) — and Exponential.qf had stale digits. The remaining four were the CustomDistribution Gompertz walkthrough, whose multi-line def used >>> where doctest needs ..., so pasting it raised IndentationError. In every case the code was right and the docs wrong, which is the reassuring direction, but a reader checking their understanding against them would have been misled. They accumulated because the doctests are not run; whether to run them in CI is a separate style decision and is left alone here. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts --- docs/changelog.rst | 28 +++++++++++++++++++ .../distributions/custom_distribution.py | 10 +++---- .../parametric/distributions/expo_weibull.py | 2 +- .../parametric/distributions/exponential.py | 2 +- .../parametric/distributions/loglogistic.py | 4 +-- .../parametric/distributions/uniform.py | 2 +- 6 files changed, 38 insertions(+), 10 deletions(-) diff --git a/docs/changelog.rst b/docs/changelog.rst index 692f77b5..0f547b94 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -4,6 +4,34 @@ Changelog v0.19.1 (unreleased) -------------------- +- **Nine wrong examples in the distribution docstrings.** Running + ``pytest --doctest-modules`` over ``distributions/`` gives 39 failures. + Thirty are the numpy-2 scalar repr (``np.float64(2.719...)`` against a + recorded ``2.719...``) and are cosmetic. Nine were not. + + Five documented outputs were simply wrong. ``Uniform.ff``'s example + called ``Uniform.sf``, and ``ExpoWeibull.cs``'s called + ``ExpoWeibull.sf`` -- in both the printed values were right for the + function being documented and wrong for the one being called, so the + example read as if the two were the same. ``LogLogistic.sf`` carried + values from some other parameterisation entirely (0.622 where the + answer is 0.988), ``LogLogistic.mean(3, 4)`` claimed ``3`` against + ``3.3322`` (the closed form is ``alpha (pi/beta) / sin(pi/beta)``), and + ``Exponential.qf`` had stale digits. + + The other four were the ``CustomDistribution`` example -- the Gompertz + walkthrough -- whose multi-line ``def`` used ``>>>`` where doctest + needs ``...``, so pasting it raised ``IndentationError``. + + In every case the code was right and the documentation was wrong, which + is the reassuring direction, but a reader checking their understanding + against these would have been misled. They accumulated precisely + because the doctests are not run. + + Whether to run doctests in CI, and what to do about the 30 cosmetic + repr differences, is tracked separately -- it is a style decision, not + a correctness one. + - **Documented why a Turnbull fit does not equal a Kaplan-Meier fit.** ``Turnbull.fit`` defaults to ``turnbull_estimator="Fleming-Harrington"`` while ``KaplanMeier.fit`` is, unsurprisingly, KM. The Turnbull EM diff --git a/surpyval/univariate/parametric/distributions/custom_distribution.py b/surpyval/univariate/parametric/distributions/custom_distribution.py index 745fbf9b..4a3daa03 100644 --- a/surpyval/univariate/parametric/distributions/custom_distribution.py +++ b/surpyval/univariate/parametric/distributions/custom_distribution.py @@ -43,16 +43,16 @@ class CustomDistribution(ParametricFitter): >>> name = 'Gompertz' >>> >>> def Hf(x, *params): - >>> return params[0] * np.exp(params[1] * x - 1) - >>> + ... return params[0] * np.exp(params[1] * x - 1) + ... >>> param_names = ['nu', 'b'] >>> bounds = ((0, None), (0, None)) >>> support = (-np.inf, np.inf) >>> Gompertz = surv.CustomDistribution( - name, Hf, param_names, bounds, support - ) + ... name, Hf, param_names, bounds, support + ... ) >>> x = np.array([1, 2, 3, 4, 5]) - >>> Gompertz.fit(x) + >>> model = Gompertz.fit(x) """ def __init__(self, name, fun, param_names, bounds, support): diff --git a/surpyval/univariate/parametric/distributions/expo_weibull.py b/surpyval/univariate/parametric/distributions/expo_weibull.py index 2cc80eee..bd7c337e 100755 --- a/surpyval/univariate/parametric/distributions/expo_weibull.py +++ b/surpyval/univariate/parametric/distributions/expo_weibull.py @@ -184,7 +184,7 @@ def cs(self, x, X, alpha, beta, mu): >>> import numpy as np >>> from surpyval import ExpoWeibull >>> x = np.array([1, 2, 3, 4, 5]) - >>> ExpoWeibull.sf(x, 1, 3, 4, 1.2) + >>> ExpoWeibull.cs(x, 1, 3, 4, 1.2) array([8.77367129e-01, 4.25451775e-01, 5.09266354e-02, 5.37452200e-04, 1.35732908e-07]) """ diff --git a/surpyval/univariate/parametric/distributions/exponential.py b/surpyval/univariate/parametric/distributions/exponential.py index 7ea07526..204c8f9d 100755 --- a/surpyval/univariate/parametric/distributions/exponential.py +++ b/surpyval/univariate/parametric/distributions/exponential.py @@ -327,7 +327,7 @@ def qf(self, p, failure_rate): >>> from surpyval import Exponential >>> p = np.array([.1, .2, .3, .4, .5]) >>> Exponential.qf(p, 3) - array([0.03512219, 0.07438118, 0.11889152, 0.17027853, 0.23104906]) + array([0.03512017, 0.07438118, 0.11889165, 0.17027521, 0.23104906]) """ return -np.log1p(-p) / failure_rate diff --git a/surpyval/univariate/parametric/distributions/loglogistic.py b/surpyval/univariate/parametric/distributions/loglogistic.py index 597b6dca..cb3c9569 100755 --- a/surpyval/univariate/parametric/distributions/loglogistic.py +++ b/surpyval/univariate/parametric/distributions/loglogistic.py @@ -57,7 +57,7 @@ def sf(self, x, alpha, beta): >>> from surpyval import LogLogistic >>> x = np.array([1, 2, 3, 4, 5]) >>> LogLogistic.sf(x, 3, 4) - array([0.62245933, 0.5621765 , 0.5 , 0.4378235 , 0.37754067]) + array([0.98780488, 0.83505155, 0.5 , 0.24035608, 0.11473088]) """ # 1 / (1 + (x/alpha)^beta): algebraically identical to the # (x/alpha)^-beta form but defined at x = 0 (sf(0) = 1) instead @@ -304,7 +304,7 @@ def mean(self, alpha, beta): -------- >>> from surpyval import LogLogistic >>> LogLogistic.mean(3, 4) - 3 + 3.332162203618775 """ if beta > 1: return (alpha * np.pi / beta) / (np.sin(np.pi / beta)) diff --git a/surpyval/univariate/parametric/distributions/uniform.py b/surpyval/univariate/parametric/distributions/uniform.py index e5fb9dd0..095ada9c 100755 --- a/surpyval/univariate/parametric/distributions/uniform.py +++ b/surpyval/univariate/parametric/distributions/uniform.py @@ -119,7 +119,7 @@ def ff(self, x, a, b): >>> import numpy as np >>> from surpyval import Uniform >>> x = np.array([1, 2, 3, 4, 5]) - >>> Uniform.sf(x, 0, 6) + >>> Uniform.ff(x, 0, 6) array([0.16666667, 0.33333333, 0.5 , 0.66666667, 0.83333333]) """ f = np.zeros_like(x) From c11cfedefa7b4a528a2de66fd7950c16c646b5c9 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 4 Aug 2026 11:54:29 +0000 Subject: [PATCH 04/13] Stop Gamma offering a probability plot it cannot draw MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Gamma.fit(x, how="MPP") now raises, joining Beta and ExpoWeibull which already declined for the same reason. A probability plot rearranges the survival function so a transform of the data falls on a straight line. Weibull gives log(-log S) = beta log x - beta log alpha: the axes do not depend on the answer. The Gamma has no such rearrangement — its CDF is the regularised incomplete gamma and the shape sits inside that special function rather than outside as an exponent. The only straight-line y-axis is the inverse incomplete gamma, which needs the shape. To draw the axis you need the answer; to get the answer you need the axis. The code broke the circle by guessing the shape from moments, drawing the plot on the guess and regressing. A wrong guess means a wrong axis, points that are no longer straight on it, and a line fitted through a curve — a confident wrong estimate rather than an error. An offset made it worse, since the shift distorts the low-x end hardest and that is where the shape information is. plot() is untouched: it transforms with the fitted parameters, so the axis is correct by the time it is drawn. MLE, MSE and MOM are unchanged, offset included. Deleting the 118-line Gamma.mpp override takes the rr="x" mis-inversion and the censored-data LinAlgError from #257 with it, by making both unreachable. The three tests that asserted Gamma MPP recovery now assert the refusal and keep their offset-recovery coverage under MLE, and the test_fit.py MPP sweep gates on supports_mpp rather than a hardcoded list so a new distribution will not need it edited. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts --- docs/changelog.rst | 33 ++++ .../parametric/test_distribution_fixes.py | 19 ++- .../tests/univariate/parametric/test_fit.py | 5 +- .../parametric/test_offset_divergence.py | 23 ++- .../parametric/distributions/gamma.py | 141 +++--------------- 5 files changed, 90 insertions(+), 131 deletions(-) diff --git a/docs/changelog.rst b/docs/changelog.rst index 0f547b94..4a643cf4 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -4,6 +4,39 @@ Changelog v0.19.1 (unreleased) -------------------- +- **``Gamma`` no longer offers probability plotting as a fit method.** + ``Gamma.fit(x, how="MPP")`` now raises, joining ``Beta`` and + ``ExpoWeibull``, which already declined for the same reason. + + A probability plot works by rearranging the survival function so some + transform of the data falls on a straight line. For a Weibull, + ``log(-log S) = beta log x - beta log alpha`` — the axes do not depend + on the answer, so you can draw them before knowing anything. The + Gamma has no such rearrangement: its CDF is the regularised incomplete + gamma function and the shape sits *inside* that special function + rather than outside as an exponent. The only straight-line y-axis is + the inverse incomplete gamma, which needs the shape. To draw the axis + you need the answer; to get the answer you need the axis. + + The code broke the circle by guessing the shape from moments, drawing + the plot on that guess, and regressing. When the guess was off, the + axis was the wrong axis, the points were no longer straight on it, and + the regression fitted a line through a curve — returning a confident + wrong estimate rather than an error. An offset made it worse: the + shift distorts the low-``x`` end hardest, which is exactly where the + shape information lives. + + ``plot()`` is unaffected. It transforms with the *fitted* parameters, + so by the time the plot is drawn the axis is the right one — the + probability plot of an MLE-fitted Gamma remains a valid diagnostic. + Fitting is unchanged for MLE (the default), MSE and MOM. + + The 118-line ``Gamma.mpp`` override is deleted with it, which removes + the ``rr="x"`` mis-inversion and the censored-data ``LinAlgError`` from + #257 by making both paths unreachable. The MPP sweep in ``test_fit.py`` + now gates on each distribution's ``supports_mpp`` flag instead of a + hardcoded exclusion list, so it stays correct without editing. + - **Nine wrong examples in the distribution docstrings.** Running ``pytest --doctest-modules`` over ``distributions/`` gives 39 failures. Thirty are the numpy-2 scalar repr (``np.float64(2.719...)`` against a diff --git a/surpyval/tests/univariate/parametric/test_distribution_fixes.py b/surpyval/tests/univariate/parametric/test_distribution_fixes.py index 54580272..b8ffe64c 100644 --- a/surpyval/tests/univariate/parametric/test_distribution_fixes.py +++ b/surpyval/tests/univariate/parametric/test_distribution_fixes.py @@ -44,18 +44,31 @@ def test_exponential_offset_mpp_rr_x_inversion(): def test_gamma_censored_mpp_rr_x_does_not_crash(): + # Was a LinAlgError: rr="x" regressed the filtered y against the raw + # x whenever censoring filtered any point (#257). Gamma declines MPP + # entirely now (#158), so the crashing path is unreachable -- the + # refusal arrives before any regression is attempted. np.random.seed(1) x = Gamma.random(300, 3, 2) c = (x > 2.5).astype(int) - m = Gamma.fit(np.minimum(x, 2.5), c=c, how="MPP", rr="x") - assert np.all(np.isfinite(m.params)) + with pytest.raises(ValueError, match="does not work with probability"): + Gamma.fit(np.minimum(x, 2.5), c=c, how="MPP", rr="x") @pytest.mark.parametrize("rr", ["x", "y"]) def test_gamma_offset_mpp_recovers_offset(rr): + # Gamma no longer offers MPP at all (#158): the probability plot's + # own y-axis is the inverse incomplete gamma, which needs the shape + # being estimated, so the fit regressed against an axis built from a + # guess and returned a confident wrong answer. It now refuses. np.random.seed(2) x = Gamma.random(300, 3, 2) + 10 - m = Gamma.fit(x, how="MPP", rr=rr, offset=True) + with pytest.raises(ValueError, match="does not work with probability"): + Gamma.fit(x, how="MPP", rr=rr, offset=True) + + # The offset recovery this test existed to protect still holds under + # the estimators Gamma does support. + m = Gamma.fit(x, offset=True) assert m.gamma == pytest.approx(10.0, abs=1.5) assert m.params[0] == pytest.approx(3.0, rel=0.5) diff --git a/surpyval/tests/univariate/parametric/test_fit.py b/surpyval/tests/univariate/parametric/test_fit.py index cc7e5706..b2d63f20 100644 --- a/surpyval/tests/univariate/parametric/test_fit.py +++ b/surpyval/tests/univariate/parametric/test_fit.py @@ -246,7 +246,10 @@ def test_mle_convergence_small(dist, random_parameters, kind): "dist,random_parameters,rr", generate_mpp_test_cases(), ids=idfunc ) def test_mpp(dist, random_parameters, rr): - if dist not in [Beta, ExpoWeibull]: + # Gate on the distribution's own flag rather than a hardcoded + # list: Beta, ExpoWeibull and Gamma all lack a linearising + # probability plot, and a new one should not need this edited. + if dist.supports_mpp: for n in FIT_SIZES: test_params = [] tol = 0.025 diff --git a/surpyval/tests/univariate/parametric/test_offset_divergence.py b/surpyval/tests/univariate/parametric/test_offset_divergence.py index fe14c938..18284712 100644 --- a/surpyval/tests/univariate/parametric/test_offset_divergence.py +++ b/surpyval/tests/univariate/parametric/test_offset_divergence.py @@ -13,6 +13,7 @@ """ import numpy as np +import pytest from scipy.stats import wasserstein_distance from surpyval import Gamma, Rayleigh @@ -64,19 +65,29 @@ def _summary(dist, true_params, how, seed=0): } -def test_mpp_offset_gamma_parameters_recovered(): +def test_mpp_offset_gamma_is_refused(): """MPP offset on Gamma used to land on an absurd parameter tuple (huge shape, gamma far off) that merely mimicked the true distribution. This was the divergence this module existed to pin - down; #257 fixed the initialiser (multi-started shape search) and - the rr="x" inversion, so the *parameters* are now recovered too.""" - s = _summary(Gamma, (3.0, 2.0), how="MPP") + down. #257 fixed the initialiser and the rr="x" inversion, which + recovered the parameters, but the underlying circularity remained: + the plot's y-axis is the inverse incomplete gamma, so it has to be + drawn from a guess at the very shape being estimated. Gamma now + declines MPP outright (#158) rather than returning a confident + answer off a wrong axis.""" + x = Gamma.random(500, 3.0, 2.0) + TRUE_GAMMA + with pytest.raises(ValueError, match="does not work with probability"): + Gamma.fit(x, how="MPP", offset=True) + + +def test_mle_offset_gamma_parameters_recovered(): + """The recovery the MPP test above was really protecting, under an + estimator Gamma supports.""" + s = _summary(Gamma, (3.0, 2.0), how="MLE") - # The parameters are now close to the truth. assert abs(s["fit"].gamma - TRUE_GAMMA) < 1.0, s["fit"].gamma assert s["max_param_rel_err"] < 0.20 # within 20% - # ...and the distribution remains essentially exact. assert s["mean_rel_err"] < 0.01 # mean within 1% assert s["median_rel_err"] < 0.03 # median within 3% assert s["std_rel_err"] < 0.10 # spread within 10% diff --git a/surpyval/univariate/parametric/distributions/gamma.py b/surpyval/univariate/parametric/distributions/gamma.py index 51a32c8e..5be16717 100755 --- a/surpyval/univariate/parametric/distributions/gamma.py +++ b/surpyval/univariate/parametric/distributions/gamma.py @@ -1,13 +1,8 @@ -import warnings - from autograd.scipy.special import gamma as agamma from autograd.scipy.special import gammaln as agammaln -from scipy.optimize import minimize from scipy.special import digamma, gammaincinv -from scipy.stats import pearsonr from surpyval import np -from surpyval.univariate.nonparametric import plotting_positions from surpyval.univariate.parametric.parametric_fitter import ParametricFitter from surpyval.utils.autograd_gamma_compat import gammainc as agammainc from surpyval.utils.autograd_gamma_compat import gammainccln as agammainccln @@ -35,6 +30,26 @@ def __init__(self, name): param_map={"alpha": 0, "beta": 1}, plot_x_scale="linear", ) + # The Gamma has no linearising probability plot, for the same + # reason as the Beta above it: the CDF is the regularised + # incomplete gamma function, and the shape sits *inside* that + # special function rather than outside it as an exponent. The + # only straight-line y-axis is the inverse incomplete gamma, + # which needs the shape -- so to draw the axis you need the + # answer, and to get the answer you need the axis. + # + # MPP broke the circle by guessing the shape from moments, + # drawing the plot on that guess and regressing. When the guess + # is off the axis is the wrong axis, the points are no longer + # straight on it, and the regression fits a line through a + # curve -- returning a confident, wrong estimate rather than an + # error. An offset makes it worse: the shift distorts the low-x + # end hardest, which is exactly where the shape information is. + # + # Fit by MLE (the default), MSE or MOM instead. ``plot()`` still + # works, because it transforms with the *fitted* parameters, so + # the axis is the right one by the time it is drawn. + self.supports_mpp = False @staticmethod def _moment_estimate(x): @@ -490,121 +505,5 @@ def mpp_inv_y_transform(self, y, *params): def mpp_x_transform(self, x, gamma=0): return x - gamma - def mpp( - self, - x, - c=None, - n=None, - t=None, - heuristic="Nelson-Aalen", - rr="y", - on_d_is_0=False, - offset=False, - ): - # Forward the truncation windows (previously dropped, #280). - x_pp, r, d, F = plotting_positions( - x, c=c, n=n, t=t, heuristic=heuristic - ) - - results = {} - - if on_d_is_0: - pass - else: - F = F[d > 0] - x_pp = x_pp[d > 0] - - if (F == 1).any(): - mask = F != 1 - warnings.warn( - "Some heuristic values for CDF = 1 have been " - "encountered in plotting points and have been " - "ignored.", - stacklevel=2, - ) - F = F[mask] - x_pp = x_pp[mask] - - init = self._parameter_initialiser(x_pp, c, n) - - mask = np.isfinite(F) - if not mask.all(): - warnings.warn( - "Some Infinite values encountered in plotting " - "points and have been ignored.", - stacklevel=2, - ) - F = F[mask] - x_pp = x_pp[mask] - - if offset: - - def fun(a): - return -pearsonr(x_pp, self.mpp_y_transform(F, a, 1.0))[0] - - # The moment-based init is computed from the *unshifted* data, - # which diverges for strongly offset data (tiny relative spread - # -> huge alpha) and strands the correlation search at a bad - # local optimum. Try several starts and keep the best (#257). - starts = {float(init[0]), 0.5, 1.0, 2.0, 5.0} - best = None - for a0 in starts: - res = minimize(fun, [a0], bounds=((1e-8, None),)) - if best is None or res.fun < best.fun: - best = res - alpha = best.x[0] - - y_pp = self.mpp_y_transform(F, alpha) - - if rr == "y": - # y = beta * (x - gamma): slope is beta, intercept is - # -beta * gamma. - params = np.polyfit(x_pp, y_pp, 1) - beta = params[0] - gamma = -params[1] / beta - elif rr == "x": - # x = y / beta + gamma: slope is 1/beta and the intercept - # IS gamma (#257). - params = np.polyfit(y_pp, x_pp, 1) - beta = 1.0 / params[0] - gamma = params[1] - - results["gamma"] = gamma - results["params"] = np.array([alpha, beta]) - - return results - else: - if rr == "y": - x_pp = x_pp[:, np.newaxis] - - def fun(alpha): - y_pp = self.mpp_y_transform(F, alpha, 1.0) - return np.linalg.lstsq(x_pp, y_pp)[1] - - res = minimize(fun, init[0], bounds=((0, None),)) - alpha = res.x[0] - y_pp = self.mpp_y_transform(F, alpha, 1.0) - beta, residuals, _, _ = np.linalg.lstsq(x_pp, y_pp) - beta = beta[0] - else: - # Regress against the same filtered plotting positions used - # everywhere else — the raw ``x`` argument has a different - # length whenever any point was filtered (#257). - - def fun(a): - y = self.mpp_y_transform(F, a, 1.0)[:, np.newaxis] - return np.linalg.lstsq(y, x_pp)[1] - - res = minimize(fun, init[0], bounds=((0, None),)) - alpha = res.x[0] - beta = np.linalg.lstsq( - self.mpp_y_transform(F, alpha, 1.0)[:, np.newaxis], x_pp - )[0][0] - beta = 1.0 / beta - - results["params"] = np.array([alpha, beta]) - - return results - Gamma: ParametricFitter = Gamma_("Gamma") From e949f33bfb48b3faa1b955478c3c5dd5226e9c7f Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 4 Aug 2026 12:13:40 +0000 Subject: [PATCH 05/13] Record the real output in the distribution docstring examples MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit pytest --doctest-modules over distributions/ is now green: 139 examples, no failures, down from 39 failing. Most were the numpy 2 scalar repr. Weibull.mean(3, 4) prints np.float64(2.7192074311664314); the docstring recorded the bare float numpy 1 used to print. The examples now carry the wrapper because that is what a reader sees at their own prompt. The alternative was np.set_printoptions(legacy="1.25") in a fixture, which keeps the docstrings prettier by showing people output their session will not produce — prettier, but not true. Four qf examples exceeded 79 columns once the real output was recorded, numpy having rewrapped the arrays differently from the hand-wrapping. Rather than reflow them into something numpy would never emit, those examples now take fewer probabilities, so what is printed is exactly what that input produces. Two scalar examples had drifted in the last digit and are re-recorded against a direct run. Every rewrite was gated on the numbers agreeing to 1e-12 first, so a genuine mismatch could not be blessed by the sweep; nothing failed that gate. Two files needed hand correction afterwards where mean and moment share a value and the text search matched the wrong occurrence. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts --- docs/changelog.rst | 24 +++++++++++++++++++ .../parametric/distributions/beta.py | 4 ++-- .../parametric/distributions/beta4.py | 2 +- .../parametric/distributions/binomial.py | 14 +++++------ .../parametric/distributions/exponential.py | 4 ++-- .../parametric/distributions/gamma.py | 4 ++-- .../parametric/distributions/gumbel.py | 7 +++--- .../parametric/distributions/gumbel_lev.py | 2 +- .../parametric/distributions/logistic.py | 6 ++--- .../parametric/distributions/loglogistic.py | 4 ++-- .../parametric/distributions/lognormal.py | 11 ++++----- .../parametric/distributions/normal.py | 8 +++---- .../parametric/distributions/rayleigh.py | 4 ++-- .../parametric/distributions/uniform.py | 4 ++-- .../parametric/distributions/weibull.py | 4 ++-- 15 files changed, 62 insertions(+), 40 deletions(-) diff --git a/docs/changelog.rst b/docs/changelog.rst index 4a643cf4..0794a3c1 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -4,6 +4,30 @@ Changelog v0.19.1 (unreleased) -------------------- +- **The distribution docstring examples now show what you actually see.** + ``pytest --doctest-modules`` over ``distributions/`` is green: 139 + examples, no failures. Previously 39 failed. + + Most were the numpy 2 scalar repr. ``Weibull.mean(3, 4)`` prints + ``np.float64(2.7192074311664314)``, where the docstring recorded the + bare ``2.7192074311664314`` that numpy 1 used to print. The examples + now record the wrapper, because that is what appears at a prompt -- + the alternative was ``np.set_printoptions(legacy="1.25")`` in a test + fixture, which would have kept the docstrings prettier by showing + readers something their own session will not produce. + + Four ``qf`` examples printed wider than the 79-column limit once the + real output was recorded, numpy having rewrapped the arrays. Rather + than hand-wrap them into something numpy would not emit, those + examples take fewer probabilities: what is shown is exactly what that + input produces. + + Two scalar examples had drifted in the last digit, and are re-recorded. + + With the examples now true, ``--doctest-modules`` is worth running in + CI -- that decision is separate, and is what would stop this + recurring. + - **``Gamma`` no longer offers probability plotting as a fit method.** ``Gamma.fit(x, how="MPP")`` now raises, joining ``Beta`` and ``ExpoWeibull``, which already declined for the same reason. diff --git a/surpyval/univariate/parametric/distributions/beta.py b/surpyval/univariate/parametric/distributions/beta.py index 670cf66e..b0ccca59 100755 --- a/surpyval/univariate/parametric/distributions/beta.py +++ b/surpyval/univariate/parametric/distributions/beta.py @@ -333,7 +333,7 @@ def moment(self, n, alpha, beta): -------- >>> from surpyval import Beta >>> Beta.moment(2, 3, 4) - 0.21428571428571427 + np.float64(0.2142857142857143) """ return np.exp(abetaln(n + alpha, beta) - abetaln(alpha, beta)) @@ -369,7 +369,7 @@ def entropy(self, alpha, beta): -------- >>> from surpyval import Beta >>> Beta.entropy(3, 4) - -0.3443445622221013 + np.float64(-0.3443445622221013) """ return ( abetaln(alpha, beta) diff --git a/surpyval/univariate/parametric/distributions/beta4.py b/surpyval/univariate/parametric/distributions/beta4.py index 815d5a66..3b2a6f89 100644 --- a/surpyval/univariate/parametric/distributions/beta4.py +++ b/surpyval/univariate/parametric/distributions/beta4.py @@ -411,7 +411,7 @@ def moment(self, m, alpha, beta, a, b): -------- >>> from surpyval import Beta4 >>> Beta4.moment(1, 3, 4, 2, 3) - 2.4285714285714284 + np.float64(2.428571428571429) """ scale = b - a total = 0.0 diff --git a/surpyval/univariate/parametric/distributions/binomial.py b/surpyval/univariate/parametric/distributions/binomial.py index 5038e049..ee6b8532 100644 --- a/surpyval/univariate/parametric/distributions/binomial.py +++ b/surpyval/univariate/parametric/distributions/binomial.py @@ -64,7 +64,7 @@ def df(self, x, n, p): -------- >>> from surpyval import Binomial >>> Binomial.df(2, 5, 0.3) - 0.3086999999999998 + np.float64(0.3086999999999998) """ return binom.pmf(x, n, p) @@ -97,7 +97,7 @@ def ff(self, x, n, p): -------- >>> from surpyval import Binomial >>> Binomial.ff(2, 5, 0.3) - 0.83692 + np.float64(0.83692) """ return binom.cdf(x, n, p) @@ -129,7 +129,7 @@ def sf(self, x, n, p): -------- >>> from surpyval import Binomial >>> Binomial.sf(2, 5, 0.3) - 0.16308 + np.float64(0.16308) """ return binom.sf(x, n, p) @@ -215,7 +215,7 @@ def qf(self, q, n, p): -------- >>> from surpyval import Binomial >>> Binomial.qf(0.5, 5, 0.3) - 1.0 + np.float64(1.0) """ return binom.ppf(q, n, p) @@ -289,7 +289,7 @@ def moment(self, m, n, p): -------- >>> from surpyval import Binomial >>> Binomial.moment(1, 5, 0.3) - 1.5 + np.float64(1.5) """ return binom.moment(m, n, p) @@ -302,7 +302,7 @@ def entropy(self, n, p): -------- >>> from surpyval import Binomial >>> Binomial.entropy(5, 0.3) - 1.413614855283445 + np.float64(1.413614855283445) """ return binom.entropy(n, p) @@ -414,7 +414,7 @@ def from_params(self, params): >>> from surpyval import Binomial >>> model = Binomial.from_params([5, 0.3]) >>> model.mean() - 1.5 + np.float64(1.5) """ params = np.atleast_1d(np.asarray(params, dtype=float)) diff --git a/surpyval/univariate/parametric/distributions/exponential.py b/surpyval/univariate/parametric/distributions/exponential.py index 204c8f9d..f760ebcd 100755 --- a/surpyval/univariate/parametric/distributions/exponential.py +++ b/surpyval/univariate/parametric/distributions/exponential.py @@ -386,7 +386,7 @@ def moment(self, n, failure_rate): -------- >>> from surpyval import Exponential >>> Exponential.moment(2, 3) - 0.2222222222222222 + np.float64(0.2222222222222222) """ return factorial(n) / (failure_rate**n) @@ -414,7 +414,7 @@ def entropy(self, failure_rate): -------- >>> from surpyval import Exponential >>> Exponential.entropy(3) - -0.09861228866810978 + np.float64(-0.09861228866810978) """ return 1 - np.log(failure_rate) diff --git a/surpyval/univariate/parametric/distributions/gamma.py b/surpyval/univariate/parametric/distributions/gamma.py index 5be16717..2f7450d1 100755 --- a/surpyval/univariate/parametric/distributions/gamma.py +++ b/surpyval/univariate/parametric/distributions/gamma.py @@ -411,7 +411,7 @@ def moment(self, n, alpha, beta): -------- >>> from surpyval import Gamma >>> Gamma.moment(3, 3, 4) - 0.9375 + np.float64(0.9375) """ return agamma(n + alpha) / (beta**n * agamma(alpha)) @@ -445,7 +445,7 @@ def entropy(self, alpha, beta): -------- >>> from surpyval import Gamma >>> Gamma.entropy(3, 4) - 0.46128414924312033 + np.float64(0.46128414924312033) """ return ( alpha diff --git a/surpyval/univariate/parametric/distributions/gumbel.py b/surpyval/univariate/parametric/distributions/gumbel.py index efd8aceb..4dc5bbd7 100755 --- a/surpyval/univariate/parametric/distributions/gumbel.py +++ b/surpyval/univariate/parametric/distributions/gumbel.py @@ -243,10 +243,9 @@ def qf(self, p, mu, sigma): -------- >>> import numpy as np >>> from surpyval import Gumbel - >>> p = np.array([.1, .2, .3, .4, .5]) + >>> p = np.array([0.1, 0.3, 0.5]) >>> Gumbel.qf(p, 3, 2) - array([-1.50073465e+00, 1.20026481e-04, 9.38139134e-01, 1.65654602e+00, - 2.26697416e+00]) + array([-1.50073465, 0.93813913, 2.26697416]) """ return mu + sigma * (np.log(-np.log1p(-p))) @@ -325,7 +324,7 @@ def entropy(self, mu, sigma): -------- >>> from surpyval import Gumbel >>> Gumbel.entropy(3, 2) - 2.270362845461478 + np.float64(2.270362845461478) """ return np.log(sigma) + euler_gamma + 1 diff --git a/surpyval/univariate/parametric/distributions/gumbel_lev.py b/surpyval/univariate/parametric/distributions/gumbel_lev.py index 352cfc76..8b5f3639 100644 --- a/surpyval/univariate/parametric/distributions/gumbel_lev.py +++ b/surpyval/univariate/parametric/distributions/gumbel_lev.py @@ -335,7 +335,7 @@ def entropy(self, mu, sigma): -------- >>> from surpyval import GumbelLEV >>> GumbelLEV.entropy(3, 2) - 2.270362845461478 + np.float64(2.270362845461478) """ return np.log(sigma) + euler_gamma + 1 diff --git a/surpyval/univariate/parametric/distributions/logistic.py b/surpyval/univariate/parametric/distributions/logistic.py index 88924a63..53cdd424 100755 --- a/surpyval/univariate/parametric/distributions/logistic.py +++ b/surpyval/univariate/parametric/distributions/logistic.py @@ -226,9 +226,9 @@ def qf(self, p, mu, sigma): -------- >>> import numpy as np >>> from surpyval import Logistic - >>> p = np.array([.1, .2, .3, .4, .5]) + >>> p = np.array([0.1, 0.2, 0.3, 0.4]) >>> Logistic.qf(p, 3, 4) - array([-5.78889831, -2.54517744, -0.38919144, 1.37813957, 3. ]) + array([-5.78889831, -2.54517744, -0.38919144, 1.37813957]) """ return mu + sigma * (np.log(p) - np.log1p(-p)) @@ -311,7 +311,7 @@ def entropy(self, mu, sigma): -------- >>> from surpyval import Logistic >>> Logistic.entropy(3, 4) - 3.386294361119891 + np.float64(3.386294361119891) """ return np.log(sigma) + 2 diff --git a/surpyval/univariate/parametric/distributions/loglogistic.py b/surpyval/univariate/parametric/distributions/loglogistic.py index cb3c9569..54b0b215 100755 --- a/surpyval/univariate/parametric/distributions/loglogistic.py +++ b/surpyval/univariate/parametric/distributions/loglogistic.py @@ -304,7 +304,7 @@ def mean(self, alpha, beta): -------- >>> from surpyval import LogLogistic >>> LogLogistic.mean(3, 4) - 3.332162203618775 + np.float64(3.332162203618775) """ if beta > 1: return (alpha * np.pi / beta) / (np.sin(np.pi / beta)) @@ -382,7 +382,7 @@ def entropy(self, alpha, beta): -------- >>> from surpyval import LogLogistic >>> LogLogistic.entropy(3, 4) - 1.7123179275482192 + np.float64(1.7123179275482192) """ return np.log(alpha / beta) + 2 diff --git a/surpyval/univariate/parametric/distributions/lognormal.py b/surpyval/univariate/parametric/distributions/lognormal.py index 36ff05d4..571059a2 100755 --- a/surpyval/univariate/parametric/distributions/lognormal.py +++ b/surpyval/univariate/parametric/distributions/lognormal.py @@ -301,10 +301,9 @@ def qf(self, p, mu, sigma): -------- >>> import numpy as np >>> from surpyval import LogNormal - >>> p = np.array([.1, .2, .3, .4, .5]) + >>> p = np.array([0.1, 0.2, 0.3, 0.4]) >>> LogNormal.qf(p, 3, 4) - array([ 0.11928899, 0.69316658, 2.46550819, 7.29078766, - 20.08553692]) + array([0.11928899, 0.69316658, 2.46550819, 7.29078766]) """ return np.exp(scipy_norm.ppf(p, mu, sigma)) @@ -334,7 +333,7 @@ def mean(self, mu, sigma): -------- >>> from surpyval import LogNormal >>> LogNormal.mean(3, 4) - 59874.14171519782 + np.float64(59874.14171519782) """ return np.exp(mu + (sigma**2) / 2) @@ -366,7 +365,7 @@ def moment(self, n, mu, sigma): -------- >>> from surpyval import LogNormal >>> LogNormal.moment(2, 3, 4) - 3.1855931757113756e+16 + np.float64(3.1855931757113756e+16) """ return np.exp(n * mu + (n**2 * sigma**2) / 2) @@ -396,7 +395,7 @@ def entropy(self, mu, sigma): -------- >>> from surpyval import LogNormal >>> LogNormal.entropy(3, 4) - 5.805232894324563 + np.float64(5.805232894324563) """ return mu + 0.5 * np.log(2 * np.pi * np.e * sigma**2) diff --git a/surpyval/univariate/parametric/distributions/normal.py b/surpyval/univariate/parametric/distributions/normal.py index 3afcd536..8b34429c 100755 --- a/surpyval/univariate/parametric/distributions/normal.py +++ b/surpyval/univariate/parametric/distributions/normal.py @@ -308,9 +308,9 @@ def qf(self, p, mu, sigma): -------- >>> import numpy as np >>> from surpyval import Normal - >>> p = np.array([.1, .2, .3, .4, .5]) + >>> p = np.array([0.1, 0.2, 0.3, 0.4]) >>> Normal.qf(p, 3, 4) - array([-2.12620626, -0.36648493, 0.90239795, 1.98661159, 3. ]) + array([-2.12620626, -0.36648493, 0.90239795, 1.98661159]) """ return scipy_norm.ppf(p, mu, sigma) @@ -372,7 +372,7 @@ def moment(self, n, mu, sigma): -------- >>> from surpyval import Normal >>> Normal.moment(2, 3, 4) - 25.0 + np.float64(25.0) """ return scipy_norm.moment(n, mu, sigma) @@ -402,7 +402,7 @@ def entropy(self, mu, sigma): -------- >>> from surpyval import Normal >>> Normal.entropy(3, 4) - 2.8052328943245635 + np.float64(2.8052328943245635) """ return 0.5 * np.log(2 * np.pi * np.e * sigma**2) diff --git a/surpyval/univariate/parametric/distributions/rayleigh.py b/surpyval/univariate/parametric/distributions/rayleigh.py index 98e72f54..719ac981 100644 --- a/surpyval/univariate/parametric/distributions/rayleigh.py +++ b/surpyval/univariate/parametric/distributions/rayleigh.py @@ -299,7 +299,7 @@ def mean(self, sigma): -------- >>> from surpyval import Rayleigh >>> Rayleigh.mean(3) - 3.7599424119465006 + np.float64(3.7599424119465006) """ return sigma * np.sqrt(np.pi / 2) @@ -329,7 +329,7 @@ def moment(self, n, sigma): -------- >>> from surpyval import Rayleigh >>> Rayleigh.moment(2, 3) - 18.0 + np.float64(18.0) """ return (sigma**n) * (2 ** (n / 2)) * gamma_func(1 + n / 2) diff --git a/surpyval/univariate/parametric/distributions/uniform.py b/surpyval/univariate/parametric/distributions/uniform.py index 095ada9c..608eb615 100755 --- a/surpyval/univariate/parametric/distributions/uniform.py +++ b/surpyval/univariate/parametric/distributions/uniform.py @@ -341,7 +341,7 @@ def moment(self, n, a, b): -------- >>> from surpyval import Uniform >>> Uniform.moment(2, 0, 6) - 12.0 + np.float64(12.0) """ if n == 0: return 1 @@ -377,7 +377,7 @@ def entropy(self, a, b): -------- >>> from surpyval import Uniform >>> Uniform.entropy(0, 6) - 1.791759469228055 + np.float64(1.791759469228055) """ return np.log(b - a) diff --git a/surpyval/univariate/parametric/distributions/weibull.py b/surpyval/univariate/parametric/distributions/weibull.py index 59e8ac5a..f66a6bd0 100755 --- a/surpyval/univariate/parametric/distributions/weibull.py +++ b/surpyval/univariate/parametric/distributions/weibull.py @@ -302,7 +302,7 @@ def mean(self, alpha, beta): -------- >>> from surpyval import Weibull >>> Weibull.mean(3, 4) - 2.7192074311664314 + np.float64(2.7192074311664314) """ return alpha * gamma_func(1 + 1.0 / beta) @@ -334,7 +334,7 @@ def moment(self, n, alpha, beta): -------- >>> from surpyval import Weibull >>> Weibull.moment(2, 3, 4) - 7.976042329074821 + np.float64(7.976042329074821) """ return alpha**n * gamma_func(1 + n / beta) From 7dd534eb72b8ee4cd84eaad4821fbe04038d8ca3 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 4 Aug 2026 12:49:20 +0000 Subject: [PATCH 06/13] Run the docstring examples in CI A docstring example is a promise about what the library prints, and it is the first thing a user or a coding agent reaches for -- help() is faster than opening the docs. Nothing was checking it, so it drifted. Turning on --doctest-modules over the package gave 59 failing tests; this fixes all of them and adds the step to the deployment workflow. Beyond the cosmetic drift (numpy 2 scalar reprs, optimiser output from two rewrites ago), the run found: - Twelve examples that could not run at all. Six regression docstrings (PH, AH, PO, AFT, AcceleratedLife, Frailty) were sketches -- `model = PH(Weibull).fit(x, Z=covariates, c=c)` with none of x, covariates or c ever defined. Four used >>> on the continuation lines of a multi-line call, so pasting them raised SyntaxError. plotting_positions imported from a module that moved several releases ago. ParametricFitter.fit demonstrated how='MPP' on interval-censored input, which now correctly requires the Turnbull heuristic and raises without it. - The five ParametricRegressionModel prediction examples (sf, ff, df, hf, Hf) had been copied from the univariate Parametric class and never adapted: they built a Weibull.from_params([10, 3]) and called it with no covariates, documenting a signature the method does not have. - Parametric.var() claimed 11.229 for a Weibull(10, 3). The variance is 10.533; the code was right. - Several examples fitted unseeded random data and then recorded specific digits. They now seed. Parametric.hf and Parametric.Hf returned a 0-d array for scalar input where sf, ff, df and qf returned a numpy scalar, as did cs -- against their own Returns sections, which promise "the scalar value ... if a scalar was passed". np.where does not collapse a 0-d result; [()] does, and is a no-op on a real array. Two doctest flags are set in pyproject.toml. NORMALIZE_WHITESPACE, because numpy picks its own line breaks and column padding for an array and both move with the widest element. ELLIPSIS, so an example ending in a fit can write 529.05371... rather than all seventeen digits: the trailing digits of an optimiser's output are not part of what the docstring promises, and they move with the BLAS and the platform. Array reprs are left exact -- numpy already prints only eight significant digits there. 229 examples, all passing. Closes #158. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts --- .github/workflows/actions.yml | 17 ++++ docs/changelog.rst | 69 ++++++++++++-- pyproject.toml | 19 ++++ surpyval/degradation/degradation_analysis.py | 4 +- surpyval/recurrent/parametric/cox_lewis.py | 16 ++-- surpyval/recurrent/parametric/crow_amsaa.py | 16 ++-- surpyval/recurrent/parametric/duane.py | 16 ++-- surpyval/recurrent/parametric/hpp.py | 13 ++- .../regression/proportional_intensity.py | 6 +- .../renewal/generalized_one_renewal.py | 30 +++--- .../recurrent/renewal/generalized_renewal.py | 26 +++--- surpyval/univariate/information_criteria.py | 8 +- .../univariate/nonparametric/nonparametric.py | 2 +- .../nonparametric/plotting_positions.py | 2 +- .../univariate/nonparametric/success_run.py | 2 +- .../univariate/parametric/mixture_model.py | 8 +- surpyval/univariate/parametric/parametric.py | 39 +++++--- .../parametric/parametric_fitter.py | 24 ++--- .../accelerated_failure_time/aft_fitter.py | 8 +- .../accelerated_life/accelerated_life.py | 8 +- .../regression/additive_hazards/__init__.py | 8 +- .../additive_hazards_fitter.py | 9 +- .../univariate/regression/frailty/__init__.py | 11 ++- .../regression/parametric_regression_model.py | 91 +++++++++++-------- .../proportional_hazards/__init__.py | 8 +- .../proportional_hazards_fitter.py | 32 +++---- .../proportional_odds_fitter.py | 8 +- .../univariate/regression/regression_data.py | 19 +++- .../univariate/regression/tvc_schedule.py | 2 + surpyval/utils/__init__.py | 20 ++-- 30 files changed, 361 insertions(+), 180 deletions(-) diff --git a/.github/workflows/actions.yml b/.github/workflows/actions.yml index 65be7595..5f2bcd6e 100644 --- a/.github/workflows/actions.yml +++ b/.github/workflows/actions.yml @@ -82,6 +82,23 @@ jobs: python -m pytest -n auto --cov=surpyval --cov-report= --ignore=surpyval/tests/alpha --run-ml + # Execute the ``>>>`` examples in the docstrings and compare their + # printed output. This is what a user (or an agent) sees from + # ``help(Weibull.fit)``, so it is a documented promise like any + # other, and it went stale silently until it was checked. Kept as + # its own step -- separate from the suite above -- so a failure + # reads as "the docs drifted", not "a test broke", and so the + # doctest collection does not disturb coverage or xdist. + # + # ``surpyval/tests`` is excluded because the test modules have no + # user-facing examples, and ``surpyval/alpha`` because it is not + # part of the release contract (same reason the suite skips it). + # Option flags are set in pyproject.toml. + - name: doctests + run: + python -m pytest --doctest-modules surpyval + --ignore=surpyval/tests --ignore=surpyval/alpha + - name: coverage run: | coverage report diff --git a/docs/changelog.rst b/docs/changelog.rst index 0794a3c1..280e01cd 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -4,6 +4,66 @@ Changelog v0.19.1 (unreleased) -------------------- +- **CI now runs the docstring examples.** ``pytest --doctest-modules`` + over the package is a new step in the deployment workflow, and every + one of the 229 docstring examples passes. It was 59 failing tests + when the flag was first turned on. + + A docstring example is a promise about what the library prints, and + it is the one users and coding agents reach for first -- + ``help(Weibull.fit)`` is faster than opening the docs. Nothing was + checking it, so it drifted: examples recorded the output of an + optimiser two rewrites ago, of numpy 1's scalar repr, of a module + that has since moved. + + What the run found, beyond the cosmetic drift: + + - Twelve examples could not run at all. Six regression docstrings + (``PH``, ``AH``, ``PO``, ``AFT``, ``AcceleratedLife``, ``Frailty``) + were sketches -- ``model = PH(Weibull).fit(x, Z=covariates, c=c)`` + with ``x``, ``covariates`` and ``c`` never defined. Four more used + ``>>>`` on the continuation lines of a multi-line call, so pasting + them raised ``SyntaxError``. ``plotting_positions`` imported from + ``surpyval.nonparametric``, which moved to + ``surpyval.univariate.nonparametric`` several releases ago. + ``ParametricFitter.fit`` demonstrated ``how='MPP'`` on + interval-censored input, which now (correctly) requires the Turnbull + heuristic and raises without it. All are now runnable, with data. + + - The five ``ParametricRegressionModel`` prediction examples + (``sf``, ``ff``, ``df``, ``hf``, ``Hf``) had been copied from the + univariate ``Parametric`` class and never adapted: they built a + ``Weibull.from_params([10, 3])`` and called it with no covariates at + all, documenting a signature the method does not have. They now fit + a ``WeibullPH`` and pass ``Z``. + + - ``Parametric.var()`` claimed 11.229 for a Weibull(10, 3). The + variance is 10.533 (``100 Gamma(5/3) - (10 Gamma(4/3))^2``); the + code was right. + + - Several examples fitted unseeded random data and then recorded + specific digits, which cannot be reproducible. They now seed. + + ``Parametric.hf`` and ``Parametric.Hf`` returned a 0-d array + (``array(0.012)``) for scalar input where ``sf``, ``ff``, ``df`` and + ``qf`` all returned a numpy scalar, and ``cs`` did the same; their + own ``Returns`` sections promised "the scalar value ... if a scalar + was passed". That is now true. The 0-d array came from ``np.where``, + which does not collapse. + + Two doctest option flags are set in ``pyproject.toml``. + ``NORMALIZE_WHITESPACE``, because numpy picks its own line breaks and + column padding for an array and both move with the widest element -- + without it an example is only correct at the exact wrapping it was + captured at. ``ELLIPSIS``, so an example that ends in a fit can write + ``529.05371...`` rather than all seventeen digits: the trailing digits + of an optimiser's output are not part of what the docstring is + promising, and they move with the BLAS and the platform. Array reprs + are left exact -- numpy already prints only eight significant digits + there. + + This closes #158. + - **The distribution docstring examples now show what you actually see.** ``pytest --doctest-modules`` over ``distributions/`` is green: 139 examples, no failures. Previously 39 failed. @@ -25,8 +85,7 @@ v0.19.1 (unreleased) Two scalar examples had drifted in the last digit, and are re-recorded. With the examples now true, ``--doctest-modules`` is worth running in - CI -- that decision is separate, and is what would stop this - recurring. + CI, which is what stops this recurring; it is turned on above. - **``Gamma`` no longer offers probability plotting as a fit method.** ``Gamma.fit(x, how="MPP")`` now raises, joining ``Beta`` and @@ -83,11 +142,7 @@ v0.19.1 (unreleased) In every case the code was right and the documentation was wrong, which is the reassuring direction, but a reader checking their understanding against these would have been misled. They accumulated precisely - because the doctests are not run. - - Whether to run doctests in CI, and what to do about the 30 cosmetic - repr differences, is tracked separately -- it is a style decision, not - a correctness one. + because the doctests were not run, which is addressed above. - **Documented why a Turnbull fit does not equal a Kaplan-Meier fit.** ``Turnbull.fit`` defaults to ``turnbull_estimator="Fleming-Harrington"`` diff --git a/pyproject.toml b/pyproject.toml index 98f6d8ae..bc4849d8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -49,6 +49,25 @@ where = ["."] [tool.setuptools.package-data] "*" = ["datasets/*.csv", "py.typed"] +# pytest +[tool.pytest.ini_options] +# Applied to the ``--doctest-modules`` run (see .github/workflows/actions.yml), +# which executes the ``>>>`` examples in the docstrings and compares their +# printed output exactly. +# +# NORMALIZE_WHITESPACE: numpy chooses its own line breaks and column +# padding when it reprs an array, and both change with the width of the +# widest element. Without this an example is only correct for the exact +# terminal wrapping it was captured at, and a docstring wrapped to 79 +# columns by hand -- as most of these are -- can never match. +# +# ELLIPSIS: lets an example that ends in a fit write ``529.05...`` +# instead of all seventeen digits. The last few digits of an optimiser's +# output are not part of the promise the docstring is making, and they +# move with the BLAS, the platform and any change to the fitting +# ladder. +doctest_optionflags = "NORMALIZE_WHITESPACE ELLIPSIS" + # For Black pre-commit hook [tool.black] line-length = 79 diff --git a/surpyval/degradation/degradation_analysis.py b/surpyval/degradation/degradation_analysis.py index 2652d4bf..f6370e29 100644 --- a/surpyval/degradation/degradation_analysis.py +++ b/surpyval/degradation/degradation_analysis.py @@ -1193,8 +1193,8 @@ class DegradationAnalysis_: Censored Units : 0 Life Distribution : Weibull Parameters : - alpha: 441.4780882117898 - beta: 6.987078993008337 + alpha: 441.47809... + beta: 6.9870788... >>> model.pseudo_failure_times array([451.61290323, 500. , 318.18181818, 378.37837838]) """ diff --git a/surpyval/recurrent/parametric/cox_lewis.py b/surpyval/recurrent/parametric/cox_lewis.py index f1b3270a..140a7ad2 100644 --- a/surpyval/recurrent/parametric/cox_lewis.py +++ b/surpyval/recurrent/parametric/cox_lewis.py @@ -27,19 +27,19 @@ class CoxLewis(NHPPFitter): Process : Cox-Lewis Fitted by : MLE Parameters : - alpha: 0.3848528712360503 - beta: 0.1939447728437042 + alpha: 0.38481273... + beta: 0.19396672... >>> model.cif([1, 2, 3, 4, 5, 6]) - array([ 1.6215655 , 3.59019342, 5.98016527, 8.88166096, 12.40416155, - 16.68058024]) + array([ 1.62151879, 3.59013322, 5.98014113, 8.88174429, 12.40445268, + 16.6812175 ]) >>> >>> model.iif([1, 2, 3, 4, 5, 6]) - array([1.78389227, 2.16569736, 2.62921991, 3.19194983, 3.87512041, - 4.70450946]) + array([1.78385983, 2.16570551, 2.62928751, 3.19210196, 3.87539016, + 4.70494021]) >>> >>> model.inv_cif([1, 2, 3, 4, 5, 6]) - array([0.63923607, 1.20789182, 1.72001505, 2.18583659, 2.61303902, - 3.00753845]) + array([0.63925589, 1.20792021, 1.72004461, 2.18586234, 2.61305756, + 3.00754742]) """ def __init__(self): diff --git a/surpyval/recurrent/parametric/crow_amsaa.py b/surpyval/recurrent/parametric/crow_amsaa.py index ef66cf05..c427c654 100644 --- a/surpyval/recurrent/parametric/crow_amsaa.py +++ b/surpyval/recurrent/parametric/crow_amsaa.py @@ -27,19 +27,19 @@ class CrowAMSAA(NHPPFitter): Process : Crow-AMSAA Fitted by : MLE Parameters : - alpha: 3129.2801848331596 - beta: 1.239258986741094 + alpha: 913.84662... + beta: 1.4781707... >>> model.cif([1, 2, 3, 4, 5, 6]) - array([4.65842366e-05, 1.09974784e-04, 1.81767316e-04, 2.59625443e-04, - 3.42329130e-04, 4.29111278e-04]) + array([4.20072057e-05, 1.17030084e-04, 2.13103439e-04, 3.26040266e-04, + 4.53440995e-04, 5.93696079e-04]) >>> >>> model.iif([1, 2, 3, 4, 5, 6]) - array([5.77299348e-05, 6.81436206e-05, 7.50855947e-05, 8.04357921e-05, - 8.48468915e-05, 8.86300027e-05]) + array([6.20938211e-05, 8.64952211e-05, 1.05001087e-04, 1.20485793e-04, + 1.34052640e-04, 1.46264026e-04]) >>> >>> model.inv_cif([1, 2, 3, 4, 5, 6]) - array([ 3129.2801836 , 5474.64210552, 7593.6351178 , 9577.82762328, - 11467.45338053, 13284.98316248]) + array([ 913.84662107, 1460.57434899, 1921.54912724, 2334.39329941, + 2714.78099355, 3071.15581638]) """ def __init__(self): diff --git a/surpyval/recurrent/parametric/duane.py b/surpyval/recurrent/parametric/duane.py index 71f7963f..23bb1f2e 100644 --- a/surpyval/recurrent/parametric/duane.py +++ b/surpyval/recurrent/parametric/duane.py @@ -27,19 +27,19 @@ class Duane(NHPPFitter): Process : Duane Fitted by : MLE Parameters : - alpha: 1.2392945732132952 - b: 4.6568641229556424e-05 + alpha: 1.4782020... + b: 4.199455086392048e-05 >>> model.cif([1, 2, 3, 4, 5, 6]) - array([4.65686412e-05, 1.09940677e-04, 1.81713565e-04, 2.59551323e-04, - 3.42234115e-04, 4.28994958e-04]) + array([4.19945509e-05, 1.16997373e-04, 2.13046585e-04, 3.25956224e-04, + 4.53327287e-04, 5.93550595e-04]) >>> >>> model.iif([1, 2, 3, 4, 5, 6]) - array([5.77122644e-05, 6.81244421e-05, 7.50655449e-05, 8.04151364e-05, - 8.48257764e-05, 8.86085206e-05]) + array([6.20764328e-05, 8.64728804e-05, 1.04975302e-04, 1.20457293e-04, + 1.34021869e-04, 1.46231288e-04]) >>> >>> model.inv_cif([1, 2, 3, 4, 5, 6]) - array([ 3129.4028404 , 5474.76881015, 7593.73955973, 9577.89554535, - 11467.47544342, 13284.95262974]) + array([ 913.90063856, 1460.64614435, 1921.63239303, 2334.48481047, + 2714.87871658, 3071.25832646]) """ def __init__(self): diff --git a/surpyval/recurrent/parametric/hpp.py b/surpyval/recurrent/parametric/hpp.py index 4bbc1ea9..e5c7d78d 100644 --- a/surpyval/recurrent/parametric/hpp.py +++ b/surpyval/recurrent/parametric/hpp.py @@ -34,18 +34,17 @@ class HPP(CountingProcess): Process : Homogeneous Poisson Process Fitted by : MLE Parameters : - lambda: 0.000498450145693719 + lambda: 0.0023047023... >>> model.cif([1, 2, 3, 4, 5, 6]) - array([0.00049845, 0.0009969 , 0.00149535, 0.0019938 , 0.00249225, - 0.0029907 ]) + array([0.0023047 , 0.0046094 , 0.00691411, 0.00921881, 0.01152351, + 0.01382821]) >>> >>> model.iif([1, 2, 3, 4, 5, 6]) - array([0.00049845, 0.00049845, 0.00049845, 0.00049845, 0.00049845, - 0.00049845]) + array([0.0023047, 0.0023047, 0.0023047, 0.0023047, 0.0023047, 0.0023047]) >>> >>> model.inv_cif([1, 2, 3, 4, 5, 6]) - array([ 2006.21869336, 4012.43738672, 6018.65608009, 8024.87477345, - 10031.09346681, 12037.31216017]) + array([ 433.89551258, 867.79102516, 1301.68653774, 1735.58205032, + 2169.4775629 , 2603.37307548]) """ def __init__(self): diff --git a/surpyval/recurrent/regression/proportional_intensity.py b/surpyval/recurrent/regression/proportional_intensity.py index 3844f2fa..4a4d4035 100644 --- a/surpyval/recurrent/regression/proportional_intensity.py +++ b/surpyval/recurrent/regression/proportional_intensity.py @@ -37,10 +37,10 @@ class ProportionalIntensityModel( >>> c = data['arrest'].values >>> Z = data[["fin", "age", "race", "wexp", "mar", "paro", "prio"]].values >>> model = ProportionalIntensityNHPP.fit(x, Z, c, dist=CrowAMSAA) - >>> type(model) - surpyval.recurrent.regression.proportional_intensity.ProportionalIntensityModel + >>> type(model).__name__ + 'ProportionalIntensityModel' >>> model.cif([1, 2, 3], Z.mean(axis=0)) - array([0.00625402, 0.04304137, 0.13302238]) + array([8.84210972e-07, 2.79074784e-05, 2.10220821e-04]) """ def __repr__(self): diff --git a/surpyval/recurrent/renewal/generalized_one_renewal.py b/surpyval/recurrent/renewal/generalized_one_renewal.py index b54e148d..60b8b303 100644 --- a/surpyval/recurrent/renewal/generalized_one_renewal.py +++ b/surpyval/recurrent/renewal/generalized_one_renewal.py @@ -60,16 +60,16 @@ class GeneralizedOneRenewal(RenewalFitMixin): ========================= Distribution : Weibull Fitted by : MLE - Restoration Factor : -0.1730179893443181 + Restoration Factor : -0.17301846... Parameters : - alpha: 1.3919016662855024 - beta: 5.008872636271443 + alpha: 1.3919045... + beta: 5.0088611... >>> >>> np.random.seed(0) >>> np_model = model.count_terminated_simulation(len(x), 5000) >>> np_model.mcf(np.array([1, 2, 3, 4, 5, 6])) - array([0.1696 , 1.181 , 2.287 , 3.6696 , 5.58237921, - 8.54474127]) + array([0.1696 , 1.181 , 2.287 , 3.6694 , 5.58237925, + 8.54474531]) """ @staticmethod @@ -208,10 +208,10 @@ def fit_from_recurrent_data(self, data, dist=Weibull, init=None): ========================= Distribution : Weibull Fitted by : MLE - Restoration Factor : 0.4270960618530103 + Restoration Factor : 0.34027890... Parameters : - alpha: 1.3494830373118245 - beta: 2.7838386997223212 + alpha: 1.4115217... + beta: 3.5499343... """ self._check_dist_eligible(dist) validate_renewal_censoring(data.c, type(self).__name__) @@ -290,10 +290,10 @@ def fit(self, x, i=None, c=None, n=None, dist=Weibull, init=None): ========================= Distribution : Weibull Fitted by : MLE - Restoration Factor : 0.4270960618530103 + Restoration Factor : 0.34027890... Parameters : - alpha: 1.3494830373118245 - beta: 2.7838386997223212 + alpha: 1.4115217... + beta: 3.5499343... """ data = handle_xicn(x, i, c, n) return self.fit_from_recurrent_data(data, dist=dist, init=init) @@ -325,10 +325,10 @@ def fit_from_parameters(self, params, q, dist=Weibull): >>> from surpyval.recurrent import GeneralizedOneRenewal >>> >>> model = GeneralizedOneRenewal.fit_from_parameters( - [10, 2], - 0.2, - dist=Weibull - ) + ... [10, 2], + ... 0.2, + ... dist=Weibull + ... ) >>> model G1 Renewal SurPyval Model ========================= diff --git a/surpyval/recurrent/renewal/generalized_renewal.py b/surpyval/recurrent/renewal/generalized_renewal.py index 71210fa7..a72dadf6 100644 --- a/surpyval/recurrent/renewal/generalized_renewal.py +++ b/surpyval/recurrent/renewal/generalized_renewal.py @@ -60,10 +60,10 @@ class GeneralizedRenewal(RenewalFitMixin): Distribution : Weibull Fitted by : MLE Kijima Type : i - Restoration Factor : 0.1573211400037486 + Restoration Factor : 0.15732122... Parameters : - alpha: 1.261338468404201 - beta: 8.93900788677076 + alpha: 1.2613379... + beta: 8.9390232... >>> >>> np.random.seed(0) >>> np_model = model.count_terminated_simulation(len(x), 5000) @@ -259,10 +259,10 @@ def fit_from_recurrent_data( Distribution : Weibull Fitted by : MLE Kijima Type : i - Restoration Factor : 1.594694243423234e-11 + Restoration Factor : 1.3316262291443964e-16 Parameters : - alpha: 2.399029078569064 - beta: 2.753920439616154 + alpha: 2.3990296... + beta: 2.7539200... """ validate_renewal_censoring(data.c, type(self).__name__) reject_left_truncation(data, type(self).__name__) @@ -344,10 +344,10 @@ def fit( Distribution : Weibull Fitted by : MLE Kijima Type : i - Restoration Factor : 1.594694243423234e-11 + Restoration Factor : 1.3316262291443964e-16 Parameters : - alpha: 2.399029078569064 - beta: 2.753920439616154 + alpha: 2.3990296... + beta: 2.7539200... """ data = handle_xicn(x, i, c, n) return self.fit_from_recurrent_data(data, dist, kijima, init=init) @@ -381,10 +381,10 @@ def fit_from_parameters(self, params, q, kijima="i", dist=Weibull): >>> from surpyval.recurrent import GeneralizedRenewal >>> >>> model = GeneralizedRenewal.fit_from_parameters( - [10, 2], - 0.2, - dist=Normal - ) + ... [10, 2], + ... 0.2, + ... dist=Normal + ... ) >>> model Generalized Renewal SurPyval Model ================================== diff --git a/surpyval/univariate/information_criteria.py b/surpyval/univariate/information_criteria.py index 7db588f8..f9bd852e 100644 --- a/surpyval/univariate/information_criteria.py +++ b/surpyval/univariate/information_criteria.py @@ -59,7 +59,7 @@ def neg_ll(self) -> float: >>> x = Weibull.random(100, 10, 3) >>> model = Weibull.fit(x) >>> model.neg_ll() - 262.52685642385734 + 262.52685... """ if getattr(self, "data", None) is None: raise ValueError("Must have been fit with data") @@ -88,7 +88,7 @@ def bic(self) -> float: >>> x = Weibull.random(100, 10, 3) >>> model = Weibull.fit(x) >>> model.bic() - 534.2640532196908 + np.float64(534.26405...) References ---------- @@ -124,7 +124,7 @@ def aic(self) -> float: >>> x = Weibull.random(100, 10, 3) >>> model = Weibull.fit(x) >>> model.aic() - 529.0537128477147 + 529.05371... """ if hasattr(self, "_aic"): return self._aic @@ -152,7 +152,7 @@ def aic_c(self) -> float: >>> x = Weibull.random(100, 10, 3) >>> model = Weibull.fit(x) >>> model.aic_c() - 529.1774241879209 + np.float64(529.17742...) """ if hasattr(self, "_aic_c"): return self._aic_c diff --git a/surpyval/univariate/nonparametric/nonparametric.py b/surpyval/univariate/nonparametric/nonparametric.py index b9f41ad1..fadbeed5 100755 --- a/surpyval/univariate/nonparametric/nonparametric.py +++ b/surpyval/univariate/nonparametric/nonparametric.py @@ -686,7 +686,7 @@ def mean(self, tau: float | None = None) -> float: >>> x = np.array([1, 2, 3, 4, 5]) >>> model = KaplanMeier.fit(x) >>> model.mean() - 3.0 + 3.0000000000000004 """ if np.min(self.x) < 0: raise ValueError( diff --git a/surpyval/univariate/nonparametric/plotting_positions.py b/surpyval/univariate/nonparametric/plotting_positions.py index 254ac29f..976ea04e 100755 --- a/surpyval/univariate/nonparametric/plotting_positions.py +++ b/surpyval/univariate/nonparametric/plotting_positions.py @@ -97,7 +97,7 @@ def plotting_positions( Examples -------- - >>> from surpyval.nonparametric import plotting_positions + >>> from surpyval.univariate.nonparametric import plotting_positions >>> import numpy as np >>> x = np.array([1, 2, 3, 4, 5, 6, 7, 8]) >>> x, r, d, F = plotting_positions(x, heuristic="Filliben") diff --git a/surpyval/univariate/nonparametric/success_run.py b/surpyval/univariate/nonparametric/success_run.py index a115b70c..5d27b6f3 100755 --- a/surpyval/univariate/nonparametric/success_run.py +++ b/surpyval/univariate/nonparametric/success_run.py @@ -32,7 +32,7 @@ def success_run(n, confidence=None, alpha=None): >>> from surpyval import success_run >>> success_run(10) - 0.7411344491069477 + np.float64(0.7411344491069477) """ if confidence and alpha: raise ValueError("Only one of confidence or alpha can be specified") diff --git a/surpyval/univariate/parametric/mixture_model.py b/surpyval/univariate/parametric/mixture_model.py index 6615d657..4d36919d 100755 --- a/surpyval/univariate/parametric/mixture_model.py +++ b/surpyval/univariate/parametric/mixture_model.py @@ -297,11 +297,11 @@ def fit( Sub-Distributions : 2 Fitted by : EM Weights : - 0.6094710980384728, - 0.39052890196152723 + 0.61848918..., + 0.38151081... Parameters : - alpha: [ 5.8855232 17.23187124] - beta: [ 2.04051304 11.01565277] + alpha: [ 6.32508961 17.37701969] + beta: [ 1.83105154 12.01392721] """ data = SurpyvalData(x=x, c=c, n=n, t=t, tl=tl, tr=tr, xl=xl, xr=xr) diff --git a/surpyval/univariate/parametric/parametric.py b/surpyval/univariate/parametric/parametric.py index 9e7a25c6..a0b1dcb7 100755 --- a/surpyval/univariate/parametric/parametric.py +++ b/surpyval/univariate/parametric/parametric.py @@ -549,7 +549,7 @@ def sf(self, x: npt.ArrayLike) -> npt.NDArray: >>> from surpyval import Weibull >>> model = Weibull.from_params([10, 3]) >>> model.sf(2) - 0.9920319148370607 + np.float64(0.9920319148370607) >>> model.sf([1, 2, 3, 4, 5]) array([0.9990005 , 0.99203191, 0.97336124, 0.938005 , 0.8824969 ]) """ @@ -591,7 +591,7 @@ def ff(self, x: npt.ArrayLike) -> npt.NDArray: >>> from surpyval import Weibull >>> model = Weibull.from_params([10, 3]) >>> model.ff(2) - 0.007968085162939342 + np.float64(0.007968085162939372) >>> model.ff([1, 2, 3, 4, 5]) array([0.0009995 , 0.00796809, 0.02663876, 0.061995 , 0.1175031 ]) """ @@ -632,7 +632,7 @@ def df(self, x: npt.ArrayLike) -> npt.NDArray: >>> from surpyval import Weibull >>> model = Weibull.from_params([10, 3]) >>> model.df(2) - 0.01190438297804473 + np.float64(0.01190438297804473) >>> model.df([1, 2, 3, 4, 5]) array([0.002997 , 0.01190438, 0.02628075, 0.04502424, 0.06618727]) """ @@ -679,7 +679,7 @@ def hf(self, x: npt.ArrayLike) -> npt.NDArray: >>> from surpyval import Weibull >>> model = Weibull.from_params([10, 3]) >>> model.hf(2) - 0.012000000000000002 + np.float64(0.012000000000000002) >>> model.hf([1, 2, 3, 4, 5]) array([0.003, 0.012, 0.027, 0.048, 0.075]) """ @@ -687,7 +687,12 @@ def hf(self, x: npt.ArrayLike) -> npt.NDArray: if (self.p == 1) and (self.f0 == 0): xg = x - self.gamma # type: ignore[operator] s0 = getattr(self.dist, "support", (-np.inf, np.inf))[0] - return np.where(xg < s0, 0.0, self.dist.hf(xg, *self.params)) + out = np.where(xg < s0, 0.0, self.dist.hf(xg, *self.params)) + # ``np.where`` hands back a 0-d array for scalar input, so a + # scalar argument used to come out as ``array(0.012)`` while + # every sibling method returned a numpy scalar. ``[()]`` is a + # no-op on a real array and unwraps the 0-d case. + return out[()] else: return self.df(x) / self.sf(x) @@ -719,7 +724,7 @@ def Hf(self, x: npt.ArrayLike) -> npt.NDArray: >>> from surpyval import Weibull >>> model = Weibull.from_params([10, 3]) >>> model.Hf(2) - 0.008000000000000002 + np.float64(0.008000000000000002) >>> model.Hf([1, 2, 3, 4, 5]) array([0.001, 0.008, 0.027, 0.064, 0.125]) """ @@ -728,7 +733,8 @@ def Hf(self, x: npt.ArrayLike) -> npt.NDArray: if (self.p == 1) and (self.f0 == 0): xg = x - self.gamma # type: ignore[operator] s0 = getattr(self.dist, "support", (-np.inf, np.inf))[0] - return np.where(xg < s0, 0.0, self.dist.Hf(xg, *self.params)) + out = np.where(xg < s0, 0.0, self.dist.Hf(xg, *self.params)) + return out[()] else: return -np.log(self.sf(x)) @@ -757,7 +763,7 @@ def qf(self, p: npt.ArrayLike) -> npt.NDArray: >>> from surpyval import Weibull >>> model = Weibull.from_params([10, 3]) >>> model.qf(0.2) - 6.06542793124108 + np.float64(6.06542793124108) >>> model.qf([.1, .2, .3, .4, .5]) array([4.72308719, 6.06542793, 7.09181722, 7.99387877, 8.84997045]) @@ -823,14 +829,14 @@ def cs(self, x: npt.ArrayLike, X: npt.ArrayLike) -> npt.NDArray: >>> from surpyval import Weibull >>> model = Weibull.from_params([10, 3]) >>> model.cs(11, 10) - 0.00025840046151723767 + np.float64(0.00025840046151723767) """ x = np.asarray(x) X = np.asarray(X) Xg = X - self.gamma # type: ignore[operator] cs = np.array(self.dist.cs(x, Xg, *self.params)) cs[cs > 1.0] = 1 - return cs + return cs[()] def random( self, @@ -972,7 +978,7 @@ def mean(self) -> float: >>> from surpyval import Weibull >>> model = Weibull.from_params([10, 3]) >>> model.mean() - 8.929795115692489 + np.float64(8.929795115692489) """ if not hasattr(self, "_mean"): # Defective mean: the zero-inflated mass f0 sits at 0 and @@ -998,7 +1004,7 @@ def var(self) -> float: >>> from surpyval import Weibull >>> model = Weibull.from_params([10, 3]) >>> model.var() - 11.229... + np.float64(10.533288486847923) """ m1 = self.dist._moment(1, *self.params) m2 = self.dist._moment(2, *self.params) @@ -1071,7 +1077,7 @@ def entropy(self) -> float: >>> from surpyval import Normal >>> model = Normal.from_params([10, 3]) >>> model.entropy() - 2.5175508218727822 + np.float64(2.5175508218727822) Notes ----- @@ -1556,16 +1562,19 @@ def plot( Returns ------- - plot : list - list of a matplotlib plot object + plot : matplotlib.axes.Axes + the axes the probability plot was drawn onto Examples -------- + >>> import numpy as np >>> from surpyval import Weibull + >>> np.random.seed(1) >>> x = Weibull.random(100, 10, 3) >>> model = Weibull.fit(x) >>> model.plot() + """ if ax is None: ax = plt.gcf().gca() diff --git a/surpyval/univariate/parametric/parametric_fitter.py b/surpyval/univariate/parametric/parametric_fitter.py index 43d2842f..f01cde16 100755 --- a/surpyval/univariate/parametric/parametric_fitter.py +++ b/surpyval/univariate/parametric/parametric_fitter.py @@ -753,6 +753,7 @@ def fit( -------- >>> from surpyval import Weibull >>> import numpy as np + >>> np.random.seed(1) >>> x = Weibull.random(100, 10, 4) >>> model = Weibull.fit(x) >>> print(model) @@ -761,8 +762,8 @@ def fit( Distribution : Weibull Fitted by : MLE Parameters : - alpha: 10.551521182640098 - beta: 3.792549834495306 + alpha: 9.8150187... + beta: 3.7987404... >>> Weibull.fit(x, how='MPS', fixed={'alpha' : 10}) Parametric SurPyval Model ========================= @@ -770,15 +771,16 @@ def fit( Fitted by : MPS Parameters : alpha: 10.0 - beta: 3.4314657446866836 - >>> Weibull.fit(xl=x-1, xr=x+1, how='MPP') + beta: 3.6707965... + >>> Weibull.fit(xl=np.floor(x), xr=np.ceil(x), how='MPP', + ... heuristic='Turnbull') Parametric SurPyval Model ========================= Distribution : Weibull Fitted by : MPP Parameters : - alpha: 9.943092756713078 - beta: 8.613016934518258 + alpha: 9.9501683... + beta: 3.2119714... >>> c = np.zeros_like(x) >>> c[x > 13] = 1 >>> x[x > 13] = 13 @@ -790,8 +792,8 @@ def fit( Distribution : Weibull Fitted by : MLE Parameters : - alpha: 10.363725328793413 - beta: 4.9886821457305865 + alpha: 9.8935844... + beta: 3.7868860... """ surv_data = SurpyvalData( @@ -892,10 +894,10 @@ def fit_from_df( ========================= Distribution : Weibull Fitted by : MLE - Offset (gamma) : 39.76562962867477 + Offset (gamma) : 39.765577... Parameters : - alpha: 7.141925216146524 - beta: 2.6204524040137844 + alpha: 7.1419836... + beta: 2.6204759... """ if not isinstance(df, pd.DataFrame): diff --git a/surpyval/univariate/regression/accelerated_failure_time/aft_fitter.py b/surpyval/univariate/regression/accelerated_failure_time/aft_fitter.py index 13d594e3..ae87b354 100644 --- a/surpyval/univariate/regression/accelerated_failure_time/aft_fitter.py +++ b/surpyval/univariate/regression/accelerated_failure_time/aft_fitter.py @@ -144,8 +144,14 @@ def AFT(distribution): Examples -------- + >>> import numpy as np >>> from surpyval import Weibull >>> from surpyval import AFT - >>> model = AFT(Weibull).fit(x, Z=covariates, c=c) + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> model = AFT(Weibull).fit(x, Z=Z) + >>> model.params.round(3) + array([9.629, 1.751, 0.473]) """ return AFTFitter(distribution) diff --git a/surpyval/univariate/regression/accelerated_life/accelerated_life.py b/surpyval/univariate/regression/accelerated_life/accelerated_life.py index 3fa5fc39..4626f273 100644 --- a/surpyval/univariate/regression/accelerated_life/accelerated_life.py +++ b/surpyval/univariate/regression/accelerated_life/accelerated_life.py @@ -38,9 +38,15 @@ def AcceleratedLife(distribution, life_model): Examples -------- + >>> import numpy as np >>> from surpyval import Weibull >>> from surpyval import AcceleratedLife, Power - >>> model = AcceleratedLife(Weibull, Power).fit(x, c=c, Z=stress) + >>> np.random.seed(1) + >>> stress = np.repeat([20.0, 30.0, 40.0], 40).reshape(-1, 1) + >>> x = Weibull.random(120, 10, 3) * (100.0 / stress[:, 0]) + >>> model = AcceleratedLife(Weibull, Power).fit(x, Z=stress) + >>> model.params.round(3) + array([ 1. , 2.831, 558.686, -0.828]) """ if distribution.name not in _LIFE_PARAM_MAP: supported = list(_LIFE_PARAM_MAP.keys()) diff --git a/surpyval/univariate/regression/additive_hazards/__init__.py b/surpyval/univariate/regression/additive_hazards/__init__.py index da90dd75..9a840da1 100644 --- a/surpyval/univariate/regression/additive_hazards/__init__.py +++ b/surpyval/univariate/regression/additive_hazards/__init__.py @@ -32,8 +32,14 @@ def AH(distribution): Examples -------- + >>> import numpy as np >>> from surpyval import Weibull, AH - >>> model = AH(Weibull).fit(x, Z=covariates, c=c) + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> model = AH(Weibull).fit(x, Z=Z) + >>> model.params.round(3) + array([9.332, 1.851, 0.086]) """ return AdditiveHazardsFitter.create(distribution) diff --git a/surpyval/univariate/regression/additive_hazards/additive_hazards_fitter.py b/surpyval/univariate/regression/additive_hazards/additive_hazards_fitter.py index 59c7c13e..b2115f4a 100644 --- a/surpyval/univariate/regression/additive_hazards/additive_hazards_fitter.py +++ b/surpyval/univariate/regression/additive_hazards/additive_hazards_fitter.py @@ -200,8 +200,15 @@ def fit( Examples -------- - >>> from surpyval import WeibullAH + >>> import numpy as np + >>> from surpyval import Weibull, WeibullAH + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> c = np.zeros(100) >>> model = WeibullAH.fit(x=x, Z=Z, c=c) + >>> model.params.round(3) + array([9.332, 1.851, 0.086]) """ data, prep = prepare_regression_fit( self, diff --git a/surpyval/univariate/regression/frailty/__init__.py b/surpyval/univariate/regression/frailty/__init__.py index c09dfeaa..724df8bd 100644 --- a/surpyval/univariate/regression/frailty/__init__.py +++ b/surpyval/univariate/regression/frailty/__init__.py @@ -34,8 +34,17 @@ def Frailty(distribution, family="gamma"): Examples -------- + >>> import numpy as np >>> from surpyval import Frailty, Weibull - >>> model = Frailty(Weibull).fit(x, Z=Z, c=c, groups=unit_id) + >>> np.random.seed(1) + >>> unit_id = np.repeat(np.arange(20), 5) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> model = Frailty(Weibull).fit(x, Z=Z, groups=unit_id) + >>> model.beta.round(3) + array([0.829]) + >>> model.n_groups + 20 """ return FrailtyFitter.create(distribution, family) diff --git a/surpyval/univariate/regression/parametric_regression_model.py b/surpyval/univariate/regression/parametric_regression_model.py index a36c8664..efa9d8e7 100644 --- a/surpyval/univariate/regression/parametric_regression_model.py +++ b/surpyval/univariate/regression/parametric_regression_model.py @@ -430,12 +430,14 @@ def sf( Examples -------- - >>> from surpyval import Weibull - >>> model = Weibull.from_params([10, 3]) - >>> model.sf(2) - 0.9920319148370607 - >>> model.sf([1, 2, 3, 4, 5]) - array([0.9990005 , 0.99203191, 0.97336124, 0.938005 , 0.8824969 ]) + >>> import numpy as np + >>> from surpyval import Weibull, WeibullPH + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> model = WeibullPH.fit(x, Z) + >>> model.sf([1, 2, 3], [[0], [0], [1]]).round(4) + array([0.9812, 0.9382, 0.7429]) """ if isinstance(x, list): x = np.array(x) @@ -669,12 +671,14 @@ def ff( Examples -------- - >>> from surpyval import Weibull - >>> model = Weibull.from_params([10, 3]) - >>> model.ff(2) - 0.007968085162939342 - >>> model.ff([1, 2, 3, 4, 5]) - array([0.0009995 , 0.00796809, 0.02663876, 0.061995 , 0.1175031 ]) + >>> import numpy as np + >>> from surpyval import Weibull, WeibullPH + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> model = WeibullPH.fit(x, Z) + >>> model.ff([1, 2, 3], [[0], [0], [1]]).round(4) + array([0.0188, 0.0618, 0.2571]) """ if isinstance(x, list): x = np.array(x) @@ -708,12 +712,14 @@ def df( Examples -------- - >>> from surpyval import Weibull - >>> model = Weibull.from_params([10, 3]) - >>> model.df(2) - 0.01190438297804473 - >>> model.df([1, 2, 3, 4, 5]) - array([0.002997 , 0.01190438, 0.02628075, 0.04502424, 0.06618727]) + >>> import numpy as np + >>> from surpyval import Weibull, WeibullPH + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> model = WeibullPH.fit(x, Z) + >>> model.df([1, 2, 3], [[0], [0], [1]]).round(4) + array([0.0326, 0.0524, 0.1289]) """ if isinstance(x, list): x = np.array(x) @@ -748,12 +754,14 @@ def hf( Examples -------- - >>> from surpyval import Weibull - >>> model = Weibull.from_params([10, 3]) - >>> model.hf(2) - 0.012000000000000002 - >>> model.hf([1, 2, 3, 4, 5]) - array([0.003, 0.012, 0.027, 0.048, 0.075]) + >>> import numpy as np + >>> from surpyval import Weibull, WeibullPH + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> model = WeibullPH.fit(x, Z) + >>> model.hf([1, 2, 3], [[0], [0], [1]]).round(4) + array([0.0332, 0.0559, 0.1735]) """ if isinstance(x, list): x = np.array(x) @@ -789,12 +797,14 @@ def Hf( Examples -------- - >>> from surpyval import Weibull - >>> model = Weibull.from_params([10, 3]) - >>> model.Hf(2) - 0.008000000000000002 - >>> model.Hf([1, 2, 3, 4, 5]) - array([0.001, 0.008, 0.027, 0.064, 0.125]) + >>> import numpy as np + >>> from surpyval import Weibull, WeibullPH + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> model = WeibullPH.fit(x, Z) + >>> model.Hf([1, 2, 3], [[0], [0], [1]]).round(4) + array([0.0189, 0.0638, 0.2972]) """ if isinstance(x, list): x = np.array(x) @@ -826,15 +836,22 @@ def random( Examples -------- - >>> from surpyval import Weibull - >>> model = Weibull.from_params([10, 3]) - >>> np.random.seed(1) - >>> model.random(1) - array([8.14127103]) - >>> from surpyval import WeibullPH + >>> import numpy as np + >>> from surpyval import Weibull, WeibullPH >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) >>> model = WeibullPH.fit(x, Z) - >>> x_rand, Z_rand = model.random(10, Z[:1]) + >>> np.random.seed(1) + >>> x_rand, Z_rand = model.random(5, Z[:1]) + >>> x_rand.round(3) + array([ 8.919, 5.095, 33.929, 10.666, 13.97 ]) + >>> Z_rand + array([[0.], + [0.], + [0.], + [0.], + [0.]]) """ # Dispatch to the regression fitter's own covariate-aware sampler # (#261): the previous implementation ignored ``Z`` entirely and diff --git a/surpyval/univariate/regression/proportional_hazards/__init__.py b/surpyval/univariate/regression/proportional_hazards/__init__.py index ae74e0fa..add16b8c 100644 --- a/surpyval/univariate/regression/proportional_hazards/__init__.py +++ b/surpyval/univariate/regression/proportional_hazards/__init__.py @@ -31,9 +31,15 @@ def PH(distribution): Examples -------- + >>> import numpy as np >>> from surpyval import Weibull >>> from surpyval import PH - >>> model = PH(Weibull).fit(x, Z=covariates, c=c) + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> model = PH(Weibull).fit(x, Z=Z) + >>> model.params.round(3) + array([9.629, 1.751, 0.829]) """ return ProportionalHazardsFitter.create(distribution) diff --git a/surpyval/univariate/regression/proportional_hazards/proportional_hazards_fitter.py b/surpyval/univariate/regression/proportional_hazards/proportional_hazards_fitter.py index f0f17f64..d62f1ac5 100644 --- a/surpyval/univariate/regression/proportional_hazards/proportional_hazards_fitter.py +++ b/surpyval/univariate/regression/proportional_hazards/proportional_hazards_fitter.py @@ -188,17 +188,13 @@ def fit( >>> from surpyval import WeibullPH >>> from surpyval.datasets import load_tires_data - >>> from autograd import numpy as anp - >>> import numpy as np - >>> >>> data = load_tires_data() - >>> >>> x = data['Survival'].values >>> c = data['Censoring'].values >>> Z = data[[ - 'Wedge gauge', 'Interbelt gauge', 'Peel force', - 'Wedge gauge×peel force' - ]].values + ... 'Wedge gauge', 'Interbelt gauge', 'Peel force', + ... 'Wedge gauge×peel force' + ... ]].values >>> model = WeibullPH.fit(x=x, Z=Z, c=c) >>> model Parametric Regression SurPyval Model @@ -208,13 +204,13 @@ def fit( Regression Model : Log Linear [e^(beta'Z)] Fitted by : MLE Distribution : - alpha: 0.24255054642143947 - beta: 16.057791674515805 + alpha: 0.24255136... + beta: 16.057785... Regression Model : - beta_0: -9.165062641226692 - beta_1: -7.998599877425742 - beta_2: -27.503283340963034 - beta_3: 18.38550143851751 + beta_0: -9.1650627... + beta_1: -7.9985730... + beta_2: -27.503185... + beta_3: 18.385445... >>> model = WeibullPH.fit(x=x, Z=Z, c=c, fixed={"beta": 15}) >>> model Parametric Regression SurPyval Model @@ -224,13 +220,13 @@ def fit( Regression Model : Log Linear [e^(beta'Z)] Fitted by : MLE Distribution : - alpha: 0.23772915681951018 + alpha: 0.23772966... beta: 15.0 Regression Model : - beta_0: -8.628333861229965 - beta_1: -7.617541980158942 - beta_2: -25.952407717383302 - beta_3: 17.270173771235655 + beta_0: -8.6283269... + beta_1: -7.6175293... + beta_2: -25.952367... + beta_3: 17.270148... """ data, prep = prepare_regression_fit( self, diff --git a/surpyval/univariate/regression/proportional_odds/proportional_odds_fitter.py b/surpyval/univariate/regression/proportional_odds/proportional_odds_fitter.py index e6074230..d0304bbd 100644 --- a/surpyval/univariate/regression/proportional_odds/proportional_odds_fitter.py +++ b/surpyval/univariate/regression/proportional_odds/proportional_odds_fitter.py @@ -200,8 +200,14 @@ def PO(distribution): Examples -------- + >>> import numpy as np >>> from surpyval import Logistic >>> from surpyval import PO - >>> model = PO(Logistic).fit(x, Z=covariates, c=c) + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Logistic.random(100, 10, 2) + 2.0 * Z[:, 0] + >>> model = PO(Logistic).fit(x, Z=Z) + >>> model.params.round(3) + array([9.708, 2.337, 0.918]) """ return ProportionalOddsFitter(distribution) diff --git a/surpyval/univariate/regression/regression_data.py b/surpyval/univariate/regression/regression_data.py index 259390ec..b85bf809 100644 --- a/surpyval/univariate/regression/regression_data.py +++ b/surpyval/univariate/regression/regression_data.py @@ -307,11 +307,26 @@ def fit_from_df( Examples -------- - >>> from surpyval import WeibullPH + >>> import numpy as np + >>> import pandas as pd + >>> from surpyval import Weibull, WeibullPH + >>> np.random.seed(1) + >>> age = np.random.uniform(20, 60, 100) + >>> weight = np.random.uniform(50, 100, 100) + >>> time = Weibull.random(100, 10, 2) * np.exp(-0.02 * (age - 40)) + >>> df = pd.DataFrame({ + ... "time": time, + ... "age": age, + ... "weight": weight, + ... "censored": np.zeros(100, dtype=int), + ... }) >>> model = WeibullPH.fit_from_df( ... df, x_col="time", Z_cols=["age", "weight"], c_col="censored" ... ) - >>> model.sf([10, 20], df[["age", "weight"]]) + >>> model.feature_names + ['age', 'weight'] + >>> model.sf([10, 20], df[["age", "weight"]].head(2)).round(4) + array([0.4757, 0.0024]) """ Z, feature_names, model_spec = design_matrix_from_df( df, Z_cols, formula diff --git a/surpyval/univariate/regression/tvc_schedule.py b/surpyval/univariate/regression/tvc_schedule.py index f2589a4c..25ae5d8b 100644 --- a/surpyval/univariate/regression/tvc_schedule.py +++ b/surpyval/univariate/regression/tvc_schedule.py @@ -448,7 +448,9 @@ def from_expression(cls, expr, horizon, resolution=1.0, t0=0.0): Examples -------- >>> StepSchedule.from_expression("0.9 if t % 24 < 8 else 0.3", 96) + StepSchedule(step, 9 segment(s), p=1) >>> StepSchedule.from_expression("0.3 * 2 ** floor(t / 1000)", 5000) + StepSchedule(step, 6 segment(s), p=1) """ exprs = [expr] if isinstance(expr, str) else list(expr) if len(exprs) == 0: diff --git a/surpyval/utils/__init__.py b/surpyval/utils/__init__.py index 5f6064f9..604ae09a 100755 --- a/surpyval/utils/__init__.py +++ b/surpyval/utils/__init__.py @@ -302,7 +302,7 @@ def xrd_handler(x, r, d): >>> r array([5, 4, 3, 2, 1]) >>> d - array([1, 1, 1, 1, 1])) + array([1, 1, 1, 1, 1]) """ try: @@ -824,7 +824,7 @@ def xcnt_to_xrd(x, c=None, n=None, t=None, **kwargs): >>> n = np.array([1, 1, 1, 1, 1]) >>> x, r, d = xcnt_to_xrd(x, c, n) >>> x - array([1, 2, 3, 4, 5]) + array([1., 2., 3., 4., 5.]) >>> r array([5, 4, 3, 2, 1]) >>> d @@ -917,8 +917,8 @@ def xrd_to_xcnt(x, r, d): >>> r = np.array([5, 4, 3, 2, 1]) >>> d = np.array([1, 0, 0, 1, 1]) >>> x, c, n, t = xrd_to_xcnt(x, r, d) - >>> x, c, n, t - array([1, 2, 3, 4, 5]) + >>> x + array([1., 2., 3., 4., 5.]) >>> c array([0, 1, 1, 0, 0]) >>> n @@ -928,7 +928,7 @@ def xrd_to_xcnt(x, r, d): [-inf, inf], [-inf, inf], [-inf, inf], - [-inf, inf]])) + [-inf, inf]]) """ n_f = np.copy(d) x_f = np.copy(x) @@ -1001,7 +1001,7 @@ def fsli_to_xcnt(f=None, s=None, l=None, i=None): >>> i = [] >>> x, c, n, t = fsli_to_xcnt(f, s, l, i) >>> x - array([1, 2, 3, 4, 5]) + array([1., 2., 3., 4., 5.]) >>> c array([0, 1, 1, 0, 0]) >>> n @@ -1476,8 +1476,12 @@ def fs_to_xrd(f, s): >>> f = [1, 4, 5] >>> s = [2, 3] >>> x, r, d = fs_to_xrd(f, s) - >>> x, r, d - (array([1, 2, 3, 4, 5]), array([5, 4, 3, 2, 1]), array([1, 0, 0, 1, 1])) + >>> x + array([1., 2., 3., 4., 5.]) + >>> r + array([5, 4, 3, 2, 1]) + >>> d + array([1, 0, 0, 1, 1]) """ x, c, n, _ = fs_to_xcnt(f, s) return xcnt_to_xrd(x, c, n) From 88553fa0809cd4c60e154a0d7376126ea18636e6 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 4 Aug 2026 21:25:08 +0000 Subject: [PATCH 07/13] Compare the numbers in doctest output as numbers CI failed the doctest step on 3.12 and 3.13 while passing on 3.11. The diagnosis was the assumption made in the previous commit: that numpy's eight-significant-digit array repr was tight enough to pin. It is not. The Duane example lands on b = 4.1995e-05 under 3.11 and 4.2032e-05 under 3.12 -- a third-significant-figure difference in an optimiser whose fit is poorly conditioned on ten events. Sixteen of the 229 examples disagree somewhere in their digits between those Pythons. Trimming each documented number back to the digits that agree everywhere would make the docstring show something the reader's own session will not produce, which is what these examples exist to avoid. So the examples record the real output, in full, and the comparison changes instead. conftest.py patches doctest.OutputChecker.check_output with a fallback that runs only after the ordinary text comparison has failed. It fires when the two outputs are identical apart from their numeric literals -- same words, same brackets, same integer-versus-float shape, so "1" never matches "1." and a dtype change is still a failure -- and then compares the numbers pairwise. rel_tol is 1e-3, set by the loosest genuine disagreement observed with no margin beyond it; abs_tol is 1e-12 for a restoration factor whose true value is zero and which surfaces as 1e-16 with whatever mantissa the optimiser stopped on. Patched on the base class rather than installed as a checker: pytest builds its own LiteralsOutputChecker subclass and calls up to this method, so the override survives both plain doctest and pytest without depending on pytest internals. What this forgives is a value drifting inside the tolerance. What it still catches is every defect the sweep found: a stale value from another parameterisation, the wrong function being called, the wrong shape, an exception, a missing import. test_doctest_checker.py pins both halves, using the real output pairs observed on different Pythons -- so the tolerance cannot be widened without a test saying why. ELLIPSIS is dropped from the option flags; it was only there to support the trimmed numbers, and leaving it on invites reintroducing them. Verified against real 3.11, 3.12 and 3.13 interpreters locally: 229 doctests pass on each, and the same 3.12 run shows 16 text-level differences with the fallback disabled. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts --- .github/workflows/actions.yml | 3 +- conftest.py | 92 ++++++++++++++ docs/changelog.rst | 43 +++++-- pyproject.toml | 21 ++-- surpyval/degradation/degradation_analysis.py | 4 +- surpyval/recurrent/parametric/cox_lewis.py | 4 +- surpyval/recurrent/parametric/crow_amsaa.py | 4 +- surpyval/recurrent/parametric/duane.py | 2 +- surpyval/recurrent/parametric/hpp.py | 2 +- .../renewal/generalized_one_renewal.py | 18 +-- .../recurrent/renewal/generalized_renewal.py | 14 +-- surpyval/tests/test_doctest_checker.py | 119 ++++++++++++++++++ surpyval/univariate/information_criteria.py | 8 +- .../univariate/parametric/mixture_model.py | 4 +- .../parametric/parametric_fitter.py | 20 +-- .../proportional_hazards_fitter.py | 24 ++-- 16 files changed, 307 insertions(+), 75 deletions(-) create mode 100644 surpyval/tests/test_doctest_checker.py diff --git a/.github/workflows/actions.yml b/.github/workflows/actions.yml index 5f2bcd6e..4fc6629c 100644 --- a/.github/workflows/actions.yml +++ b/.github/workflows/actions.yml @@ -93,7 +93,8 @@ jobs: # ``surpyval/tests`` is excluded because the test modules have no # user-facing examples, and ``surpyval/alpha`` because it is not # part of the release contract (same reason the suite skips it). - # Option flags are set in pyproject.toml. + # Option flags are set in pyproject.toml; conftest.py adds the + # numeric comparison that lets the examples record real output. - name: doctests run: python -m pytest --doctest-modules surpyval diff --git a/conftest.py b/conftest.py index 1df740b0..4e8babc3 100644 --- a/conftest.py +++ b/conftest.py @@ -31,8 +31,100 @@ the run dies on "unrecognized arguments". """ +import doctest +import math +import re + import pytest +# --------------------------------------------------------------------------- +# Numeric comparison for the ``--doctest-modules`` run +# --------------------------------------------------------------------------- +# doctest compares printed output as text. That is the wrong test for a +# library whose examples end in an optimiser: the same fit lands on +# ``b = 4.1995e-05`` under one Python and ``4.2032e-05`` under the next, +# and numpy prints eight significant digits either way, so a byte-exact +# comparison fails on a difference no reader would call a difference. +# +# The alternative -- trimming every documented number to the digits that +# happen to agree everywhere -- makes the docstring show something the +# user's own session will not produce, which is the thing these examples +# exist to avoid. So the examples record the real output, in full, and +# the numbers in it are compared as numbers. +# +# The fallback only runs after the ordinary text comparison has failed, +# and only fires when the two outputs are identical apart from their +# numeric literals -- same words, same brackets, same integer-vs-float +# shape ("1" never matches "1.", which is a dtype change worth +# failing on). What it forgives is the value drifting inside a +# tolerance. What it still catches is everything that actually went +# wrong when this was first switched on: a stale value from another +# parameterisation, a different function being called, the wrong array +# shape, an exception, a missing import. +# +# RTOL is set by the loosest genuine disagreement between supported +# Pythons -- the Duane example above, at 9e-4 -- with no margin beyond +# that. ATOL exists for the one other case, a restoration factor whose +# true value is zero and which surfaces as 1e-16 with whatever sign and +# mantissa the optimiser stopped on; relative tolerance is meaningless +# there. +RTOL = 1e-3 +ATOL = 1e-12 + +_NUMBER = re.compile(r"[-+]?(?:\d+\.\d*|\.\d+|\d+)(?:[eE][-+]?\d+)?") + + +def _skeleton(text: str) -> str: + """The text with each number replaced by its *kind*. + + Integers and floats get different placeholders so that a change in + dtype -- ``array([1, 2])`` becoming ``array([1., 2.])`` -- is still + a failure rather than two numbers that happen to be equal. + + Whitespace is dropped entirely. numpy pads an array's columns to its + widest element, so shortening one number moves the spaces around + every other: ``[ 6.32508961 17.37701969]`` against + ``[ 6.3250866 17.377018 ]``. Those spaces carry no meaning the + numeric comparison below has not already made. + """ + + def mark(match: re.Match) -> str: + token = match.group(0) + return "~f" if ("." in token or "e" in token or "E" in token) else "~i" + + return "".join(_NUMBER.sub(mark, text).split()) + + +def _numerically_equal(want: str, got: str) -> bool: + if _skeleton(want) != _skeleton(got): + return False + wants = _NUMBER.findall(want) + gots = _NUMBER.findall(got) + if not wants or len(wants) != len(gots): + return False + return all( + math.isclose(float(w), float(g), rel_tol=RTOL, abs_tol=ATOL) + for w, g in zip(wants, gots) + ) + + +_text_check_output = doctest.OutputChecker.check_output + + +def _check_output(self, want, got, optionflags): + if _text_check_output(self, want, got, optionflags): + return True + return _numerically_equal(want, got) + + +# Patched on the base class rather than installed as a checker: pytest +# builds its own ``LiteralsOutputChecker`` subclass and calls up to this +# method, so overriding here survives both plain ``doctest`` and pytest, +# and does not depend on pytest's internals. +_patched = _check_output # type: ignore[assignment] +doctest.OutputChecker.check_output = _patched # type: ignore[method-assign] + + OPT_IN = { "ml": ( "--run-ml", diff --git a/docs/changelog.rst b/docs/changelog.rst index 280e01cd..663670e8 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -51,16 +51,39 @@ v0.19.1 (unreleased) was passed". That is now true. The 0-d array came from ``np.where``, which does not collapse. - Two doctest option flags are set in ``pyproject.toml``. - ``NORMALIZE_WHITESPACE``, because numpy picks its own line breaks and - column padding for an array and both move with the widest element -- - without it an example is only correct at the exact wrapping it was - captured at. ``ELLIPSIS``, so an example that ends in a fit can write - ``529.05371...`` rather than all seventeen digits: the trailing digits - of an optimiser's output are not part of what the docstring is - promising, and they move with the BLAS and the platform. Array reprs - are left exact -- numpy already prints only eight significant digits - there. + **The numbers in the examples are compared as numbers.** doctest + compares printed output as text, which is the wrong test for a library + whose examples end in an optimiser: the same ``Duane`` fit lands on + ``b = 4.1995e-05`` under Python 3.11 and ``4.2032e-05`` under 3.12, + and numpy prints eight significant digits either way. Sixteen of the + 229 examples disagree between those two Pythons somewhere in their + digits. + + The obvious workaround -- trimming each documented number back to the + digits that agree everywhere -- makes the docstring show something the + reader's own session will not produce, which is precisely what these + examples exist to avoid. So the examples record the real output, in + full, and ``conftest.py`` installs a fallback comparison that runs + only after the ordinary text comparison has failed. It fires when the + two outputs are identical apart from their numeric literals -- same + words, same brackets, same integer-versus-float shape, so ``1`` never + matches ``1.`` and a dtype change is still a failure -- and then + compares the numbers with ``rel_tol=1e-3``, set by the loosest genuine + disagreement between supported Pythons with no margin beyond it, and + ``abs_tol=1e-12`` for a restoration factor whose true value is zero + and which surfaces as ``1e-16`` with whatever mantissa the optimiser + stopped on. + + What that forgives is a value drifting inside the tolerance. What it + still catches is every defect listed above: a stale value from another + parameterisation, the wrong function being called, the wrong shape, an + exception, a missing import. ``surpyval/tests/test_doctest_checker.py`` + pins both halves of that, using the real output pairs observed on + different Pythons. + + ``NORMALIZE_WHITESPACE`` is set in ``pyproject.toml`` for the same + reason: numpy picks its own line breaks and column padding for an + array and both move with the width of the widest element. This closes #158. diff --git a/pyproject.toml b/pyproject.toml index bc4849d8..c1eb7d4d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -53,20 +53,17 @@ where = ["."] [tool.pytest.ini_options] # Applied to the ``--doctest-modules`` run (see .github/workflows/actions.yml), # which executes the ``>>>`` examples in the docstrings and compares their -# printed output exactly. +# printed output. # -# NORMALIZE_WHITESPACE: numpy chooses its own line breaks and column -# padding when it reprs an array, and both change with the width of the -# widest element. Without this an example is only correct for the exact -# terminal wrapping it was captured at, and a docstring wrapped to 79 -# columns by hand -- as most of these are -- can never match. +# numpy chooses its own line breaks and column padding when it reprs an +# array, and both change with the width of the widest element. Without +# this an example is only correct for the exact terminal wrapping it was +# captured at, and a docstring wrapped to 79 columns by hand -- as most +# of these are -- can never match. # -# ELLIPSIS: lets an example that ends in a fit write ``529.05...`` -# instead of all seventeen digits. The last few digits of an optimiser's -# output are not part of the promise the docstring is making, and they -# move with the BLAS, the platform and any change to the fitting -# ladder. -doctest_optionflags = "NORMALIZE_WHITESPACE ELLIPSIS" +# The other half of the comparison is in conftest.py, which compares the +# numbers in the output as numbers rather than as text. +doctest_optionflags = "NORMALIZE_WHITESPACE" # For Black pre-commit hook [tool.black] diff --git a/surpyval/degradation/degradation_analysis.py b/surpyval/degradation/degradation_analysis.py index f6370e29..0197fd5c 100644 --- a/surpyval/degradation/degradation_analysis.py +++ b/surpyval/degradation/degradation_analysis.py @@ -1193,8 +1193,8 @@ class DegradationAnalysis_: Censored Units : 0 Life Distribution : Weibull Parameters : - alpha: 441.47809... - beta: 6.9870788... + alpha: 441.47809611105606 + beta: 6.987078889297555 >>> model.pseudo_failure_times array([451.61290323, 500. , 318.18181818, 378.37837838]) """ diff --git a/surpyval/recurrent/parametric/cox_lewis.py b/surpyval/recurrent/parametric/cox_lewis.py index 140a7ad2..e66f9cdf 100644 --- a/surpyval/recurrent/parametric/cox_lewis.py +++ b/surpyval/recurrent/parametric/cox_lewis.py @@ -27,8 +27,8 @@ class CoxLewis(NHPPFitter): Process : Cox-Lewis Fitted by : MLE Parameters : - alpha: 0.38481273... - beta: 0.19396672... + alpha: 0.384812737762836 + beta: 0.19396672109211047 >>> model.cif([1, 2, 3, 4, 5, 6]) array([ 1.62151879, 3.59013322, 5.98014113, 8.88174429, 12.40445268, 16.6812175 ]) diff --git a/surpyval/recurrent/parametric/crow_amsaa.py b/surpyval/recurrent/parametric/crow_amsaa.py index c427c654..2bd8eac7 100644 --- a/surpyval/recurrent/parametric/crow_amsaa.py +++ b/surpyval/recurrent/parametric/crow_amsaa.py @@ -27,8 +27,8 @@ class CrowAMSAA(NHPPFitter): Process : Crow-AMSAA Fitted by : MLE Parameters : - alpha: 913.84662... - beta: 1.4781707... + alpha: 913.8466210685444 + beta: 1.4781707110680866 >>> model.cif([1, 2, 3, 4, 5, 6]) array([4.20072057e-05, 1.17030084e-04, 2.13103439e-04, 3.26040266e-04, 4.53440995e-04, 5.93696079e-04]) diff --git a/surpyval/recurrent/parametric/duane.py b/surpyval/recurrent/parametric/duane.py index 23bb1f2e..0974d1ca 100644 --- a/surpyval/recurrent/parametric/duane.py +++ b/surpyval/recurrent/parametric/duane.py @@ -27,7 +27,7 @@ class Duane(NHPPFitter): Process : Duane Fitted by : MLE Parameters : - alpha: 1.4782020... + alpha: 1.478202089169939 b: 4.199455086392048e-05 >>> model.cif([1, 2, 3, 4, 5, 6]) array([4.19945509e-05, 1.16997373e-04, 2.13046585e-04, 3.25956224e-04, diff --git a/surpyval/recurrent/parametric/hpp.py b/surpyval/recurrent/parametric/hpp.py index e5c7d78d..43cbdf62 100644 --- a/surpyval/recurrent/parametric/hpp.py +++ b/surpyval/recurrent/parametric/hpp.py @@ -34,7 +34,7 @@ class HPP(CountingProcess): Process : Homogeneous Poisson Process Fitted by : MLE Parameters : - lambda: 0.0023047023... + lambda: 0.0023047023327236213 >>> model.cif([1, 2, 3, 4, 5, 6]) array([0.0023047 , 0.0046094 , 0.00691411, 0.00921881, 0.01152351, 0.01382821]) diff --git a/surpyval/recurrent/renewal/generalized_one_renewal.py b/surpyval/recurrent/renewal/generalized_one_renewal.py index 60b8b303..a9d1937c 100644 --- a/surpyval/recurrent/renewal/generalized_one_renewal.py +++ b/surpyval/recurrent/renewal/generalized_one_renewal.py @@ -60,10 +60,10 @@ class GeneralizedOneRenewal(RenewalFitMixin): ========================= Distribution : Weibull Fitted by : MLE - Restoration Factor : -0.17301846... + Restoration Factor : -0.1730184624683848 Parameters : - alpha: 1.3919045... - beta: 5.0088611... + alpha: 1.3919045968817332 + beta: 5.008861189641614 >>> >>> np.random.seed(0) >>> np_model = model.count_terminated_simulation(len(x), 5000) @@ -208,10 +208,10 @@ def fit_from_recurrent_data(self, data, dist=Weibull, init=None): ========================= Distribution : Weibull Fitted by : MLE - Restoration Factor : 0.34027890... + Restoration Factor : 0.3402789091696592 Parameters : - alpha: 1.4115217... - beta: 3.5499343... + alpha: 1.4115217370254167 + beta: 3.5499343659245564 """ self._check_dist_eligible(dist) validate_renewal_censoring(data.c, type(self).__name__) @@ -290,10 +290,10 @@ def fit(self, x, i=None, c=None, n=None, dist=Weibull, init=None): ========================= Distribution : Weibull Fitted by : MLE - Restoration Factor : 0.34027890... + Restoration Factor : 0.3402789091696592 Parameters : - alpha: 1.4115217... - beta: 3.5499343... + alpha: 1.4115217370254167 + beta: 3.5499343659245564 """ data = handle_xicn(x, i, c, n) return self.fit_from_recurrent_data(data, dist=dist, init=init) diff --git a/surpyval/recurrent/renewal/generalized_renewal.py b/surpyval/recurrent/renewal/generalized_renewal.py index a72dadf6..eaa40bcf 100644 --- a/surpyval/recurrent/renewal/generalized_renewal.py +++ b/surpyval/recurrent/renewal/generalized_renewal.py @@ -60,10 +60,10 @@ class GeneralizedRenewal(RenewalFitMixin): Distribution : Weibull Fitted by : MLE Kijima Type : i - Restoration Factor : 0.15732122... + Restoration Factor : 0.15732122999163628 Parameters : - alpha: 1.2613379... - beta: 8.9390232... + alpha: 1.261337933121844 + beta: 8.93902321971521 >>> >>> np.random.seed(0) >>> np_model = model.count_terminated_simulation(len(x), 5000) @@ -261,8 +261,8 @@ def fit_from_recurrent_data( Kijima Type : i Restoration Factor : 1.3316262291443964e-16 Parameters : - alpha: 2.3990296... - beta: 2.7539200... + alpha: 2.399029668688425 + beta: 2.753920042066547 """ validate_renewal_censoring(data.c, type(self).__name__) reject_left_truncation(data, type(self).__name__) @@ -346,8 +346,8 @@ def fit( Kijima Type : i Restoration Factor : 1.3316262291443964e-16 Parameters : - alpha: 2.3990296... - beta: 2.7539200... + alpha: 2.399029668688425 + beta: 2.753920042066547 """ data = handle_xicn(x, i, c, n) return self.fit_from_recurrent_data(data, dist, kijima, init=init) diff --git a/surpyval/tests/test_doctest_checker.py b/surpyval/tests/test_doctest_checker.py new file mode 100644 index 00000000..0c9eaa78 --- /dev/null +++ b/surpyval/tests/test_doctest_checker.py @@ -0,0 +1,119 @@ +"""The numeric fallback used by the ``--doctest-modules`` run. + +``conftest._numerically_equal`` decides whether two blocks of doctest +output that differ as text are the same as numbers. It is the mechanism +that lets the docstring examples record real, untrimmed output while +still passing on every supported Python, so it needs its own tests: too +strict and the doctest step fails on a toolchain difference, too loose +and it stops catching the documentation drift it was added to catch. + +The "tolerated" cases below are not invented. Each is a real pair of +outputs seen for the same example on two different Pythons in CI. +""" + +import pytest + +from conftest import _numerically_equal + +# (documented output, output seen on another Python) +TOLERATED = [ + pytest.param( + " alpha: 913.84662107\n beta: 1.4781707110680866\n", + " alpha: 913.8468395959314\n beta: 1.4781709287931744\n", + id="crow-amsaa-fit", + ), + pytest.param( + # The loosest genuine disagreement found: 9e-4 relative on ``b``. + " alpha: 1.478202089169939\n b: 4.199455086392048e-05\n", + " alpha: 1.4781024854527343\n b: 4.203204430005199e-05\n", + id="duane-fit", + ), + pytest.param( + " beta_2 : -0.02147901865302258\n", + " beta_2 : -0.021479544898597686\n", + id="hpp-proportional-intensity-coefficient", + ), + pytest.param( + # A restoration factor whose true value is zero; the optimiser + # stops on whatever mantissa it stops on. Only ``abs_tol`` can + # rescue this one. + "Restoration Factor : 1.3316262291443964e-16\n", + "Restoration Factor : 1.1551809284521243e-16\n", + id="numerical-zero", + ), + pytest.param( + "1.8227536487527594\n", + "1.822753648752769\n", + id="scipy-special-last-digits", + ), + pytest.param( + # numpy re-pads the columns when an element gets shorter, so the + # closing bracket moves. + " alpha: [ 6.32508961 17.37701969]\n", + " alpha: [ 6.3250866 17.377018 ]\n", + id="array-repadded", + ), +] + +# Every one of these is a real defect this sweep found in the docstrings, +# or the class of defect it found. None may be forgiven. +REJECTED = [ + pytest.param( + "3\n", + "3.332162203618775\n", + id="value-from-another-parameterisation", + ), + pytest.param( + "11.229\n", "10.533288486847923\n", id="documented-variance-was-wrong" + ), + pytest.param( + "array([1, 2, 3, 4, 5])\n", + "array([1., 2., 3., 4., 5.])\n", + id="dtype-changed", + ), + pytest.param( + "array([0.83333333, 0.66666667])\n", + "array([0.16666667, 0.33333333])\n", + id="example-called-the-wrong-function", + ), + pytest.param( + "array([1.0, 2.0])\n", "array([1.0, 2.0, 3.0])\n", id="shape-changed" + ), + pytest.param("alpha: 1.0\n", "beta: 1.0\n", id="different-label"), + pytest.param("array([1.0])\n", "3.0\n", id="scalar-against-array"), + pytest.param("", "1.0\n", id="expected-nothing"), + pytest.param("1.0\n", "", id="produced-nothing"), + pytest.param( + # 2e-3 relative, twice the tolerance. + "8.929795115692489\n", + "8.947654505923874\n", + id="drifted-past-the-tolerance", + ), +] + + +@pytest.mark.parametrize("want, got", TOLERATED) +def test_the_same_number_under_a_different_toolchain_is_accepted(want, got): + assert _numerically_equal(want, got) + + +@pytest.mark.parametrize("want, got", REJECTED) +def test_a_real_difference_is_still_a_failure(want, got): + assert not _numerically_equal(want, got) + + +def test_output_with_no_numbers_is_not_silently_accepted(): + # Nothing to compare numerically, so the fallback must decline and + # leave the verdict to the ordinary text comparison. + assert not _numerically_equal("Fitted by : MLE\n", "Fitted by : MPS\n") + assert not _numerically_equal("Fitted by : MLE\n", "Fitted by : MLE\n") + + +def test_the_checker_is_installed_on_the_doctest_base_class(): + # The examples' precision depends on this patch being live for the + # whole doctest run, including under pytest's own subclass. + import doctest + + checker = doctest.OutputChecker() + assert checker.check_output("1.0000000\n", "1.0000001\n", 0) + assert not checker.check_output("1.0\n", "2.0\n", 0) diff --git a/surpyval/univariate/information_criteria.py b/surpyval/univariate/information_criteria.py index f9bd852e..d45960c3 100644 --- a/surpyval/univariate/information_criteria.py +++ b/surpyval/univariate/information_criteria.py @@ -59,7 +59,7 @@ def neg_ll(self) -> float: >>> x = Weibull.random(100, 10, 3) >>> model = Weibull.fit(x) >>> model.neg_ll() - 262.52685... + 262.52685642390634 """ if getattr(self, "data", None) is None: raise ValueError("Must have been fit with data") @@ -88,7 +88,7 @@ def bic(self) -> float: >>> x = Weibull.random(100, 10, 3) >>> model = Weibull.fit(x) >>> model.bic() - np.float64(534.26405...) + np.float64(534.2640532197888) References ---------- @@ -124,7 +124,7 @@ def aic(self) -> float: >>> x = Weibull.random(100, 10, 3) >>> model = Weibull.fit(x) >>> model.aic() - 529.05371... + 529.0537128478127 """ if hasattr(self, "_aic"): return self._aic @@ -152,7 +152,7 @@ def aic_c(self) -> float: >>> x = Weibull.random(100, 10, 3) >>> model = Weibull.fit(x) >>> model.aic_c() - np.float64(529.17742...) + np.float64(529.1774241880189) """ if hasattr(self, "_aic_c"): return self._aic_c diff --git a/surpyval/univariate/parametric/mixture_model.py b/surpyval/univariate/parametric/mixture_model.py index 4d36919d..2a6f80ea 100755 --- a/surpyval/univariate/parametric/mixture_model.py +++ b/surpyval/univariate/parametric/mixture_model.py @@ -297,8 +297,8 @@ def fit( Sub-Distributions : 2 Fitted by : EM Weights : - 0.61848918..., - 0.38151081... + 0.6184891886499861, + 0.381510811350014 Parameters : alpha: [ 6.32508961 17.37701969] beta: [ 1.83105154 12.01392721] diff --git a/surpyval/univariate/parametric/parametric_fitter.py b/surpyval/univariate/parametric/parametric_fitter.py index f01cde16..3c52b9bc 100755 --- a/surpyval/univariate/parametric/parametric_fitter.py +++ b/surpyval/univariate/parametric/parametric_fitter.py @@ -762,8 +762,8 @@ def fit( Distribution : Weibull Fitted by : MLE Parameters : - alpha: 9.8150187... - beta: 3.7987404... + alpha: 9.815018791049368 + beta: 3.798740470368033 >>> Weibull.fit(x, how='MPS', fixed={'alpha' : 10}) Parametric SurPyval Model ========================= @@ -771,7 +771,7 @@ def fit( Fitted by : MPS Parameters : alpha: 10.0 - beta: 3.6707965... + beta: 3.670796510564323 >>> Weibull.fit(xl=np.floor(x), xr=np.ceil(x), how='MPP', ... heuristic='Turnbull') Parametric SurPyval Model @@ -779,8 +779,8 @@ def fit( Distribution : Weibull Fitted by : MPP Parameters : - alpha: 9.9501683... - beta: 3.2119714... + alpha: 9.950168329892755 + beta: 3.211971411540382 >>> c = np.zeros_like(x) >>> c[x > 13] = 1 >>> x[x > 13] = 13 @@ -792,8 +792,8 @@ def fit( Distribution : Weibull Fitted by : MLE Parameters : - alpha: 9.8935844... - beta: 3.7868860... + alpha: 9.893584496413128 + beta: 3.78688602908912 """ surv_data = SurpyvalData( @@ -894,10 +894,10 @@ def fit_from_df( ========================= Distribution : Weibull Fitted by : MLE - Offset (gamma) : 39.765577... + Offset (gamma) : 39.76557772434183 Parameters : - alpha: 7.1419836... - beta: 2.6204759... + alpha: 7.141983615103902 + beta: 2.62047590823775 """ if not isinstance(df, pd.DataFrame): diff --git a/surpyval/univariate/regression/proportional_hazards/proportional_hazards_fitter.py b/surpyval/univariate/regression/proportional_hazards/proportional_hazards_fitter.py index d62f1ac5..8497d0bc 100644 --- a/surpyval/univariate/regression/proportional_hazards/proportional_hazards_fitter.py +++ b/surpyval/univariate/regression/proportional_hazards/proportional_hazards_fitter.py @@ -204,13 +204,13 @@ def fit( Regression Model : Log Linear [e^(beta'Z)] Fitted by : MLE Distribution : - alpha: 0.24255136... - beta: 16.057785... + alpha: 0.2425513627560218 + beta: 16.057785182711932 Regression Model : - beta_0: -9.1650627... - beta_1: -7.9985730... - beta_2: -27.503185... - beta_3: 18.385445... + beta_0: -9.165062726518311 + beta_1: -7.998573055929788 + beta_2: -27.50318580568538 + beta_3: 18.385445332039488 >>> model = WeibullPH.fit(x=x, Z=Z, c=c, fixed={"beta": 15}) >>> model Parametric Regression SurPyval Model @@ -220,13 +220,13 @@ def fit( Regression Model : Log Linear [e^(beta'Z)] Fitted by : MLE Distribution : - alpha: 0.23772966... - beta: 15.0 + alpha: 0.237729668424067 + beta: 15.0 Regression Model : - beta_0: -8.6283269... - beta_1: -7.6175293... - beta_2: -25.952367... - beta_3: 17.270148... + beta_0: -8.62832691738283 + beta_1: -7.617529362323243 + beta_2: -25.952367249502934 + beta_3: 17.270148387391387 """ data, prep = prepare_regression_fit( self, From 96483bcac3de33ee4734fc7be70c2cfe21be9a7b Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 4 Aug 2026 22:02:56 +0000 Subject: [PATCH 08/13] Strip before the numeric doctest comparison The doctest step failed on all three Pythons, on one example: ProportionalIntensityHPP, whose numbers were well inside tolerance. The fallback never got to look at them. Its expected output carries markers, which the text comparison substitutes for empty lines before matching; the fallback did not, so the two skeletons differed on the marker text and it declined before comparing a single number. Why local verification missed it: the fallback only runs once the text comparison has failed, and on this machine that example's fit matches the recorded value exactly. Running against real 3.11, 3.12 and 3.13 interpreters exercised the fallback on the sixteen examples that drift *here* -- and this was not one of them. Passing on three Pythons said nothing about the examples whose numbers agreed on all three. So --doctest-force-numeric is added alongside the fix: it routes every example whose output contains a number through the numeric comparison regardless of whether the text matched, exercising the fallback against all 229 rather than today's accidental few. Outputs with no numbers keep the text comparison; there is nothing in them to compare. CI runs the doctest step a second time under the flag, which costs fifteen seconds and turns "the fallback handles the examples that happen to drift on this machine" into "the fallback handles the examples". The regression test uses the ProportionalIntensityHPP block whole, markers and all, rather than one extracted coefficient -- an extracted line would have gone on passing, since what broke was not a number. Verified on 3.11, 3.12 and 3.13: 229 doctests pass in both modes on each. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts --- .github/workflows/actions.yml | 12 +++++++ conftest.py | 47 +++++++++++++++++++++++- docs/changelog.rst | 10 ++++++ surpyval/tests/test_doctest_checker.py | 50 ++++++++++++++++++++++++-- 4 files changed, 116 insertions(+), 3 deletions(-) diff --git a/.github/workflows/actions.yml b/.github/workflows/actions.yml index 4fc6629c..eee43ccc 100644 --- a/.github/workflows/actions.yml +++ b/.github/workflows/actions.yml @@ -100,6 +100,18 @@ jobs: python -m pytest --doctest-modules surpyval --ignore=surpyval/tests --ignore=surpyval/alpha + # The numeric comparison above only runs for an example whose + # output has actually drifted -- a handful on any one machine, and + # a different handful on each. That leaves a gap in it invisible + # until CI hits the one example that needed it. This second pass + # routes every example through it, so the fallback is checked + # against all of them. Fifteen seconds. + - name: doctests (numeric comparison forced) + run: + python -m pytest --doctest-modules surpyval + --ignore=surpyval/tests --ignore=surpyval/alpha + --doctest-force-numeric + - name: coverage run: | coverage report diff --git a/conftest.py b/conftest.py index 4e8babc3..227d5e0b 100644 --- a/conftest.py +++ b/conftest.py @@ -1,4 +1,8 @@ -"""Opt-in gating for the slow parts of the test suite. +"""Suite-wide fixtures: doctest number comparison, and opt-in gating. + +The first half of this file makes the ``--doctest-modules`` run compare +the numbers in an example's output as numbers rather than as text; see +the comment above ``RTOL``. The rest is the opt-in gating below. Two groups are skipped unless asked for, because both are expensive and neither guards a regression that the default run would miss quickly: @@ -72,6 +76,7 @@ ATOL = 1e-12 _NUMBER = re.compile(r"[-+]?(?:\d+\.\d*|\.\d+|\d+)(?:[eE][-+]?\d+)?") +_BLANKLINE = re.compile(r"(?m)^%s\s*?$" % re.escape(doctest.BLANKLINE_MARKER)) def _skeleton(text: str) -> str: @@ -96,6 +101,12 @@ def mark(match: re.Match) -> str: def _numerically_equal(want: str, got: str) -> bool: + # ```` stands for an empty line in the expected output. + # The text comparison substitutes it before matching, so this one has + # to as well, or a model repr with a blank line in it can never reach + # the numeric comparison at all. + want = _BLANKLINE.sub("", want) + if _skeleton(want) != _skeleton(got): return False wants = _NUMBER.findall(want) @@ -125,6 +136,26 @@ def _check_output(self, want, got, optionflags): doctest.OutputChecker.check_output = _patched # type: ignore[method-assign] +def _forced_check_output(self, want, got, optionflags): + """As above, but the numeric path is the *only* path. + + The fallback normally runs only when an example's output has + actually drifted, which on any one machine is a handful of them. A + gap in it -- the ```` markers it did not strip, say -- + therefore stays invisible locally and surfaces in CI, on whichever + Python happens to compute a different last digit. + + Under ``--doctest-force-numeric`` every example whose output + contains a number is compared numerically instead, so the fallback + is exercised against all 229 of them rather than against today's + accidental few. Outputs with no numbers keep the text comparison; + there is nothing in them for this to compare. + """ + if not _NUMBER.search(want): + return _text_check_output(self, want, got, optionflags) + return _numerically_equal(want, got) + + OPT_IN = { "ml": ( "--run-ml", @@ -147,6 +178,16 @@ def pytest_addoption(parser): default=False, help=f"run the {description} (skipped by default)", ) + parser.addoption( + "--doctest-force-numeric", + action="store_true", + default=False, + help=( + "compare every doctest example's numbers numerically, not " + "only those whose text has drifted; exercises the fallback " + "against all of them" + ), + ) def pytest_configure(config): @@ -154,6 +195,10 @@ def pytest_configure(config): config.addinivalue_line( "markers", f"{mark}: {description}; opt in with {flag}" ) + if config.getoption("--doctest-force-numeric"): + doctest.OutputChecker.check_output = ( # type: ignore[method-assign] + _forced_check_output # type: ignore[assignment] + ) def pytest_collection_modifyitems(config, items): diff --git a/docs/changelog.rst b/docs/changelog.rst index 663670e8..4c2cdb79 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -81,6 +81,16 @@ v0.19.1 (unreleased) pins both halves of that, using the real output pairs observed on different Pythons. + The fallback only runs when an example has actually drifted, which on + any one machine is a handful of them -- and a different handful on + each. A gap in it is therefore invisible locally and surfaces in CI, + on whichever Python computed a different last digit. So the doctest + step runs twice: once normally, and once under + ``--doctest-force-numeric``, which routes every example whose output + contains a number through the numeric comparison. Fifteen seconds, and + the fallback is exercised against all 229 examples rather than + today's accidental few. + ``NORMALIZE_WHITESPACE`` is set in ``pyproject.toml`` for the same reason: numpy picks its own line breaks and column padding for an array and both move with the width of the widest element. diff --git a/surpyval/tests/test_doctest_checker.py b/surpyval/tests/test_doctest_checker.py index 0c9eaa78..9193ff3f 100644 --- a/surpyval/tests/test_doctest_checker.py +++ b/surpyval/tests/test_doctest_checker.py @@ -11,10 +11,41 @@ outputs seen for the same example on two different Pythons in CI. """ +import doctest + import pytest from conftest import _numerically_equal +# The ``ProportionalIntensityHPP`` example, verbatim: the documented +# block on the left, the block CI produced on the right. Kept whole +# rather than reduced to one coefficient because the thing that broke +# here was not a number -- it was the ```` markers, which the +# text comparison strips from the expected output before matching and +# which the numeric fallback has to strip too. A test on an extracted +# line would have gone on passing. +HPP_DOCUMENTED = """\ +Base Rate Parameters: + lambda : 0.012395105741757225 + +Covariate Coefficients: + beta_0 : 0.06397367067847898 + beta_1 : 0.011491178797116433 + beta_2 : -0.02147901865302258 + +""" + +HPP_OBSERVED = """\ +Base Rate Parameters: + lambda : 0.012395109943236718 + +Covariate Coefficients: + beta_0 : 0.06397361219411789 + beta_1 : 0.011491197469342556 + beta_2 : -0.021479544898597686 + +""" + # (documented output, output seen on another Python) TOLERATED = [ pytest.param( @@ -53,6 +84,7 @@ " alpha: [ 6.3250866 17.377018 ]\n", id="array-repadded", ), + pytest.param(HPP_DOCUMENTED, HPP_OBSERVED, id="repr-with-blank-lines"), ] # Every one of these is a real defect this sweep found in the docstrings, @@ -112,8 +144,22 @@ def test_output_with_no_numbers_is_not_silently_accepted(): def test_the_checker_is_installed_on_the_doctest_base_class(): # The examples' precision depends on this patch being live for the # whole doctest run, including under pytest's own subclass. - import doctest - checker = doctest.OutputChecker() assert checker.check_output("1.0000000\n", "1.0000001\n", 0) assert not checker.check_output("1.0\n", "2.0\n", 0) + + +def test_the_whole_comparison_path_handles_blank_lines(): + # Not the fallback in isolation but the method doctest actually + # calls, with the flags the doctest step actually runs under. The + # ```` handling lives in the text comparison, so only + # this route proves the two halves agree about it. + checker = doctest.OutputChecker() + assert checker.check_output( + HPP_DOCUMENTED, HPP_OBSERVED, doctest.NORMALIZE_WHITESPACE + ) + assert not checker.check_output( + HPP_DOCUMENTED, + HPP_OBSERVED.replace("0.0123951", "0.0198765"), + doctest.NORMALIZE_WHITESPACE, + ) From 42c98047728f5af07ce7a48b8eefbe1ba952358b Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 4 Aug 2026 23:41:14 +0000 Subject: [PATCH 09/13] Make the docs toolchain a `docs` extra `pip install -e ".[docs]"` now installs everything needed to build the documentation, alongside the `tests` extra that was already there. docs/requirements.txt is removed rather than kept alongside the extra: two copies of the same pinned toolchain is the arrangement that drifts, and the repo already settled this pattern for tests, where requirements_dev.txt is `-e .[tests]` plus tools rather than a second pin list. Read the Docs installs the extra directly via extra_requirements, which is their documented form for exactly this; Contributing.rst loses a step. The pins are carried over unchanged, including the ipykernel==6.31.0 cap and the reason for it. matplotlib is not repeated in the extra -- it is a runtime dependency of the package, installed alongside. Verified by a complete `sphinx -b html` in a clean 3.12 virtualenv built only from `pip install ".[docs]"`. Part of #141. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts --- .readthedocs.yaml | 3 ++- docs/Contributing.rst | 3 +-- docs/changelog.rst | 16 ++++++++++++++++ docs/requirements.txt | 17 ----------------- pyproject.toml | 23 +++++++++++++++++++++++ 5 files changed, 42 insertions(+), 20 deletions(-) delete mode 100644 docs/requirements.txt diff --git a/.readthedocs.yaml b/.readthedocs.yaml index 3b8c8153..c014b91d 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -34,4 +34,5 @@ python: install: - method: pip path: . - - requirements: docs/requirements.txt \ No newline at end of file + extra_requirements: + - docs \ No newline at end of file diff --git a/docs/Contributing.rst b/docs/Contributing.rst index e39cc43a..e2dfdb45 100644 --- a/docs/Contributing.rst +++ b/docs/Contributing.rst @@ -41,8 +41,7 @@ To build the documentation locally: .. code-block:: bash - pip install -e . - pip install -r docs/requirements.txt + pip install -e ".[docs]" sphinx-build -b html docs docs/_build/html When writing documentation, prefer ``.. jupyter-execute::`` over static diff --git a/docs/changelog.rst b/docs/changelog.rst index 4c2cdb79..e526d7ae 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -4,6 +4,22 @@ Changelog v0.19.1 (unreleased) -------------------- +- **The documentation toolchain is a ``docs`` extra.** + ``pip install -e ".[docs]"`` now installs everything needed to build + the documentation, alongside the ``tests`` extra that was already + there. ``docs/requirements.txt`` is gone: its pins moved into + ``pyproject.toml`` verbatim, and Read the Docs installs the extra + directly via ``extra_requirements``. Keeping both would have meant two + copies of the same pinned toolchain, which is the arrangement that + drifts. + + The pins are unchanged, including the ``ipykernel==6.31.0`` cap and + the reason for it -- jupyter-sphinx notebook execution dies against + the ipykernel 7 line. ``matplotlib`` is not repeated in the extra; it + is a runtime dependency of the package, which is installed alongside. + + Part of #141. + - **CI now runs the docstring examples.** ``pytest --doctest-modules`` over the package is a new step in the deployment workflow, and every one of the 229 docstring examples passes. It was 59 failing tests diff --git a/docs/requirements.txt b/docs/requirements.txt deleted file mode 100644 index 22d611ee..00000000 --- a/docs/requirements.txt +++ /dev/null @@ -1,17 +0,0 @@ -# Documentation build toolchain (used by ReadTheDocs and local builds). -# -# This is a coherent, fully pinned set validated by a complete docs -# build (every jupyter-execute cell runs) in a clean virtualenv on top -# of `pip install .`. Bump versions together and re-validate with: -# python -m sphinx -b html docs docs/_build/html -# -# ipykernel is capped below 7: jupyter-sphinx notebook execution hangs -# or crashes (zmq "Socket operation on non-socket") against the -# ipykernel 7 line. -matplotlib>=3.10 -sphinx==8.2.3 -sphinx_rtd_theme==3.1.0 -sphinx-copybutton==0.5.2 -sphinx-notfound-page==1.1.0 -jupyter-sphinx==0.5.3 -ipykernel==6.31.0 diff --git a/pyproject.toml b/pyproject.toml index c1eb7d4d..67545eea 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,6 +39,29 @@ tests = [ # (MongoDB compatibility); pymongo provides the bson package. "pymongo>=4", ] +docs = [ + # The documentation build toolchain, used by Read the Docs and by a + # local ``sphinx-build``. Fully pinned rather than bounded: the docs + # execute every ``.. jupyter-execute::`` cell during the build, so a + # toolchain change can break the build in ways a version range would + # let through silently. Bump these together and re-validate with a + # complete build: + # + # pip install -e ".[docs]" + # python -m sphinx -b html docs docs/_build/html + # + # matplotlib is not repeated here -- it is a runtime dependency of + # the package itself, which the docs install alongside these. + "sphinx==8.2.3", + "sphinx_rtd_theme==3.1.0", + "sphinx-copybutton==0.5.2", + "sphinx-notfound-page==1.1.0", + "jupyter-sphinx==0.5.3", + # Capped below 7: jupyter-sphinx notebook execution hangs or dies + # with a zmq "Socket operation on non-socket" against the ipykernel + # 7 line. + "ipykernel==6.31.0", +] [project.urls] Homepage = "https://github.com/derrynknife/SurPyval" From 03c6106b734ce065e62ac6a4f0fc90e8fa4b3131 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 4 Aug 2026 23:41:34 +0000 Subject: [PATCH 10/13] Fix the docs build broken by the Gamma MPP removal The offset-threshold section of "Parametric SurPyval Modelling" ran a jupyter-execute cell looping over ['MPP', 'MOM', 'MSE', 'MPS', 'MLE'] for a shifted Gamma. Gamma.fit(how="MPP") now raises, so that cell raised, and since documentation cells execute during the build the whole build failed. Nothing caught it. CI does not build the documentation, and Read the Docs builds only master and tags, so this would have surfaced as a failed hosted build at the next release rather than on the pull request that caused it. It was found by running a build to validate the docs extra. The prose around the cell had gone stale in the same way: it described the multi-start probability-plotting search that the removal deleted, and quoted an MPP tolerance from test_offset_divergence.py that no longer exists. It now explains why the Gamma has no probability plot -- the shape sits inside the regularised incomplete gamma rather than outside as an exponent, so the only straight-line axis is the inverse incomplete gamma, which needs the shape being estimated -- and notes that plot() is unaffected because by then the parameters are known. Verified by a complete build: succeeded, with only the 18 pre-existing warnings (duplicate changelog labels, the rtd-theme deprecation, and the ProportionalIntensityNHPP_ autodoc imports). Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts --- docs/Parametric SurPyval Modelling.rst | 16 +++++++++++++--- docs/changelog.rst | 24 ++++++++++++++++++++++++ 2 files changed, 37 insertions(+), 3 deletions(-) diff --git a/docs/Parametric SurPyval Modelling.rst b/docs/Parametric SurPyval Modelling.rst index 895632c7..26b3f134 100644 --- a/docs/Parametric SurPyval Modelling.rst +++ b/docs/Parametric SurPyval Modelling.rst @@ -393,19 +393,29 @@ In practice that sensitivity is the estimator's problem to solve, not yours. Shi x = surv.Gamma.random(10_000, 3.0, 2.0) + 10.0 print('Truth : gamma=10.000, alpha=3.000, beta=2.000') - for how in ['MPP', 'MOM', 'MSE', 'MPS', 'MLE']: + for how in ['MOM', 'MSE', 'MPS', 'MLE']: m = surv.Gamma.fit(x, offset=True, how=how) print('{:6s}: gamma={:.3f}, alpha={:.3f}, beta={:.3f}'.format( how, m.gamma, *m.params)) -This was not always so. Earlier versions could land on an absurd tuple - a negative ``gamma`` with a shape parameter inflated by two orders of magnitude - which was nonetheless an acceptable *distribution*, precisely because of the flat trade-off described above. Two separate causes were at work, and both are now fixed. The probability-plotting search was stranded at a poor local optimum by a single starting shape, so it now tries several and keeps the best correlation. The moment-based initialisers took their moments from the *unshifted* data, where the offset dominates every moment and the shape estimate explodes - 649 for a true shape of 3 - so they now estimate the threshold first and read the remaining parameters off ``x - gamma``. +``MPP`` is absent from that list because the Gamma does not offer it. A +probability plot needs a straight-line y-axis that can be drawn *before* +the parameters are known; the Gamma's CDF is the regularised incomplete +gamma function, with the shape inside the special function rather than +outside as an exponent, so the only such axis is the inverse incomplete +gamma — which needs the shape. To draw the axis you need the answer. +``Gamma.fit(x, how="MPP")`` raises rather than guessing a shape to draw +the axis with; ``Gamma.plot()`` is unaffected, since by then the fitted +parameters are in hand. + +This was not always so. Earlier versions could land on an absurd tuple - a negative ``gamma`` with a shape parameter inflated by two orders of magnitude - which was nonetheless an acceptable *distribution*, precisely because of the flat trade-off described above. The moment-based initialisers took their moments from the *unshifted* data, where the offset dominates every moment and the shape estimate explodes - 649 for a true shape of 3 - so they now estimate the threshold first and read the remaining parameters off ``x - gamma``. The underlying caution still stands, though, and it is worth keeping in mind for your own data: - **Judge an offset fit by what it predicts, not only by the printed parameters.** Plot it against the non-parametric estimate, or compare the survival function, quantiles, mean and variance. Two parameter tuples that look very different can imply nearly the same distribution. - **If you need ``gamma`` itself to be meaningful** - you are interpreting it as a guaranteed minimum life, say - prefer ``MLE``, which remains the most accurate on the parameters, and treat a single point estimate of a threshold with care regardless of method. -``test_offset_divergence.py`` in the test suite pins this down with measured KL and Wasserstein distances alongside parameter tolerances: ``MLE`` is held to 5% on every parameter, ``MOM`` to 10% and ``MPP`` to 20%, with the implied distributions essentially identical in all three cases. +``test_offset_divergence.py`` in the test suite pins this down with measured KL and Wasserstein distances alongside parameter tolerances: ``MLE`` is held to 5% on every parameter across the offsettable distributions, and ``MOM`` to 10%, with the implied distributions essentially identical either way. Fixing parameters ----------------- diff --git a/docs/changelog.rst b/docs/changelog.rst index e526d7ae..c928d521 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -4,6 +4,30 @@ Changelog v0.19.1 (unreleased) -------------------- +- **Fixed a documentation build broken by the Gamma MPP removal.** The + offset-threshold section of *Parametric SurPyval Modelling* ran a + ``jupyter-execute`` cell looping over + ``['MPP', 'MOM', 'MSE', 'MPS', 'MLE']`` for a shifted Gamma. Since + ``Gamma.fit(how="MPP")`` now raises, that cell raised, and because + documentation cells are executed during the build the whole build + failed. + + Nothing caught it: continuous integration does not build the + documentation, and Read the Docs builds only ``master`` and tags, so + it would have surfaced as a failed hosted build at the next release + rather than on the change that caused it. It was found by running a + build to validate the ``docs`` extra below. + + The prose around the cell had gone stale the same way -- it described + the multi-start probability-plotting search that the removal deleted, + and quoted an ``MPP`` tolerance from ``test_offset_divergence.py`` + that no longer exists. It now explains why the Gamma has no + probability plot at all: the shape sits inside the regularised + incomplete gamma rather than outside as an exponent, so the only + straight-line axis is the inverse incomplete gamma, which needs the + very shape being estimated. ``Gamma.plot()`` is unaffected, since by + then the parameters are known. + - **The documentation toolchain is a ``docs`` extra.** ``pip install -e ".[docs]"`` now installs everything needed to build the documentation, alongside the ``tests`` extra that was already From d4cd31b58de510cd1e52671b58b93218c4b30217 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 5 Aug 2026 00:00:06 +0000 Subject: [PATCH 11/13] Build the documentation in CI on the release pull request The docs execute every `.. jupyter-execute::` cell as they build, so they are a second test suite that exercises the public API for real -- and one that a change touching no documentation file can break. Removing Gamma's probability-plot fitting did exactly that, and nothing noticed, because Read the Docs builds only master and tags: it would have surfaced as a broken hosted build after the release rather than on the pull request that caused it. The job is conditioned on `github.base_ref == 'master'`, which is set only for pull_request events, so it runs on the develop -> master release pull request and nowhere else. Deliberately not on pushes to master: Read the Docs rebuilds there anyway, and by then the gate has nothing left to gate. It matches .readthedocs.yaml rather than the test jobs -- Python 3.12, the package installed via its own `docs` extra -- because the point is to reproduce the hosted build, and it uploads the rendered HTML as an artifact for review on the release PR. Not built with -W. There are 18 pre-existing warnings, mostly duplicate labels from autosectionlabel meeting the changelog's repeated section headings; clearing those and then failing on warning here and in .readthedocs.yaml together is a separate change, and turning it on before then would fail every release. The residual gap is deliberate and now documented: a break introduced on a pull request into develop is caught when the release is prepared, not when it lands. Building on every pull request would cost minutes on each, and a path filter would not have helped here -- the change that broke the build was in gamma.py, not under docs/. Contributing.rst claimed the documentation build already ran once per pull request. It did not run in CI at all. It now describes what runs where, and names the trade-off. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts --- .github/workflows/actions.yml | 51 +++++++++++++++++++++++++++++++++++ docs/Contributing.rst | 17 +++++++++--- docs/changelog.rst | 30 +++++++++++++++++++++ 3 files changed, 94 insertions(+), 4 deletions(-) diff --git a/.github/workflows/actions.yml b/.github/workflows/actions.yml index eee43ccc..e2acf853 100644 --- a/.github/workflows/actions.yml +++ b/.github/workflows/actions.yml @@ -123,3 +123,54 @@ jobs: with: name: coverage-html-report path: htmlcov/ + + # Build the documentation on the release pull request only. + # + # The docs execute every ``.. jupyter-execute::`` cell during the + # build, so they are a second test suite that runs the public API for + # real -- and one that can be broken by a change that touches no + # documentation file at all. Removing Gamma's probability-plot fitting + # did exactly that: a cell looping over the fit methods started + # raising, and nothing noticed, because Read the Docs builds only + # ``master`` and tags. The failure would have appeared as a broken + # hosted build after the release rather than on the change that caused + # it. + # + # ``github.base_ref`` is set only for pull_request events, so this + # runs on the develop -> master release pull request and nowhere else. + # It is deliberately not run on pushes to master: Read the Docs + # rebuilds there anyway, and by that point the gate has nothing left + # to gate. + # + # The point is to reproduce the hosted build, so this matches + # .readthedocs.yaml rather than the test jobs above: Python 3.12, and + # the package installed with its own ``docs`` extra. + docs: + if: github.event_name == 'pull_request' && github.base_ref == 'master' + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + + - name: Set-up Python + uses: actions/setup-python@v5 + with: + python-version: '3.12' + + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install ".[docs]" + + # Not -W: there are pre-existing warnings (duplicate changelog + # labels from autosectionlabel, an rtd-theme deprecation, three + # autodoc imports). Clearing those is worth doing, and then + # failing on warning here and in .readthedocs.yaml together -- + # turning it on before then would just fail every release. + - name: sphinx build + run: python -m sphinx -b html docs docs/_build/html + + - name: Upload documentation artifact + uses: actions/upload-artifact@v4 + with: + name: docs-html + path: docs/_build/html/ diff --git a/docs/Contributing.rst b/docs/Contributing.rst index e2dfdb45..75bb412f 100644 --- a/docs/Contributing.rst +++ b/docs/Contributing.rst @@ -22,10 +22,19 @@ documentation build from running on every change: Continuous integration (``.github/workflows/actions.yml``) therefore runs on **pull requests into develop or master** and on **pushes to master**, rather -than on every push to every branch. Read the Docs is configured to build -``master`` and tags only. The net effect is that the full test suite and the -documentation build run once per pull request and once per release, instead of -once per intermediate commit. +than on every push to every branch. The lint and test jobs run on both; the +documentation build runs on the release pull request into ``master`` only, +where it reproduces the hosted build. Read the Docs itself is configured to +build ``master`` and tags only. + +The documentation build is gated at the release rather than on every pull +request because it executes every code cell in the documentation, which takes +around three minutes from cold rather than the seconds a lint job costs. The +trade-off is that a change which breaks a +documentation example is caught when the release is prepared rather than when +it is merged into ``develop`` -- so if you change the behaviour of a public +function, it is worth building the docs locally before opening the pull +request. Documentation ------------- diff --git a/docs/changelog.rst b/docs/changelog.rst index c928d521..0d85cbbd 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -4,6 +4,36 @@ Changelog v0.19.1 (unreleased) -------------------- +- **The documentation build runs in CI on the release pull request.** + The docs execute every ``.. jupyter-execute::`` cell as they build, so + they are a second test suite that exercises the public API for real -- + and one that a change touching no documentation file at all can break, + as the Gamma entry below did. Read the Docs builds only ``master`` and + tags, so until now that break would have surfaced as a failed hosted + build *after* a release. + + The new job is conditioned on ``github.base_ref == 'master'``, which + is set only for pull requests, so it runs on the ``develop`` -> + ``master`` release pull request and nowhere else. It is not run on + pushes to ``master`` either: Read the Docs rebuilds there anyway, and + by then the gate has nothing left to gate. It matches + ``.readthedocs.yaml`` rather than the test jobs -- Python 3.12, the + package installed via its own ``docs`` extra -- because its purpose is + to reproduce the hosted build, and it uploads the rendered HTML as an + artifact. + + It does not build with ``-W``; the build currently emits 18 warnings, + mostly duplicate labels from ``autosectionlabel`` meeting the + changelog's repeated section headings. Clearing those and then failing + on warning here and in ``.readthedocs.yaml`` together is worth doing + separately. + + The residual gap is deliberate: a documentation break introduced on a + pull request into ``develop`` is caught when the release is prepared, + not when it lands. Building on every pull request would cost minutes + on each, and a path filter would not have helped here -- the change + that broke the build was in ``gamma.py``, not under ``docs/``. + - **Fixed a documentation build broken by the Gamma MPP removal.** The offset-threshold section of *Parametric SurPyval Modelling* ran a ``jupyter-execute`` cell looping over From 0aee893522c34840006435dfb277d0b064ef65cc Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 5 Aug 2026 00:08:21 +0000 Subject: [PATCH 12/13] Run the test suite on the release pull request only Pull requests into `develop` now run lint alone -- about a minute, against the nine the suite takes across three interpreters. The suite still runs in full on the release pull request into `master` and on pushes to `master`. The reason is the edit-review loop. With a single maintainer running the suite locally before pushing, the pull-request run was mostly confirming what was already known, while being the slowest part of working on the package. What this gives up is real and is written down rather than glossed: a failure that appears on only one interpreter is now found when the release is prepared, with a release's worth of commits to search rather than one. That is not hypothetical -- the doctest numeric comparison landed green on 3.11 and failed on 3.12 and 3.13, and it was the pull-request run that caught it. Contributing.rst gains a table of which jobs run on which event, the timings that motivate the split, and what to run locally to compensate: the suite across more than one interpreter when touching numerics, and a docs build when changing the behaviour of a public function. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts --- .github/workflows/actions.yml | 24 +++++++++++++++++ docs/Contributing.rst | 49 +++++++++++++++++++++++++---------- docs/changelog.rst | 19 ++++++++++++++ 3 files changed, 79 insertions(+), 13 deletions(-) diff --git a/.github/workflows/actions.yml b/.github/workflows/actions.yml index e2acf853..159bce14 100644 --- a/.github/workflows/actions.yml +++ b/.github/workflows/actions.yml @@ -6,6 +6,11 @@ run-name: SurPyval CI # master (i.e. release merges + tags). Feature work accumulates on `develop` # via PRs; day-to-day pushes to feature branches no longer each trigger the # full suite. See docs/Contributing.rst for the branching / release flow. +# +# The jobs below are not all gated the same way. Lint runs on everything; +# the test suite and the documentation build run only on the release pull +# request into `master`, where the time they cost is worth paying. Each +# job carries the reasoning for its own condition. on: pull_request: branches: [develop, master] @@ -40,7 +45,26 @@ jobs: - name: black run: black --check $SRC + # The test suite runs on the release pull request into `master`, and on + # pushes to `master`, but not on pull requests into `develop`. + # + # It is roughly nine minutes across the three interpreters, against + # about one for lint, and paying that on every feature pull request + # made the edit-review loop the slowest part of working on the + # package. Development is single-maintainer and the suite is run + # locally before pushing, so the pull-request run was mostly + # confirming what was already known. + # + # What that gives up is real, and worth naming: a failure that only + # appears on one interpreter is now found when the release is + # prepared, with a release's worth of commits to search rather than + # one. That is not hypothetical -- the doctest comparison landed + # green on 3.11 and broke on 3.12 and 3.13, and it was the + # pull-request run that caught it. Run the suite locally across more + # than one interpreter when touching numerics, or open the release + # pull request early and let it sit. surpyval_ci: + if: github.event_name != 'pull_request' || github.base_ref == 'master' runs-on: ubuntu-latest strategy: fail-fast: false diff --git a/docs/Contributing.rst b/docs/Contributing.rst index 75bb412f..0da94ffc 100644 --- a/docs/Contributing.rst +++ b/docs/Contributing.rst @@ -22,19 +22,42 @@ documentation build from running on every change: Continuous integration (``.github/workflows/actions.yml``) therefore runs on **pull requests into develop or master** and on **pushes to master**, rather -than on every push to every branch. The lint and test jobs run on both; the -documentation build runs on the release pull request into ``master`` only, -where it reproduces the hosted build. Read the Docs itself is configured to -build ``master`` and tags only. - -The documentation build is gated at the release rather than on every pull -request because it executes every code cell in the documentation, which takes -around three minutes from cold rather than the seconds a lint job costs. The -trade-off is that a change which breaks a -documentation example is caught when the release is prepared rather than when -it is merged into ``develop`` -- so if you change the behaviour of a public -function, it is worth building the docs locally before opening the pull -request. +than on every push to every branch. Not every job runs on every event: + +.. list-table:: + :header-rows: 1 + :widths: 30 70 + + * - Event + - Jobs + * - Pull request into ``develop`` + - lint only (about a minute) + * - Pull request into ``master`` (the release) + - lint, the test suite across three interpreters, and the + documentation build (about ten minutes) + * - Push to ``master`` / tag + - lint and the test suite; Read the Docs rebuilds the hosted + documentation + +The test suite and the documentation build are both gated at the release +rather than on every pull request because of what they cost: the suite is +roughly nine minutes across the three interpreters and the documentation build +around three from cold, against about one for lint. Paying that on every +feature pull request made the edit-review loop the slowest part of working on +the package, and with a single maintainer running the suite locally before +pushing, the pull-request run was mostly confirming what was already known. + +The trade-off is real and worth understanding before you rely on it. A failure +that appears on only one interpreter, or a change that breaks a documentation +example, is now found when the release pull request is opened -- with a +release's worth of commits to search through rather than one. So: + +* Run the suite locally before pushing, and across more than one interpreter + when you have touched anything numerical. +* Build the documentation locally when you change the behaviour of a public + function, since documentation cells call the real API. +* On a long-running branch, open the release pull request early and let it sit, + so the full run has somewhere to fail before the release itself. Documentation ------------- diff --git a/docs/changelog.rst b/docs/changelog.rst index 0d85cbbd..4ba1cf87 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -4,6 +4,25 @@ Changelog v0.19.1 (unreleased) -------------------- +- **The test suite runs on the release pull request, not on every one.** + Pull requests into ``develop`` now run lint only, about a minute + against the nine the suite takes across three interpreters. The suite + still runs in full on the release pull request into ``master`` and on + pushes to ``master``. + + The reason is the edit-review loop: with a single maintainer running + the suite locally before pushing, the pull-request run was mostly + confirming what was already known, and it was the slowest part of + working on the package. + + What this gives up is stated rather than glossed: a failure that + appears on only one interpreter is now found when the release is + prepared, with a release's worth of commits to search rather than one. + That is not hypothetical -- the doctest numeric comparison two entries + below landed green on 3.11 and failed on 3.12 and 3.13, and it was the + pull-request run that caught it. ``Contributing.rst`` now says which + jobs run on which event, and what to run locally to compensate. + - **The documentation build runs in CI on the release pull request.** The docs execute every ``.. jupyter-execute::`` cell as they build, so they are a second test suite that exercises the public API for real -- From e8594da7de98425cb2b911e1b498ff2245494af6 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 5 Aug 2026 00:45:46 +0000 Subject: [PATCH 13/13] Add a local multi-interpreter check script With the test suite no longer running on pull requests into `develop`, this is the other half of that trade: one command runs the suite and both doctest passes on 3.11, 3.12 and 3.13, and refuses to say "passed" unless all of them did. The failure it guards against is not hypothetical. The doctest numeric comparison passed on 3.11 -- the interpreter it was written on -- and failed on 3.12 and 3.13, because an optimiser landed on a different last digit. Nothing short of running the other interpreters finds that. Environments live in a git-ignored .venvs/ and are reused, so only the first run pays for the installs. uv is used when available and it falls back to venv and pip when not; an interpreter that is not installed is reported rather than fatal. Lint is deliberately absent -- it runs on every pull request already, so it is not what this is for. The command list is a copy of the workflow's, with a comment saying so: if the two drift this stops being a preview of CI and becomes its own thing that can pass while CI fails. Contributing.rst names it as the compensation for the CI split rather than leaving "run it locally across interpreters" as advice with no mechanism. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts --- .gitignore | 3 + docs/Contributing.rst | 15 ++- docs/changelog.rst | 14 +++ scripts/check_all_pythons.py | 195 +++++++++++++++++++++++++++++++++++ 4 files changed, 226 insertions(+), 1 deletion(-) create mode 100644 scripts/check_all_pythons.py diff --git a/.gitignore b/.gitignore index 9a75e36c..f28a9153 100755 --- a/.gitignore +++ b/.gitignore @@ -24,3 +24,6 @@ __pycache__/ publish.sh + +# Local multi-interpreter check environments (scripts/check_all_pythons.py) +.venvs/ diff --git a/docs/Contributing.rst b/docs/Contributing.rst index 0da94ffc..3fede8b8 100644 --- a/docs/Contributing.rst +++ b/docs/Contributing.rst @@ -53,7 +53,20 @@ example, is now found when the release pull request is opened -- with a release's worth of commits to search through rather than one. So: * Run the suite locally before pushing, and across more than one interpreter - when you have touched anything numerical. + when you have touched anything numerical. ``scripts/check_all_pythons.py`` + does exactly that -- it runs what continuous integration would have run, on + 3.11, 3.12 and 3.13: + + .. code-block:: bash + + python scripts/check_all_pythons.py # all three + python scripts/check_all_pythons.py 3.12 # just one + python scripts/check_all_pythons.py --skip-install # reuse as-is + + It keeps its environments in ``.venvs/`` (git-ignored) and reuses them, so + only the first run pays for the installs. It uses ``uv`` when that is + available and falls back to ``venv`` and ``pip`` when it is not. + * Build the documentation locally when you change the behaviour of a public function, since documentation cells call the real API. * On a long-running branch, open the release pull request early and let it sit, diff --git a/docs/changelog.rst b/docs/changelog.rst index 4ba1cf87..468aeb9e 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -4,6 +4,20 @@ Changelog v0.19.1 (unreleased) -------------------- +- **``scripts/check_all_pythons.py`` runs the CI checks locally on every + supported interpreter.** With the suite no longer running on pull + requests into ``develop`` (below), this is the other half of the + trade: one command runs the test suite and both doctest passes on + 3.11, 3.12 and 3.13, and refuses to say "passed" unless all of them + did. + + It keeps its environments in a git-ignored ``.venvs/`` and reuses + them, so only the first run pays for the installs; it uses ``uv`` + when available and falls back to ``venv`` and ``pip`` when not, and + reports an interpreter that is not installed rather than failing on + it. The command list is deliberately a copy of the workflow's, so + what it runs is what CI would have run. + - **The test suite runs on the release pull request, not on every one.** Pull requests into ``develop`` now run lint only, about a minute against the nine the suite takes across three interpreters. The suite diff --git a/scripts/check_all_pythons.py b/scripts/check_all_pythons.py new file mode 100644 index 00000000..5f052f14 --- /dev/null +++ b/scripts/check_all_pythons.py @@ -0,0 +1,195 @@ +#!/usr/bin/env python3 +"""Run the deployment checks locally, across every supported Python. + +Continuous integration runs the test suite only on the release pull +request into ``master`` (see docs/Contributing.rst). That keeps the +edit-review loop fast, at the cost of finding an interpreter-specific +failure at release time rather than when it lands. This script is the +other half of that trade: it runs what CI would have run, on all three +interpreters, before you push. + +It exists because the failure it guards against is not hypothetical. +The doctest numeric comparison added in 0.19.1 passed on 3.11 -- the +interpreter it was written on -- and failed on 3.12 and 3.13, because an +optimiser landed on a different last digit. Nothing short of actually +running the other interpreters would have found it. + + python scripts/check_all_pythons.py # 3.11, 3.12, 3.13 + python scripts/check_all_pythons.py 3.12 # just one + python scripts/check_all_pythons.py --skip-install # reuse as-is + +Environments live in ``.venvs/py3.X`` (git-ignored) and are reused +between runs, so only the first is slow. ``uv`` is used when it is +installed, because it makes the install step take seconds; otherwise +this falls back to ``venv`` and ``pip``, which works the same and takes +longer. + +Exits non-zero if any check on any interpreter fails. +""" + +from __future__ import annotations + +import argparse +import shutil +import subprocess +import sys +import time +from pathlib import Path + +ROOT = Path(__file__).resolve().parent.parent +VENVS = ROOT / ".venvs" +DEFAULT_VERSIONS = ["3.11", "3.12", "3.13"] + +# What CI runs, minus lint (which runs on every pull request anyway, so +# it is not what this script is for). Kept in the same order and with +# the same arguments as .github/workflows/actions.yml -- if they drift, +# this stops being a preview of CI and becomes its own thing. +CHECKS: list[tuple[str, list[str]]] = [ + ( + "test suite", + [ + "-m", + "pytest", + "-n", + "auto", + "-q", + "--ignore=surpyval/tests/alpha", + "--run-ml", + ], + ), + ( + "doctests", + [ + "-m", + "pytest", + "--doctest-modules", + "surpyval", + "-q", + "--ignore=surpyval/tests", + "--ignore=surpyval/alpha", + ], + ), + ( + "doctests (numeric forced)", + [ + "-m", + "pytest", + "--doctest-modules", + "surpyval", + "-q", + "--ignore=surpyval/tests", + "--ignore=surpyval/alpha", + "--doctest-force-numeric", + ], + ), +] + + +def run(cmd: list[str], **kwargs) -> int: + """Run a command, streaming its output, and return its exit status.""" + return subprocess.call(cmd, cwd=ROOT, **kwargs) + + +def interpreter_for(venv: Path) -> Path: + bin_dir = "Scripts" if sys.platform == "win32" else "bin" + exe = "python.exe" if sys.platform == "win32" else "python" + return venv / bin_dir / exe + + +def ensure_env(version: str, skip_install: bool) -> Path | None: + """Create (or reuse) the environment for ``version``. + + Returns the interpreter path, or None if the interpreter is not + available on this machine -- which is a thing to report, not a + thing to crash on. + """ + venv = VENVS / f"py{version}" + python = interpreter_for(venv) + uv = shutil.which("uv") + + if not python.exists(): + print(f"\n[{version}] creating {venv.relative_to(ROOT)}") + if uv: + created = run([uv, "venv", "--python", version, str(venv)]) + else: + base = shutil.which(f"python{version}") + if base is None: + print(f"[{version}] no python{version} on PATH -- skipping") + return None + created = run([base, "-m", "venv", str(venv)]) + if created != 0 or not python.exists(): + print(f"[{version}] could not create an environment -- skipping") + return None + + if not skip_install: + print(f"[{version}] installing") + # pytest-xdist is not in the tests extra; it is a development + # tool, and the suite's ``-n auto`` needs it. + target = ["-e", ".[tests]", "pytest-xdist"] + if uv: + installed = run( + [uv, "pip", "install", "-q", "--python", str(python), *target] + ) + else: + installed = run( + [str(python), "-m", "pip", "install", "-q", *target] + ) + if installed != 0: + print(f"[{version}] install failed -- skipping") + return None + + return python + + +def main() -> int: + parser = argparse.ArgumentParser( + description="Run the CI checks locally on every supported Python." + ) + parser.add_argument( + "versions", + nargs="*", + default=DEFAULT_VERSIONS, + help="Python versions to check (default: %s)" + % " ".join(DEFAULT_VERSIONS), + ) + parser.add_argument( + "--skip-install", + action="store_true", + help="reuse the environments as they are, without reinstalling", + ) + args = parser.parse_args() + + VENVS.mkdir(exist_ok=True) + results: list[tuple[str, str, str, float]] = [] + + for version in args.versions: + python = ensure_env(version, args.skip_install) + if python is None: + results.append((version, "environment", "UNAVAILABLE", 0.0)) + continue + for name, cmd in CHECKS: + print(f"\n[{version}] {name}") + started = time.monotonic() + status = run([str(python), *cmd]) + elapsed = time.monotonic() - started + results.append( + (version, name, "ok" if status == 0 else "FAILED", elapsed) + ) + + print("\n" + "=" * 60) + failed = False + for version, name, status, elapsed in results: + if status != "ok": + failed = True + print(f"{version:6s} {name:28s} {status:12s} {elapsed:6.1f}s") + print("=" * 60) + + if failed: + print("\nSomething failed. Do not push.") + return 1 + print("\nAll checks passed on every interpreter.") + return 0 + + +if __name__ == "__main__": + sys.exit(main())