fix(#124): make GLS convergence invariant to observable amplitude

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

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

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
 

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):

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: