Skip to content

Inverse-distance weights use the actual distance; p defaults to 1 (#427)#448

Open
lmoresi wants to merge 1 commit into
developmentfrom
bugfix/idw-actual-distance
Open

Inverse-distance weights use the actual distance; p defaults to 1 (#427)#448
lmoresi wants to merge 1 commit into
developmentfrom
bugfix/idw-actual-distance

Conversation

@lmoresi

@lmoresi lmoresi commented Jul 27, 2026

Copy link
Copy Markdown
Member

Fixes #427.

The defect

KDTree.query / find_closest_n_points return squared distances, and the weighting used them directly:

weights = 1 / np.power(epsilon + distance_n, p)     # distance_n is r²

So the decay was r^(-2p), not the r^(-p) the argument is named and documented for. At the old default p=2 the scheme was inverse fourth power — close to nearest-neighbour — not the "inverse distance (squared)" the docstring claimed.

The retained old_* implementations took sqrt explicitly, so this changed silently at a rewrite. Nothing caught it: no test pinned the exponent, and the only accuracy test used rtol=5e-2, loose enough that 1/r² and 1/r⁴ both pass.

The change (maintainer decision)

  1. Take the square root, so p applies to the actual distance.
  2. Default p 2 → 1, so the default is genuinely inverse distance.

epsilon stays 1e-12 but now floors r rather than , so its scale is the one it reads as — it was an effective 1e-6 length floor.

Verification

Two sources at distance 1 and 2 from the target; interpolating a field that is 1 at the near point and 0 at the far one recovers the weight ratio directly:

p measured w₁/w₂ expected 2^p
1 2.0000 2
2 4.0000 4
3 8.0000 8

The exponent had no test at all before; it does now, including one pinning the default at p=1.

How much it actually changes

Linear field through the swarm-proxy stencil (nnn=dim+1), max relative error:

p=1 p=2
2D 8.79e-2 8.77e-2
3D 2.05e-1 2.04e-1

The choice of p barely matters. What mattered was that p meant something other than its name — and that since #430, order=0 is the fallback rather than the primary path, so this is a correctness-of-contract fix rather than an accuracy one.

Call sites

Stopped hard-coding p=2: SwarmVariable.rbf_interpolate lets the default apply, MeshVariable.rbf_interpolate's own p default follows to 1, and the nnn=1 checkpoint call drops a positional p it never used (nnn=1 short-circuits before any weighting).

Test gap closed

test_interpolation_matrix_agrees_with_the_value_path now sweeps nnn and p as well as order. It previously pinned order only and hard-coded p=2 on the value side — so it compared two different weightings the moment the default changed, and duly failed here. That gap was flagged in this series' own adversarial review; this is it biting.

Regression

  • level_1 tier_a/tier_b: 567 passed, 0 failed
  • level_2 tier_a/tier_b: 319 passed, 0 failed

Underworld development team with AI support from Claude Code

The kd-tree returns SQUARED distances, and the weighting used them
directly:

    weights = 1 / (epsilon + distance_n)**p     # distance_n is r^2

so the decay was r^(-2p), not the r^(-p) the argument is named and
documented for. At the old default p=2 the scheme was inverse FOURTH
power -- close to nearest-neighbour, not the "inverse distance (squared)"
the docstring claimed. The retained old_* implementations took sqrt
explicitly, so the semantics changed silently at a rewrite and nothing
caught it: no test pinned the exponent, and the only accuracy test used
rtol=5e-2, loose enough that 1/r^2 and 1/r^4 both pass.

Two changes, maintainer decision:

  - take the square root, so p applies to the actual distance;
  - default p 2 -> 1, so the default is genuinely inverse distance.

epsilon stays 1e-12 but now floors r rather than r^2, so its scale is the
one it reads as (it was an effective 1e-6 length floor).

Pinned by test: two sources at distance 1 and 2 give a weight ratio of
exactly 2**p for p in (1, 2, 3), and the default gives 2. The exponent
had no test at all before.

Measured effect on a linear field through the swarm-proxy stencil
(nnn=dim+1): 2D 8.79e-2 at p=1 vs 8.77e-2 at p=2; 3D 2.05e-1 vs 2.04e-1.
The choice of p barely matters -- what mattered was that p meant
something other than its name, and that order=0 is now the fallback
rather than the primary path.

Call sites updated to stop hard-coding p=2: SwarmVariable.rbf_interpolate
now lets the default apply, MeshVariable.rbf_interpolate's own p default
follows to 1, and the nnn=1 checkpoint call drops a positional p it never
used.

test_interpolation_matrix_agrees_with_the_value_path now sweeps nnn and p
as well as order. It previously pinned order only and hard-coded p=2 on
the value side, so it compared two different weightings the moment the
default changed -- a gap flagged in this PR series' own review.

level_1 567 passed, level_2 319 passed.

Underworld development team with AI support from Claude Code
Copilot AI review requested due to automatic review settings July 27, 2026 12:21

Copilot AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

@lmoresi

lmoresi commented Jul 27, 2026

Copy link
Copy Markdown
Member Author

Adversarial review

This changes numbers for every existing order=0 caller, and I have not measured most of them

The weighting goes from an effective 1/r⁴ to 1/r. That is a large change in stencil sharpness — the old scheme was close to nearest-neighbour, the new one is markedly smoother. I measured it on one thing: a linear field through the proxy stencil, where it is a wash (8.79e-2 vs 8.77e-2).

Consumers still on order=0 that I did not measure:

  • MeshVariable.rbf_interpolate — the RBF rung of the uw.function.evaluate point-location fallback ladder. Points that fail location now get a smoother extrapolant. Still bounded (positive weights, convex combination), so no contract is broken, but the values move.
  • adaptivity.py:903 mesh2mesh transfer.
  • The order=1 degenerate fallback, which uses these weights at the handful of points where the affine fit fails.

The suites pass, but "no test asserts a number that moved" is not the same as "nothing moved". A caller relying on the old near-nearest-neighbour sharpness will see different answers with no deprecation signal.

The measurement I chose flatters the change

A linear field is close to the worst possible probe for a weighting-exponent change: inverse distance smears a gradient regardless of p, so the error is dominated by the scheme's inability to reproduce linears, not by the exponent. That the two agree to three significant figures says little.

A field with curvature, or a discontinuity, would separate them properly — 1/r⁴ tracks a jump much more sharply than 1/r. I did not run one. The honest statement is "p barely matters for the linear-field metric I used", not "p barely matters".

epsilon is now arguably in the wrong place

epsilon = 1e-12 on the actual distance is a 1e-12 length floor. For non-dimensional coordinates on an O(1) domain that is fine. For a kd-tree built in metres over a planetary domain it is still fine; for one built in units where the domain is O(1e-12) it is not. It was previously an effective 1e-6 floor, which was wrong differently. Neither version scales epsilon to the data, and nothing tests coincident points.

The default change is a second decision riding on a bugfix

Taking the square root is unambiguously a bug fix. Moving p from 2 to 1 is a taste decision about default smoothing, and it is bundled into the same commit. They are separable: the fix alone (p=2 on actual distance → 1/r²) would have been the conservative landing. Bundled, a reviewer cannot approve one without the other, and a future bisect cannot separate them.

Maintainer chose p=1 explicitly, so this is recorded rather than contested — but it should be visible that the PR does two things.

What holds up

  • The exponent now has a direct test, on a two-point geometry where the expected ratio is exact arithmetic rather than a tolerance. That test would have caught the original regression the moment it was introduced.
  • The consistency test's p/nnn sweep — flagged as missing in this series' own earlier review — failed on this change and was fixed rather than loosened. That is the gap closing as intended.
  • The call sites stopped hard-coding p, so the default is now expressed in exactly one place.

@lmoresi

lmoresi commented Jul 27, 2026

Copy link
Copy Markdown
Member Author

Closing the gap I flagged: measured on curved and discontinuous fields

The review above said the linear-field metric was the worst possible probe for an exponent change and that I had not run a better one. Running it now. p=4 on the actual distance reproduces the old effective 1/r⁴, so the last column is the previous behaviour.

Max error, nnn=dim+1, 4000 sources / 400 targets:

p=1 (new default) p=2 p=4 (old behaviour)
2D linear 2.79e-2 2.61e-2 2.57e-2
2D curved (sin 6x + |x|²) 5.73e-2 5.86e-2 6.11e-2
2D jump 1.00 1.00 1.00
3D linear 1.29e-1 1.30e-1 1.33e-1
3D curved 1.69e-1 1.66e-1 1.77e-1
3D jump 8.30e-1 8.90e-1 9.60e-1

So the conclusion survives contact with a better probe, and mildly favours the change:

  • Curved fields: p=1 is best in 2D and within 2% of best in 3D; the old 1/r⁴ is worst in both. Sharper weighting does not help a smooth field — it just narrows the effective stencil toward one point.
  • Discontinuity: this is where I expected 1/r⁴ to win, and it does not. In 2D all three saturate at 1.0 (a target between two materials gets a value near 0.5 either way). In 3D p=1 is clearly best, 8.30e-1 against 9.60e-1.
  • Linear: the old behaviour is marginally better in 2D (2.57e-2 vs 2.79e-2) and worse in 3D. This is the flattering-metric case the review called out — the differences are noise against a scheme that cannot reproduce linears at all.

That removes the substantive objection. The remaining ones stand: this still changes numbers for order=0 consumers I have not individually exercised (evaluate's fallback rung, adaptivity mesh2mesh, the order=1 degenerate fallback), epsilon is still an unscaled absolute length floor with no coincident-point test, and the bugfix and the default change are still bundled in one commit.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants