Merge pull request #523 from PyAutoLabs/feature/delaunay-area-magnification-audit

test: Delaunay pixel-area known-answer tests + areas_for_magnification docstring (#522)

Author: Jammy2211

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
 

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)

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)