From 46aaac67b282f3464ece7836e68e75c6d01e06f8 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 4 Sep 2026 22:38:50 +0000 Subject: [PATCH] fix: Delaunay areas_for_magnification returns barycentric dual areas (#524) MeshGeometryDelaunay.areas_for_magnification returned scipy Voronoi cell areas (unbounded cells zeroed), but the Delaunay mapper is a barycentric-linear interpolant whose exact quadrature weights are the barycentric dual areas (sum of triangle_area / 3 per vertex). The phase 8 audit (#522) measured the Voronoi denominator biasing magnification by -13% to -99% on realistic meshes; the dual areas recover the identity-lens mu = 1.0 to 3e-5. - scipy_delaunay / jax_delaunay (and the Matern variants) now return the dual areas they already compute; DelaunayInterface carries them as `dual_areas`, InterpolatorDelaunay.mesh_geometry hands them to MeshGeometryDelaunay(dual_areas=), so the value is in-graph and jit-safe on the JAX path. The kNN interpolators return None and the geometry falls back to a host-side scipy triangulation only if asked. - areas_for_magnification returns the dual areas; voronoi_areas is unchanged and stays available under its own name. Docstring rewritten. - _plot_delaunay: the docstring no longer claims Gouraud shading (tripcolor defaults to flat, integral-equivalent), and the mapper's own simplices are passed as triangles= so the figure shows the mesh the inversion used. - Tests: the two phase-8 tests pinning the Voronoi semantics are flipped deliberately; new mapper-level identity test (sum(mapping_matrix @ s) * pixel_area == sum(s * areas), rel 3.4e-5, Voronoi off by 27%) and a JAX-vs-NumPy parity test. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01VVLuDrRfkZ81TBh8yVbFrN --- .../inversion/mesh/interpolator/delaunay.py | 117 ++++++++++++--- autoarray/inversion/mesh/interpolator/knn.py | 11 ++ .../inversion/mesh/mesh_geometry/delaunay.py | 91 ++++++++--- autoarray/plot/inversion.py | 33 +++- .../interpolator/test_delaunay.py | 63 +++++--- .../interpolator/test_delaunay_nn.py | 2 +- .../pixelization/mappers/test_delaunay.py | 141 ++++++++++++++++++ .../mesh_geometry/test_delaunay.py | 96 ++++++++---- 8 files changed, 461 insertions(+), 93 deletions(-) diff --git a/autoarray/inversion/mesh/interpolator/delaunay.py b/autoarray/inversion/mesh/interpolator/delaunay.py index 6c751b7d4..8d803f13c 100644 --- a/autoarray/inversion/mesh/interpolator/delaunay.py +++ b/autoarray/inversion/mesh/interpolator/delaunay.py @@ -9,7 +9,15 @@ def scipy_delaunay(points_np, query_points_np, areas_factor): - """Compute Delaunay simplices (simplices_padded) and Voronoi areas in one call.""" + """Compute the Delaunay simplices (``simplices_padded``), the query-point + mappings, the split-cross points and the barycentric dual areas in one call. + + The dual areas are returned as the sixth element: they weight the split + points here, and they are also the exact quadrature weights of the + barycentric interpolant, which + :meth:`~autoarray.inversion.mesh.mesh_geometry.delaunay.MeshGeometryDelaunay.areas_for_magnification` + consumes (PyAutoArray#524). + """ from scipy.spatial import Delaunay max_simplices = 2 * points_np.shape[0] @@ -60,7 +68,7 @@ def scipy_delaunay(points_np, query_points_np, areas_factor): delaunay_points=points_np, ) - return points, simplices_padded, mappings, split_points, splitted_mappings + return points, simplices_padded, mappings, split_points, splitted_mappings, areas # Query points are located in chunks of this size so the (chunk, N) distance @@ -287,6 +295,30 @@ def locate_chunk(q_chunk): return mappings +def _dual_areas_padded_jnp(points, simplices_padded): + """In-graph barycentric dual areas from the -1 padded simplex table. + + A masked scatter-add of ``triangle_area / 3`` into each of a triangle's + three vertices, skipping the padded (-1) rows. Equivalent to + ``barycentric_dual_area_from`` on the unpadded simplices, but written with + ``.at[].add`` so it stays inside the JIT program and differentiable with + respect to ``points``. + """ + import jax.numpy as jnp + + valid = simplices_padded[:, 0] >= 0 + s = simplices_padded.clip(min=0) + p0, p1, p2 = points[s[:, 0]], points[s[:, 1]], points[s[:, 2]] + tri_cross = (p1[:, 0] - p0[:, 0]) * (p2[:, 1] - p0[:, 1]) - ( + p1[:, 1] - p0[:, 1] + ) * (p2[:, 0] - p0[:, 0]) + contrib = jnp.where(valid, 0.5 * jnp.abs(tri_cross) / 3.0, 0.0) + areas = jnp.zeros(points.shape[0], dtype=points.dtype) + for k in range(3): + areas = areas.at[s[:, k]].add(contrib) + return areas + + def jax_delaunay(points, query_points, areas_factor=0.5): """JAX-path Delaunay construction. Only the qhull triangulation runs on the host (via ``pure_callback``); point location, dual areas and split @@ -307,16 +339,7 @@ def jax_delaunay(points, query_points, areas_factor=0.5): ) # dual areas via masked scatter-add over the padded simplices - valid = simplices_padded[:, 0] >= 0 - s = simplices_padded.clip(min=0) - p0, p1, p2 = points[s[:, 0]], points[s[:, 1]], points[s[:, 2]] - tri_cross = (p1[:, 0] - p0[:, 0]) * (p2[:, 1] - p0[:, 1]) - ( - p1[:, 1] - p0[:, 1] - ) * (p2[:, 0] - p0[:, 0]) - contrib = jnp.where(valid, 0.5 * jnp.abs(tri_cross) / 3.0, 0.0) - areas = jnp.zeros(points.shape[0], dtype=points.dtype) - for k in range(3): - areas = areas.at[s[:, k]].add(contrib) + areas = _dual_areas_padded_jnp(points, simplices_padded) split_points = split_points_from( points=points, @@ -336,7 +359,7 @@ def jax_delaunay(points, query_points, areas_factor=0.5): xp=jnp, ) - return points, simplices_padded, mappings, split_points, splitted_mappings + return points, simplices_padded, mappings, split_points, splitted_mappings, areas def barycentric_dual_area_from( @@ -441,12 +464,16 @@ def scipy_delaunay_matern(points_np, query_points_np): """ Minimal SciPy Delaunay callback for Matérn regularization. - Returns only what’s needed for mapping: + Returns only what’s needed for mapping, plus the dual areas: - points (tri.points) - simplices_padded - mappings: integer array of pixel indices for each query point, typically of shape (Q, 3), where each row gives the indices of the Delaunay mesh vertices ("pixels") associated with that query point. + - areas: the barycentric dual area of every vertex. Matérn + regularization does not need them (there are no split points), but the + magnification quadrature does, so they are returned here too + (PyAutoArray#524). """ from scipy.spatial import Delaunay @@ -472,13 +499,19 @@ def scipy_delaunay_matern(points_np, query_points_np): delaunay_points=points_np, ) - return points, simplices_padded, mappings + areas = barycentric_dual_area_from( + points, + simplices, + xp=np, + ) + + return points, simplices_padded, mappings, areas def jax_delaunay_matern(points, query_points): """JAX-path Matérn variant: qhull-only callback + JAX point location, - returning the same minimal (points, simplices_padded, mappings) contract - as ``scipy_delaunay_matern``.""" + returning the same minimal (points, simplices_padded, mappings, areas) + contract as ``scipy_delaunay_matern``.""" import jax.numpy as jnp simplices_padded, simplex_neighbors, vertex_simplex = _jax_delaunay_tables(points) @@ -492,7 +525,9 @@ def jax_delaunay_matern(points, query_points): xp=jnp, ) - return points, simplices_padded, mappings + areas = _dual_areas_padded_jnp(points, simplices_padded) + + return points, simplices_padded, mappings, areas def triangle_area_xp(c0, c1, c2, xp): @@ -589,14 +624,34 @@ def pixel_weights_delaunay_from( class DelaunayInterface: def __init__( - self, points, simplices, mappings, split_points, splitted_mappings, xp=np + self, + points, + simplices, + mappings, + split_points, + splitted_mappings, + xp=np, + dual_areas=None, ): + """ + Parameters + ---------- + dual_areas + The barycentric dual area of every mesh vertex (the sum of + ``triangle_area / 3`` over the triangles touching it), computed by + whichever Delaunay construction built this interface -- in-graph on + the JAX path, so it is safe to consume inside a ``jax.jit``. + ``None`` for interfaces built without them (e.g. the + natural-neighbour interface in ``sibson.py``), in which case the + mesh geometry falls back to a host-side computation. + """ self.points = points self.simplices = simplices self.mappings = mappings self.split_points = split_points self.splitted_mappings = splitted_mappings + self.dual_areas = dual_areas self.xp = xp @@ -653,6 +708,20 @@ def __init__( xp=xp, ) + @property + def dual_areas(self): + """ + The barycentric dual area of every mesh vertex, taken from the Delaunay + construction that already ran for the mappings -- so it is free here, + and on the JAX path it is an in-graph value safe to use inside a + ``jax.jit``. + + Subclasses that build no Delaunay tables (the kNN interpolators) + override this with ``None``, which makes the mesh geometry compute the + dual areas host-side only if something actually asks for them. + """ + return self.delaunay.dual_areas + @cached_property def mesh_geometry(self): @@ -664,6 +733,7 @@ def mesh_geometry(self): mesh=self.mesh, mesh_grid=self.mesh_grid, data_grid=self.data_grid, + dual_areas=self.dual_areas, xp=self._xp, ) @@ -700,7 +770,7 @@ def delaunay(self) -> "scipy.spatial.Delaunay": import jax.numpy as jnp - points, simplices, mappings, split_points, splitted_mappings = ( + points, simplices, mappings, split_points, splitted_mappings, areas = ( jax_delaunay( points=self.mesh_grid_xy, query_points=self.data_grid.over_sampled.array, @@ -710,7 +780,7 @@ def delaunay(self) -> "scipy.spatial.Delaunay": else: - points, simplices, mappings, split_points, splitted_mappings = ( + points, simplices, mappings, split_points, splitted_mappings, areas = ( scipy_delaunay( points_np=self.mesh_grid_xy, query_points_np=self.data_grid.over_sampled.array, @@ -724,14 +794,14 @@ def delaunay(self) -> "scipy.spatial.Delaunay": import jax.numpy as jnp - points, simplices, mappings = jax_delaunay_matern( + points, simplices, mappings, areas = jax_delaunay_matern( points=self.mesh_grid_xy, query_points=self.data_grid.over_sampled.array, ) else: - points, simplices, mappings = scipy_delaunay_matern( + points, simplices, mappings, areas = scipy_delaunay_matern( points_np=self.mesh_grid_xy, query_points_np=self.data_grid.over_sampled.array, ) @@ -746,6 +816,7 @@ def delaunay(self) -> "scipy.spatial.Delaunay": split_points=split_points, splitted_mappings=splitted_mappings, xp=self._xp, + dual_areas=areas, ) @property diff --git a/autoarray/inversion/mesh/interpolator/knn.py b/autoarray/inversion/mesh/interpolator/knn.py index 1b4774829..313bec62a 100644 --- a/autoarray/inversion/mesh/interpolator/knn.py +++ b/autoarray/inversion/mesh/interpolator/knn.py @@ -145,6 +145,17 @@ def kernel_interpolate_points(points, query_chunk, values, k, radius_scale): class InterpolatorKNearestNeighbor(InterpolatorDelaunay): + @property + def dual_areas(self): + """ + The kNN interpolators build no Delaunay tables, so there are no + in-graph dual areas to hand the mesh geometry. Returning ``None`` + keeps ``MeshGeometryDelaunay`` on its host-side fallback (which + triangulates only if the areas are actually asked for) rather than + forcing a qhull call here. + """ + return None + @cached_property def _mappings_sizes_weights(self): diff --git a/autoarray/inversion/mesh/mesh_geometry/delaunay.py b/autoarray/inversion/mesh/mesh_geometry/delaunay.py index ea16d1546..4bc332994 100644 --- a/autoarray/inversion/mesh/mesh_geometry/delaunay.py +++ b/autoarray/inversion/mesh/mesh_geometry/delaunay.py @@ -93,6 +93,39 @@ def voronoi_areas_numpy(points, qhull_options="Qbb Qc Qx Qm Q12 Pp"): class MeshGeometryDelaunay(AbstractMeshGeometry): + def __init__( + self, + mesh, + mesh_grid, + data_grid, + dual_areas=None, + xp=np, + **kwargs, + ): + """ + The geometry of a Delaunay triangulation / Voronoi mesh. + + Parameters + ---------- + dual_areas + The barycentric dual area of every mesh vertex, as computed by the + interpolator that built this geometry (see + `autoarray.inversion.mesh.interpolator.delaunay`). These are the + quadrature weights `areas_for_magnification` returns. When the + geometry is built standalone (e.g. in the tests) this is `None` and + the dual areas are computed host-side on demand instead. + """ + + super().__init__( + mesh=mesh, + mesh_grid=mesh_grid, + data_grid=data_grid, + xp=xp, + **kwargs, + ) + + self.dual_areas = dual_areas + @cached_property def mesh_grid_xy(self): """ @@ -194,23 +227,47 @@ def voronoi_areas(self): @property def areas_for_magnification(self) -> np.ndarray: """ - Returns the Voronoi cell area of every pixel in the mesh, as computed by `voronoi_areas_numpy` (a shoelace - sum over the `scipy.spatial.Voronoi` cell of each mesh point). - - Only cells that are **unbounded** in the Voronoi diagram (those `voronoi_areas_numpy` flags with the `-1` - sentinel, because their region runs to infinity and has no finite area) are set to zero. Cells that are - bounded but sit at the edge of the mesh are **kept at full size**, even though they can be far larger than - the interior cells -- an order of magnitude is routine, because a boundary cell extends out to the - circumcentres of the outermost triangles rather than being clipped to the mesh. - - These Voronoi areas are **not** the barycentric dual areas used by the Delaunay interpolator - (`barycentric_dual_area_from` in `autoarray.inversion.mesh.interpolator.delaunay`), which assign each vertex - the sum of `triangle_area / 3` over the triangles touching it. The dual areas tile the convex hull of the - mesh exactly and integrate the piecewise-linear reconstruction exactly; these Voronoi areas do neither, so - `sum(reconstruction * areas_for_magnification)` is not the integral of the reconstructed source. + Returns the **barycentric dual area** of every pixel in the mesh: the sum of `triangle_area / 3` over the + Delaunay triangles that touch that vertex. + + These are the exact quadrature weights of the mesh's piecewise-linear (barycentric) interpolant. The + Delaunay mapper reconstructs the source as the linear interpolant through the vertex values `s_i`, and for + that interpolant + + integral over the hull of f = sum_i s_i * dual_i + + holds exactly, not approximately -- the dual areas tile the convex hull of the mesh exactly (they sum to the + hull area), so `sum(reconstruction * areas_for_magnification)` is the integral of the reconstructed source + over the mesh. + + The value is taken from the interpolator, which computes it **in-graph** as part of the same Delaunay + construction that builds the mappings (`scipy_delaunay` / `jax_delaunay` and their Matérn variants). It is + therefore free here, and on the JAX path it is a traced value -- so this property is safe to evaluate inside + a `jax.jit` (e.g. PyAutoLens' per-sample latent evaluation). When the geometry is constructed standalone, + without an interpolator (`dual_areas is None`), the dual areas are computed host-side from a + `scipy.spatial.Delaunay` of the mesh grid, which gives the identical answer. + + This is **not** `voronoi_areas`, which sums the shoelace area of each point's `scipy.spatial.Voronoi` cell. + Those are correct arithmetic but the wrong quantity for this use: a Voronoi cell that is bounded but sits at + the edge of the mesh extends out to the circumcentres of the outermost triangles instead of being clipped to + the hull, so it can be orders of magnitude larger than that vertex's dual area. Weighting the reconstruction + by them biased magnification by -13% to -99% across the audited configurations (PyAutoArray#522); the fix is + PyAutoArray#524. `voronoi_areas` itself is unchanged and remains available for the geometric uses that + genuinely want a Voronoi cell. """ - areas = self.voronoi_areas + if self.dual_areas is not None: + return self.dual_areas - areas[areas == -1] = 0.0 + import scipy.spatial + + from autoarray.inversion.mesh.interpolator.delaunay import ( + barycentric_dual_area_from, + ) - return areas + mesh_grid_xy = np.asarray(self.mesh_grid_xy) + + return barycentric_dual_area_from( + mesh_grid_xy, + scipy.spatial.Delaunay(mesh_grid_xy).simplices, + xp=np, + ) diff --git a/autoarray/plot/inversion.py b/autoarray/plot/inversion.py index a99f465b6..1078b91ba 100644 --- a/autoarray/plot/inversion.py +++ b/autoarray/plot/inversion.py @@ -269,8 +269,20 @@ def _plot_rectangular( def _plot_delaunay(ax, pixel_values, mapper, norm, colormap, extent, is_subplot=False): """Render a Delaunay or KNN pixelization reconstruction onto *ax*. - Uses ``ax.tripcolor`` with Gouraud shading so that the reconstructed - flux is interpolated smoothly across the triangulated source-plane mesh. + Uses ``ax.tripcolor`` with its default *flat* shading: each triangle is + painted with the mean of its three vertex values. That is not the + piecewise-linear reconstruction pointwise, but its area-weighted integral + is identical to it by linearity -- the mean of the three vertices times the + triangle area is exactly the integral of the linear interpolant over that + triangle, which is the same sum the barycentric dual areas + (``mesh_geometry.areas_for_magnification``) accumulate per vertex. + + The triangulation is taken from the mapper's own Delaunay tables + (``mapper.interpolator.delaunay.simplices``, with the ``-1`` padded rows + dropped) so the figure shows the mesh the inversion actually used rather + than a triangulation matplotlib re-derives. Mappers with no such tables + (kNN, mocks) fall back to matplotlib's own triangulation. + A colorbar is attached after rendering. Parameters @@ -316,7 +328,22 @@ def _plot_delaunay(ax, pixel_values, mapper, norm, colormap, extent, is_subplot= else: vals = pixel_values - tc = ax.tripcolor(x, y, vals, cmap=colormap, norm=norm) + # The mesh grid is (y,x) and tripcolor is given x=column 1, y=column 0; the + # simplices index the same points in the same order, so they pass through + # unchanged. + triangles = None + + try: + simplices = np.asarray(mapper.interpolator.delaunay.simplices) + if simplices.ndim == 2 and simplices.shape[1] == 3: + triangles = simplices[simplices[:, 0] >= 0] + except AttributeError: + triangles = None + + if triangles is not None and triangles.shape[0] > 0: + tc = ax.tripcolor(x, y, triangles, vals, cmap=colormap, norm=norm) + else: + tc = ax.tripcolor(x, y, vals, cmap=colormap, norm=norm) from autoarray.plot.utils import _apply_colorbar _apply_colorbar(tc, ax, is_subplot=is_subplot) diff --git a/test_autoarray/inversion/pixelization/interpolator/test_delaunay.py b/test_autoarray/inversion/pixelization/interpolator/test_delaunay.py index 03928948a..f09324976 100644 --- a/test_autoarray/inversion/pixelization/interpolator/test_delaunay.py +++ b/test_autoarray/inversion/pixelization/interpolator/test_delaunay.py @@ -75,11 +75,15 @@ def test__scipy_delaunay__split(grid_2d_sub_1_7x7): # * `barycentric_dual_area_from` (here) -- sum of triangle_area / 3 over the # triangles touching each vertex. These tile the convex hull exactly and # integrate the piecewise-linear interpolant exactly. -# * `MeshGeometryDelaunay.areas_for_magnification` -- scipy Voronoi cell -# areas with only the *unbounded* cells zeroed. These do NOT tile the hull -# and do NOT integrate the interpolant. +# * `MeshGeometryDelaunay.areas_for_magnification` -- since PyAutoArray#524 +# these ARE the dual areas (taken straight from the interpolator, which +# computes them in-graph). They used to be scipy Voronoi cell areas with +# only the *unbounded* cells zeroed, which do NOT tile the hull and do NOT +# integrate the interpolant; that quantity is still available, unchanged, +# as `voronoi_areas` / `voronoi_areas_numpy`. # -# The tests below pin both facts, including the size of the divergence. +# The tests below pin both facts, including the size of the divergence between +# the dual areas and the Voronoi areas. # ---------------------------------------------------------------------------- # jax is an `[optional]` extra and is absent on the NumPy-only matrix env, so @@ -107,9 +111,13 @@ def test__barycentric_dual_area__single_triangle(): def test__barycentric_dual_area__sums_to_convex_hull_area(): """ The dual areas partition the convex hull exactly, so they sum to the hull - area. The Voronoi areas behind `areas_for_magnification` do not -- on this - configuration they overshoot the hull by ~29%, because the bounded boundary - cells extend well outside the hull and are kept. + area. The Voronoi areas -- which `areas_for_magnification` returned before + PyAutoArray#524 and which `voronoi_areas_numpy` still returns -- do not: on + this configuration they overshoot the hull by ~29%, because the bounded + boundary cells extend well outside the hull and are kept. + + Since #524 a `MeshGeometryDelaunay` built on the same points returns the + dual areas from `areas_for_magnification`, which is also asserted here. """ import scipy.spatial @@ -135,6 +143,18 @@ def test__barycentric_dual_area__sums_to_convex_hull_area(): ) assert voronoi.sum() != pytest.approx(hull_area, rel=1.0e-2) + # ... and the mesh geometry now hands out the dual areas, not those Voronoi + # areas. `mesh_grid_xy` may swap the (y,x) columns to (x,y); dual areas are + # invariant under that swap, so comparing against the dual areas of the same + # `points` array is exact. + mesh_geometry = aa.MeshGeometryDelaunay( + mesh=aa.mesh.Delaunay(pixels=40), + mesh_grid=aa.Grid2DIrregular(points), + data_grid=aa.Grid2D.uniform(shape_native=(7, 7), pixel_scales=1.0), + ) + + assert mesh_geometry.areas_for_magnification == pytest.approx(dual, rel=1.0e-10) + def test__linear_field_integral__dual_areas_exact(): """ @@ -155,9 +175,9 @@ def test__linear_field_integral__dual_areas_exact(): p1 = points[simplices[:, 1]] p2 = points[simplices[:, 2]] - cross = (p1[:, 0] - p0[:, 0]) * (p2[:, 1] - p0[:, 1]) - ( - p1[:, 1] - p0[:, 1] - ) * (p2[:, 0] - p0[:, 0]) + cross = (p1[:, 0] - p0[:, 0]) * (p2[:, 1] - p0[:, 1]) - (p1[:, 1] - p0[:, 1]) * ( + p2[:, 0] - p0[:, 0] + ) tri_area = 0.5 * np.abs(cross) @@ -174,9 +194,9 @@ def test__linear_field_integral__dual_areas_exact(): voronoi_integral = (f * voronoi).sum() - assert abs(voronoi_integral - exact) / exact > 0.01, ( - f"Voronoi-weighted integral {voronoi_integral} vs exact {exact}" - ) + assert ( + abs(voronoi_integral - exact) / exact > 0.01 + ), f"Voronoi-weighted integral {voronoi_integral} vs exact {exact}" @requires_jax @@ -186,10 +206,12 @@ def test__barycentric_dual_area__numpy_matches_jax_in_graph(): (a masked scatter-add over the padded simplices) rather than calling `barycentric_dual_area_from`. The two must agree. - The dual areas are not returned by either path; they enter it as the split - point weights (`areas_factor * sqrt(areas)`), so the split points -- which - both paths do return -- are the observable that pins them. Comparing them - exercises the real library path rather than a hand-rolled copy of it. + Since PyAutoArray#524 both paths return the dual areas as the sixth element + of their tuple (the mesh geometry consumes them as the magnification + quadrature weights), so they are compared directly. The split points -- + which carry the dual areas as `areas_factor * sqrt(areas)` -- are still + compared too, since they exercise the same values through the real library + path. """ import jax.numpy as jnp @@ -198,10 +220,13 @@ def test__barycentric_dual_area__numpy_matches_jax_in_graph(): points = rng.random((25, 2)) query_points = rng.random((30, 2)) - _, simplices_np, _, split_np, _ = scipy_delaunay(points, query_points, 0.5) - _, simplices_jx, _, split_jx, _ = jax_delaunay( + _, simplices_np, _, split_np, _, areas_np = scipy_delaunay( + points, query_points, 0.5 + ) + _, simplices_jx, _, split_jx, _, areas_jx = jax_delaunay( jnp.asarray(points), jnp.asarray(query_points), 0.5 ) assert np.array_equal(np.asarray(simplices_np), np.asarray(simplices_jx)) assert np.asarray(split_jx) == pytest.approx(np.asarray(split_np), abs=1.0e-10) + assert np.asarray(areas_jx) == pytest.approx(np.asarray(areas_np), abs=1.0e-10) diff --git a/test_autoarray/inversion/pixelization/interpolator/test_delaunay_nn.py b/test_autoarray/inversion/pixelization/interpolator/test_delaunay_nn.py index 696221ff1..5fe83ae53 100644 --- a/test_autoarray/inversion/pixelization/interpolator/test_delaunay_nn.py +++ b/test_autoarray/inversion/pixelization/interpolator/test_delaunay_nn.py @@ -88,7 +88,7 @@ def test__smooth_source_mapping_is_numerically_close_to_delaunay(): query_y, query_x = np.meshgrid(query_axis, query_axis, indexing="ij") query = np.stack([query_y.ravel(), query_x.ravel()], axis=1) - _, _, delaunay_mappings, _, _ = scipy_delaunay(points, query, areas_factor=0.5) + _, _, delaunay_mappings, _, _, _ = scipy_delaunay(points, query, areas_factor=0.5) delaunay_weights = pixel_weights_delaunay_from( data_grid=query, mesh_grid=points, diff --git a/test_autoarray/inversion/pixelization/mappers/test_delaunay.py b/test_autoarray/inversion/pixelization/mappers/test_delaunay.py index aed75a938..7bf1723d5 100644 --- a/test_autoarray/inversion/pixelization/mappers/test_delaunay.py +++ b/test_autoarray/inversion/pixelization/mappers/test_delaunay.py @@ -1,4 +1,7 @@ +import importlib.util + import numpy as np +import pytest import scipy.spatial import autoarray as aa @@ -95,3 +98,141 @@ def test__pix_indexes_for_sub_slim_index__delaunay_mesh__matches_util_and_expect assert ( mapper.pix_sizes_for_sub_slim_index == np.array([1, 1, 3, 1, 1, 1, 1, 1, 1]) ).all() + + +# ---------------------------------------------------------------------------- +# Magnification quadrature (PyAutoArray#524). +# +# `mesh_geometry.areas_for_magnification` returns the barycentric dual areas, +# which are the exact quadrature weights of the mapper's barycentric-linear +# reconstruction. The identity test below pins that against the mapping matrix +# itself -- the mapper's own definition of how flux is spread over the mesh. +# ---------------------------------------------------------------------------- + +# jax is an `[optional]` extra and is absent on the NumPy-only matrix env, so +# the JAX parity test skips rather than fails there (same convention as +# test_knn_barycentric.py). +requires_jax = pytest.mark.skipif( + importlib.util.find_spec("jax") is None, + reason="requires jax (installed via the [optional] extras; absent on the NumPy-only matrix env)", +) + + +def _magnification_identity_setup(): + """A Delaunay mapper whose mesh hull is pinned to the data footprint. + + The data grid's pixel centres span +/- 0.95 at a pixel scale of 0.1, so its + footprint is exactly [-1, 1] ** 2 (area 4). The mesh's outer ring sits on + +/- 1.0, so every data pixel centre is inside the hull (no out-of-hull + nearest-vertex fallback) and the hull area is exactly the footprint area -- + the two integrals below therefore cover the same region. + + Interior mesh points are jittered by <= 0.3 of the lattice spacing with a + seeded rng, so the triangulation is irregular rather than a regular lattice. + """ + data_grid = aa.Grid2D.uniform( + shape_native=(20, 20), pixel_scales=0.1, over_sample_size=1 + ) + + n = 8 + axis = np.linspace(-1.0, 1.0, n) + spacing = axis[1] - axis[0] + + mesh_y, mesh_x = np.meshgrid(axis, axis, indexing="ij") + mesh_points = np.stack([mesh_y.ravel(), mesh_x.ravel()], axis=1) + + rng = np.random.default_rng(42) + + interior = (np.abs(mesh_points[:, 0]) < 1.0 - 1.0e-9) & ( + np.abs(mesh_points[:, 1]) < 1.0 - 1.0e-9 + ) + mesh_points[interior] += ( + rng.uniform(-0.3, 0.3, size=(int(interior.sum()), 2)) * spacing + ) + + reconstruction = rng.uniform(0.5, 2.0, size=mesh_points.shape[0]) + + return data_grid, mesh_points, reconstruction + + +def test__areas_for_magnification__integrates_reconstruction_like_mapping_matrix(): + """ + The mapping matrix and the dual areas are two routes to the same integral of + the reconstruction over the source plane: + + * `(mapping_matrix @ s).sum() * pixel_scale ** 2` -- a Riemann sum of the + interpolant over the data pixels, using the mapper's own weights. + * `(s * areas_for_magnification).sum()` -- the exact quadrature of the same + piecewise-linear interpolant over the mesh hull. + + With the hull pinned to the data footprint they agree to the half-pixel + border error (~3e-5 here). The Voronoi areas -- what + `areas_for_magnification` returned before PyAutoArray#524 -- get the same + integral wrong by ~27%, which is asserted too so the flipped semantics are + documented by a number. + """ + import scipy.spatial + + data_grid, mesh_points, reconstruction = _magnification_identity_setup() + + mesh = aa.mesh.Delaunay(pixels=mesh_points.shape[0]) + + interpolator = mesh.interpolator_from( + source_plane_data_grid=data_grid, + source_plane_mesh_grid=aa.Grid2DIrregular(mesh_points), + ) + + mapper = aa.Mapper(interpolator=interpolator) + + # no data pixel took the out-of-hull nearest-vertex fallback + assert (np.asarray(mapper.interpolator.delaunay.mappings)[:, 1] >= 0).all() + + areas = np.asarray(mapper.mesh_geometry.areas_for_magnification) + + hull_area = scipy.spatial.ConvexHull(mesh_points).volume + + assert areas.sum() == pytest.approx(hull_area, rel=1.0e-10) + assert hull_area == pytest.approx(4.0, rel=1.0e-10) + + f_map = (np.asarray(mapper.mapping_matrix) @ reconstruction).sum() * 0.1**2 + f_dual = (reconstruction * areas).sum() + + assert f_map == pytest.approx(f_dual, rel=1.0e-4) + + voronoi = np.asarray(mapper.mesh_geometry.voronoi_areas) + voronoi = np.where(voronoi == -1.0, 0.0, voronoi) + + f_voronoi = (reconstruction * voronoi).sum() + + assert abs(f_voronoi - f_map) / f_map > 1.0e-2, ( + f"Voronoi-weighted integral {f_voronoi} vs mapping-matrix integral " + f"{f_map}; the two area definitions are not interchangeable." + ) + + +@requires_jax +def test__areas_for_magnification__jax_matches_numpy(): + """ + PyAutoLens evaluates latents inside a per-sample `jax.jit`, so + `areas_for_magnification` must come from the interpolator's in-graph dual + areas rather than a host-side scipy call. Building the same mesh geometry on + both backends and comparing pins that. + """ + import jax.numpy as jnp + + data_grid, mesh_points, _ = _magnification_identity_setup() + + mesh = aa.mesh.Delaunay(pixels=mesh_points.shape[0]) + + numpy_areas = mesh.interpolator_from( + source_plane_data_grid=data_grid, + source_plane_mesh_grid=aa.Grid2DIrregular(mesh_points), + ).mesh_geometry.areas_for_magnification + + jax_areas = mesh.interpolator_from( + source_plane_data_grid=data_grid, + source_plane_mesh_grid=aa.Grid2DIrregular(mesh_points), + xp=jnp, + ).mesh_geometry.areas_for_magnification + + assert np.asarray(jax_areas) == pytest.approx(np.asarray(numpy_areas), abs=1.0e-10) diff --git a/test_autoarray/inversion/pixelization/mesh_geometry/test_delaunay.py b/test_autoarray/inversion/pixelization/mesh_geometry/test_delaunay.py index 1c1a14d0e..b5d553d34 100644 --- a/test_autoarray/inversion/pixelization/mesh_geometry/test_delaunay.py +++ b/test_autoarray/inversion/pixelization/mesh_geometry/test_delaunay.py @@ -59,13 +59,21 @@ def test__voronoi_areas_via_delaunay_from(grid_2d_sub_1_7x7): def test__areas_for_magnification__uniform_lattice(grid_2d_sub_1_7x7): """ - Known-answer test: on an n x n unit-spacing lattice every *interior* Voronoi - cell is the unit square around its point (area exactly 1.0), and every point - on the convex hull has an unbounded Voronoi region, which - `areas_for_magnification` zeroes. - - The total is therefore (n - 2) ** 2, not n ** 2 -- the summed Delaunay - "source-plane area" is strictly smaller than the region the mesh covers. + Known-answer test on an n x n unit-spacing lattice, FLIPPED by PyAutoArray#524. + + `areas_for_magnification` now returns the *barycentric dual areas* (the sum of + triangle_area / 3 over the triangles touching each vertex), which tile the convex + hull of the mesh exactly. The total is therefore the hull area (n - 1) ** 2, and + every boundary point carries a real, non-zero area -- where the old Voronoi-based + definition zeroed the whole hull (every hull point has an unbounded Voronoi cell) + and summed to (n - 2) ** 2. + + Individual cell values are deliberately not pinned: on a perfectly regular lattice + the dual area of a vertex depends on which diagonal qhull picks for each square + (interior values alternate 2/3 and 4/3, corners are 1/6 or 1/3), so only the + diagonal-independent facts are asserted here. The elementwise identity against + `barycentric_dual_area_from` is pinned by + `test__areas_for_magnification__equals_barycentric_dual_area` below. (scipy's `Qbb Qc Qx Qm Q12 Pp` options handle the perfectly regular lattice without degenerate output, so no jitter of the interior points is needed.) @@ -85,30 +93,46 @@ def test__areas_for_magnification__uniform_lattice(grid_2d_sub_1_7x7): areas = mesh.areas_for_magnification.reshape(n, n) - interior = areas[1:-1, 1:-1] - assert interior == pytest.approx(np.ones((n - 2, n - 2)), 1.0e-8) + # the dual areas tile the hull exactly + assert areas.sum() == pytest.approx(float((n - 1) ** 2), 1.0e-8) + + # ... which is NOT the old, Voronoi-based total + assert areas.sum() != pytest.approx(float((n - 2) ** 2), 1.0e-2) + + # every cell, boundary included, now carries a real area + assert (areas > 0.0).all() - assert areas[0, :] == pytest.approx(np.zeros(n), 1.0e-8) - assert areas[-1, :] == pytest.approx(np.zeros(n), 1.0e-8) - assert areas[:, 0] == pytest.approx(np.zeros(n), 1.0e-8) - assert areas[:, -1] == pytest.approx(np.zeros(n), 1.0e-8) + # every dual area is a sum of (unit-square-half) / 3 terms, so it is bounded by + # the 1 (corner) to 8 (interior) triangles a lattice vertex can touch + assert areas.min() >= 0.5 / 3.0 - 1.0e-12 + assert areas.max() <= 8.0 * 0.5 / 3.0 + 1.0e-12 - assert areas.sum() == pytest.approx(float((n - 2) ** 2), 1.0e-8) + # the hull boundary is exactly the ring the old definition threw away + boundary = areas.sum() - areas[1:-1, 1:-1].sum() + assert boundary > 0.0 -def test__areas_for_magnification__bounded_boundary_cells_are_kept(grid_2d_sub_1_7x7): +def test__areas_for_magnification__equals_barycentric_dual_area(grid_2d_sub_1_7x7): """ - Pins the CURRENT semantics of `areas_for_magnification`: only cells whose - Voronoi region is *unbounded* (the `-1` sentinel) are zeroed. A cell that is - bounded but sits at the edge of the mesh is kept at full size, even when it - is an order of magnitude larger than the interior cells -- here index 3 keeps - an area of ~29.8 against index 1's ~1.4. - - This is the bias candidate flagged by the euclid-dr1-prep phase 8 audit: a - magnification denominator built from these areas is inflated by the huge - bounded boundary cells. A later fix should flip this test *deliberately*, - not silently. + FLIPPED by PyAutoArray#524 (was + `test__areas_for_magnification__bounded_boundary_cells_are_kept`, which pinned the + Voronoi semantics the euclid-dr1-prep phase 8 audit flagged as the bias source and + which a later fix was told to flip deliberately). + + `areas_for_magnification` is now the barycentric dual area of every vertex, so on + this 6-point configuration it equals `barycentric_dual_area_from` elementwise and + sums to the convex-hull area. + + In particular index 3 is no longer the ~29.8 bounded-but-huge Voronoi cell, and + index 4 -- whose Voronoi region is unbounded, and which the old definition zeroed + outright -- now carries a real area. """ + import scipy.spatial + + from autoarray.inversion.mesh.interpolator.delaunay import ( + barycentric_dual_area_from, + ) + mesh_grid = aa.Grid2DIrregular( [[0.0, 0.0], [1.1, 0.6], [2.1, 0.1], [0.4, 1.1], [1.1, 7.1], [2.1, 1.1]] ) @@ -121,11 +145,23 @@ def test__areas_for_magnification__bounded_boundary_cells_are_kept(grid_2d_sub_1 areas = mesh.areas_for_magnification - assert areas[1] == pytest.approx(1.39137102, 1.0e-4) - # bounded, huge, and kept: - assert areas[3] == pytest.approx(29.836324, 1.0e-4) - # unbounded (-1 in `voronoi_areas`), and therefore zeroed: - assert areas[4] == pytest.approx(0.0, 1.0e-8) + mesh_grid_xy = np.asarray(mesh.mesh_grid_xy) + + expected = barycentric_dual_area_from( + mesh_grid_xy, + scipy.spatial.Delaunay(mesh_grid_xy).simplices, + xp=np, + ) + + assert areas == pytest.approx(expected, rel=1.0e-10) + + hull_area = scipy.spatial.ConvexHull(mesh_grid_xy).volume + + assert areas.sum() == pytest.approx(hull_area, rel=1.0e-10) + + # the two indexes the old, Voronoi-based definition got wrong + assert areas[3] != pytest.approx(29.836324, 1.0e-2) + assert areas[4] > 0.0 def test__areas_for_magnification__repeat_calls_agree(grid_2d_sub_1_7x7):