From da4c6cacc82b95b1a0e89b886c8f7138909d5d25 Mon Sep 17 00:00:00 2001
From: Marco Hostettler <marcohost33@gmail.com>
Date: Tue, 1 Sep 2026 20:18:22 +0200
Subject: [PATCH 1/8] fix(#124): make GLS optimizer residuals amplitude-scale
 invariant

---
 src/liouscope/fitting/gls.py | 15 ++++++++++++++-
 1 file changed, 14 insertions(+), 1 deletion(-)

diff --git a/src/liouscope/fitting/gls.py b/src/liouscope/fitting/gls.py
index 73b54f4..50ffffa 100644
--- a/src/liouscope/fitting/gls.py
+++ b/src/liouscope/fitting/gls.py
@@ -141,12 +141,25 @@ def fit_gls_ar1(
             degenerate=True,
         )
 
+    # Issue #124: the mathematical least-squares optimum is unchanged when an
+    # observable is multiplied by a positive constant, but SciPy's numerical
+    # termination is not. In particular, TRF declares success when its scaled
+    # gradient falls below ``gtol``; a curve such as ``1e-40*exp(-1.3*t)`` can
+    # therefore return the INITIAL RATE as a "converged" estimate after one
+    # evaluation. Divide only the residuals presented to the optimiser by the
+    # curve's own finite, non-zero scale. This multiplies the objective by one
+    # positive constant, so the minimiser is identical, while the numerical
+    # problem becomes amplitude-scale invariant. Raw residuals, AR(1) rho,
+    # sigma and likelihood below remain in the caller's original data units.
+    fit_scale = y_scale
+
     rho = 0.0
     success = True
     for _ in range(n_iters):
+
         def residual(params: np.ndarray, rho_local: float = rho) -> np.ndarray:
             y_hat = model(t, params)
-            r = y - y_hat
+            r = (y - y_hat) / fit_scale
             return _whiten(r, rho_local)
 
         ls_kwargs: dict[str, object] = {"max_nfev": max_nfev}

From 8a694a79029f71a2d702f0d46786cbbb8941e1b1 Mon Sep 17 00:00:00 2001
From: Marco Hostettler <marcohost33@gmail.com>
Date: Tue, 1 Sep 2026 20:18:44 +0200
Subject: [PATCH 2/8] test(#124): discriminate tiny-amplitude seed false
 convergence

---
 tests/test_gls_amplitude_scale.py | 54 +++++++++++++++++++++++++++++++
 1 file changed, 54 insertions(+)
 create mode 100644 tests/test_gls_amplitude_scale.py

diff --git a/tests/test_gls_amplitude_scale.py b/tests/test_gls_amplitude_scale.py
new file mode 100644
index 0000000..1269fe0
--- /dev/null
+++ b/tests/test_gls_amplitude_scale.py
@@ -0,0 +1,54 @@
+"""Regression tests for observable-amplitude invariance in the GLS optimiser.
+
+Issue #124: the mathematical fit is unchanged by ``y -> c*y`` for ``c > 0``,
+but feeding the raw residuals to SciPy lets the gradient termination criterion
+accept the initial seed when ``c`` is tiny.  These tests pin both directions:
+the tiny-amplitude counterexample must fit, and the ordinary-scale positive
+control must remain unchanged.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+from liouscope.fitting.gls import fit_gls_ar1
+
+_TRUE_RATE = 1.3
+_SEED_RATE = 0.2
+_T = np.linspace(0.0, 5.0, 64)
+
+
+def _fit_exponential(scale: float):
+    def model(t: np.ndarray, params: np.ndarray) -> np.ndarray:
+        return scale * np.exp(-params[0] * t)
+
+    y = scale * np.exp(-_TRUE_RATE * _T)
+    return fit_gls_ar1(
+        model,
+        _T,
+        y,
+        np.array([_SEED_RATE]),
+        bounds=(np.array([0.0]), np.array([5.0])),
+        n_iters=1,
+    )
+
+
+def test_tiny_amplitude_cannot_turn_the_seed_into_a_measurement() -> None:
+    """The exact #124 counterexample must move away from its initial seed."""
+    tiny = _fit_exponential(1.0e-40)
+
+    assert tiny.success
+    assert not tiny.degenerate
+    assert tiny.params[0] == pytest.approx(_TRUE_RATE, rel=1.0e-7, abs=1.0e-10)
+    assert abs(tiny.params[0] - _SEED_RATE) > 0.5
+
+
+def test_amplitude_rescaling_preserves_the_fitted_rate() -> None:
+    """Positive control: ordinary and tiny amplitudes represent the same fit."""
+    ordinary = _fit_exponential(1.0)
+    tiny = _fit_exponential(1.0e-40)
+
+    assert ordinary.success and tiny.success
+    assert ordinary.params[0] == pytest.approx(_TRUE_RATE, rel=1.0e-7, abs=1.0e-10)
+    assert tiny.params[0] == pytest.approx(ordinary.params[0], rel=1.0e-9, abs=1.0e-12)

From e675e619da3a03b2f192f5446363570a8a515cde Mon Sep 17 00:00:00 2001
From: Marco Hostettler <marcohost33@gmail.com>
Date: Tue, 1 Sep 2026 20:56:49 +0200
Subject: [PATCH 3/8] fix(#125): fail closed on degenerate BCa uncertainty

---
 src/liouscope/fitting/bootstrap.py | 29 +++++++++++++++++++++++++++++
 1 file changed, 29 insertions(+)

diff --git a/src/liouscope/fitting/bootstrap.py b/src/liouscope/fitting/bootstrap.py
index b936ee5..6bca6ce 100644
--- a/src/liouscope/fitting/bootstrap.py
+++ b/src/liouscope/fitting/bootstrap.py
@@ -13,6 +13,7 @@
 
 from __future__ import annotations
 
+import warnings
 from collections.abc import Callable
 
 import numpy as np
@@ -168,6 +169,14 @@ def bca_ci(
 
     Returns ``(p, 2)`` array of ``(lo, hi)`` per parameter at level
     ``1 - alpha``.
+
+    A parameter whose bootstrap distribution is degenerate, or whose adjusted
+    quantiles collapse to an exactly zero-width interval, is reported as
+    ``(NaN, NaN)`` with a warning. This follows the fail-closed semantics used
+    by SciPy's BCa implementation for degenerate bootstrap distributions:
+    arithmetic concentration is not promoted to a claim of perfect certainty
+    when the resampling model contains no measurable spread (issue #125).
+    Other non-degenerate parameters in the same fit remain reportable.
     """
     samples = np.asarray(samples, dtype=float)
     theta_hat = np.asarray(theta_hat, dtype=float)
@@ -177,6 +186,16 @@ def bca_ci(
     z_alpha_hi = ndtri(1.0 - alpha / 2.0)
     for j in range(p):
         boot_j = np.sort(samples[:, j])
+        if boot_j.size == 0 or np.all(boot_j == boot_j[0]):
+            warnings.warn(
+                f"bca_ci: bootstrap distribution for parameter {j} is "
+                "degenerate; uncertainty is unavailable rather than a "
+                "zero-width confidence claim (issue #125)",
+                RuntimeWarning,
+                stacklevel=2,
+            )
+            cis[j] = (float("nan"), float("nan"))
+            continue
         # Bias-correction z0 with half-correction for ties (Efron 1987, S3
         # audit 2026-06-04). Using a strict ``<`` only would send the
         # proportion to 0 (hence z0 -> -inf, clamped) whenever many bootstrap
@@ -214,6 +233,16 @@ def adjusted(z_alpha: float, z0_: float = z0, a_: float = a) -> float:
         # continuous-percentile estimator.
         lo = float(np.quantile(boot_j, q_lo, method="linear"))
         hi = float(np.quantile(boot_j, q_hi, method="linear"))
+        if lo == hi:
+            warnings.warn(
+                f"bca_ci: confidence interval for parameter {j} collapsed to "
+                "exactly zero width; uncertainty is unavailable rather than "
+                "perfectly known (issue #125)",
+                RuntimeWarning,
+                stacklevel=2,
+            )
+            cis[j] = (float("nan"), float("nan"))
+            continue
         cis[j] = (lo, hi)
     return cis
 

From b46daa8ca0642beb836783247ce47960cb7411d3 Mon Sep 17 00:00:00 2001
From: Marco Hostettler <marcohost33@gmail.com>
Date: Tue, 1 Sep 2026 20:57:05 +0200
Subject: [PATCH 4/8] test(#125): discriminate degenerate from ordinary BCa
 uncertainty

---
 tests/test_issue125_degenerate_bca.py | 57 +++++++++++++++++++++++++++
 1 file changed, 57 insertions(+)
 create mode 100644 tests/test_issue125_degenerate_bca.py

diff --git a/tests/test_issue125_degenerate_bca.py b/tests/test_issue125_degenerate_bca.py
new file mode 100644
index 0000000..7e2fa40
--- /dev/null
+++ b/tests/test_issue125_degenerate_bca.py
@@ -0,0 +1,57 @@
+"""Regression tests for issue #125: degenerate BCa uncertainty claims."""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+from liouscope.fitting.bootstrap import bca_ci
+
+
+def test_degenerate_bootstrap_distribution_is_unavailable_not_zero_width() -> None:
+    samples = np.ones((200, 1), dtype=float)
+    theta_hat = np.array([1.0])
+
+    with pytest.warns(RuntimeWarning, match="degenerate"):
+        ci = bca_ci(samples, theta_hat)
+
+    assert np.isnan(ci[0, 0])
+    assert np.isnan(ci[0, 1])
+
+
+def test_degenerate_parameter_does_not_withhold_an_independent_parameter() -> None:
+    x = np.linspace(-1.0, 1.0, 200)
+    samples = np.column_stack([np.ones_like(x), x])
+    theta_hat = np.array([1.0, 0.0])
+
+    with pytest.warns(RuntimeWarning, match="parameter 0"):
+        ci = bca_ci(samples, theta_hat)
+
+    assert np.all(np.isnan(ci[0]))
+    assert np.all(np.isfinite(ci[1]))
+    assert ci[1, 0] < ci[1, 1]
+
+
+def test_ordinary_nondegenerate_distribution_keeps_a_finite_interval() -> None:
+    rng = np.random.default_rng(125)
+    samples = rng.normal(loc=0.7, scale=0.03, size=(1000, 1))
+    theta_hat = np.array([0.7])
+
+    ci = bca_ci(samples, theta_hat)
+
+    assert np.all(np.isfinite(ci))
+    assert ci[0, 0] < ci[0, 1]
+    assert ci[0, 0] < theta_hat[0] < ci[0, 1]
+
+
+def test_zero_width_after_adjusted_quantiles_is_not_reported_as_certainty() -> None:
+    # Non-identical distribution, but >97.5% of its mass sits on one value.
+    # The BCa adjusted endpoints can therefore collapse to that value. The
+    # claim boundary, not merely the raw distribution, must catch the collapse.
+    samples = np.concatenate([np.ones(199), np.array([2.0])]).reshape(-1, 1)
+    theta_hat = np.array([1.0])
+
+    with pytest.warns(RuntimeWarning, match="zero width"):
+        ci = bca_ci(samples, theta_hat)
+
+    assert np.all(np.isnan(ci[0]))

From 7ad48290421a6d639c1fce8c1b29f35cacf25d97 Mon Sep 17 00:00:00 2001
From: marcohost33-maker <marcohost33@gmail.com>
Date: Fri, 11 Sep 2026 18:54:57 +0200
Subject: [PATCH 5/8] fix(#134): fall back to raw residuals where the #124
 rescaling is not representable

The #124 rescaling divides the optimiser residuals by max|y|. SciPy's
finite-difference probes step in absolute parameter units, so for a tiny
scale the rescaled residual, Jacobian or cost leaves float64; least_squares
then raised or failed and the fit was reported unsuccessful. That broke the
#147 contract (an unrepresentable MLE scale withholds only the CI, not the
fit): tests/test_issue135_likelihood_scale.py failed 2/21 after main was
merged in, and _fit_with_model("M0") regressed vs main at 1e-150 and 1e-310.

The rescaled solve now runs under np.errstate(call=...); any FP exception or
non-finite cost/fun/jac sends the iteration to the raw residuals (exactly
main's problem) with a RuntimeWarning. A finite rescaled solve that did not
converge is still a failure. Adds a Jacobian-level regression test.

Full-suite evidence run follows this commit.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0186oBqQsooNFqonGmyBd9Mp
---
 CHANGELOG.md                      |  17 +++++
 src/liouscope/fitting/gls.py      | 104 +++++++++++++++++++++++++++---
 tests/test_gls_amplitude_scale.py |  29 +++++++++
 3 files changed, 142 insertions(+), 8 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index fc8d85e..53f8fc0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -63,6 +63,23 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
     the change is not visible in `input_hash`.
 
 ### Fixed
+- **GLS optimiser residuals are divided by the curve's own scale, and fall
+  back to the raw residuals where that rescaling is not representable (issue
+  #124, PR #134).** Rescaling makes SciPy's termination amplitude-invariant
+  for a model whose amplitude is not a fitted parameter: `1e-40*exp(-1.3 t)`
+  no longer returns its seed rate as a converged fit. It does not cure a
+  FREE amplitude parameter, which still returns the seed rate at 1e-40 when
+  fitted from `p0 = [1e-40, 0.2]` without bounds.
+  SciPy's finite-difference probes, however, step in absolute parameter
+  units, so for a tiny scale the rescaled residual, Jacobian or cost can leave
+  float64 -- measured at `max|y| = 5e-324` (the #147 fixture) and, through
+  `_fit_with_model("M0")`, at 1e-150 and 1e-310, where the fit was reported
+  unsuccessful and dropped from AICc selection. Any floating-point exception or
+  non-finite result of the rescaled solve now sends that iteration to the raw
+  residuals (the pre-#124 problem) with a `RuntimeWarning` stating that the
+  #124 invariance does not hold for the curve; a finite rescaled solve that did
+  not converge is still reported unsuccessful. Raw residuals, rho, sigma and
+  the likelihood remain in the caller's data units.
 - **The trace-preservation defect overflowed on the way to a column sum that is
   exactly zero (issue #139, P2 review).** `trace_preservation_defect` assembled
   `vec(I)^H L` with an ordinary matrix product, which accumulates in ordinary
diff --git a/src/liouscope/fitting/gls.py b/src/liouscope/fitting/gls.py
index a5d9d0f..358e98c 100644
--- a/src/liouscope/fitting/gls.py
+++ b/src/liouscope/fitting/gls.py
@@ -18,7 +18,7 @@
 from dataclasses import dataclass
 
 import numpy as np
-from scipy.optimize import least_squares
+from scipy.optimize import OptimizeResult, least_squares
 
 from ..numerics.norms import scaled_log_sum_squares
 from .aicc import gaussian_log_likelihood
@@ -58,6 +58,15 @@ class GLSFitOutput:
     scale_unavailable: bool = False
 
 
+class _ScaledResidualOverflowError(ArithmeticError):
+    """A residual rescaled for the optimiser left float64 while the raw one is finite.
+
+    Private control flow of :func:`fit_gls_ar1` (PR #134). Deliberately NOT a
+    ``ValueError``/``RuntimeError``: those mean "the fit failed" there, while
+    this means only "the #124 rescaling is not representable for this curve".
+    """
+
+
 def _whiten(y: np.ndarray, rho: float) -> np.ndarray:
     if y.size < 2:
         y_copy: np.ndarray = y.copy()
@@ -191,22 +200,101 @@ def fit_gls_ar1(
     # positive constant, so the minimiser is identical, while the numerical
     # problem becomes amplitude-scale invariant. Raw residuals, AR(1) rho,
     # sigma and likelihood below remain in the caller's original data units.
+    #
+    # The rescaled problem is not always REPRESENTABLE (PR #134 x PR #147).
+    # SciPy's finite-difference probes step by ``~1.5e-8 * max(1, |x|)`` in
+    # ABSOLUTE parameter units, so the model moves by an amount unrelated to
+    # the curve's scale. Measured on the #147 fixture (``max|y| = 5e-324``):
+    # the scaled residual at ``p0`` is finite (``max = 1.0``), the first
+    # Jacobian probe gives ``1.5e-8 / 5e-324 = inf``, ``least_squares`` raises
+    # and the fit was reported unsuccessful -- a curve withheld from model
+    # selection because of its absolute amplitude, the boundary #135 removed.
+    # The overflow is a property of the rescaling, not of the fit: when a
+    # scaled residual leaves float64 while the raw one is finite, the
+    # iteration is repeated on the raw residuals, exactly the pre-#124
+    # problem, and the loss of the #124 invariance is announced rather than
+    # silent. A raw residual that is itself non-finite still fails closed.
     fit_scale = y_scale
 
+    def _residual_fn(
+        rho_local: float, scale: float
+    ) -> Callable[[np.ndarray], np.ndarray]:
+        def residual(params: np.ndarray) -> np.ndarray:
+            y_hat = model(t, params)
+            r_raw = y - y_hat
+            # errstate: the overflow is DETECTED below, so it must not also
+            # surface as a numpy RuntimeWarning (an error under this repo's
+            # ``filterwarnings = error``) before the detection can act.
+            with np.errstate(over="ignore", divide="ignore", invalid="ignore"):
+                scaled = _whiten(r_raw / scale, rho_local)
+            if not np.all(np.isfinite(scaled)) and np.all(np.isfinite(r_raw)):
+                raise _ScaledResidualOverflowError
+            return scaled
+
+        return residual
+
+    # The residual is not the only place the rescaled problem can leave
+    # float64: with a FREE amplitude parameter the Jacobian and the cost carry
+    # the same ``1/scale`` factor. Measured through ``_fit_with_model("M0")``
+    # on ``s*exp(-1.3 t)``: at s = 1e-150 and 1e-310 the rescaled optimiser
+    # overflowed inside SciPy (``dot``/``square``) and the fit was reported
+    # unsuccessful (aicc = inf), where the raw problem fits rate 1.3. So every
+    # floating-point exception raised while the RESCALED problem is solved,
+    # and any non-finite cost/residual/Jacobian it returns, counts as "not
+    # representable" and sends the iteration to the raw residuals. A finite
+    # rescaled result that merely did not converge is NOT retried: that is a
+    # genuine failure, and retrying it raw would re-admit the #124 seed.
+    def _solve_rescaled(
+        rho_local: float, x0: np.ndarray, **kw: object
+    ) -> OptimizeResult:
+        events: list[str] = []
+
+        def _record(kind: str, _flag: int) -> None:
+            events.append(kind)
+
+        with np.errstate(
+            over="call", divide="call", invalid="call", under="ignore", call=_record
+        ):
+            try:
+                res = least_squares(_residual_fn(rho_local, fit_scale), x0, **kw)
+            except (ValueError, RuntimeError):
+                if events:
+                    raise _ScaledResidualOverflowError from None
+                raise
+        if events or not (
+            np.isfinite(res.cost)
+            and np.all(np.isfinite(res.fun))
+            and np.all(np.isfinite(res.jac))
+        ):
+            raise _ScaledResidualOverflowError
+        return res
+
     rho = 0.0
     success = True
     for _ in range(n_iters):
-
-        def residual(params: np.ndarray, rho_local: float = rho) -> np.ndarray:
-            y_hat = model(t, params)
-            r = (y - y_hat) / fit_scale
-            return _whiten(r, rho_local)
-
         ls_kwargs: dict[str, object] = {"max_nfev": max_nfev}
         if bounds is not None:
             ls_kwargs["bounds"] = bounds
         try:
-            result = least_squares(residual, p, **ls_kwargs)
+            try:
+                if fit_scale == 1.0:
+                    result = least_squares(
+                        _residual_fn(rho, fit_scale), p, **ls_kwargs
+                    )
+                else:
+                    result = _solve_rescaled(rho, p, **ls_kwargs)
+            except _ScaledResidualOverflowError:
+                warnings.warn(
+                    "fit_gls_ar1: residuals rescaled by the curve's own scale "
+                    f"({fit_scale:.3e}) are not representable as float64 at the "
+                    "optimiser's probe points; fitting the unscaled residuals "
+                    "instead, so the amplitude-scale invariance of issue #124 "
+                    "does not hold for this curve (PR #134).",
+                    RuntimeWarning,
+                    stacklevel=2,
+                )
+                fit_scale = 1.0
+                result = least_squares(_residual_fn(rho, fit_scale), p, **ls_kwargs)
             p = result.x
             success = result.success
         except (ValueError, RuntimeError):
diff --git a/tests/test_gls_amplitude_scale.py b/tests/test_gls_amplitude_scale.py
index 1269fe0..923b7fc 100644
--- a/tests/test_gls_amplitude_scale.py
+++ b/tests/test_gls_amplitude_scale.py
@@ -52,3 +52,32 @@ def test_amplitude_rescaling_preserves_the_fitted_rate() -> None:
     assert ordinary.success and tiny.success
     assert ordinary.params[0] == pytest.approx(_TRUE_RATE, rel=1.0e-7, abs=1.0e-10)
     assert tiny.params[0] == pytest.approx(ordinary.params[0], rel=1.0e-9, abs=1.0e-12)
+
+
+@pytest.mark.parametrize("scale", [1.0e-150, 1.0e-310])
+def test_unrepresentable_rescaling_falls_back_instead_of_dropping_the_fit(
+    scale: float,
+) -> None:
+    """PR #134: the #124 rescaling must not withhold a fit the raw problem makes.
+
+    With a FREE amplitude parameter the rescaled Jacobian and cost carry the
+    same ``1/scale`` factor as the residual. Measured on the merge of #124 into
+    main: ``_fit_with_model("M0")`` overflowed inside SciPy at these scales and
+    returned ``success=False`` / ``aicc=inf``, while main fitted rate 1.3. The
+    Jacobian-level overflow is the case the residual-level check cannot see, so
+    it is pinned here separately from the #147 subnormal fixture.
+    """
+    import warnings
+
+    from liouscope.diagnostics import relaxation as relaxation_mod
+
+    y = scale * np.exp(-_TRUE_RATE * _T)
+    with warnings.catch_warnings(record=True) as caught:
+        warnings.simplefilter("always")
+        fit_result, _ = relaxation_mod._fit_with_model("M0", _T, y)
+    assert any(
+        "not representable as float64 at the" in str(w.message) for w in caught
+    ), [str(w.message) for w in caught]
+    assert fit_result.success
+    assert np.isfinite(fit_result.aicc)
+    assert fit_result.params[1] == pytest.approx(_TRUE_RATE, rel=1.0e-7)

From 9a623e8e1c1a23505f25466ff108d96a23fd748c Mon Sep 17 00:00:00 2001
From: marcohost33-maker <marcohost33@gmail.com>
Date: Fri, 11 Sep 2026 19:38:57 +0200
Subject: [PATCH 6/8] fix(#134): keep the raw path verbatim and count only
 rescaling FP events

Equalita NO-GO for 7ad4829:
1. The residual detector also ran on the RAW residuals; for max|y| >= ~1e155
   their whitening overflows although they are finite, and the private
   _ScaledResidualOverflowError escaped fit_gls_ar1 and _fit_with_model
   (main: success=False). The raw path is now main's residual verbatim.
2. Three detector parts were unpinned. The explicit residual check is gone
   (the division now counts as an FP event itself); one discriminating test
   each for the non-finite arm, exception-without-event, and the
   no-retry-on-nonconverged rule. Own mutation run: 7/7 mutants red.
3. The model is evaluated under the caller's FP policy, so a benign model
   event (sinc 0/0 at t=0) no longer forces the fallback: 1e-40 now fits
   rate 1.3003772913207243 (2.0: 1.3003772913320941) instead of the seed.
4. CHANGELOG corrected.

Full-suite evidence run follows this commit (~13 min).

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0186oBqQsooNFqonGmyBd9Mp
---
 CHANGELOG.md                      |  17 ++--
 src/liouscope/fitting/gls.py      |  54 +++++++++----
 tests/test_gls_amplitude_scale.py | 124 ++++++++++++++++++++++++++++++
 3 files changed, 173 insertions(+), 22 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 53f8fc0..6c529ce 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -74,11 +74,18 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
   units, so for a tiny scale the rescaled residual, Jacobian or cost can leave
   float64 -- measured at `max|y| = 5e-324` (the #147 fixture) and, through
   `_fit_with_model("M0")`, at 1e-150 and 1e-310, where the fit was reported
-  unsuccessful and dropped from AICc selection. Any floating-point exception or
-  non-finite result of the rescaled solve now sends that iteration to the raw
-  residuals (the pre-#124 problem) with a `RuntimeWarning` stating that the
-  #124 invariance does not hold for the curve; a finite rescaled solve that did
-  not converge is still reported unsuccessful. Raw residuals, rho, sigma and
+  unsuccessful and dropped from AICc selection. A floating-point exception
+  raised by the rescaled arithmetic itself (the division, the whitening, or
+  SciPy's own computation on the rescaled values -- NOT the evaluation of the
+  model, which keeps the caller's floating-point policy) or a non-finite
+  cost/residual/Jacobian of the rescaled solve now sends that iteration to the
+  raw residuals, evaluated exactly as before #124, with a `RuntimeWarning`
+  stating that the #124 invariance does not hold for the curve. A finite
+  rescaled solve that did not converge, and an exception with no
+  floating-point event, are still reported unsuccessful. For `max|y| >= ~1e150`
+  the fallback fires as well and the fit fails closed exactly as on main
+  (`success=False`), without leaking the private exception that an earlier
+  revision of this fix raised there (PR #134 review). Raw residuals, rho, sigma and
   the likelihood remain in the caller's data units.
 - **The trace-preservation defect overflowed on the way to a column sum that is
   exactly zero (issue #139, P2 review).** `trace_preservation_defect` assembled
diff --git a/src/liouscope/fitting/gls.py b/src/liouscope/fitting/gls.py
index 358e98c..de8742a 100644
--- a/src/liouscope/fitting/gls.py
+++ b/src/liouscope/fitting/gls.py
@@ -16,6 +16,7 @@
 import warnings
 from collections.abc import Callable
 from dataclasses import dataclass
+from typing import Any
 
 import numpy as np
 from scipy.optimize import OptimizeResult, least_squares
@@ -216,20 +217,35 @@ def fit_gls_ar1(
     # silent. A raw residual that is itself non-finite still fails closed.
     fit_scale = y_scale
 
-    def _residual_fn(
-        rho_local: float, scale: float
+    def _raw_residual_fn(rho_local: float) -> Callable[[np.ndarray], np.ndarray]:
+        # The unscaled problem, verbatim as before #124: no detection, no
+        # errstate. A non-finite raw residual reaches ``least_squares``, which
+        # raises ``ValueError`` and the fit fails closed (PR #134 review: the
+        # detector used to run here too and leaked a private exception from
+        # the public API for ``max|y| >= ~1e155``).
+        def residual(params: np.ndarray) -> np.ndarray:
+            return _whiten(y - model(t, params), rho_local)
+
+        return residual
+
+    def _rescaled_residual_fn(
+        rho_local: float,
+        scale: float,
+        caller_err: dict[str, Any],
+        caller_call: Any,
     ) -> Callable[[np.ndarray], np.ndarray]:
+        # Runs inside the counting errstate of ``_solve_rescaled``. The MODEL
+        # and the raw difference are evaluated under the CALLER's floating-point
+        # policy instead: an event there (e.g. a sinc factor's 0/0 at t = 0,
+        # finite after ``np.where``) belongs to the raw problem, not to the
+        # rescaling, and must not trigger the fallback (PR #134 review: it did,
+        # and returned the seed rate as a converged fit at 1e-40).
         def residual(params: np.ndarray) -> np.ndarray:
-            y_hat = model(t, params)
-            r_raw = y - y_hat
-            # errstate: the overflow is DETECTED below, so it must not also
-            # surface as a numpy RuntimeWarning (an error under this repo's
-            # ``filterwarnings = error``) before the detection can act.
-            with np.errstate(over="ignore", divide="ignore", invalid="ignore"):
-                scaled = _whiten(r_raw / scale, rho_local)
-            if not np.all(np.isfinite(scaled)) and np.all(np.isfinite(r_raw)):
-                raise _ScaledResidualOverflowError
-            return scaled
+            with np.errstate(call=caller_call, **caller_err):
+                r_raw = y - model(t, params)
+            # Counted: an overflow of the division or the whitening IS the
+            # rescaled problem leaving float64.
+            return _whiten(r_raw / scale, rho_local)
 
         return residual
 
@@ -252,11 +268,17 @@ def _solve_rescaled(
         def _record(kind: str, _flag: int) -> None:
             events.append(kind)
 
+        caller_err = dict(np.geterr())
+        caller_call = np.geterrcall()
         with np.errstate(
             over="call", divide="call", invalid="call", under="ignore", call=_record
         ):
             try:
-                res = least_squares(_residual_fn(rho_local, fit_scale), x0, **kw)
+                res = least_squares(
+                    _rescaled_residual_fn(rho_local, fit_scale, caller_err, caller_call),
+                    x0,
+                    **kw,
+                )
             except (ValueError, RuntimeError):
                 if events:
                     raise _ScaledResidualOverflowError from None
@@ -278,9 +300,7 @@ def _record(kind: str, _flag: int) -> None:
         try:
             try:
                 if fit_scale == 1.0:
-                    result = least_squares(
-                        _residual_fn(rho, fit_scale), p, **ls_kwargs
-                    )
+                    result = least_squares(_raw_residual_fn(rho), p, **ls_kwargs)
                 else:
                     result = _solve_rescaled(rho, p, **ls_kwargs)
             except _ScaledResidualOverflowError:
@@ -294,7 +314,7 @@ def _record(kind: str, _flag: int) -> None:
                     stacklevel=2,
                 )
                 fit_scale = 1.0
-                result = least_squares(_residual_fn(rho, fit_scale), p, **ls_kwargs)
+                result = least_squares(_raw_residual_fn(rho), p, **ls_kwargs)
             p = result.x
             success = result.success
         except (ValueError, RuntimeError):
diff --git a/tests/test_gls_amplitude_scale.py b/tests/test_gls_amplitude_scale.py
index 923b7fc..99df7ad 100644
--- a/tests/test_gls_amplitude_scale.py
+++ b/tests/test_gls_amplitude_scale.py
@@ -81,3 +81,127 @@ def test_unrepresentable_rescaling_falls_back_instead_of_dropping_the_fit(
     assert fit_result.success
     assert np.isfinite(fit_result.aicc)
     assert fit_result.params[1] == pytest.approx(_TRUE_RATE, rel=1.0e-7)
+
+
+# --------------------------------------------------------------------------
+# PR #134 review (Equalita NO-GO for 7ad4829): one discriminating test per
+# part of the representability detector.
+# --------------------------------------------------------------------------
+
+_FALLBACK = "not representable as float64 at the"
+_NOISE = np.random.default_rng(20260911).standard_normal(_T.size)
+
+
+def _free_amplitude(t: np.ndarray, params: np.ndarray) -> np.ndarray:
+    return params[0] * np.exp(-params[1] * t)
+
+
+def _recorded(fn, *args, **kwargs):  # type: ignore[no-untyped-def]
+    import warnings
+
+    with warnings.catch_warnings(record=True) as caught:
+        warnings.simplefilter("always")
+        out = fn(*args, **kwargs)
+    return out, [str(w.message) for w in caught]
+
+
+@pytest.mark.parametrize("scale", [1.0e200, 1.0e300])
+def test_large_amplitude_fails_closed_without_leaking_the_private_exception(
+    scale: float,
+) -> None:
+    """The raw path is main's problem verbatim: ``success=False``, no raise.
+
+    Measured on 7ad4829: the residual detector also ran on the RAW residuals,
+    whose whitening overflows here although they are finite, and its private
+    ``ArithmeticError`` escaped both public entry points (20/20 cases).
+    """
+    from liouscope.diagnostics import relaxation as relaxation_mod
+
+    y = scale * np.exp(-_TRUE_RATE * _T)
+    fit, _ = _recorded(fit_gls_ar1, _free_amplitude, _T, y, np.array([scale, 0.2]))
+    assert not fit.success
+    (fit_result, _), _ = _recorded(
+        relaxation_mod._fit_with_model, "M0", _T, y * (1.0 + 1.0e-3 * _NOISE)
+    )
+    assert not fit_result.success
+    assert np.isinf(fit_result.aicc)
+
+
+def test_benign_model_float_event_does_not_trigger_the_fallback() -> None:
+    """A 0/0 inside the MODEL belongs to the raw problem, not to the rescaling.
+
+    ``sin(t)/t`` at ``t = 0`` raises 'invalid' and is replaced by ``np.where``;
+    the model output is finite. When that event was counted, every fit of this
+    model fell back and returned the seed rate 0.2 as converged at 1e-40.
+    """
+
+    def run(scale: float) -> tuple[float, list[str]]:
+        def model(t: np.ndarray, params: np.ndarray) -> np.ndarray:
+            return scale * np.exp(-params[0] * t) * np.where(t > 0, np.sin(t) / t, 1.0)
+
+        with np.errstate(invalid="ignore"):  # building the data is not under test
+            y = model(_T, np.array([_TRUE_RATE])) * (1.0 + 1.0e-3 * _NOISE)
+        fit, messages = _recorded(
+            fit_gls_ar1, model, _T, y, np.array([_SEED_RATE]),
+            bounds=(np.array([0.0]), np.array([5.0])), n_iters=1,
+        )
+        assert fit.success
+        return float(fit.params[0]), messages
+
+    ordinary, _ = run(2.0)
+    tiny, messages = run(1.0e-40)
+    assert not any(_FALLBACK in m for m in messages), messages
+    assert tiny == pytest.approx(ordinary, rel=1.0e-9)
+    assert abs(tiny - _SEED_RATE) > 0.5
+
+
+def test_nonfinite_rescaled_result_without_float_events_falls_back(
+    monkeypatch: pytest.MonkeyPatch,
+) -> None:
+    """The non-finite arm: a NaN Jacobian set by assignment raises no FP event."""
+    import liouscope.fitting.gls as gls_mod
+
+    real = gls_mod.least_squares
+    calls = {"n": 0}
+
+    def first_call_poisoned(*args, **kwargs):  # type: ignore[no-untyped-def]
+        res = real(*args, **kwargs)
+        calls["n"] += 1
+        if calls["n"] == 1:
+            res.jac = np.full_like(res.jac, np.nan)
+        return res
+
+    monkeypatch.setattr(gls_mod, "least_squares", first_call_poisoned)
+    y = 1.0e-40 * np.exp(-_TRUE_RATE * _T) * (1.0 + 1.0e-3 * _NOISE)
+    fit, messages = _recorded(
+        fit_gls_ar1, lambda t, p: 1.0e-40 * np.exp(-p[0] * t), _T, y,
+        np.array([_SEED_RATE]), bounds=(np.array([0.0]), np.array([5.0])), n_iters=1,
+    )
+    assert any(_FALLBACK in m for m in messages), messages
+    assert calls["n"] == 2
+
+
+def test_exception_without_float_events_is_not_rerouted() -> None:
+    """A model error in the rescaled solve is a failed fit, not a scale problem."""
+
+    def failing(t: np.ndarray, params: np.ndarray) -> np.ndarray:
+        if params[0] != _SEED_RATE:
+            raise RuntimeError("model refuses this rate")
+        return 1.0e-40 * np.exp(-params[0] * t)
+
+    y = 1.0e-40 * np.exp(-_TRUE_RATE * _T) * (1.0 + 1.0e-3 * _NOISE)
+    fit, messages = _recorded(fit_gls_ar1, failing, _T, y, np.array([_SEED_RATE]), n_iters=1)
+    assert not fit.success
+    assert not any(_FALLBACK in m for m in messages), messages
+
+
+def test_finite_nonconverged_rescaled_solve_is_not_retried_raw() -> None:
+    """Retrying a finite non-converged solve raw would re-admit the #124 seed."""
+    y = 1.0e-40 * np.exp(-_TRUE_RATE * _T) * (1.0 + 1.0e-3 * _NOISE)
+    fit, messages = _recorded(
+        fit_gls_ar1, lambda t, p: 1.0e-40 * np.exp(-p[0] * t), _T, y,
+        np.array([_SEED_RATE]), bounds=(np.array([0.0]), np.array([5.0])),
+        n_iters=1, max_nfev=1,
+    )
+    assert not fit.success
+    assert not any(_FALLBACK in m for m in messages), messages

From c218f0944a74b3b58386246c06a5ad94b7af9cf5 Mon Sep 17 00:00:00 2001
From: marcohost33-maker <marcohost33@gmail.com>
Date: Fri, 11 Sep 2026 20:16:10 +0200
Subject: [PATCH 7/8] fix(#134 round 3): scale-safe AR(1) rho, CITATION entry,
 single import form

1. Codex P1 (gls.py:248): rho between Cochrane-Orcutt iterations came from
   ar1_correlation on UNSCALED residuals; its dot products underflow to 0
   below ~1e-162 (corrected rho collapses to the floor 1/(n-3)) and overflow
   to NaN above ~1e154. Measured with n_iters=3 on one AR(1) curve: rho
   0.012987 / 0.44637 / NaN (success=False) at 1e-170 / 1e0 / 1e160.
   ar1_correlation now normalises by an exact power of two (frexp/ldexp);
   rho is bit-identical where the old arithmetic stayed in range (200/200
   random series 1e-100..1e100) and 0.44636512337/0.44636512373/
   0.44636512351 at the three scales after the fix.
2. Codex P1 (CHANGELOG.md:68): CITATION.cff "Also pending for the next cut"
   block for the GLS rescaling, fallback and rho normalisation (DoD #5).
3. CodeQL: test_gls_amplitude_scale.py no longer imports liouscope.fitting.gls
   both as a module and via from-import.

Full-suite evidence run follows this commit (~9-13 min).

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0186oBqQsooNFqonGmyBd9Mp
---
 CHANGELOG.md                      |  9 +++++
 CITATION.cff                      | 16 +++++++++
 src/liouscope/fitting/neff.py     | 17 ++++++++-
 tests/test_gls_amplitude_scale.py | 58 +++++++++++++++++++++++++++++--
 4 files changed, 96 insertions(+), 4 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6c529ce..7533b5e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -87,6 +87,15 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
   (`success=False`), without leaking the private exception that an earlier
   revision of this fix raised there (PR #134 review). Raw residuals, rho, sigma and
   the likelihood remain in the caller's data units.
+- **The AR(1) lag-1 autocorrelation is formed from residuals normalised by an
+  exact power of two (PR #134, round-3 review).** `ar1_correlation` took its dot
+  products on unscaled residuals, which underflow to 0 below ~1e-162 -- the
+  corrected rho then collapses to its floor `1/(n-3)` -- and overflow to NaN
+  above ~1e154. Measured through `fit_gls_ar1` with `n_iters=3` on one AR(1)
+  correlated curve: rho 0.012987 at 1e-170, 0.44637 at 1e0 and NaN with
+  `success=False` at 1e160. Scaling by `2**k` is exact, so rho is bit-identical
+  wherever the old computation stayed in range (200/200 random series between
+  1e-100 and 1e100).
 - **The trace-preservation defect overflowed on the way to a column sum that is
   exactly zero (issue #139, P2 review).** `trace_preservation_defect` assembled
   `vec(I)^H L` with an ordinary matrix product, which accumulates in ordinary
diff --git a/CITATION.cff b/CITATION.cff
index eecf07a..8d0e720 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -140,6 +140,22 @@ abstract: >-
 # MANIFEST_SCHEMA is deliberately unchanged -- `trace_defect` is a report
 # field, not a run-manifest field, so the manifest contract is untouched and
 # `input_hash` does not move.
+# Also pending for the next cut: the GLS AR(1) fit presents the optimiser with
+# residuals divided by the curve's own scale max|y| (issue #124), so SciPy's
+# termination no longer turns the seed rate into a "converged" estimate for a
+# tiny-amplitude curve whose amplitude is not a fitted parameter; where that
+# rescaled problem leaves float64, the iteration falls back, with a
+# RuntimeWarning, to the unscaled residuals exactly as before (PR #134). The
+# lag-1 autocorrelation behind the AR(1) whitening is formed from residuals
+# normalised by an exact power of two, so it no longer underflows to its floor
+# below ~1e-162 or overflows to NaN above ~1e154; where the old computation
+# stayed in range it is bit-identical. This CHANGES NUMERICAL RESULTS and is a
+# METHODOLOGY correction: fitted rates, rho, AICc and the selected model can
+# differ for curves far from unit scale. It must NOT be described as making the
+# fit amplitude-invariant in general: a FREE amplitude parameter still returns
+# its seed rate at 1e-40, and the fallback keeps the unscaled problem's own
+# limits. MANIFEST_SCHEMA is deliberately unchanged -- none of these quantities
+# is a run-manifest field, so `input_hash` does not move.
 keywords:
   - open quantum systems
   - Lindblad
diff --git a/src/liouscope/fitting/neff.py b/src/liouscope/fitting/neff.py
index b888b59..5b932f2 100644
--- a/src/liouscope/fitting/neff.py
+++ b/src/liouscope/fitting/neff.py
@@ -15,6 +15,7 @@
 
 from __future__ import annotations
 
+import math
 import warnings
 
 import numpy as np
@@ -100,9 +101,23 @@ def ar1_correlation(residuals: np.ndarray) -> float:
     coefficient is estimator-convention dependent, so treat it as indicative
     rather than a pinned analytic identity.
     """
-    x = np.asarray(residuals, dtype=float) - float(np.mean(residuals))
+    x = np.asarray(residuals, dtype=float)
     if x.size < 2:
         return 0.0
+    # rho is a RATIO of two quadratic forms, so it is invariant under a common
+    # rescaling of the residuals -- but the dot products that form it are not:
+    # residuals below ~1e-162 square to 0 (rho_hat = 0, and the corrected
+    # estimate collapses to its floor 1/(n-3)), and above ~1e154 they overflow
+    # to NaN (PR #134 review, measured through fit_gls_ar1 at 1e-170 / 1e160).
+    # Normalise by the exact power of two nearest the peak magnitude. Scaling
+    # by 2**k is exact in binary floating point whenever nothing under- or
+    # overflows, so every mean, product and sum below is the old one times an
+    # exact power of two and rho_hat is BIT-IDENTICAL wherever the unnormalised
+    # computation stayed in range -- the audit-pinned formula does not move.
+    peak = float(np.max(np.abs(x)))
+    if math.isfinite(peak) and peak > 0.0:
+        x = np.ldexp(x, -math.frexp(peak)[1])
+    x = x - float(np.mean(x))
     num = float(np.dot(x[:-1], x[1:]))
     den = float(np.dot(x, x))
     if den == 0.0:
diff --git a/tests/test_gls_amplitude_scale.py b/tests/test_gls_amplitude_scale.py
index 99df7ad..82c536c 100644
--- a/tests/test_gls_amplitude_scale.py
+++ b/tests/test_gls_amplitude_scale.py
@@ -159,9 +159,9 @@ def test_nonfinite_rescaled_result_without_float_events_falls_back(
     monkeypatch: pytest.MonkeyPatch,
 ) -> None:
     """The non-finite arm: a NaN Jacobian set by assignment raises no FP event."""
-    import liouscope.fitting.gls as gls_mod
+    from scipy.optimize import least_squares as real
+
 
-    real = gls_mod.least_squares
     calls = {"n": 0}
 
     def first_call_poisoned(*args, **kwargs):  # type: ignore[no-untyped-def]
@@ -171,7 +171,7 @@ def first_call_poisoned(*args, **kwargs):  # type: ignore[no-untyped-def]
             res.jac = np.full_like(res.jac, np.nan)
         return res
 
-    monkeypatch.setattr(gls_mod, "least_squares", first_call_poisoned)
+    monkeypatch.setattr("liouscope.fitting.gls.least_squares", first_call_poisoned)
     y = 1.0e-40 * np.exp(-_TRUE_RATE * _T) * (1.0 + 1.0e-3 * _NOISE)
     fit, messages = _recorded(
         fit_gls_ar1, lambda t, p: 1.0e-40 * np.exp(-p[0] * t), _T, y,
@@ -205,3 +205,55 @@ def test_finite_nonconverged_rescaled_solve_is_not_retried_raw() -> None:
     )
     assert not fit.success
     assert not any(_FALLBACK in m for m in messages), messages
+
+
+# --------------------------------------------------------------------------
+# PR #134 round 3 (Codex P1): the AR(1) rho between Cochrane-Orcutt iterations
+# was estimated from UNSCALED residuals, whose dot products underflow to 0
+# below ~1e-162 and overflow to NaN above ~1e154.
+# --------------------------------------------------------------------------
+
+
+def _ar1_curve(scale: float) -> tuple[np.ndarray, np.ndarray]:
+    t = np.linspace(0.0, 5.0, 80)
+    rng = np.random.default_rng(20260911)
+    e = np.zeros(t.size)
+    nu = rng.standard_normal(t.size)
+    for i in range(t.size):
+        e[i] = (0.6 * e[i - 1] if i else 0.0) + nu[i]
+    return t, scale * (np.exp(-_TRUE_RATE * t) + 1.0e-3 * e)
+
+
+@pytest.mark.parametrize("scale", [1.0e-170, 1.0e160])
+def test_ar1_rho_between_iterations_is_amplitude_invariant(scale: float) -> None:
+    """Amplitude-equivalent correlated curves give the same rho, rate, success.
+
+    Measured on 9a623e8: rho 0.012987 (= the corrected floor 1/(n-3), i.e. an
+    underflowed rho_hat of 0) at 1e-170, 0.44637 at 1e0, and NaN with
+    ``success=False`` at 1e160.
+    """
+
+    def fit(s: float):  # type: ignore[no-untyped-def]
+        t, y = _ar1_curve(s)
+        out, _ = _recorded(
+            fit_gls_ar1, lambda tt, p: s * p[0] * np.exp(-p[1] * tt), t, y,
+            np.array([1.0, 0.8]), n_iters=3,
+        )
+        return out
+
+    ref = fit(1.0)
+    other = fit(scale)
+    assert ref.success and np.isfinite(ref.rho_ar1) and ref.rho_ar1 > 0.3
+    assert other.success == ref.success
+    assert other.rho_ar1 == pytest.approx(ref.rho_ar1, rel=1.0e-8)
+    assert other.params[1] == pytest.approx(ref.params[1], rel=1.0e-8)
+
+
+def test_ar1_correlation_is_bit_identical_under_power_of_two_rescaling() -> None:
+    """The normalisation must not move rho where the old arithmetic was in range."""
+    from liouscope.fitting.neff import ar1_correlation
+
+    _, y = _ar1_curve(1.0)
+    base = ar1_correlation(y)
+    for k in (-900, -540, 500, 1000):
+        assert ar1_correlation(np.ldexp(y, k)) == base, k

From 28c7db2c401331ed6bcb78040a0c73a128689ada Mon Sep 17 00:00:00 2001
From: marcohost33-maker <marcohost33@gmail.com>
Date: Fri, 11 Sep 2026 20:33:35 +0200
Subject: [PATCH 8/8] fix(#134 round 3b): pin the caller's FP policy,
 filterable fallback warning

Equalita follow-up on 9a623e8:
(a) Mutation E3 (model under the detector's errstate) stayed green. New
    test: the model raises 0/0 only at the seed, i.e. inside the solve, so
    the post-fit evaluation cannot mask it; with np.errstate(invalid="raise")
    the fit must raise FloatingPointError, as it does on main. Red under E3
    and under "model counted", green with the fix and against main's gls.py.
(b) CHANGELOG: the fallback at max|y| >= ~1e150 applies to a FREE amplitude
    only (fixed amplitude: no fallback, rate 1.3 at 1e150/1e200/1e300), and
    the warning is new relative to main. Free amplitude at 1e-40 still
    returns the seed 0.2 after the rho change (claim re-measured).
(c) The fallback warning is AmplitudeRescalingFallbackWarning, a
    RuntimeWarning subclass, emitted per fit, so -W error users can filter
    it by class. Test red when the category is reverted to RuntimeWarning.

Full-suite evidence run follows this commit (~7-13 min).

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0186oBqQsooNFqonGmyBd9Mp
---
 CHANGELOG.md                      | 13 +++++---
 src/liouscope/fitting/gls.py      | 15 +++++++++-
 tests/test_gls_amplitude_scale.py | 50 +++++++++++++++++++++++++++++++
 3 files changed, 73 insertions(+), 5 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7533b5e..6d40ff1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -82,10 +82,15 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
   raw residuals, evaluated exactly as before #124, with a `RuntimeWarning`
   stating that the #124 invariance does not hold for the curve. A finite
   rescaled solve that did not converge, and an exception with no
-  floating-point event, are still reported unsuccessful. For `max|y| >= ~1e150`
-  the fallback fires as well and the fit fails closed exactly as on main
-  (`success=False`), without leaking the private exception that an earlier
-  revision of this fix raised there (PR #134 review). Raw residuals, rho, sigma and
+  floating-point event, are still reported unsuccessful. For a FREE amplitude
+  parameter at `max|y| >= ~1e150` the fallback fires as well (measured at
+  1e150, 1e200 and 1e300; a fixed-amplitude model does not fall back there and
+  fits rate 1.3) and the fit fails closed exactly as on main (`success=False`),
+  without leaking the private exception that an earlier revision of this fix
+  raised there (PR #134 review). Unlike main, such a fit now emits the fallback
+  warning. It has the dedicated category
+  `liouscope.fitting.gls.AmplitudeRescalingFallbackWarning`, a `RuntimeWarning`
+  subclass, so a caller running with `-W error` can filter it by class. Raw residuals, rho, sigma and
   the likelihood remain in the caller's data units.
 - **The AR(1) lag-1 autocorrelation is formed from residuals normalised by an
   exact power of two (PR #134, round-3 review).** `ar1_correlation` took its dot
diff --git a/src/liouscope/fitting/gls.py b/src/liouscope/fitting/gls.py
index de8742a..38aedbc 100644
--- a/src/liouscope/fitting/gls.py
+++ b/src/liouscope/fitting/gls.py
@@ -68,6 +68,19 @@ class _ScaledResidualOverflowError(ArithmeticError):
     """
 
 
+class AmplitudeRescalingFallbackWarning(RuntimeWarning):
+    """The #124 residual rescaling was not representable; the fit ran unscaled.
+
+    A dedicated ``RuntimeWarning`` subclass (PR #134 review): for a free
+    amplitude parameter at ``max|y| >= ~1e150`` the fallback fires where the
+    pre-#124 code was silent, with a result identical to it. A caller running
+    with ``-W error`` can filter exactly this notice by class instead of all
+    RuntimeWarnings. It is emitted per fit, not once per process: which curve
+    lost the #124 invariance is audit information, and a process-global
+    "once" registry would make it depend on call order.
+    """
+
+
 def _whiten(y: np.ndarray, rho: float) -> np.ndarray:
     if y.size < 2:
         y_copy: np.ndarray = y.copy()
@@ -310,7 +323,7 @@ def _record(kind: str, _flag: int) -> None:
                     "optimiser's probe points; fitting the unscaled residuals "
                     "instead, so the amplitude-scale invariance of issue #124 "
                     "does not hold for this curve (PR #134).",
-                    RuntimeWarning,
+                    AmplitudeRescalingFallbackWarning,
                     stacklevel=2,
                 )
                 fit_scale = 1.0
diff --git a/tests/test_gls_amplitude_scale.py b/tests/test_gls_amplitude_scale.py
index 82c536c..555ea3b 100644
--- a/tests/test_gls_amplitude_scale.py
+++ b/tests/test_gls_amplitude_scale.py
@@ -257,3 +257,53 @@ def test_ar1_correlation_is_bit_identical_under_power_of_two_rescaling() -> None
     base = ar1_correlation(y)
     for k in (-900, -540, 500, 1000):
         assert ar1_correlation(np.ldexp(y, k)) == base, k
+
+
+# --------------------------------------------------------------------------
+# PR #134 round 3 (Equalita follow-up on 9a623e8)
+# --------------------------------------------------------------------------
+
+
+def _sinc_model(t: np.ndarray, params: np.ndarray) -> np.ndarray:
+    # Same values at every parameter; the 0/0 at t = 0 is only RAISED at the
+    # seed, i.e. inside the optimiser's first evaluation. The post-fit model
+    # evaluation (at ~1.3) is event-free, so it cannot mask what the solve did.
+    if params[0] == _SEED_RATE:
+        sinc = np.where(t > 0, np.sin(t) / t, 1.0)
+    else:
+        sinc = np.where(t > 0, np.sin(t) / np.where(t > 0, t, 1.0), 1.0)
+    return 1.0e-40 * np.exp(-params[0] * t) * sinc
+
+
+def test_model_runs_under_the_callers_floating_point_policy() -> None:
+    """The detector's errstate must not replace the caller's for the MODEL.
+
+    A caller that asks for ``invalid="raise"`` gets the model's 0/0 during the
+    solve as a ``FloatingPointError``, exactly as on main. Evaluating the
+    model under the detector's own errstate (ignore or call) swallows it.
+    """
+    y = _sinc_model(_T, np.array([_TRUE_RATE]))  # event-free away from the seed
+    with np.errstate(invalid="raise"):
+        with pytest.raises(FloatingPointError):
+            fit_gls_ar1(
+                _sinc_model, _T, y, np.array([_SEED_RATE]),
+                bounds=(np.array([0.0]), np.array([5.0])), n_iters=1,
+            )
+
+
+def test_fallback_warning_has_a_dedicated_filterable_category() -> None:
+    """``-W error`` users can silence exactly the fallback notice by class."""
+    import warnings
+
+    from liouscope.fitting.gls import AmplitudeRescalingFallbackWarning
+
+    assert issubclass(AmplitudeRescalingFallbackWarning, RuntimeWarning)
+    scale = 1.0e200
+    y = scale * np.exp(-_TRUE_RATE * _T)
+    with warnings.catch_warnings(record=True) as caught:
+        warnings.simplefilter("always")
+        fit = fit_gls_ar1(_free_amplitude, _T, y, np.array([scale, 0.2]))
+    fallback = [w for w in caught if _FALLBACK in str(w.message)]
+    assert fallback, [str(w.message) for w in caught]
+    assert all(w.category is AmplitudeRescalingFallbackWarning for w in fallback)
+    assert not fit.success