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

docs/developer/subsystems/interpolation.md

@@ -25,6 +25,11 @@ That is what `order` selects on
 
 $$ w_j = \frac{d_j^{-p}}{\sum_k d_k^{-p}} $$
 
+$d_j$ is the **actual** distance to the neighbour, and `p` defaults to 1 — so
+the default really is inverse distance. (Until #427 the kd-tree's squared
+distances were used directly, making the decay $r^{-2p}$: the documented
+default of `p=2` was in fact $1/r^4$.)
+
 Weights are positive and sum to one. Consequences, both of them important:
 
 - A **constant** field is reproduced exactly.

src/underworld3/ckdtree.pyx

@@ -434,7 +434,7 @@ cdef class KDTree:
             coords,
             data,
             nnn = None,
-            p=2,
+            p=1,
             verbose = False,
             order = 0,
             monotone = False,
@@ -476,7 +476,7 @@ cdef class KDTree:
             nearest-neighbour values without distance weighting.
         p : int, optional
             Power index for distance weighting: ``weight = 1/distance^p``
-            (default 2). Used by ``order=0`` only.
+            (default 1, i.e. inverse distance). Used by ``order=0`` only.
         verbose : bool, optional
             Print progress messages (default False).
         order : int, optional
@@ -738,13 +738,15 @@ cdef class KDTree:
                 )
             return closest_n, np.ones(closest_n.shape), None, degenerate
 
-        # can decompose weighting vecotrs as IDW is a linear relationship
-        # build normalise weight vectors and multiply that with known data
-        # TODO(BUG): issue #427 — `distance_n` holds SQUARED distances, so the
-        # decay is r^(-2p), not the documented r^(-p), and `epsilon` floors r
-        # at ~1e-6 rather than 1e-12.
+        # Inverse distance weighting: w = 1 / (eps + r)^p, normalised.
+        # `find_closest_n_points` returns SQUARED distances, so the square root
+        # is taken here -- without it the decay is r^(-2p) rather than the
+        # r^(-p) the argument names (issue #427). epsilon is a length floor for
+        # a target sitting exactly on a source point; it is on r, not r^2, so
+        # its scale is the one it reads as.
         epsilon = 1e-12
-        weights = 1 / np.power(epsilon + distance_n[:], p)
+        distance = np.sqrt(distance_n[:])
+        weights = 1 / np.power(epsilon + distance, p)
         n_weights = (weights.T / np.sum(weights, axis=1)).T
 
         if order == 0:
@@ -791,7 +793,7 @@ cdef class KDTree:
 
         return closest_n, linear_weights, wide, degenerate
 
-    def interpolation_matrix(self, coords, nnn=None, p=2, order=0):
+    def interpolation_matrix(self, coords, nnn=None, p=1, order=0):
         """Sparse operator mapping values on the KD-tree points to ``coords``.
 
         ``T @ data`` is exactly what :meth:`rbf_interpolator_local` returns for
@@ -895,8 +897,8 @@ cdef class KDTree:
         nnn : int
             The number of neighbour points to sample from. If `1`, no distance averaging is done.
         p : int
-            The power index to calculate weights, i.e., pow(distance, -p).
-            Used by ``order=0`` only.
+            The power index to calculate weights, i.e., pow(distance, -p),
+            on the actual distance. Used by ``order=0`` only.
         verbose : bool
             Print when mapping occurs
         order : int, optional

src/underworld3/discretisation/discretisation_mesh_variables.py

@@ -886,7 +886,7 @@ def _get_kdtree(self):
 
         return self._kdtree
 
-    def rbf_interpolate(self, new_coords, nnn=None, p=2, verbose=False,
+    def rbf_interpolate(self, new_coords, nnn=None, p=1, verbose=False,
                         order=0, monotone=False):
         """Interpolate variable data to new coordinates using RBF.
 
@@ -907,7 +907,8 @@ def rbf_interpolate(self, new_coords, nnn=None, p=2, verbose=False,
         nnn : int, optional
             Number of nearest neighbours (default: 4 for 3D, 3 for 2D).
         p : float, optional
-            Power parameter for inverse distance weighting (default: 2).
+            Power parameter for inverse distance weighting on the actual
+            distance (default: 1, i.e. inverse distance).
         verbose : bool, optional
             Print progress information.
         order : int, optional
@@ -1352,7 +1353,7 @@ def read_timestep(
             # ``nnn=1`` — exact match for round-trip reads, sensible
             # nearest-neighbour fallback for cross-mesh reads.
             result.array[:, 0, :] = kdt.rbf_interpolator_local(
-                local_query, landed_D, 1, 2, verbose
+                local_query, landed_D, nnn=1, verbose=verbose
             )
         elif local_query.shape[0] > 0:
             # No saved data landed on this rank — leave query payload zero

src/underworld3/swarm.py

@@ -1518,7 +1518,7 @@ def rbf_interpolate(self, new_coords, verbose=False, nnn=None, order=1,
         D = raw_data.copy()
         kdt = self.swarm._get_kdtree()
         values = kdt.rbf_interpolator_local(
-            new_coords, D, nnn, 2, verbose, order=order, monotone=monotone
+            new_coords, D, nnn, verbose=verbose, order=order, monotone=monotone
         )
 
         return values

tests/test_0102_kdtree_linear_rbf.py

@@ -424,16 +424,23 @@ def test_kdtree_neighbours_match_brute_force_on_random_data(dim):
 # --------------------------------------------------------------------------
 @pytest.mark.parametrize("dim", [2, 3])
 @pytest.mark.parametrize("order", [0, 1])
-def test_interpolation_matrix_agrees_with_the_value_path(dim, order):
-    """`T @ data` must be what the value API returns, or they will drift."""
+@pytest.mark.parametrize("nnn", [None, 8])
+@pytest.mark.parametrize("p", [1, 2])
+def test_interpolation_matrix_agrees_with_the_value_path(dim, order, nnn, p):
+    """`T @ data` must be what the value API returns, or they will drift.
+
+    Sweeps `nnn` and `p` as well as `order`: the earlier version pinned only
+    `order`, and hard-coded `p=2` on the value side, so it silently compared
+    two different weightings as soon as the default `p` changed.
+    """
     rng = np.random.default_rng(808 + dim)
     source = rng.random((500, dim))
     target = 0.1 + 0.8 * rng.random((60, dim))
     values = rng.standard_normal((source.shape[0], 2))
 
     kdt = uw.kdtree.KDTree(source)
-    T = kdt.interpolation_matrix(target, order=order)
-    direct = kdt.rbf_interpolator_local(target, values, None, 2, False, order=order)
+    T = kdt.interpolation_matrix(target, nnn=nnn, p=p, order=order)
+    direct = kdt.rbf_interpolator_local(target, values, nnn, p, False, order=order)
 
     assert T.shape == (target.shape[0], source.shape[0])
     assert np.abs(T @ values - direct).max() < 1.0e-12
@@ -465,6 +472,44 @@ def test_both_entry_points_reject_nnn_1_with_order_1(dim):
     assert kdt.interpolation_matrix(target, nnn=1).shape == (15, 200)
 
 
+@pytest.mark.parametrize("p", [1, 2, 3])
+def test_inverse_distance_decays_as_the_named_power(p):
+    """`p` must apply to the distance, not its square (issue #427).
+
+    The kd-tree returns squared distances, and the weighting used them
+    directly, so the decay was r^(-2p) while the argument was documented as
+    r^(-p). Nothing pinned the exponent, so it went unnoticed through a
+    rewrite.
+
+    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 returns the near
+    point's normalised weight, so w1/w2 is recoverable and must be 2**p.
+    """
+    source = np.array([[1.0, 0.0], [2.0, 0.0]])
+    target = np.array([[0.0, 0.0]])
+    data = np.array([[1.0], [0.0]])
+
+    kdt = uw.kdtree.KDTree(source)
+    near = kdt.rbf_interpolator_local(target, data, 2, p, False)[0, 0]
+    ratio = near / (1.0 - near)
+
+    assert np.isclose(ratio, 2.0 ** p, rtol=1e-10), (
+        f"weight ratio {ratio:.4f} implies decay r^-{np.log2(ratio):.2f}, "
+        f"expected r^-{p}"
+    )
+
+
+def test_inverse_distance_default_is_inverse_distance():
+    """The default is p=1 — genuinely inverse distance, not inverse square."""
+    source = np.array([[1.0, 0.0], [2.0, 0.0]])
+    target = np.array([[0.0, 0.0]])
+    data = np.array([[1.0], [0.0]])
+
+    kdt = uw.kdtree.KDTree(source)
+    near = kdt.rbf_interpolator_local(target, data, 2)[0, 0]
+    assert np.isclose(near / (1.0 - near), 2.0, rtol=1e-10)
+
+
 def test_stencil_larger_than_the_cloud_reports_what_went_wrong():
     """The old message named a function the caller never invoked."""
     rng = np.random.default_rng(17)