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