Uh oh!
There was an error while loading. Please reload this page.
fix(regularisers): correct diag_bending/diag_all corner cross-term sign error - #49
fix(regularisers): correct diag_bending/diag_all corner cross-term sign error#49balbasty wants to merge 1 commit into
Conversation
The boundary-corrected diagonal for diag_bending (field) and
diag_bending/diag_all (flow) computed the corner weight as
w*(fx0*fy0 + fx1*fy0 + fx1*fy0 + fx1*fy1) -- fx1*fy0 counted twice and
fx0*fy1 dropped. The correct expansion of (fx0+fx1)*(fy0+fy1) is
fx0*fy0 + fx1*fy0 + fx0*fy1 + fx1*fy1, matching what matvec_bending/
matvec_all already compute a few lines up in the same files.
The bug cancels at fully-symmetric corners (both axes flip the same way)
and only manifests where the two axes' one-sided boundary signs differ,
so it was invisible to the existing interior-only diag tests. It
corrupts the shipped field_diag/flow_diag output at boundary voxels
under sign-flipping bounds (Zero/DST1/DST2) whenever bending is active.
Fixes 6 functions across field/{2,3}d.h and flow/{2,3}d.h (12 wrong
terms total): field/2d diag_bending (1), field/3d diag_bending (3),
flow/2d diag_bending + diag_all (1+1), flow/3d diag_bending + diag_all
(3+3).
Closes#48
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015xcQBMEdA5eQtTsx2wB3AC…er (#56) * chore(deps): bump teeny submodule to main (5c46bd4) The pin sat at dcd591c, the tip of teeny's stale `claude/fastfields-teeny-refactor-js42id` integration branch, which is 125 commits behind teeny's `main`. That is now the branch this project tracks, so point `.gitmodules` at `main` too rather than at a branch that has not moved since the refactor started. Brings in, among 125 commits: - dispatch_values, product-form runtime->static dispatch (teeny#465) - the forward+backward scan_ sweep idiom docs (teeny#466) - scan_/scan (teeny#254), unfold (#256), subsample (#258), peel_zip, index_select, in-place maximum_/minimum_, sqdist/dist, host-side atomic fetch_add_/sub_ (#257) - the breaking rename take_along -> slice_along (teeny#423) The rename is the only breaking change that could reach us; a tree-wide grep for `take_along` across all fastfields repos finds no call sites (only two prose mentions of unrelated names in TEENIFICATION_REVIEW.md), so no downstream fixups are needed. Verified: `make -C fastfields-cpu-lib clean test` with CXX=clang++ and CXX=g++, before and after the bump. All 11 suites PASS with byte-identical check counts (2352/4622/704/5092/326/465/854/186/4130/191/4577). * refactor(reg_field): one N-D tap-table engine for the field regulariser `regularisers/field/{1,2,3}d.h` hand-expanded the same stencil roughly 30 times: D in {1,2,3} x {absolute, membrane, bending} x {plain, RLS, JRLS} x {matvec, diag, kernel}. That transcription has already shipped several bugs (kernels#39/#40/#41/#48/#49), and it is where the flow port would have had to repeat the exercise again. Replace all three files with `regularisers/field/nd.h`, one engine generic in D, built on a new shared primitive: * `stap.h` -- per-axis boundary-folded tap tables. A corner tap is the componentwise combination of two axis taps, valid because the boundary fold is separable per axis; that separability is what makes an N-D-generic engine possible at all. Energy-agnostic, so #50 phase 2's flow/Lame engine builds on the same header. * `bounds.h` -- `bound::dyn<B>` ported verbatim from `main` so every boundary access goes through it rather than `utils<B>`, preserving the FF_STATIC_BOUND_* runtime-fallback mechanism (kernels#42), plus the `supports_bending` and `index_stays_inbounds` predicates. * `cuda_switch.h` -- `FF_INLINE`. Without it clang leaves `stencil_matvec` out of line and calls it once per channel, which would have been a real regression against the old straight-line bodies. RLS and JRLS are now one parameterisation (`wmode`: an optional weight accessor, per-channel or shared) rather than six hand-written copies, and `wmode::none` swaps in an empty type so the plain stencils carry no trace of it. Two deliberate behaviour changes, both from #50's settled decisions: * `diag_*` returns the EXACT matrix diagonal at every voxel (decision 1): `w0 + sum_t w_t*((folded onto centre ? sgn_t : 0) - 1)`. The old `w0 - sum_t w_t*sgn_t` agrees in the interior and differs at every boundary voxel under every condition, not only in #41's diag_bending corner case. * The bending self-adjointness predicate lands here; the throw belongs at the host dispatch entry (decision 2), in a companion cpu-lib PR. and three bug fixes that generic code cannot avoid making: * the RLS weight map is no longer read out of bounds (kernels#39) -- demonstrated under ASan; * the RLS kernel-table rescaling now runs for every D (it sized itself from the STATIC channel count, so it silently never ran for 2D-bending / 3D-membrane / 3D-bending); * the bending-RLS first-order coefficient no longer drops one corner weight and double-counts another for D >= 2, which had made the old operator asymmetric by 0.6%-9% even under DCT2/DFT. kernels#40 is deliberately NOT fixed -- it lives in the weight table, not the stencil -- and is flagged at its construction site. Verified with a standalone old-vs-new element-identity probe over D x 8 bounds x C x every voxel x every variant (218,136 comparisons, 0 failures, plain matvec/kernel/tables BIT-IDENTICAL, worst deviation elsewhere 2.7e-15), an explicit assembled-matrix self-adjointness survey, and `make -C fastfields-cpu-lib test` on clang++ and g++ -- all 11 suites PASS with unchanged check counts. Closes #55. Part of #50. * docs(bounds): correct supports_bending's factual claims per review The self-adjointness rationale asserted reach-1 energies (absolute, membrane) are exact under all eight boundary conditions, and that DST1 breaks reach-2 (bending) via a mis-landed ±2 fold. Both were measured false during fastfields-kernels#56's review: membrane is NOT self-adjoint under DCT1 (0.29-0.47 relative asymmetry, same mechanism as bending's DCT1 failure), and DST1 is EXACTLY self-adjoint for field bending at every D (its ±2 fold lands back on the centre voxel, not on a mismatched one). supports_bending()'s return value is unchanged in this commit -- nothing calls it yet (that lands with the corresponding cpu-lib dispatch-entry check), so widening/narrowing its actual rejection set is a follow-up to fastfields-kernels#50's Decision 2, not this fix. This commit only stops the comment from documenting something false. --------- Co-authored-by: Claude <noreply@anthropic.com>
balbasty
commented
Aug 18, 2026
Closing this — not because the analysis was wrong. It wasn't: I re-derived the defect independently and this PR's diagnosis is correct in every particular. Why it's being closedThis PR comes from the The analysis was sound — credit where it's dueVerified independently rather than taken on trust:
Where the fix actually landedThe equivalent fix is now on The pin chain has been cascaded so the fix is live on both backends: Regression coverageThe On Nothing here has been dropped on the floor. This PR's branch has not been modified. Generated by Claude Code |
Agent: claude-fastfields-to-teeny
Closes#48
The bug
The boundary-corrected diagonal in
diag_bending(field) anddiag_bending/diag_all(flow) computed each corner weight asfx1*fy0is counted twice andfx0*fy1is dropped. The correct expansion of(fx0+fx1)*(fy0+fy1)iswhich is exactly what
matvec_bending/matvec_allalready compute a few lines up in the same files, e.g.regularisers/field/2d.h:+ w11 * (get(x0+y0, fx0*fy0) + get(x1+y0, fx1*fy0) + get(x0+y1, fx0*fy1) + get(x1+y1, fx1*fy1))The error cancels at fully-symmetric corners (both axes flipping the same way) and only shows up where the two axes' one-sided boundary signs differ — which is why the existing interior-only diag tests never saw it. It corrupts the shipped
field_diag/flow_diagpublic output at boundary voxels under sign-flipping bounds (Zero / DST1 / DST2) whenever bending is active, and hence also the Jacobi preconditioner built from that diagonal.Sites fixed — 12 terms, matching the issue exactly
regularisers/field/2d.hdiag_bendingw11xy)regularisers/field/3d.hdiag_bendingw110xy,w101xz,w011yz)regularisers/flow/2d.hdiag_bending,diag_allregularisers/flow/3d.hdiag_bending,diag_allTotal 12 wrong terms across 6 functions — the exact count the issue predicted, and the exact same per-file line counts as the already-verified fix on
claude/fastfields-teeny-refactor-js42id(5ed785a, ex-PR #32). Confirmed by grep that no remaining doubled-term expression exists anywhere underregularisers/.This is a re-application, not a cherry-pick:
mainhas since converted these fromstaticmethods usingImpl::qualification to instance methods on abound::dyn<B>-holding object (bound_utils_x.sign(...)), per kernels #42 / thebound::type::Dynamicwork. The corrected expressions are textually identical to the verified fix; only their surroundings differ.Verification
Built
fastfields-cpu-lib@origin/main(6aaa159) against this branch via the dev-tree symlinks (cpu-lib/impl -> cpu-impl@b6eae48 -> kernels@this),make test CXX=clang++:test_reg_field4610 checks / 0 fail,test_reg_flow5859 / 0No existing test changes value, i.e. nothing on
maincurrently covers the buggy boundary corners.Regression test — fastfields/fastfields-cpu-lib#53
Since
fastfields-kernelshas no test suite of its own, the regression test for this fix lives downstream: fastfields/fastfields-cpu-lib#53 carries the three boundary-symmetry tests (test-only, no production code). Verified revert/restore cycle against that branch:test_reg_fieldtest_reg_flow15295c2— pre-fix9b6517f— this PR's headAll 20 pre-fix failures show a uniform ~8.0 discrepancy against a
1e-5tolerance, so it is a decisive detector. The invariant needs no reference implementation: on a square domain with the same boundary condition on every axis (DST2, bending active), the operator's diagonal must be symmetric under axis swap, sodiag(0,j,c) == diag(j,0,c)— and fordiag_allwithshears == div,diag(0,j,c) == diag(j,0,1-c)under simultaneous axis + channel swap.Merge this PR first. cpu-lib#53 stays red until this lands and the
cpu-impl -> kernels/cpu-lib -> implpins are bumped; that ordering is spelled out in its body.