Local RBF interpolator with a linear-reproduction guarantee

[lmoresi]

Why

UW3's local scattered-point interpolator was inverse-distance (Shepard) weighting — KDTree.rbf_interpolator_local. Its weights are positive and sum to one, so it reproduces a constant exactly but not a linear field.

This is general-purpose numerical machinery and the aim is simply that it be accurate. A local interpolant that cannot reproduce a uniform gradient smears every field that has one, everywhere it is used — proxy variables, mesh-to-mesh transfer, evaluation at scattered points. No tuning downstream recovers what it discarded.

The codebase already had a linear-exact scheme — custom_mg.rbf_prolongation, polyharmonic with an affine tail — but global and dense (nnz/row == n_coarse). So UW3 had sparse-but-constant-only and linear-exact-but-dense, and nothing that was both. This adds the missing one.

Baseline, measured before writing any of it

The probe is an exactly linear field, because it lies inside both the P1 and the P2 proxy space: the discretisation contributes nothing and every bit of the measured error is transfer error.

Sweeping cell size, proxy degree, particle density and nnn (raw output in ~/+Simulations/linear-rbf-proxy/step0_shepard_baseline.txt), inverse distance gives 1e-3 to 2e-2 relative error, and:

• it falls only first order with mesh refinement and with particle density;
• it does not fall at all with stencil size — at nnn=12 it is unchanged to two significant figures, because the error is set by the asymmetry of the neighbourhood, not its size;
• boundary nodes are ~6x worse than interior ones.

So the accuracy loss is real, dominant, and not reachable by tuning.

The scheme

order= selects the polynomial reproduction order:

	weights	reproduces	bounded?
order=0 (unchanged default)	1/dᵖ, normalised	constants	yes — convex combination
order=1	polyharmonic r² log r + affine tail	constants and linears	no — signed weights

For each target, solve the small saddle system on its own kNN stencil in a local scaled frame (the kernel is not scale-invariant, so this is required for conditioning, not cosmetic):

[ A   P ] [ w ]   [ φ(|x* - x_j|) ]
[ Pᵀ  0 ] [ λ ] = [ 1,  x*        ]

The lower block is the reproduction constraint Pᵀw = [1, x*]. Cost is one (nnn+d+1)³ dense solve per point, batched and chunked; the result is sparse at nnn non-zeros per row.

Measured

Swarm proxy against the analytic field, at the proxy's own nodes:

	2D order=0	2D order=1	3D order=0	3D order=1
linear, max	4.6e-3	5.1e-16	4.8e-3	5.5e-16
quadratic, max	9.2e-3	1.0e-4	8.5e-3	2.4e-4
quadratic, rms	2.1e-3	2.1e-5	2.0e-3	2.3e-5

(cellSize=1/8, fill_param=4, proxy degree 1; nnn = dim+1 vs 2*(dim+1).)

The linear column is the clean statement of the guarantee. The quadratic column is what makes it a general accuracy improvement rather than a special case: ~90x better in 2D, ~36x in 3D on max error, ~100x on rms. Full sweep in ~/+Simulations/linear-rbf-proxy/step1_shepard_vs_linear.txt.

Other validation, all measured:

• partition of unity |Σw − 1| ≤ 4.4e-16
• linear reproduction |Σ w_j x_j − x*| ≤ 5e-15
• sparsity exactly nnn non-zeros per row
• convergence: on a quadratic field, halving the point spacing cuts order=1 error by >3x while order=0 barely moves
• order=0 positional calls are bit-identical (np.array_equal) to before

What changes behaviour

Only the swarm proxy path: SwarmVariable.rbf_interpolate / _rbf_to_meshVar default to order=1, nnn=2*(dim+1).

Deliberately left on order=0:

• MeshVariable.rbf_interpolate — it is the RBF rung of the uw.function.evaluate point-location ladder, whose contract in design/point-location-capability.md is "bounded, topology-free, honest". Signed weights would break a documented promise. Also covers adaptivity.py mesh2mesh.
• read_timestep (both) — nnn=1 is an exact checkpoint round-trip.
• IndexSwarmVariable — level-set fractions must stay in [0,1]; it hand-rolls its own IDW and is untouched.

Traps handled

Degenerate stencils. Pre-screen the affine block by SVD, solve only the healthy ones, then validate the answer (both reproduction identities + finiteness) rather than guessing a condition number. Failures fall back to inverse distance with one counted warning — #424's lesson is that a silent geometric fallback surfaces months later. Tested with collinear 2D, coplanar 3D, collapsed stencil, duplicate sources, and a mixed batch where one singular member must not take the batched solve down with it. Never NaN.

nnn guard. order=1 requires nnn >= dim+2 and the error says why: at dim+1 the affine block is square, the RBF block is inert, and the scheme is bare barycentric interpolation on the neighbour simplex — singular whenever those points are collinear/coplanar. Today's proxy default of dim+1 was precisely that fragile case.

Limiting — and a defect fixed in review (d395424). The first version of monotone clipped the result to the min/max of its stencil's source values. That is wrong for a linear-exact scheme, and wrong subtly: a target outside the convex hull of its own nnn neighbours has a value outside their range even for an exactly linear field, and proxy nodes are routinely in that position — not only at domain boundaries. So the clip could not distinguish legitimate extrapolation from ringing and fired on the linear part, returning a linear field to ~4e-3 (inverse-distance level).

Replaced with the slope-limiter discipline — never limit the linear reconstruction, limit only the correction on top of it: fit an affine trend to the stencil by least squares, keep it at the target untouched, and bound the remaining RBF correction to the non-affine residual the stencil actually shows. Measured, the limiter is now a no-op on any field the scheme reproduces exactly (moves the answer by ~1e-15). It still bites on rough data — on sin(6x)+|x|² it moved the answer by 2e-3 (2D) / 5e-2 (3D), and in 3D that costs accuracy (9.0e-2 → 1.3e-1, still well under inverse distance at 2.7e-1), the ordinary limiter trade.

It bounds new oscillation relative to the local trend, not absolute range — a quantity with hard physical bounds still needs its own clip. Off by default.

Empty ranks. A rank with fewer particles than the affine tail needs drops to order=0 automatically. The collective structure of the refresh is unchanged: every rank still performs the same single read-then-write per proxy (SWARM-07 / #313), and the starved-rank hold-current policy is untouched.

Cost

~5x per proxy refresh, and refresh runs at every solve entry:

	2D, 24k particles / 4.9k nodes	3D, 160k particles / 12.7k nodes
order=0, nnn=dim+1 (old default)	4.5 ms	21.7 ms
order=0, nnn=2(dim+1)	5.3 ms	27.2 ms
order=1, nnn=2(dim+1) (new default)	23.9 ms	111.5 ms

Widening the stencil is ~20% of that; the dense per-point solves are the rest. Small next to a 3D Stokes solve, but not free for proxy-heavy models with cheap solves. Obvious untaken lever: cache the weights while particle positions are unchanged (they are geometry-only, and Swarm._get_kdtree() is already invalidated on exactly the right events).

Parallel

Proxy refresh is rank-local with no halo exchange (SWARM-15). A linear-exact stencil reproduces a linear field exactly from any neighbourhood, one-sided or not — so for linear fields the seam error is zero and the proxy becomes np-independent. This reduces np-dependence; it does not remove it for fields with curvature. The halo exchange in SWARM-15 is still the real fix, and the parallel test says so rather than claiming more.

tests/parallel/test_0776_... passes at np=2 and np=4, numbered in the 07xx band so scripts/test.sh and test_levels.sh actually collect it.

Regression

	passed	failed
feature/linear-rbf, level_1 tier_a+b	494	0
feature/linear-rbf, level_2 tier_a+b	312	0
development baseline, level_1 tier_a+b	464	2

The counts reconcile: 464 + 28 new serial tests + the 2 that fail on the baseline = 494. Those 2 are test_0103_jit_rampable_constants (from #416); they pass in the freshly-built worktree, so it is a stale JIT cache in the main-repo checkout, unrelated to this branch. No test changed status.

Lagrangian / semi-Lagrangian specifically, since the proxy is on the DDt critical path: test_0857, test_1052, test_1053, test_1054, test_1100, test_1120 — 30 passed.

Not claimed

Whether this helps the free-surface topography noise recorded in the planning file is an open question, not the motivation. No free-surface case was run here; everything measured is interpolation accuracy plus the regression suite. The Crameri Case-2 comparison would settle it and that item stays open.

Defects found in passing — filed, not fixed (Charter §9)

• Swarm proxy refresh is broken whenever reference quantities are set (dimensional coords into a non-dimensional kd-tree) #426 — swarm proxy refresh is broken whenever reference quantities are set: dimensional meshVar.coords into a non-dimensional kd-tree. Three sites. Also _create_proxy_variable drops units=.
• rbf_interpolator_local: inverse-distance weights use SQUARED distances, so the documented p semantics are wrong (effective exponent is 2p) #427 — inverse-distance weights use squared distances, so the decay is r^(-2p) not the documented r^(-p); epsilon floors r at ~1e-6.
• Dead parameters and dead code in the RBF interpolation path (and a test silently discarding its nnn) #428 — dead meth/rubbish parameters, dead _rbf_reduce_to_meshVar/_nn_proxy, and test_0505 silently passing its nnn into meth.
• custom_mg.rbf_prolongation: no conditioning path on an indefinite saddle solve, and its only test asserts convergence rather than accuracy #429 — custom_mg.rbf_prolongation has no LinAlgError/conditioning path and its only test asserts convergence, not accuracy.

# TODO(BUG): markers naming each issue are at the exact locations.

Out of scope

custom_mg.rbf_prolongation is not rewritten, and #425's nested topological transfer is untouched. An assessment for #425 is posted on that issue.

Docs

New docs/developer/subsystems/interpolation.md — there was no interpolation or kd-tree subsystem doc. Registered in the explicit toctree and the authority table in docs/developer/index.md, cross-linked from data-access.md.

Underworld development team with AI support from Claude Code

[lmoresi]

Adversarial review

Reviewing my own PR against the strongest objections I can make to it. Two of these are things I had not measured when I opened it; both are now measured and both are genuine costs.

1. It is 5x more expensive per proxy refresh, and I had not measured that

Proxy refresh runs at every solve entry (Mesh.update_lvec -> Swarm._sync_before_assembly) and inside every DDt step, so a constant-factor slowdown there is not free. Measured, simplex box, proxy_degree=2, fill_param=3:

	2D, 24k particles / 4.9k nodes	3D, 160k particles / 12.7k nodes
order=0, nnn=dim+1 (old default)	4.5 ms	21.7 ms
order=0, nnn=2(dim+1)	5.3 ms	27.2 ms
order=1, nnn=2(dim+1) (new default)	23.9 ms	111.5 ms

Decomposed: widening the stencil costs ~20%, and the dense per-point solves cost the remaining ~4.5x. So it is the solves, not the kNN query.

In absolute terms 112 ms is small next to a 3D Stokes solve, and the accuracy buys ~2 orders of magnitude. But I should not have opened the PR claiming a default change without this number, and anyone running a proxy-heavy model with cheap solves (pure advection, tracer-only) will feel it. If that turns out to matter, the obvious lever is caching the weights while the particle positions are unchanged — the weights are geometry-only, and Swarm._get_kdtree() is already invalidated on exactly the events that would invalidate them. Not done here; deliberately out of scope.

2. The motivating claim is not actually demonstrated

The external planning item this is aimed at attributes residual free-surface topography noise to the P1-proxy/RBF particle->mesh mapping. I did not run a free-surface case. Everything measured here is static interpolation accuracy on analytic fields plus the existing regression suite.

So the honest status is: this removes a demonstrated, dominant source of transfer error that is a plausible cause of that noise. Whether topography noise actually drops to the mesh-level-set level is unverified, and the Crameri Case-2 comparison is the test that would settle it. Nobody should read this PR as closing that item.

3. monotone=True is much weaker than it sounds, and I nearly shipped it as the safe option

Clamping to the stencil range and linear exactness are in direct conflict: wherever the target lies outside its neighbours' convex hull, the exact linear value is outside the stencil range, so the clamp throws it away. Measured on the 2D case, order=1 + monotone on a linear field gives 4.0e-3 max error — straight back to inverse-distance level — while interior nodes stay exact.

This is called out in the docstring, the PR body and the subsystem doc, because the intuitive reading ("bounded and accurate") is wrong. It is False by default for that reason.

The uncomfortable corollary: order=1 really does give up the maximum principle, and the swarm proxy feeds Lagrangian/Lagrangian_Swarm psi_star history. The regression suite (30 Lagrangian/SL tests) is clean, but the suite is not a proof that no model overshoots somewhere it matters.

4. Things I checked that could have been wrong, and were not

• Did the evaluate fallback rung silently change? No. _function.pyx:1371 and :1523 both resolve to MeshVariable, whose default stays order=0, so the "bounded, topology-free" contract in design/point-location-capability.md is intact. Same for adaptivity.py:903.
• Is order=0 really bit-identical? Asserted with np.array_equal, not allclose.
• Does a singular member kill a batched solve? Yes, np.linalg.solve raises for the whole batch and does not say which row — which is why the rank pre-screen exists, with a pinv retry behind it and a reproduction check in front of the result. Tested with a deliberately mixed batch.
• Does the degeneracy branch break collective symmetry? No. It lives entirely inside the numeric kernel; _rbf_to_meshVar still performs exactly one read-then-write per rank per proxy (SWARM-07 / fix(swarm): Track-0 stale-cache and migration-semantics fixes (BF-02/05/06/07/08 + #286) #313).
• Memory. First version budgeted only the stacked saddle matrices and missed the (nnn, nnn, dim) offset temporary, which is the larger one in 3D — real peak was ~3x the stated budget. Fixed in bedfef2.
• Coupling. _normalise_monotone imported from function/ on every call, including the default path, putting a low-level module behind the evaluator's import graph for no benefit. Now returns before the import. Also bedfef2.

5. Where I would push back on a reviewer

If asked to make order=1 the default for MeshVariable.rbf_interpolate too, for consistency: no. The asymmetry is deliberate and documented. That method is a fallback for points the FE machinery could not locate, and its value is that it cannot invent an out-of-range answer. Consistency of defaults is worth less than keeping a stated contract.

6. Weakest remaining part

test_0102's convergence test asserts order=1 error falls by >3x under halving the point spacing, which is a loose proxy for "at least O(h^2)". It is a threshold chosen to be robust on a small grid rather than a measured rate, and it would not catch a scheme that had silently degraded to O(h^1.8). The exactness and partition-of-unity tests are tight (1e-12) and carry the real weight; the convergence one is a smoke test and should be read as such.

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[lmoresi]

Correction to the review above — the limiter was a defect, not a trade-off

Maintainer review made a point that turns out to be decisive: if it hurts a linear field, it is clearly a problem — how can that overshoot through internal curvature?

That is exactly right, and it identifies a design error rather than a caveat. My §3 above described the behaviour and then rationalised it. It should have been treated as a bug.

The mechanism. The clamp took its bound from the stencil's own source values. But a target that lies outside the convex hull of its nnn nearest neighbours has a value outside their range even for an exactly linear field — that is what extrapolation is. And proxy nodes are routinely in that position, not only at domain boundaries: the nnn nearest particles to an interior node can easily all sit on one side of it. So the clamp had no way to distinguish legitimate extrapolation from ringing, and clipped the linear part. There is no "internal curvature" involved at all; the overshoot it was reacting to was not overshoot.

The fix (d395424) is the standard slope-limiter discipline: never limit the linear reconstruction, limit only the correction on top of it.

1. least-squares affine fit over the stencil (rbf_stencil.affine_trend)
2. keep that trend at the target untouched
3. bound the remaining RBF correction to the non-affine residual the stencil actually exhibits

For an exactly linear field every residual is zero, so the bound is [0, 0], the correction is forced to zero, and the result is the exact affine value. Linear reproduction survives the limiter by construction, not by tolerance.

Measured:

field	limiter off	limiter on	limiter moved the answer by
linear, 2D	1.3e-15	8.9e-16	1.3e-15 (was 4e-3)
linear, 3D	3.6e-15	2.7e-15	1.8e-15
sin(6x)+abs(x)^2, 2D	5.3e-3	5.5e-3	2.0e-3
sin(6x)+abs(x)^2, 3D	9.0e-2	1.3e-1	5.1e-2

So it is now a no-op on anything the scheme reproduces exactly, and still limits rough data. The 3D row is the ordinary limiter trade — accuracy given up for boundedness where the correction genuinely exceeds the observed residual — and it is still well below inverse distance (2.7e-1) on that field.

What it promises is now stated precisely, because the previous wording overpromised: it bounds new oscillation relative to the local trend, not absolute range. A quantity with hard physical bounds (a fraction in [0,1]) still needs its own clip on top. Two tests that asserted the old, wrong guarantee have been replaced by ones asserting the new one, including that the limiter leaves a linear field alone.

Also corrected: the framing

The PR body has been rewritten. Motivation is general interpolation accuracy — a local interpolant that cannot reproduce a uniform gradient smears every field that has one, everywhere it is used. I had over-anchored on the free-surface topography item from the planning file; that is now stated as an open question this does not answer, not as the reason for the work. Nothing about the measurements changed, only what they are claimed to be for.

Re-ran level_1 tier_a+b after the limiter change: 494 passed, 0 failed.

[lmoresi]

This is a comment from the point of view of a user of this tool (for mesh -> mesh transfers)

1. The raw-array path has a footgun the value path does not

KDTree.rbf_interpolator_local(..., order=1) falls back to IDW on a degenerate stencil, counted and warned. linear_exact_weights returns zeros plus a mask and leaves substitution to the caller.

That asymmetry is fine in itself, but the raw path is the one anyone building an operator must use — you cannot get a reusable matrix out of the value API. So the path with the sharp edge is the path serious consumers are pushed onto, and a zero row is not a small inaccuracy: for a transfer operator it silently breaks partition of unity.

Either give the matrix path the same automatic fallback, or make the docstring say plainly: "returns zero rows for degenerate stencils; the caller MUST substitute — see degenerate". One sentence.

2. Test with a RANDOM field, not a linear one

Your measured |sum w_j x_j - x*| <= 5e-15 is a linear-reproduction check. Necessary, but it cannot detect wrong-neighbour selection — any local averaging of nearby values reproduces a linear function exactly.

I watched an equivalent test pass in both directions this session, on code that was and was not correct. A random nodal field with an independently computed reference is the discriminating test. Worth adding to test_0102 while you are in there.

3. Warning !! Do not use uw.function.evaluate as your deciding call on accuracy.

There are still cases where this returns wrong values at points lying exactly on cell edges/faces (#432 — a recurrence). Swarm proxy points and mesh-transfer targets land on cell boundaries routinely, so this is live for your validation, not hypothetical.

4. Three things only you can state, and downstream cannot discover

• Tie-breaking determinism. If two sources are equidistant from a target, is KDTree.query deterministic, and identical regardless of point ordering or partition? Adapted meshes in this repo are held to bit-confluence; a transfer built on a non-deterministic tie would break that.

• Purity. Confirm the function is a pure function of its arrays with no hidden global state — that is what makes it valid to call on an all-gathered cloud in parallel.

• No column guarantee. Row-wise kNN gives nnn non-zeros per row and says nothing about columns. Not yours to fix, but one sentence in the docs prevents a downstream consumer assuming it. (For MG a zero column makes PtAP singular — that is 3D: shape relaxation can break the custom-P FMG build (zero-column transfer) — silent GAMG fallback #424, and this does not fix it.)

5. Make the DEFAULT nnn safe, not just the minimum

You enforce nnn >= dim+2. But at exactly dim+1 the affine block becomes square and the RBF block inert, and dim+2 is only one clear of that. The recommendation in the notes from this PR is 2*(dim+1).

Enforcing a minimum is not the same as defaulting well. If the default is inherited from the IDW path (nnn=4), then in 3D it is below your own enforced minimum and will raise. Please make the default dimension-dependent and set it to the recommended value, so the safe thing happens without the caller knowing any of this.

6. Two documentation facts that will bite silently

• .coords_nd, never .coords — building a tree from dimensional coordinates under an active units model fails (Swarm proxy refresh is broken whenever reference quantities are set (dimensional coords into a non-dimensional kd-tree) #426). Invisible in a unit-less test.
• Output carries no units, and cost is ~5x IDW — fine for a one-time operator build, think twice in a hot loop. Both worth stating at the point of use.

Everything else in my previous comment is for whoever wires this into multigrid, and is recorded on #425 instead.

[lmoresi]

Deployability pass — three gaps closed, one split off

Making this safe to sweep into development for swarms. Each item was measured, not assumed.

1. Blocker: proxied swarm variables did not work under units at all — fixed in #434

Verified reproducing before any change:

proxy refresh under units: FAILED -> ValueError: KD-tree was built with
dimensionless coordinates, but query coordinates have units 'meter'.

The swarm kd-tree is built from non-dimensional particle coordinates; the refresh queried it with MeshVariable.coords, which dimensionalises once reference quantities are set. Three live sites, IndexSwarmVariable included.

This is a pre-existing bug, not caused by this branch, so it went to its own branch as #434 (bugfix/swarm-proxy-coords-nd) against development — cherry-pickable, and it does not need this feature to land. Regression test written first and shown to fail. #430 should merge after #434.

The other half of #426 — _create_proxy_variable not forwarding units= — is deliberately left open: attaching units there makes var.sym unit-aware and changes symbolic behaviour everywhere proxies are composed, so it needs its own design.

2. Isolated spikes at degenerate points — fixed here (7b6e398)

Running the new default at scale exposed something the accuracy sweeps could not. 3D, 190k particles, 29k proxy nodes, fill_param=2: 2 stencils in 28824 failed the affine-rank test, and those two nodes carried the entire max error — 2.2e-3 while every other node sat at 1e-15.

That shape of error is worse than its size suggests: a single non-exact node among exact ones is a spike, and spikes damage derivatives and integrals far more than a smooth error of equal magnitude.

Degenerate stencils are now retried on a wider neighbourhood (4× nnn, capped at the cloud size) before falling back. Only the failed points are re-queried.

	before	after
3D cell=1/16 fill=2, np=1	2.213e-03	8.882e-16
3D cell=1/16 fill=2, np=4	2.452e-03	1.599e-15

A cloud that is coplanar over a wide region genuinely cannot be rescued by widening, and there the behaviour is unchanged and correct: bounded inverse-distance values plus a loud warning. Both paths now have tests — one asserting recovery, one asserting the honest fallback fires.

3. IndexSwarmVariable (materials) — measured, deliberately NOT changed

It never used rbf_interpolate at all; it hand-rolls inverse distance (nnn=5, radius_s=0.5). Verified: _meshVar is None, and neither rbf_interpolate nor rbf_interpolator_local appears in _update_proxy_variables.

The structural reason not to port it: a material indicator is piecewise constant. Away from an interface both schemes are exact, because both reproduce constants — linear exactness adds nothing. At the interface the field is discontinuous, so no polynomial-reproducing scheme is exact either; signed weights only add overshoot where the data jumps.

Measured on a straight interface at x = 0.5 (exactly representable):

scheme	interface error, median	level-set range
inverse distance, nnn=5	5.2e-3 – 1.1e-2	[0, 1] exactly
order=1, nnn=6	6.9e-3 – 7.8e-3	[0, 1] exactly
order=1, nnn=8	5.2e-3 – 7.9e-3	[-0.038, 1.038]

A wash on accuracy — better at coarse resolution, equal or worse at fine, with the ordering flipping — and nnn=8 breaks the [0,1] bound by 3.8%. Partition of unity survives either way (all indices share one weight set), but a negative material fraction is still wrong and constitutive_models.py consumes these directly.

So the swarm story is deliberately split: smooth fields take order=1 for a two-orders-of-magnitude gain; the material jump keeps inverse distance because there is nothing to win and a bound to lose. Documented with the numbers in subsystems/interpolation.md, including the caveat that the interface metric's maximum is set by node spacing rather than the weights, so only the median is informative.

4. Scale and parallel

Up to 190k particles / 31k nodes, 2D and 3D, np=1 and np=4:

• linear field exact everywhere — worst case 1.6e-15 across all cases
• np=1 and np=4 agree to round-off
• degenerate fallback after retry: zero warnings on every case tested

Status

I consider the interpolator deployable for swarms once #434 lands first. Remaining known limits, all documented rather than hidden: ~5× refresh cost with no weight caching, the SWARM-15 rank-local seam for fields with curvature, and the monotone default left off pending a model that actually rings.

Regression after all of the above: level_1 496 passed, level_2 312 passed, np=2 parallel 6 passed, 0 failures.

Artifacts: ~/+Simulations/linear-rbf-proxy/ (step2_index_swarm_levelset.txt, step3_scale_and_parallel.py).

[lmoresi]

Acting on the multigrid-transfer session's review

Six points, from the first non-swarm consumer. All six were right; two were live bugs rather than documentation gaps. Landed in 44ab4fb.

1. The raw-array path had a footgun the value path did not — fixed properly

Correct, and worse than stated. linear_exact_weights returns zeros plus a mask for degenerate rows, and a zero row silently interpolates to zero — not merely inaccurately. Since the raw kernel was the only route to an operator, the sharp edge faced exactly the serious consumers. My own documented recipe for building T did not handle the mask.

Rather than document the asymmetry, I removed it: KDTree.interpolation_matrix(coords, nnn=None, p=2, order=0) returns the sparse operator with the value path's full handling — widen the stencil, then inverse distance with a warning. Rows are never empty.

Both paths are now built on one _local_stencil, so the operator and the values it produces cannot drift, and a test asserts T @ data == rbf_interpolator_local(...) for both orders and both dimensions. The raw kernel keeps a Warnings block pointing at the operator API.

2. Test with a random field, not a linear one — correct

The reasoning is exactly right: |Σ w_j x_j − x*| is evaluated against whatever points were handed to the solver, so it holds even if the tree returned the wrong neighbours — and a linear field then interpolates exactly anyway. My whole validation was blind to neighbour selection.

Added: a brute-force kNN + dense-solve oracle that shares no code with the library, run on random values in 2D and 3D. Plus the control that makes it meaningful — asserting that deliberately shuffled stencils do change the answer. Without that the test proves nothing.

3. uw.function.evaluate as oracle (#432) — checked, not affected, now documented

Verified: none of my tests use evaluate; they compare against analytic fields directly. But the warning is worth carrying, because proxy nodes sit on cell vertices and edges routinely — so it is now an explicit "do not validate this way" note in the subsystem doc, pointing at #432.

4. Three things only I can state — now stated and tested

• Determinism/purity: the weights are a pure function of geometry — no global state, no RNG, no accumulation. Repeated calls, and separate KDTree instances over equal point sets, are bit-identical.
• Tie-breaking: with several exactly-equidistant sources, which is selected is unspecified, but the choice is stable across calls and instances. Tested with a symmetric ring where every candidate ties.
• No column guarantee: every row has nnn non-zeros; nothing guarantees every column is non-empty, so a Galerkin PᵀAP can still be singular (3D: shape relaxation can break the custom-P FMG build (zero-column transfer) — silent GAMG fallback #424). Stated in the interpolation_matrix docstring and the docs, flagged as a property of any row-wise kNN construction and not something linear-exactness fixes.

5. Dimension-dependent default nnn — this was a live bug

Confirmed before changing anything:

3D default nnn + order=1: RAISES -> A linear-exact stencil needs at least
dim + 2 = 5 neighbours, got nnn = 4.

rbf_interpolator_local(coords, data, order=1) raised on its own default in 3D. Enforcing a minimum is not defaulting well, as you put it. nnn=None now resolves to 2*(dim+1) at order=1 and keeps 4 at order=0, so back-compatibility is intact and both are exercised in 2D and 3D by test.

6. Silent-failure docs facts — added

A dedicated section: .coords_nd never .coords (#426), no units on output, ~5× the order=0 cost, and the evaluate warning from (3).

One I broke and you would have caught

The refactor swapped query(k=1) for find_closest_n_points, which added a stencil axis to the nnn=1 checkpoint round-trip path. Only a snapshot test noticed — could not broadcast (30,1,1) into (30,1). Fixed, and the shape contract now has its own test instead of relying on a distant consumer.

Regression

level_1 511 passed, level_2 312 passed, np=2 parallel 4 passed, np=4 scale sweep unchanged at 1.6e-15 worst case. 0 failures.

On the #425 items

Agreed and understood as not-a-demand: column repair, sweeping nnn, the iteration-count bar against GAMG, and topological-over-geometric where a parent map exists all belong to whoever wires this into multigrid. The one thing I have moved toward them is that interpolation_matrix now exists, so that work starts from an operator with sane degenerate handling rather than from the raw kernel.

[lmoresi]

Adversarial review of the post-review commits (retrospective)

The adversarial review on this PR covered the state at d3954245. Three substantial commits landed after it — the deployability pass (7b6e3981), the operator refactor (44ab4fbb) and the merge from development — and none of them got one. Running it now against the merged state, 4d043aeb.

One real defect, filed as #443

interpolation_matrix(nnn=1, order=1) silently returns a nearest-neighbour operator where the value path raises:

kdt.rbf_interpolator_local(tgt, f, 1, 2, False, order=1)   # ValueError
kdt.interpolation_matrix(tgt, nnn=1, order=1)              # no error, nnz/row = 1

The caller asked for a linear-exact operator and got a constants-only one, with nothing recording the downgrade. This is exactly the failure mode the rest of the PR was written to avoid.

It also falsifies a claim I made in the PR body: that the two paths share one core and "cannot drift apart". They share _local_stencil, but _local_stencil early-returns for nnn == 1 before the order handling, and the guard lives in rbf_interpolator_local_from_kdtree only. Sharing a core is not the same as sharing the guards, and I asserted the stronger thing.

The consistency test parametrises order but not nnn, and never exercises nnn=1, so it could not have caught this.

The refactor was riskier than the diff suggests

44ab4fbb restructured a working, tested path to extract _local_stencil. It broke the nnn=1 checkpoint round-trip by swapping query(k=1) (which reshapes to (n,)) for find_closest_n_points (which does not) — caught only by a snapshot test, two subsystems away, as could not broadcast (30,1,1) into (30,1).

That is a near miss, not a clean catch. Had the snapshot format tests not existed, a silent shape change would have reached read_timestep. The shape contract now has its own test, but the general lesson stands: the refactor was justified on DRY grounds and I ran it on a path whose only coverage was incidental.

Claims that are weaker than they read

• "~5× cost" is one hardware, one mesh family, proxy_degree=2, fill_param=3. It is not a scaling study, and the retry path added in 7b6e3981 was not re-benchmarked — the cost of widening is bounded by how many stencils fail, which was 2 in 28824 on the case measured and could be much worse on a strongly graded mesh.
• "np=1 and np=4 agree to round-off" is for a linear field, which is the case where rank-locality provably does not matter. For a curved field the seam error is unmeasured — SWARM-15's own migration plan asks for exactly that measurement and it still has not been made.
• The IndexSwarmVariable decision rests on an interface metric that groups unstructured-simplex nodes into rows by y. Its maximum is identical across all schemes because it is set by node spacing, which I noted — but that also means the metric is crude, and the conclusion leans on the median plus the [0,1] violation. The bound violation is solid; the accuracy wash is softer than the table makes it look.

Smaller things

• interpolation_matrix with nnn > self.n raises RuntimeError: Error in rbf_interpolator_local_from_kdtree ..., naming a function the caller never called.
• The degenerate-fallback warning uses stacklevel=3, tuned for the value path's depth; from interpolation_matrix it is off by a frame.
• interpolation_matrix returns a scipy CSR. Charter §11 is explicit that serial libraries are usually a poor choice; it is the right currency for custom_mg (which is serial-only gated) but the method carries no note that it is rank-local, which an MPI consumer would want stated.
• The monotone path re-gathers np.asarray(self.points)[closest_n], a second large gather of data already materialised upstream.

What holds up

• The measurements are reproducible from the scripts in ~/+Simulations/linear-rbf-proxy/, and the negative results (limiter cost in 3D, IndexSwarmVariable wash, degenerate fallback fraction) are reported rather than buried.
• order=0 bit-identity is asserted with array_equal, not allclose, and survived the refactor.
• The random-data oracle with a shuffled-stencil control is the one test in the file that can actually detect wrong-neighbour selection; adding it was the single most valuable review outcome.