From e62fc29e0f9a59e3a1cebe7252f89d7b08c818ca Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 28 Oct 2026 20:42:51 +0000 Subject: [PATCH 1/6] fix: RectangularUniform areas_for_magnification use its own uniform cells The uniform mesh's geometry routed edges_transformed / areas through the adaptive meshes' kernel-CDF guard-ring convention, which does not describe the uniform mapper (bilinear between overlay pixel centres, no CDF, no guard ring). Magnifications from it were biased by -77% to +259% (PyAutoArray#659). The uniform interpolator now builds its geometry with transform="uniform", whose edges are the uniform cells centred on every overlay node (the same cells its imshow plot draws). xp-aware; jit/grad checked. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_016nhQ8NkVvJhhqNSibeGcBQ --- .../mesh/interpolator/rectangular_uniform.py | 1 + .../inversion/mesh/mesh_geometry/abstract.py | 7 +++-- .../mesh/mesh_geometry/rectangular.py | 30 +++++++++++++++++++ 3 files changed, 35 insertions(+), 3 deletions(-) diff --git a/autoarray/inversion/mesh/interpolator/rectangular_uniform.py b/autoarray/inversion/mesh/interpolator/rectangular_uniform.py index 009e5392b..a7abd65c2 100644 --- a/autoarray/inversion/mesh/interpolator/rectangular_uniform.py +++ b/autoarray/inversion/mesh/interpolator/rectangular_uniform.py @@ -158,6 +158,7 @@ def mesh_geometry(self): mesh_grid=self.mesh_grid, data_grid=self.data_grid, mesh_weight_map=self.mesh_weight_map, + transform="uniform", xp=self._xp, ) diff --git a/autoarray/inversion/mesh/mesh_geometry/abstract.py b/autoarray/inversion/mesh/mesh_geometry/abstract.py index 958f27a02..cd7665c44 100644 --- a/autoarray/inversion/mesh/mesh_geometry/abstract.py +++ b/autoarray/inversion/mesh/mesh_geometry/abstract.py @@ -24,9 +24,10 @@ def __init__( self.kernel_bandwidth = kernel_bandwidth self.kernel_knots = kernel_knots # The rectangular lattice transform the areas / edges route through: - # "kernel" (the RTU meshes and the uniform mesh's fallback) or "rank" - # (the Bilinear meshes) — must match the mapper's transform so the - # geometry is consistent with where the mapper scatters its flux. + # "kernel" (the RTU meshes), "rank" (the Bilinear meshes) or "uniform" + # (RectangularUniform: no CDF, cells centred on the overlay nodes) — + # must match the mapper's transform so the geometry is consistent + # with where the mapper scatters its flux. self.transform = transform self._use_jax = xp is not np diff --git a/autoarray/inversion/mesh/mesh_geometry/rectangular.py b/autoarray/inversion/mesh/mesh_geometry/rectangular.py index 1fb4775a6..0ca2227ab 100644 --- a/autoarray/inversion/mesh/mesh_geometry/rectangular.py +++ b/autoarray/inversion/mesh/mesh_geometry/rectangular.py @@ -474,7 +474,18 @@ def edges_transformed(self): mapper scatters its flux (issue #372). Guard-node edges fall outside the CDF domain and clamp to the data span in the inverse interp, squashing the border guard cells to the region they actually cover. + + ``RectangularUniform`` (``transform="uniform"``) is the exception: its + mapper interpolates bilinearly between the overlay grid's pixel + centres with no CDF and no guard ring, so its edges are the uniform + cell boundaries half a node spacing either side of every node (the + same cells its ``imshow`` plot draws). Routing it through the kernel + CDF and guard-ring convention instead biased its + ``areas_for_magnification`` by -77% to +259% across the audited + configurations (PyAutoArray#659). """ + if self.transform == "uniform": + return self._edges_uniform # Node-midpoint edges in unit CDF space (see docstring); the guard # denominator keeps the degenerate n <= 3 meshes finite (their @@ -515,6 +526,25 @@ def edges_transformed(self): xp=self._xp, ) + @property + def _edges_uniform(self): + """ + The cell edges of the uniform rectangular mesh, shape ``(n + 1, 2)`` + packing (y_edges descending, x_edges ascending): every overlay-grid + node sits at the centre of its cell, so the cells tile the overlay's + bounding box with one node spacing per cell. + """ + xp = self._xp + grid = self.mesh_grid.array + n_y, n_x = self.shape_native + y_max = xp.max(grid[:, 0]) + x_min = xp.min(grid[:, 1]) + dy = (y_max - xp.min(grid[:, 0])) / (n_y - 1) + dx = (xp.max(grid[:, 1]) - x_min) / (n_x - 1) + edges_y = y_max + 0.5 * dy - xp.arange(n_y + 1) * dy + edges_x = x_min - 0.5 * dx + xp.arange(n_x + 1) * dx + return xp.stack([edges_y, edges_x]).T + @property def neighbors(self) -> Neighbors: """ From 114204b4a332bb0e76a703e7b7cc18ba138351bc Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 28 Oct 2026 20:42:51 +0000 Subject: [PATCH 2/6] test: per-mesh areas_for_magnification correctness (identity-lens quadrature) For all nine meshes: areas sum to the analytic covered domain, the identity-lens magnification of a reconstructed Gaussian against the mapper's exact quadrature, and pinned guard-cell / hull-vertex / zeroed-ring semantics. Approximate meshes carry measured, documented tolerances (PyAutoArray#659). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_016nhQ8NkVvJhhqNSibeGcBQ --- .../test_areas_for_magnification.py | 288 ++++++++++++++++++ 1 file changed, 288 insertions(+) create mode 100644 test_autoarray/inversion/pixelization/mesh_geometry/test_areas_for_magnification.py 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 new file mode 100644 index 000000000..d8394bb97 --- /dev/null +++ b/test_autoarray/inversion/pixelization/mesh_geometry/test_areas_for_magnification.py @@ -0,0 +1,288 @@ +""" +Per-mesh correctness of `mesh_geometry.areas_for_magnification` (PyAutoArray#659). + +Every pixelized magnification (PyAutoLens `analysis/latent.py`, the workspace +`source_science.py` scripts) divides the lensed image-plane flux by +`sum(reconstruction * areas_for_magnification)`, so the areas must be the +quadrature weights of the reconstructed source surface. + +The tests use an identity lens: the source-plane data grid is a uniform image grid, +so every data sub-pixel covers a known source-plane area. For any reconstruction +`s` the mapper's own interpolant then integrates to + + Q_k = sum_j mapping_matrix_over_sampled[j, k] * sub_pixel_area + integral of the reconstructed source = sum_k s_k * Q_k + +and `sum_k s_k * A_k / sum_k s_k * Q_k` is the identity-lens magnification, which +must be 1. The checks are: + +(a) the areas sum to the analytic area of the domain the mesh covers — the data + grid's bounding box for the adaptive rectangular meshes (the CDF maps that box + onto the unit square), the overlay box for `RectangularUniform`, and the convex + hull of the vertices for the Delaunay family (barycentric dual areas tile it); +(b) for a smooth Gaussian source reconstructed on the mesh, the identity-lens + magnification is 1 and the integrated reconstruction matches the analytic + integral, within a per-mesh tolerance; +(c) the guard / boundary-cell semantics are pinned. + +Tolerances for meshes whose areas are an approximation rather than the exact +quadrature of their interpolant are set at the measured error (see +`MAGNIFICATION_TOLERANCE`) and documented there; they are not loosened to make a +test pass. +""" + +import numpy as np +import pytest +import scipy.spatial + +import autoarray as aa + +SHAPE = 24 +PIXEL_SCALE = 0.2 # image box is [-2.4, 2.4]^2 +SUB_SIZE = 2 +MESH_SHAPE = (12, 12) + +GAUSS_CENTRE = (0.2, -0.1) +GAUSS_SIGMA = 0.5 + + +# Absolute tolerance on |mu_identity - 1| (the mapper's exact quadrature of the +# reconstruction over its area-weighted sum) for every mesh, at this module's 24x24 +# data grid and 12x12-equivalent mesh, a deliberately coarse mesh (Gaussian sigma ~ +# 1.25 mesh spacings) on which approximation errors are largest. Values were measured +# with this exact configuration (PyAutoArray#659) and set ~1.5-2x above the +# measurement, quoted in brackets. +# +# Exact quadrature (tight): +# - Delaunay: barycentric dual areas are the exact quadrature of the linear +# interpolant (PyAutoArray#524) [-1.1e-4]. +# - RectangularUniform: uniform cells centred on the overlay nodes, the trapezoid +# weights of its bilinear interpolant [-1.2e-3]. Before PyAutoArray#659 its areas +# came from the adaptive meshes' kernel-CDF guard-ring geometry, which does not +# describe the uniform mapper at all, and this error was +55%. +# - RectangularRTUAdaptDensity: node-midpoint cells of a near-linear kernel CDF +# [+2.5e-3]. +# +# Documented approximations (flagged in PyAutoArray#659 for a ruling): +# - RectangularBilinearAdaptDensity / RectangularBilinearAdaptImage / +# RectangularRTUAdaptImage: `areas_for_magnification` is each node's +# node-midpoint cell pushed through the inverse CDF. That is not the integral of +# the node's bilinear hat once the CDF is non-linear across a cell (the rank +# CDF's staircase on tied coordinates; a steep adapt-image CDF). The exact +# per-axis hat quadrature brings these to ~1e-3 in a prototype +# [+3.1e-2, +2.4e-1, +8.0e-2; ~1-3% for all three at 20x20 and finer]. +# - 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]. +MAGNIFICATION_TOLERANCE = { + "RectangularUniform": 2.5e-3, + "RectangularBilinearAdaptDensity": 5.0e-2, + "RectangularBilinearAdaptImage": 3.5e-1, + "RectangularRTUAdaptDensity": 5.0e-3, + "RectangularRTUAdaptImage": 1.2e-1, + "Delaunay": 3.0e-4, + "DelaunayNN": 1.0e-2, + "KNearestNeighbor": 3.0e-2, + "KNNBarycentric": 7.0e-2, +} + +RECTANGULAR_ADAPTIVE = [ + "RectangularBilinearAdaptDensity", + "RectangularBilinearAdaptImage", + "RectangularRTUAdaptDensity", + "RectangularRTUAdaptImage", +] +DELAUNAY_FAMILY = ["Delaunay", "DelaunayNN", "KNearestNeighbor", "KNNBarycentric"] + + +def _gaussian(grid): + return np.exp( + -((grid[:, 0] - GAUSS_CENTRE[0]) ** 2 + (grid[:, 1] - GAUSS_CENTRE[1]) ** 2) + / (2.0 * GAUSS_SIGMA**2) + ) + + +def _vertices(): + """ + A jittered 12x12 lattice spanning the image box exactly (the outer ring is not + jittered), so the Delaunay hull is the whole image box. + """ + rng = np.random.default_rng(1) + t = np.linspace(0.0, 1.0, MESH_SHAPE[0]) + yy, xx = np.meshgrid(t, t, indexing="ij") + nodes = np.stack([yy.ravel(), xx.ravel()], axis=1) + interior = ((nodes > 0.0) & (nodes < 1.0)).all(axis=1) + jitter = rng.uniform(-0.3, 0.3, nodes.shape) / (MESH_SHAPE[0] - 1) + nodes = nodes + jitter * interior[:, None] + half = 0.5 * SHAPE * PIXEL_SCALE + return -half + 2.0 * half * nodes + + +def _mesh_from(name): + if name.startswith("Rectangular"): + return getattr(aa.mesh, name)(shape=MESH_SHAPE) + return getattr(aa.mesh, name)(pixels=MESH_SHAPE[0] * MESH_SHAPE[1]) + + +_CACHE = {} + + +def _setup(name): + """ + The mapper, its areas, its exact quadrature `Q` and a Gaussian reconstruction, + cached per mesh so the parametrized tests share one build. + """ + if name in _CACHE: + return _CACHE[name] + + grid = aa.Grid2D.uniform( + shape_native=(SHAPE, SHAPE), + pixel_scales=PIXEL_SCALE, + over_sample_size=SUB_SIZE, + ) + adapt_data = aa.Array2D(values=_gaussian(grid.array) + 0.01, mask=grid.mask) + + interpolator = _mesh_from(name).interpolator_from( + source_plane_data_grid=grid, + source_plane_mesh_grid=aa.Grid2DIrregular(_vertices()), + adapt_data=adapt_data, + ) + mapper = aa.Mapper(interpolator=interpolator) + + areas = np.asarray(mapper.mesh_geometry.areas_for_magnification) + + mapping_matrix = np.asarray(mapper.mapping_matrix_over_sampled) + sub_area = (PIXEL_SCALE / SUB_SIZE) ** 2 + quadrature = mapping_matrix.sum(axis=0) * sub_area + + # Noise-free least-squares reconstruction of the Gaussian at the data sub-pixels + # (a vanishing ridge only fixes the unmapped guard cells at zero). + data = _gaussian(np.asarray(grid.over_sampled.array)) + curvature = mapping_matrix.T @ mapping_matrix + ridge = 1.0e-8 * np.mean(np.diag(curvature)) * np.eye(curvature.shape[0]) + reconstruction = np.linalg.solve(curvature + ridge, mapping_matrix.T @ data) + + _CACHE[name] = (grid, mapper, areas, quadrature, reconstruction) + return _CACHE[name] + + +@pytest.mark.parametrize("name", list(MAGNIFICATION_TOLERANCE)) +def test__areas_sum_to_the_covered_domain(name): + grid, mapper, areas, _, _ = _setup(name) + + assert np.all(np.isfinite(areas)) + assert np.all(areas >= 0.0) + assert areas.shape == (mapper.params,) + + if name in DELAUNAY_FAMILY: + vertices = np.asarray(mapper.mesh_geometry.mesh_grid_xy) + domain = scipy.spatial.ConvexHull(vertices).volume + rel = 1.0e-10 + elif name == "RectangularUniform": + # the overlay grid spans the over-sampled data plus a 1e-8 buffer + domain = np.prod(np.ptp(np.asarray(grid.over_sampled.array), axis=0)) + rel = 1.0e-6 + else: + # the adaptive CDF maps the data grid's bounding box onto the unit square + domain = np.prod(np.ptp(np.asarray(grid.array), axis=0)) + rel = 1.0e-10 + + assert areas.sum() == pytest.approx(domain, rel=rel) + + +@pytest.mark.parametrize("name", list(MAGNIFICATION_TOLERANCE)) +def test__identity_lens_magnification_of_a_gaussian_source(name): + _, _, areas, quadrature, reconstruction = _setup(name) + + source_flux = np.sum(reconstruction * areas) + image_flux = np.sum(reconstruction * quadrature) + + tol = MAGNIFICATION_TOLERANCE[name] + + assert image_flux / source_flux == pytest.approx(1.0, abs=tol) + + # the integrated reconstruction against the analytic Gaussian integral (the + # source sits well inside the image box, so its truncation is negligible) + analytic = 2.0 * np.pi * GAUSS_SIGMA**2 + assert source_flux / analytic == pytest.approx(1.0, abs=tol) + + +@pytest.mark.parametrize("name", RECTANGULAR_ADAPTIVE) +def test__rectangular_adaptive__guard_cells_have_zero_area_and_no_flux(name): + """ + The adaptive rectangular mapper is node-based with a guard ring in unit CDF + space: rows 0 and 1 (above U_y = 1) and columns 0 and n - 1 (outside + [0, 1] in U_x) never receive mapping weight, and their node-midpoint edges clamp + to the data span so their area is exactly zero. Every other cell is live and + carries a positive area — including the bottom row (U_y = 0), which is live + although it sits in the `zeroed_pixels` edge ring. + """ + _, mapper, areas, quadrature, _ = _setup(name) + + n_y, n_x = MESH_SHAPE + mapped = quadrature.reshape(n_y, n_x) > 0.0 + area = areas.reshape(n_y, n_x) + + guard = np.zeros((n_y, n_x), dtype=bool) + guard[:2, :] = True + guard[:, 0] = True + guard[:, -1] = True + + np.testing.assert_array_equal(mapped, ~guard) + np.testing.assert_array_equal(area[guard], 0.0) + assert np.all(area[~guard] > 0.0) + + +def test__rectangular_uniform__every_cell_is_one_node_spacing(): + """ + `RectangularUniform` interpolates between the overlay grid's pixel centres with + no CDF and no guard ring: every cell, edge cells included, is one node spacing + on a side and receives mapping weight. + """ + _, mapper, areas, quadrature, _ = _setup("RectangularUniform") + + mesh_grid = np.asarray(mapper.mesh_geometry.mesh_grid.array) + spacing = np.ptp(mesh_grid, axis=0) / (np.array(MESH_SHAPE) - 1) + + np.testing.assert_allclose(areas, np.prod(spacing), rtol=1.0e-10) + assert np.all(quadrature > 0.0) + + edges = np.asarray(mapper.mesh_geometry.edges_transformed) + # y edges descend from the top of the overlay box, x edges ascend + np.testing.assert_allclose(np.diff(edges[:, 0]), -spacing[0], rtol=1.0e-10) + np.testing.assert_allclose(np.diff(edges[:, 1]), spacing[1], rtol=1.0e-10) + assert edges[0, 0] == pytest.approx(mesh_grid[:, 0].max() + 0.5 * spacing[0]) + assert edges[0, 1] == pytest.approx(mesh_grid[:, 1].min() - 0.5 * spacing[1]) + + +@pytest.mark.parametrize("name", DELAUNAY_FAMILY) +def test__delaunay_family__hull_vertices_keep_their_area_and_zeroing_does_not_change_it( + name, +): + """ + Every vertex, hull vertices included, carries a positive barycentric dual area + (the old Voronoi definition zeroed unbounded hull cells). The areas are a + property of the vertices alone: zeroing the last vertices (`zeroed_pixels`, the + appended edge ring) fixes their reconstruction to zero but does not change any + area. + """ + grid, mapper, areas, _, _ = _setup(name) + + assert np.all(areas > 0.0) + + mesh = getattr(aa.mesh, name)( + pixels=MESH_SHAPE[0] * MESH_SHAPE[1] - 10, zeroed_pixels=10 + ) + zeroed_mapper = aa.Mapper( + interpolator=mesh.interpolator_from( + source_plane_data_grid=grid, + source_plane_mesh_grid=aa.Grid2DIrregular(_vertices()), + ) + ) + + np.testing.assert_allclose( + np.asarray(zeroed_mapper.mesh_geometry.areas_for_magnification), + areas, + rtol=1.0e-12, + ) From 850f7b3ce8582e3decab720d85e3abe8a71435b3 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 28 Oct 2026 20:43:15 +0000 Subject: [PATCH 3/6] test: classify RectangularRTUAdaptDensity areas as a documented approximation Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_016nhQ8NkVvJhhqNSibeGcBQ --- .../mesh_geometry/test_areas_for_magnification.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) 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 d8394bb97..2d8ce3caa 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 @@ -60,10 +60,12 @@ # weights of its bilinear interpolant [-1.2e-3]. Before PyAutoArray#659 its areas # came from the adaptive meshes' kernel-CDF guard-ring geometry, which does not # describe the uniform mapper at all, and this error was +55%. -# - RectangularRTUAdaptDensity: node-midpoint cells of a near-linear kernel CDF -# [+2.5e-3]. # # Documented approximations (flagged in PyAutoArray#659 for a ruling): +# - RectangularRTUAdaptDensity: node-midpoint cells, accurate here because a +# density-adapted kernel CDF of an unlensed grid is near-linear [+2.5e-3]; the +# lensed workspace check (autolens_workspace_test +# scripts/imaging/mesh_magnification.py) measures -1.8%. # - RectangularBilinearAdaptDensity / RectangularBilinearAdaptImage / # RectangularRTUAdaptImage: `areas_for_magnification` is each node's # node-midpoint cell pushed through the inverse CDF. That is not the integral of From 16f22be6d13918fcbfd527c3d23182d01abf6af2 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 28 Oct 2026 20:49:13 +0000 Subject: [PATCH 4/6] docs: accurate RectangularUniform area/interpolation wording (review) Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_016nhQ8NkVvJhhqNSibeGcBQ --- .../mesh/mesh/rectangular_uniform.py | 7 +++--- .../mesh/mesh_geometry/rectangular.py | 24 ++++++++++++------- .../test_areas_for_magnification.py | 14 +++++++---- 3 files changed, 29 insertions(+), 16 deletions(-) diff --git a/autoarray/inversion/mesh/mesh/rectangular_uniform.py b/autoarray/inversion/mesh/mesh/rectangular_uniform.py index fad20eeda..a8a9c0b9d 100644 --- a/autoarray/inversion/mesh/mesh/rectangular_uniform.py +++ b/autoarray/inversion/mesh/mesh/rectangular_uniform.py @@ -15,9 +15,10 @@ class RectangularUniform(RectangularRTUAdaptDensity): - Indices increase left-to-right across rows and top-to-bottom between rows. - Each source-plane coordinate is associated with the rectangular pixel - in which it lies. No interpolation is performed — every coordinate - contributes fully to a single pixel. + Each source-plane coordinate is interpolated bilinearly between the four + surrounding pixel centres of the uniform overlay grid (extrapolating + linearly beyond the outermost centres), so it contributes to up to four + pixels. Uniform behaviour ----------------- diff --git a/autoarray/inversion/mesh/mesh_geometry/rectangular.py b/autoarray/inversion/mesh/mesh_geometry/rectangular.py index 0ca2227ab..d764fa879 100644 --- a/autoarray/inversion/mesh/mesh_geometry/rectangular.py +++ b/autoarray/inversion/mesh/mesh_geometry/rectangular.py @@ -434,11 +434,14 @@ def origin(self) -> Tuple[float, float]: @property def areas_transformed(self): """ - The source-plane area of every pixel in the rectangular pixelization, - transformed through the kernel-density CDF the mapper adapts with. - - Cell boundaries are the mapper's node-midpoint edges, including - clamped guard cells. Reusing ``edges_transformed`` also keeps the + The source-plane area of every pixel in the rectangular pixelization: + the outer product of the cell spacings of ``edges_transformed``. + + For the adaptive meshes the cells are transformed through the CDF the + mapper adapts with (rank or kernel) and their boundaries are the + mapper's node-midpoint edges, including clamped guard cells. For + ``RectangularUniform`` (``transform="uniform"``) they are the uniform + cells centred on the overlay nodes (see ``edges_transformed``). Reusing ``edges_transformed`` also keeps the slim data grid aligned with its adapt-image weights. """ edges = self.edges_transformed @@ -461,8 +464,9 @@ def areas_for_magnification(self): @property def edges_transformed(self): """ - The source-plane cell edges of every mesh pixel, transformed through - the adaptive CDF — one more edge than pixels per axis, shape + The source-plane cell edges of every mesh pixel — transformed through + the adaptive CDF for the adaptive meshes, uniform for + ``RectangularUniform`` — one more edge than pixels per axis, shape ``(n + 1, 2)`` packing (y_edges, x_edges). The interpolation mapper is node-based: reconstruction value @@ -479,7 +483,11 @@ def edges_transformed(self): mapper interpolates bilinearly between the overlay grid's pixel centres with no CDF and no guard ring, so its edges are the uniform cell boundaries half a node spacing either side of every node (the - same cells its ``imshow`` plot draws). Routing it through the kernel + same cells its ``imshow`` plot draws). These cells tile the overlay + box exactly, so the area sum is exact; per node they approximate the + quadrature of the clipped + extrapolated bilinear interpolant only + near the boundary ring, which is negligible for a source away from the + mesh edge. Routing it through the kernel CDF and guard-ring convention instead biased its ``areas_for_magnification`` by -77% to +259% across the audited configurations (PyAutoArray#659). 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 2d8ce3caa..00547be78 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 @@ -53,13 +53,17 @@ # with this exact configuration (PyAutoArray#659) and set ~1.5-2x above the # measurement, quoted in brackets. # -# Exact quadrature (tight): +# Tight: # - Delaunay: barycentric dual areas are the exact quadrature of the linear # interpolant (PyAutoArray#524) [-1.1e-4]. -# - RectangularUniform: uniform cells centred on the overlay nodes, the trapezoid -# weights of its bilinear interpolant [-1.2e-3]. Before PyAutoArray#659 its areas -# came from the adaptive meshes' kernel-CDF guard-ring geometry, which does not -# describe the uniform mapper at all, and this error was +55%. +# - RectangularUniform: uniform cells centred on the overlay nodes. They tile the +# overlay box exactly (the area SUM is exact), but per node they only +# approximate the quadrature of the clipped + extrapolated bilinear interpolant +# near the boundary ring (measured Q/A ~1.74 corner, ~1.32 edge, ~0.81 second +# ring, ~1.0 interior); negligible for a source away from the mesh edge, hence +# [-1.2e-3]. Before PyAutoArray#659 its areas came from the adaptive meshes' +# kernel-CDF guard-ring geometry, which does not describe the uniform mapper at +# all, and this error was +55%. # # Documented approximations (flagged in PyAutoArray#659 for a ruling): # - RectangularRTUAdaptDensity: node-midpoint cells, accurate here because a From 32a9073eb1cffe3346338678ea097cedb55a9805 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Thu, 29 Oct 2026 20:36:13 +0000 Subject: [PATCH 5/6] fix: exact hat-quadrature areas_for_magnification on adaptive rectangular meshes The adaptive rectangular meshes' areas_for_magnification pushed each node's node-midpoint cell through the inverse CDF. That is not the integral of the node's bilinear basis function once the CDF bends across a cell (the rank CDF's staircase on tied coordinates, a steep adapt-image CDF), and biased magnifications by up to +14% lensed / +11% unlensed. New adaptive_rectangular_hat_areas_from integrates each node's separable basis hat_r(g_y(y)) * hat_c(g_x(x)) exactly per axis, through the same forward transform the mapper uses, over the over-sampled data's bounding box. JAX-jit/grad safe; numpy and JAX agree to 1e-14. areas_transformed and edges_transformed (plotting) are unchanged. Lensed workspace check vs simulator truth: BilinearAdaptImage +14.43% -> -0.57%, RTUAdaptImage +8.56% -> +0.08%; identity BilinearAdaptDensity +11.23% -> -0.34%. Tolerances tightened to match. Refs #659 Co-Authored-By: Claude Opus 5.5 --- .../mesh/interpolator/rectangular.py | 78 +++++++++++++++++++ .../mesh/mesh_geometry/rectangular.py | 41 +++++++++- .../test_areas_for_magnification.py | 47 +++++------ .../mesh_geometry/test_rectangular.py | 14 +++- 4 files changed, 153 insertions(+), 27 deletions(-) diff --git a/autoarray/inversion/mesh/interpolator/rectangular.py b/autoarray/inversion/mesh/interpolator/rectangular.py index 20118f373..527f9eefa 100644 --- a/autoarray/inversion/mesh/interpolator/rectangular.py +++ b/autoarray/inversion/mesh/interpolator/rectangular.py @@ -478,6 +478,84 @@ def inv_full(U): return xp.abs(xp.outer(dy, dx).flatten()) +# Samples per axis of the 1D quadrature in `adaptive_rectangular_hat_areas_from`. +# Against 40001 samples on a 12x12 mesh (PyAutoArray#659) the per-node error is +# ~1e-6 on the smooth kernel CDF and ~6e-4 on the rank CDF's staircase over an +# unlensed (fully tied) grid, its worst case. +HAT_AREA_QUADRATURE_SAMPLES: int = 4097 + + +def adaptive_rectangular_hat_areas_from( + source_grid_size: int, + data_grid, + data_grid_over_sampled, + mesh_weight_map=None, + bandwidth: float = KERNEL_CDF_DEFAULT_BANDWIDTH, + n_knots: int = KERNEL_CDF_DEFAULT_KNOTS, + transform: str = "kernel", + n_samples: int = HAT_AREA_QUADRATURE_SAMPLES, + xp=np, +): + """ + The exact source-plane integral of every node's interpolation basis + function on the adaptive rectangular mesh — the quadrature weights for which + ``sum(reconstruction * areas)`` is the integral of the reconstructed source. + + The mapper of ``adaptive_rectangular_mappings_weights_via_interpolation_from`` + gives node ``(r, c)`` the separable basis function + + phi_rc(y, x) = hat_r(g_y(y)) * hat_c(g_x(x)), + g(.) = (source_grid_size - 3) * U(.) + 1, + + where ``U`` is the forward CDF transform and ``hat_k(g) = max(0, 1 - |g - k|)`` + in index space (rows flipped, ``row = n - i``). Its integral is therefore + the outer product of two 1D integrals, ``q_y(r) * q_x(c)``, each evaluated + here by trapezoid quadrature of ``hat_k(g(t))`` along the axis through the + SAME forward transform the mapper uses. Pushing the node-midpoint cell edges + through the inverse CDF instead (``areas_transformed``) is only correct where + the CDF is linear across a cell; on a steep adapt-image CDF or the rank CDF's + tied-coordinate staircase it was biased by up to ~20% per node. + + The domain is the bounding box of the over-sampled data grid — the region the + mapper actually samples — so the areas sum to that box's area. Guard nodes + (rows 0 and 1, columns 0 and n - 1) are never reached by ``g`` and get zero. + Differentiable under JAX (no sorts on the kernel path; the rank path's sort + is the same one the mapper already uses). + """ + mu = data_grid.mean(axis=0) + scale = data_grid.std(axis=0).min() + source_grid_scaled = (data_grid - mu) / scale + + transform_func, _ = _transforms_from( + source_grid_scaled, + mesh_pixels=source_grid_size, + mesh_weight_map=mesh_weight_map, + bandwidth=bandwidth, + n_knots=n_knots, + transform=transform, + xp=xp, + ) + + lo = xp.min(data_grid_over_sampled, axis=0) + hi = xp.max(data_grid_over_sampled, axis=0) + u = xp.linspace(0.0, 1.0, n_samples) + samples = lo[None, :] + u[:, None] * (hi - lo)[None, :] + + g = (source_grid_size - 3) * transform_func((samples - mu) / scale) + 1 + + # The mapper's `flatten` puts y index-space node k on row n - k and x node k + # on column k, so row r holds y node n - r (row 0 is node n, never reached). + n = source_grid_size + nodes = xp.stack([n - xp.arange(n), xp.arange(n)], axis=0) + hats = xp.clip(1.0 - xp.abs(g[:, :, None] - nodes[None, :, :]), 0.0, None) + + # trapezoid rule along each axis: (n_samples, 2, n) -> (2, n) + dt = ((hi - lo) / (n_samples - 1))[:, None] + q = dt * (xp.sum(hats, axis=0) - 0.5 * (hats[0] + hats[-1])) + + return xp.outer(q[0], q[1]).ravel() + + def adaptive_rectangular_mappings_weights_via_interpolation_from( source_grid_size: int, data_grid, diff --git a/autoarray/inversion/mesh/mesh_geometry/rectangular.py b/autoarray/inversion/mesh/mesh_geometry/rectangular.py index d764fa879..f323c0a33 100644 --- a/autoarray/inversion/mesh/mesh_geometry/rectangular.py +++ b/autoarray/inversion/mesh/mesh_geometry/rectangular.py @@ -452,14 +452,49 @@ def areas_transformed(self): @property def areas_for_magnification(self): """ - The area of every pixel in the rectangular pixelization. + The quadrature weight of every pixel in the rectangular pixelization: + ``sum(reconstruction * areas_for_magnification)`` is the integral of + the reconstructed source. + + For the adaptive meshes this is the exact integral of each node's + bilinear basis function through the mapper's forward CDF (see + ``adaptive_rectangular_hat_areas_from``), not the node-midpoint cell + area of ``areas_transformed``, which was biased by up to ~20% per node + wherever the CDF bends across a cell (PyAutoArray#659). For + ``RectangularUniform`` it is the uniform cell area. Returns ------- ndarray - The area of every pixel in the rectangular pixelization. + The quadrature weight of every pixel in the rectangular pixelization. """ - return self.areas_transformed + if self.transform == "uniform": + return self.areas_transformed + + from autoarray.inversion.mesh.interpolator.rectangular import ( + KERNEL_CDF_DEFAULT_BANDWIDTH, + KERNEL_CDF_DEFAULT_KNOTS, + adaptive_rectangular_hat_areas_from, + ) + + return adaptive_rectangular_hat_areas_from( + source_grid_size=self.shape_native[0], + data_grid=self.data_grid.array, + data_grid_over_sampled=self.data_grid.over_sampled.array, + mesh_weight_map=self.mesh_weight_map, + bandwidth=( + self.kernel_bandwidth + if self.kernel_bandwidth is not None + else KERNEL_CDF_DEFAULT_BANDWIDTH + ), + n_knots=( + self.kernel_knots + if self.kernel_knots is not None + else KERNEL_CDF_DEFAULT_KNOTS + ), + transform=self.transform, + xp=self._xp, + ) @property def edges_transformed(self): 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 00547be78..37c8620ca 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 @@ -16,10 +16,11 @@ and `sum_k s_k * A_k / sum_k s_k * Q_k` is the identity-lens magnification, which must be 1. The checks are: -(a) the areas sum to the analytic area of the domain the mesh covers — the data - grid's bounding box for the adaptive rectangular meshes (the CDF maps that box - onto the unit square), the overlay box for `RectangularUniform`, and the convex - hull of the vertices for the Delaunay family (barycentric dual areas tile it); +(a) the areas sum to the analytic area of the domain the mesh covers — the + over-sampled data grid's bounding box for the rectangular meshes (the adaptive + meshes integrate their basis functions over it; `RectangularUniform`'s overlay + grid spans it), and the convex hull of the vertices for the Delaunay family + (barycentric dual areas tile it); (b) for a smooth Gaussian source reconstructed on the mesh, the identity-lens magnification is 1 and the integrated reconstruction matches the analytic integral, within a per-mesh tolerance; @@ -65,28 +66,27 @@ # kernel-CDF guard-ring geometry, which does not describe the uniform mapper at # all, and this error was +55%. # -# Documented approximations (flagged in PyAutoArray#659 for a ruling): -# - RectangularRTUAdaptDensity: node-midpoint cells, accurate here because a -# density-adapted kernel CDF of an unlensed grid is near-linear [+2.5e-3]; the -# lensed workspace check (autolens_workspace_test -# scripts/imaging/mesh_magnification.py) measures -1.8%. -# - RectangularBilinearAdaptDensity / RectangularBilinearAdaptImage / -# RectangularRTUAdaptImage: `areas_for_magnification` is each node's -# node-midpoint cell pushed through the inverse CDF. That is not the integral of -# the node's bilinear hat once the CDF is non-linear across a cell (the rank -# CDF's staircase on tied coordinates; a steep adapt-image CDF). The exact -# per-axis hat quadrature brings these to ~1e-3 in a prototype -# [+3.1e-2, +2.4e-1, +8.0e-2; ~1-3% for all three at 20x20 and finer]. +# - The adaptive rectangular meshes: the exact per-axis integral of each node's +# bilinear hat through the mapper's forward CDF +# (`adaptive_rectangular_hat_areas_from`). The residual is this coarse grid's +# sub-pixel quadrature `Q`, not the areas [BilinearAdaptDensity -3.1e-4, +# BilinearAdaptImage -4.2e-3, RTUAdaptDensity +2.3e-3, RTUAdaptImage -5.3e-3]. +# The node-midpoint cells they replaced (PyAutoArray#659) were off by +# [+3.1e-2, +2.4e-1, +2.5e-3, +8.0e-2]: the inverse CDF of a cell's edges is +# 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]. MAGNIFICATION_TOLERANCE = { "RectangularUniform": 2.5e-3, - "RectangularBilinearAdaptDensity": 5.0e-2, - "RectangularBilinearAdaptImage": 3.5e-1, + "RectangularBilinearAdaptDensity": 1.0e-3, + "RectangularBilinearAdaptImage": 8.0e-3, "RectangularRTUAdaptDensity": 5.0e-3, - "RectangularRTUAdaptImage": 1.2e-1, + "RectangularRTUAdaptImage": 1.0e-2, "Delaunay": 3.0e-4, "DelaunayNN": 1.0e-2, "KNearestNeighbor": 3.0e-2, @@ -190,8 +190,9 @@ def test__areas_sum_to_the_covered_domain(name): domain = np.prod(np.ptp(np.asarray(grid.over_sampled.array), axis=0)) rel = 1.0e-6 else: - # the adaptive CDF maps the data grid's bounding box onto the unit square - domain = np.prod(np.ptp(np.asarray(grid.array), axis=0)) + # the hats are integrated over the over-sampled data's bounding box and + # partition unity on it + domain = np.prod(np.ptp(np.asarray(grid.over_sampled.array), axis=0)) rel = 1.0e-10 assert areas.sum() == pytest.approx(domain, rel=rel) @@ -219,8 +220,8 @@ def test__rectangular_adaptive__guard_cells_have_zero_area_and_no_flux(name): """ The adaptive rectangular mapper is node-based with a guard ring in unit CDF space: rows 0 and 1 (above U_y = 1) and columns 0 and n - 1 (outside - [0, 1] in U_x) never receive mapping weight, and their node-midpoint edges clamp - to the data span so their area is exactly zero. Every other cell is live and + [0, 1] in U_x) never receive mapping weight, and their hats are never reached + by the transformed coordinate so their area is exactly zero. Every other cell is live and carries a positive area — including the bottom row (U_y = 0), which is live although it sits in the `zeroed_pixels` edge ring. """ diff --git a/test_autoarray/inversion/pixelization/mesh_geometry/test_rectangular.py b/test_autoarray/inversion/pixelization/mesh_geometry/test_rectangular.py index 075f6167a..345c51973 100644 --- a/test_autoarray/inversion/pixelization/mesh_geometry/test_rectangular.py +++ b/test_autoarray/inversion/pixelization/mesh_geometry/test_rectangular.py @@ -294,4 +294,16 @@ def test__areas_transformed__adapt_image_matches_cell_edges(over_sample_size): assert np.all(np.isfinite(areas)) assert np.all(areas >= 0.0) np.testing.assert_allclose(areas, expected, rtol=1e-12, atol=1e-12) - np.testing.assert_array_equal(geometry.areas_for_magnification, areas) + + # `areas_for_magnification` is the exact hat quadrature, not these cells + # (PyAutoArray#659): positive exactly where the mapper puts flux (the cells + # above give this coarse mesh's live bottom row zero area), summing to the + # over-sampled data box. + hat_areas = np.asarray(geometry.areas_for_magnification) + mapped = np.asarray(mapper.mapping_matrix_over_sampled).sum(axis=0) > 0.0 + assert hat_areas.shape == (36,) + assert np.all(hat_areas >= 0.0) + np.testing.assert_array_equal(hat_areas > 0.0, mapped) + assert hat_areas.sum() == pytest.approx( + np.prod(np.ptp(np.asarray(grid.over_sampled.array), axis=0)), rel=1e-10 + ) From ac3b54b36d70c0f34546e056fc9291e5d9ce076e Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Thu, 29 Oct 2026 21:03:06 +0000 Subject: [PATCH 6/6] test: skip KNN mesh areas tests without jax (unittest-nojax leg) Co-Authored-By: Claude Opus 5.5 --- .../mesh_geometry/test_areas_for_magnification.py | 4 ++++ 1 file changed, 4 insertions(+) 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 37c8620ca..23db88d96 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 @@ -139,6 +139,10 @@ def _setup(name): The mapper, its areas, its exact quadrature `Q` and a Gaussian reconstruction, cached per mesh so the parametrized tests share one build. """ + if name in ("KNearestNeighbor", "KNNBarycentric"): + # the KNN interpolators import jax (unittest-nojax leg) + pytest.importorskip("jax") + if name in _CACHE: return _CACHE[name]