diff --git a/.github/workflows/pypi.yml b/.github/workflows/pypi.yml index 330d4ba..32033c3 100644 --- a/.github/workflows/pypi.yml +++ b/.github/workflows/pypi.yml @@ -71,6 +71,7 @@ jobs: url: https://pypi.org/p/liouscope permissions: id-token: write + attestations: write contents: read steps: - uses: actions/checkout@11bd71901bbe5b1630ceea73d27597364c9af683 # v4.2.2 @@ -103,6 +104,18 @@ jobs: test "$SOURCE_VERSION" = "$EXPECTED_VERSION" test "$HEAD_SHA" = "$TAG_SHA" test "$EVENT_SHA" = "$TAG_SHA" + - name: Generate build provenance attestation + uses: actions/attest-build-provenance@0f67c3f4856b2e3261c31976d6725780e5e4c373 # v4.1.1 + with: + subject-path: "dist/*" + - name: Verify build provenance attestation + env: + GH_TOKEN: ${{ github.token }} + run: | + set -euo pipefail + for artifact in dist/*; do + gh attestation verify "$artifact" --repo "${{ github.repository }}" + done - name: Publish via Trusted Publishing uses: pypa/gh-action-pypi-publish@cef221092ed1bacb1cc03d23a2d87d1d172e277b # v1.14.0 with: diff --git a/CHANGELOG.md b/CHANGELOG.md index 9968fbc..9804c9c 100644 --- a/CHANGELOG.md +++ b/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. - **D16 was published from the spectrum for which D1/D3/D4 had just been withheld (PR #121, round-23 review, finding 13).** Where the zero-mode certificate is applicable but unresolved, the spectral layer reports D1, D3, diff --git a/CITATION.cff b/CITATION.cff index 18b5fac..202b5eb 100644 --- a/CITATION.cff +++ b/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 diff --git a/docs/GOVERNANCE_AUDIT_2026-09-01.md b/docs/GOVERNANCE_AUDIT_2026-09-01.md new file mode 100644 index 0000000..d107ea6 --- /dev/null +++ b/docs/GOVERNANCE_AUDIT_2026-09-01.md @@ -0,0 +1,48 @@ +# LiouScope governance audit — 2026-09-01 + +Repository: `marcohost33-maker/Liouscope` +Observed `main`: `386041b4072a776a985beb34af5815fb344c74f9` + +## Scope + +This is a point-in-time repository-governance evidence lock. It does not amend the historical scientific or packaging claims of `RELEASE_AUDIT_v0.5.0.md` and does not certify settings that the available integration cannot read. + +## Required checks observed on `main` + +GitHub branch metadata reports `main` as protected, enforcement level `everyone`, with these required status checks: + +- `test (ubuntu-latest, 3.10)` +- `test (ubuntu-latest, 3.11)` +- `test (ubuntu-latest, 3.12)` +- `test (ubuntu-latest, 3.13)` +- `test (ubuntu-latest, 3.14)` +- `qutip-cross-check (3.11)` +- `qutip-cross-check (3.12)` +- `quality contract` + +## Negative-control proof + +Draft PR #141 intentionally adds one inert `workflow_dispatch`-only workflow containing a mutable `actions/checkout@main` reference. The `Quality Contract` workflow completed with `failure`; its `quality contract` job failed specifically at `Check workflow hardening`. + +This demonstrates that the aggregate required check receives and propagates a repository-policy failure. The control workflow is deliberately non-mergeable test material and must be removed after the evidence is captured. + +## Positive / non-workflow proof + +This audit and the accompanying `QUALITY_WORKFLOW_OS.md` update are documentation-only. Their PR is used to verify that the required `quality contract` is emitted for a non-workflow change and reaches a terminal result rather than remaining indefinitely `Expected`. + +## Explicitly unverified protection settings + +The integration can read the branch summary but receives HTTP 403 for the detailed branch-protection endpoint. Therefore this audit does not claim independently verified values for: + +- required approving-review count; +- force-push allowance; +- branch deletion allowance; +- administrator bypass; +- linear-history or signed-commit requirements; +- conversation-resolution requirements. + +These settings remain a repository-settings verification item if a release or governance claim depends on them. + +## Verdict + +`quality contract` is observed as a required status check and its fail-closed path is load-bearing. Issue #133 can be closed once the documentation-only positive-control PR is green and the negative-control PR has been closed/reset without merging its deliberate violation. diff --git a/docs/QUALITY_WORKFLOW_OS.md b/docs/QUALITY_WORKFLOW_OS.md index 94b0e3b..5d3025b 100644 --- a/docs/QUALITY_WORKFLOW_OS.md +++ b/docs/QUALITY_WORKFLOW_OS.md @@ -103,6 +103,28 @@ Quality claims must include at least one success metric and one counter-metric: | release frequency | failed release recovery time | | coverage | escaped defect rate | +### Current enforced `main` contract — observed 2026-09-01 + +This is an **observed governance snapshot**, not a timeless claim. At `main@386041b4072a776a985beb34af5815fb344c74f9`, GitHub's branch metadata reports `main` as protected with required-status-check enforcement level `everyone` and the following required checks: + +- `test (ubuntu-latest, 3.10)` +- `test (ubuntu-latest, 3.11)` +- `test (ubuntu-latest, 3.12)` +- `test (ubuntu-latest, 3.13)` +- `test (ubuntu-latest, 3.14)` +- `qutip-cross-check (3.11)` +- `qutip-cross-check (3.12)` +- `quality contract` + +The `quality contract` is therefore merge-boundary evidence, not merely advisory workflow output, for this observed configuration. + +The enforcement path is tested in both directions: + +- ordinary PR heads have produced `quality contract = success`; +- controlled negative PR #141 added an inert, `workflow_dispatch`-only workflow with an intentionally mutable `actions/checkout@main` reference; its `quality contract` failed specifically at `Check workflow hardening`. The violating workflow is test evidence only and must never be merged. + +A connector-access limitation remains explicit: the repository branch endpoint exposes the required-check set above, but the integration receives HTTP 403 for the full branch-protection detail endpoint. This audit therefore does **not** infer or certify review-count, force-push, deletion, administrator-bypass, or other protection settings that were not independently observable. Those settings must be verified in GitHub repository settings before any release process depends on them. + ## 3. Repository Quality Delta Score (RQDS v0.2) `RQDS = 0.20*CI_Reliability + 0.20*Security_Posture + 0.15*Evidence_Coverage + 0.15*Maintainability + 0.15*Delivery_Stability + 0.10*Observability + 0.05*Cost_Discipline - Penalty` diff --git a/src/liouscope/_types.py b/src/liouscope/_types.py index f689393..67ddbf2 100644 --- a/src/liouscope/_types.py +++ b/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) diff --git a/src/liouscope/diagnostics/relaxation.py b/src/liouscope/diagnostics/relaxation.py index d8433bc..666dddd 100644 --- a/src/liouscope/diagnostics/relaxation.py +++ b/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 diff --git a/src/liouscope/fitting/aicc.py b/src/liouscope/fitting/aicc.py index 2130a6d..f078164 100644 --- a/src/liouscope/fitting/aicc.py +++ b/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: diff --git a/src/liouscope/fitting/bootstrap.py b/src/liouscope/fitting/bootstrap.py index b936ee5..1ec84e1 100644 --- a/src/liouscope/fitting/bootstrap.py +++ b/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 diff --git a/src/liouscope/fitting/gls.py b/src/liouscope/fitting/gls.py index 73b54f4..7440643 100644 --- a/src/liouscope/fitting/gls.py +++ b/src/liouscope/fitting/gls.py @@ -12,6 +12,7 @@ from __future__ import annotations +import math import warnings from collections.abc import Callable from dataclasses import dataclass @@ -19,6 +20,7 @@ import numpy as np from scipy.optimize import least_squares +from ..numerics.norms import scaled_log_sum_squares from .aicc import gaussian_log_likelihood from .models import saturation_watch from .neff import _AR1_SMALL_N, ar1_correlation_corrected @@ -41,6 +43,19 @@ class GLSFitOutput: #: reports a fit that ran and ended on a magnitude plateau; here there was #: nothing to fit. Implies ``success`` is False and ``params`` is NaN. degenerate: bool = False + #: True when the residual Gaussian scale has no finite usable MLE for model + #: selection (issue #135). Distinct from ``degenerate`` above: the curve + #: may carry variation and the optimiser may have run, but exact-zero RSS + #: (or an unrepresentable positive MLE scale) cannot support AICc/CI claims. + likelihood_degenerate: bool = False + #: True when the profile likelihood IS computable in log space but the + #: positive MLE scale ``exp(log_sigma)`` is not representable as a float64 + #: (round-1 review of PR #147). ``sigma`` is then NaN while + #: ``log_likelihood`` stays finite: the fit remains a valid AICc candidate, + #: and only the scale-dependent evidence -- the parametric bootstrap, hence + #: the CI -- is withheld. Distinct from ``likelihood_degenerate``, where the + #: likelihood itself has no finite value. + scale_unavailable: bool = False def _whiten(y: np.ndarray, rho: float) -> np.ndarray: @@ -192,20 +207,83 @@ def residual(params: np.ndarray, rho_local: float = rho) -> np.ndarray: ) whitened = _whiten(residuals_raw, rho) n = whitened.size - sigma = float(np.sqrt(max(np.dot(whitened, whitened) / max(n, 1), 1.0e-30))) + log_rss = scaled_log_sum_squares(whitened) + if log_rss == float("-inf") or not math.isfinite(log_rss): + warnings.warn( + "fit_gls_ar1: residual Gaussian scale has no finite positive MLE " + "for model selection; likelihood/AICc/CI evidence is unavailable " + "(issue #135).", + RuntimeWarning, + stacklevel=2, + ) + return GLSFitOutput( + params=p, + residuals=residuals_raw, + rho_ar1=rho, + sigma=float("nan"), + log_likelihood=float("nan"), + success=False, + saturated=tuple(sorted(fired)), + likelihood_degenerate=True, + ) + + log_sigma = 0.5 * (log_rss - math.log(n)) + try: + sigma = float(math.exp(log_sigma)) + except OverflowError: + sigma = float("inf") + scale_unavailable = not math.isfinite(sigma) or sigma <= 0.0 + if scale_unavailable: + # ROUND-1 REVIEW (PR #147). This branch used to return + # ``success=False, log_likelihood=nan, likelihood_degenerate=True``, + # which withheld the whole FIT because a number the likelihood never + # needs could not be materialised. ``log_rss`` is finite here -- the + # gate above already refused the case where it is not -- so the profile + # likelihood is computable directly in log space by + # ``gaussian_log_likelihood``, which never forms ``sigma``. Only + # ``exp(log_sigma)`` left the float64 range. + # + # Measured on the reviewer's construction: one minimum-subnormal + # residual (5e-324) in an otherwise-zero series of 64 points gives + # ``log_rss = -1488.88``, ``log_sigma = -746.52`` and a perfectly + # finite profile log-likelihood of ``+47686.44`` -- while + # ``exp(-746.52)`` underflows to ``0.0`` because the smallest + # subnormal is ``exp(-744.44)``. Marking the fit unsuccessful there + # makes ``aicc()`` return ``inf`` in ``_fit_one`` and the model drops + # out of selection, which reintroduces exactly the ABSOLUTE SCALE + # BOUNDARY into model selection that issue #135 removed: the same + # curve in different rate units is or is not a candidate. + # + # So the fit stays selectable and only the scale-dependent evidence is + # withheld. ``sigma`` is NaN, never 0.0 or inf: it is consumed by + # ``_ar1_resample`` as the standard deviation of the innovation, where + # 0.0 would silently generate a bootstrap of IDENTICAL replicates -- + # a zero-width CI, which is the failure mode of an uncertainty + # pipeline, not a wide one. ``parametric_bootstrap`` refuses on the + # dedicated flag rather than on ``success``, so "no interval" and "no + # estimate" stay distinguishable. + warnings.warn( + "fit_gls_ar1: the positive residual MLE scale is not representable " + f"as float64 (log sigma = {log_sigma:.6g}); the log-space " + "likelihood and AICc remain available, but bootstrap/CI evidence " + "is withheld (issue #135, PR #147 review).", + RuntimeWarning, + stacklevel=2, + ) + # Prais-Winsten exact AR(1) likelihood: _whiten keeps observation 0 # (scaled by sqrt(1-rho^2)), so the transform has log-Jacobian - # 0.5*log(1-rho^2). Omitting it makes the reported value a hybrid of the - # exact and conditional likelihoods and biases cross-model AICc (each - # model fits its own rho) toward under-fitting high-rho models. + # 0.5*log(1-rho^2). The profile likelihood is evaluated directly from + # log(RSS), not by squaring ``sigma`` or materialising RSS. jac = 0.5 * float(np.log(max(1.0 - rho * rho, 1.0e-12))) - log_lik = gaussian_log_likelihood(whitened, sigma=sigma) + jac + log_lik = gaussian_log_likelihood(whitened) + jac return GLSFitOutput( params=p, residuals=residuals_raw, rho_ar1=rho, - sigma=sigma, + sigma=float("nan") if scale_unavailable else sigma, log_likelihood=log_lik, success=success, saturated=tuple(sorted(fired)), + scale_unavailable=scale_unavailable, ) diff --git a/src/liouscope/numerics/norms.py b/src/liouscope/numerics/norms.py index 067d801..712973f 100644 --- a/src/liouscope/numerics/norms.py +++ b/src/liouscope/numerics/norms.py @@ -73,6 +73,40 @@ def scaled_euclidean_norm(values: np.ndarray) -> float: return float(np.ldexp(scaled_norm, exponent)) +def scaled_log_sum_squares(values: np.ndarray) -> float: + """Return ``log(sum(abs(values)**2))`` without spurious under/overflow. + + The return value is ``-inf`` for an exact all-zero input, ``nan`` when any + component is NaN, and ``inf`` when any component is infinite. For finite, + non-zero float64 input the logarithm remains finite even when the true sum + of squares (or its square root) is outside the representable float64 range. + + This is the likelihood-facing companion to :func:`scaled_euclidean_norm`: + both use the same exact power-of-two scaling contract, but this function + never reconstructs RSS in ordinary floating-point units. + """ + arr = np.asarray(values) + if arr.size == 0: + return float("-inf") + + real = np.asarray(np.real(arr), dtype=float) + imag = np.asarray(np.imag(arr), dtype=float) + if np.any(np.isnan(real)) or np.any(np.isnan(imag)): + return float("nan") + if np.any(np.isinf(real)) or np.any(np.isinf(imag)): + return float("inf") + + scaled = _finite_component_scale(arr) + if scaled is None: + return float("-inf") + scaled_real, scaled_imag, exponent = scaled + sumsq = float( + np.sum(scaled_real * scaled_real, dtype=float) + + np.sum(scaled_imag * scaled_imag, dtype=float) + ) + return float(math.log(sumsq) + 2.0 * exponent * math.log(2.0)) + + def scaled_cancellation_ratio(values: np.ndarray) -> float: """Return ``abs(sum(values)) / sum(abs(values))`` scale-safely. diff --git a/tests/test_issue135_likelihood_scale.py b/tests/test_issue135_likelihood_scale.py new file mode 100644 index 0000000..8231541 --- /dev/null +++ b/tests/test_issue135_likelihood_scale.py @@ -0,0 +1,417 @@ +"""Issue #135: scale-safe Gaussian likelihood and degenerate RSS semantics.""" + +from __future__ import annotations + +import math +import warnings + +import numpy as np +import pytest + +from liouscope.diagnostics import relaxation as relaxation_mod +from liouscope.fitting.aicc import aicc, choose_model, gaussian_log_likelihood +from liouscope.fitting.bootstrap import parametric_bootstrap +from liouscope.fitting.gls import GLSFitOutput, fit_gls_ar1 +from liouscope.fitting.models import M0 +from liouscope.io.export import _to_jsonable +from liouscope.numerics.norms import scaled_log_sum_squares + + +def test_scaled_log_rss_spans_underflow_and_overflow_regimes(): + tiny = np.array([1.0e-320, -2.0e-320]) + huge = np.array([1.0e308, -1.0e308]) + log_tiny = scaled_log_sum_squares(tiny) + log_huge = scaled_log_sum_squares(huge) + assert np.isfinite(log_tiny) + assert np.isfinite(log_huge) + assert log_tiny == pytest.approx(math.log(5.0) + 2.0 * math.log(1.0e-320), rel=2e-5) + assert log_huge == pytest.approx(math.log(2.0) + 2.0 * math.log(1.0e308), rel=1e-14) + assert scaled_log_sum_squares(np.zeros(4)) == float("-inf") + + +def test_profile_loglikelihood_is_scale_covariant_and_delta_invariant(): + x = np.linspace(0.2, 1.2, 64) + r0 = 0.7 + x + 0.03 * np.sin(5.0 * x) + r1 = 0.9 + 1.1 * x - 0.02 * np.cos(3.0 * x) + reference_delta = None + reference_winner = None + for scale in (1.0e-150, 1.0e-40, 1.0, 1.0e40, 1.0e150): + ll0 = gaussian_log_likelihood(scale * r0) + ll1 = gaussian_log_likelihood(scale * r1) + assert np.isfinite(ll0) and np.isfinite(ll1) + delta = ll1 - ll0 + scores = { + "M0": aicc(ll0, k=2, n_eff=64.0), + "M1": aicc(ll1, k=3, n_eff=64.0), + } + winner = choose_model(scores) + if reference_delta is None: + reference_delta = delta + reference_winner = winner + else: + assert delta == pytest.approx(reference_delta, rel=1e-11, abs=1e-10) + assert winner == reference_winner + + +def test_profile_loglikelihood_handles_true_rss_above_float_range(): + ll = gaussian_log_likelihood(np.array([1.0e308, -1.0e308])) + assert np.isfinite(ll) + + +def test_zero_rss_unknown_sigma_is_unavailable_but_known_sigma_is_valid(): + residuals = np.zeros(8) + assert np.isnan(gaussian_log_likelihood(residuals)) + ll = gaussian_log_likelihood(residuals, sigma=2.0) + expected = -0.5 * residuals.size * math.log(2.0 * math.pi * 4.0) + assert ll == pytest.approx(expected, rel=1e-14, abs=1e-14) + + +def test_gls_exact_fit_marks_likelihood_degenerate_and_bootstrap_refuses(): + t = np.linspace(0.0, 4.0, 48) + p0 = np.array([1.25, 0.6]) + y = M0(t, p0) + with pytest.warns(RuntimeWarning, match="likelihood/AICc/CI evidence is unavailable"): + fit = fit_gls_ar1(M0, t, y, p0, n_iters=1) + assert not fit.success + assert fit.likelihood_degenerate + assert not fit.degenerate + assert np.isnan(fit.sigma) + assert np.isnan(fit.log_likelihood) + # The repository treats unexpected warnings as errors. Here the warning is + # part of the intended public contract: the base fit first reports why its + # likelihood evidence is unusable, then bootstrap refuses the non-estimate. + with pytest.warns(RuntimeWarning, match="likelihood/AICc/CI evidence is unavailable"): + with pytest.raises(RuntimeError, match="likelihood scale is degenerate"): + parametric_bootstrap(M0, t, y, p0, B=4) + + +def test_ordinary_noisy_gls_positive_control_is_not_likelihood_degenerate(rng): + t = np.linspace(0.0, 4.0, 80) + p = np.array([1.25, 0.6]) + y = M0(t, p) + 1.0e-3 * rng.standard_normal(t.size) + fit = fit_gls_ar1(M0, t, y, p, n_iters=1) + assert fit.success + assert not fit.likelihood_degenerate + assert np.isfinite(fit.sigma) and fit.sigma > 0.0 + assert np.isfinite(fit.log_likelihood) + + +def test_likelihood_degenerate_state_reaches_fitresult_and_is_nonselectable(monkeypatch): + t = np.linspace(0.0, 1.0, 16) + y = np.exp(-t) + fake = GLSFitOutput( + params=np.array([1.0, 1.0]), + residuals=np.zeros_like(y), + rho_ar1=0.0, + sigma=float("nan"), + log_likelihood=float("nan"), + success=False, + likelihood_degenerate=True, + ) + monkeypatch.setattr(relaxation_mod, "fit_gls_ar1", lambda *args, **kwargs: fake) + fit_result, _ = relaxation_mod._fit_with_model("M0", t, y) + assert not fit_result.success + assert fit_result.likelihood_degenerate + assert np.isinf(fit_result.aicc) + + +# -------------------------------------------------------------------------- +# PR #147 round-1 review: an unrepresentable MLE scale withheld the whole FIT, +# not just the scale-dependent evidence. +# -------------------------------------------------------------------------- + +#: Smallest positive float64. One such residual among an otherwise-zero series +#: keeps ``log_rss`` finite while ``exp(0.5 * (log_rss - log n))`` underflows: +#: the profile likelihood is computable, its scale is not. +_MIN_SUBNORMAL = 5.0e-324 + + +def _constant_model(t: np.ndarray, p: np.ndarray) -> np.ndarray: + return np.full(t.shape, float(p[0])) + + +def _underflowing_scale_case() -> tuple[np.ndarray, np.ndarray]: + n = 64 + t = np.linspace(0.0, 5.0, n) + y = np.zeros(n) + y[7] = _MIN_SUBNORMAL + return t, y + + +def test_the_underflowing_case_really_has_a_finite_log_likelihood() -> None: + """Positive control: the fixture must exercise the finding, not a NaN RSS. + + Without this the two assertions below could pass on a case where nothing + was computable in the first place. + """ + _, y = _underflowing_scale_case() + log_rss = scaled_log_sum_squares(y) + assert math.isfinite(log_rss), log_rss + log_sigma = 0.5 * (log_rss - math.log(y.size)) + assert math.isfinite(log_sigma), log_sigma + assert math.exp(log_sigma) == 0.0, ( + "the fixture no longer underflows; the finding cannot be reached" + ) + assert math.isfinite(gaussian_log_likelihood(y)) + + +def test_unrepresentable_scale_keeps_the_fit_selectable() -> None: + """The fit stays an AICc candidate; only ``sigma`` is withheld. + + Measured before the repair: ``success=False``, ``log_likelihood=nan`` and + ``likelihood_degenerate=True`` -- so ``_fit_with_model`` scored the model + ``inf`` and dropped it from selection. That reintroduces an ABSOLUTE scale + boundary into model selection, which is the boundary issue #135 removed: + the same curve in different rate units is or is not a candidate. + """ + t, y = _underflowing_scale_case() + # Recorded rather than ``pytest.warns``: a missing warning must fail this + # test at an ASSERTION, so that removing the guard is attributable. With + # ``pytest.warns`` the test dies inside pytest's own context manager, which + # a mutation run cannot tell apart from an incidental crash. + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + fit = fit_gls_ar1(_constant_model, t, y, np.array([0.0])) + assert any( + issubclass(w.category, RuntimeWarning) + and "scale is not representable" in str(w.message) + for w in caught + ), [str(w.message) for w in caught] + assert fit.success, "the fit was withheld because its SCALE could not be built" + assert math.isfinite(fit.log_likelihood), fit.log_likelihood + assert not fit.likelihood_degenerate + assert fit.scale_unavailable + assert math.isnan(fit.sigma), ( + f"sigma is {fit.sigma!r}; 0.0 would make _ar1_resample draw identical " + "replicates and produce a zero-width interval" + ) + + +def test_unrepresentable_scale_still_withholds_the_bootstrap() -> None: + """Fail-closed control: selectable must not mean an interval was invented.""" + t, y = _underflowing_scale_case() + with pytest.warns(RuntimeWarning): + with pytest.raises(Exception) as excinfo: + parametric_bootstrap( + _constant_model, t, y, np.array([0.0]), B=4, + rng=np.random.default_rng(0), + ) + assert excinfo.type is RuntimeError, excinfo.type + assert "no representable residual scale" in str(excinfo.value) + + +def test_ordinary_fit_is_untouched_by_the_scale_repair() -> None: + """Over-correction control: a normal fit keeps a finite positive sigma.""" + rng = np.random.default_rng(20260903) + t = np.linspace(0.0, 4.0, 64) + p = np.array([1.0, 0.8]) + y = M0(t, p) + 1.0e-3 * rng.standard_normal(t.size) + fit = fit_gls_ar1(M0, t, y, p, n_iters=1) + assert not fit.scale_unavailable + assert np.isfinite(fit.sigma) and fit.sigma > 0.0 + + +# -------------------------------------------------------------------------- +# PR #147 round-2 review +# -------------------------------------------------------------------------- + + +def test_scale_unavailable_reaches_fitresult() -> None: + """The reason a CI was withheld must survive into the persisted report. + + ``bca_ci_beta = (nan, nan)`` is what EVERY bootstrap or jackknife failure + produces. Without this flag on ``FitResult`` an unrepresentable residual + scale -- a fit that is otherwise successful and selectable -- was + indistinguishable from a non-converged resample once the warning stream + was gone. + """ + t, y = _underflowing_scale_case() + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + fit_result, _ = relaxation_mod._fit_with_model("M0", t, y) + assert any( + "scale is not representable" in str(w.message) for w in caught + ), [str(w.message) for w in caught] + assert fit_result.scale_unavailable + # The fit stays a candidate: the flag records why the INTERVAL is missing, + # it does not withhold the estimate (round-1 review contract). + assert fit_result.success + assert math.isfinite(fit_result.log_likelihood) + # The field is serialised, so an audit artefact carries the reason. + assert _to_jsonable(fit_result)["scale_unavailable"] is True + + +def test_ordinary_fit_reaches_fitresult_without_the_flag() -> None: + """Over-correction control: a normal fit must not be labelled scaleless.""" + rng = np.random.default_rng(20260903) + t = np.linspace(0.0, 4.0, 64) + y = M0(t, np.array([1.0, 0.8])) + 1.0e-3 * rng.standard_normal(t.size) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + fit_result, _ = relaxation_mod._fit_with_model("M0", t, y) + assert not fit_result.scale_unavailable + + +def test_explicit_sigma_likelihood_uses_the_full_representable_rss_range() -> None: + """``-0.5 * RSS`` is representable up to an RSS of ``2 * float64.max``. + + The guard compared the standardised RSS against ``float64.max`` and + returned ``-inf`` for the octave above it, although the value the formula + actually forms -- the HALVED RSS -- is perfectly finite there. Dropping a + model from selection on that boundary is the absolute-scale dependence + issue #135 exists to remove. + """ + residual = np.array([1.4e154]) + log_lik = gaussian_log_likelihood(residual, sigma=1.0) + assert math.isfinite(log_lik), log_lik + expected = -0.5 * math.log(2.0 * math.pi) - math.exp( + 2.0 * math.log(1.4e154) - math.log(2.0) + ) + assert log_lik == pytest.approx(expected, rel=1e-12) + + # POSITIVE CONTROL that the fixture really sits in the disputed octave: + # the raw standardised RSS is NOT representable, only its half is. + log_rss = scaled_log_sum_squares(residual) + log_max = math.log(np.finfo(float).max) + assert log_max < log_rss <= log_max + math.log(2.0), log_rss + + +def test_a_genuinely_unrepresentable_half_rss_still_returns_minus_inf() -> None: + """Fail-closed control: the bound moved by one octave, it did not vanish. + + The outcome is captured rather than asserted inline. Without the guard the + call raises ``OverflowError`` from ``math.exp``, and a test that dies of an + exception is indistinguishable from an incidental crash in a mutation run -- + the death must be attributable to an ASSERTION for the proof to count. + """ + residual = np.array([1.4e308]) + log_rss = scaled_log_sum_squares(residual) + assert log_rss > math.log(np.finfo(float).max) + math.log(2.0) + + outcome: object + try: + outcome = gaussian_log_likelihood(residual, sigma=1.0) + except Exception as exc: + outcome = exc + assert not isinstance(outcome, BaseException), ( + f"the likelihood raised {type(outcome).__name__} " + f"({outcome}) instead of returning -inf" + ) + assert outcome == float("-inf") + + +def test_values_inside_the_old_range_are_bit_identical() -> None: + """Widening a range must not perturb the values already inside it. + + The halving stays ``0.5 * exp(...)`` below the disputed octave, so every + likelihood the previous implementation could compute is reproduced exactly + rather than to within an ulp. + """ + for scale in (1.0e-30, 1.0, 3.7, 1.0e30, 1.0e120): + residuals = scale * np.array([0.3, -1.1, 2.0, 0.0]) + log_rss = scaled_log_sum_squares(residuals) + log_standardised = log_rss # sigma = 1 + assert log_standardised <= math.log(np.finfo(float).max) + n = residuals.size + expected = ( + -0.5 * n * math.log(2.0 * math.pi) - 0.5 * math.exp(log_standardised) + ) + assert gaussian_log_likelihood(residuals, sigma=1.0) == expected + + +# -------------------------------------------------------------------------- +# PR #147 round-3 review +# -------------------------------------------------------------------------- + + +def test_the_halved_boundary_is_decided_on_the_materialised_value() -> None: + """Round 2 moved the bound; round 3 removes the bound. + + The round-2 guard compared ``log_rss`` against ``log_max + log(2)``, two + separately rounded logarithms. On the reviewer's fixture -- ``sigma = 1`` + with a single residual ``sqrt(float_max) * sqrt(2)`` -- ``log_rss`` lands + exactly ONE ulp above that bound and the fit was dropped, although the + quantity the formula actually forms is 0.9999999999999999 of + ``float_max`` and entirely representable. + + The repair is not a wider bound. It is materialising the halved sum from + the residuals and letting float64 addition decide, so there is no rounded + threshold left to be one ulp wrong about. + """ + residual = np.array([math.sqrt(np.finfo(float).max) * math.sqrt(2.0)]) + + # PREMISE, asserted rather than assumed: the fixture is in the octave the + # round-2 bound rejected, and it misses that bound by a single ulp. + log_rss = scaled_log_sum_squares(residual) + log_max = math.log(np.finfo(float).max) + round2_bound = log_max + math.log(2.0) + assert log_rss > round2_bound, "fixture no longer trips the round-2 bound" + assert (log_rss - round2_bound) == pytest.approx( + math.ulp(log_max), rel=1e-9 + ), "the fixture no longer sits one ulp above the bound" + + log_lik = gaussian_log_likelihood(residual, sigma=1.0) + assert math.isfinite(log_lik), log_lik + + # The value is the one the formula defines, computed the direct way. + scaled = residual / 1.0 + half = float(np.sum((0.5 * scaled) * scaled)) + assert math.isfinite(half) + assert half == pytest.approx(0.9999999999999999 * np.finfo(float).max, rel=1e-15) + assert log_lik == pytest.approx(-0.5 * math.log(2.0 * math.pi) - half, rel=1e-15) + + +def test_the_log_space_route_could_not_have_reached_it() -> None: + """Why the fix is not "subtract log(2) and exponentiate". + + One ulp of a logarithm near 710 is a factor of ~1e-16 in the value, and + that is precisely the width of the decision being made here. Documented + as a live assertion because the obvious smaller repair looks correct and + is not: ``math.exp`` raises on the same input for which the direct + computation returns a finite number. + """ + residual = np.array([math.sqrt(np.finfo(float).max) * math.sqrt(2.0)]) + log_rss = scaled_log_sum_squares(residual) + + raised: BaseException | None = None + try: + math.exp(log_rss - math.log(2.0)) + except Exception as exc: + raised = exc + assert isinstance(raised, OverflowError), ( + "the log-space route no longer overflows; the round-3 rationale " + f"needs rereading (got {raised!r})" + ) + assert math.isfinite(gaussian_log_likelihood(residual, sigma=1.0)) + + +def test_a_half_rss_that_truly_overflows_is_still_dropped() -> None: + """Fail-closed control: removing the bound must not remove the refusal. + + Three shapes of genuine overflow, including one that only overflows in + the SUM -- no single term does -- which a per-term bound would have + missed. Captured and asserted rather than crashed into, so a mutation run + can attribute the death to an assertion. + """ + cases = { + "4x float_max, one residual": np.array( + [math.sqrt(np.finfo(float).max) * 2.0] + ), + "1e3x float_max, one residual": np.array( + [math.sqrt(np.finfo(float).max) * math.sqrt(1000.0)] + ), + "overflows only in the sum": np.full( + 64, math.sqrt(np.finfo(float).max) * math.sqrt(2.0) / 4.0 + ), + } + for label, residuals in cases.items(): + outcome: object + try: + outcome = gaussian_log_likelihood(residuals, sigma=1.0) + except Exception as exc: + outcome = exc + assert not isinstance(outcome, BaseException), ( + f"{label}: raised {type(outcome).__name__} ({outcome}) " + "instead of returning -inf" + ) + assert outcome == float("-inf"), label