fix(#124): make GLS convergence invariant to observable amplitude - #134
fix(#124): make GLS convergence invariant to observable amplitude#134marcohost33-maker wants to merge 10 commits into
Conversation
|
You have reached your Codex usage limits for code reviews. You can see your limits in the Codex usage dashboard. |
|
Exact-head verification is complete on |
Status 2026-09-07 — the head IS matrix-tested, but 13 commits staleCorrecting an assumption that applies to this stack generally: because Those runs exist because the now-closed exact-head verification PR #138 carried this same SHA against Staleness, measured ( It declares No issue band added: this PR targets What Marco has to decideRebase onto the current Read-only measurement. No push, no merge, no rebase, no re-run, nothing resolved. |
…nerate-bca fix(#125): fail closed on degenerate BCa intervals
… 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
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
![github-advanced-security[bot]](https://avatars.githubusercontent.com/in/57789?s=60&v=4){kind=small}
|
@codex review — please review 7ad4829..9a623e8: GLS residual scaling now falls back to the unscaled path (identical to main) when the scaled solve raises a floating-point event in /scale, whitening or SciPy arithmetic; model evaluation keeps the caller's float behaviour. Focus: amplitude invariance (#124) vs. unrepresentable-scale fits staying selectable (#147/#135), and large amplitudes (>=1e150). |
![chatgpt-codex-connector[bot]](https://avatars.githubusercontent.com/in/1144995?s=60&v=4){kind=small}
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 9a623e8e1c
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| 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) |
There was a problem hiding this comment.
When the default n_iters > 1 is used and fitted raw residuals are below roughly 1e-162 or above roughly 1e154, this normalization protects only least_squares; the subsequent ar1_correlation_corrected call still forms unscaled dot products in neff.py, which underflow to zero or overflow to NaN. The next iteration consequently whitens with a correction-floor rho or NaN, so amplitude-equivalent correlated-noise curves can produce different fitted rates or success states. Estimate rho from normalized residuals, or use scale-safe dot products, before the next iteration.
Useful? React with 👍 / 👎.
| - **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 |
There was a problem hiding this comment.
This changes the fitting methodology and can change fitted rates and model selectability, but the commit leaves CITATION.cff unchanged, so its pending-next-release methodology record omits this correction. Add a bounded description of the GLS scaling and fallback behavior to the citation metadata.
AGENTS.md reference: AGENTS.md:L121-L123
Useful? React with 👍 / 👎.
…t 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
…arning 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
|
@codex review — please re-review 9a623e8..28c7db2: AR(1) rho is now estimated from power-of-two-normalised residuals in neff.py (exact, bit-identical where the old computation stayed in range); CITATION.cff pending-methodology block; dedicated AmplitudeRescalingFallbackWarning; test pinning that user models run under the caller's float settings. Known and deliberately not in this push: estimate_neff_geyer/_autocorr has the same scaling class (N_eff = n at 1e-170/1e160). |
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 28c7db2c40
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| # 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]) |
There was a problem hiding this comment.
When callers enable strict underflow handling (for example, with np.errstate(under="raise")) and finite residuals span more than roughly 1074 binary exponents, this normalization underflows the smallest entries and raises FloatingPointError before computing rho. For example, ar1_correlation([1e150, 1e-200]) now raises under that policy, whereas the previous mean-centered calculation returns -0.5; the exception can also escape a multi-iteration GLS fit. Suppress underflow specifically around this implementation-detail scaling while retaining the caller's policy for other arithmetic.
Useful? React with 👍 / 👎.
2 participants
Summary
Closes #124 on top of PR #121 (
base=pr107-fix) as one isolated numerical-correctness slice.Defect
For a curve such as
1e-40 * exp(-1.3 t), the raw residuals and their numerical gradient are already far below SciPyleast_squarestermination scales at the starting point. TRF can therefore return the seed withsuccess=Trueafter one evaluation even though the same curve at order-one amplitude fits normally.That violates a basic invariance: multiplying an observable and its model by a positive constant must not change the fitted rate or whether a fit occurred.
Fix
The residual vector presented to
least_squaresis divided by the observed curve's own finite, non-zero magnitudemax(abs(y))after the existing #123 degeneracy guard.This is a positive constant scaling of the least-squares objective, so it does not change the mathematical minimizer. It only changes the numerical units seen by the optimizer. Raw residuals, AR(1) estimation, reported sigma, and likelihood remain in the caller's original data units.
No solver tolerances are loosened or tightened, and this deliberately does not attempt the separate time/parameter nondimensionalisation tracked in #111.
SciPy documents that
least_squaressuccess is tied toftol/xtol/gtol, with TRF's gradient criterion compared togtol; it separately describesx_scaleas a scaled-variable reformulation. The correction here addresses the residual-amplitude axis directly rather than retuning the solver.Reference: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html
Discrimination tests
tests/test_gls_amplitude_scale.pypins the exact acceptance pair from #124:exp(-1.3 t)at order-one amplitude must still converge to 1.3;1e-40 * exp(-1.3 t)must converge to the same rate and must move far away from the seed 0.2.The second assertion kills the pre-fix path, which returned the seed as a successful fit.
Scope boundary
N_eff(N_eff: Geyer's estimator is 2-4x optimistic on the uniform path #128).Change-control note
No CHANGELOG/CITATION edit is bundled into this stacked corrective slice to avoid another conflict in the already-large #121 stack. Before integration into
main, the results-changing numerical correction must be recorded once in the canonical[Unreleased]/pending-release notes together with the final CI evidence.