(fix): ND mixed-partial BC consistency — restore Clairaut symmetry#111
(fix): ND mixed-partial BC consistency — restore Clairaut symmetry#111
Conversation
The ND PreCompute build helpers `_get_effective_bc` (cubic) and `_get_effective_bc_quadratic` previously substituted `CubicFit()` / `Right(QuadraticFit())` for the user BC when building mixed partials (`p_src > 1`). This produced an axis-asymmetric tensor product and a numerical drift vs. the OnTheFly composition path (~1e-10 cubic, ~1e-5 quadratic for non-default BCs). PR #110 documented the discrepancy and relaxed four autodiff test tolerances; this commit removes the substitution at the root. Both helpers now return the user BC for every partial. By linearity and C^1 continuity of the 1D Hermite build, this makes the stored mixed partial exactly Σ_{k,l} φ_k'(x_i) ψ_l'(y_j) f[k,l], which is OnTheFly's composition formula, so PreCompute and OnTheFly become bit-equivalent (modulo a few ULPs of FP reordering noise) for every user BC. The short-grid defensive fallback (cubic: < 4 points, quadratic: < 3 points) is preserved but now emits an informative one-shot `@warn maxlog=1` via a `@noinline` helper before substituting the default. Side effects: - The cubic ND adjoint in `cubic_nd_adjoint.jl` retains its `caches`/`mixed_caches` dual-cache plumbing. After this fix the two are structurally equivalent for the common case (non-periodic, length >= 4); the redundancy is harmless and a follow-up PR can collapse them. An inline comment marks the opportunity. - The quadratic ND oneshot docstring previously carried a `!!! note "Strategy discrepancy"` explaining the ~1e-6 drift. That note is removed and replaced with a plain statement that the two strategies are bit-equivalent, with a pointer to `claudedocs/TODO/00_HIGHEST_PRIORITY_mixed_partial_bc_fix.md` for historical context. PeriodicBC is cubic-only in this codebase, so the quadratic helper does not need a Rule-2-equivalent propagation check.
New dedicated test file `test_nd_mixed_partial_bc_consistency.jl` pins the mathematical invariants that the BC fix in the preceding commit restores. Authored as the RED phase of a TDD cycle: on the pre-fix commit the math-invariant groups fail by the bug magnitudes documented in the design doc (~1e-10 cubic, ~1e-5 quadratic); after the fix they pass to `atol = 50*eps(Float64)` (~1.1e-14), which is generous enough to absorb tridiagonal / recurrence FP reordering but orders of magnitude tighter than any bug signature. Five test groups: A. PreCompute ↔ OnTheFly bit-equivalence across BCs (ZeroCurvBC, ZeroSlopeBC, CubicFit, MinCurvFit, Right(QuadraticFit)) and data functions (sin*cos, sin*y^3, sin*(y-π)) in 2D and 3D for both cubic and quadratic. B. Clairaut / axis-swap symmetry: building on `(xs, ys)` with `data` vs `(ys, xs)` with `permutedims(data)` and querying at swapped coordinates must give the same value. This is the strongest regression guard: the pre-fix PreCompute violated this by up to ~1e-5 for quadratic ZeroSlopeBC and ~1e-9 for cubic. C. Hessian off-diagonal symmetry (H[1,2] ≈ H[2,1]) for PreCompute interpolants with non-default BCs in both cubic and quadratic. Note: this assertion happens to pass on the buggy code because a single build is internally symmetric; we keep it as a forward regression guard. D. Cubic periodic-on-one-axis propagation regression-pin. Rule 2 in `_get_effective_bc` already handles this; the test asserts it continues to work after the fix. (Quadratic does not support PeriodicBC in this codebase.) E. Short-grid fallback warnings. Tests the helpers directly via `FastInterpolations._get_effective_bc` / `_get_effective_bc_quadratic` on very short grids (cubic: 3 points, quadratic: 2 points) with non-PolyFit user BCs. Asserts the `@warn ... maxlog=1` fires and the fallback BC is returned (`ZeroCurvBC()` / `MinCurvFit()`). Uses `@test_logs` for warning capture. Registered in `test/runtests.jl` next to `test_nd_oneshot_onthefly.jl`.
…d check Downstream verification that the mixed-partial BC consistency fix works end-to-end through the autodiff pipeline. PR #110 temporarily relaxed four assertions in the ND quadratic autodiff tests (two in `test_autodiff_Zygote.jl`, two in `test_autodiff_Enzyme.jl`) because Zygote/Enzyme trace the interpolant constructor (`PreCompute`) while the cross-validation seed comes from the one-shot API (now `OnTheFly` under `AutoCoeffs`), and the two paths disagreed by ~1e-6 for non-default quadratic BCs. With the BC consistency fix, PreCompute and OnTheFly produce bit-equivalent results, so the original `atol = 1.0e-10` is restored at all four sites and the TODO comment blocks are removed: - test/ext/test_autodiff_Zygote.jl - "Quadratic eval — L2 loss" : rtol=1e-4 → atol=1e-10 - "Quadratic gradient — ‖∇f‖² loss" : rtol=1e-3 → atol=1e-10 - test/ext/test_autodiff_Enzyme.jl - "Quadratic eval — L2 loss" : rtol=1e-4 → atol=1e-10 - "Quadratic gradient — ‖∇f‖² loss" : rtol=1e-3 → atol=1e-10 Also adds a new testset "AutoCoeffs: quadratic AD seed matches PreCompute exactly" in `test_nd_oneshot_onthefly.jl` that uses `==` (bit-exact equality, no tolerance) to compare OnTheFly and PreCompute for the exact call patterns autodiff uses — value query, ∂/∂x, ∂/∂y. This is a direct regression pin for the autodiff seed consistency and is stronger than the `atol = 50*eps` tolerance used in the new `test_nd_mixed_partial_bc_consistency.jl`. Verified locally: - test_autodiff_Zygote.jl passes (2m 29s) - test_autodiff_Enzyme.jl passes (68 pass, 2 broken pre-existing, 20m 49s) - test_nd_oneshot_onthefly.jl passes (1m 24s)
The full-suite run on the branch uncovered one hard failure and several pieces of documentation drift that encoded the pre-fix buggy behavior. All are fixed in this commit. **test/test_nd_coverage.jl** - `_get_effective_bc edge cases` (cubic, line 141): the test asserted `_get_effective_bc(ZeroCurvBC(), 2, grid_long) isa CubicFit` which was a characterization test on the now-removed substitution. Updated to assert `isa ZeroCurvBC` (the user BC is returned). - Cubic short-grid edge case (line 138): the fallback branch now emits an informative `@warn maxlog=1`. Wrapped the call in `@test_logs (:warn, r"Cubic ND mixed-partial build")` to assert the warning and suppress it from test output. - `_get_effective_bc_quadratic edge cases` (line 352): symmetric update — the quadratic short-grid branch also now emits `@warn maxlog=1`. Wrapped the call in `@test_logs (:warn, r"Quadratic ND mixed-partial build")`. **test/test_nd_oneshot_onthefly.jl** - "Quadratic OnTheFly with all supported BCs" (line 455): the long NOTE describing the ~1e-6 divergence as expected behavior was rewritten to describe the post-fix reality (PreCompute ≡ OnTheFly bit-equivalent for every user BC). The assertion tolerance is tightened from `rtol = 1e-5` to `atol = 50 * eps(Float64)` (~1.1e-14), matching the new regression suite. - "Quadratic BC-exact data" (line 479): the explanatory comment was rewritten as a historical note. The test is retained as a coarse smoke check; bit-identity now holds for arbitrary smooth data, not only BC-exact data, so the original motivation no longer applies. Verified locally: - test_nd_coverage.jl passes (30s) - test_nd_oneshot_onthefly.jl passes (1m 19s)
Adds a new test group F to `test_nd_mixed_partial_bc_consistency.jl` that directly invokes `_get_effective_bc` and `_get_effective_bc_quadratic` with inputs hitting every distinct return statement in both helpers. Groups A–C cover the helpers indirectly through the full ND build pipeline, but a few branches (e.g., `get_polyfit_degree(bc) > 0` short-circuit; `PeriodicBC` Rule 2 early return) are not trivially reached that way. Group F exists so `src/cubic/nd/cubic_nd_build.jl` and `src/quadratic/nd/ quadratic_nd_build.jl` land at ≥ 95% line coverage for the patch. Also consolidates group E's two `@test_logs` calls into single `match_mode = :any` blocks (warning assertion and return-value assertion together) so the short-grid fallback emits exactly one captured log per helper, with no stray stderr output from a second call outside `@test_logs`. Branches covered by group F: - Cubic `_get_effective_bc`: - Rule 1 (p_src == 1) for three BC types - Rule 2 (periodic) for long and short grids - Rule 3 left branch (`get_polyfit_degree > 0`) for both lengths - Rule 3 right branch (`length >= 4`) for ZeroSlope/ZeroCurv - Quadratic `_get_effective_bc_quadratic`: - Rule 1 for Right(QuadraticFit), MinCurvFit, ZeroCurv - Rule 2 (length >= 3) for three BC types The short-grid fallback path for both helpers is still covered by group E's `@test_logs` blocks — those calls execute every line of both the helper body and the `@noinline` warning helper.
- Restructure `_get_effective_bc_quadratic` to use explicit `if/end` for the length-≥3 branch, mirroring cubic's `_get_effective_bc` control flow so Julia coverage attribution lands on every executable line. - Add direct short-grid invocations to group F (cubic + quadratic). The group E `@test_logs` macro swallows line counters for the wrapped helper body (a known Julia coverage attribution quirk); the unwrapped direct calls give ground-truth line attribution. `maxlog=1` keeps the output silent since group E already fired the warning. Patch line coverage on this branch: 77/77 = 100.0%.
Source: - quadratic_nd_oneshot.jl: docstring lied about AutoCoeffs routing for scalar quadratic queries — said "OnTheFly" but the implementation unconditionally routes AutoCoeffs to PreCompute (intentional, for bit-exact AD seed agreement with the interpolant constructor). Updated the docstring to reflect actual behavior and explain why cubic differs. Also updated the design-doc reference to its archived location. Tests: - Move testset D (cubic periodic propagation) to its declared position between C and E, matching the lettered header comment. - Add a second 3D axis-swap permutation in group B (axes 1↔2 in addition to the existing 1↔3) for both cubic and quadratic. - Add a `!iszero` Hessian-diagonal sanity assertion to group C so a hypothetical bug zeroing the entire matrix would not silently pass. - Clarify the explanatory comment in group F: the second short-grid call is silent because there is no `@test_logs` wrapper, not because `maxlog=1` is exhausted (Test.jl bypasses the maxlog counter). All targeted tests still green: - test_nd_mixed_partial_bc_consistency.jl - test_nd_oneshot_onthefly.jl
…ion path Two related bugs in the GridIdx + coeffs interaction shipped with PR #110: 1. Scalar `interp` validated `coeffs` against the *un-promoted* method tuple, so `interp((x, y), data, (qx, GridIdx(k)); method=(CubicInterp(), PchipInterp()), coeffs=PreCompute())` rejected the call even though the GridIdx auto-promotion would have replaced PchipInterp with NoInterp, leaving a 1-D cubic problem that PreCompute fully supports. Fix: run `_promote_grididx_to_nointerp` BEFORE `_validate_nd_coeffs` so the validation sees the methods that actually drive dispatch. 2. Batch `interp!` delegating to `_interp_batch_with_grididx!` did not forward the `coeffs` kwarg, so explicit `coeffs=PreCompute()` / `coeffs=OnTheFly()` were silently dropped on the GridIdx batch path and the reduced sub-problem ran with the default AutoCoeffs fallback. Fix: add a `coeffs` kwarg to `_interp_batch_with_grididx!` and forward it at both delegation sites (early-return for fully-expanded queries and the main reduced-problem `interp!` call). Tests cover: - Scalar: GridIdx-on-Pchip + PreCompute now succeeds and matches the manually-sliced 1-D cubic call. The non-EvalValue deriv path (which disables promotion) still correctly rejects. - Batch: NoInterp axis is GridIdx-sliced → reduced 1-D Cubic problem accepts both PreCompute and OnTheFly. Differentiator: non-NoInterp GridIdx (Pchip remains in the method tuple) now correctly REJECTS `coeffs=PreCompute()` instead of silently dropping it. Both fixes verified in TDD red→green order. Touched suites all green: test_nointerp.jl, test_nd_mixed_partial_bc_consistency.jl, test_nd_oneshot_onthefly.jl, test_hetero_nd.jl, test_hetero_oneshot.jl, test_hetero_precomputed.jl.
FastInterpolations.jl BenchmarksAll benchmarks (42 total, click to expand)
This comment was automatically generated by Benchmark workflow. |
There was a problem hiding this comment.
Pull request overview
Fixes ND precomputed mixed-partial construction to apply the user-specified boundary conditions consistently for both pure and mixed partials, restoring axis symmetry (Clairaut) and aligning PreCompute() with OnTheFly() behavior. Also fixes GridIdx auto-promotion paths to correctly honor an explicit coeffs strategy.
Changes:
- Update cubic/quadratic ND mixed-partial BC selection to use the user BC for mixed partials (with a warned short-grid fallback).
- Ensure
coeffsis respected whenGridIdxtriggers method auto-promotion and in theGridIdxbatch (interp!) delegation path. - Add/adjust regression tests for PreCompute↔OnTheFly consistency, Hessian symmetry, and tightened AD tolerances.
Reviewed changes
Copilot reviewed 13 out of 13 changed files in this pull request and generated 3 comments.
Show a summary per file
| File | Description |
|---|---|
| test/test_nointerp.jl | Adds tests ensuring coeffs=PreCompute() works correctly after GridIdx auto-promotion and is forwarded through batch GridIdx paths. |
| test/test_nd_oneshot_onthefly.jl | Adds quadratic AutoCoeffs/PreCompute exact-equality pins and tightens quadratic OTF vs PreCompute expectations across BCs. |
| test/test_nd_mixed_partial_bc_consistency.jl | New regression suite for PreCompute↔OnTheFly agreement, axis-swap symmetry, Hessian symmetry, periodic propagation, and short-grid fallback warnings. |
| test/test_nd_coverage.jl | Updates BC edge-case expectations and asserts the new short-grid fallback warnings. |
| test/runtests.jl | Includes the new mixed-partial consistency regression file in the main test run. |
| test/ext/test_autodiff_Zygote.jl | Tightens tolerances now that quadratic OTF/PreCompute mismatch is addressed. |
| test/ext/test_autodiff_Enzyme.jl | Same tolerance tightening for Enzyme-based AD tests. |
| src/quadratic/nd/quadratic_nd_oneshot.jl | Changes quadratic AutoCoeffs() resolution (scalar oneshot) and updates docs to reflect new strategy behavior. |
| src/quadratic/nd/quadratic_nd_build.jl | Uses user BC for quadratic mixed partials and adds a one-shot warning on short-grid fallback to MinCurvFit(). |
| src/hetero/hetero_oneshot.jl | Reorders GridIdx promotion vs coeff strategy validation; forwards coeffs into GridIdx batch helper. |
| src/hetero/hetero_nointerp.jl | Extends _interp_batch_with_grididx! to accept/forward coeffs through its delegations. |
| src/cubic/nd/cubic_nd_build.jl | Uses user BC for cubic mixed partials (preserving periodic propagation) and adds a one-shot warning on short-grid fallback to ZeroCurvBC(). |
| src/cubic/nd/cubic_nd_adjoint.jl | Updates commentary to reflect the mixed-partial BC behavior change (no functional change). |
Comments suppressed due to low confidence (1)
src/quadratic/nd/quadratic_nd_oneshot.jl:166
- This change makes
AutoCoeffs()resolve toPreCompute()for scalarquadratic_interp(viacoeffs isa AutoCoeffs ? PreCompute() : coeffs), which is a user-visible behavior change vs the prior “scalar oneshot → OnTheFly” policy from #110. If intentional, consider calling it out explicitly in the PR description/changelog since it affects performance expectations and AD behavior for quadratic oneshots.
# Keep AutoCoeffs on the specialized PreCompute path so scalar one-shot
# calls match the interpolant constructor and AD rules.
coeffs_resolved = coeffs isa AutoCoeffs ? PreCompute() : coeffs
if coeffs_resolved isa OnTheFly
sample = @inbounds first(data)
methods = map(bc_i -> QuadraticInterp(_to_quadratic_bc(bc_i, sample)), bcs)
return _interp_nd_oneshot_onthefly(grids_typed, data, query, methods, extraps_val, searches, ops, hint)::Tr
end
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Copilot review on PR #111 flagged that docstrings/comments described PreCompute ↔ OnTheFly as "bit-equivalent" / "bit-identical" while the actual assertions use `isapprox(...; atol=50*eps(Float64))`. The wording overstates the guarantee — the residual difference is FP reordering noise in the tridiagonal/recurrence solver, ~50 ULPs at most, not zero. Reworded all four sites: - src/quadratic/nd/quadratic_nd_build.jl:233 — "bit-equivalence" → "numerical equivalence within FP noise" - src/cubic/nd/cubic_nd_build.jl:358 — short-grid warning text now says "no longer match the OnTheFly composition" instead of "not bit-equivalent" - test/test_nd_mixed_partial_bc_consistency.jl header + group-A label — "bit-equivalence" → "agreement to ~ULP / machine epsilon" - test/test_nd_oneshot_onthefly.jl:457 — "bit-equivalent results" → "agree to ~ULP tolerance" The single legitimate "bit-identical" usage at test_nd_oneshot_onthefly.jl lines 469-481 is preserved because that testset uses `==` (true bit equality) on data with `d²f/dy² = 0` everywhere, where both paths agree exactly. The "**not** bit-identical" wording in quadratic_nd_oneshot.jl:131 is also preserved — it correctly negates the property in the context of explaining why quadratic AutoCoeffs forces PreCompute. No code logic changed; comment-only commit.
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #111 +/- ##
=======================================
Coverage 96.56% 96.56%
=======================================
Files 136 136
Lines 11083 11090 +7
=======================================
+ Hits 10702 10709 +7
Misses 381 381
🚀 New features to boost your workflow:
|
* (feat): ND forwarders for pchip/cardinal/akima — unified `interp` entry Add ND-shape methods to `pchip_interp`/`cardinal_interp`/`akima_interp` (and their in-place `!` variants) that forward to `interp(grids, data; method=...)`. Covers all four ND call shapes: interpolant construction, scalar one-shot, batch one-shot, and in-place batch. Lets users keep the per-method entry point they already use in 1D instead of rewriting to `interp(...; method=PchipInterp())` by hand. `cardinal_interp` forwarders thread `tension` into `CardinalInterp(tension)`. `CubicHermiteInterp` is intentionally not forwarded — user-supplied per-axis slope arrays need a separate ND design. * test: ND forwarder coverage for pchip/cardinal/akima Smoke tests for the new `pchip_interp`/`cardinal_interp`/`akima_interp` ND forwarders. Verifies ULP-exact equivalence with the corresponding `interp(grids, data; method = ...)` calls across all four ND call shapes (build, scalar oneshot, batch oneshot, in-place batch) for 2D and 3D fixtures. Coverage matrix: - Equivalence (`===` / `==`) vs direct `interp` for every call shape - `@inferred` type stability on the scalar oneshot path - Cardinal `tension` kwarg passthrough (default vs 0.5, plus cross-check against `CardinalInterp(0.5)`) - Grid flavors: all-Vector, all-range, range×Vector, Vector×range — forwarder must preserve concrete grid tuple type into `interp` * Runic formatting * docs: address Copilot review on PR #114 — fix doc/code drift Two doc-only fixes flagged by Copilot review: 1. `src/hetero/local_hermite_nd_forward.jl`: header said "all three ND call shapes" but enumerated four (off-by-one). 2. `test/test_local_hermite_nd_forward.jl`: header claimed `===` for all four call shapes, but batch/array paths use `==` (array identity is meaningless). Clarified that scalar paths use `===`, array paths `==`. No code behavior change.
Summary
PreCompute()now applies the userbcuniformly to every partial (pure and mixed), restoring∂²f/∂x∂y == ∂²f/∂y∂x(Clairaut's identity) and makingPreCompute() ↔ OnTheFly()agree to ULP for every BC.coeffskwarg is now honored on the GridIdx auto-promotion paths (scalarinterpand batchinterp!).Background
_get_effective_bc(cubic) and_get_effective_bc_quadraticsubstitutedCubicFit()/Right(QuadraticFit())for the user BC when forming mixed partials (p_src > 1). This produced an axis-asymmetric tensor product:PreComputevsOnTheFlydrifted by ~1e-10The substitution was deliberate in the original implementation (a hidden BC override) but mathematically wrong: it broke the equivalence
between the OnTheFly composition formula and the PreCompute stored value. Returning the user BC for every partial restores the equivalence.
A short-grid defensive fallback (
length(grid) < 4cubic,< 3quadratic) remains in place. It now emits a one-shot@warn maxlog=1so users notice the substitution instead of silently losing their BC.Behavior change
For non-default BC + non-trivial data (e.g.
bc=ZeroCurvBC()on data whosed²f/dy²at the boundary is non-zero),PreCompute()results change by the bug magnitudes above. Code paths that pinned to the buggy numerics will need to regenerate baselines. Default BCs (Left(QuadraticFit()),CubicFit()) andOnTheFly()are unchanged.GridIdx coeffs fix (bonus)
Two related bugs in the
GridIdx+coeffsinteraction shipped with #110:Scalar
interp((x, y), data, (qx, GridIdx(k)); method=(CubicInterp(), PchipInterp()), coeffs=PreCompute())rejected the call even though the GridIdx auto-promotion would have replacedPchipInterpwithNoInterp, leaving a 1-D cubic problem that PreCompute supports. Fix: promote first, validate second.Batch
interp!delegating to_interp_batch_with_grididx!did not forwardcoeffs, so explicit strategies were silently dropped on the GridIdx batch path. Fix: forward the kwarg through both delegation sites.