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 <noreply@anthropic.com>

Author: Jammy2211

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):
         """

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)

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 = [