FastInterpolations 0.4.15 -> 0.4.17 changes b_res by 25% and flips the sign of the toroidal torque

[matt-pharr]

Bumping FastInterpolations from v0.4.15 to v0.4.17 — with no other change — moves the resonant area-weighted field b_res by ~25% and flips the sign of the perturbed-equilibrium toroidal torque.

This is very likely an error in FastInterpolations.

Reproduction

Single-variable test: Pkg.add(name="FastInterpolations", version="0.4.15") against an otherwise-current environment (DiffEqCallbacks 4.19.0, DoubleFloats 1.11.2, LinearAlgebra/SparseArrays 1.12.0, Julia 1.12.6). Case is examples/DIIID-like_ideal_example, n=1, run through the regression harness against the same GPEC commit on both sides.

Fixed grid (grid_type = "ldp", mpsi = 256)

quantity	0.4.17	0.4.15	diff
resonant area-weighted field b^r	[5 elem]	[5 elem]	24.62%
||resonant area-weighted field||	4.382144e-04	3.918227e-04	10.59%
Chirikov parameter	[5 elem]	[5 elem]	8.81%
island half-widths	[5 elem]	[5 elem]	3.94%
ODE steps (saved)	1481	1466	1.01%
PE toroidal torque	-5.060621e-02	-5.060621e-02	1.3e-11

The equilibrium is untouched here: q profile and pressure profile checksums are identical, and plasma volume, li1/2/3, beta_t, beta_n agree to 1e-15..1e-16. Only quantities downstream of the ODE integration move.

Auto grid (grid_type = "auto", mpsi = 0)

quantity	0.4.17	0.4.15	diff
PE toroidal torque	-5.034315e-02	+5.057643e-02	sign flip
q profile (checksum)	0fead70cb422…	71b97936ef0b…	different
pressure profile (checksum)	d553d38a470b…	958ae8fd02b7…	different
resonant area-weighted field b^r	[5 elem]	[5 elem]	4.45%
ODE steps (total)	2002	1933	3.45%

On the auto grid the two-pass measured-curvature knot placement itself changes, so the q and pressure profiles differ — the same equilibrium sampled on a different grid (plasma volume still agrees to 4e-11, singular ψ locations to 2e-14).

Notes

• Both effects are reproduced by pinning FastInterpolations alone; DiffEqCallbacks, DoubleFloats and the stdlib versions were held fixed and are not implicated.
• delta prime from the Riccati BVP path (asymptotic Frobenius basis) is unaffected (~1e-8 relative). The affected quantities — b^r, Chirikov, island half-widths — are all scalar multiples of the same resonant-surface jump computed in SingularCoupling.jl, which obtains ξ′ from a chord slope across saved ODE nodes. A shift in node placement therefore propagates into all of them.
• Numbers above are from the regression harness; the auto-grid and ldp runs used the same GPEC commit on both sides, so the only difference is the package version.

[matt-pharr]

@logan-nc @jhalpern30 @d-burg this is very concerning. How could a minor bump change things so much?

[jhalpern30]

I dug into this. Summary up front: the FastInterpolations delta is 1 ULP, and it is not the bug. The 18–25% is our own amplification in SingularCoupling.jl. Details and the controls below.

1. Which version actually moves the numbers

I ran the ldp deck (grid_type="ldp", mpsi=256, so the equilibrium is held fixed) three times, one isolated env per pinned FastInterpolations version, same GPEC commit throughout:

pair	b^r	island half-widths	Chirikov	Δ' (chord-slope)
0.4.15 → 0.4.16	18.23%	4.51%	5.87%	18.23%
0.4.16 → 0.4.17	0 — bit-identical	0	0	0
0.4.15 → 0.4.17	18.23%	4.51%	5.87%	18.23%

The whole effect enters at 15→16. 16→17 is bit-for-bit identical in a full run, so anything in the 0.4.17 release is off the hook for this case.

2. What changed in 15→16

I bisected numerically rather than by reading the ~5000-line refactor diff: three pinned scratch envs, ~530 bit-exact repr(Float64) values covering every FastInterpolations surface we touch — scalar cubic value/d1/d2/d3, vectorized eval, hint-based forward/reverse/stale sweeps and hint index traces, multi-series in-place + complex series, cumulative_integrate, the 2D rzphi bicubic (bc=(CubicFit(),PeriodicBC()), extrap=(ExtendExtrap(),WrapExtrap())) with all five mixed derivatives and 2D hints, theta wrap/OOB, nodal partials, gradient, every search and store policy.

Across the entire 0.4.15→0.4.17 range there are exactly two differences, both ~1 ULP:

• 15→16: one 1-ULP value in the multi-series scalar+hint path (1 query point out of ~200). This comes from the lean-anchor refactor — src/core/idx_stencil.jl deleted in favor of core/axis_indices.jl + core/series_lean_anchors.jl + cubic/cubic_series_payloads.jl, which moved multi-series scalar evaluation off the old _CubicAnchoredQuery/_eval_series_point! route. The weight formula (_compute_anchor_weights) and the muladd nesting are byte-identical, so this is FMA/inlining reassociation, not an arithmetic change.

  This is exactly the path our ODE drives: fmats/amats/locstab evaluated at ψ through the shared ffit_hint / spline_hint. That is the seed.

• 16→17: integrate (full-domain and limited-range) shifts ~1 ULP because @simd was added to the accumulation loops (integrate_common_1d.jl:146, integrate_fulldomain.jl). This relaxes FP-reduction associativity → multiple partial accumulators. It is deliberate and documented upstream, and it is an accuracy improvement — on a case with a known exact answer (∫(sin t + 0.4cos3t + 0.1)dt over [0,2π] = 0.6283185307179587) 0.4.16 gives …58, 0.4.17 gives …584. We consume it at Equilibrium.jl:319, but only for volume/li/beta, which is why those agree to 1e-15 and nothing downstream moves.

Everything else in the list above is bit-identical across all three versions.

3. Why 1 ULP becomes 18%

Control experiment — identical code, identical FastInterpolations, same ldp deck, changing only eulerlagrange_tolerance from 1e-10 to 1.000000001e-10 (one part in 10⁹):

resonant area-weighted field b^r    15.83 %   (surface 4: 90.58 %)
island half-widths                  10.51 %   (surface 4: 37.81 %)
Chirikov parameter                  37.81 %
Δ' (chord-slope)                     5.11 %   (surface 4: 90.58 %)

A perturbation with no physical content reproduces the same order of effect as the version bump. The equilibrium is bit-identical; the runs first diverge in ODE step acceptance (W inverse matrix was non-Hermitian ... at 19 vs 21 integration step(s)).

Mechanism, SingularCoupling.jl:18-29: _psi_bracket picks the two adjacent saved adaptive ODE nodes straddling lpsi/rpsi = psi_s ∓ 5e-4/(|n||q'|), and _hermite_cubic_val interpolates between them. Near a singular surface u varies nearly singularly, so which node pair happens to bracket that point dominates the answer — and b = singfac * xsp (line 257) is a tiny×huge near-cancellation on top. Any ULP anywhere reshuffles step acceptance → different node pair → tens of percent.

The corroboration is already in the original report: the Riccati-BVP Δ' is stable to ~1e-8 precisely because it uses the asymptotic Frobenius basis and never touches saved nodes.

4. What I think we should do

1. Don't pin 0.4.15. It would hide the symptom. 1 ULP in one evaluation path is within any reasonable numerical contract, and the next unrelated change (a compiler bump, a loop reorder, a thread-count change) will reproduce this.
2. Fix the diagnostic. Take ξ′ from the singular-surface asymptotic solution the way the Riccati Δ' path already does, or force saved ODE nodes at lpsi/rpsi so the bracket is deterministic rather than whatever the controller happened to accept.
3. Adjust the regression thresholds meanwhile. b^r, chirikov_parameter, island_half_width and the chord-slope delta_prime are currently at noise_threshold = 1e-8 in diiid_n1.toml. Their real reproducibility floor is tens of percent, so they are presently reporting noise as regressions. Either raise the thresholds to reflect that or drop them as tracked quantities until (2) lands.

5. Still open — probably a separate issue

I did not reproduce or explain the auto-grid toroidal-torque sign flip (−5.034e-02 → +5.058e-02). Note the magnitude is preserved to 0.5%, which does not look like cancellation noise — it looks like a discrete sign or ordering selection flipping. I tested only the ldp grid, deliberately, to hold the equilibrium fixed.

The auto grid also adds a second, discrete amplifier that the ldp grid doesn't have: _density_from_curvature! in GridRefinement.jl has a d4n < noise_floor && continue and a clamp(…, H_TARGET_MIN, H_TARGET_MAX), so one ULP in the curvature measure can add or drop a whole knot and change the entire pass-2 grid. That is consistent with the q/pressure checksums differing there. Worth its own issue rather than folding into this one.

Reproduction

Versions are already in the depot as ~/.julia/packages/FastInterpolations/{DBh2Y,66hkg,i7S5w} = 0.4.15/0.4.16/0.4.17. For the full-run comparison, make one scratch env per version so the repo Manifest.toml is untouched:

using Pkg
Pkg.activate("/tmp/gpec15"); Pkg.develop(path="/path/to/JPEC")
Pkg.add(name="FastInterpolations", version="0.4.15")

then run the same deck under each and diff perturbed_equilibrium/singular_coupling/* in gpec.h5. The API-level bit-exact probe is a standalone script; happy to push it somewhere if it's useful, though note integrate(x, y; method) and coeffs= on ND only exist in 0.4.17, so a shared probe has to avoid them.

[jhalpern30]

@matt-pharr I had Claude look into this, see above. It seems to not be an issue with FastInterpolations, but rather an extreme sensitivity of the resonant values to bit-wise changes in splines. Was this sensitivity not in the Fortran? I thought I remember Nik saying at one point you had done a whole scan on the jump width calculation

[matt-pharr]

@jhalpern30 @logan-nc Looking into this, for some reason the decision was made to use the slope of the two xi points nearest to the singularity, instead of using the derivative calculated by the ODE solver? That seems to be the reason for this. Going back to using the ode solver derivative (like fortran) eliminates the problem. I will open a pull request shortly.

[logan-nc]

there was a period of time where the Julia derivative solution wasn't trustworthy. maybe it was from that? I think the the sing spot normalization was also in flux for a while. it should all be good now? so matching Fortran is a good idea 😅