diff --git a/autoarray/inversion/mesh/mesh_geometry/delaunay.py b/autoarray/inversion/mesh/mesh_geometry/delaunay.py index 093fa6172..ea16d1546 100644 --- a/autoarray/inversion/mesh/mesh_geometry/delaunay.py +++ b/autoarray/inversion/mesh/mesh_geometry/delaunay.py @@ -194,11 +194,20 @@ def voronoi_areas(self): @property def areas_for_magnification(self) -> np.ndarray: """ - Returns the area of every Voronoi pixel in the Voronoi mesh. - - Pixels at boundaries can sometimes have large unrealistic areas, which can impact the magnification - calculation. This method therefore sets their areas to zero so they do not impact the magnification - calculation. + 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. """ areas = self.voronoi_areas diff --git a/test_autoarray/inversion/pixelization/interpolator/test_delaunay.py b/test_autoarray/inversion/pixelization/interpolator/test_delaunay.py index 4d3bbe559..03928948a 100644 --- a/test_autoarray/inversion/pixelization/interpolator/test_delaunay.py +++ b/test_autoarray/inversion/pixelization/interpolator/test_delaunay.py @@ -1,7 +1,15 @@ +import importlib.util + import numpy as np import pytest import autoarray as aa +from autoarray.inversion.mesh.interpolator.delaunay import ( + barycentric_dual_area_from, + jax_delaunay, + scipy_delaunay, +) +from autoarray.inversion.mesh.mesh_geometry.delaunay import voronoi_areas_numpy def test__scipy_delaunay__simplices(grid_2d_sub_1_7x7): @@ -57,3 +65,143 @@ def test__scipy_delaunay__split(grid_2d_sub_1_7x7): assert mesh_grid.delaunay.splitted_mappings[-1, :] == pytest.approx( [5, 1, 2], 1.0e-4 ) + + +# ---------------------------------------------------------------------------- +# Barycentric dual areas (euclid-dr1-prep phase 8 audit). +# +# Two different "areas" exist for a Delaunay mesh on two different code paths: +# +# * `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. +# +# The tests below pin both facts, including the size of the divergence. +# ---------------------------------------------------------------------------- + +# jax is an `[optional]` extra and is absent on the NumPy-only matrix env, so +# the in-graph 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 test__barycentric_dual_area__single_triangle(): + + points = np.array([[0.0, 0.0], [0.0, 4.0], [3.0, 0.0]]) + simplices = np.array([[0, 1, 2]]) + + area = 0.5 * 3.0 * 4.0 + + dual = barycentric_dual_area_from(points, simplices, xp=np) + + assert dual == pytest.approx(np.full(3, area / 3.0), 1.0e-10) + assert dual.sum() == pytest.approx(area, 1.0e-10) + + +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. + """ + import scipy.spatial + + points = np.random.default_rng(1).random((40, 2)) + + simplices = scipy.spatial.Delaunay(points).simplices + + dual = barycentric_dual_area_from(points, simplices, xp=np) + + hull_area = scipy.spatial.ConvexHull(points).volume + + assert dual.sum() == pytest.approx(hull_area, rel=1.0e-10) + + voronoi = voronoi_areas_numpy(points) + voronoi = np.where(voronoi == -1.0, 0.0, voronoi) + + ratio = voronoi.sum() / hull_area + + assert ratio == pytest.approx(1.2868, 1.0e-3), ( + f"Voronoi areas (unbounded cells zeroed) sum to {voronoi.sum()} against a " + f"convex-hull area of {hull_area} (ratio {ratio}); the two area " + f"definitions are not interchangeable." + ) + assert voronoi.sum() != pytest.approx(hull_area, rel=1.0e-2) + + +def test__linear_field_integral__dual_areas_exact(): + """ + For a field that is linear over the mesh, the piecewise-linear interpolant + is the field itself, so its integral over the hull is exactly + `sum(f_i * dual_area_i)`. The Voronoi areas get the same integral wrong by + tens of percent. + """ + import scipy.spatial + + points = np.random.default_rng(1).random((40, 2)) + + simplices = scipy.spatial.Delaunay(points).simplices + + f = 0.3 + 0.7 * points[:, 0] - 0.2 * points[:, 1] + + p0 = points[simplices[:, 0]] + 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]) + + tri_area = 0.5 * np.abs(cross) + + # exact: on each triangle a linear field integrates to (mean of its three + # vertex values) x (triangle area) + exact = (f[simplices].mean(axis=1) * tri_area).sum() + + dual = barycentric_dual_area_from(points, simplices, xp=np) + + assert (f * dual).sum() == pytest.approx(exact, rel=1.0e-10) + + voronoi = voronoi_areas_numpy(points) + voronoi = np.where(voronoi == -1.0, 0.0, voronoi) + + voronoi_integral = (f * voronoi).sum() + + assert abs(voronoi_integral - exact) / exact > 0.01, ( + f"Voronoi-weighted integral {voronoi_integral} vs exact {exact}" + ) + + +@requires_jax +def test__barycentric_dual_area__numpy_matches_jax_in_graph(): + """ + `jax_delaunay` carries its own in-graph copy of the dual-area computation + (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. + """ + import jax.numpy as jnp + + rng = np.random.default_rng(3) + + 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( + 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) diff --git a/test_autoarray/inversion/pixelization/mesh_geometry/test_delaunay.py b/test_autoarray/inversion/pixelization/mesh_geometry/test_delaunay.py index ca00035ca..1c1a14d0e 100644 --- a/test_autoarray/inversion/pixelization/mesh_geometry/test_delaunay.py +++ b/test_autoarray/inversion/pixelization/mesh_geometry/test_delaunay.py @@ -55,3 +55,100 @@ def test__voronoi_areas_via_delaunay_from(grid_2d_sub_1_7x7): assert voronoi_areas[1] == pytest.approx(1.39137102, 1.0e-4) assert voronoi_areas[3] == pytest.approx(29.836324, 1.0e-4) assert voronoi_areas[4] == pytest.approx(-1.0, 1.0e-4) + + +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. + + (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.) + """ + n = 5 + + y, x = np.meshgrid( + np.arange(n, dtype=float), np.arange(n, dtype=float), indexing="ij" + ) + points = np.stack([y.ravel(), x.ravel()], axis=1) + + mesh = aa.MeshGeometryDelaunay( + mesh=aa.mesh.Delaunay(pixels=n * n), + mesh_grid=aa.Grid2DIrregular(points), + data_grid=grid_2d_sub_1_7x7.over_sampled, + ) + + 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) + + 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) + + assert areas.sum() == pytest.approx(float((n - 2) ** 2), 1.0e-8) + + +def test__areas_for_magnification__bounded_boundary_cells_are_kept(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. + """ + 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]] + ) + + mesh = aa.MeshGeometryDelaunay( + mesh=aa.mesh.Delaunay(pixels=6), + mesh_grid=mesh_grid, + data_grid=grid_2d_sub_1_7x7.over_sampled, + ) + + 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) + + +def test__areas_for_magnification__repeat_calls_agree(grid_2d_sub_1_7x7): + """ + `areas_for_magnification` mutates the array it gets from `voronoi_areas` + (`areas[areas == -1] = 0.0`). `voronoi_areas` is an uncached property that + recomputes each call, so the `-1` sentinel must never be written back into + any state shared between calls. + """ + 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]] + ) + + mesh = aa.MeshGeometryDelaunay( + mesh=aa.mesh.Delaunay(pixels=6), + mesh_grid=mesh_grid, + data_grid=grid_2d_sub_1_7x7.over_sampled, + ) + + first = mesh.areas_for_magnification + second = mesh.areas_for_magnification + + assert first == pytest.approx(second, 1.0e-8) + + # the sentinel is still there afterwards -- nothing was written back + assert mesh.voronoi_areas[4] == pytest.approx(-1.0, 1.0e-8)