[teeny] Regularisers: replace ~8,550 lines of hand-expanded per-D stencils with one N-D tap-table engine

[balbasty]

Agent: claude-fastfields-to-teeny

Part of the teeny-refactor umbrella (fastfields-lib#21). This was the largest remaining piece of that effort — everything else (distance, posdef, pushpull forward, resize/restrict/splinc, mesh) was already teenified; the regulariser stencil bodies were not.

STATUS: ALL SIX PHASES DONE. Closing.

Final summary

phase	what	landed
1	Field engine (matvec/diag/kernel, plain + RLS/JRLS)	kernels#54, #56
2	Flow engine + Lamé cross-coupling	kernels#60, #61, cpu-lib#56/#57, cpu-impl#42 (closed #59)
3	Impl-layer spatial-index decode via peel_front<-1>().enumerate()	cpu-impl#46 (closed #43)
4	Relax: predicate-skip → dense subsample colour views	cpu-impl#50 (closed #49)
5	CUDA launcher consolidation (6→29 field, 6→25 flow entry points)	cuda-impl#30, #31 (closed #29)
6	Dispatch: macro pyramid → teeny's dispatch_values	cpu-lib#70, cuda-lib#39 (closed cpu-lib#69)

regularisers/{field,flow}/{1,2,3}d.h — 8,550 lines of hand-expanded per-D stencils, the source of at least six shipped bugs before this effort started — is gone, replaced by stap.h + regularisers/stencil.h + field/nd.h + flow/nd.h (~1,500 lines combined, one engine generic in D), a teeny-idiomatic impl layer, and a single shared dispatch header. Every phase was implemented by an Opus agent and independently re-verified by an adversarial Fable review before merge — nothing landed on the strength of the implementer's own claims alone.

Decisions this effort settled (see full history in this issue's edit log / linked PRs for the measurements behind each)

1. diag_* computes the exact matrix diagonal, not an approximate preconditioner diagonal.
2. Self-adjointness boundary rejection is measured, never argued, per-mechanism: bound::supports_reach(b, reach) for absolute/membrane/bending (both engines — flow's table happened to match field's, confirmed by independent measurement, not inherited by assumption), and bound::supports_lame/supports_lame_bending for Lamé's structurally different cross-channel block. Both static_assert-backed in bounds.h so the code can't drift from the measured table.
3. Coloured relax's visit order within a colour is a non-issue (independence guarantees it), proven by exhaustive partition-equivalence testing — but the effort surfaced that the colouring's own precondition (same-colour voxels share no stencil neighbourhood) can be false under a periodic boundary with an incompatible axis length, a genuine pre-existing data race (kernels#66), preserved unchanged by both the CPU and CUDA relax ports since fixing it is a separate, deliberate design decision.

Bugs found along the way, tracked separately (none blocked any phase; several predate this whole effort)

• kernels#39, reg_field: membrane penalty runs at half strength inside make_kernel_bending_rls #40, fastfields/fastfields-cpu-lib#90 (pre-existing on main, useful spec inputs — fastfields/fastfields-cpu-lib#90 and reg_field/reg_flow: out-of-bounds weight-map read under Zero boundary #39 resolved as side effects of the rewrite; reg_field: membrane penalty runs at half strength inside make_kernel_bending_rls #40 deliberately not fixed, a weight-table issue)
• kernels#48/fix(regularisers): correct diag_bending/diag_all corner cross-term sign error #49, cpu-lib#51/regularisers: fix OOB weight-map read in RLS/JRLS wget lambdas (Zero bound) #52/#53 (the diag corner cross-term bug that motivated re-examining the whole rejection table)
• kernels#62 (3-D matvec_all perf regression, characterized, not fixed)
• kernels#66 (DFT-boundary colouring data race — the most safety-relevant finding of the whole effort, confirmed by two independent measurement methods)
• cpu-impl#51 (bending RLS/JRLS colour-loop under-sweep, dates to 2025-10-15 — predates this effort entirely)
• cpu-lib#66 (field_relax never exported, whole family untested)
• kernels#67 (RLS weight-neighbourhood heap-buffer-overflow, same family as reg_field/reg_flow: out-of-bounds weight-map read under Zero boundary #39)
• kernels#71 (CUDA regulariser build entirely blocked by a missing device annotation — fixed, merged as fix(bounds): mark index_stays_inbounds CUHOSTDEV #72)

Deliberately left open, not part of this umbrella's scope

• fastfields-cuda-lib dispatch wiring for phase 5's new 54-entry cuda-impl surface (only the pre-existing 6 matvec/diag entries are exported and dispatch-migrated; kernel_*/relax_*/RLS/JRLS/Lamé/all are implemented and nvcc-verified at the cuda-impl header layer but not yet wired through cuda-lib's ABI or given the Decision 2 boundary-rejection check).
• The bugs listed above under "found along the way" — each is its own tracked issue, to be picked up independently.

Thank you for the close engagement on the friction points throughout — particularly catching that the relax colouring's correctness precondition, not just its visit order, needed scrutiny (kernels#66), and the repeated insistence on independent re-measurement over trusting a prior analysis. Both materially changed what shipped.

[balbasty]

Agent: claude-fastfields-to-teeny

Phase 2 (flow) — Decision 2's flow half, measured

Decision 2 says flow's Lamé cross-coupling term needs its own measurement and
must not inherit field's set in either direction. It doesn't, and the reason is
sharper than "it might differ": the Lamé cross block is not a reach question
at all, so no value of supports_reach could have been the right answer.

Measured by two independent methods — assemble A column by column (matvec on
unit vectors) and take max|A − Aᵀ|/max|A|; and, with no matrix assembled,
|⟨Av,w⟩ − ⟨v,Aw⟩| over random v, w — on several grids, D = 1..3, isotropic
and anisotropic voxels, and against both the old hand-expanded engine and
the new one (they agree, so these are properties of the discretisation, not of
the rewrite):

bound	absolute	membrane	bending	lamé 1d / nd	lamé+bending 1d / nd
Zero	ok	ok	ok	ok / ok	ok / ok
Replicate	ok	ok	REJECT	ok / REJECT	REJECT / REJECT
DCT1	ok	REJECT	REJECT	REJECT / REJECT	REJECT / REJECT
DCT2	ok	ok	ok	ok / ok	ok / ok
DST1	ok	ok	ok	ok / REJECT	ok / REJECT
DST2	ok	ok	ok	ok / ok	ok / ok
DFT	ok	ok	ok	ok / ok	ok / ok
NoCheck	ok	ok	ok	ok / ok	ok / ok

Every ok is a measured 0. Three findings:

1. Flow's membrane and bending DO share field's table — measured, not
   assumed, and they match cell for cell.
2. The Lamé cross block's reject set is {Replicate, DCT1, DST1}, which is
   neither reach set. Replicate is exact at reach 1 on the same axis but not
   here (its clamp is idempotent, so a corner whose two taps both clamp onto the
   centre has no partner entry); DST1 is exact at reach 2 for field's bending
   — that was fix(bounds): correct the self-adjointness rejection set; add the reach-1 predicate #58's correction — but not here. Either inheritance would have
   been wrong, in opposite directions.
3. The 1d / nd split is real: no axis pair at D = 1 means no cross block
   and no extra condition, so Replicate + lamé is accepted at D = 1 and
   rejected at D ≥ 2. Hence supports_lame(b, ndim) takes ndim.

bounds.h's transpose() comment was wrong about Replicate

It argued the accepted set — "exactly self-adjoint for the half-sample family
(DCT2↔DST2) and for DFT / Zero / Replicate / NoCheck" — and Replicate
measures 0.086–0.11. Corrected to point at the measured predicate. Same shape of
error as this issue's two earlier Decision 2 corrections; arguing a fold's
involutivity has now been wrong three times for three different reasons, which
seems worth stating plainly.

The fastfields-lib#26 transpose fix is preserved, and proven load-bearing

Not merely "the numbers still match" — mutation-tested. Removing transpose()
from the cross block flips the table exactly as pre-#26:

D=2 lamé	with transpose	without
DCT2	0	0.055
DST2	0	0.037
DST1	0.027	0

Status

Phase 2 is up for review as four PRs plus a prerequisite:

• refactor(regularisers): extract the shared N-D stencil engine into regularisers/stencil.h #60 — extract the shared per-component stencil engine into
  regularisers/stencil.h (pure move; field bit-identical, 1,528,224 checks,
  deviation 0)
• refactor(reg_flow): one N-D tap-table engine for the flow regulariser #61 — the flow engine (−3,469 lines) + the measured Lamé predicates
  (2,695,888-check element-identity probe, deviation 0 on the plain paths)
• fastfields-cpu-impl#42 — kernels pin bump
• fastfields-cpu-lib#57 — the dispatch-entry rejection (test_reg_flow
  854 → 1123)
• fastfields-cpu-lib#56 — independent: fixes cpu-lib#55, without which the
  both-compiler gate cannot be run at all

One regression found and filed rather than buried: #62, 3-D matvec_all is
1.1–1.5× slower on the new engine (everything else is at parity or up to 1.6×
faster). Bounded, characterised, correctness unaffected.

Phase 3 (impl-layer peel_front().enumerate()) is unblocked once these land.

---

Generated by Claude Code