refactor(reg_field): one N-D tap-table engine for the field regulariser

[balbasty]

Agent:claude-fastfields-to-teeny

Closesfastfields/fastfields-lib#105. Part of #50 (phase 1 of 6). Part of the teeny-refactor umbrella (fastfields-lib#21).

Stacked PR. Base is claude/bump-teeny-pin-main (#54) so the diff shows only
this change. Retarget to teeny once #54 lands.

What changed

regularisers/field/{1,2,3}d.h (4,082 lines) hand-expanded one idea about
thirty times — D ∈ {1,2,3} × {absolute, membrane, bending} × {plain, RLS, JRLS}
× {matvec, diag, kernel}. All three files are replaced by
regularisers/field/nd.h, one engine generic in D, plus one new shared
primitive.

file	Δ
regularisers/field/1d.h	−1,059
regularisers/field/2d.h	−1,351
regularisers/field/3d.h	−1,672
regularisers/field/nd.h	+1,041
stap.h	+166
bounds.h	+134
cuda_switch.h, meta.h, field.h, CLAUDE.md	+46
net	−2,705

The stencil code proper goes 4,082 → 1,207 lines (−70%), and nd.h is about 40%
comment. All 33 Kernels<Config<D,C,…>> entry points keep their exact
signatures, so cpu-impl and cuda-impl compile untouched — the impl-layer
integration is phase 3.

The primitive

stap.h — per-axis boundary-folded tap tables. A stencil of reach R turns each
axis' -R..+R taps into (off, sgn, inb), and a corner tap is just the
componentwise combination of two axis taps (off_d + off_e, sgn_d · sgn_e).
That is valid precisely because the boundary fold is separable per axis, and that
separability is the whole reason one engine can serve every D. The header is
energy-agnostic on purpose: phase 2's flow/Lamé engine builds on the same thing.

It also carries the three contractions every energy shares — sdelta (the
difference-form read), sdiag (the exact diagonal), smag (an unsigned
companion-array read) — so matvec, diag and kernel all come from one tap
enumeration. A diag/matvec disagreement, which is the kernels#48/#49 bug class,
becomes structurally impossible rather than merely unlikely.

RLS/JRLS is one parameterisation

wmode ∈ {none, split, joint} — an optional weight accessor, per-channel (RLS)
or shared and hoisted (JRLS). wmode::none swaps the weight neighbourhood for an
empty type, so every weight access sits in a discarded if constexpr branch and
the plain stencils carry no trace of it. Six hand-written bodies become one.

bound::dyn<B> throughout

bounds.h gains bound::dyn<B>ported verbatim from main (so the two
tracks merge without a conflict), and the engine writes every boundary access
through it rather than utils<B> — preserving the FF_STATIC_BOUND_*
runtime-fallback mechanism of #42 exactly as #50 requires. Also added:
bound::supports_bending (Decision 2's predicate) and
bound::index_stays_inbounds.

Deliberate behaviour changes

Decision 1 — diag_* is now the exact matrix diagonal at every voxel:

diag = w0 + Σ_t w_t · ((tap t folded onto the centre ? sgn_t : 0) − 1)

The old w0 − Σ_t w_t·sgn_t agrees in the interior. It differs at every
boundary voxel under every condition — not only in fastfields/fastfields-cpu-lib#90's diag_bending
corner case. Two concrete 1-D membrane rows (k0 = absolute, k1 = −m·v):

at x=0	true operator row	exact diag	old diag
Zero	(k0 − 2k1)·f0 + k1·f1	k0 − 2k1	k0
DCT2	(k0 − k1)·f0 + k1·f1	k0 − k1	k0 − 2k1

The new engine's diagonal is checked against a brute-force
"matvec on unit vectors" reference at every voxel of every case (see below).

Decision 2 — bending boundary validation. The predicate lands here; the
throw belongs at the host dispatch entry, in the companion cpu-lib PR.

Bugs generic code cannot avoid fixing

1. kernels#39 — the RLS weight map was read out of bounds.wget used the
   unsigned cget(wgt, off), which dereferences unconditionally; under Zero,
   NoCheck and DST1 the folded index can leave [0, n). Now gated through
   smag. Demonstrated at a Zero boundary with an exactly-sized buffer:

   === OLD ===
   ==28476==ERROR: AddressSanitizer: heap-buffer-overflow on address 0x504000000048
   READ of size 8 at 0x504000000048 thread T0
   #0 ... ff::cpu::bound::get<double, long>(double const*, long) bounds.h:373
   fastfields/fastfields-kernels#2 ... operator() old_1d.h:557
   fastfields/fastfields-kernels#3 ... matvec_membrane_rls<> old_1d.h:559
   === NEW ===
   NEW: out=0.6 (no sanitizer report above => no out-of-bounds read)
   

2. The RLS kernel-table rescaling silently never ran for most D.
   make_kernel_membrane_rls (3D) and make_kernel_bending_rls (2D and 3D)
   called get_kernelsize_*_rls() with no argument, which defaults to the
   static channel count — i.e. -1 whenever channels are dynamic, which they
   always are through the impl layer. The loop bound was negative and the tables
   were left unscaled (2× / 4× too strong). Only the 1D path ever worked. Generic
   code cannot reproduce a per-D transcription slip.

3. The bending-RLS first-order coefficient dropped one corner weight and
   double-counted another for D ≥ 2. In 2D, m21 used w13 where
   self-adjointness requires w31 (and m23 vice versa). The correct form is

   coef(±1 along d) = (w1[d] − 2·w2[d])·(Wc + W(d,s))
   − 2·w2[d]·(W(d,−s) + W2(d,s))
   − Σ_{e≠d} w11[d][e] · Σ_{t=±1} ( W(e,t) + Wcorner(d,s,e,t) )
   

   which is symmetric under the row/column swap; the old one is not. Measured
   below.

Not fixed here: kernels#40 (the membrane term inside
make_kernel_bending_rls scaled by ¼ like the bending term instead of ½). It
lives in the weight table, not the stencil, so fixing it is an independent
behaviour change; it is flagged at its construction site with a note on what the
fix needs.

Verification

1. Old-vs-new element identity

A standalone probe instantiates both engines (the old bodies extracted verbatim
from the parent commit into a parallel namespace) and compares them over
D ∈ {1,2,3} × all 8 boundary conditions × C ∈ {1,2,3} × 2 domain sizes ×
every voxel × {absolute, membrane, bending} × {plain, RLS, JRLS} ×
{matvec, diag, kernel} × the weight tables.

checks: 218136, failures: 0, skipped(documented): 124200
max relative deviation: overall 2.66e-15, on the exact-match set 0
PASSED

Identical under clang++ and g++ at -O0, -O2 and -O3, and clean under
-fsanitize=address,undefined.

The plain matvec / kernel / make_kernel* paths are bit-identical —
"max deviation on the exact-match set: 0", with the comparison tolerance set to
exactly zero. That needed the tap walk to keep the original summation order
(order-major, − before +, higher axis outermost in a corner block) and to
share one multiply across a ± pair; both are now explicit invariants with a
comment saying why. The 2.7e-15 residual is confined to the paths whose formula
genuinely changed (the exact diagonal, and RLS/JRLS).

The skips are the three documented deviations: diag_* at boundary voxels, RLS
where the old code's value came from a heap over-read, and bending-RLS for
D ≥ 2 where the old coefficient is the asymmetric one.

2. The exact diagonal, independently

At every voxel of every case the new diag_* is compared against the matrix
diagonal assembled by applying the new matvec_* to unit vectors. Included in
the count above; passes to 1e-13.

3. Self-adjointness, measured not assumed

#50's Decision 2 asserts membrane is self-adjoint under all eight conditions and
only bending fails, under Replicate / DCT1 / DST1. A second probe assembles the
operator explicitly and reports its relative asymmetry, new and old:

D = 3 (new / old relative asymmetry)
D=3 Zero | membrane 0 / 0 | bending 0 / 0 | bend-RLS 2.3e-17 / 0.093
D=3 Replicate | membrane 0 / 0 | bending 0.06 / 0.06 | bend-RLS 0.075 / 0.075
D=3 DCT1 | membrane 0.29 / 0.29 | bending 0.42 / 0.42 | bend-RLS 0.45 / 0.45
D=3 DCT2 | membrane 0 / 0 | bending 0 / 0 | bend-RLS 6.2e-17 / 0.012
D=3 DST1 | membrane 0 / 0 | bending 0 / 0 | bend-RLS 2.3e-17 / 0.011
D=3 DST2 | membrane 0 / 0 | bending 0 / 0 | bend-RLS 3.1e-17 / 0.0066
D=3 DFT | membrane 0 / 0 | bending 0 / 0 | bend-RLS 2.2e-17 / 0.012
D=3 NoCheck | membrane 0 / 0 | bending 0 / 0 | bend-RLS 2.3e-17 / 0.093

Three things fall out, and two of them contradict the issue's premise — see
"Findings for review" below. The one that does not: the new bending-RLS operator
is symmetric to ~1e-17 wherever plain bending is, while the old one is off by
0.6%–9% even under DCT2/DFT. That is independent confirmation that the new
first-order coefficient (bug 3 above) is the right one — no reference
implementation needed, just the requirement that the operator be a symmetric
matrix.

4. Codegen — no lost compile-time folding

-O3, one wrapper per translation unit, whole-.text measurement:

case	new bytes / fp-mul / calls	old bytes / fp-mul / calls
1D membrane	761 / 2 / 0	751 / 2 / 0
1D bending	1197 / 3 / 0	1182 / 3 / 0
2D membrane	1250 / 3 / 0	1240 / 3 / 0
2D bending	3120 / 6 / 0	2100 / 6 / 1
3D membrane	1952 / 4 / 0	1486 / 4 / 1
3D bending	2742 / 9 / 0	3852 / 10 / 1

Floating-point multiply counts match exactly — one per stencil weight — which is
only possible if D, the tap loops and the if constexpr variant selection all
fold away. No residual runtime loop, and no residual call.

That last column is a fix, not a freebie: before adding FF_INLINE, clang left
stencil_matvec out of line and called it once per channel.inline is a
linkage keyword; where inlining is a codegen requirement rather than a hint the
code now says so, and cuda_switch.h gained the macro next to CUDEV.

5. The suite

make -C fastfields-cpu-lib clean test, both compilers:

test_distance 2352 | test_distance_mesh 4622 | test_distance_spline 704
test_posdef 5092 | test_pushpull 326 | test_reg_field 465 | test_reg_flow 854
test_reg_op 186 | test_resize 4130 | test_restrict 191 | test_splinc 4577

failures: 0, PASSED on every line — identical to before the change.

CUDA is compile-only in this environment (no nvcc, no GPU in CI); cuda-impl's
reg_field.h calls the same unchanged signatures, and phase 5 covers that path.

Findings for review

Two measurements disagree with #50's Decision 2 as written. I implemented the
decision as specified — it is settled, and it errs conservative — but they should
be on the record:

1. DCT1 is not self-adjoint at reach 1 either. The issue asks to confirm that
   absolute/membrane are exactly self-adjoint under every condition and only
   bending has a problem. They are not: membrane under DCT1 measures 0.29–0.47
   relative asymmetry (whole-sample reflection folds x−1 onto x+1 at x=0, so
   A[0][1] picks up the fold while A[1][0] does not). This is identical in the
   old code, so it is pre-existing, not introduced here — but it means the
   exclusion set is arguably too narrow, and DCT1 should perhaps be rejected for
   membrane as well.

2. DST1 bending measures exactly self-adjoint (0 asymmetry, every D), so its
   exclusion is over-restrictive for the field energy. The bounds.h
transpose() argument that motivates excluding DST1 is about the Lamé
   cross-coupling block, which is phase 2 — it may well still be needed there.

Neither changes what this PR does. Both are worth a decision before phase 2 fixes
the exclusion set in stone.

Judgment calls made while implementing

• bound::dyn<B> did not exist on this branch.[teeny] Regularisers: replace ~8,550 lines of hand-expanded per-D stencils with one N-D tap-table engine #50 says to write every
  boundary access through it, so it is ported verbatim from main. The
  FF_STATIC_BOUND_* macros and BoundVec that drive it from the dispatch
  layers are not pulled in — they are bounds: implement bound::type::Dynamic runtime fallback + build-time static/dynamic policy #42's own concern and phase 6's.
• diag_* uses unit strides to build its tap table. It never reads data, so
  it only needs to know "did this tap fold onto the centre", and a real stride of
  0 (a broadcast axis) would make that test lie.
• A corner's fold test is per axis, not off_d + off_e == 0 — two offsets
  can cancel without either folding.
• The plain path is not expressed as RLS with weights ≡ 1. It would be the
  prettier factoring, but kernels#40 makes the two genuinely different today, so
  unifying them would have smuggled that fix in. wmode::none is a separate
  if constexpr arm instead, and the reason is written down at reg_field: membrane penalty runs at half strength inside make_kernel_bending_rls #40's site.

---

Generated by Claude Code

[claude]

Agent:claude-fastfields-to-teeny

About the red ❌ on test-via-cpu-lib

That check is not reporting on this diff. The workflow hardcodes ref: main when it checks out the cpu-lib harness, then pins only the nested kernels submodule to the commit under test — so any PR based on the teeny branch is built by the pre-teeny-refactor cpu-lib and cpu-impl, which cannot compile teeny-branch kernels at all.

Running the workflow's own recipe against teeny's tip e579add (which predates this PR) fails identically:

./impl/kernels/distance/spline.h:3:10: fatal error: 'teeny/teeny.h' file not found

Filed as fastfields/fastfields-lib#104, deliberately not fixed here (unrelated concern, own branch).

The gate that was actually run

make -C fastfields-cpu-lib clean test against the dev-tree checkout, on both compilers, before and after:

test_distance 2352 | test_distance_mesh 4622 | test_distance_spline 704
test_posdef 5092 | test_pushpull 326 | test_reg_field 465
test_reg_flow 854 | test_reg_op 186 | test_resize 4130
test_restrict 191 | test_splinc 4577

failures: 0, PASSED on every line, exit=0 — 23,499 checks, identical counts to the pre-change baseline.

---

Generated by Claude Code

[balbasty]

Agent:claude-fastfields-to-teeny

Independent adversarial review — measurements reproduced from first principles

I re-derived the self-adjointness survey without reusing this PR's probe: a fresh harness assembles the operator matrix column-by-column (A e_j for every basis vector, per the kernels#43 method) against both this branch's engine and the parent commit's hand-expanded code (compiled side-by-side in a renamed namespace). Grids 7 / 5×6 / 4×5×6, voxel size (1, 1.2, 0.8), max|A−Aᵀ|/max|A|:

D	bound	membrane new/old	bending new/old
1	DCT1	0.500 / 0.500	0.500 / 0.500
2	DCT1	0.295 / 0.295	0.425 / 0.425
3	DCT1	0.240 / 0.240	0.376 / 0.376
1–3	Replicate	0 / 0	0.049–0.167
1–3	DST1	0 / 0	exactly 0 / 0
1–3	Zero, DCT2, DST2, DFT, NoCheck	0 / 0	0 / 0

Both contested findings confirmed. (1) DCT1 breaks self-adjointness at reach 1 (mechanism: at x=0 the −1 tap folds onto +1, so A[0][1] = 2·w1 while A[1][0] = w1 — verifiable by hand). Identical in the old code, so pre-existing. (2) DST1 bending is exactly self-adjoint for the field energy, every D, pure and mixed (its ±2 fold lands back on the centre — a diagonal entry — and its ±1 fold hits the sign-0 phantom node; nothing off-diagonal is ever unmatched).

One required fix before merge: the new supports_bending doc-comment in bounds.h states "every condition here is involutive at reach 1, so [absolute/membrane] are self-adjoint under all eight" and that DCT1/DST1's "reflected ±2 tap lands on a different voxel than the reverse fold does." Both statements are contradicted by this PR's own Findings-for-review section and by my measurements (DCT1 fails at reach 1; DST1's fold is harmless). The code is fine; the in-tree rationale is not, and it is exactly the text a future reader will trust when revisiting the exclusion set.

Also independently verified, all pass:

• Bit-exactness: plain matvec/kernel/tables old-vs-new identical at every voxel (boundary voxels included), all 8 bounds × D 1–3 × both energies — 4,288/4,288 comparisons equal to tolerance 0. Membrane-RLS/JRLS and 1D bending-RLS agree to ≤2.3e-15.
• Exact diagonal: new diag_* matches the assembled matrix diagonal to 5.7e-14 across every D × bound × energy (old diag off by up to 221 absolute at boundaries, as expected under Decision 1).
• kernels#39: reproduced the heap over-read in the oldmatvec_membrane_rls under ASan with an exactly-sized buffer (bounds.h:318 via wget); new engine clean on the same input.
• RLS rescale bug: reproduced with a dynamic-C instantiation — old 3D membrane-RLS table 2× too strong (−1.3 vs −0.65), old 2D/3D bending-RLS 4× (−18.1 vs −4.525); 1D and 2D-membrane unaffected, exactly as described.
• Bending-RLS coefficient fix: old assembled operator asymmetric 2.1–4.1% for D≥2 even under Zero/DCT2/DST2/DFT/NoCheck; new ≤9e-17. The symmetry requirement alone confirms the new coefficient.
• Codegen: -O3 3D-bending-DCT2, dynamic C: clang 0 calls / 9 double-multiplies (old: 1 call, 10), g++ 9 multiplies, no residual per-tap loop. FF_INLINE is doing its job.
• bound::dyn: byte-identical to main's (diffed directly); no FF_STATIC_BOUND_*/BoundVec leakage beyond a comment.
• Full cpu-lib suite green on clang++ and g++, counts matching the PR body.

Recommendation posted with the review verdict on the companion cpu-lib PR: this PR is sound; fix the bounds.h comment, then merge (after #54, retargeted to teeny). The exclusion-set decision (add DCT1+membrane, drop DST1 for field bending?) belongs on #50 before cpu-lib#54 lands, not here.

---

Generated by Claude Code

---

Generated by Claude Code

[claude]

Agent:claude-fastfields-to-teeny

Correction to this PR's verification claim (self-reported)

This PR's body says the cpu-lib suite passed on bothclang++ and g++. The g++ half of that was not measuring anything, and I want it on the record rather than buried.

make clean in fastfields-cpu-lib only removes build/*.o and build/*.so — notbuild/testobj/*.o or the build/test_* binaries. Dependency tracking means a header change still forces a rebuild, so the clang++ numbers were genuine. But make clean && make test CXX=g++ right after a clang run found every binary up to date and simply re-ran the clang-built ones, printing no compile lines and exiting 0. And make test CXX=g++ cannot actually work here anyway — CXXFLAGS hardcodes -ferror-limit=1, which gcc rejects.

Filed as fastfields-cpu-lib#55.

Re-verified properly

From-scratch rm -rf build rebuilds, one per compiler, with the clang-only flags dropped for the gcc run — 30 real compile invocations each, confirmed by counting them in the logs:

suite	clang++	g++
test_distance	2352	2352
test_distance_mesh	4622	4622
test_distance_spline	704	704
test_posdef	5092	5092
test_pushpull	326	326
test_reg_field	525	525
test_reg_flow	854	854
test_reg_op	186	186
test_resize	4130	4130
test_restrict	191	191
test_splinc	4577	4577

failures: 0 / PASSED throughout, exit=0 both. So the conclusion this PR reported is unchanged — the engine is clean on both compilers — but it is now actually measured on both.

What was not affected: every g++ claim in this PR that came from a direct compiler invocation rather than through the cpu-lib Makefile — the element-identity probe (run at -O0/-O2/-O3 on both), the self-adjointness survey, the ASan/UBSan run (built with g++), and the codegen comparison table (measured separately per compiler). Those were real.

---

Generated by Claude Code