diff --git a/CHANGELOG.md b/CHANGELOG.md index fc8d85e..6d40ff1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -63,6 +63,44 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). the change is not visible in `input_hash`. ### Fixed +- **GLS optimiser residuals are divided by the curve's own scale, and fall + back to the raw residuals where that rescaling is not representable (issue + #124, PR #134).** Rescaling makes SciPy's termination amplitude-invariant + for a model whose amplitude is not a fitted parameter: `1e-40*exp(-1.3 t)` + no longer returns its seed rate as a converged fit. It does not cure a + FREE amplitude parameter, which still returns the seed rate at 1e-40 when + fitted from `p0 = [1e-40, 0.2]` without bounds. + SciPy's finite-difference probes, however, step in absolute parameter + units, so for a tiny scale the rescaled residual, Jacobian or cost can leave + float64 -- measured at `max|y| = 5e-324` (the #147 fixture) and, through + `_fit_with_model("M0")`, at 1e-150 and 1e-310, where the fit was reported + unsuccessful and dropped from AICc selection. A floating-point exception + raised by the rescaled arithmetic itself (the division, the whitening, or + SciPy's own computation on the rescaled values -- NOT the evaluation of the + model, which keeps the caller's floating-point policy) or a non-finite + cost/residual/Jacobian of the rescaled solve now sends that iteration to the + raw residuals, evaluated exactly as before #124, with a `RuntimeWarning` + stating that the #124 invariance does not hold for the curve. A finite + rescaled solve that did not converge, and an exception with no + floating-point event, are still reported unsuccessful. For a FREE amplitude + parameter at `max|y| >= ~1e150` the fallback fires as well (measured at + 1e150, 1e200 and 1e300; a fixed-amplitude model does not fall back there and + fits rate 1.3) and the fit fails closed exactly as on main (`success=False`), + without leaking the private exception that an earlier revision of this fix + raised there (PR #134 review). Unlike main, such a fit now emits the fallback + warning. It has the dedicated category + `liouscope.fitting.gls.AmplitudeRescalingFallbackWarning`, a `RuntimeWarning` + subclass, so a caller running with `-W error` can filter it by class. Raw residuals, rho, sigma and + the likelihood remain in the caller's data units. +- **The AR(1) lag-1 autocorrelation is formed from residuals normalised by an + exact power of two (PR #134, round-3 review).** `ar1_correlation` took its dot + products on unscaled residuals, which underflow to 0 below ~1e-162 -- the + corrected rho then collapses to its floor `1/(n-3)` -- and overflow to NaN + above ~1e154. Measured through `fit_gls_ar1` with `n_iters=3` on one AR(1) + correlated curve: rho 0.012987 at 1e-170, 0.44637 at 1e0 and NaN with + `success=False` at 1e160. Scaling by `2**k` is exact, so rho is bit-identical + wherever the old computation stayed in range (200/200 random series between + 1e-100 and 1e100). - **The trace-preservation defect overflowed on the way to a column sum that is exactly zero (issue #139, P2 review).** `trace_preservation_defect` assembled `vec(I)^H L` with an ordinary matrix product, which accumulates in ordinary diff --git a/CITATION.cff b/CITATION.cff index eecf07a..8d0e720 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -140,6 +140,22 @@ abstract: >- # MANIFEST_SCHEMA is deliberately unchanged -- `trace_defect` is a report # field, not a run-manifest field, so the manifest contract is untouched and # `input_hash` does not move. +# Also pending for the next cut: the GLS AR(1) fit presents the optimiser with +# residuals divided by the curve's own scale max|y| (issue #124), so SciPy's +# termination no longer turns the seed rate into a "converged" estimate for a +# tiny-amplitude curve whose amplitude is not a fitted parameter; where that +# rescaled problem leaves float64, the iteration falls back, with a +# RuntimeWarning, to the unscaled residuals exactly as before (PR #134). The +# lag-1 autocorrelation behind the AR(1) whitening is formed from residuals +# normalised by an exact power of two, so it no longer underflows to its floor +# below ~1e-162 or overflows to NaN above ~1e154; where the old computation +# stayed in range it is bit-identical. This CHANGES NUMERICAL RESULTS and is a +# METHODOLOGY correction: fitted rates, rho, AICc and the selected model can +# differ for curves far from unit scale. It must NOT be described as making the +# fit amplitude-invariant in general: a FREE amplitude parameter still returns +# its seed rate at 1e-40, and the fallback keeps the unscaled problem's own +# limits. MANIFEST_SCHEMA is deliberately unchanged -- none of these quantities +# is a run-manifest field, so `input_hash` does not move. keywords: - open quantum systems - Lindblad diff --git a/src/liouscope/fitting/bootstrap.py b/src/liouscope/fitting/bootstrap.py index 1ec84e1..548c1d8 100644 --- a/src/liouscope/fitting/bootstrap.py +++ b/src/liouscope/fitting/bootstrap.py @@ -13,6 +13,7 @@ from __future__ import annotations +import warnings from collections.abc import Callable import numpy as np @@ -185,6 +186,14 @@ def bca_ci( Returns ``(p, 2)`` array of ``(lo, hi)`` per parameter at level ``1 - alpha``. + + A parameter whose bootstrap distribution is degenerate, or whose adjusted + quantiles collapse to an exactly zero-width interval, is reported as + ``(NaN, NaN)`` with a warning. This follows the fail-closed semantics used + by SciPy's BCa implementation for degenerate bootstrap distributions: + arithmetic concentration is not promoted to a claim of perfect certainty + when the resampling model contains no measurable spread (issue #125). + Other non-degenerate parameters in the same fit remain reportable. """ samples = np.asarray(samples, dtype=float) theta_hat = np.asarray(theta_hat, dtype=float) @@ -194,6 +203,16 @@ def bca_ci( z_alpha_hi = ndtri(1.0 - alpha / 2.0) for j in range(p): boot_j = np.sort(samples[:, j]) + if boot_j.size == 0 or np.all(boot_j == boot_j[0]): + warnings.warn( + f"bca_ci: bootstrap distribution for parameter {j} is " + "degenerate; uncertainty is unavailable rather than a " + "zero-width confidence claim (issue #125)", + RuntimeWarning, + stacklevel=2, + ) + cis[j] = (float("nan"), float("nan")) + continue # Bias-correction z0 with half-correction for ties (Efron 1987, S3 # audit 2026-06-04). Using a strict ``<`` only would send the # proportion to 0 (hence z0 -> -inf, clamped) whenever many bootstrap @@ -231,6 +250,16 @@ def adjusted(z_alpha: float, z0_: float = z0, a_: float = a) -> float: # continuous-percentile estimator. lo = float(np.quantile(boot_j, q_lo, method="linear")) hi = float(np.quantile(boot_j, q_hi, method="linear")) + if lo == hi: + warnings.warn( + f"bca_ci: confidence interval for parameter {j} collapsed to " + "exactly zero width; uncertainty is unavailable rather than " + "perfectly known (issue #125)", + RuntimeWarning, + stacklevel=2, + ) + cis[j] = (float("nan"), float("nan")) + continue cis[j] = (lo, hi) return cis diff --git a/src/liouscope/fitting/gls.py b/src/liouscope/fitting/gls.py index 135419b..38aedbc 100644 --- a/src/liouscope/fitting/gls.py +++ b/src/liouscope/fitting/gls.py @@ -16,9 +16,10 @@ import warnings from collections.abc import Callable from dataclasses import dataclass +from typing import Any import numpy as np -from scipy.optimize import least_squares +from scipy.optimize import OptimizeResult, least_squares from ..numerics.norms import scaled_log_sum_squares from .aicc import gaussian_log_likelihood @@ -58,6 +59,28 @@ class GLSFitOutput: scale_unavailable: bool = False +class _ScaledResidualOverflowError(ArithmeticError): + """A residual rescaled for the optimiser left float64 while the raw one is finite. + + Private control flow of :func:`fit_gls_ar1` (PR #134). Deliberately NOT a + ``ValueError``/``RuntimeError``: those mean "the fit failed" there, while + this means only "the #124 rescaling is not representable for this curve". + """ + + +class AmplitudeRescalingFallbackWarning(RuntimeWarning): + """The #124 residual rescaling was not representable; the fit ran unscaled. + + A dedicated ``RuntimeWarning`` subclass (PR #134 review): for a free + amplitude parameter at ``max|y| >= ~1e150`` the fallback fires where the + pre-#124 code was silent, with a result identical to it. A caller running + with ``-W error`` can filter exactly this notice by class instead of all + RuntimeWarnings. It is emitted per fit, not once per process: which curve + lost the #124 invariance is audit information, and a process-global + "once" registry would make it depend on call order. + """ + + def _whiten(y: np.ndarray, rho: float) -> np.ndarray: if y.size < 2: y_copy: np.ndarray = y.copy() @@ -181,19 +204,130 @@ def fit_gls_ar1( degenerate=True, ) + # Issue #124: the mathematical least-squares optimum is unchanged when an + # observable is multiplied by a positive constant, but SciPy's numerical + # termination is not. In particular, TRF declares success when its scaled + # gradient falls below ``gtol``; a curve such as ``1e-40*exp(-1.3*t)`` can + # therefore return the INITIAL RATE as a "converged" estimate after one + # evaluation. Divide only the residuals presented to the optimiser by the + # curve's own finite, non-zero scale. This multiplies the objective by one + # positive constant, so the minimiser is identical, while the numerical + # problem becomes amplitude-scale invariant. Raw residuals, AR(1) rho, + # sigma and likelihood below remain in the caller's original data units. + # + # The rescaled problem is not always REPRESENTABLE (PR #134 x PR #147). + # SciPy's finite-difference probes step by ``~1.5e-8 * max(1, |x|)`` in + # ABSOLUTE parameter units, so the model moves by an amount unrelated to + # the curve's scale. Measured on the #147 fixture (``max|y| = 5e-324``): + # the scaled residual at ``p0`` is finite (``max = 1.0``), the first + # Jacobian probe gives ``1.5e-8 / 5e-324 = inf``, ``least_squares`` raises + # and the fit was reported unsuccessful -- a curve withheld from model + # selection because of its absolute amplitude, the boundary #135 removed. + # The overflow is a property of the rescaling, not of the fit: when a + # scaled residual leaves float64 while the raw one is finite, the + # iteration is repeated on the raw residuals, exactly the pre-#124 + # problem, and the loss of the #124 invariance is announced rather than + # silent. A raw residual that is itself non-finite still fails closed. + fit_scale = y_scale + + def _raw_residual_fn(rho_local: float) -> Callable[[np.ndarray], np.ndarray]: + # The unscaled problem, verbatim as before #124: no detection, no + # errstate. A non-finite raw residual reaches ``least_squares``, which + # raises ``ValueError`` and the fit fails closed (PR #134 review: the + # detector used to run here too and leaked a private exception from + # the public API for ``max|y| >= ~1e155``). + def residual(params: np.ndarray) -> np.ndarray: + return _whiten(y - model(t, params), rho_local) + + return residual + + def _rescaled_residual_fn( + rho_local: float, + scale: float, + caller_err: dict[str, Any], + caller_call: Any, + ) -> Callable[[np.ndarray], np.ndarray]: + # Runs inside the counting errstate of ``_solve_rescaled``. The MODEL + # and the raw difference are evaluated under the CALLER's floating-point + # policy instead: an event there (e.g. a sinc factor's 0/0 at t = 0, + # finite after ``np.where``) belongs to the raw problem, not to the + # rescaling, and must not trigger the fallback (PR #134 review: it did, + # and returned the seed rate as a converged fit at 1e-40). + def residual(params: np.ndarray) -> np.ndarray: + with np.errstate(call=caller_call, **caller_err): + r_raw = y - model(t, params) + # Counted: an overflow of the division or the whitening IS the + # rescaled problem leaving float64. + return _whiten(r_raw / scale, rho_local) + + return residual + + # The residual is not the only place the rescaled problem can leave + # float64: with a FREE amplitude parameter the Jacobian and the cost carry + # the same ``1/scale`` factor. Measured through ``_fit_with_model("M0")`` + # on ``s*exp(-1.3 t)``: at s = 1e-150 and 1e-310 the rescaled optimiser + # overflowed inside SciPy (``dot``/``square``) and the fit was reported + # unsuccessful (aicc = inf), where the raw problem fits rate 1.3. So every + # floating-point exception raised while the RESCALED problem is solved, + # and any non-finite cost/residual/Jacobian it returns, counts as "not + # representable" and sends the iteration to the raw residuals. A finite + # rescaled result that merely did not converge is NOT retried: that is a + # genuine failure, and retrying it raw would re-admit the #124 seed. + def _solve_rescaled( + rho_local: float, x0: np.ndarray, **kw: object + ) -> OptimizeResult: + events: list[str] = [] + + def _record(kind: str, _flag: int) -> None: + events.append(kind) + + caller_err = dict(np.geterr()) + caller_call = np.geterrcall() + with np.errstate( + over="call", divide="call", invalid="call", under="ignore", call=_record + ): + try: + res = least_squares( + _rescaled_residual_fn(rho_local, fit_scale, caller_err, caller_call), + x0, + **kw, + ) + except (ValueError, RuntimeError): + if events: + raise _ScaledResidualOverflowError from None + raise + if events or not ( + np.isfinite(res.cost) + and np.all(np.isfinite(res.fun)) + and np.all(np.isfinite(res.jac)) + ): + raise _ScaledResidualOverflowError + return res + rho = 0.0 success = True for _ in range(n_iters): - def residual(params: np.ndarray, rho_local: float = rho) -> np.ndarray: - y_hat = model(t, params) - r = y - y_hat - return _whiten(r, rho_local) - ls_kwargs: dict[str, object] = {"max_nfev": max_nfev} if bounds is not None: ls_kwargs["bounds"] = bounds try: - result = least_squares(residual, p, **ls_kwargs) + try: + if fit_scale == 1.0: + result = least_squares(_raw_residual_fn(rho), p, **ls_kwargs) + else: + result = _solve_rescaled(rho, p, **ls_kwargs) + except _ScaledResidualOverflowError: + warnings.warn( + "fit_gls_ar1: residuals rescaled by the curve's own scale " + f"({fit_scale:.3e}) are not representable as float64 at the " + "optimiser's probe points; fitting the unscaled residuals " + "instead, so the amplitude-scale invariance of issue #124 " + "does not hold for this curve (PR #134).", + AmplitudeRescalingFallbackWarning, + stacklevel=2, + ) + fit_scale = 1.0 + result = least_squares(_raw_residual_fn(rho), p, **ls_kwargs) p = result.x success = result.success except (ValueError, RuntimeError): diff --git a/src/liouscope/fitting/neff.py b/src/liouscope/fitting/neff.py index b888b59..5b932f2 100644 --- a/src/liouscope/fitting/neff.py +++ b/src/liouscope/fitting/neff.py @@ -15,6 +15,7 @@ from __future__ import annotations +import math import warnings import numpy as np @@ -100,9 +101,23 @@ def ar1_correlation(residuals: np.ndarray) -> float: coefficient is estimator-convention dependent, so treat it as indicative rather than a pinned analytic identity. """ - x = np.asarray(residuals, dtype=float) - float(np.mean(residuals)) + x = np.asarray(residuals, dtype=float) if x.size < 2: return 0.0 + # rho is a RATIO of two quadratic forms, so it is invariant under a common + # rescaling of the residuals -- but the dot products that form it are not: + # residuals below ~1e-162 square to 0 (rho_hat = 0, and the corrected + # estimate collapses to its floor 1/(n-3)), and above ~1e154 they overflow + # to NaN (PR #134 review, measured through fit_gls_ar1 at 1e-170 / 1e160). + # Normalise by the exact power of two nearest the peak magnitude. Scaling + # by 2**k is exact in binary floating point whenever nothing under- or + # overflows, so every mean, product and sum below is the old one times an + # exact power of two and rho_hat is BIT-IDENTICAL wherever the unnormalised + # computation stayed in range -- the audit-pinned formula does not move. + peak = float(np.max(np.abs(x))) + if math.isfinite(peak) and peak > 0.0: + x = np.ldexp(x, -math.frexp(peak)[1]) + x = x - float(np.mean(x)) num = float(np.dot(x[:-1], x[1:])) den = float(np.dot(x, x)) if den == 0.0: diff --git a/tests/test_gls_amplitude_scale.py b/tests/test_gls_amplitude_scale.py new file mode 100644 index 0000000..555ea3b --- /dev/null +++ b/tests/test_gls_amplitude_scale.py @@ -0,0 +1,309 @@ +"""Regression tests for observable-amplitude invariance in the GLS optimiser. + +Issue #124: the mathematical fit is unchanged by ``y -> c*y`` for ``c > 0``, +but feeding the raw residuals to SciPy lets the gradient termination criterion +accept the initial seed when ``c`` is tiny. These tests pin both directions: +the tiny-amplitude counterexample must fit, and the ordinary-scale positive +control must remain unchanged. +""" + +from __future__ import annotations + +import numpy as np +import pytest + +from liouscope.fitting.gls import fit_gls_ar1 + +_TRUE_RATE = 1.3 +_SEED_RATE = 0.2 +_T = np.linspace(0.0, 5.0, 64) + + +def _fit_exponential(scale: float): + def model(t: np.ndarray, params: np.ndarray) -> np.ndarray: + return scale * np.exp(-params[0] * t) + + y = scale * np.exp(-_TRUE_RATE * _T) + return fit_gls_ar1( + model, + _T, + y, + np.array([_SEED_RATE]), + bounds=(np.array([0.0]), np.array([5.0])), + n_iters=1, + ) + + +def test_tiny_amplitude_cannot_turn_the_seed_into_a_measurement() -> None: + """The exact #124 counterexample must move away from its initial seed.""" + tiny = _fit_exponential(1.0e-40) + + assert tiny.success + assert not tiny.degenerate + assert tiny.params[0] == pytest.approx(_TRUE_RATE, rel=1.0e-7, abs=1.0e-10) + assert abs(tiny.params[0] - _SEED_RATE) > 0.5 + + +def test_amplitude_rescaling_preserves_the_fitted_rate() -> None: + """Positive control: ordinary and tiny amplitudes represent the same fit.""" + ordinary = _fit_exponential(1.0) + tiny = _fit_exponential(1.0e-40) + + assert ordinary.success and tiny.success + assert ordinary.params[0] == pytest.approx(_TRUE_RATE, rel=1.0e-7, abs=1.0e-10) + assert tiny.params[0] == pytest.approx(ordinary.params[0], rel=1.0e-9, abs=1.0e-12) + + +@pytest.mark.parametrize("scale", [1.0e-150, 1.0e-310]) +def test_unrepresentable_rescaling_falls_back_instead_of_dropping_the_fit( + scale: float, +) -> None: + """PR #134: the #124 rescaling must not withhold a fit the raw problem makes. + + With a FREE amplitude parameter the rescaled Jacobian and cost carry the + same ``1/scale`` factor as the residual. Measured on the merge of #124 into + main: ``_fit_with_model("M0")`` overflowed inside SciPy at these scales and + returned ``success=False`` / ``aicc=inf``, while main fitted rate 1.3. The + Jacobian-level overflow is the case the residual-level check cannot see, so + it is pinned here separately from the #147 subnormal fixture. + """ + import warnings + + from liouscope.diagnostics import relaxation as relaxation_mod + + y = scale * np.exp(-_TRUE_RATE * _T) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + fit_result, _ = relaxation_mod._fit_with_model("M0", _T, y) + assert any( + "not representable as float64 at the" in str(w.message) for w in caught + ), [str(w.message) for w in caught] + assert fit_result.success + assert np.isfinite(fit_result.aicc) + assert fit_result.params[1] == pytest.approx(_TRUE_RATE, rel=1.0e-7) + + +# -------------------------------------------------------------------------- +# PR #134 review (Equalita NO-GO for 7ad4829): one discriminating test per +# part of the representability detector. +# -------------------------------------------------------------------------- + +_FALLBACK = "not representable as float64 at the" +_NOISE = np.random.default_rng(20260911).standard_normal(_T.size) + + +def _free_amplitude(t: np.ndarray, params: np.ndarray) -> np.ndarray: + return params[0] * np.exp(-params[1] * t) + + +def _recorded(fn, *args, **kwargs): # type: ignore[no-untyped-def] + import warnings + + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + out = fn(*args, **kwargs) + return out, [str(w.message) for w in caught] + + +@pytest.mark.parametrize("scale", [1.0e200, 1.0e300]) +def test_large_amplitude_fails_closed_without_leaking_the_private_exception( + scale: float, +) -> None: + """The raw path is main's problem verbatim: ``success=False``, no raise. + + Measured on 7ad4829: the residual detector also ran on the RAW residuals, + whose whitening overflows here although they are finite, and its private + ``ArithmeticError`` escaped both public entry points (20/20 cases). + """ + from liouscope.diagnostics import relaxation as relaxation_mod + + y = scale * np.exp(-_TRUE_RATE * _T) + fit, _ = _recorded(fit_gls_ar1, _free_amplitude, _T, y, np.array([scale, 0.2])) + assert not fit.success + (fit_result, _), _ = _recorded( + relaxation_mod._fit_with_model, "M0", _T, y * (1.0 + 1.0e-3 * _NOISE) + ) + assert not fit_result.success + assert np.isinf(fit_result.aicc) + + +def test_benign_model_float_event_does_not_trigger_the_fallback() -> None: + """A 0/0 inside the MODEL belongs to the raw problem, not to the rescaling. + + ``sin(t)/t`` at ``t = 0`` raises 'invalid' and is replaced by ``np.where``; + the model output is finite. When that event was counted, every fit of this + model fell back and returned the seed rate 0.2 as converged at 1e-40. + """ + + def run(scale: float) -> tuple[float, list[str]]: + def model(t: np.ndarray, params: np.ndarray) -> np.ndarray: + return scale * np.exp(-params[0] * t) * np.where(t > 0, np.sin(t) / t, 1.0) + + with np.errstate(invalid="ignore"): # building the data is not under test + y = model(_T, np.array([_TRUE_RATE])) * (1.0 + 1.0e-3 * _NOISE) + fit, messages = _recorded( + fit_gls_ar1, model, _T, y, np.array([_SEED_RATE]), + bounds=(np.array([0.0]), np.array([5.0])), n_iters=1, + ) + assert fit.success + return float(fit.params[0]), messages + + ordinary, _ = run(2.0) + tiny, messages = run(1.0e-40) + assert not any(_FALLBACK in m for m in messages), messages + assert tiny == pytest.approx(ordinary, rel=1.0e-9) + assert abs(tiny - _SEED_RATE) > 0.5 + + +def test_nonfinite_rescaled_result_without_float_events_falls_back( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """The non-finite arm: a NaN Jacobian set by assignment raises no FP event.""" + from scipy.optimize import least_squares as real + + + calls = {"n": 0} + + def first_call_poisoned(*args, **kwargs): # type: ignore[no-untyped-def] + res = real(*args, **kwargs) + calls["n"] += 1 + if calls["n"] == 1: + res.jac = np.full_like(res.jac, np.nan) + return res + + monkeypatch.setattr("liouscope.fitting.gls.least_squares", first_call_poisoned) + y = 1.0e-40 * np.exp(-_TRUE_RATE * _T) * (1.0 + 1.0e-3 * _NOISE) + fit, messages = _recorded( + fit_gls_ar1, lambda t, p: 1.0e-40 * np.exp(-p[0] * t), _T, y, + np.array([_SEED_RATE]), bounds=(np.array([0.0]), np.array([5.0])), n_iters=1, + ) + assert any(_FALLBACK in m for m in messages), messages + assert calls["n"] == 2 + + +def test_exception_without_float_events_is_not_rerouted() -> None: + """A model error in the rescaled solve is a failed fit, not a scale problem.""" + + def failing(t: np.ndarray, params: np.ndarray) -> np.ndarray: + if params[0] != _SEED_RATE: + raise RuntimeError("model refuses this rate") + return 1.0e-40 * np.exp(-params[0] * t) + + y = 1.0e-40 * np.exp(-_TRUE_RATE * _T) * (1.0 + 1.0e-3 * _NOISE) + fit, messages = _recorded(fit_gls_ar1, failing, _T, y, np.array([_SEED_RATE]), n_iters=1) + assert not fit.success + assert not any(_FALLBACK in m for m in messages), messages + + +def test_finite_nonconverged_rescaled_solve_is_not_retried_raw() -> None: + """Retrying a finite non-converged solve raw would re-admit the #124 seed.""" + y = 1.0e-40 * np.exp(-_TRUE_RATE * _T) * (1.0 + 1.0e-3 * _NOISE) + fit, messages = _recorded( + fit_gls_ar1, lambda t, p: 1.0e-40 * np.exp(-p[0] * t), _T, y, + np.array([_SEED_RATE]), bounds=(np.array([0.0]), np.array([5.0])), + n_iters=1, max_nfev=1, + ) + assert not fit.success + assert not any(_FALLBACK in m for m in messages), messages + + +# -------------------------------------------------------------------------- +# PR #134 round 3 (Codex P1): the AR(1) rho between Cochrane-Orcutt iterations +# was estimated from UNSCALED residuals, whose dot products underflow to 0 +# below ~1e-162 and overflow to NaN above ~1e154. +# -------------------------------------------------------------------------- + + +def _ar1_curve(scale: float) -> tuple[np.ndarray, np.ndarray]: + t = np.linspace(0.0, 5.0, 80) + rng = np.random.default_rng(20260911) + e = np.zeros(t.size) + nu = rng.standard_normal(t.size) + for i in range(t.size): + e[i] = (0.6 * e[i - 1] if i else 0.0) + nu[i] + return t, scale * (np.exp(-_TRUE_RATE * t) + 1.0e-3 * e) + + +@pytest.mark.parametrize("scale", [1.0e-170, 1.0e160]) +def test_ar1_rho_between_iterations_is_amplitude_invariant(scale: float) -> None: + """Amplitude-equivalent correlated curves give the same rho, rate, success. + + Measured on 9a623e8: rho 0.012987 (= the corrected floor 1/(n-3), i.e. an + underflowed rho_hat of 0) at 1e-170, 0.44637 at 1e0, and NaN with + ``success=False`` at 1e160. + """ + + def fit(s: float): # type: ignore[no-untyped-def] + t, y = _ar1_curve(s) + out, _ = _recorded( + fit_gls_ar1, lambda tt, p: s * p[0] * np.exp(-p[1] * tt), t, y, + np.array([1.0, 0.8]), n_iters=3, + ) + return out + + ref = fit(1.0) + other = fit(scale) + assert ref.success and np.isfinite(ref.rho_ar1) and ref.rho_ar1 > 0.3 + assert other.success == ref.success + assert other.rho_ar1 == pytest.approx(ref.rho_ar1, rel=1.0e-8) + assert other.params[1] == pytest.approx(ref.params[1], rel=1.0e-8) + + +def test_ar1_correlation_is_bit_identical_under_power_of_two_rescaling() -> None: + """The normalisation must not move rho where the old arithmetic was in range.""" + from liouscope.fitting.neff import ar1_correlation + + _, y = _ar1_curve(1.0) + base = ar1_correlation(y) + for k in (-900, -540, 500, 1000): + assert ar1_correlation(np.ldexp(y, k)) == base, k + + +# -------------------------------------------------------------------------- +# PR #134 round 3 (Equalita follow-up on 9a623e8) +# -------------------------------------------------------------------------- + + +def _sinc_model(t: np.ndarray, params: np.ndarray) -> np.ndarray: + # Same values at every parameter; the 0/0 at t = 0 is only RAISED at the + # seed, i.e. inside the optimiser's first evaluation. The post-fit model + # evaluation (at ~1.3) is event-free, so it cannot mask what the solve did. + if params[0] == _SEED_RATE: + sinc = np.where(t > 0, np.sin(t) / t, 1.0) + else: + sinc = np.where(t > 0, np.sin(t) / np.where(t > 0, t, 1.0), 1.0) + return 1.0e-40 * np.exp(-params[0] * t) * sinc + + +def test_model_runs_under_the_callers_floating_point_policy() -> None: + """The detector's errstate must not replace the caller's for the MODEL. + + A caller that asks for ``invalid="raise"`` gets the model's 0/0 during the + solve as a ``FloatingPointError``, exactly as on main. Evaluating the + model under the detector's own errstate (ignore or call) swallows it. + """ + y = _sinc_model(_T, np.array([_TRUE_RATE])) # event-free away from the seed + with np.errstate(invalid="raise"): + with pytest.raises(FloatingPointError): + fit_gls_ar1( + _sinc_model, _T, y, np.array([_SEED_RATE]), + bounds=(np.array([0.0]), np.array([5.0])), n_iters=1, + ) + + +def test_fallback_warning_has_a_dedicated_filterable_category() -> None: + """``-W error`` users can silence exactly the fallback notice by class.""" + import warnings + + from liouscope.fitting.gls import AmplitudeRescalingFallbackWarning + + assert issubclass(AmplitudeRescalingFallbackWarning, RuntimeWarning) + scale = 1.0e200 + y = scale * np.exp(-_TRUE_RATE * _T) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + fit = fit_gls_ar1(_free_amplitude, _T, y, np.array([scale, 0.2])) + fallback = [w for w in caught if _FALLBACK in str(w.message)] + assert fallback, [str(w.message) for w in caught] + assert all(w.category is AmplitudeRescalingFallbackWarning for w in fallback) + assert not fit.success diff --git a/tests/test_issue125_degenerate_bca.py b/tests/test_issue125_degenerate_bca.py new file mode 100644 index 0000000..7e2fa40 --- /dev/null +++ b/tests/test_issue125_degenerate_bca.py @@ -0,0 +1,57 @@ +"""Regression tests for issue #125: degenerate BCa uncertainty claims.""" + +from __future__ import annotations + +import numpy as np +import pytest + +from liouscope.fitting.bootstrap import bca_ci + + +def test_degenerate_bootstrap_distribution_is_unavailable_not_zero_width() -> None: + samples = np.ones((200, 1), dtype=float) + theta_hat = np.array([1.0]) + + with pytest.warns(RuntimeWarning, match="degenerate"): + ci = bca_ci(samples, theta_hat) + + assert np.isnan(ci[0, 0]) + assert np.isnan(ci[0, 1]) + + +def test_degenerate_parameter_does_not_withhold_an_independent_parameter() -> None: + x = np.linspace(-1.0, 1.0, 200) + samples = np.column_stack([np.ones_like(x), x]) + theta_hat = np.array([1.0, 0.0]) + + with pytest.warns(RuntimeWarning, match="parameter 0"): + ci = bca_ci(samples, theta_hat) + + assert np.all(np.isnan(ci[0])) + assert np.all(np.isfinite(ci[1])) + assert ci[1, 0] < ci[1, 1] + + +def test_ordinary_nondegenerate_distribution_keeps_a_finite_interval() -> None: + rng = np.random.default_rng(125) + samples = rng.normal(loc=0.7, scale=0.03, size=(1000, 1)) + theta_hat = np.array([0.7]) + + ci = bca_ci(samples, theta_hat) + + assert np.all(np.isfinite(ci)) + assert ci[0, 0] < ci[0, 1] + assert ci[0, 0] < theta_hat[0] < ci[0, 1] + + +def test_zero_width_after_adjusted_quantiles_is_not_reported_as_certainty() -> None: + # Non-identical distribution, but >97.5% of its mass sits on one value. + # The BCa adjusted endpoints can therefore collapse to that value. The + # claim boundary, not merely the raw distribution, must catch the collapse. + samples = np.concatenate([np.ones(199), np.array([2.0])]).reshape(-1, 1) + theta_hat = np.array([1.0]) + + with pytest.warns(RuntimeWarning, match="zero width"): + ci = bca_ci(samples, theta_hat) + + assert np.all(np.isnan(ci[0]))