From a4f43df7d8c44e723b66c65671dadd826e617cb3 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 27 Jul 2026 10:45:50 +1000 Subject: [PATCH 1/5] Local RBF interpolator with a linear-reproduction guarantee UW3's local scattered-point interpolator was inverse-distance (Shepard) weighting. Its weights are positive and sum to one, so it reproduces a constant exactly but not a linear field: anything with a gradient is smeared, and the error does not vanish as the stencil tightens. That interpolator refreshes the swarm proxy mesh variables, which is what integration and derivative evaluation read. Adds a second scheme, selected by order=: order=0 inverse distance, unchanged and still the KDTree default order=1 polyharmonic (r^2 log r) with an affine tail, so both sum(w) = 1 and sum(w_j x_j) = x* hold by construction The saddle system is solved per target point on its own kNN stencil in a local scaled frame, batched and chunked. Cost is one small dense solve per point; the result is still sparse at nnn non-zeros per row. The codebase already had a linear-exact scheme in custom_mg.rbf_prolongation, but global and dense; this is the sparse equivalent. Measured, swarm proxy against an analytic field (a linear field lies exactly inside both the P1 and P2 proxy space, so all of the error is particle-to-node transfer): 2D, cellSize 1/8, fill 4, degree 1 linear max 4.6e-3 -> 5.1e-16 quadratic max 9.2e-3 -> 1.0e-4 rms 2.1e-3 -> 2.1e-5 3D, cellSize 1/8, fill 4, degree 1 linear max 4.8e-3 -> 5.5e-16 quadratic max 8.5e-3 -> 2.4e-4 rms 2.0e-3 -> 2.3e-5 Widening the order=0 stencil does not close the gap: at nnn=12 the inverse-distance error is unchanged to two significant figures, because it is set by the asymmetry of the neighbourhood rather than its size. Swarm proxy defaults therefore change to order=1, nnn=2*(dim+1). Left deliberately on order=0: MeshVariable.rbf_interpolate, because it is the RBF rung of the evaluate point-location ladder whose documented contract is that it is bounded; read_timestep, which wants an exact nnn=1 round-trip; and IndexSwarmVariable, whose level sets must stay in [0,1]. Degenerate stencils (collinear in 2D, coplanar in 3D) are pre-screened by the rank of the affine block, the computed weights are then validated against the two reproduction identities, and anything failing falls back to inverse distance with a single counted warning. Never NaN, never silent. order=1 requires nnn >= dim+2 and says why if it is not: at dim+1 the RBF block is inert and the scheme is bare barycentric interpolation on the neighbour simplex. monotone= clamps to the stencil range for consumers that need bounded values, reusing the evaluator's vocabulary. Note it discards linear exactness wherever the target sits outside its neighbours' convex hull, so it is not a free safety net. Defects found in passing and filed rather than fixed here: #426 (proxy refresh broken under an active units model), #427 (inverse-distance p semantics), #428 (dead parameters and dead code), #429 (custom_mg RBF prolongation conditioning). Underworld development team with AI support from Claude Code --- docs/developer/index.md | 2 + docs/developer/subsystems/data-access.md | 4 + docs/developer/subsystems/interpolation.md | 187 +++++++++++++ src/underworld3/ckdtree.pyx | 121 ++++++++- .../discretisation_mesh_variables.py | 23 +- src/underworld3/swarm.py | 57 +++- src/underworld3/utilities/rbf_stencil.py | 201 ++++++++++++++ .../test_0776_linear_rbf_proxy_parallel.py | 101 +++++++ tests/test_0102_kdtree_linear_rbf.py | 254 ++++++++++++++++++ tests/test_0506_swarm_proxy_linear_exact.py | 120 +++++++++ 10 files changed, 1046 insertions(+), 24 deletions(-) create mode 100644 docs/developer/subsystems/interpolation.md create mode 100644 src/underworld3/utilities/rbf_stencil.py create mode 100644 tests/parallel/test_0776_linear_rbf_proxy_parallel.py create mode 100644 tests/test_0102_kdtree_linear_rbf.py create mode 100644 tests/test_0506_swarm_proxy_linear_exact.py diff --git a/docs/developer/index.md b/docs/developer/index.md index eed57fbf4..a94685c67 100644 --- a/docs/developer/index.md +++ b/docs/developer/index.md @@ -35,6 +35,7 @@ same topic are reference or historical material subordinate to the governing doc |-------|--------------------| | Coding style & API conventions | [UW3 Style Charter](UW3_STYLE_CHARTER.md) (detailed reference: [Style Guide](UW3_Style_and_Patterns_Guide.md)) | | Data access | [subsystems/data-access.md](subsystems/data-access.md) (internals reference: [NDArray System](UW3_Developers_NDArrays.md)) | +| Local scattered-point interpolation | [subsystems/interpolation.md](subsystems/interpolation.md) | | Units | [design/UNITS_SIMPLIFIED_DESIGN_2025-11.md](design/UNITS_SIMPLIFIED_DESIGN_2025-11.md) | | Testing tiers | [TESTING-RELIABILITY-SYSTEM.md](TESTING-RELIABILITY-SYSTEM.md) | | Branching & releases | [guides/branching-strategy.md](guides/branching-strategy.md) | @@ -188,6 +189,7 @@ subsystems/constitutive-models-theory subsystems/constitutive-models-anisotropy subsystems/swarm-system subsystems/data-access +subsystems/interpolation subsystems/expressions-functions subsystems/containers subsystems/checkpointing-system diff --git a/docs/developer/subsystems/data-access.md b/docs/developer/subsystems/data-access.md index cd1ef76a0..764d57315 100644 --- a/docs/developer/subsystems/data-access.md +++ b/docs/developer/subsystems/data-access.md @@ -103,6 +103,10 @@ symbolic_expr = swarm_var.sym # Triggers RBF interpolation if stale This pattern prevents circular dependencies when callbacks would trigger nested PETSc field access. +The interpolation itself — which fields the proxy transfer reproduces exactly, +how to choose `nnn`, and when boundedness matters more than accuracy — is +covered in [Local Scattered-Point Interpolation](interpolation.md). + ### Vector Availability Variables set `_available=True` by default, ensuring solvers can access vectors without modification. Lazy initialization creates vectors on first access: diff --git a/docs/developer/subsystems/interpolation.md b/docs/developer/subsystems/interpolation.md new file mode 100644 index 000000000..4d40e535d --- /dev/null +++ b/docs/developer/subsystems/interpolation.md @@ -0,0 +1,187 @@ +--- +title: "Local Scattered-Point Interpolation" +--- + +## Purpose + +Underworld needs to move field values between point sets that have no shared +mesh topology: particles to proxy mesh nodes, nodes to arbitrary evaluation +points, one mesh to another after adaptation. The tool for that is a *local* +interpolant — for each target point, take its `nnn` nearest source points from +a kd-tree and form a weighted average. + +The only question that matters is **which fields the weights reproduce +exactly**. That is what `order` selects on +{py:meth}`underworld3.ckdtree.KDTree.rbf_interpolator_local`. + +## The two schemes + +### `order=0` — inverse distance (Shepard) + +$$ w_j = \frac{d_j^{-p}}{\sum_k d_k^{-p}} $$ + +Weights are positive and sum to one. Consequences, both of them important: + +- A **constant** field is reproduced exactly. +- The result is a **convex combination** of the neighbouring values, so it can + never overshoot them. The interpolant is intrinsically bounded. +- A field with a **gradient is smeared**, and the error does *not* vanish as + the stencil tightens. Adding neighbours does not help; refining the point + spacing helps only at first order. + +### `order=1` — polyharmonic RBF with an affine tail + +For a target $x^*$ with neighbours $x_j$, solve the small saddle-point system + +$$ +\begin{bmatrix} A & P \\ P^{T} & 0 \end{bmatrix} +\begin{bmatrix} w \\ \lambda \end{bmatrix} += +\begin{bmatrix} \varphi(|x^* - x_j|) \\ 1,\; x^* \end{bmatrix} +$$ + +with $A_{ij} = \varphi(|x_i - x_j|)$, the polyharmonic (thin-plate) kernel +$\varphi(r) = r^2 \log r$, and the affine tail $P = [1,\; x_j]$. + +The lower block *is* the reproduction constraint $P^{T} w = [1,\; x^*]$, so + +$$ \sum_j w_j = 1, \qquad \sum_j w_j x_j = x^* $$ + +hold to round-off by construction: **constants and linear fields are exact**. + +Weights are signed, so the interpolant can overshoot the neighbouring values — +the price of giving up the convex combination. + +Cost is one dense $(nnn + d + 1)^3$ solve per target point, and the result is +still sparse: `nnn` non-zeros per row. The weights depend only on geometry, so +one solve serves every data component. + +## Measured + +Relative error of a swarm proxy variable against the analytic field, at the +proxy's own nodal coordinates. A **linear** field lies exactly inside both the +P1 and P2 proxy spaces, so the discretisation contributes nothing and every bit +of the measured error is particle-to-node transfer error. + +2D simplex box, `cellSize=1/8`, `fill_param=4`, proxy degree 1: + +| field | `order=0`, `nnn=3` | `order=1`, `nnn=6` | +|---|---|---| +| linear, max | 4.6e-3 | **5.1e-16** | +| quadratic, max | 9.2e-3 | 1.0e-4 | +| quadratic, rms | 2.1e-3 | 2.1e-5 | + +3D simplex box, `cellSize=1/8`, `fill_param=4`, proxy degree 1: + +| field | `order=0`, `nnn=4` | `order=1`, `nnn=8` | +|---|---|---| +| linear, max | 4.8e-3 | **5.5e-16** | +| quadratic, max | 8.5e-3 | 2.4e-4 | +| quadratic, rms | 2.0e-3 | 2.3e-5 | + +Two things to read off. First, the linear-field column is the clean statement: +one scheme is exact, the other is not. Second — and this is the part that +actually matters for physics — the improvement carries over to a field with +curvature, by roughly two orders of magnitude in rms. + +Widening the `order=0` stencil does **not** close the gap. At `nnn=12` the +inverse-distance error is unchanged to two significant figures, because the +error is set by the asymmetry of the neighbourhood, not by how many points are +in it. + +## Choosing `nnn` + +`order=1` requires **`nnn >= dim + 2`** and raises `ValueError` below that. + +The reason is worth knowing. At exactly `nnn == dim + 1` the affine block $P$ +is square, so the constraint $P^{T} w = [1, x^*]$ alone determines $w$: the RBF +block is inert and the scheme collapses to bare **barycentric interpolation on +the neighbour simplex**. That is singular whenever those `dim + 1` points are +collinear (2D) or coplanar (3D) — common near boundaries and on graded meshes. + +The swarm proxy default is `nnn = 2 * (dim + 1)` (6 in 2D, 8 in 3D), which +leaves enough slack that degenerate neighbourhoods are rare. + +## Degenerate stencils + +If the neighbours cannot support an affine fit, the saddle system is singular. +This is not hypothetical: it happens at boundaries and in graded meshes. + +The implementation (`underworld3.utilities.rbf_stencil.linear_exact_weights`): + +1. **Pre-screens** the affine block by SVD — degenerate if + $\sigma_{\min} < 10^{-8}\,\sigma_{\max}$. +2. **Solves** only the healthy stencils, batched and chunked over targets. +3. **Validates the answer** rather than guessing a condition number: it checks + the two reproduction identities and finiteness on the computed weights, and + demotes anything that fails. +4. **Falls back** to the inverse-distance weights for those points — less + accurate there, but finite and bounded. It never returns NaN. +5. **Warns once** with the count and fraction that fell back. A silent + geometric fallback is the failure mode behind issue #424; a per-point + warning would be unusable, a single counted one is not. + +## Boundedness: `monotone` + +`monotone=True` (or `"clamp"`, matching the vocabulary of +{py:func}`underworld3.function.evaluate`) bounds each result to the min/max of +its own stencil's source values, restoring the boundedness that `order=1` gives +up. + +```{warning} +Clamping and linear exactness are in genuine conflict. Where the target lies +outside the convex hull of its neighbours the exact linear value is *outside* +the stencil range, so the clamp discards it. + +Measured on the 2D case above, `order=1` with `monotone=True` gives a max error +of 4.0e-3 on a linear field — back to the inverse-distance level — while +keeping the interior nodes exact. Use it when a consumer genuinely requires +bounded values, not as a default safety net. +``` + +Consumers that do depend on boundedness, and are therefore deliberately left on +`order=0`: + +- The **RBF rung of the point-location fallback ladder** in + `uw.function.evaluate` — `docs/developer/design/point-location-capability.md` + states its contract as "bounded, topology-free, honest". + `MeshVariable.rbf_interpolate` therefore keeps `order=0` as its default. +- `IndexSwarmVariable` level sets, which carry material fractions in $[0,1]$ and + are divided through by their sum in `constitutive_models.py`. That path + hand-rolls its own inverse-distance weighting and is untouched. + +## Where each is used + +| Caller | Default | Why | +|---|---|---| +| `SwarmVariable.rbf_interpolate` / `_rbf_to_meshVar` (proxy refresh) | `order=1`, `nnn=2*(dim+1)` | Accuracy; the proxy feeds integration, derivatives and the Lagrangian time derivative | +| `MeshVariable.rbf_interpolate` | `order=0` | Documented bounded contract (see above) | +| `SwarmVariable.read_timestep`, `MeshVariable.read_timestep` | `nnn=1` | Exact checkpoint round-trip, no averaging wanted | +| `IndexSwarmVariable._update_proxy_variables` | its own IDW | Level-set fractions must stay in $[0,1]$ | + +A rank holding fewer particles than the affine tail needs drops to `order=0` +automatically rather than failing the refresh. + +## Parallel behaviour + +Proxy refresh is **rank-local**: the kd-tree indexes only this rank's +particles and there is no halo exchange (SWARM-15, see +`docs/developer/design/SWARM_MODERNIZATION_DESIGN_2026-07.md` §4). A proxy node +near a partition seam therefore gathers from a one-sided neighbourhood. + +Linear exactness improves this but does not fix it. 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 is np-independent +(pinned by `tests/parallel/test_0776_linear_rbf_proxy_parallel.py`). For a +field with curvature a one-sided stencil still differs from a centred one, so +np-dependence remains. The halo exchange in SWARM-15 is still the real fix. + +## Related + +- `docs/developer/subsystems/data-access.md` — lazy proxy updates +- `docs/developer/design/SWARM_MODERNIZATION_DESIGN_2026-07.md` §4 — rank-local seams +- `docs/developer/design/point-location-capability.md` — the evaluate fallback ladder +- `utilities/custom_mg.py` — a *global, dense* polyharmonic prolongation for + multigrid transfers. Same kernel and tail, but every coarse point enters every + row, so `PᵀAP` is dense. The local scheme documented here is the sparse + equivalent. diff --git a/src/underworld3/ckdtree.pyx b/src/underworld3/ckdtree.pyx index d8eb90b90..733d30bc1 100644 --- a/src/underworld3/ckdtree.pyx +++ b/src/underworld3/ckdtree.pyx @@ -34,6 +34,28 @@ def total_constructed(): return _total_constructed +def _normalise_monotone(monotone): + """Resolve the ``monotone`` argument to a bool, using UW3's one vocabulary. + + The canonical spelling lives with the evaluator + (``function.functions_unit_system._normalize_monotone``) and is reused here + so there is a single definition of what the word accepts. Only the + ``"clamp"`` mode has meaning for a local stencil: ``"pick"`` re-evaluates + out-of-bounds points through the FE path, which a kd-tree knows nothing + about. + """ + from underworld3.function.functions_unit_system import _normalize_monotone + + mode = _normalize_monotone(monotone) + if mode == "pick": + raise ValueError( + "monotone='pick' has no meaning for a local kd-tree stencil — it " + "re-evaluates out-of-bounds points through the finite-element path. " + "Use monotone='clamp' here, or uw.function.evaluate(..., monotone='pick')." + ) + return mode == "clamp" + + cdef class KDTree: """ Unit-aware KD-Tree for spatial indexing and queries. @@ -408,14 +430,29 @@ cdef class KDTree: nnn = 4, p=2, verbose = False, + order = 0, + monotone = False, ): """ - Interpolate data to target coordinates using inverse distance weighting. + Interpolate data from the KD-tree points to arbitrary target coordinates. + + Two local schemes are available, selected by ``order`` — the highest + degree of polynomial the weights reproduce exactly: - This is a convenience wrapper around :meth:`rbf_interpolator_local_from_kdtree`. - It performs radial basis function (RBF) interpolation using inverse distance - weighting to map known data values from the KD-tree points to arbitrary - target coordinates. + ``order=0`` (default) + Inverse-distance (Shepard) weighting. Weights are positive and sum + to one, so a **constant** field is reproduced exactly and the result + is a convex combination of the neighbouring values — bounded, but a + field with a gradient is smeared and the error does not vanish as + the stencil tightens. + + ``order=1`` + Polyharmonic RBF with an affine tail, so **constant and linear** + fields are reproduced exactly (:math:`\\sum_j w_j = 1` and + :math:`\\sum_j w_j x_j = x^*`). Weights may be negative, so the + result can overshoot the neighbouring values; see ``monotone``. + + Both are local: ``nnn`` non-zero weights per target point. Parameters ---------- @@ -428,10 +465,20 @@ cdef class KDTree: nnn : int, optional Number of nearest neighbours to use for interpolation (default 4). If 1, returns raw nearest-neighbour values without distance weighting. + ``order=1`` requires ``nnn >= dim + 2``. p : int, optional - Power index for distance weighting: ``weight = 1/distance^p`` (default 2). + Power index for distance weighting: ``weight = 1/distance^p`` + (default 2). Used by ``order=0`` only. verbose : bool, optional Print progress messages (default False). + order : int, optional + Polynomial reproduction order, 0 (default) or 1. + monotone : bool or str, optional + Bound each interpolated value to the min/max of its own stencil's + source values. ``False`` (default) or ``True`` / ``"clamp"``. + Restores boundedness at ``order=1``, but note that clamping + necessarily discards linear exactness wherever the target lies + outside the convex hull of its neighbours (boundary points). Returns ------- @@ -444,7 +491,7 @@ cdef class KDTree: query : Find nearest neighbours without interpolation. """ return self.rbf_interpolator_local_from_kdtree( - coords, data, nnn, p, verbose, + coords, data, nnn, p, verbose, order, monotone, ) def old_rbf_interpolator_local_from_kdtree(self, @@ -576,9 +623,10 @@ cdef class KDTree: return Values - def rbf_interpolator_local_from_kdtree(self, coords, data, nnn, p, verbose): + def rbf_interpolator_local_from_kdtree(self, coords, data, nnn, p, verbose, + order=0, monotone=False): """ - Performs an inverse distance (squared) mapping of data to the target `coords`. + Map data held on the KD-tree points onto the target `coords`. This method is unit-aware: if the KD-tree was built with unit-aware coordinates, it will automatically convert query coordinates to match before interpolation. @@ -596,15 +644,29 @@ 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) + The power index to calculate weights, i.e., pow(distance, -p). + Used by ``order=0`` only. verbose : bool Print when mapping occurs + order : int, optional + Polynomial reproduction order: 0 for inverse-distance (constants + exact), 1 for polyharmonic + affine tail (constants and linears + exact). See :meth:`rbf_interpolator_local`. + monotone : bool or str, optional + ``False``, ``True`` or ``"clamp"`` — bound each result to the + min/max of its own stencil's source values. Returns ------- ndarray Interpolated data values at target coordinates """ + if order not in (0, 1): + raise ValueError( + f"order must be 0 (inverse distance) or 1 (linear-exact), got {order!r}." + ) + monotone_clamp = _normalise_monotone(monotone) + # Convert coordinates to match tree's coordinate system coords_converted = self._convert_coords_to_tree_units(coords) @@ -638,19 +700,58 @@ cdef class KDTree: if nnn == 1: # only use nearest neighbour raw data + if order == 1: + raise ValueError( + "order=1 needs at least dim + 2 neighbours to determine the " + f"affine tail; nnn=1 selects the raw nearest-neighbour path." + ) return data[closest_n] # 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 (query + # defaults to sqr_dists=True), so the decay is r^(-2p), not the + # documented r^(-p), and `epsilon` floors r at ~1e-6 rather than 1e-12. epsilon = 1e-12 weights = 1 / np.power(epsilon + distance_n[:], p) n_weights = (weights.T / np.sum(weights, axis=1)).T kdata = data[closest_n[:]] + if order == 1: + from underworld3.utilities.rbf_stencil import linear_exact_weights + + stencil_coords = np.asarray(self.points)[closest_n[:]] + linear_weights, degenerate = linear_exact_weights( + coords_converted, stencil_coords + ) + # A stencil that cannot support an affine fit (collinear in 2D, + # coplanar in 3D) keeps the inverse-distance weights: less accurate + # there, but finite and bounded. Report it once — a silent + # geometric fallback is the failure mode of issue #424. + linear_weights[degenerate] = n_weights[degenerate] + n_weights = linear_weights + + n_degenerate = int(degenerate.sum()) + if n_degenerate: + import warnings + + warnings.warn( + f"rbf_interpolator_local(order=1): {n_degenerate} of " + f"{degenerate.size} stencils could not support an affine fit " + "(collinear/coplanar neighbours) and fell back to inverse-" + "distance weighting. Increasing nnn usually removes them.", + stacklevel=2, + ) + # magic with einstein summation power vals = np.einsum("sdc,sd->sc", kdata, n_weights) # print(valz) + if monotone_clamp: + # Bound to the stencil's own range. This is what restores the + # convex-combination property that order=1 gives up. + vals = np.clip(vals, kdata.min(axis=1), kdata.max(axis=1)) + if verbose and uw.mpi.rank == 0: print(f"Mapping values ... finished", flush=True) diff --git a/src/underworld3/discretisation/discretisation_mesh_variables.py b/src/underworld3/discretisation/discretisation_mesh_variables.py index 1a71d2c45..3e6bdabcc 100644 --- a/src/underworld3/discretisation/discretisation_mesh_variables.py +++ b/src/underworld3/discretisation/discretisation_mesh_variables.py @@ -886,24 +886,41 @@ def _get_kdtree(self): return self._kdtree - def rbf_interpolate(self, new_coords, meth=0, p=2, verbose=False, nnn=None, rubbish=None): + def rbf_interpolate(self, new_coords, meth=0, p=2, verbose=False, nnn=None, + rubbish=None, order=0, monotone=False): """Interpolate variable data to new coordinates using RBF. Uses inverse distance weighting with k-nearest neighbors to interpolate values from mesh nodes to arbitrary coordinates. + The default ``order=0`` is deliberate and differs from the swarm proxy + path: this method is the RBF rung of the point-location fallback ladder + in :func:`underworld3.function.evaluate`, whose documented contract is + that it is *bounded*. Inverse-distance weights are a convex combination + and so cannot overshoot; ``order=1`` weights can. Pass ``order=1`` + explicitly where linear exactness matters more than boundedness. + Parameters ---------- new_coords : numpy.ndarray Target coordinates of shape ``(n_points, dim)``. meth : int, optional Interpolation method (reserved, currently unused). + TODO(BUG): issue #428 — ``meth`` and ``rubbish`` are dead + parameters, and ``tests/test_0505_rbf_swarm_mesh.py`` passes its + ``nnn`` into ``meth`` positionally, so that test silently + discards it. p : float, optional Power parameter for inverse distance weighting (default: 2). verbose : bool, optional Print progress information. nnn : int, optional Number of nearest neighbors (default: 4 for 3D, 3 for 2D). + order : int, optional + Polynomial reproduction order, 0 (default, bounded) or 1 + (constants and linears exact; requires ``nnn >= dim + 2``). + monotone : bool or str, optional + Bound each value to the min/max of its own stencil. Returns ------- @@ -928,7 +945,9 @@ def rbf_interpolate(self, new_coords, meth=0, p=2, verbose=False, nnn=None, rubb # Use cached KDTree for interpolation kdt = self._get_kdtree() - values = kdt.rbf_interpolator_local(new_coords, D, nnn, p=p, verbose=verbose) + values = kdt.rbf_interpolator_local( + new_coords, D, nnn, p=p, verbose=verbose, order=order, monotone=monotone + ) return values diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py index da80015ed..65a67840d 100644 --- a/src/underworld3/swarm.py +++ b/src/underworld3/swarm.py @@ -1140,9 +1140,15 @@ def _update_proxy_if_stale(self): return # Maybe rbf_interpolate for this one and meshVar is a special case - def _rbf_to_meshVar(self, meshVar, nnn=None, verbose=False): + def _rbf_to_meshVar(self, meshVar, nnn=None, verbose=False, order=1, + monotone=False): """ - Here is how it works: for each particle, create a distance-weighted average on the node data + Refresh a proxy mesh variable from the particles. + + Each proxy node gathers from its ``nnn`` nearest particles. The + default weights reproduce linear fields exactly (``order=1``), so a + field with a uniform gradient transfers without smearing; ``nnn`` and + ``order`` are resolved in :meth:`rbf_interpolate`. Todo: caching the k-d trees etc for the proxy-mesh-variable nodal points Todo: some form of global fall-back for when there are no particles on a processor @@ -1151,9 +1157,6 @@ def _rbf_to_meshVar(self, meshVar, nnn=None, verbose=False): # Mapping to the coordinates of the variable from the # particle coords - if nnn is None: - nnn = self.swarm.mesh.dim + 1 - if meshVar.mesh != self.swarm.mesh: # If this is our own proxy variable and mesh has changed, recreate it if hasattr(self, "_meshVar") and meshVar is self._meshVar: @@ -1163,6 +1166,11 @@ def _rbf_to_meshVar(self, meshVar, nnn=None, verbose=False): else: raise RuntimeError("Cannot map a swarm to a different mesh") + # TODO(BUG): issue #426 — MeshVariable.coords is dimensionalised when + # the model has reference quantities active, but the swarm kd-tree is + # built from non-dimensional particle coordinates, so this raises from + # ckdtree._convert_coords_to_tree_units. Should be `coords_nd`. + # Same defect at _update_proxy_variables (IndexSwarmVariable). new_coords = meshVar.coords # Starved-rank guard (SWARM-07): with <= 1 local particles there is @@ -1190,7 +1198,9 @@ def _rbf_to_meshVar(self, meshVar, nnn=None, verbose=False): ) Values = current_values else: - Values = self.rbf_interpolate(new_coords, verbose=verbose, nnn=nnn) + Values = self.rbf_interpolate( + new_coords, verbose=verbose, nnn=nnn, order=order, monotone=monotone + ) meshVar.data[...] = Values[...] @@ -1479,12 +1489,16 @@ def _object_viewer(self): display(self.data), return - def rbf_interpolate(self, new_coords, verbose=False, nnn=None): + def rbf_interpolate(self, new_coords, verbose=False, nnn=None, order=1, + monotone=False): """ Radial basis function interpolation of particle data to arbitrary points. - Uses inverse-distance weighting to interpolate particle values - to new coordinate locations. + By default this reproduces constant *and linear* fields exactly + (``order=1``): a polyharmonic kernel with an affine tail over the + ``nnn`` nearest particles. Inverse-distance weighting (``order=0``) + reproduces only constants, so any field with a gradient is smeared by + an error that does not vanish as the particles crowd together. Parameters ---------- @@ -1493,7 +1507,18 @@ def rbf_interpolate(self, new_coords, verbose=False, nnn=None): verbose : bool, default=False Print diagnostic information during interpolation. nnn : int, optional - Number of nearest neighbors to use. Defaults to ``mesh.dim + 1``. + Number of nearest neighbours to use. Defaults to + ``2 * (mesh.dim + 1)`` — comfortably above the ``dim + 2`` that the + affine tail needs, so that near-degenerate particle neighbourhoods + do not have to fall back. + order : int, default=1 + Polynomial reproduction order: 1 (constants and linears exact) or + 0 (constants only, inverse distance). Drops to 0 automatically on a + rank holding too few particles to determine the affine tail. + monotone : bool or str, default=False + Bound each value to the min/max of its own particle stencil. Note + this discards linear exactness wherever the target sits outside the + convex hull of its neighbours. Returns ------- @@ -1528,15 +1553,23 @@ def rbf_interpolate(self, new_coords, verbose=False, nnn=None): return np.zeros((new_coords.shape[0], data_size[1])) if nnn is None: - nnn = self.swarm.mesh.dim + 1 + nnn = 2 * (self.swarm.mesh.dim + 1) if nnn > data_size[0]: nnn = data_size[0] + # A rank holding fewer particles than the affine tail needs cannot + # support a linear fit at all. Inverse distance still gives a sensible + # answer there, so degrade rather than fail the whole refresh. + if order == 1 and nnn < self.swarm.mesh.dim + 2: + order = 0 + # Use direct PETSc access to avoid callback circular dependency D = raw_data.copy() kdt = self.swarm._get_kdtree() - values = kdt.rbf_interpolator_local(new_coords, D, nnn, 2, verbose) + values = kdt.rbf_interpolator_local( + new_coords, D, nnn, 2, verbose, order=order, monotone=monotone + ) return values diff --git a/src/underworld3/utilities/rbf_stencil.py b/src/underworld3/utilities/rbf_stencil.py new file mode 100644 index 000000000..1eaa00974 --- /dev/null +++ b/src/underworld3/utilities/rbf_stencil.py @@ -0,0 +1,201 @@ +"""Local RBF-FD stencil weights that reproduce linear fields exactly. + +Inverse-distance (Shepard) weights are positive and sum to one, so they +reproduce a *constant* exactly but not a linear field: any field with a +gradient is smeared, and the error does not vanish as the stencil tightens. + +The weights built here do reproduce linears. For a target point +:math:`x^*` and its :math:`m` nearest source points :math:`x_j`, solve the +small saddle-point system + +.. math:: + + \\begin{bmatrix} A & P \\\\ P^{T} & 0 \\end{bmatrix} + \\begin{bmatrix} w \\\\ \\lambda \\end{bmatrix} + = \\begin{bmatrix} \\varphi(|x^* - x_j|) \\\\ 1,\\; x^* \\end{bmatrix} + +with :math:`A_{ij} = \\varphi(|x_i - x_j|)`, the polyharmonic (thin-plate) +kernel :math:`\\varphi(r) = r^2 \\log r`, and the affine tail +:math:`P = [1,\\; x_j]`. The lower block *is* the reproduction constraint +:math:`P^{T} w = [1,\\; x^*]`, so + +.. math:: + + \\sum_j w_j = 1, \\qquad \\sum_j w_j x_j = x^* + +hold to round-off by construction — constants and linears are exact. The +weights depend only on geometry, so one solve serves every data component. + +Cost is one dense :math:`(m + d + 1)^3` solve per target point, and the +result is sparse: :math:`m` non-zeros per row. + +This module is deliberately pure NumPy with no PETSc dependency; it is the +numeric core behind ``KDTree.rbf_interpolator_local(..., order=1)``. +""" + +import numpy as np + +# The kernel is not scale-invariant: under r -> s.r it gains an r^2 log(s) +# term, which is quadratic and therefore NOT absorbed by the affine tail. So +# every stencil is solved in a local frame (target at the origin, farthest +# neighbour at unit radius); conditioning is then independent of the absolute +# stencil size, and both tolerances below are dimensionless. +_RANK_TOL = 1.0e-8 # smallest/largest singular value of the affine block +_REPRODUCTION_TOL = 1.0e-6 # residual of the two reproduction identities + +# Peak working set per chunk of the stacked saddle matrices. Sized for memory, +# not speed: the full (N, m+d+1, m+d+1) stack for a million targets would be +# gigabytes. +_CHUNK_BYTES = 64 * 1024 * 1024 + + +def _polyharmonic(r): + """:math:`\\varphi(r) = r^2 \\log r`, with :math:`\\varphi(0) = 0`.""" + # The clip only keeps log() finite at coincident points; r**2 drives the + # product to zero there regardless of the value substituted. + safe = np.where(r > 0.0, r, 1.0) + return np.where(r > 0.0, r ** 2 * np.log(safe), 0.0) + + +def _solve_stencils(matrices, rhs): + """Batched saddle-point solve, tolerant of an exactly singular member. + + ``np.linalg.solve`` raises for the whole batch if any one matrix is + singular, and gives no way to tell which. The pre-screen in + :func:`linear_exact_weights` removes the common case; this catches the + residue by falling back to the pseudo-inverse for that chunk, which never + raises. The post-solve reproduction check then rejects anything the + pseudo-inverse merely papered over. + """ + try: + return np.linalg.solve(matrices, rhs[..., None])[..., 0] + except np.linalg.LinAlgError: + return (np.linalg.pinv(matrices) @ rhs[..., None])[..., 0] + + +def linear_exact_weights(target_coords, neighbour_coords): + """Stencil weights that reproduce constant and linear fields exactly. + + Parameters + ---------- + target_coords : numpy.ndarray + Points to interpolate to, shape ``(n_targets, dim)``. + neighbour_coords : numpy.ndarray + Source points of each target's stencil, shape + ``(n_targets, nnn, dim)``. Normally the ``nnn`` nearest neighbours + from a kd-tree query. + + Returns + ------- + weights : numpy.ndarray + Shape ``(n_targets, nnn)``. Interpolate with + ``(weights[:, :, None] * data[stencil]).sum(axis=1)``. + degenerate : numpy.ndarray + Boolean, shape ``(n_targets,)``. True where the stencil could not + support an affine fit — collinear neighbours in 2D, coplanar in 3D, + or a stencil that collapsed onto its target. **Rows flagged here are + returned as zeros**; the caller must substitute its own fallback + (Shepard weights over the same neighbours are the natural choice). + + Notes + ----- + ``nnn`` must be at least ``dim + 2``. At exactly ``dim + 1`` the affine + block is square, the constraint alone determines the weights, and the + scheme degenerates to bare barycentric interpolation on the neighbour + simplex — which is singular whenever those points are collinear or + coplanar, a common situation near boundaries and in graded meshes. + """ + target_coords = np.ascontiguousarray(target_coords, dtype=np.float64) + neighbour_coords = np.ascontiguousarray(neighbour_coords, dtype=np.float64) + + n_targets, nnn, dim = neighbour_coords.shape + if target_coords.shape != (n_targets, dim): + raise ValueError( + f"target_coords has shape {target_coords.shape}, expected " + f"({n_targets}, {dim}) to match neighbour_coords " + f"{neighbour_coords.shape}." + ) + if nnn < dim + 2: + raise ValueError( + f"A linear-exact stencil needs at least dim + 2 = {dim + 2} " + f"neighbours, got nnn = {nnn}. At dim + 1 the affine tail is " + "exactly determined and the scheme reduces to barycentric " + "interpolation on the neighbour simplex, which is singular for " + "collinear (2D) or coplanar (3D) neighbours." + ) + + weights = np.zeros((n_targets, nnn), dtype=np.float64) + degenerate = np.zeros(n_targets, dtype=bool) + + size = nnn + dim + 1 + rows_per_chunk = max(1, _CHUNK_BYTES // (size * size * 8)) + + for start in range(0, n_targets, rows_per_chunk): + stop = min(start + rows_per_chunk, n_targets) + _weights_for_chunk( + target_coords[start:stop], + neighbour_coords[start:stop], + weights[start:stop], + degenerate[start:stop], + ) + + return weights, degenerate + + +def _weights_for_chunk(targets, neighbours, weights_out, degenerate_out): + """Fill one chunk of ``weights``/``degenerate`` in place.""" + n_chunk, nnn, dim = neighbours.shape + + # Local frame: target at the origin, farthest neighbour at unit radius. + offsets = neighbours - targets[:, None, :] + radius = np.linalg.norm(offsets, axis=2).max(axis=1) + collapsed = radius <= 0.0 # every neighbour sits on the target + scale = np.where(collapsed, 1.0, radius) + y = offsets / scale[:, None, None] + + # Affine block, and its rank as the primary degeneracy test. + P = np.concatenate([np.ones((n_chunk, nnn, 1)), y], axis=2) + singular_values = np.linalg.svd(P, compute_uv=False) + rank_deficient = singular_values[:, -1] <= _RANK_TOL * singular_values[:, 0] + + healthy = ~(collapsed | rank_deficient) + degenerate_out[...] = ~healthy + if not healthy.any(): + return + + y_h = y[healthy] + P_h = P[healthy] + n_h = y_h.shape[0] + size = nnn + dim + 1 + + pair_distance = np.linalg.norm(y_h[:, :, None, :] - y_h[:, None, :, :], axis=3) + + matrices = np.zeros((n_h, size, size), dtype=np.float64) + matrices[:, :nnn, :nnn] = _polyharmonic(pair_distance) + matrices[:, :nnn, nnn:] = P_h + matrices[:, nnn:, :nnn] = np.transpose(P_h, (0, 2, 1)) + + # RHS: kernel from the target (the origin) to each neighbour, then the + # polynomial basis evaluated at the origin, [1, 0, ..., 0]. + rhs = np.zeros((n_h, size), dtype=np.float64) + rhs[:, :nnn] = _polyharmonic(np.linalg.norm(y_h, axis=2)) + rhs[:, nnn] = 1.0 + + w = _solve_stencils(matrices, rhs)[:, :nnn] + + # Validate the answer rather than guessing a condition number: these two + # identities are what the whole construction is for. + partition_error = np.abs(w.sum(axis=1) - 1.0) + linear_error = np.linalg.norm( + np.einsum("nm,nmd->nd", w, y_h), axis=1 + ) + reproduced = ( + np.isfinite(w).all(axis=1) + & (partition_error <= _REPRODUCTION_TOL) + & (linear_error <= _REPRODUCTION_TOL) + ) + + accepted = np.zeros(n_chunk, dtype=bool) + accepted[healthy] = reproduced + weights_out[accepted] = w[reproduced] + degenerate_out[...] = ~accepted diff --git a/tests/parallel/test_0776_linear_rbf_proxy_parallel.py b/tests/parallel/test_0776_linear_rbf_proxy_parallel.py new file mode 100644 index 000000000..a82430e2f --- /dev/null +++ b/tests/parallel/test_0776_linear_rbf_proxy_parallel.py @@ -0,0 +1,101 @@ +"""Linear-exact swarm proxy transfer under MPI. + +Proxy refresh is strictly rank-local: the kd-tree indexes only this rank's +particles, and there is no halo exchange (SWARM-15, recorded in +``docs/developer/design/SWARM_MODERNIZATION_DESIGN_2026-07.md`` §4). A proxy +node near a partition seam therefore gathers from a one-sided neighbourhood, +which is why proxy values have historically been np-dependent. + +A linear-exact stencil reproduces a linear field exactly from *any* +neighbourhood, one-sided or not. So for a linear field the seam error is +zero, and the proxy becomes np-independent. That is what this test pins. + +It does **not** claim the seam problem is fixed in general: for a field with +curvature the one-sided stencil still differs from a centred one, so np +dependence remains. The test asserts the linear case only, which is the part +that is genuinely exact. + +Run with: + mpirun -n 2 python -m pytest --with-mpi tests/parallel/test_0776_linear_rbf_proxy_parallel.py + mpirun -n 4 python -m pytest --with-mpi tests/parallel/test_0776_linear_rbf_proxy_parallel.py +""" + +import numpy as np +import pytest + +import underworld3 as uw +from underworld3.meshing import UnstructuredSimplexBox + +pytestmark = [ + pytest.mark.level_1, + pytest.mark.tier_b, + pytest.mark.mpi(min_size=2), + pytest.mark.timeout(120), +] + + +def _linear(coords): + return 0.5 + coords @ np.arange(1, coords.shape[1] + 1, dtype=float) + + +def _build(): + mesh = UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=1.0 / 12.0 + ) + swarm = uw.swarm.Swarm(mesh) + var = swarm.add_variable(name="f", size=1, proxy_degree=1) + swarm.populate(fill_param=4) + return mesh, swarm, var + + +@pytest.mark.mpi(min_size=2) +def test_proxy_is_linear_exact_on_every_rank(): + """Every rank's owned proxy nodes carry the analytic linear field.""" + mesh, swarm, var = _build() + + proxy = var._meshVar + var.data[:, 0] = _linear(swarm._particle_coordinates.data) + var._rbf_to_meshVar(proxy) + + node_coords = np.asarray(proxy.coords) + expected = _linear(node_coords) + local_error = ( + np.abs(np.asarray(proxy.data[:, 0]) - expected).max() / np.abs(expected).max() + if node_coords.shape[0] + else 0.0 + ) + + # Collective: a rank-local assertion would let one bad rank pass silently + # while the others reported success (the house failure mode). + global_error = uw.mpi.comm.allreduce(local_error, uw.MPI.MAX) + + assert global_error < 1.0e-12, ( + f"rank-local linear-exact transfer left a global max error of " + f"{global_error:.3e}; seam nodes should still be exact for a linear field" + ) + + del swarm + del mesh + + +@pytest.mark.mpi(min_size=2) +def test_proxy_constant_is_exact_and_collective_refresh_completes(): + """A constant is exact under either scheme; this checks the refresh itself + stays collective (every rank performs the same read-then-write).""" + mesh, swarm, var = _build() + + proxy = var._meshVar + var.data[:, 0] = 2.75 + var._rbf_to_meshVar(proxy) + + local_error = ( + np.abs(np.asarray(proxy.data[:, 0]) - 2.75).max() + if proxy.coords.shape[0] + else 0.0 + ) + global_error = uw.mpi.comm.allreduce(local_error, uw.MPI.MAX) + + assert global_error < 1.0e-12, f"constant field error {global_error:.3e}" + + del swarm + del mesh diff --git a/tests/test_0102_kdtree_linear_rbf.py b/tests/test_0102_kdtree_linear_rbf.py new file mode 100644 index 000000000..2f4ecfc08 --- /dev/null +++ b/tests/test_0102_kdtree_linear_rbf.py @@ -0,0 +1,254 @@ +"""Local RBF stencil weights that reproduce linear fields exactly. + +``KDTree.rbf_interpolator_local(..., order=1)`` builds polyharmonic +(:math:`r^2 \\log r`) weights with an affine tail, so both + + sum_j w_j = 1 (constants exact) + sum_j w_j x_j = x* (linears exact) + +hold by construction. ``order=0`` is inverse-distance weighting, which gets +the first identity but not the second — the property this whole module +exists to add. + +The weights themselves are tested through +``underworld3.utilities.rbf_stencil.linear_exact_weights``, because the +partition-of-unity and sparsity claims are about the *weights* and the +value-returning KDTree API cannot expose them. +""" + +import numpy as np +import pytest + +import underworld3 as uw +from underworld3.utilities.rbf_stencil import linear_exact_weights + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_b] + + +def _stencils(source, target, nnn): + """Brute-force kNN, so the test does not depend on the kd-tree it checks.""" + d = np.linalg.norm(target[:, None, :] - source[None, :, :], axis=2) + idx = np.argsort(d, axis=1)[:, :nnn] + return idx, source[idx] + + +def _linear(coords): + return 0.5 + coords @ np.arange(1, coords.shape[1] + 1, dtype=float) + + +# -------------------------------------------------------------------------- +# The two reproduction identities, at the weights level +# -------------------------------------------------------------------------- +@pytest.mark.parametrize("dim", [2, 3]) +@pytest.mark.parametrize("nnn_rule", ["minimum", "default", "wide"]) +def test_weights_reproduce_constants_and_linears(dim, nnn_rule): + nnn = {"minimum": dim + 2, "default": 2 * (dim + 1), "wide": 14}[nnn_rule] + rng = np.random.default_rng(20260727 + dim) + source = rng.random((500, dim)) + target = rng.random((80, dim)) + + idx, stencil = _stencils(source, target, nnn) + weights, degenerate = linear_exact_weights(target, stencil) + + assert not degenerate.any(), ( + f"{degenerate.sum()} random-cloud stencils were rejected as degenerate" + ) + + partition = np.abs(weights.sum(axis=1) - 1.0).max() + assert partition < 1.0e-12, f"partition of unity violated by {partition:.3e}" + + reproduced = (weights[:, :, None] * stencil).sum(axis=1) + linear = np.abs(reproduced - target).max() + assert linear < 1.0e-12, f"linear reproduction violated by {linear:.3e}" + + values = _linear(source) + interpolated = (weights * values[idx]).sum(axis=1) + error = np.abs(interpolated - _linear(target)).max() + assert error < 1.0e-12, f"linear field interpolated to {error:.3e}, expected round-off" + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_weights_are_sparse(dim): + """Exactly nnn non-zeros per row: this is the point of a *local* scheme.""" + nnn = 2 * (dim + 1) + rng = np.random.default_rng(99) + source = rng.random((300, dim)) + target = rng.random((40, dim)) + + _, stencil = _stencils(source, target, nnn) + weights, _ = linear_exact_weights(target, stencil) + + assert weights.shape == (40, nnn) + assert np.count_nonzero(weights, axis=1).max() <= nnn + + +# -------------------------------------------------------------------------- +# Degeneracy: the affine block loses rank near boundaries and on graded meshes +# -------------------------------------------------------------------------- +def test_collinear_stencil_is_flagged_not_nan(): + stencil = np.zeros((1, 5, 2)) + stencil[0, :, 0] = np.linspace(0.0, 1.0, 5) + stencil[0, :, 1] = 0.3 + target = np.array([[0.5, 0.7]]) + + weights, degenerate = linear_exact_weights(target, stencil) + + assert degenerate[0], "collinear 2D stencil should be flagged degenerate" + assert np.isfinite(weights).all(), "degenerate stencil must not produce NaN" + + +def test_coplanar_stencil_is_flagged_not_nan(): + stencil = np.zeros((1, 6, 3)) + stencil[0, :, 0] = np.linspace(0.0, 1.0, 6) + stencil[0, :, 1] = np.linspace(0.0, 0.5, 6) + stencil[0, :, 2] = 0.2 + target = np.array([[0.4, 0.2, 0.9]]) + + weights, degenerate = linear_exact_weights(target, stencil) + + assert degenerate[0], "coplanar 3D stencil should be flagged degenerate" + assert np.isfinite(weights).all(), "degenerate stencil must not produce NaN" + + +def test_degenerate_and_healthy_stencils_in_one_batch(): + """A singular member must not take the whole batched solve down with it.""" + healthy = np.array([[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]]) + collinear = np.array([[[0.0, 0.0], [0.25, 0.0], [0.5, 0.0], [1.0, 0.0]]]) + stencil = np.concatenate([healthy, collinear]) + target = np.array([[0.3, 0.3], [0.3, 0.3]]) + + weights, degenerate = linear_exact_weights(target, stencil) + + assert not degenerate[0] and degenerate[1] + assert np.isfinite(weights).all() + assert abs(weights[0].sum() - 1.0) < 1.0e-12 + + +def test_stencil_collapsed_onto_target_is_flagged(): + weights, degenerate = linear_exact_weights(np.zeros((1, 2)), np.zeros((1, 4, 2))) + assert degenerate[0] + assert np.isfinite(weights).all() + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_nnn_below_dim_plus_two_is_rejected(dim): + """At dim+1 the affine tail is exactly determined and the scheme reduces to + bare barycentric interpolation on the neighbour simplex — singular whenever + those points are collinear/coplanar, which is common.""" + with pytest.raises(ValueError, match=f"dim \\+ 2 = {dim + 2}"): + linear_exact_weights(np.zeros((1, dim)), np.zeros((1, dim + 1, dim))) + + +# -------------------------------------------------------------------------- +# Through the KDTree API +# -------------------------------------------------------------------------- +@pytest.mark.parametrize("dim", [2, 3]) +def test_kdtree_order1_beats_order0_on_a_linear_field(dim): + nnn = 2 * (dim + 1) + rng = np.random.default_rng(1234 + dim) + source = rng.random((2000, dim)) + target = rng.random((150, dim)) + data = _linear(source)[:, None] + expected = _linear(target) + + kdt = uw.kdtree.KDTree(source) + shepard = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=0) + exact = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=1) + + exact_error = np.abs(exact[:, 0] - expected).max() + shepard_error = np.abs(shepard[:, 0] - expected).max() + + assert exact_error < 1.0e-12, f"order=1 gave {exact_error:.3e} on a linear field" + assert shepard_error > 1.0e-4, ( + "order=0 is expected to smear a linear field; if this fails the " + f"baseline has changed (got {shepard_error:.3e})" + ) + + +def test_kdtree_order0_is_unchanged_by_the_new_arguments(): + """Back-compatibility: the positional call must be bit-identical.""" + rng = np.random.default_rng(5) + source = rng.random((400, 2)) + target = rng.random((50, 2)) + data = _linear(source)[:, None] + + kdt = uw.kdtree.KDTree(source) + positional = kdt.rbf_interpolator_local(target, data, 4, 2, False) + explicit = kdt.rbf_interpolator_local(target, data, 4, 2, False, order=0, + monotone=False) + + assert np.array_equal(positional, explicit) + + +def test_kdtree_monotone_clamps_to_the_stencil_range(): + """order=1 weights are signed and can overshoot; monotone bounds them.""" + dim = 2 + nnn = 2 * (dim + 1) + rng = np.random.default_rng(77) + source = rng.random((600, dim)) + target = rng.random((200, dim)) + data = _linear(source)[:, None] + + kdt = uw.kdtree.KDTree(source) + unbounded = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=1) + clamped = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=1, + monotone=True) + + lo, hi = data.min(), data.max() + assert clamped.min() >= lo - 1.0e-12 and clamped.max() <= hi + 1.0e-12 + # Sanity: the clamp is doing something on this configuration, otherwise the + # test would pass vacuously. + assert np.abs(clamped - unbounded).max() > 0.0 + + +def test_kdtree_rejects_bad_order_and_monotone_mode(): + rng = np.random.default_rng(11) + source = rng.random((200, 2)) + target = rng.random((10, 2)) + data = _linear(source)[:, None] + kdt = uw.kdtree.KDTree(source) + + with pytest.raises(ValueError, match="order must be 0"): + kdt.rbf_interpolator_local(target, data, 6, 2, False, order=2) + + with pytest.raises(ValueError, match="no meaning for a local kd-tree"): + kdt.rbf_interpolator_local(target, data, 6, 2, False, order=1, + monotone="pick") + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_error_falls_with_point_spacing_faster_than_inverse_distance(dim): + """On a quadratic field, order=1 should converge under refinement and + order=0 should barely move — inverse-distance error is set by the stencil + geometry, not by how close the points are.""" + + def quadratic(coords): + return 0.5 + (coords ** 2).sum(axis=1) + + nnn = 2 * (dim + 1) + errors = {} + for n_per_side in (10, 20): + axis = (np.arange(n_per_side) + 0.5) / n_per_side + grid = np.meshgrid(*([axis] * dim), indexing="ij") + source = np.stack([g.ravel() for g in grid], axis=1) + rng = np.random.default_rng(3) + target = 0.2 + 0.6 * rng.random((100, dim)) + data = quadratic(source)[:, None] + expected = quadratic(target) + + kdt = uw.kdtree.KDTree(source) + for order in (0, 1): + got = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=order) + errors[(order, n_per_side)] = np.abs(got[:, 0] - expected).max() + + order1_ratio = errors[(1, 10)] / errors[(1, 20)] + order0_ratio = errors[(0, 10)] / errors[(0, 20)] + + assert order1_ratio > 3.0, ( + f"order=1 should converge at least ~O(h^2) under halving the spacing, " + f"got a factor of {order1_ratio:.2f}" + ) + assert order1_ratio > order0_ratio, ( + f"order=1 convergence ({order1_ratio:.2f}x) should beat order=0 " + f"({order0_ratio:.2f}x)" + ) diff --git a/tests/test_0506_swarm_proxy_linear_exact.py b/tests/test_0506_swarm_proxy_linear_exact.py new file mode 100644 index 000000000..e2d6196bf --- /dev/null +++ b/tests/test_0506_swarm_proxy_linear_exact.py @@ -0,0 +1,120 @@ +"""Swarm proxy variables transfer a linear field without smearing it. + +A proxy mesh variable is refreshed from the particles by local RBF +interpolation. That transfer used to be inverse-distance (Shepard) weighting, +which reproduces constants but not gradients; the proxy default is now +``order=1``, which reproduces both. + +Why a *linear* field is the right probe: it lies exactly inside both the P1 +and the P2 proxy space, so the finite element discretisation contributes +nothing at all to the measured error. Everything left is particle -> node +transfer error, which makes the assertion below unambiguous. + +Measured on this configuration when the default was changed (2026-07-27): +inverse distance gave a relative max error of ~5e-3 on the linear field and +~9e-3 on a quadratic one; the linear-exact scheme gives round-off and ~1e-4 +respectively. +""" + +import numpy as np +import pytest + +import underworld3 as uw +from underworld3.meshing import UnstructuredSimplexBox + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_b] + + +def _linear(coords): + return 0.5 + coords @ np.arange(1, coords.shape[1] + 1, dtype=float) + + +def _quadratic(coords): + return 0.5 + (coords ** 2).sum(axis=1) + + +@pytest.fixture +def proxied_swarm(): + mesh = UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=1.0 / 8.0 + ) + swarm = uw.swarm.Swarm(mesh) + var = swarm.add_variable(name="f", size=1, proxy_degree=1) + swarm.populate(fill_param=4) + yield swarm, var + del swarm + del mesh + + +def _proxy_error(var, field, **refresh): + """Refresh the proxy from an analytic field and return its relative error.""" + proxy = var._meshVar + particle_coords = var.swarm._particle_coordinates.data.copy() + node_coords = np.asarray(proxy.coords) + + var.data[:, 0] = field(particle_coords) + var._rbf_to_meshVar(proxy, **refresh) + + expected = field(node_coords) + return np.abs(np.asarray(proxy.data[:, 0]) - expected).max() / np.abs(expected).max() + + +def test_proxy_default_reproduces_a_linear_field(proxied_swarm): + """The shipped default must be linear-exact — no arguments passed.""" + _, var = proxied_swarm + error = _proxy_error(var, _linear) + assert error < 1.0e-12, ( + f"proxy of an exactly linear field carries a relative error of " + f"{error:.3e}; the default transfer is not linear-exact" + ) + + +def test_inverse_distance_smears_the_same_linear_field(proxied_swarm): + """The control arm: without order=1 the error is real and much larger. + + This is what stops the test above passing for the wrong reason (e.g. a + tolerance that any scheme would meet on this mesh). + """ + _, var = proxied_swarm + exact = _proxy_error(var, _linear, order=1) + shepard = _proxy_error(var, _linear, order=0, nnn=3) + + assert shepard > 1.0e-4, ( + f"inverse distance is expected to smear a linear field, got {shepard:.3e}" + ) + assert shepard > 1.0e6 * max(exact, 1.0e-16) + + +def test_proxy_default_improves_a_quadratic_field(proxied_swarm): + """Linear exactness is not a trick that only helps linear fields.""" + _, var = proxied_swarm + exact = _proxy_error(var, _quadratic, order=1) + shepard = _proxy_error(var, _quadratic, order=0, nnn=3) + + assert exact < shepard / 10.0, ( + f"order=1 gave {exact:.3e} vs inverse distance {shepard:.3e} on a " + "quadratic field; expected at least a 10x improvement" + ) + + +def test_proxy_still_reproduces_a_constant_exactly(proxied_swarm): + """Constants were already exact under inverse distance — do not regress.""" + _, var = proxied_swarm + error = _proxy_error(var, lambda coords: np.full(coords.shape[0], 3.25)) + assert error < 1.0e-12, f"constant field carries error {error:.3e}" + + +def test_proxy_monotone_bounds_values_to_the_particle_range(proxied_swarm): + """order=1 weights are signed, so bounded quantities can opt back in.""" + _, var = proxied_swarm + proxy = var._meshVar + particle_coords = var.swarm._particle_coordinates.data.copy() + + var.data[:, 0] = _linear(particle_coords) + lo, hi = var.data[:, 0].min(), var.data[:, 0].max() + + var._rbf_to_meshVar(proxy, monotone=True) + values = np.asarray(proxy.data[:, 0]) + + assert values.min() >= lo - 1.0e-12 + assert values.max() <= hi + 1.0e-12 From bedfef261b5d60767f264b209748f3dac77a484b Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 27 Jul 2026 10:48:59 +1000 Subject: [PATCH 2/5] Account for the pairwise-offset temporary in chunk sizing; skip the evaluator import when no limiter is asked for The chunk budget counted only the stacked saddle matrices and missed the (nnn, nnn, dim) offset array the kernel distances are formed from, which is the larger of the two in 3D. _normalise_monotone now returns before importing from function/, so the default path does not couple a low-level module to the evaluator. Underworld development team with AI support from Claude Code --- src/underworld3/ckdtree.pyx | 6 ++++++ src/underworld3/utilities/rbf_stencil.py | 7 ++++++- 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/src/underworld3/ckdtree.pyx b/src/underworld3/ckdtree.pyx index 733d30bc1..800c5dcc4 100644 --- a/src/underworld3/ckdtree.pyx +++ b/src/underworld3/ckdtree.pyx @@ -44,6 +44,12 @@ def _normalise_monotone(monotone): out-of-bounds points through the FE path, which a kd-tree knows nothing about. """ + # Early out before the import: the overwhelmingly common call has no + # limiter, and ckdtree is a low-level module that should not take a + # dependency on the evaluator just to be told "no". + if monotone is False or monotone is None: + return False + from underworld3.function.functions_unit_system import _normalize_monotone mode = _normalize_monotone(monotone) diff --git a/src/underworld3/utilities/rbf_stencil.py b/src/underworld3/utilities/rbf_stencil.py index 1eaa00974..1d7b0dc97 100644 --- a/src/underworld3/utilities/rbf_stencil.py +++ b/src/underworld3/utilities/rbf_stencil.py @@ -127,8 +127,13 @@ def linear_exact_weights(target_coords, neighbour_coords): weights = np.zeros((n_targets, nnn), dtype=np.float64) degenerate = np.zeros(n_targets, dtype=bool) + # Budget the two largest per-row temporaries: the stacked saddle matrices + # (size x size) and the pairwise offset array (nnn x nnn x dim) that the + # kernel distances are formed from. The latter is easy to overlook and is + # the bigger of the two in 3D. size = nnn + dim + 1 - rows_per_chunk = max(1, _CHUNK_BYTES // (size * size * 8)) + bytes_per_row = 8 * (size * size + nnn * nnn * dim) + rows_per_chunk = max(1, _CHUNK_BYTES // bytes_per_row) for start in range(0, n_targets, rows_per_chunk): stop = min(start + rows_per_chunk, n_targets) From d395424584435a8996b192d87fc35f15a53bee23 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 27 Jul 2026 11:07:30 +1000 Subject: [PATCH 3/5] Limit the RBF correction, not the linear reconstruction The monotone limiter clipped the result to the min/max of its stencil's source values. That is wrong for a linear-exact scheme, and wrong in a way that is easy to miss: a target outside the convex hull of its own nnn neighbours has a value outside their range even for an exactly linear field. 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 -- measured, it returned a linear field to ~4e-3, i.e. inverse-distance level. A limiter that destroys the guarantee the scheme exists to provide is a defect, not a documented trade-off. Replaced with the slope-limiter discipline: never limit the linear reconstruction, limit only the correction on top of it. 1. least-squares affine fit over the stencil (utilities/rbf_stencil.py, affine_trend) 2. keep that trend at the target untouched 3. bound the remaining RBF correction to the range of 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, was 4e-3). It still bites on rough data -- on sin(6x) + |x|^2 it moved the answer by 2e-3 in 2D and 5e-2 in 3D, and in 3D that costs accuracy (9.0e-2 -> 1.3e-1, still well below inverse distance at 2.7e-1), which is the ordinary limiter trade. What it promises is now stated precisely: it bounds new oscillation relative to the local trend, not absolute range. A quantity with hard physical bounds still needs its own clip. The two tests that asserted the old, wrong guarantee are replaced by ones asserting the new one, including that the limiter leaves a linear field alone. Also reframes the motivation throughout: this is general interpolation accuracy. Whether it helps the free-surface topography noise is an open question, not the reason for the work. Underworld development team with AI support from Claude Code --- docs/developer/subsystems/interpolation.md | 67 ++++++++++++++------ src/underworld3/ckdtree.pyx | 41 ++++++++++--- src/underworld3/swarm.py | 6 +- src/underworld3/utilities/rbf_stencil.py | 44 +++++++++++++ tests/test_0102_kdtree_linear_rbf.py | 68 +++++++++++++++++---- tests/test_0506_swarm_proxy_linear_exact.py | 26 ++++---- 6 files changed, 196 insertions(+), 56 deletions(-) diff --git a/docs/developer/subsystems/interpolation.md b/docs/developer/subsystems/interpolation.md index 4d40e535d..e9ee17d6d 100644 --- a/docs/developer/subsystems/interpolation.md +++ b/docs/developer/subsystems/interpolation.md @@ -10,8 +10,13 @@ points, one mesh to another after adaptation. The tool for that is a *local* interpolant — for each target point, take its `nnn` nearest source points from a kd-tree and form a weighted average. -The only question that matters is **which fields the weights reproduce -exactly**. That is what `order` selects on +This is general-purpose numerical machinery, and the goal is straightforward: +**the interpolation should be accurate**. The single question that determines +that is **which fields the weights reproduce exactly**. A scheme that cannot +reproduce a uniform gradient will smear every field that has one, everywhere it +is used, and no amount of tuning elsewhere recovers what it threw away. + +That is what `order` selects on {py:meth}`underworld3.ckdtree.KDTree.rbf_interpolator_local`. ## The two schemes @@ -80,9 +85,9 @@ of the measured error is particle-to-node transfer error. | quadratic, rms | 2.0e-3 | 2.3e-5 | Two things to read off. First, the linear-field column is the clean statement: -one scheme is exact, the other is not. Second — and this is the part that -actually matters for physics — the improvement carries over to a field with -curvature, by roughly two orders of magnitude in rms. +one scheme is exact, the other is not. Second — and this is the part that makes +it a general accuracy improvement rather than a special case — the gain carries +over to a field with curvature, by roughly two orders of magnitude in rms. Widening the `order=0` stencil does **not** close the gap. At `nnn=12` the inverse-distance error is unchanged to two significant figures, because the @@ -121,26 +126,48 @@ The implementation (`underworld3.utilities.rbf_stencil.linear_exact_weights`): geometric fallback is the failure mode behind issue #424; a per-point warning would be unusable, a single counted one is not. -## Boundedness: `monotone` +## Limiting: `monotone` + +`order=1` weights are signed, so the interpolant can oscillate where the data +is rough. `monotone=True` (or `"clamp"`) bounds that. -`monotone=True` (or `"clamp"`, matching the vocabulary of -{py:func}`underworld3.function.evaluate`) bounds each result to the min/max of -its own stencil's source values, restoring the boundedness that `order=1` gives -up. +The important part is *what* it bounds. The obvious limiter — clip the result +to the min/max of the stencil's source values — is **wrong for a linear-exact +scheme**, and wrong in a way that is easy to miss: ```{warning} -Clamping and linear exactness are in genuine conflict. Where the target lies -outside the convex hull of its neighbours the exact linear value is *outside* -the stencil range, so the clamp discards it. - -Measured on the 2D case above, `order=1` with `monotone=True` gives a max error -of 4.0e-3 on a linear field — back to the inverse-distance level — while -keeping the interior nodes exact. Use it when a consumer genuinely requires -bounded values, not as a default safety net. +A target that lies outside the convex hull of its own `nnn` neighbours has a +value outside their range **even for an exactly linear field**. Proxy nodes are +routinely in that position, and not only at domain boundaries. A raw min/max +clip therefore cannot distinguish legitimate extrapolation from ringing, and +fires on the linear part — destroying the one guarantee the scheme exists to +provide. ``` -Consumers that do depend on boundedness, and are therefore deliberately left on -`order=0`: +So the limiter here follows the slope-limiter discipline instead: **never limit +the linear reconstruction, limit only the correction on top of it.** + +1. Fit an affine function to the stencil data by least squares. +2. Keep that trend at the target untouched. +3. Bound the remaining RBF correction to the range of the non-affine residual + the stencil actually exhibits. + +Consequences, both measured: + +- On any field the scheme already reproduces exactly, the limiter is a **no-op** + (it moves the answer by ~1e-15). Linear reproduction survives it. +- On a rough field it does bite — and, like any limiter, it trades accuracy for + boundedness where the correction genuinely exceeds the observed residual. On + a `sin(6x) + |x|²` test field it moved the answer by 2e-3 (2D) and 5e-2 (3D); + in 3D that made the max error slightly worse (9.0e-2 → 1.3e-1), still well + below inverse distance (2.7e-1). + +Note what this does and does not promise: it bounds *new oscillation relative +to the local trend*, not absolute range. A quantity that must stay inside hard +physical bounds (a fraction in $[0,1]$) needs its own clip on top. + +Consumers that depend on absolute boundedness, and are therefore deliberately +left on `order=0`: - The **RBF rung of the point-location fallback ladder** in `uw.function.evaluate` — `docs/developer/design/point-location-capability.md` diff --git a/src/underworld3/ckdtree.pyx b/src/underworld3/ckdtree.pyx index 800c5dcc4..c65ff6a77 100644 --- a/src/underworld3/ckdtree.pyx +++ b/src/underworld3/ckdtree.pyx @@ -480,11 +480,14 @@ cdef class KDTree: order : int, optional Polynomial reproduction order, 0 (default) or 1. monotone : bool or str, optional - Bound each interpolated value to the min/max of its own stencil's - source values. ``False`` (default) or ``True`` / ``"clamp"``. - Restores boundedness at ``order=1``, but note that clamping - necessarily discards linear exactness wherever the target lies - outside the convex hull of its neighbours (boundary points). + ``False`` (default) or ``True`` / ``"clamp"``. Limits the + **non-affine part** of the interpolant: the local least-squares + affine trend is preserved exactly, and only the RBF correction on + top of it is bounded by the correction actually present in the + stencil. This is the slope-limiter discipline — the linear + reconstruction is never clipped — so linear reproduction survives + the limiter and it is a no-op on any field the scheme already + reproduces exactly. It bounds new oscillation, not absolute range. Returns ------- @@ -754,9 +757,31 @@ cdef class KDTree: # print(valz) if monotone_clamp: - # Bound to the stencil's own range. This is what restores the - # convex-combination property that order=1 gives up. - vals = np.clip(vals, kdata.min(axis=1), kdata.max(axis=1)) + if order == 0: + # Already a convex combination; the clip is exact-arithmetic + # redundant and only guards round-off. + vals = np.clip(vals, kdata.min(axis=1), kdata.max(axis=1)) + else: + # Limit the CORRECTION, never the linear part -- the same + # discipline as a slope limiter in a second-order FV scheme. + # + # Clipping the total against the stencil's raw min/max would + # be wrong here: a target outside the convex hull of its own + # neighbours has a value outside their range even for an + # exactly linear field, so the naive clip cannot distinguish + # legitimate extrapolation from ringing and destroys the + # reproduction guarantee. + from underworld3.utilities.rbf_stencil import affine_trend + + trend_at_target, trend_at_stencil = affine_trend( + coords_converted, stencil_coords, kdata + ) + residual = kdata - trend_at_stencil + vals = trend_at_target + np.clip( + vals - trend_at_target, + residual.min(axis=1), + residual.max(axis=1), + ) if verbose and uw.mpi.rank == 0: print(f"Mapping values ... finished", flush=True) diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py index 65a67840d..4b4a3e53c 100644 --- a/src/underworld3/swarm.py +++ b/src/underworld3/swarm.py @@ -1516,9 +1516,9 @@ def rbf_interpolate(self, new_coords, verbose=False, nnn=None, order=1, 0 (constants only, inverse distance). Drops to 0 automatically on a rank holding too few particles to determine the affine tail. monotone : bool or str, default=False - Bound each value to the min/max of its own particle stencil. Note - this discards linear exactness wherever the target sits outside the - convex hull of its neighbours. + Limit the non-affine part of the interpolant to the non-affine + variation present in the particle stencil. The local linear trend + is preserved, so this does not cost linear exactness. Returns ------- diff --git a/src/underworld3/utilities/rbf_stencil.py b/src/underworld3/utilities/rbf_stencil.py index 1d7b0dc97..14895c46a 100644 --- a/src/underworld3/utilities/rbf_stencil.py +++ b/src/underworld3/utilities/rbf_stencil.py @@ -49,6 +49,50 @@ _CHUNK_BYTES = 64 * 1024 * 1024 +def affine_trend(target_coords, neighbour_coords, values): + """Least-squares affine fit of ``values`` over each stencil. + + Returns the fit evaluated at the target and at each neighbour, which is + what a limiter needs in order to bound the *non-affine* part of an + interpolant without touching its linear part. + + Parameters + ---------- + target_coords : numpy.ndarray + Shape ``(n_targets, dim)``. + neighbour_coords : numpy.ndarray + Shape ``(n_targets, nnn, dim)``. + values : numpy.ndarray + Stencil data, shape ``(n_targets, nnn, n_components)``. + + Returns + ------- + at_target : numpy.ndarray + Shape ``(n_targets, n_components)``. + at_neighbours : numpy.ndarray + Shape ``(n_targets, nnn, n_components)``. + """ + target_coords = np.ascontiguousarray(target_coords, dtype=np.float64) + neighbour_coords = np.ascontiguousarray(neighbour_coords, dtype=np.float64) + + n_targets, nnn, _ = neighbour_coords.shape + + # Local frame again, so the constant column and the coordinate columns are + # comparably scaled and the fit is not dominated by the offset. + offsets = neighbour_coords - target_coords[:, None, :] + radius = np.linalg.norm(offsets, axis=2).max(axis=1) + y = offsets / np.where(radius > 0.0, radius, 1.0)[:, None, None] + + P = np.concatenate([np.ones((n_targets, nnn, 1)), y], axis=2) + # pinv rather than solve: a rank-deficient stencil gives the minimum-norm + # fit instead of raising, and those points are limited conservatively. + coefficients = np.linalg.pinv(P) @ values + + # The target sits at the origin of the local frame, so its affine value is + # the constant coefficient. + return coefficients[:, 0, :], P @ coefficients + + def _polyharmonic(r): """:math:`\\varphi(r) = r^2 \\log r`, with :math:`\\varphi(0) = 0`.""" # The clip only keeps log() finite at coincident points; r**2 drives the diff --git a/tests/test_0102_kdtree_linear_rbf.py b/tests/test_0102_kdtree_linear_rbf.py index 2f4ecfc08..619fcfbd4 100644 --- a/tests/test_0102_kdtree_linear_rbf.py +++ b/tests/test_0102_kdtree_linear_rbf.py @@ -180,25 +180,69 @@ def test_kdtree_order0_is_unchanged_by_the_new_arguments(): assert np.array_equal(positional, explicit) -def test_kdtree_monotone_clamps_to_the_stencil_range(): - """order=1 weights are signed and can overshoot; monotone bounds them.""" +@pytest.mark.parametrize("dim", [2, 3]) +def test_monotone_does_not_touch_a_linear_field(dim): + """The limiter must not cost the guarantee the scheme exists to provide. + + A naive clamp against the stencil's raw min/max does exactly that: a + target outside the convex hull of its own neighbours has a value outside + their range even for an exactly linear field, so such a clamp cannot tell + legitimate extrapolation from ringing. The limiter therefore bounds the + non-affine part only, leaving the linear reconstruction alone. + """ + nnn = 2 * (dim + 1) + rng = np.random.default_rng(77 + dim) + source = rng.random((1500, dim)) + target = rng.random((250, dim)) + data = _linear(source)[:, None] + expected = _linear(target) + + kdt = uw.kdtree.KDTree(source) + unlimited = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=1) + limited = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=1, + monotone=True) + + assert np.abs(limited[:, 0] - expected).max() < 1.0e-12, ( + "the limiter destroyed linear reproduction" + ) + assert np.abs(limited - unlimited).max() < 1.0e-12, ( + "the limiter should be a no-op on a field the scheme reproduces exactly" + ) + + +def test_monotone_bounds_the_correction_on_a_curved_field(): + """It must still do something: bound the RBF correction to the non-affine + variation actually present in the stencil.""" dim = 2 nnn = 2 * (dim + 1) - rng = np.random.default_rng(77) - source = rng.random((600, dim)) + rng = np.random.default_rng(4242) + source = rng.random((800, dim)) target = rng.random((200, dim)) - data = _linear(source)[:, None] + values = 0.5 + (source ** 2).sum(axis=1) + np.sin(6.0 * source[:, 0]) + data = values[:, None] kdt = uw.kdtree.KDTree(source) - unbounded = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=1) - clamped = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=1, + unlimited = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=1) + limited = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=1, monotone=True) - lo, hi = data.min(), data.max() - assert clamped.min() >= lo - 1.0e-12 and clamped.max() <= hi + 1.0e-12 - # Sanity: the clamp is doing something on this configuration, otherwise the - # test would pass vacuously. - assert np.abs(clamped - unbounded).max() > 0.0 + assert np.abs(limited - unlimited).max() > 1.0e-6, ( + "the limiter had no effect on a curved field, so it is not limiting" + ) + + # The stated guarantee: the deviation from the local affine trend never + # exceeds the deviation the stencil itself shows. + from underworld3.utilities.rbf_stencil import affine_trend + + _, stencil = _stencils(source, target, nnn) + idx, _ = _stencils(source, target, nnn) + stencil_values = data[idx] + at_target, at_stencil = affine_trend(target, stencil, stencil_values) + residual = stencil_values - at_stencil + correction = limited - at_target + + assert (correction <= residual.max(axis=1) + 1.0e-12).all() + assert (correction >= residual.min(axis=1) - 1.0e-12).all() def test_kdtree_rejects_bad_order_and_monotone_mode(): diff --git a/tests/test_0506_swarm_proxy_linear_exact.py b/tests/test_0506_swarm_proxy_linear_exact.py index e2d6196bf..b2ad40bb1 100644 --- a/tests/test_0506_swarm_proxy_linear_exact.py +++ b/tests/test_0506_swarm_proxy_linear_exact.py @@ -104,17 +104,17 @@ def test_proxy_still_reproduces_a_constant_exactly(proxied_swarm): assert error < 1.0e-12, f"constant field carries error {error:.3e}" -def test_proxy_monotone_bounds_values_to_the_particle_range(proxied_swarm): - """order=1 weights are signed, so bounded quantities can opt back in.""" - _, var = proxied_swarm - proxy = var._meshVar - particle_coords = var.swarm._particle_coordinates.data.copy() - - var.data[:, 0] = _linear(particle_coords) - lo, hi = var.data[:, 0].min(), var.data[:, 0].max() +def test_proxy_monotone_keeps_the_linear_field_exact(proxied_swarm): + """Turning the limiter on must not cost the proxy its exactness. - var._rbf_to_meshVar(proxy, monotone=True) - values = np.asarray(proxy.data[:, 0]) - - assert values.min() >= lo - 1.0e-12 - assert values.max() <= hi + 1.0e-12 + Proxy nodes routinely sit outside the convex hull of their own particle + stencil — not only at the domain boundary — so a limiter that clipped + against the stencil's raw min/max would fire on the linear part. This one + limits only the non-affine correction. + """ + _, var = proxied_swarm + limited = _proxy_error(var, _linear, monotone=True) + assert limited < 1.0e-12, ( + f"proxy with the limiter on carries {limited:.3e} on an exactly linear " + "field; the limiter is clipping the linear reconstruction" + ) From 7b6e3981af459713881256b7a2b3c33dca7e9c28 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 27 Jul 2026 14:23:33 +1000 Subject: [PATCH 4/5] Widen the stencil before surrendering a degenerate point; settle IndexSwarmVariable Two gaps closed after scale testing. 1. Degenerate stencils are now RETRIED on a wider neighbourhood (4x nnn, capped at the cloud size) before falling back to inverse distance. Found by running the new default at scale: a 3D low-density case (190k particles, 29k proxy nodes, fill_param=2) had 2 stencils in 28824 fail the affine-rank test, and those two nodes carried the entire max error -- 2.2e-3 while every other node was at 1e-15. That is the damaging shape of error: an isolated non-exact node sitting among exact ones is a SPIKE, which corrupts derivatives and integrals far more than a smooth error of the same size. A degenerate neighbourhood is usually a local accident of point placement, so reaching further escapes it. Only the failed points are re-queried, so cost tracks how many actually failed. 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 cannot be rescued by widening, and there the honest outcome is unchanged: bounded inverse-distance values plus a loud warning. Both paths now have tests. 2. IndexSwarmVariable stays on inverse distance, measured rather than assumed. A material indicator is piecewise constant: away from an interface both schemes reproduce it exactly, and at the interface it is discontinuous so no polynomial reproduction is exact either. Measured on a straight interface at x=0.5, the interface-position error is a wash (order=1 better at coarse resolution, equal or worse at fine, ordering flipping between cases) while nnn=8 pushed level sets to [-0.038, 1.038]. Partition of unity survives either way, 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, the material jump does not. Also notes the libcpp `bool` cimport trap in ckdtree.pyx, which makes `dtype=bool` a compile error in that module. Regression: level_1 496 passed, level_2 312 passed, np=2 parallel 6 passed, np=1/np=4 proxy agreement to round-off at 190k particles. Underworld development team with AI support from Claude Code --- docs/developer/subsystems/interpolation.md | 43 +++++++++++ src/underworld3/ckdtree.pyx | 88 ++++++++++++++++++---- tests/test_0102_kdtree_linear_rbf.py | 73 ++++++++++++++++++ 3 files changed, 188 insertions(+), 16 deletions(-) diff --git a/docs/developer/subsystems/interpolation.md b/docs/developer/subsystems/interpolation.md index e9ee17d6d..a3a912c19 100644 --- a/docs/developer/subsystems/interpolation.md +++ b/docs/developer/subsystems/interpolation.md @@ -166,6 +166,49 @@ Note what this does and does not promise: it bounds *new oscillation relative to the local trend*, not absolute range. A quantity that must stay inside hard physical bounds (a fraction in $[0,1]$) needs its own clip on top. +### Material level sets stay on `order=0` — measured, not assumed + +`IndexSwarmVariable` builds one level-set MeshVariable per material index and +keeps its own inverse-distance weighting. It was tested against `order=1` and +**deliberately not changed**. + +The reason is structural. A material indicator is **piecewise constant**, not +smooth. Away from an interface both schemes reproduce it exactly, because both +reproduce constants — linear exactness has nothing to add. At the interface the +field is *discontinuous*, so no polynomial-reproducing scheme is exact either; +signed weights simply add overshoot where the data has a jump. + +Measured on a straight interface at `x = 0.5` (exactly representable, so any +displacement of the recovered 0.5 contour is scheme error): + +| 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]`** | + +The accuracy result is a wash — `order=1` is better at the coarse resolution +and equal or worse at the fine one, with the ordering flipping between cases — +while `nnn=8` violates the `[0, 1]` bound by ~3.8%. + +Partition of unity survives either way (all indices share one weight set, and +the indicator flags sum to one per particle, so the level sets sum to +`Σ w_j = 1` regardless of sign). But a *negative* material fraction is still +physically wrong, and `constitutive_models.py` consumes these directly. + +```{note} +The interface metric groups nodes into rows by `y` and interpolates the 0.5 +crossing, which is crude on an unstructured simplex mesh — the *maximum* error +is identical across all schemes because it is set by node spacing, not by the +weights. Only the median is informative, and it is the median that shows no +consistent gain. +``` + +So the swarm story is deliberately split: the plain `SwarmVariable` proxy takes +`order=1` because its fields are smooth and the gain is two orders of magnitude; +`IndexSwarmVariable` keeps inverse distance because its field is a jump, where +there is nothing to gain and a bound to lose. + Consumers that depend on absolute boundedness, and are therefore deliberately left on `order=0`: diff --git a/src/underworld3/ckdtree.pyx b/src/underworld3/ckdtree.pyx index c65ff6a77..8abc9b931 100644 --- a/src/underworld3/ckdtree.pyx +++ b/src/underworld3/ckdtree.pyx @@ -632,6 +632,47 @@ cdef class KDTree: return Values + def _retry_degenerate_stencils(self, coords_converted, data, tree_points, + vals, degenerate, nnn): + """Re-solve the failed points on a wider stencil. + + A neighbourhood that cannot support an affine fit is usually a local + accident of point placement, so reaching further normally escapes it. + Only the failed points are re-queried, so the cost is set by how many + actually failed (measured: a handful in tens of thousands). + + Returns the updated values and the mask of points that failed even on + the wider stencil, which keep their inverse-distance values. + """ + from underworld3.utilities.rbf_stencil import linear_exact_weights + + wide = min(4 * nnn, self.n) + if wide <= nnn: + return vals, degenerate + + failed = np.flatnonzero(degenerate) + targets = np.ascontiguousarray(coords_converted[failed], dtype=np.float64) + + # The raw query: coords_converted is already in tree units, so the + # unit-aware `query` would reject it when the tree carries units. + indices, _ = self.find_closest_n_points(wide, targets) + if np.any(indices >= self.n): + return vals, degenerate + + weights, still_degenerate = linear_exact_weights( + targets, tree_points[indices] + ) + + recovered = ~still_degenerate + if recovered.any(): + vals[failed[recovered]] = np.einsum( + "sdc,sd->sc", data[indices[recovered]], weights[recovered] + ) + + remaining = degenerate.copy() + remaining[failed] = still_degenerate + return vals, remaining + def rbf_interpolator_local_from_kdtree(self, coords, data, nnn, p, verbose, order=0, monotone=False): """ @@ -707,6 +748,10 @@ cdef class KDTree: # print(f"kd-tree diagnostics: d.shape - {distance_n.shape}, c.shape - {closest_n.shape}") print(f"Mapping values with nnn - {nnn} & p {p} ... start", flush=True) + # np.bool_, not bool: this module cimports the C++ `bool` from libcpp, + # which shadows the Python builtin and will not compile here. + degenerate = np.zeros(coords_converted.shape[0], dtype=np.bool_) + if nnn == 1: # only use nearest neighbour raw data if order == 1: @@ -729,33 +774,44 @@ cdef class KDTree: if order == 1: from underworld3.utilities.rbf_stencil import linear_exact_weights - stencil_coords = np.asarray(self.points)[closest_n[:]] + tree_points = np.asarray(self.points) + stencil_coords = tree_points[closest_n[:]] linear_weights, degenerate = linear_exact_weights( coords_converted, stencil_coords ) # A stencil that cannot support an affine fit (collinear in 2D, - # coplanar in 3D) keeps the inverse-distance weights: less accurate - # there, but finite and bounded. Report it once — a silent - # geometric fallback is the failure mode of issue #424. + # coplanar in 3D) keeps the inverse-distance weights for now: finite + # and bounded, but not linear-exact. linear_weights[degenerate] = n_weights[degenerate] n_weights = linear_weights - n_degenerate = int(degenerate.sum()) - if n_degenerate: - import warnings - - warnings.warn( - f"rbf_interpolator_local(order=1): {n_degenerate} of " - f"{degenerate.size} stencils could not support an affine fit " - "(collinear/coplanar neighbours) and fell back to inverse-" - "distance weighting. Increasing nnn usually removes them.", - stacklevel=2, - ) - # magic with einstein summation power vals = np.einsum("sdc,sd->sc", kdata, n_weights) # print(valz) + if order == 1 and degenerate.any(): # noqa: F821 - set in the order==1 block + # Widen the stencil for the few points that failed, rather than + # accepting inverse distance there. A degenerate neighbourhood is + # usually a local accident of point placement, and reaching further + # fixes it -- worth doing because a single non-exact node sitting + # among exact ones is an isolated SPIKE, which is far more damaging + # to derivatives and integrals than a smooth error of the same size. + vals, degenerate = self._retry_degenerate_stencils( + coords_converted, data, tree_points, vals, degenerate, nnn + ) + + n_degenerate = int(degenerate.sum()) + if n_degenerate: + import warnings + + warnings.warn( + f"rbf_interpolator_local(order=1): {n_degenerate} of " + f"{degenerate.size} stencils could not support an affine fit " + "(collinear/coplanar neighbours) even after widening, and fell " + "back to inverse-distance weighting.", + stacklevel=2, + ) + if monotone_clamp: if order == 0: # Already a convex combination; the clip is exact-arithmetic diff --git a/tests/test_0102_kdtree_linear_rbf.py b/tests/test_0102_kdtree_linear_rbf.py index 619fcfbd4..19f6b88e4 100644 --- a/tests/test_0102_kdtree_linear_rbf.py +++ b/tests/test_0102_kdtree_linear_rbf.py @@ -130,6 +130,79 @@ def test_stencil_collapsed_onto_target_is_flagged(): assert np.isfinite(weights).all() +def test_degenerate_stencil_is_retried_on_a_wider_neighbourhood(): + """A locally-coplanar neighbourhood must be escaped, not surrendered to. + + Most source points lie on the plane z = 0, with a sparse off-plane set + further away. A target near the plane has all its nearest neighbours in + it, so the affine block is rank-deficient and the z-gradient of a linear + field is simply not representable from that stencil. Widening reaches the + off-plane points and recovers exactness. + + This matters because a single non-exact node sitting among exact ones is + an isolated spike, which damages derivatives far more than a smooth error + of the same magnitude. Measured on a 3D low-density swarm before the + retry existed: 2.2e-3 max error from 2 nodes in 28824. + """ + # A SMALL coplanar patch -- a local accident, which is what widening is + # for. (A globally coplanar cloud is a different situation, covered by + # test_locally_unrecoverable_cloud_falls_back_loudly.) + axis = np.linspace(0.45, 0.55, 5) + gx, gy = np.meshgrid(axis, axis, indexing="ij") + in_plane = np.stack([gx.ravel(), gy.ravel(), np.zeros(gx.size)], axis=1) + + rng = np.random.default_rng(19) + off_plane = rng.random((80, 3)) + + source = np.vstack([in_plane, off_plane]) + target = np.array([[0.5, 0.5, 0.002], [0.49, 0.51, -0.001]]) + + nnn = 8 + _, stencil = _stencils(source, target, nnn) + _, degenerate = linear_exact_weights(target, stencil) + assert degenerate.all(), ( + "test setup no longer produces coplanar stencils, so it is not " + "exercising the retry path" + ) + + # A z-dependent linear field cannot be recovered from a coplanar stencil, + # so exactness here proves the wider neighbourhood was actually used. + data = _linear(source)[:, None] + kdt = uw.kdtree.KDTree(source) + got = kdt.rbf_interpolator_local(target, data, nnn, 2, False, order=1) + + error = np.abs(got[:, 0] - _linear(target)).max() + assert error < 1.0e-12, ( + f"degenerate stencils were not recovered by widening: error {error:.3e}" + ) + + +def test_locally_unrecoverable_cloud_falls_back_loudly(): + """When widening cannot help, say so rather than pretending. + + A source cloud that is coplanar over a wide region is genuinely unable to + determine an affine fit in the third direction. Widening the stencil is + futile, and the honest outcome is a bounded inverse-distance answer plus a + warning -- not a silent loss of the guarantee (the failure mode of #424). + """ + axis = np.linspace(0.0, 1.0, 22) + gx, gy = np.meshgrid(axis, axis, indexing="ij") + source = np.stack([gx.ravel(), gy.ravel(), np.zeros(gx.size)], axis=1) + target = np.array([[0.5, 0.5, 0.02]]) + + data = _linear(source)[:, None] + kdt = uw.kdtree.KDTree(source) + + with pytest.warns(UserWarning, match="could not support an affine fit"): + got = kdt.rbf_interpolator_local(target, data, 8, 2, False, order=1) + + assert np.isfinite(got).all(), "an unrecoverable stencil must not give NaN" + lo, hi = data.min(), data.max() + assert lo - 1e-12 <= got[0, 0] <= hi + 1e-12, ( + "the inverse-distance fallback should still be bounded" + ) + + @pytest.mark.parametrize("dim", [2, 3]) def test_nnn_below_dim_plus_two_is_rejected(dim): """At dim+1 the affine tail is exactly determined and the scheme reduces to From 44ab4fbbc9420db3f1e348d9289bd8b1c9aebd53 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 27 Jul 2026 14:45:58 +1000 Subject: [PATCH 5/5] Add the operator form; make nnn default per dimension; test with random data Acting on review from the multigrid-transfer session, which is the first non-swarm consumer. 1. `KDTree.interpolation_matrix(coords, nnn=None, p=2, order=0)` returns the sparse transfer operator directly. Previously the only route to an operator was the raw `linear_exact_weights` kernel -- which returns ZEROS plus a mask for degenerate rows, because it has no kd-tree and so can neither widen nor fall back. That put the sharp edge in front of exactly the serious consumers, and a zero row silently interpolates to zero, which is worse than an inaccurate answer. The operator path shares the value path's handling: widen, 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; a test asserts `T @ data` equals the value API for both orders and both dimensions. 2. The default nnn now depends on order AND dimension. It was fixed at 4, which is below the enforced minimum of dim+2 at order=1 in 3D -- so `rbf_interpolator_local(coords, data, order=1)` raised on its own default. Enforcing a minimum is not the same as defaulting well. order=0 keeps 4. 3. Tests now use RANDOM data against an independent implementation. Linear reproduction cannot detect wrong-neighbour selection: sum_j w_j x_j = x* is checked against whatever points were handed to the solver, so it holds even if the tree returned the wrong ones, and a linear field then interpolates exactly anyway. Added a brute-force kNN + dense-solve oracle that shares no code with the library, plus a control asserting that shuffled stencils DO change the answer -- without which the test proves nothing. 4. Stated the three contract facts a consumer cannot infer: repeated calls and separate trees over equal points are bit-identical; tied kNN distances resolve deterministically (which point wins is unspecified, that it is stable is not); and row-wise kNN gives NO column guarantee, so a Galerkin PtAP can still be singular (#424) -- not this kernel's to fix, but a consumer should not have to discover it. 5. Documented the silent-failure facts: .coords_nd never .coords (#426), no units on output, ~5x cost, and do not validate against uw.function.evaluate (#432 -- P1 values on cell edges in 3D, and proxy nodes sit there routinely). Caught by the refactor and now pinned: swapping query(k=1) for find_closest_n_points added a stencil axis to the nnn=1 checkpoint path. Only a snapshot test noticed; it has its own test now. level_1 511 passed, level_2 312 passed, np=2 parallel 4 passed, np=4 scale sweep unchanged at 1.6e-15 worst case. Underworld development team with AI support from Claude Code --- docs/developer/subsystems/interpolation.md | 65 ++++- src/underworld3/ckdtree.pyx | 306 +++++++++++++++------ src/underworld3/utilities/rbf_stencil.py | 16 +- tests/test_0102_kdtree_linear_rbf.py | 215 +++++++++++++++ 4 files changed, 518 insertions(+), 84 deletions(-) diff --git a/docs/developer/subsystems/interpolation.md b/docs/developer/subsystems/interpolation.md index a3a912c19..6dc77eaed 100644 --- a/docs/developer/subsystems/interpolation.md +++ b/docs/developer/subsystems/interpolation.md @@ -94,6 +94,39 @@ inverse-distance error is unchanged to two significant figures, because the error is set by the asymmetry of the neighbourhood, not by how many points are in it. +## Getting the operator, not just the values + +When the same transfer is applied more than once — or when the operator *is* +the product, as for a multigrid prolongation — build it once: + +```python +T = kdt.interpolation_matrix(target_coords, order=1) # scipy CSR +values = T @ data # identical to the value API +``` + +The weights depend only on geometry, so one build serves every field and every +component; the value API re-solves each call. `T @ data` is asserted equal to +`rbf_interpolator_local(...)` in the test suite so the two cannot drift. + +```{warning} +`rbf_stencil.linear_exact_weights` is the raw kernel underneath, and it returns +**zeros** for degenerate rows plus a mask. It has no kd-tree, so it can neither +widen a stencil nor fall back. A caller that ignores the mask silently +interpolates to zero. Use `interpolation_matrix` unless you specifically need +the bare weights — its rows are never empty. +``` + +### Row-wise construction carries no column guarantee + +Every row has `nnn` non-zeros. **Nothing guarantees every column is non-empty.** +A source point that is not among the `nnn` nearest neighbours of any target +produces an empty column, and a consumer forming a Galerkin coarse operator +$P^{T} A P$ then gets a singular matrix — the failure mode of issue #424. + +This is a property of *any* row-wise kNN construction, not of this kernel, and +it is not fixed by the scheme being linear-exact. Consumers that need full +column rank must check for and repair empty columns themselves. + ## Choosing `nnn` `order=1` requires **`nnn >= dim + 2`** and raises `ValueError` below that. @@ -105,7 +138,37 @@ the neighbour simplex**. That is singular whenever those `dim + 1` points are collinear (2D) or coplanar (3D) — common near boundaries and on graded meshes. The swarm proxy default is `nnn = 2 * (dim + 1)` (6 in 2D, 8 in 3D), which -leaves enough slack that degenerate neighbourhoods are rare. +leaves enough slack that degenerate neighbourhoods are rare. That is also what +`nnn=None` resolves to at `order=1`; `order=0` keeps its historical default of +4. The default has to depend on `order` *and* dimension — a fixed default of 4 +is below the enforced minimum in 3D and would simply raise. + +## Things that fail silently if you get them wrong + +- **Build kd-trees from `.coords_nd`, never `.coords`.** `MeshVariable.coords` + dimensionalises once the model has reference quantities set, while particle + coordinates and `MeshVariable._get_kdtree` are non-dimensional. Mixing them + raises from `_convert_coords_to_tree_units` (issue #426). +- **The output carries no units.** Coordinates are converted into the tree's + frame, but `data` passes through untouched and the return is a plain array. + Re-attach units at the boundary if the caller needs them. +- **`order=1` costs about 5x `order=0`** — a dense $(nnn+d+1)^3$ solve per + target point. Irrelevant if you build an operator once; think twice in a hot + loop. +- **Do not validate against `uw.function.evaluate`.** It returns wrong values + for P1 fields at points lying exactly on cell edges in 3D (issue #432), and + proxy nodes sit on cell boundaries routinely. Compare against an analytic + field, or against an independent implementation, instead. + +## Determinism + +The weights are a pure function of the geometry: no global state, no random +numbers, no accumulation across calls. Repeated calls, and separate `KDTree` +instances built from equal point sets, return bit-identical results. + +Where several source points are exactly equidistant from a target, *which* of +them the kNN search selects is unspecified — but the choice is deterministic, +so results reproduce run to run. Both properties are pinned by tests. ## Degenerate stencils diff --git a/src/underworld3/ckdtree.pyx b/src/underworld3/ckdtree.pyx index 8abc9b931..c52f6be50 100644 --- a/src/underworld3/ckdtree.pyx +++ b/src/underworld3/ckdtree.pyx @@ -433,7 +433,7 @@ cdef class KDTree: def rbf_interpolator_local(self, coords, data, - nnn = 4, + nnn = None, p=2, verbose = False, order = 0, @@ -469,9 +469,11 @@ cdef class KDTree: Known data values at KD-tree points. Shape should be ``(n_points,)`` or ``(n_points, n_components)``. nnn : int, optional - Number of nearest neighbours to use for interpolation (default 4). - If 1, returns raw nearest-neighbour values without distance weighting. - ``order=1`` requires ``nnn >= dim + 2``. + Number of nearest neighbours to use. Defaults to 4 for ``order=0`` + and ``2 * (dim + 1)`` for ``order=1`` — the default has to depend + on both, because ``order=1`` requires ``nnn >= dim + 2`` and a + fixed default of 4 would raise in 3D. If 1, returns raw + 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. @@ -632,8 +634,15 @@ cdef class KDTree: return Values - def _retry_degenerate_stencils(self, coords_converted, data, tree_points, - vals, degenerate, nnn): + def _resolve_nnn(self, nnn, order): + """Default stencil size. ``order=1`` needs at least ``dim + 2``, so a + fixed default cannot serve both schemes (it would raise in 3D).""" + if nnn is not None: + return nnn + return 2 * (self.ndim + 1) if order == 1 else 4 + + def _retry_degenerate_stencils(self, coords_converted, tree_points, + degenerate, nnn): """Re-solve the failed points on a wider stencil. A neighbourhood that cannot support an affine fit is usually a local @@ -641,14 +650,15 @@ cdef class KDTree: Only the failed points are re-queried, so the cost is set by how many actually failed (measured: a handful in tens of thousands). - Returns the updated values and the mask of points that failed even on - the wider stencil, which keep their inverse-distance values. + Returns ``(rows, indices, weights, remaining)`` — the recovered rows and + their wider stencils, plus the mask of points that failed even after + widening and so keep their inverse-distance weights. """ from underworld3.utilities.rbf_stencil import linear_exact_weights wide = min(4 * nnn, self.n) if wide <= nnn: - return vals, degenerate + return None, None, None, degenerate failed = np.flatnonzero(degenerate) targets = np.ascontiguousarray(coords_converted[failed], dtype=np.float64) @@ -657,21 +667,197 @@ cdef class KDTree: # unit-aware `query` would reject it when the tree carries units. indices, _ = self.find_closest_n_points(wide, targets) if np.any(indices >= self.n): - return vals, degenerate + return None, None, None, degenerate weights, still_degenerate = linear_exact_weights( targets, tree_points[indices] ) + remaining = degenerate.copy() + remaining[failed] = still_degenerate + recovered = ~still_degenerate - if recovered.any(): - vals[failed[recovered]] = np.einsum( - "sdc,sd->sc", data[indices[recovered]], weights[recovered] + if not recovered.any(): + return None, None, None, remaining + + return (failed[recovered], indices[recovered], weights[recovered], + remaining) + + def _local_stencil(self, coords_converted, nnn, p, order): + """Neighbour indices and weights for every target point. + + The single source of the interpolation weights: both + :meth:`rbf_interpolator_local` and :meth:`interpolation_matrix` are + built on this, so the operator and the values it produces cannot drift + apart. + + Returns + ------- + indices, weights : ndarray + Shape ``(n_targets, nnn)``. The primary stencil. + wide : tuple or None + ``(rows, indices, weights)`` for targets re-solved on a wider + stencil after their primary one proved degenerate. Their primary + row is superseded and must be discarded by the caller. + degenerate : ndarray + Bool, shape ``(n_targets,)``. Targets still degenerate after + widening; these carry inverse-distance weights, which are bounded + but not linear-exact. + """ + # find_closest_n_points returns (indices, dist_sqr) -- the reverse of + # query()'s (dist, indices). It is used here rather than query() + # because coords_converted is already in tree units, and query() would + # convert a second time. + closest_n, distance_n = self.find_closest_n_points( + nnn, np.ascontiguousarray(coords_converted, dtype=np.float64) + ) + + # valid indices are 0..n-1; the empty-tree sentinel (0 with n=0) + # must trip this guard, so the comparison is >= (issue #399). + if np.any(closest_n >= self.n): + raise RuntimeError( + "Error in rbf_interpolator_local_from_kdtree - a nearest neighbour wasn't found" ) - remaining = degenerate.copy() - remaining[failed] = still_degenerate - return vals, remaining + # np.bool_, not bool: this module cimports the C++ `bool` from libcpp, + # which shadows the Python builtin and will not compile here. + degenerate = np.zeros(coords_converted.shape[0], dtype=np.bool_) + + if nnn == 1: + 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. + epsilon = 1e-12 + weights = 1 / np.power(epsilon + distance_n[:], p) + n_weights = (weights.T / np.sum(weights, axis=1)).T + + if order == 0: + return closest_n, n_weights, None, degenerate + + from underworld3.utilities.rbf_stencil import linear_exact_weights + + tree_points = np.asarray(self.points) + linear_weights, degenerate = linear_exact_weights( + coords_converted, tree_points[closest_n[:]] + ) + # A stencil that cannot support an affine fit (collinear in 2D, + # coplanar in 3D) keeps the inverse-distance weights: finite and + # bounded, but not linear-exact. + linear_weights[degenerate] = n_weights[degenerate] + + wide = None + if degenerate.any(): + # Widen before surrendering. A single non-exact node sitting among + # exact ones is an isolated SPIKE, far more damaging to derivatives + # and integrals than a smooth error of the same size. + rows, wide_indices, wide_weights, degenerate = ( + self._retry_degenerate_stencils( + coords_converted, tree_points, degenerate, nnn + ) + ) + if rows is not None: + wide = (rows, wide_indices, wide_weights) + + n_degenerate = int(degenerate.sum()) + if n_degenerate: + import warnings + + warnings.warn( + f"rbf_interpolator_local(order=1): {n_degenerate} of " + f"{degenerate.size} stencils could not support an affine fit " + "(collinear/coplanar neighbours) even after widening, and fell " + "back to inverse-distance weighting.", + stacklevel=3, + ) + + return closest_n, linear_weights, wide, degenerate + + def interpolation_matrix(self, coords, nnn=None, p=2, order=0): + """Sparse operator mapping values on the KD-tree points to ``coords``. + + ``T @ data`` is exactly what :meth:`rbf_interpolator_local` returns for + the same arguments — both are built from the same weights. + + Use this instead of the value API when the same transfer is applied + more than once, or when the operator itself is the product (a multigrid + prolongation, for example). The weights depend only on geometry, so one + build serves every field and every component; the value API re-solves + each time. + + Parameters + ---------- + coords : array-like + Target coordinates, shape ``(n_coords, dim)``. + nnn : int, optional + Neighbours per target. Defaults to 4 for ``order=0`` and + ``2 * (dim + 1)`` for ``order=1``. + p : int, optional + Inverse-distance power; used by ``order=0`` only. + order : int, optional + Polynomial reproduction order, 0 or 1. See + :meth:`rbf_interpolator_local`. + + Returns + ------- + scipy.sparse.csr_matrix + Shape ``(n_coords, self.n)``. Rows carry ``nnn`` non-zeros, except + those re-solved on a wider stencil. + + Notes + ----- + Degenerate stencils are handled exactly as in the value path — widened, + then fall back to inverse-distance weights with a warning. Rows are + never empty. + + **Row-wise construction gives no column guarantee.** Every row has + ``nnn`` non-zeros, but a source point that is not among the neighbours + of any target produces an empty column. Consumers that need full column + rank -- a Galerkin coarse operator :math:`P^{T} A P`, for instance -- + must check for and repair empty columns themselves (see issue #424). + """ + import scipy.sparse as sp + + if order not in (0, 1): + raise ValueError( + f"order must be 0 (inverse distance) or 1 (linear-exact), got {order!r}." + ) + coords_converted = self._convert_coords_to_tree_units(coords) + if coords_converted.shape[1] != self.ndim: + raise RuntimeError( + f"Interpolation coordinates dimensionality " + f"({coords_converted.shape[1]}) is different to kD-tree " + f"dimensionality ({self.ndim})." + ) + nnn = self._resolve_nnn(nnn, order) + + indices, weights, wide, _ = self._local_stencil( + coords_converted, nnn, p, order + ) + + n_targets = coords_converted.shape[0] + keep = np.ones(n_targets, dtype=np.bool_) + if wide is not None: + keep[wide[0]] = False + + rows = np.repeat(np.flatnonzero(keep), nnn) + cols = indices[keep].ravel() + vals = weights[keep].ravel() + + if wide is not None: + wide_rows, wide_indices, wide_weights = wide + rows = np.concatenate( + [rows, np.repeat(wide_rows, wide_indices.shape[1])] + ) + cols = np.concatenate([cols, wide_indices.ravel()]) + vals = np.concatenate([vals, wide_weights.ravel()]) + + return sp.csr_matrix( + (vals, (rows, cols.astype(np.int64))), shape=(n_targets, self.n) + ) def rbf_interpolator_local_from_kdtree(self, coords, data, nnn, p, verbose, order=0, monotone=False): @@ -729,29 +915,11 @@ cdef class KDTree: f"Data does not match kd-tree size array ({data.shape[0]} v ({self.n}))" ) - coords_contiguous = np.ascontiguousarray(coords_converted) - # query nnn points to the coords - # distance_n is a list of distance to the nearest neighbours for all coords_contiguous - # closest_n is the index of the neighbours from ncoords for all coords_contiguous - # Note: query() returns sqr_dists=True by default, and we use the converted coords - distance_n, closest_n = self.query(coords, k=nnn) - - # valid indices are 0..n-1; the empty-tree sentinel (0 with n=0) - # must trip this guard, so the comparison is >= (issue #399). - if np.any(closest_n >= self.n): - raise RuntimeError( - "Error in rbf_interpolator_local_from_kdtree - a nearest neighbour wasn't found" - ) + nnn = self._resolve_nnn(nnn, order) if verbose and uw.mpi.rank == 0: - # For Debugging - # print(f"kd-tree diagnostics: d.shape - {distance_n.shape}, c.shape - {closest_n.shape}") print(f"Mapping values with nnn - {nnn} & p {p} ... start", flush=True) - # np.bool_, not bool: this module cimports the C++ `bool` from libcpp, - # which shadows the Python builtin and will not compile here. - degenerate = np.zeros(coords_converted.shape[0], dtype=np.bool_) - if nnn == 1: # only use nearest neighbour raw data if order == 1: @@ -759,57 +927,31 @@ cdef class KDTree: "order=1 needs at least dim + 2 neighbours to determine the " f"affine tail; nnn=1 selects the raw nearest-neighbour path." ) + closest_n, _ = self.find_closest_n_points( + 1, np.ascontiguousarray(coords_converted, dtype=np.float64) + ) + # (n, 1) -> (n,): the nearest-neighbour path returns data rows + # directly, so the stencil axis must not survive into the result. + # query(k=1) used to do this reshape for us. + closest_n = closest_n.reshape(-1) + if np.any(closest_n >= self.n): + raise RuntimeError( + "Error in rbf_interpolator_local_from_kdtree - a nearest neighbour wasn't found" + ) return data[closest_n] - # 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 (query - # defaults to sqr_dists=True), so the decay is r^(-2p), not the - # documented r^(-p), and `epsilon` floors r at ~1e-6 rather than 1e-12. - epsilon = 1e-12 - weights = 1 / np.power(epsilon + distance_n[:], p) - n_weights = (weights.T / np.sum(weights, axis=1)).T + closest_n, n_weights, wide, degenerate = self._local_stencil( + coords_converted, nnn, p, order + ) kdata = data[closest_n[:]] - if order == 1: - from underworld3.utilities.rbf_stencil import linear_exact_weights - - tree_points = np.asarray(self.points) - stencil_coords = tree_points[closest_n[:]] - linear_weights, degenerate = linear_exact_weights( - coords_converted, stencil_coords - ) - # A stencil that cannot support an affine fit (collinear in 2D, - # coplanar in 3D) keeps the inverse-distance weights for now: finite - # and bounded, but not linear-exact. - linear_weights[degenerate] = n_weights[degenerate] - n_weights = linear_weights - # magic with einstein summation power vals = np.einsum("sdc,sd->sc", kdata, n_weights) - # print(valz) - - if order == 1 and degenerate.any(): # noqa: F821 - set in the order==1 block - # Widen the stencil for the few points that failed, rather than - # accepting inverse distance there. A degenerate neighbourhood is - # usually a local accident of point placement, and reaching further - # fixes it -- worth doing because a single non-exact node sitting - # among exact ones is an isolated SPIKE, which is far more damaging - # to derivatives and integrals than a smooth error of the same size. - vals, degenerate = self._retry_degenerate_stencils( - coords_converted, data, tree_points, vals, degenerate, nnn - ) - - n_degenerate = int(degenerate.sum()) - if n_degenerate: - import warnings - warnings.warn( - f"rbf_interpolator_local(order=1): {n_degenerate} of " - f"{degenerate.size} stencils could not support an affine fit " - "(collinear/coplanar neighbours) even after widening, and fell " - "back to inverse-distance weighting.", - stacklevel=2, + if wide is not None: + wide_rows, wide_indices, wide_weights = wide + vals[wide_rows] = np.einsum( + "sdc,sd->sc", data[wide_indices], wide_weights ) if monotone_clamp: @@ -829,8 +971,12 @@ cdef class KDTree: # reproduction guarantee. from underworld3.utilities.rbf_stencil import affine_trend + # The bound comes from the primary stencil even for the few + # rows re-solved on a wider one: it is still a legitimate local + # bound, and the affine trend -- the part being preserved -- is + # what matters here. trend_at_target, trend_at_stencil = affine_trend( - coords_converted, stencil_coords, kdata + coords_converted, np.asarray(self.points)[closest_n[:]], kdata ) residual = kdata - trend_at_stencil vals = trend_at_target + np.clip( diff --git a/src/underworld3/utilities/rbf_stencil.py b/src/underworld3/utilities/rbf_stencil.py index 14895c46a..20f4acf4e 100644 --- a/src/underworld3/utilities/rbf_stencil.py +++ b/src/underworld3/utilities/rbf_stencil.py @@ -137,9 +137,19 @@ def linear_exact_weights(target_coords, neighbour_coords): degenerate : numpy.ndarray Boolean, shape ``(n_targets,)``. True where the stencil could not support an affine fit — collinear neighbours in 2D, coplanar in 3D, - or a stencil that collapsed onto its target. **Rows flagged here are - returned as zeros**; the caller must substitute its own fallback - (Shepard weights over the same neighbours are the natural choice). + or a stencil that collapsed onto its target. + + Warnings + -------- + **Rows flagged ``degenerate`` are returned as ZEROS.** A caller that + ignores the mask gets a silent interpolation to zero, which is worse than + an inaccurate answer. This is a raw kernel: it has no kd-tree, so it can + neither widen the stencil nor fall back to inverse distance. + + Prefer :meth:`underworld3.ckdtree.KDTree.interpolation_matrix` unless you + genuinely need the bare weights. It builds the same weights, then widens + degenerate stencils and falls back to inverse distance with a warning, so + its rows are never empty. Notes ----- diff --git a/tests/test_0102_kdtree_linear_rbf.py b/tests/test_0102_kdtree_linear_rbf.py index 19f6b88e4..652861a87 100644 --- a/tests/test_0102_kdtree_linear_rbf.py +++ b/tests/test_0102_kdtree_linear_rbf.py @@ -333,6 +333,221 @@ def test_kdtree_rejects_bad_order_and_monotone_mode(): monotone="pick") +# -------------------------------------------------------------------------- +# Oracle tests on a RANDOM field. +# +# Linear reproduction cannot detect wrong-neighbour selection: the identity +# sum_j w_j x_j = x* is evaluated against whatever points were handed to the +# solver, so it holds even if the kd-tree returned the wrong neighbours, and +# a linear field then interpolates exactly anyway. Only a field whose values +# vary independently of position -- i.e. random data -- makes the choice of +# neighbour observable. +# -------------------------------------------------------------------------- +def _reference_interpolate(source, target, values, nnn): + """An independent implementation, written the obvious slow way. + + Brute-force kNN, then one dense saddle solve per target point, with no + local rescaling. Deliberately does not share code with the library. + """ + out = np.empty(target.shape[0]) + dim = source.shape[1] + for i, x in enumerate(target): + d = np.linalg.norm(source - x, axis=1) + idx = np.argsort(d)[:nnn] + pts = source[idx] + + def phi(r): + r = np.where(r == 0.0, 1.0e-300, r) + return r ** 2 * np.log(r) + + rr = np.linalg.norm(pts[:, None, :] - pts[None, :, :], axis=2) + P = np.hstack([np.ones((nnn, 1)), pts]) + M = np.block([[phi(rr), P], [P.T, np.zeros((dim + 1, dim + 1))]]) + rhs = np.concatenate([phi(np.linalg.norm(pts - x, axis=1)), [1.0], x]) + w = np.linalg.solve(M, rhs)[:nnn] + out[i] = w @ values[idx] + return out + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_matches_an_independent_implementation_on_random_data(dim): + """Random values, so wrong-neighbour selection cannot hide.""" + nnn = 2 * (dim + 1) + rng = np.random.default_rng(31337 + dim) + source = rng.random((400, dim)) + target = 0.15 + 0.7 * rng.random((40, dim)) + values = rng.standard_normal(source.shape[0]) + + kdt = uw.kdtree.KDTree(source) + got = kdt.rbf_interpolator_local(target, values[:, None], nnn, 2, False, order=1) + want = _reference_interpolate(source, target, values, nnn) + + error = np.abs(got[:, 0] - want).max() + assert error < 1.0e-8, ( + f"disagrees with an independent implementation by {error:.3e} on " + "random data — check neighbour selection, not just reproduction" + ) + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_kdtree_neighbours_match_brute_force_on_random_data(dim): + """The value must change if the stencil is wrong, and it does not.""" + nnn = 2 * (dim + 1) + rng = np.random.default_rng(555 + dim) + source = rng.random((600, dim)) + target = 0.2 + 0.6 * rng.random((50, dim)) + values = rng.standard_normal(source.shape[0]) + + kdt = uw.kdtree.KDTree(source) + got = kdt.rbf_interpolator_local(target, values[:, None], nnn, 2, False, order=1) + + idx, stencil = _stencils(source, target, nnn) + weights, degenerate = linear_exact_weights(target, stencil) + assert not degenerate.any() + want = (weights * values[idx]).sum(axis=1) + + assert np.abs(got[:, 0] - want).max() < 1.0e-10 + + # Control: a deliberately wrong stencil must give a different answer, or + # the test above proves nothing. + rolled = np.roll(idx, 1, axis=0) + wrong_weights, _ = linear_exact_weights(target, source[rolled]) + wrong = (wrong_weights * values[rolled]).sum(axis=1) + assert np.abs(wrong - want).max() > 1.0e-3, ( + "shuffling the stencils did not change the result, so this test " + "cannot detect wrong-neighbour selection" + ) + + +# -------------------------------------------------------------------------- +# The operator form +# -------------------------------------------------------------------------- +@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.""" + 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) + + assert T.shape == (target.shape[0], source.shape[0]) + assert np.abs(T @ values - direct).max() < 1.0e-12 + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_interpolation_matrix_rows_are_never_empty(dim): + """The raw-weights helper zeroes degenerate rows; the operator must not. + + A zero row would silently interpolate to zero, which is worse than the + inverse-distance fallback it replaces. + """ + axis = np.linspace(0.0, 1.0, 12) + grids = np.meshgrid(*([axis] * dim), indexing="ij") + source = np.stack([g.ravel() for g in grids], axis=1) + rng = np.random.default_rng(2) + target = rng.random((80, dim)) + + kdt = uw.kdtree.KDTree(source) + T = kdt.interpolation_matrix(target, order=1) + + per_row = np.diff(T.indptr) + assert per_row.min() > 0, "empty row in the interpolation operator" + row_sums = np.asarray(T.sum(axis=1)).ravel() + assert np.abs(row_sums - 1.0).max() < 1.0e-12 + + +# -------------------------------------------------------------------------- +# Contract: determinism and purity +# -------------------------------------------------------------------------- +def test_repeated_calls_are_bit_identical(): + """No global state, no RNG, no accumulation: the weights are a pure + function of the geometry.""" + rng = np.random.default_rng(64) + source = rng.random((300, 3)) + target = rng.random((40, 3)) + values = rng.standard_normal((300, 1)) + + kdt = uw.kdtree.KDTree(source) + first = kdt.rbf_interpolator_local(target, values, 8, 2, False, order=1) + second = kdt.rbf_interpolator_local(target, values, 8, 2, False, order=1) + assert np.array_equal(first, second) + + other = uw.kdtree.KDTree(source.copy()) + assert np.array_equal( + other.rbf_interpolator_local(target, values, 8, 2, False, order=1), first + ) + + +def test_equidistant_neighbours_resolve_deterministically(): + """Ties in the kNN search must not make the result run-dependent. + + Which of several equidistant source points is chosen is not specified — + but it must be the same choice every time, or results stop reproducing. + """ + # A target at the centre of a symmetric ring: all ring points are exactly + # equidistant, so the nnn selection is a tie. + angles = np.arange(12) * (2.0 * np.pi / 12.0) + ring = np.stack([np.cos(angles), np.sin(angles)], axis=1) + source = np.vstack([ring, 2.0 * ring]) + target = np.array([[0.0, 0.0]]) + values = np.arange(source.shape[0], dtype=float)[:, None] + + kdt = uw.kdtree.KDTree(source) + results = [ + kdt.rbf_interpolator_local(target, values, 6, 2, False, order=1) + for _ in range(5) + ] + for r in results[1:]: + assert np.array_equal(r, results[0]), "tied kNN selection is not deterministic" + + fresh = uw.kdtree.KDTree(source.copy()) + assert np.array_equal( + fresh.rbf_interpolator_local(target, values, 6, 2, False, order=1), results[0] + ) + + +def test_nearest_neighbour_path_keeps_its_shape(): + """`nnn=1` returns source rows directly, with no stencil axis. + + This is the checkpoint round-trip path (`read_timestep` uses `nnn=1`). + A refactor that swapped `query(k=1)` for `find_closest_n_points` silently + added an axis here, and only a snapshot test noticed. + """ + rng = np.random.default_rng(12) + source = rng.random((200, 2)) + target = rng.random((30, 2)) + data = rng.standard_normal((200, 1)) + + kdt = uw.kdtree.KDTree(source) + got = kdt.rbf_interpolator_local(target, data, 1, 2, False) + + assert got.shape == (30, 1), f"nnn=1 returned shape {got.shape}, expected (30, 1)" + + # It really is the nearest neighbour's value, not an average. + d = np.linalg.norm(target[:, None, :] - source[None, :, :], axis=2) + assert np.array_equal(got[:, 0], data[np.argmin(d, axis=1), 0]) + + +@pytest.mark.parametrize("order", [0, 1]) +def test_default_nnn_works_in_both_dimensions(order): + """A default that raises is not a default (the fixed nnn=4 did, at + order=1 in 3D, where dim + 2 = 5).""" + for dim in (2, 3): + rng = np.random.default_rng(9) + source = rng.random((300, dim)) + target = rng.random((20, dim)) + data = _linear(source)[:, None] + kdt = uw.kdtree.KDTree(source) + got = kdt.rbf_interpolator_local(target, data, order=order) + assert got.shape == (20, 1) + assert np.isfinite(got).all() + + @pytest.mark.parametrize("dim", [2, 3]) def test_error_falls_with_point_spacing_faster_than_inverse_distance(dim): """On a quadratic field, order=1 should converge under refinement and