fix(#130): make TP evidence norms underflow/overflow-safe

CHANGELOG.md

@@ -6,7 +6,82 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased]
 
+### Changed
+- **The Gaussian likelihood behind AICc is evaluated in log-RSS space, and an
+  exact-zero RSS is now an explicit abstention (issue #135).** This is a
+  METHODOLOGY change with user-visible consequences: it can reorder AICc,
+  change which of M0..M3b is selected, and therefore change the reported decay
+  rate for a given run. Three parts:
+  - `gaussian_log_likelihood` evaluates the profile likelihood at the Gaussian
+    MLE directly from `log(RSS)` via the new
+    `numerics.norms.scaled_log_sum_squares`, instead of forming RSS in ordinary
+    units. Residual series whose true sum of squares lies outside the float64
+    range — a rescaled curve is enough — previously produced `0` or `inf` RSS
+    and an unusable likelihood. The value is now scale-covariant, and
+    likelihood *differences* between models are invariant under a common
+    rescaling of the residuals, which is the quantity AICc actually consumes.
+  - An exact-zero RSS has no finite interior MLE for the positive scale
+    parameter, so the profile likelihood is `NaN` rather than an invented
+    absolute epsilon variance. `fit_gls_ar1` reports this as
+    `likelihood_degenerate=True` with `success=False`, `aicc()` then scores the
+    model `inf`, and `parametric_bootstrap` refuses: a perfect fit yields no
+    interval rather than a zero-width one. Passing an explicit positive `sigma`
+    still gives a finite likelihood for zero residuals, which is the
+    well-posed case.
+  - Runs whose reported model or rate changes are those where the previous
+    absolute-RSS path had over- or underflowed, or where the winner was chosen
+    against a floored likelihood. Re-run any archived analysis whose selected
+    model matters; the run manifest does not record per-model likelihoods, so
+    the change is not visible in `input_hash`.
+
 ### Fixed
+- **The reason a confidence interval was withheld did not reach the report (PR
+  #147, round-2 external review).** When the residual MLE scale is not
+  representable as float64 the fit stays a valid AICc candidate and only its
+  interval is withheld — but `_fit_with_model` copied `likelihood_degenerate`
+  and dropped `scale_unavailable`, so the persisted report showed only
+  `bca_ci_beta = (nan, nan)`. That is the value ANY bootstrap or jackknife
+  failure produces, and the distinguishing information lived in a
+  `RuntimeWarning` an artefact does not keep. `FitResult` carries
+  `scale_unavailable` now and `_fit_with_model` copies it; `io.export`
+  serialises `FitResult` field-wise, so the reason travels into dumped
+  reports.
+- **The one-half factor was applied after the overflow guard (PR #147, round-2
+  external review).** For the explicit-`sigma` path the standardised RSS never
+  appears on its own: it enters the log-likelihood as `-0.5 * RSS`, so the
+  representable range reaches `2 * float64.max`. The guard compared against
+  `log(float64.max)` and returned `-inf` for the octave above it — measured:
+  `sigma = 1` with a single residual near `1.4e154` has RSS ~`1.96e308` and a
+  finite log-likelihood of ~`-9.8e307`, and the model was dropped from
+  selection on that arithmetic boundary rather than on the data. The bound now
+  carries `+ log(2)` and the exponential is taken after subtracting `log(2)`
+  in the octave that needs it; below that octave the arithmetic is unchanged,
+  so no previously computable likelihood moves by even one ulp.
+- **A fit was withheld entirely because a number its likelihood never uses
+  could not be materialised (PR #147, round-1 review).** When the whitened
+  residuals have a finite, non-zero RMS below the smallest positive float64,
+  `log_rss` and `log_sigma` are both finite and the profile likelihood is
+  computable — but `exp(log_sigma)` underflows to `0.0`, and `fit_gls_ar1`
+  then returned `success=False`, `log_likelihood=NaN` and
+  `likelihood_degenerate=True`. `_fit_with_model` scored that model `inf` and
+  dropped it from selection, which reintroduces exactly the ABSOLUTE SCALE
+  BOUNDARY into model selection that the change above removed: the same curve
+  in different rate units either is or is not a candidate. Measured on one
+  minimum-subnormal residual (`5e-324`) among 64 otherwise-zero points:
+  `log_rss = -1488.88`, `log_sigma = -746.52`, profile log-likelihood
+  `+47686.44`, `exp(log_sigma) = 0.0`.
+
+  The fit now stays selectable with its finite log-space likelihood, and only
+  the scale-dependent evidence is withheld. `GLSFitOutput` gains
+  `scale_unavailable` (additive, default `False`); `sigma` is `NaN` there,
+  deliberately not `0.0`, because `_ar1_resample` consumes it as the
+  innovation standard deviation and `0.0` would generate identical replicates
+  — a zero-width confidence interval, which is the failure mode of an
+  uncertainty pipeline, not a conservative one. `parametric_bootstrap` refuses
+  on the new flag rather than on `success`, so "no interval" and "no estimate"
+  stay distinguishable, and `compute_relaxation_layer` reports the CI as `NaN`
+  through its existing handler. Three mutations, one per guard line, are
+  proven to discriminate.
 - **A degenerate spectrum let one eigenvector certify two modes (PR #121,
   round-26 external review, P1).** `certified_nonzero_modes` picked the
   eigenpair for an in-band candidate with `argmin(|ref - lambda|)`, a lookup by

CITATION.cff

@@ -84,6 +84,23 @@ abstract: >-
 # likewise scale-relative rather than absolute (issue #109), which closes a
 # fail-open path that admitted non-GKSL generators at small ||H||; this changes
 # which inputs are ACCEPTED, not the numerics of accepted ones.
+# Also pending for the next cut: the Gaussian likelihood behind AICc is
+# evaluated in log-RSS space rather than by forming the residual sum of squares
+# in ordinary units, and an exact-zero RSS is an explicit abstention rather than
+# a floored variance (issue #135). This CHANGES NUMERICAL RESULTS and is a
+# METHODOLOGY correction, not an addition: AICc ordering, the selected model
+# M0..M3b and therefore the reported decay rate can differ for runs whose
+# residual scale previously over- or underflowed, and a perfect fit now yields
+# no confidence interval instead of a zero-width one. The likelihood value is
+# scale-covariant and likelihood DIFFERENCES between models are invariant under
+# a common rescaling of the residuals, which is the quantity AICc consumes. It
+# must NOT be described as making model selection scale-free in general: only
+# the likelihood path was corrected, and the residual-scale parameter itself
+# still has a float64 representation limit -- a fit whose MLE scale underflows
+# keeps its log-space likelihood and stays selectable, but its bootstrap
+# interval is withheld (round-1 review of PR #147). Archived analyses whose
+# selected model matters should be re-run; per-model likelihoods are not run-
+# manifest fields, so this change is not visible in input_hash.
 keywords:
   - open quantum systems
   - Lindblad

src/liouscope/_types.py

@@ -187,6 +187,16 @@ class FitResult:
     n_eff: float
     residual_ar1_rho: float
     success: bool
+    likelihood_degenerate: bool = False
+    #: True when the fit is a valid model-selection candidate -- its log-space
+    #: likelihood and AICc are finite -- but the positive MLE residual scale is
+    #: not representable as float64, so the parametric bootstrap and therefore
+    #: the confidence interval were withheld (issue #135, PR #147 round-2
+    #: review). Without this field the persisted report showed only
+    #: ``bca_ci_beta = (nan, nan)``, which is the same value ANY bootstrap or
+    #: jackknife failure produces: once the warning stream is gone, the reason
+    #: the interval is missing was unrecoverable from the artefact.
+    scale_unavailable: bool = False
 
 
 @dataclass(frozen=True, slots=True, kw_only=True)

src/liouscope/diagnostics/relaxation.py

@@ -215,6 +215,16 @@ def _fit_with_model(
         n_eff=n_eff,
         residual_ar1_rho=fit.rho_ar1,
         success=fit.success,
+        likelihood_degenerate=fit.likelihood_degenerate,
+        # ROUND-2 REVIEW (PR #147). This conversion carried
+        # ``likelihood_degenerate`` and dropped ``scale_unavailable``, so the
+        # newly supported underflow case -- ``success=True`` with a withheld
+        # ``sigma`` -- reached the report indistinguishable from an ordinary
+        # fit. Its bootstrap then fails and ``compute_relaxation_layer``
+        # records ``bca_ci_beta = (nan, nan)``, exactly what a jackknife
+        # failure or a non-converged resample produces. The reason lived only
+        # in a RuntimeWarning, which a persisted artefact does not keep.
+        scale_unavailable=fit.scale_unavailable,
     )
     return fit_result, fit.params
 

src/liouscope/fitting/aicc.py

@@ -11,10 +11,13 @@
 
 from __future__ import annotations
 
+import math
 from collections.abc import Mapping
 
 import numpy as np
 
+from ..numerics.norms import scaled_log_sum_squares
+
 
 def aicc(log_likelihood: float, k: int, n_eff: float) -> float:
     """Return ``AICc`` for given log-likelihood, parameter count, and N_eff."""
@@ -29,18 +32,109 @@ def gaussian_log_likelihood(
     *,
     sigma: float | None = None,
 ) -> float:
-    """Gaussian log-likelihood for ``y - y_hat``.
+    """Gaussian log-likelihood for ``y - y_hat`` without absolute RSS floors.
 
-    Uses MLE sigma if ``sigma`` is omitted.
+    When ``sigma`` is omitted, evaluate the profile likelihood at the Gaussian
+    MLE ``sigma_hat**2 = RSS / n`` directly in log-RSS space. Exact zero RSS has
+    no finite interior MLE for the positive scale parameter and therefore
+    returns NaN (model-selection likelihood unavailable) instead of inventing
+    an absolute epsilon variance. With an explicitly supplied finite positive
+    ``sigma``, zero residuals remain a valid finite likelihood.
     """
     residuals = np.asarray(residuals, dtype=float)
     n = residuals.size
-    rss = float(np.dot(residuals, residuals))
+    if n == 0:
+        return float("nan")
+
+    log_rss = scaled_log_sum_squares(residuals)
+    if math.isnan(log_rss):
+        return float("nan")
+
+    log_2pi = math.log(2.0 * math.pi)
     if sigma is None:
-        sigma_sq = max(rss / n, 1.0e-30)
+        if log_rss == float("-inf"):
+            return float("nan")
+        if log_rss == float("inf"):
+            return float("-inf")
+        return float(-0.5 * n * (log_2pi + 1.0 + log_rss - math.log(n)))
+
+    sigma = float(sigma)
+    if not math.isfinite(sigma) or sigma <= 0.0:
+        return float("nan")
+    log_sigma = math.log(sigma)
+    if log_rss == float("-inf"):
+        half_standardised_rss = 0.0
+    elif log_rss == float("inf"):
+        return float("-inf")
     else:
-        sigma_sq = sigma * sigma
-    return float(-0.5 * n * (np.log(2.0 * np.pi * sigma_sq) + rss / (n * sigma_sq)))
+        log_standardised_rss = log_rss - 2.0 * log_sigma
+        # ROUND-2 REVIEW (PR #147). The standardised RSS never appears on its
+        # own: it enters the result as ``-0.5 * RSS``. Its representable range
+        # therefore reaches ``2 * float64.max``, and comparing against
+        # ``log(float64.max)`` refused an entire octave of perfectly finite
+        # likelihoods. Measured: ``sigma = 1`` with a single residual near
+        # 1.4e154 has RSS ~1.96e308 and a log-likelihood of ~-9.8e307, and
+        # this guard returned ``-inf`` -- which drops the model from AICc
+        # selection on an arithmetic accident rather than on the data, the
+        # same absolute-scale boundary issue #135 exists to remove.
+        #
+        # The one-half factor is now applied BEFORE the decision: the bound
+        # carries ``+ log(2)`` and, in the octave that only the halved value
+        # can represent, the exponential is taken after subtracting
+        # ``log(2)``. Below that octave the arithmetic is left exactly as it
+        # was (``0.5 * exp(...)``), so no previously computed likelihood
+        # changes by even one ulp -- widening a range must not perturb the
+        # values already inside it.
+        # ROUND-3 REVIEW (PR #147). The round-2 bound was still a decision
+        # made in LOG SPACE, and at this magnitude log space cannot make it.
+        # Measured with ``sigma = 1`` and a single residual
+        # ``sqrt(float_max) * sqrt(2)``: ``log_rss`` comes out at
+        # 710.475860073944 while ``log_max + log(2)`` rounds to
+        # 710.4758600739439 -- ONE ulp lower -- so the guard returned
+        # ``-inf`` although the mathematically relevant half-RSS is
+        # 1.7976931348623155e+308, a perfectly finite 0.9999999999999999 of
+        # ``float_max``. Nor can the halved value be recovered by subtracting
+        # ``log(2)`` and exponentiating: one ulp of a logarithm near 710 is a
+        # factor of ~1e-16 in the value, which straddles the overflow edge
+        # exactly where this test is made, and ``math.exp`` raises
+        # ``OverflowError`` on the same input (measured).
+        #
+        # So the halved quantity is MATERIALISED from the residuals instead
+        # of inferred from their logarithm. Halving each term BEFORE
+        # accumulating keeps every partial sum inside float64 whenever the
+        # true half-RSS is inside it, and lets the float64 addition itself
+        # decide the boundary -- there is no rounded bound left to be one ulp
+        # wrong about. An overflow here is now a measurement rather than a
+        # prediction: ``inf`` means the half-RSS genuinely does not fit
+        # (verified for 4x, 1e3x float_max, and for 64 residuals whose halved
+        # sum overflows), and only then is the fit dropped.
+        #
+        # The lower octave is untouched, so no previously computed likelihood
+        # changes by even one ulp -- the same promise round 2 made.
+        log_max = math.log(np.finfo(float).max)
+        if log_standardised_rss > log_max:
+            scaled = residuals / sigma
+            with np.errstate(over="ignore", invalid="ignore"):
+                half_standardised_rss = float(np.sum((0.5 * scaled) * scaled))
+            # Tested for OVERFLOW specifically, not for non-finiteness.
+            # The mutation run reported a blanket isfinite() check BLIND, and
+            # measuring why exposed a latent state collapse: for an infinite
+            # half-RSS the check is provably equivalent (the subtraction
+            # below already yields -inf), and its only remaining effect would
+            # be to turn a NaN into -inf -- converting "model-selection
+            # likelihood unavailable", which this function returns NaN for by
+            # contract, into "definitively dropped". Those are different
+            # states. NaN is unreachable here (a NaN residual is caught by
+            # the isnan(log_rss) gate above, and every term 0.5*s*s is
+            # non-negative), so this is guarding an invariant rather than a
+            # path -- which is the reason to state it precisely.
+            if half_standardised_rss == float("inf"):
+                return float("-inf")
+        else:
+            half_standardised_rss = 0.5 * math.exp(log_standardised_rss)
+    return float(
+        -0.5 * n * log_2pi - n * log_sigma - half_standardised_rss
+    )
 
 
 def choose_model(aiccs: Mapping[str, float]) -> str:

src/liouscope/fitting/bootstrap.py

@@ -52,6 +52,21 @@ def parametric_bootstrap(
     if rng is None:
         rng = np.random.default_rng(0)
     base = fit_gls_ar1(model, t, y, p0, bounds=bounds)
+    if base.scale_unavailable:
+        # Round-1 review (PR #147). The base fit is a legitimate AICc
+        # candidate -- its log-space likelihood is finite -- but its innovation
+        # scale is not representable, so there is nothing to draw replicates
+        # from. ``sigma`` is NaN there by construction, and ``rng.normal(0, nan)``
+        # returns NaN rather than raising, so every replicate would be fitted to
+        # NaN data and the resulting interval would be an artefact. Refusing
+        # here routes the case into the caller's existing handler, which reports
+        # the CI as NaN -- "fit uncertainty UNKNOWN", which is the true state.
+        raise RuntimeError(
+            "parametric_bootstrap: the base fit has no representable residual "
+            "scale (issue #135), so replicates cannot be simulated; the fit "
+            "itself remains a valid model-selection candidate, only its "
+            "interval is unavailable"
+        )
     if not base.success:
         # Round-17 review (PR #121). Every replicate is simulated AROUND
         # ``theta_hat``; if the base fit ended on the model's magnitude
@@ -65,6 +80,8 @@ def parametric_bootstrap(
             + (f" (saturated: {', '.join(base.saturated)})" if base.saturated else "")
             + (" (the curve carries no resolvable variation, issue #123)"
                if base.degenerate else "")
+            + (" (the residual likelihood scale is degenerate, issue #135)"
+               if base.likelihood_degenerate else "")
             + "; a bootstrap around a non-estimate has no meaning"
         )
     theta_hat = base.params