fix(#122): repair the two-scale relaxation window instead of disclosing it

CITATION.cff

@@ -84,6 +84,56 @@ 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 DEFAULT relaxation time grid spans a fixed
+# number of e-foldings of the slowest mode, [0, 10 / Delta], instead of the
+# absolute window [0, 10] (PR #115). This CHANGES NUMERICAL RESULTS -- for every
+# run that does not pass an explicit t_grid and whose gap is not Delta = 1 --
+# and must be described as a correction, not an addition: the previously
+# reported rates were measured on a window whose relation to the dynamics
+# depended on the caller's unit of time, so on an amplitude-damped qubit under
+# the pure rescale L -> cL the fitted beta_D drifted 22%, the D17 linear rate
+# was wrong by a factor ~430 at c = 1e-6, and four different A-classes were
+# reported for one system. It must NOT be described as making the fitted
+# relaxation rates invariant under an arbitrary change of rate units: the
+# least-squares solver's own convergence controls remain unit-dependent (open in
+# #111), and the mechanism verdict still gates on the rate-dimensioned
+# henrici_eta threshold (open in #101). The claim is bounded to what was
+# measured -- beta_D and beta_D_linear track the rescaling to <=1.2e-3 relative
+# over c in [1e-6, 1e6] with a stable AICc winner, and the A-class is stable
+# over c in [1e-6, 1] -- on the two qubit systems tested.
+# Also pending for the next cut: on a system whose rates are separated by more
+# than about eightfold the default window is now TWO-SCALE rather than uniform,
+# and the GLS residual model switches with it from discrete AR(1) to
+# continuous-time CAR(1) (per-step exp(-theta dt_k)). This CHANGES NUMERICAL
+# RESULTS for exactly those runs and is a correction, not an addition: on two
+# independent damped qubits with rates 1e-6 and 1 the uniform window's winning
+# M2 fit reported a second rate of 2.17e-05 for a component that was never
+# sampled (true value 1), and the D17 linear rate sat 0.58 relative from the
+# certified gap; the two-scale window recovers 1.10 and 0.021 respectively.
+# It must NOT be described as resolving arbitrary multiscale dynamics: the two
+# segments resolve the FASTEST and the SLOWEST mode, and an intermediate
+# timescale can still fall between the coarse late samples -- measured 0.0019
+# samples per e-folding for the middle mode of a 1e-6 / 1e-3 / 1 system, where
+# the UnderResolvedTransientWarning correctly still fires. Nor is the
+# small-sample bias correction of the AR(1) path carried over to CAR(1): every
+# way of anchoring it to a single step of a multi-decade grid was measured to
+# destroy the estimator's discrimination, so it is disclosed rather than
+# applied.
+# Also pending for the next cut: where the zero-mode certificate is applicable
+# but NOT resolved, the spectral layer now withholds D3 (oscillating_gap) and
+# D4 (spectral_spread) as NaN and reports has_complex_pairs as None, under the
+# same predicate that already withheld D1 (round-17 review, PR #127). This
+# CHANGES REPORTED RESULTS for exactly those runs and is a correction, not an
+# addition: the three values were previously read off the candidate spectrum
+# that the layer's own warning declares unreliable (measured on the issue-#113
+# stiff fixture: D3 = 0.0 and has_complex_pairs = False, i.e. an assertion that
+# no oscillatory separation exists, taken from a spectrum whose slow modes lie
+# below the eigensolver's backward error). Consumers see the A8 oscillatory
+# rung become UNEVALUABLE instead of NOT_SUPPORTED for those runs. It must NOT
+# be described as making the spectral layer reliable on stiff generators: the
+# unresolved case is still unresolved, only no longer reported as measured.
+# D2/D2b are unaffected -- they are computed from the operator and the steady
+# state rather than from the candidate spectrum.
 # 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
@@ -101,6 +151,18 @@ abstract: >-
 # 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.
+# Also pending for the next cut: the Hermiticity validation of the Hamiltonian
+# measures its defect against the scale of the GENERATOR, not of H alone (PR
+# #127): max(gauge-fixed max|H0|, max|sum_k gamma_k L0_k^dag L0_k| / 2), read in
+# the canonical Lindblad gauge (traceless jump operators, compensated H0). As
+# with issue #109 this changes which inputs are ACCEPTED, not the numerics of
+# accepted ones: a numerically pure-gauge Hamiltonian with a dissipator is now
+# accepted, while an identity-dominated H such as [[1e308, 1], [0, 1e308]], a
+# null dissipator c*I and the Lindblad gauge L -> L + c*I no longer excuse a
+# Hermiticity defect. Without jump operators the verdict is main's gauge-fixed
+# one, and an H whose gauge-fixed scale is not finite is refused. It must NOT be
+# described as settling whether a large physical dissipation may excuse a
+# defect of the coherent part: that question is open (cross-family review).
 keywords:
   - open quantum systems
   - Lindblad

docs/explanation/layers-and-taxonomy.md

@@ -175,3 +175,118 @@ not influence any verdict yet: the switch requires the preregistered
 calibration study and independent physics review specified in issue #101
 (slice C), including gapless-normal negative controls, before any threshold
 is chosen.
+
+## The relaxation window is measured in the system's own relaxation time
+
+Everything the relaxation layer reports — D5, D6, D7, the M0..M3b AICc
+comparison, `beta_D`, its BCa interval and the D17 gap-rate check — is fitted
+on a time grid. A decay rate has dimension `1/time`, so an *absolute* default
+window would be an unstated claim about the caller's unit of time.
+
+When `t_grid` is omitted, `diagnose()` therefore spans
+
+```text
+t in [0, RELAXATION_HORIZON / Delta],   80 uniform samples
+```
+
+with `Delta` the D1 gap and `RELAXATION_HORIZON = 10` — a fixed number of
+e-foldings of the slowest mode, which is the only window carried along by the
+rescaling `L → cL`. The fitted rates track that rescaling to ≤1.2e-3 relative
+over twelve decades of rate units
+(`tests/test_relaxation_grid_scale.py`). At `Delta = 1` the grid is
+bit-identical to the historical `linspace(0.0, 10.0, 80)`; when no decay scale
+is resolved (`Delta <= 0`) that historical window is used, since there is then
+no timescale to scale by.
+
+The grid is **uniform whenever a uniform grid suffices** — which is whenever
+the spread between the slowest and the fastest mode stays under about
+`(n_points − 1) / horizon = 7.9`. Past that the window switches to a two-scale
+grid (below), and the residual model switches with it.
+
+Which window produced a given run is recorded on the report, so it never has
+to be inferred — including the grid itself, which is the abscissa the exported
+D5/D6/D7 curves are sampled on:
+
+```python
+report.relaxation.t_grid_source
+#   "caller" | "gap_scaled" | "gap_scaled_multiscale" | "legacy_fixed"
+report.relaxation.t_grid_span
+report.relaxation.t_grid         # the sampling, not just its extent
+report.relaxation.residual_model
+#   what the fits were ACTUALLY whitened with, not what the grid asked for:
+#   "ar1"               uniform grid, discrete AR(1)
+#   "car1"              non-uniform grid, every fit whitened continuous-time
+#   "car1_fallback_ar1" CAR(1) theta failed on every fit -> AR(1) fallback
+#   "car1_mixed"        some fits CAR(1), some fallen back
+#   "car1_unavailable"  non-uniform grid and no fit succeeded
+```
+
+### What one uniform window cannot do — and what replaced it
+
+The two requirements pull against each other. The window must reach `~1/Δ` to
+see the slowest mode relax; the step must stay below `~1/r` to see a mode at
+rate `r` at all. Eighty uniform samples over ten e-foldings give
+
+```text
+samples_per_fast_efolding = min over modes of  1 / (r · blind_r)
+```
+
+where `blind_r` is the largest interval the grid leaves unsampled **while that
+mode still has amplitude** — i.e. the largest gap that starts before `1/r`,
+counting the lead-in `[0, t[0]]` as a gap. The lead-in matters because
+`diagnose()` accepts any non-negative start: on `linspace(100, 101, 101)` a
+rate-1 mode is sampled a hundred times per e-folding by its step and is still
+long gone by the first sample. On a uniform grid starting at zero this reduces
+to `1 / (r_max · dt) ≈ 7.9 · Δ / r_max`, so roughly an **eightfold** spread of
+timescales is the most one uniform grid can straddle.
+
+Widening the window does not help: an absolute window resolves the fast mode
+and misses the relaxation entirely, which is worse for the quantity this layer
+reports (measured `beta_D_linear` `3.5e4` relative from the true gap, against
+`0.58` for the uniform gap-scaled window).
+
+What *does* help is a non-uniform grid — and the reason this layer long
+declined to use one turned out to be false. The objection was that the GLS
+layer "whitens with a single AR(1) coefficient, which presumes a constant
+sample interval". That is a property of the **discrete parametrisation**, not
+of the noise. The stationary continuous-time process (Ornstein–Uhlenbeck,
+equivalently CAR(1)) has `Corr(t, t+d) = exp(−θ·d)` for any `d`, so on an
+arbitrary grid one whitens with the per-step `a_k = exp(−θ·dt_k)` and rescales
+by `sqrt(1 − a_k²)` to keep the result homoskedastic. Measured on a two-scale
+grid with exact OU noise, the median `|lag-1 autocorrelation|` of the whitened
+residuals is `0.374` with one constant `ρ` — taken at the median step, that
+scheme's best case — against `0.073` with the per-step coefficient. On a
+uniform grid both schemes give `0.073`, which is what shows the contrast comes
+from the grid rather than from the comparison.
+
+So the default window is now **repaired**, not merely disclosed. When the
+uniform grid cannot resolve the fastest mode (`max(−Re λ)`, taken from the
+spectrum, never guessed), half the points cover `[0, horizon/r_max]` and the
+rest carry the window out to `horizon/Δ`; `liouscope.fitting.car1` supplies the
+whitening, the exact `N_eff = n² / Σ_jk exp(−θ|t_j − t_k|)` and the
+exact-transition bootstrap resampler that go with it. On the reviewer's case
+(rates `1e-6` and `1`) the AICc winner M2 now recovers the fast rate as `1.10`
+against a true `1.0`, where the uniform window's M2 reported `2.17e-05` for it,
+and the D17 linear rate lands `0.021` from the gap instead of `0.58`.
+
+The disclosure remains for what the two-scale grid still cannot reach:
+
+- a **caller-supplied** `t_grid` that does not resolve its own system;
+- an **intermediate** timescale. The two segments resolve the fastest and the
+  slowest mode; a mode between them can fall entirely between the coarse late
+  samples. On three damped qubits at `1e-6`, `1e-3` and `1` the middle mode is
+  sampled `0.0019` times per e-folding and the warning fires, naming that mode.
+
+Below one sample per e-folding the layer emits an
+`UnderResolvedTransientWarning` and records
+`report.relaxation.samples_per_fast_efolding`. The reported rates then describe
+the dynamics the window *does* resolve, and a caller who needs the missing
+component must supply a `t_grid` covering it — reading the resulting rates as
+describing *that* window.
+
+This closes the time-grid unit dependence only. The `henrici_eta > 1.0` gate
+above is unaffected: `henrici_eta` is rate-dimensioned, so on an
+amplitude-damped qubit it equals the rescaling factor `c` exactly and still
+flips A5 → A10 between `c = 1` and `c = 3` whatever the grid. A third,
+independent dependence sits in the least-squares solver's own convergence
+controls (issue #111). Neither is asserted away by the grid work.

src/liouscope/_diagnostics.py

@@ -68,7 +68,20 @@ def diagnose(
     rho_steady_state
         Optional pre-computed steady state.
     t_grid
-        Time grid for the relaxation layer.
+        Time grid for the relaxation layer. When omitted the window is derived
+        from the system's own slowest relaxation time — ``[0, 10 / Delta]``
+        with ``Delta`` the D1 gap, sampled uniformly at 80 points (see
+        :func:`liouscope.diagnostics.relaxation.default_relaxation_grid`). The
+        default is therefore invariant under a pure change of rate units
+        ``L -> cL``; an absolute default would make every fitted rate, the AICc
+        model choice and the reported A-class depend on the caller's unit of
+        time. When the spectrum is spread too widely for one uniform grid to
+        resolve both ends (beyond roughly eightfold), the window becomes
+        two-scale and the residual model switches from AR(1) to CAR(1) with it.
+        Which window was used is recorded on
+        ``report.relaxation.t_grid_source`` (``"caller"`` / ``"gap_scaled"`` /
+        ``"gap_scaled_multiscale"`` / ``"legacy_fixed"``), ``.t_grid_span``,
+        ``.t_grid`` and ``.residual_model``.
     include_mpemba
         Compute D19/D20.
     bootstrap_B
@@ -162,6 +175,12 @@ def diagnose(
         rho_initial=rho_initial,
         rho_steady_state=rho_steady_state,
         t_grid=t_grid,
+        # D1 is already computed above; forwarding it keeps the DEFAULT
+        # relaxation window tied to the system's own relaxation time instead of
+        # an absolute one (see relaxation.default_relaxation_grid). Without
+        # this, a pure change of rate units L -> cL moved beta_D by >20% and
+        # changed the reported A-class, on identical physics.
+        gap=spectral.gap,
         bootstrap_B=bootstrap_B,
         seed=resolved_seed,
     )

src/liouscope/_types.py

@@ -187,6 +187,17 @@ class FitResult:
     n_eff: float
     residual_ar1_rho: float
     success: bool
+    # Which residual model the GLS layer actually used, and hence what
+    # ``residual_ar1_rho`` and ``n_eff`` mean:
+    #   NaN    -- discrete AR(1) on a uniform grid; ``residual_ar1_rho`` is the
+    #             fitted lag-1 correlation, ``n_eff`` the Geyer IPS estimate.
+    #   finite -- continuous-time CAR(1) on a NON-uniform grid; the whitening
+    #             used ``exp(-theta dt_k)`` per step, ``residual_ar1_rho`` is
+    #             that correlation at the mean step (reported for continuity
+    #             only), and ``n_eff`` is the exact CAR(1) value
+    #             ``n^2 / sum_jk exp(-theta |t_j - t_k|)``.
+    # Additive + defaulted, so older callers and serialised fits stay valid.
+    residual_theta_car1: float = float("nan")
     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
@@ -223,6 +234,42 @@ class RelaxationResult:
     # uses. Additive + defaulted so older callers / serialised reports stay valid.
     beta_D_linear: float = float("nan")
     linear_fit_model: str = "none"
+    # Provenance of the time grid the whole layer was fitted on. Every fitted
+    # quantity above is conditional on that window, so recording HOW it was
+    # chosen is part of the audit trail rather than a convenience:
+    #   "caller"       -- the caller supplied t_grid explicitly
+    #   "gap_scaled"   -- default window [0, HORIZON / Delta], rate-unit invariant
+    #   "legacy_fixed" -- no usable decay scale (Delta <= 0); absolute fallback
+    # Additive + defaulted, so older callers and serialised reports stay valid.
+    t_grid_source: str = "caller"
+    t_grid_span: float = float("nan")
+    # The grid ITSELF, not merely its span. The three curves above are y-values
+    # sampled on it, and a span alone does not identify the sampling: [0, 1, 10]
+    # and [0, 9, 10] share a span of 10 while describing materially different
+    # trajectories. Without this the exported report carries ordinates with no
+    # abscissa, so a downstream consumer cannot re-fit, re-plot or audit the
+    # rates it reports. Additive + defaulted like the fields above.
+    t_grid: np.ndarray | None = None
+    # How finely this grid samples the FASTEST decaying mode, as samples per
+    # e-folding. Below ``relaxation.MIN_SAMPLES_PER_FAST_EFOLD`` that mode was
+    # stepped over rather than measured, so the reported rates describe only
+    # the slow dynamics the window resolves and an
+    # ``UnderResolvedTransientWarning`` is emitted. ``inf`` when nothing decays.
+    samples_per_fast_efolding: float = float("nan")
+    # Residual model the reported M0..M3b hierarchy was ACTUALLY whitened
+    # with -- read off the fits, not off the grid (PR #127 review):
+    #   "ar1"                uniform grid, historical discrete AR(1);
+    #   "car1"               non-uniform grid, every fit whitened with the
+    #                        continuous-time exp(-theta dt_k) per step;
+    #   "car1_fallback_ar1"  non-uniform grid, but CAR(1) theta estimation
+    #                        failed on every fit (degenerate residuals), so
+    #                        ``fit_gls_ar1`` fell back to discrete AR(1);
+    #   "car1_mixed"         non-uniform grid, some fits CAR(1), some fallen
+    #                        back -- one label cannot cover the hierarchy;
+    #   "car1_unavailable"   non-uniform grid and no fit succeeded at all.
+    # The per-fit value is ``FitResult.residual_theta_car1``. Additive +
+    # defaulted.
+    residual_model: str = "ar1"
 
 
 @dataclass(frozen=True, slots=True, kw_only=True)