fix: exact mesh magnification areas (RectangularUniform + adaptive rectangular) with per-mesh tests (arc phase 4a)

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,

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,
         )
 

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

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
 

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
@@ -449,20 +452,56 @@ 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):
         """
-        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
@@ -474,7 +513,22 @@ 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). 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).
         """
+        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 +569,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:
         """