From 0bc03d4a2bf83abfc3ffa292261605e350b821f5 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Thu, 29 Oct 2026 22:01:35 +0000 Subject: [PATCH] feat: exact kNN mesh magnification areas by quadrature of the interpolant KNearestNeighbor and KNNBarycentric used the Delaunay barycentric dual areas, which integrate only the linear Delaunay interpolant exactly. Their own interpolants differ, biasing magnification by up to +6.10% lensed and -5.11% unlensed against simulator truth. New MeshGeometryKNearestNeighbor (built by the kNN interpolators) computes areas_for_magnification as the integral of each vertex's basis function over the convex hull: the interpolator's own weights at the centres of a 256^2 grid, masked to the hull, summed per vertex and rescaled to the hull area. Host-side and cached on first access, like the Delaunay fallback, so fits never pay for it. Identity-lens error: KNN 1.6e-2 -> 1.3e-5, KNNBarycentric 4.4e-2 -> 1.8e-3. The partial-block neighbour bug (#609) is deliberately not fixed here; the quadrature integrates whatever neighbours the interpolator selects. Closes #665 Co-Authored-By: Claude Opus 5.5 --- autoarray/inversion/mesh/interpolator/knn.py | 19 +++++ .../inversion/mesh/mesh_geometry/delaunay.py | 84 +++++++++++++++++++ .../test_areas_for_magnification.py | 21 +++-- 3 files changed, 117 insertions(+), 7 deletions(-) diff --git a/autoarray/inversion/mesh/interpolator/knn.py b/autoarray/inversion/mesh/interpolator/knn.py index 9180c798a..2d4967ff8 100644 --- a/autoarray/inversion/mesh/interpolator/knn.py +++ b/autoarray/inversion/mesh/interpolator/knn.py @@ -150,6 +150,25 @@ class InterpolatorKNearestNeighbor(InterpolatorDelaunay): # Excluding them would change pixel signals broadly and add a jump in the likelihood (PyAutoArray#627). exclude_zero_weight_mappings = False + @cached_property + def mesh_geometry(self): + """ + The kNN geometry: Delaunay vertex geometry whose + `areas_for_magnification` integrates this interpolator's own basis + functions (see `MeshGeometryKNearestNeighbor`). + """ + from autoarray.inversion.mesh.mesh_geometry.delaunay import ( + MeshGeometryKNearestNeighbor, + ) + + return MeshGeometryKNearestNeighbor( + mesh=self.mesh, + mesh_grid=self.mesh_grid, + data_grid=self.data_grid, + interpolator_cls=type(self), + xp=self._xp, + ) + @property def dual_areas(self): """ diff --git a/autoarray/inversion/mesh/mesh_geometry/delaunay.py b/autoarray/inversion/mesh/mesh_geometry/delaunay.py index 4bc332994..ad3786bf3 100644 --- a/autoarray/inversion/mesh/mesh_geometry/delaunay.py +++ b/autoarray/inversion/mesh/mesh_geometry/delaunay.py @@ -271,3 +271,87 @@ def areas_for_magnification(self) -> np.ndarray: scipy.spatial.Delaunay(mesh_grid_xy).simplices, xp=np, ) + + +# Quadrature cells per axis for `MeshGeometryKNearestNeighbor.areas_for_magnification`. +# On a 144-vertex mesh the identity-lens magnification moved by < 3e-5 between +# 200, 400 and 800 cells per axis (PyAutoArray#665); 256 costs ~1 s there. +KNN_AREA_QUADRATURE_PIXELS: int = 256 + + +class MeshGeometryKNearestNeighbor(MeshGeometryDelaunay): + def __init__(self, mesh, mesh_grid, data_grid, interpolator_cls, xp=np, **kwargs): + """ + The geometry of the kNN meshes (`KNearestNeighbor`, `KNNBarycentric`): + Delaunay vertex geometry, but magnification areas that integrate the + kNN interpolant itself. + + Parameters + ---------- + interpolator_cls + The kNN interpolator class that built this geometry; its basis + functions are evaluated on a quadrature grid to compute + `areas_for_magnification`. + """ + super().__init__( + mesh=mesh, + mesh_grid=mesh_grid, + data_grid=data_grid, + dual_areas=None, + xp=xp, + **kwargs, + ) + + self.interpolator_cls = interpolator_cls + + @cached_property + def areas_for_magnification(self) -> np.ndarray: + """ + The integral of every vertex's interpolation basis function over the + convex hull of the mesh, so `sum(reconstruction * areas_for_magnification)` + is the integral of the reconstructed source. + + The Delaunay barycentric dual areas are these integrals only for the linear + barycentric interpolant. The Wendland-kNN and 3-nearest barycentric + interpolants are not linear on the Delaunay triangles, and weighting by the + dual areas biased their magnification by up to ~6% (PyAutoArray#665). + + The integral is a midpoint quadrature: the interpolator's own mapping + weights are evaluated at the centres of a uniform + `KNN_AREA_QUADRATURE_PIXELS`^2 grid over the vertices' bounding box, the + centres outside the convex hull are dropped, and the weights times the + cell area are summed per vertex. The total is rescaled to the hull area, + which removes the cells straddling the hull edge from the sum (the same + domain the Delaunay dual areas tile). + + Computed host-side on first access, like the Delaunay dual-area fallback. + """ + import scipy.spatial + + nodes = np.asarray(self.mesh_grid_xy) + lo = nodes.min(axis=0) + hi = nodes.max(axis=0) + + n = KNN_AREA_QUADRATURE_PIXELS + u = (np.arange(n) + 0.5) / n + grid_0, grid_1 = np.meshgrid( + lo[0] + u * (hi[0] - lo[0]), lo[1] + u * (hi[1] - lo[1]), indexing="ij" + ) + points = np.stack([grid_0.ravel(), grid_1.ravel()], axis=1) + points = points[scipy.spatial.Delaunay(nodes).find_simplex(points) >= 0] + + interpolator = self.interpolator_cls( + mesh=self.mesh, + mesh_grid=self.mesh_grid, + data_grid=points, + xp=np, + ) + mappings, _, weights = interpolator._mappings_sizes_weights + + areas = np.bincount( + np.asarray(mappings).ravel(), + weights=np.asarray(weights).ravel(), + minlength=nodes.shape[0], + ) + + return areas * scipy.spatial.ConvexHull(nodes).volume / np.sum(areas) diff --git a/test_autoarray/inversion/pixelization/mesh_geometry/test_areas_for_magnification.py b/test_autoarray/inversion/pixelization/mesh_geometry/test_areas_for_magnification.py index 23db88d96..8db639fdd 100644 --- a/test_autoarray/inversion/pixelization/mesh_geometry/test_areas_for_magnification.py +++ b/test_autoarray/inversion/pixelization/mesh_geometry/test_areas_for_magnification.py @@ -76,11 +76,18 @@ # not the integral of the node's hat once the CDF bends across the cell (the # rank CDF's staircase on tied coordinates; a steep adapt-image CDF). # -# Documented approximations (flagged in PyAutoArray#659; phase 4b): -# - DelaunayNN / KNearestNeighbor / KNNBarycentric inherit the Delaunay dual -# areas, which are exact only for the linear interpolant; for the Sibson, -# Wendland-kNN and 3-nearest-barycentric interpolants they are an approximation -# [-4.3e-3, +1.6e-2, -4.4e-2]. +# - KNearestNeighbor / KNNBarycentric: the integral of each vertex's kNN basis +# function over the convex hull, by midpoint quadrature of the interpolator's +# own weights (`MeshGeometryKNearestNeighbor`, PyAutoArray#665) +# [+1.3e-5, +1.8e-3; KNNBarycentric's clipped weights are discontinuous, so +# the coarse sub-pixel `Q` is the limit there]. The Delaunay dual areas they +# used before measured [+1.6e-2, -4.4e-2]. The quadrature integrates whatever +# neighbours the interpolator selects, so it stays exact with or without the +# open partial-block neighbour bug PyAutoArray#609. +# +# Documented approximation: +# - DelaunayNN inherits the Delaunay dual areas, exact only for the linear +# interpolant; for its Sibson interpolant they are an approximation [-4.3e-3]. MAGNIFICATION_TOLERANCE = { "RectangularUniform": 2.5e-3, "RectangularBilinearAdaptDensity": 1.0e-3, @@ -89,8 +96,8 @@ "RectangularRTUAdaptImage": 1.0e-2, "Delaunay": 3.0e-4, "DelaunayNN": 1.0e-2, - "KNearestNeighbor": 3.0e-2, - "KNNBarycentric": 7.0e-2, + "KNearestNeighbor": 1.0e-4, + "KNNBarycentric": 4.0e-3, } RECTANGULAR_ADAPTIVE = [