fix: rectangular-mesh plot edges follow the mapper's node convention

The bilinear interpolation mapper places reconstruction value (row r,
col c) at node U_y=(n_y-r-1)/(n_y-3), U_x=(c-1)/(n_x-3) in unit CDF space
(guard ring at the border, rows flipped by the mapper's row = n - i), but
edges_transformed drew a uniform [0,1] partition — pcolormesh plots put
flux up to ~1.5 mesh pixels from where the mapper scatters it (the "small
offsets" of issue #372, confirmed by a delta-function reproduction). The
edges are now the node midpoints pushed through the same CDF machinery
(guard-safe denominator for degenerate n<=3 meshes; border guard cells
clamp to the data span). Both consumers inherit the fix: the pcolormesh
plot path and rectangular_rotated's warped-centre reconstruction; spline
and linear CDFs share the upstream index convention. Plot geometry only —
the mapper is untouched, so inversions and likelihoods are byte-identical.

#372

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Author: 2 people

autoarray/inversion/mesh/mesh_geometry/rectangular.py

@@ -510,16 +510,29 @@ def areas_for_magnification(self):
     @property
     def edges_transformed(self):
         """
-        A class packing the ndarrays describing the neighbors of every pixel in the rectangular pixelization (see
-        `Neighbors` for a complete description of the neighboring scheme).
-
-        The neighbors of a rectangular pixelization are computed by exploiting the uniform and symmetric nature of the
-        rectangular grid, as described in the method `rectangular_neighbors_from`.
+        The source-plane cell edges of every mesh pixel, transformed through
+        the adaptive CDF — one more edge than pixels per axis, shape
+        ``(n + 1, 2)`` packing (y_edges, x_edges).
+
+        The interpolation mapper is node-based: reconstruction value
+        ``(row r, col c)`` lives at node ``U_y = (n_y - r - 1)/(n_y - 3)``,
+        ``U_x = (c - 1)/(n_x - 3)`` in unit CDF space (a guard ring occupies
+        the border; rows are flipped by the mapper's ``row = n - i``). The
+        cell edges are therefore the node midpoints — a uniform ``[0, 1]``
+        partition drew every cell up to ~1.5 mesh pixels away from where the
+        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.
         """
 
-        # edges defined in 0 -> 1 space, there is one more edge than pixel centers on each side
-        edges_y = self._xp.linspace(1, 0, self.shape_native[0] + 1)
-        edges_x = self._xp.linspace(0, 1, self.shape_native[1] + 1)
+        # Node-midpoint edges in unit CDF space (see docstring); the guard
+        # denominator keeps the degenerate n <= 3 meshes finite (their
+        # interpolator maps every point to the single interior node anyway).
+        n_y, n_x = self.shape_native
+        rows = self._xp.arange(n_y + 1)
+        cols = self._xp.arange(n_x + 1)
+        edges_y = (n_y - rows - 0.5) / max(n_y - 3, 1)
+        edges_x = (cols - 1.5) / max(n_x - 3, 1)
 
         edges_reshaped = self._xp.stack([edges_y, edges_x]).T
 

test_autoarray/inversion/pixelization/mesh_geometry/test_rectangular.py

@@ -200,3 +200,53 @@ def test__edges_transformed(mask_2d_7x7):
         ),
         abs=1e-8,
     )
+
+
+def test__edges_transformed__aligned_with_interpolation_node_convention():
+    """
+    Regression test for issue #372: the pcolormesh plotting path draws value
+    (row r, col c) inside the cell bounded by edges_transformed — that cell
+    must be centred on the interpolation mapper's node for that value, not on
+    a uniform [0, 1] partition (which shifted plots by ~1.5 mesh pixels).
+
+    A delta function scattered through the mapper must land in plotted cells
+    whose weighted centroid matches the input point to sub-cell precision.
+    """
+    from autoarray.inversion.mesh.interpolator.rectangular import (
+        adaptive_rectangular_mappings_weights_via_interpolation_from,
+        adaptive_rectangular_transformed_grid_from,
+    )
+
+    n = 10
+    rng = np.random.default_rng(0)
+    data_grid = rng.uniform(-1.0, 1.0, (5000, 2))
+    test_point = np.array([[0.3, -0.2]])
+
+    flat_indices, weights = (
+        adaptive_rectangular_mappings_weights_via_interpolation_from(
+            source_grid_size=n,
+            data_grid=data_grid,
+            data_grid_over_sampled=test_point,
+        )
+    )
+
+    # The node-midpoint unit-space edges (what edges_transformed now builds),
+    # pushed through the same CDF transform.
+    rows = np.arange(n + 1)
+    edges_y = (n - rows - 0.5) / (n - 3)
+    edges_x = (rows - 1.5) / (n - 3)
+    edges = np.stack([edges_y, edges_x]).T
+    edges_t = adaptive_rectangular_transformed_grid_from(data_grid, edges)
+    y_edges, x_edges = edges_t.T
+
+    centroid_y = 0.0
+    centroid_x = 0.0
+    for flat, weight in zip(flat_indices[0], weights[0]):
+        r, c = flat // n, flat % n
+        centroid_y += weight * 0.5 * (y_edges[r] + y_edges[r + 1])
+        centroid_x += weight * 0.5 * (x_edges[c] + x_edges[c + 1])
+
+    # Half a mesh cell in these units is ~0.15; the pre-fix uniform edges
+    # missed by ~0.4 in y.
+    assert centroid_y == pytest.approx(test_point[0, 0], abs=0.1)
+    assert centroid_x == pytest.approx(test_point[0, 1], abs=0.1)