diff --git a/CHANGELOG.md b/CHANGELOG.md index 9a6a9029..9c67dd79 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. - **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 diff --git a/CITATION.cff b/CITATION.cff index 18b5faca..202b5eb2 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/src/liouscope/_types.py b/src/liouscope/_types.py index f6893937..67ddbf22 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 d8433bc7..666ddddb 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 2130a6de..f0781649 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 b936ee56..1ec84e17 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 83bbeb87..135419b1 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: @@ -217,20 +232,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/linalg.py b/src/liouscope/numerics/linalg.py index bffbffca..56cfb143 100644 --- a/src/liouscope/numerics/linalg.py +++ b/src/liouscope/numerics/linalg.py @@ -29,6 +29,7 @@ ZERO_MODE_APOSTERIORI_MARGIN, ZERO_MODE_EPS_FACTOR, ) +from .norms import scaled_cancellation_ratio, scaled_euclidean_norm @dataclass(frozen=True, slots=True) @@ -299,23 +300,27 @@ def operator_zero_tolerance( def trace_preservation_defect(L_super: np.ndarray) -> tuple[float, float]: - """Return ``(||vec(I)^H L||_2, ||L||_F)`` for a vectorised generator. + """Return robust ``(||vec(I)^H L||_2, ||L||_F)`` for a generator. Every GKSL generator is trace preserving, which in column-stacking convention is the exact statement ``vec(I)^H L = 0``. The quotient of the two returned numbers is therefore a dimensionless, unit-invariant measure of how far ``L`` is from being a legal generator. - Both norms square their entries before summing, so a generator expressed - in small enough -- but still perfectly normal -- rate units loses its - reference scale to underflow: for ``diag([0, -1e-200, -2e-200, -3e-200])`` - every square falls below ``5e-324`` and the pair comes back ``(0.0, 0.0)``. - The applicability gate then reads ``0 > tp_rtol * tiny``, which is False, - so a demonstrably non-trace-preserving operator becomes certificate- - applicable -- while the SAME operator at ``1e-150`` is correctly refused. - A pure change of rate unit must not decide whether an operator is a legal - generator, so the norms are recomputed on a scaled copy when the reference - scale has been lost (round-23 review, PR #121). + Issue #130: both quantities are evidence, so their arithmetic must be + correct independently of the policy that consumes them. Ordinary + ``np.linalg.norm`` squares before summing and can therefore return the + false sentinels ``0.0`` (underflow) or ``inf`` (overflow) even when the + mathematical norm of finite input is representable. Both norms now use + :func:`scaled_euclidean_norm`, an xLASSQ-style power-of-two scaled + sum-of-squares primitive. It returns a finite norm whenever the true + float64 result is representable and returns infinity only when the input or + the true norm is genuinely non-representable. + + This deliberately separates measurement from certificate policy. A very + large but finite operator must not be rejected merely because a naive norm + overflowed; structural guards such as ``band_discriminates`` decide whether + the resulting zero-mode evidence actually discriminates. """ L_super = np.asarray(L_super) n = L_super.shape[0] @@ -323,95 +328,46 @@ def trace_preservation_defect(L_super: np.ndarray) -> tuple[float, float]: if d * d != n: return float("nan"), float("nan") vec_i = np.eye(d, dtype=complex).reshape(-1, order="F") - tp_row = vec_i.conj() @ L_super - defect = float(np.linalg.norm(tp_row)) - fro = float(np.linalg.norm(L_super, ord="fro")) - # The rescue is deliberately ONE-DIRECTIONAL: it repairs underflow and - # leaves overflow alone. ``fro = inf`` is not an accident to be worked - # around but the input to the round-21 refusal of a non-finite reference - # scale -- the guard that closed the round-20 counterexample (cancelling - # +-1e308 entries, spectrum {1,2,3,4}, certified as entirely stationary). - # Computing that Frobenius norm by scaling would make it finite again at - # ~1.4e308, readmit the operator, and silently reopen a hole whose repair - # CI had to prove across five interpreter versions. Underflow and overflow - # look symmetric here and are not: only one of them has a documented - # decision behind it. - if fro == 0.0 and L_super.size: - # ``max(|Re|, |Im|)`` rather than ``max|.|``: the modulus of a finite - # complex entry can itself overflow, and this quantity cannot. It is - # within sqrt(2) of the true maximum magnitude, which is all that is - # needed -- it only selects the exponent to shift by. - m = float(np.max(np.maximum(np.abs(np.real(L_super)), np.abs(np.imag(L_super))))) - if m > 0.0: - # Scaled by a POWER OF TWO, not by ``m`` itself. Two reasons, the - # second of which cost a regression before it was measured: - # - # 1. ``ldexp`` shifts the exponent and leaves every mantissa bit - # alone, so the round trip is EXACT -- dividing by ``m`` and - # multiplying back introduces two roundings into a number whose - # whole purpose is to be compared against a tolerance. - # 2. ``L / m`` is COMPLEX division, and for a subnormal divisor - # NumPy's algorithm forms an intermediate reciprocal that - # overflows: measured for ``diag([0, -1e-310, -2e-310, - # -3e-310])``, ``L / m`` came back ``[nan, -inf, -inf, -inf]`` - # and both norms became NaN. The round-21 guard then refused the - # operator -- fail-closed, but on a value nothing had measured, - # and a healthy GKSL generator at that scale was refused with it - # (``applicable`` went True -> False for valid input). Repairing - # a fail-open must not install a fail-closed defect in its place. - # - # ``frexp`` returns ``m = mantissa * 2**exp`` with the mantissa in - # [0.5, 1), so shifting by ``-exp`` puts every entry in a range - # where neither the squares below nor the reconstruction can leave - # the normal range. A genuinely zero operator has ``m == 0`` and - # keeps ``(0.0, 0.0)``, which is the correct answer, not a lost one. - exp = int(np.frexp(m)[1]) - scaled = np.ldexp(np.real(L_super), -exp) + 1j * np.ldexp( - np.imag(L_super), -exp - ) - defect = float( - np.ldexp(float(np.linalg.norm(vec_i.conj() @ scaled)), exp) - ) - fro = float(np.ldexp(float(np.linalg.norm(scaled, ord="fro")), exp)) - # ROUND-24 REVIEW (PR #121), finding B2. The rescue above is triggered by - # the DENOMINATOR alone, and the two norms do not share a scale. The - # NUMERATOR can underflow on its own while the operator around it is - # perfectly ordinary: add a representable ``1e-200`` trace-preservation - # violation to an O(1) generator and ``vec(I)^H L`` is - # ``[1e-200, 0, 0, 0]``, whose squares fall below ``5e-324`` -- so - # ``np.linalg.norm`` returns ``0.0`` while ``fro`` stays at 1.62 and the - # branch above never fires. With ``tp_rtol = 0`` (or any tolerance small - # enough) the gate then reads ``0.0 > 0.0``, which is False, and BOTH - # certificate APIs report ``applicable=True, trace_defect=0.0`` for an - # operator that is demonstrably not trace preserving. Measured before this - # repair. - # - # The condition is deliberately not ``defect == 0.0`` alone: a genuine GKSL - # generator has ``vec(I)^H L`` identically zero, and ``0.0`` is then the - # CORRECT answer rather than a lost one. What distinguishes the two is - # whether the source expression is itself zero, so that is what is asked -- - # "did a nonzero quantity collapse", not "is the result zero". - # - # The repair stays ONE-DIRECTIONAL in the same sense as the one above: it - # is reachable only from an exact zero, so no overflowing quantity is made - # finite and the round-21 refusal of a non-finite reference scale is - # untouched. - if defect == 0.0 and np.any(tp_row != 0.0): - m_row = float( - np.max(np.maximum(np.abs(np.real(tp_row)), np.abs(np.imag(tp_row)))) - ) - # ``m_row > 0`` follows from the ``np.any`` above: a vector with a - # nonzero entry has a nonzero largest component. Scaling by a power of - # two selected from that component -- not by the OPERATOR's scale, - # which is precisely the quantity this row does not share. - exp_row = int(np.frexp(m_row)[1]) - row_scaled = np.ldexp(np.real(tp_row), -exp_row) + 1j * np.ldexp( - np.imag(tp_row), -exp_row - ) - defect = float( - np.ldexp(float(np.linalg.norm(row_scaled)), exp_row) + defect_vector = vec_i.conj() @ L_super + return scaled_euclidean_norm(defect_vector), scaled_euclidean_norm(L_super) + + +def trace_preservation_componentwise_error(L_super: np.ndarray) -> float: + """Return the worst componentwise backward error of ``vec(I)^H L = 0``. + + The global ratio ``||vec(I)^H L|| / ||L||`` is necessary but insufficient + as an applicability test: entries of ``L`` that do not participate in a + particular trace equation can make the global operator norm enormous and + thereby hide an order-one violation in another column. Trace preservation + is instead a family of scalar equations, one for every superoperator + column ``j``: + + ``sum_i L[(i, i), j] = 0``. + + For each equation we therefore compute + ``abs(sum terms) / sum(abs(terms))`` and return the largest ratio. This is + the componentwise relative backward-error reading: every equation is + normalised only by the represented coefficients that can actually cancel + in that equation. :func:`scaled_cancellation_ratio` evaluates the ratio in + a power-of-two scaled domain, so the verdict is stable for subnormal and + near-overflow rate units. An exact all-zero equation contributes zero. + """ + arr = np.asarray(L_super) + if arr.ndim != 2 or arr.shape[0] != arr.shape[1] or arr.size == 0: + return float("nan") + if not np.all(np.isfinite(arr)): + return float("nan") + n = arr.shape[0] + d = int(round(np.sqrt(n))) + if d * d != n: + return float("nan") + trace_rows = np.arange(0, n, d + 1, dtype=int) + return float( + max( + (scaled_cancellation_ratio(arr[trace_rows, j]) for j in range(n)), + default=0.0, ) - return defect, fro + ) def underflow_safe_norm(v: np.ndarray) -> float: @@ -1123,6 +1079,7 @@ def certified_eigvals( L_super = require_finite_square_2d(L_super, name="L_super") L_c = np.asarray(L_super, dtype=complex) tp_defect, fro = trace_preservation_defect(L_c) + tp_componentwise = trace_preservation_componentwise_error(L_c) norm2 = float(np.linalg.norm(L_c, 2)) if L_c.size else 0.0 bound = rtol * float(np.finfo(float).eps) * norm2 @@ -1151,6 +1108,13 @@ def certified_eigvals( # RELATIVE to the operator's own scale is not a statement that # can be made about infinity, so it is refused here. or not np.isfinite(fro) + # ISSUE #130. The normwise ratio above can be diluted by enormous + # entries that do not participate in a violated trace equation. The + # componentwise gate normalises each column only by the diagonal-output + # terms that can cancel in that equation, then takes the worst column. + # Both views must pass before the exact zero-mode theorem is applicable. + or not np.isfinite(tp_componentwise) + or tp_componentwise > tp_rtol # ROUND-23 REVIEW (PR #121), same finding one line further down. # ``max(fro, tiny)`` floors the REFERENCE SCALE at a constant that # is larger than the operator whenever the operator is subnormal, @@ -1394,6 +1358,7 @@ def certified_eig( L_super = require_finite_square_2d(L_super, name="L_super") L_c = np.asarray(L_super, dtype=complex) tp_defect, fro = trace_preservation_defect(L_c) + tp_componentwise = trace_preservation_componentwise_error(L_c) norm2 = float(np.linalg.norm(L_c, 2)) if L_c.size else 0.0 bound = rtol * float(np.finfo(float).eps) * norm2 @@ -1422,6 +1387,10 @@ def certified_eig( # RELATIVE to the operator's own scale is not a statement that # can be made about infinity, so it is refused here. or not np.isfinite(fro) + # ISSUE #130. Keep the eigenvector ladder on the exact same + # applicability contract as the eigenvalue-only ladder above. + or not np.isfinite(tp_componentwise) + or tp_componentwise > tp_rtol # ROUND-23 REVIEW (PR #121), same finding one line further down. # ``max(fro, tiny)`` floors the REFERENCE SCALE at a constant that # is larger than the operator whenever the operator is subnormal, diff --git a/src/liouscope/numerics/norms.py b/src/liouscope/numerics/norms.py new file mode 100644 index 00000000..712973f9 --- /dev/null +++ b/src/liouscope/numerics/norms.py @@ -0,0 +1,146 @@ +"""Overflow/underflow-safe norm and cancellation primitives. + +The Euclidean norm implementation follows the same numerical contract as +LAPACK's xLASSQ family: scale values before squaring, accumulate in the scaled +domain, then restore the scale. LiouScope uses a power-of-two scale so the +rescaling itself does not introduce a reciprocal overflow for subnormal complex +values and does not round mantissas merely to choose numerical units. +""" + +from __future__ import annotations + +import math + +import numpy as np + + +def _finite_component_scale(values: np.ndarray) -> tuple[np.ndarray, np.ndarray, int] | None: + """Return real/imag parts shifted so their largest component is O(1). + + ``None`` denotes an exact all-zero array. NaN/inf are deliberately left to + public callers because their desired sentinel semantics differ by operation. + """ + arr = np.asarray(values) + real = np.asarray(np.real(arr), dtype=float) + imag = np.asarray(np.imag(arr), dtype=float) + component_max = float( + max( + float(np.max(np.abs(real))) if real.size else 0.0, + float(np.max(np.abs(imag))) if imag.size else 0.0, + ) + ) + if component_max == 0.0: + return None + exponent = int(np.frexp(component_max)[1]) + return np.ldexp(real, -exponent), np.ldexp(imag, -exponent), exponent + + +def scaled_euclidean_norm(values: np.ndarray) -> float: + """Return ``sqrt(sum(abs(values)**2))`` without spurious under/overflow. + + This is the vector 2-norm and, when ``values`` is a matrix, its Frobenius + norm. For finite input it returns a finite number whenever the true norm is + representable as ``float64``. It returns ``inf`` only when an infinite + input component or the mathematical norm itself is not representable, and + propagates ``NaN`` rather than fabricating evidence from corrupted input. + + Complex values are treated as their real and imaginary components, so the + calculation is mathematically ``sqrt(sum(re**2 + im**2))``. We avoid + computing ``abs(z)`` before selecting the scale because a finite complex + value near the top of the floating-point range can overflow in the modulus. + """ + arr = np.asarray(values) + if arr.size == 0: + return 0.0 + + 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 0.0 + scaled_real, scaled_imag, exponent = scaled + sumsq = float( + np.sum(scaled_real * scaled_real, dtype=float) + + np.sum(scaled_imag * scaled_imag, dtype=float) + ) + scaled_norm = float(np.sqrt(sumsq)) + with np.errstate(over="ignore", under="ignore"): + 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. + + This is a componentwise relative backward-error primitive for a scalar + linear equation assembled by cancellation. It is in ``[0, 1]`` for finite + input, is exactly scale invariant under representable non-zero rescaling, + and returns zero for an all-zero equation. + + Power-of-two scaling prevents both overflow and subnormal underflow before + summation. ``math.fsum`` is used separately on real and imaginary parts so + the numerator measures the represented equation rather than ordinary + left-to-right summation noise. The denominator is accumulated from stable + ``hypot`` magnitudes in the same scaled domain, so unrelated large entries + outside the equation cannot dilute its defect. + """ + arr = np.asarray(values) + if arr.size == 0: + return 0.0 + real = np.asarray(np.real(arr), dtype=float) + imag = np.asarray(np.imag(arr), dtype=float) + if np.any(~np.isfinite(real)) or np.any(~np.isfinite(imag)): + return float("nan") + + scaled = _finite_component_scale(arr) + if scaled is None: + return 0.0 + scaled_real, scaled_imag, _ = scaled + sum_real = math.fsum(float(x) for x in scaled_real.ravel()) + sum_imag = math.fsum(float(x) for x in scaled_imag.ravel()) + numerator = math.hypot(sum_real, sum_imag) + denominator = math.fsum( + math.hypot(float(re), float(im)) + for re, im in zip(scaled_real.ravel(), scaled_imag.ravel(), strict=True) + ) + if denominator == 0.0: + return 0.0 + return float(min(1.0, numerator / denominator)) diff --git a/tests/test_issue130_componentwise_tp.py b/tests/test_issue130_componentwise_tp.py new file mode 100644 index 00000000..8df6ff80 --- /dev/null +++ b/tests/test_issue130_componentwise_tp.py @@ -0,0 +1,62 @@ +"""Issue #130: componentwise trace-preservation backward-error regressions.""" + +from __future__ import annotations + +import numpy as np +import pytest + +from liouscope.numerics.linalg import trace_preservation_componentwise_error +from liouscope.numerics.norms import scaled_cancellation_ratio + + +def test_scaled_cancellation_ratio_exact_cancellation_at_extreme_scales() -> None: + for scale in (1.0e-300, 1.0, 1.0e300): + values = np.array([scale, -scale], dtype=complex) + assert scaled_cancellation_ratio(values) == 0.0 + + +def test_scaled_cancellation_ratio_detects_uncancelled_equation() -> None: + for scale in (1.0e-300, 1.0, 1.0e300): + values = np.array([scale, 0.0], dtype=complex) + assert scaled_cancellation_ratio(values) == pytest.approx(1.0) + + +def test_scaled_cancellation_ratio_is_scale_invariant_for_partial_cancellation() -> None: + reference = np.array([2.0 + 3.0j, -2.0 - 2.0j, 0.25 - 0.5j]) + expected = scaled_cancellation_ratio(reference) + assert 0.0 < expected < 1.0 + for scale in (1.0e-250, 1.0e-100, 1.0, 1.0e100, 1.0e250): + got = scaled_cancellation_ratio(reference * scale) + assert got == pytest.approx(expected, rel=2.0e-15, abs=0.0) + + +def test_tp_componentwise_error_ignores_unrelated_huge_entries() -> None: + # d=2 => trace-output rows are vectorised matrix rows 0 and 3. + op = np.zeros((4, 4), dtype=complex) + op[1, 0] = 1.0e300 # huge but irrelevant to the trace equation + op[0, 2] = 1.0 + assert trace_preservation_componentwise_error(op) == pytest.approx(1.0) + + +def test_tp_componentwise_error_accepts_exact_trace_cancellation_at_huge_scale() -> None: + op = np.zeros((4, 4), dtype=complex) + op[0, 1] = 1.0e300 + op[3, 1] = -1.0e300 + assert trace_preservation_componentwise_error(op) == 0.0 + + +def test_tp_componentwise_error_is_rate_unit_invariant() -> None: + op = np.zeros((4, 4), dtype=complex) + op[0, 0] = 3.0 + op[3, 0] = -2.0 + expected = 1.0 / 5.0 + for scale in (1.0e-300, 1.0e-200, 1.0, 1.0e200, 1.0e300): + got = trace_preservation_componentwise_error(op * scale) + assert got == pytest.approx(expected, rel=2.0e-15, abs=0.0) + + +def test_tp_componentwise_error_fails_closed_on_bad_shape_and_nonfinite() -> None: + assert np.isnan(trace_preservation_componentwise_error(np.ones((3, 3)))) + bad = np.zeros((4, 4), dtype=complex) + bad[0, 0] = np.nan + assert np.isnan(trace_preservation_componentwise_error(bad)) diff --git a/tests/test_issue135_likelihood_scale.py b/tests/test_issue135_likelihood_scale.py new file mode 100644 index 00000000..8231541a --- /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 diff --git a/tests/test_pr121_review_round20.py b/tests/test_pr121_review_round20.py index 8ed36204..d09c3489 100644 --- a/tests/test_pr121_review_round20.py +++ b/tests/test_pr121_review_round20.py @@ -139,11 +139,12 @@ def test_a_resolved_certificate_still_publishes_d3_d4_and_the_flag(monkeypatch) def _overflowing_norm_operator() -> np.ndarray: """The reviewer's counterexample, entry for entry. - Finite 4x4, cancelling ``+-1e308`` in the ``vec(I)^H L`` combination, so - the absolute trace defect is ``sqrt(17) = 4.12`` while ``||L||_F`` - overflows to ``inf``. Eigenvalues are exactly ``{1, 2, 3, 4}``: column 0 - and column 2 each carry a single diagonal entry, and the remaining 2x2 on - indices ``{1, 3}`` is lower triangular. + Finite 4x4, cancelling ``+-1e308`` in one trace equation. A naive + sum-of-squares Frobenius norm overflows, but the TRUE norm is the finite + representable value ``sqrt(2)*1e308``. The absolute trace defect is + ``sqrt(17) = 4.12`` and the componentwise trace-equation error is 1: other + enormous entries are not allowed to dilute the violated columns. + Eigenvalues are exactly ``{1, 2, 3, 4}``. """ L = np.diag([1.0, 2.0, 3.0, 4.0]).astype(complex) L[0, 1] = 1.0e308 @@ -152,14 +153,13 @@ def _overflowing_norm_operator() -> np.ndarray: def _nilpotent_wide_band_operator() -> np.ndarray: - """The same failure WITHOUT any overflow, to show the class is general. + """The same band failure WITHOUT a trace-preservation failure. Strictly upper triangular with ``1e10`` entries confined to row 1, so ``vec(I)^H L = 0`` EXACTLY (trace preserving, defect 0.0), ``||L||_F`` is a perfectly ordinary ``1.4e10``, and every eigenvalue is exactly zero. The band ``rtol * eps * ||L||_2`` is then ``3.1e-3`` and swallows all four - -- a nilpotent generator certified as a four-dimensional stationary - manifold, on finite, well-scaled arithmetic. + -- a nilpotent generator whose band has discriminated nothing. """ L = np.zeros((4, 4), dtype=complex) L[1, 2] = 1.0e10 @@ -167,34 +167,36 @@ def _nilpotent_wide_band_operator() -> np.ndarray: return L -def test_the_overflowing_operator_is_the_state_the_reviewer_described() -> None: - """Pins the preconditions, so a numpy change cannot hollow out the test.""" - from liouscope.numerics.linalg import trace_preservation_defect +def test_large_operator_measurement_and_local_tp_error_are_both_real() -> None: + """Pin the corrected arithmetic and the independent policy precondition.""" + from liouscope.numerics.linalg import ( + trace_preservation_componentwise_error, + trace_preservation_defect, + ) L = _overflowing_norm_operator() - assert np.all(np.isfinite(L)), "the input must be FINITE for this to be the bug" - with np.errstate(over="ignore"): - defect, fro = trace_preservation_defect(L) + assert np.all(np.isfinite(L)), "the input must be FINITE for this to be the case" + defect, fro = trace_preservation_defect(L) + local_error = trace_preservation_componentwise_error(L) + assert defect == pytest.approx(np.sqrt(17.0)), "a real, order-one trace defect" - assert not np.isfinite(fro), ( - "the whole failure is that ``defect > tp_rtol * fro`` compares " - "against infinity; if ||L||_F stops overflowing the test is hollow" + assert np.isfinite(fro), "a representable norm must not be turned into infinity" + assert fro == pytest.approx(np.sqrt(2.0) * 1.0e308, rel=2.0e-15) + assert local_error == pytest.approx(1.0), ( + "at least one trace equation has no cancellation at all; unrelated " + "1e308 entries must not make that violation look relatively small" ) assert np.allclose(np.sort(np.linalg.eigvals(L).real), [1.0, 2.0, 3.0, 4.0]) @pytest.mark.parametrize( - "make", [_overflowing_norm_operator, _nilpotent_wide_band_operator] + ("make", "expected_applicable"), + [(_overflowing_norm_operator, False), (_nilpotent_wide_band_operator, True)], ) -def test_a_band_that_swallows_the_spectrum_does_not_certify(make) -> None: - """THE regression, on both ladders. - - Before this guard the certificate read ``certified=True, resolved=True`` - with ``zero_mode_count = 4`` -- every eigenvalue declared stationary, - which for a d=2 system asserts that the entire operator space is a steady - state. ``resolved=True`` is the damaging half: it tells the spectral layer - that D1/D3/D4 may be published. - """ +def test_a_non_tp_input_or_a_nondiscriminating_band_does_not_certify( + make, expected_applicable: bool +) -> None: + """Both ladders fail closed, but the two counterexamples fail differently.""" from liouscope.numerics.linalg import certified_eig, certified_eigvals L = make() @@ -203,64 +205,35 @@ def test_a_band_that_swallows_the_spectrum_does_not_certify(make) -> None: _, cert_vals = certified_eigvals(L) _, cert_vec = certified_eig(L) - # ROUND 21. Both operators must be refused, but for DIFFERENT reasons, and - # this test names which -- asserting only the shared verdict would let - # either mechanism rot unnoticed behind the other. - # - # * finite scale (the nilpotent twin): the certificate legitimately - # applies, and the BAND GUARD catches it -- the band swallowed the whole - # spectrum, so it discriminated nothing. - # * non-finite scale (the overflowing operator): there is no measurable - # scale to be trace preserving *relative to*. ``defect <= tp_rtol * fro`` - # against an infinite ``fro`` admits everything, so the band drawn from - # it was never grounded in a measurement. Which modes it then swallows - # is decided by the LAPACK build -- 4 of 4 on Windows, 2 of 4 on the CI - # runners, and the 2-of-4 case slipped past a guard that only asks - # whether the band took *everything*. The refusal now happens one level - # earlier, where it is platform-independent. - with np.errstate(over="ignore"): - finite_scale = bool(np.isfinite(np.linalg.norm(L, ord="fro"))) - for cert in (cert_vals, cert_vec): assert cert.certified is False assert cert.resolved is False - if finite_scale: - assert cert.applicable is True - assert cert.zero_mode_count == 0 - else: - assert cert.applicable is False, ( - "an operator whose own scale overflows cannot be judged " - "trace preserving relative to that scale" - ) - - -def test_a_non_finite_scale_is_refused_before_any_band_is_drawn() -> None: - """ROUND 21: the fail-closed must key on the SCALE, not on the outcome. - - ``trace_preservation_defect`` returns ``(defect, ||L||_F)`` and the - applicability test is ``defect <= tp_rtol * ||L||_F``. A non-finite - ``defect`` was already refused; a non-finite ``||L||_F`` was not, and that - asymmetry is the whole bug -- a comparison against infinity is vacuously - true, so every operator passed and the band that followed was drawn from a - scale nobody measured. - - This pins the invariant directly rather than through a downstream verdict, - so it cannot be satisfied by accident the way the band guard was. + assert cert.applicable is expected_applicable + assert cert.zero_mode_count == 0 + + +def test_componentwise_tp_gate_refuses_local_defect_before_zero_mode_claim() -> None: + """The #130 regression: finite global norm must not create a fail-open. + + The true Frobenius norm is ~1.4e308, so a GLOBAL normwise quotient makes + ``sqrt(17)/||L||`` microscopic. Trace preservation, however, consists of + separate cancellation equations. Several columns here have order-one + equation error, which invalidates the exact-zero-mode theorem regardless + of unrelated entries elsewhere in the operator. """ from liouscope.numerics.linalg import ( certified_eigvals, + trace_preservation_componentwise_error, trace_preservation_defect, ) L = _overflowing_norm_operator() - with np.errstate(over="ignore"): - defect, fro = trace_preservation_defect(L) + defect, fro = trace_preservation_defect(L) + local_error = trace_preservation_componentwise_error(L) - assert np.isfinite(defect), "the defect itself is finite -- that is the trap" - assert not np.isfinite(fro), ( - "if the Frobenius norm stops overflowing this test is hollow and the " - "asymmetry it guards no longer exists" - ) + assert np.isfinite(defect) + assert np.isfinite(fro) + assert local_error == pytest.approx(1.0) with np.errstate(over="ignore"), warnings.catch_warnings(): warnings.simplefilter("ignore") diff --git a/tests/test_pr121_review_round23.py b/tests/test_pr121_review_round23.py index 4f1c67e5..bc48aa28 100644 --- a/tests/test_pr121_review_round23.py +++ b/tests/test_pr121_review_round23.py @@ -20,10 +20,11 @@ sit below the eigensolver's backward error. Every guard here is paired with a positive control: a guard that refuses -everything discriminates nothing and must not be able to pass this file. The -overflow twin of finding 12 has its own test for the opposite reason -- it -pins a refusal that must NOT be repaired, because the round-21 guard depends -on it. +everything discriminates nothing and must not be able to pass this file. Issue +#130 extends finding 12 symmetrically: mathematically representable large norms +must stay finite, while trace-preservation applicability also uses a separate +componentwise backward-error check so unrelated huge entries cannot dilute a +violated trace equation. """ from __future__ import annotations @@ -39,6 +40,7 @@ from liouscope.numerics.linalg import ( certified_eig, certified_eigvals, + trace_preservation_componentwise_error, trace_preservation_defect, ) @@ -57,7 +59,7 @@ def _non_tp_operator() -> np.ndarray: # --------------------------------------------------------------------------- -# Finding 12: use underflow-safe norms for the TP gate. +# Finding 12: use underflow/overflow-safe norms for the TP gate. # --------------------------------------------------------------------------- @@ -110,35 +112,34 @@ def test_a_healthy_generator_in_the_same_rate_units_still_certifies() -> None: assert cert.certified, f"a GKSL generator at scale {scale:g} was not certified" -def test_the_overflow_direction_is_deliberately_not_rescued() -> None: - """The round-21 refusal must SURVIVE the round-23 repair. +def test_representable_large_norm_is_measured_then_rejected_structurally() -> None: + """Correct measurement must not create a fail-open applicability claim. - The two directions look symmetric and are not. ``||L||_F = inf`` for the - round-20 counterexample (cancelling +-1e308 entries, spectrum {1,2,3,4}) is - the input to the round-21 guard, which refuses a non-finite reference - scale. Its true Frobenius norm is about 1.4e308 and therefore perfectly - representable, so a scaled computation would hand back a finite number, - readmit the operator and reopen a hole that took five interpreter versions - of CI to close. This test goes red the moment someone makes the rescue - symmetric. + The round-20 counterexample has two cancelling ``+-1e308`` entries. Its + true Frobenius norm is about ``sqrt(2)*1e308`` and is therefore finite in + float64, even though a naive sum of squares overflows. That finite global + scale cannot prove trace preservation: several trace equations have no + cancellation at all, so their componentwise backward error is 1. The exact + zero-mode theorem must therefore be inapplicable before the band is used. """ - # The round-20 operator, entry for entry (see - # tests/test_pr121_review_round20.py::_overflowing_norm_operator). op = np.diag([1.0, 2.0, 3.0, 4.0]).astype(complex) op[0, 1] = 1.0e308 op[3, 1] = -1.0e308 assert np.all(np.isfinite(op)), "the input must be FINITE for this to be the case" - with np.errstate(over="ignore"): - defect, fro = trace_preservation_defect(op) + + defect, fro = trace_preservation_defect(op) + local_error = trace_preservation_componentwise_error(op) assert defect == pytest.approx(np.sqrt(17.0)), "a real, order-one trace defect" - assert not np.isfinite(fro), ( - "the overflowing reference scale was made finite; the round-21 guard " - "that refuses a non-finite scale now has nothing to refuse" - ) + assert np.isfinite(fro), "a representable Frobenius norm must not overflow" + assert fro == pytest.approx(np.sqrt(2.0) * 1.0e308, rel=2.0e-15) + assert local_error == pytest.approx(1.0) + with np.errstate(over="ignore"), warnings.catch_warnings(): warnings.simplefilter("ignore") cert = certified_eigvals(op)[-1] - assert not cert.applicable and not cert.certified + assert not cert.applicable + assert not cert.certified + assert not cert.resolved def test_the_repair_reaches_into_the_subnormal_range() -> None: diff --git a/tests/test_pr121_review_round24.py b/tests/test_pr121_review_round24.py index c6d8bd66..fc0e5001 100644 --- a/tests/test_pr121_review_round24.py +++ b/tests/test_pr121_review_round24.py @@ -228,18 +228,26 @@ def test_a_genuine_generator_still_reports_an_exactly_zero_defect() -> None: assert cert.applicable is True -def test_the_one_directional_overflow_policy_is_untouched() -> None: - """The round-21 decision must survive the round-24 repair. - - ``fro = inf`` is the documented INPUT to the non-finite-reference-scale - refusal that closed the round-20 counterexample. The numerator repair is - reachable only from an exact zero, so it can neither observe nor undo an - overflow -- this pins that, so a later widening of the trigger cannot - quietly readmit the operator. +def test_the_round_20_counterexample_is_still_refused() -> None: + """The round-21 DECISION survives; its former MECHANISM does not. + + Round 21 refused the round-20 counterexample because ``||L||_F`` came back + ``inf``. That sentinel was never the physics -- it was an artefact of + squaring before summing. The true Frobenius norm of this operator is + ``sqrt(2) * 1e308 ~ 1.414e308``, comfortably inside float64, verified + against exact rational arithmetic: ``np.linalg.norm`` was simply WRONG + here, and ``scaled_euclidean_norm`` (issue #130) reports the representable + value to within one ulp. + + So this test pins the property that had to be protected -- the operator is + REFUSED -- instead of the accident that used to deliver it. The refusal now + rests on ``trace_preservation_componentwise_error == 1``: several trace + equations have no cancellation at all, so no global scale can certify + trace preservation. That is a stronger argument than an overflow, because + it does not evaporate when the arithmetic improves. """ # The round-20 counterexample, entry for entry: cancelling +-1e308 in the - # ``vec(I)^H L`` combination, so the defect is sqrt(17) while ||L||_F - # overflows. + # ``vec(I)^H L`` combination, so the defect is sqrt(17). L = np.diag([1.0, 2.0, 3.0, 4.0]).astype(complex) L[0, 1] = 1.0e308 L[3, 1] = -1.0e308 @@ -247,9 +255,18 @@ def test_the_one_directional_overflow_policy_is_untouched() -> None: with np.errstate(over="ignore"): defect, fro = trace_preservation_defect(L) assert np.isfinite(defect), "the defect itself is finite -- that is the trap" - assert not np.isfinite(fro), ( - "if the Frobenius norm stops overflowing, the round-21 refusal this " - "test guards no longer has an input and the test is hollow" + assert defect == pytest.approx(np.sqrt(17.0)) + # The measurement is now CORRECT rather than saturated. + assert np.isfinite(fro), "the true norm is representable; inf was a defect" + assert fro == pytest.approx(np.sqrt(2.0) * 1.0e308, rel=2.0e-15) + # THE PROTECTED PROPERTY. If this ever flips to True, the round-20 hole is + # open again -- regardless of how the norms are computed. + with np.errstate(over="ignore"), warnings.catch_warnings(): + warnings.simplefilter("ignore") + cert = certified_eigvals(L)[-1] + assert cert.applicable is False, ( + "the round-20 counterexample must stay inapplicable; it is refused by " + "the componentwise backward error, not by an overflowing norm" ) diff --git a/tests/test_scaled_norms.py b/tests/test_scaled_norms.py new file mode 100644 index 00000000..90eba316 --- /dev/null +++ b/tests/test_scaled_norms.py @@ -0,0 +1,68 @@ +"""Numerical-contract tests for the scaled sum-of-squares primitive (#130).""" + +from __future__ import annotations + +import numpy as np +import pytest + +from liouscope.numerics.norms import scaled_euclidean_norm + + +def test_tiny_finite_values_do_not_underflow_to_zero() -> None: + values = np.array([1.0e-200, -2.0e-200, 3.0e-200]) + got = scaled_euclidean_norm(values) + + assert got > 0.0 + # Compare after dividing by the known input scale: pytest.approx's default + # absolute tolerance would otherwise dwarf a 1e-200 expected value and let + # a wildly inaccurate non-zero result pass. + assert got / 1.0e-200 == pytest.approx(np.sqrt(14.0), rel=2.0e-15) + + +def test_subnormal_values_keep_a_nonzero_norm() -> None: + values = np.array([1.0e-320, -2.0e-320, 3.0e-320]) + got = scaled_euclidean_norm(values) + + assert got > 0.0 + # Subnormal quantisation limits the attainable relative accuracy here, so + # compare the dimensionless ratio rather than using an absolute tolerance. + assert got / 1.0e-320 == pytest.approx(np.sqrt(14.0), rel=2.0e-4) + + +def test_large_finite_values_stay_finite_when_true_norm_is_representable() -> None: + values = np.array([1.0e308, -1.0e308, 1.0e308]) + got = scaled_euclidean_norm(values) + + assert np.isfinite(got) + assert got == pytest.approx(np.sqrt(3.0) * 1.0e308, rel=2.0e-15) + + +def test_nonrepresentable_true_norm_fails_as_infinity_not_a_false_finite_value() -> None: + values = np.array([1.0e308, -1.0e308, 1.0e308, -1.0e308]) + assert scaled_euclidean_norm(values) == float("inf") + + +def test_complex_components_follow_the_same_scaled_contract() -> None: + got = scaled_euclidean_norm(np.array([1.0e308 + 1.0e308j])) + + assert np.isfinite(got) + assert got == pytest.approx(np.sqrt(2.0) * 1.0e308, rel=2.0e-15) + + +@pytest.mark.parametrize("scale", [1.0e-250, 1.0e-100, 1.0, 1.0e100, 1.0e250]) +def test_norm_is_scale_equivariant_over_representable_ranges(scale: float) -> None: + base = np.array([1.0, -2.0, 3.0, -4.0]) * 1.0e-10 + reference = scaled_euclidean_norm(base) + got = scaled_euclidean_norm(scale * base) + expected = abs(scale) * reference + + assert got > 0.0 + # Compare the ratio to one so tiny expected values cannot pass through + # pytest.approx's default absolute tolerance. + assert got / expected == pytest.approx(1.0, rel=3.0e-15) + + +def test_zero_nan_and_inf_semantics_are_explicit() -> None: + assert scaled_euclidean_norm(np.zeros(4)) == 0.0 + assert np.isnan(scaled_euclidean_norm(np.array([1.0, np.nan]))) + assert scaled_euclidean_norm(np.array([1.0, np.inf])) == float("inf")