rbf_interpolator_local: inverse-distance weights use SQUARED distances, so the documented p semantics are wrong (effective exponent is 2p)

[lmoresi]

Found while adding the linear-exact local RBF interpolator (feature/linear-rbf).

The defect

KDTree.rbf_interpolator_local_from_kdtree in src/underworld3/ckdtree.pyx:

distance_n, closest_n = self.query(coords, k=nnn)   # sqr_dists=True by DEFAULT
...
epsilon = 1e-12
weights = 1 / np.power(epsilon + distance_n[:], p)

KDTree.query returns squared distances unless sqr_dists=False is passed, and it is not. So the weight is

w = 1 / (eps + r^2)^p

i.e. the effective decay is r^(-2p). At the documented default p=2 the scheme is inverse fourth power, not inverse square — much sharper, and close to nearest-neighbour.

The docstring says:

p : int — The power index to calculate weights, i.e., pow(distance, -p)

The retained old_rbf_interpolator_local_from_kdtree took np.sqrt(distance_n[:, j]) explicitly and used a fixed 1/r, so the semantics changed silently when the method was rewritten.

Second defect, same expression

epsilon = 1e-12 is added to a squared distance, so the coincident-point regularisation floor is at r ~ 1e-6 in length units, not 1e-6^2. The old code used epsilon = 1e-9 on the actual distance.

Third, cosmetic

coords_contiguous = np.ascontiguousarray(coords_converted)   # never used
...
distance_n, closest_n = self.query(coords, k=nnn)            # redoes the conversion

coords_contiguous is dead, and unit-aware input is converted twice.

Why it has not been noticed

tests/test_0505_rbf_swarm_mesh.py is the only accuracy test and uses rtol=5e-2, loose enough that 1/r^2 vs 1/r^4 does not show. It also never varies p — every parametrised case is p=2.

Deciding the fix

This is a behaviour change either way, so it needs a decision rather than a quiet correction:

• Fix the code (take sqrt, or halve the exponent) — matches the documentation and the pre-rewrite behaviour, but moves every existing order=0 result.
• Fix the documentation — states what the code actually does, costs nothing, but leaves p meaning something non-obvious.

Nothing in feature/linear-rbf changes this; the order=0 path there is bit-identical to today.

[lmoresi]

Fixed in #448.