fix: adaptive rectangular meshes report their true pixel centres

autoarray/inversion/mappers/abstract.py

@@ -140,7 +140,11 @@ def source_plane_data_grid(self):
 
     @property
     def source_plane_mesh_grid(self):
-        return self.interpolator.mesh_grid
+        return self.interpolator.source_plane_mesh_grid
+
+    @property
+    def regularization_mesh_grid(self):
+        return self.interpolator.regularization_mesh_grid
 
     @property
     def adapt_data(self):

autoarray/inversion/mesh/interpolator/abstract.py

@@ -30,6 +30,24 @@ def _xp(self):
             return jnp
         return np
 
+    @property
+    def source_plane_mesh_grid(self):
+        """
+        The source-plane (y,x) position of every mesh pixel, in the order of the reconstruction. For most meshes
+        this is the `mesh_grid` the interpolator was built from; the adaptive rectangular interpolator overrides it.
+        """
+        return self.mesh_grid
+
+    @property
+    def regularization_mesh_grid(self):
+        """
+        The (y,x) positions kernel regularizations measure pixel-to-pixel distances on. For most meshes these are
+        the source-plane pixel positions. For the adaptive rectangular meshes they are the regular lattice
+        `mesh_grid`, so distances are measured in the mesh's stretched (CDF) coordinates, see
+        `InterpolatorRectangular.source_plane_mesh_grid`.
+        """
+        return self.mesh_grid
+
     @property
     def _mappings_sizes_weights(self):
         raise NotImplementedError(

autoarray/inversion/mesh/interpolator/rectangular.py

@@ -478,6 +478,73 @@ def inv_full(U):
     return xp.abs(xp.outer(dy, dx).flatten())
 
 
+def adaptive_rectangular_node_centres_from(
+    source_grid_shape,
+    data_grid,
+    mesh_weight_map=None,
+    bandwidth: float = KERNEL_CDF_DEFAULT_BANDWIDTH,
+    n_knots: int = KERNEL_CDF_DEFAULT_KNOTS,
+    transform: str = "kernel",
+    xp=np,
+):
+    """
+    The source-plane (y,x) centre of every node of the adaptive rectangular mesh,
+    in the mapper's pixel order — the position where node ``(r, c)``'s
+    interpolation basis function peaks, so that reconstruction value ``k`` lives
+    at centre ``k``.
+
+    The mapper of ``adaptive_rectangular_mappings_weights_via_interpolation_from``
+    places the source plane in index space via ``g = (source_grid_size - 3) * U + 1``
+    with ``U`` the forward CDF transform, and its ``flatten`` puts y index-space
+    node ``k`` on row ``n - k`` and x node ``k`` on column ``k``. Node ``(r, c)``
+    is therefore at unit-square coordinate
+
+        U_y = (n - r - 1) / (n - 3),    U_x = (c - 1) / (n - 3),
+
+    mapped back through the SAME inverse CDF. The transform is separable, so each
+    axis is inverted once and the centres are their outer grid (row-major from
+    the top-left, like the uniform overlay this replaces).
+
+    Guard nodes (rows 0 and 1, columns 0 and ``n - 1``) have ``U`` outside
+    ``[0, 1]``; the inverse CDF clamps them to the edge of the data span, so they
+    sit on the border of the region the mapper samples (they carry zero
+    interpolation weight and zero area).
+
+    The mapper indexes a square mesh; for a non-square ``source_grid_shape`` each
+    axis uses its own pixel count, so there is still one centre per pixel.
+    """
+    n_y, n_x = source_grid_shape
+
+    mu = data_grid.mean(axis=0)
+    scale = data_grid.std(axis=0).min()
+    source_grid_scaled = (data_grid - mu) / scale
+
+    _, inv_transform = _transforms_from(
+        source_grid_scaled,
+        mesh_pixels=n_y,
+        mesh_weight_map=mesh_weight_map,
+        bandwidth=bandwidth,
+        n_knots=n_knots,
+        transform=transform,
+        xp=xp,
+    )
+
+    # the transform is per axis, so invert each axis's node coordinates in its own
+    # column (the other column is unused)
+    rows = xp.arange(n_y)
+    cols = xp.arange(n_x)
+    unit_y = (n_y - rows - 1) / (n_y - 3)
+    unit_x = (cols - 1) / (n_x - 3)
+
+    y_nodes = inv_transform(xp.stack([unit_y, xp.zeros(n_y)], axis=1))[:, 0]
+    x_nodes = inv_transform(xp.stack([xp.zeros(n_x), unit_x], axis=1))[:, 1]
+
+    y = xp.repeat(y_nodes * scale + mu[0], n_x)
+    x = xp.tile(x_nodes * scale + mu[1], n_y)
+
+    return xp.stack([y, x], axis=1)
+
+
 # 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
@@ -765,6 +832,23 @@ def mesh_geometry(self):
             xp=self._xp,
         )
 
+    @cached_property
+    def source_plane_mesh_grid(self):
+        """
+        The source-plane (y,x) centre of every mesh pixel, where the interpolant takes each reconstruction value
+        (``mesh_geometry.node_centres``). Use these to place, plot, export or interpolate the reconstruction.
+
+        ``mesh_grid`` is the uniform lattice before the adaptive CDF is applied. It is kept for the mesh
+        geometry's extent and is the ``regularization_mesh_grid``: kernel regularizations (``MaternKernel``,
+        ``GaussianKernel``, ...) and ``MaternSparse`` measure distances on it, i.e. in the mesh's stretched
+        coordinates, so a kernel's ``scale`` shrinks in the source plane where the mesh packs its pixels.
+        Distances on the true pixel centres gave a lower Bayesian evidence on a test lens, so the lattice is the
+        deliberate choice (PyAutoArray#667).
+        """
+        from autoarray.structures.grids.irregular_2d import Grid2DIrregular
+
+        return Grid2DIrregular(values=self.mesh_geometry.node_centres)
+
     @cached_property
     def _mappings_sizes_weights(self):
         mappings, weights = (

autoarray/inversion/mesh/mesh_geometry/rectangular.py

@@ -496,6 +496,47 @@ def areas_for_magnification(self):
             xp=self._xp,
         )
 
+    @property
+    def node_centres(self):
+        """
+        The source-plane (y,x) centre of every pixel, in the order of the
+        reconstruction: the position at which the mapper's interpolant takes
+        each pixel's value.
+
+        For the adaptive meshes these are the nodes of the regular lattice in
+        CDF space mapped back through the mapper's inverse CDF (see
+        ``adaptive_rectangular_node_centres_from``), so they bunch up where the
+        mesh adapts. The uniform overlay ``mesh_grid`` is the lattice before
+        that mapping and is not where the values live. For
+        ``RectangularUniform`` the two are the same.
+        """
+        if self.transform == "uniform":
+            return self.mesh_grid.array
+
+        from autoarray.inversion.mesh.interpolator.rectangular import (
+            KERNEL_CDF_DEFAULT_BANDWIDTH,
+            KERNEL_CDF_DEFAULT_KNOTS,
+            adaptive_rectangular_node_centres_from,
+        )
+
+        return adaptive_rectangular_node_centres_from(
+            source_grid_shape=self.shape_native,
+            data_grid=self.data_grid.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):
         """

autoarray/inversion/mock/mock_mapper.py

@@ -51,6 +51,12 @@ def source_plane_mesh_grid(self):
             return super().source_plane_mesh_grid
         return self._source_plane_mesh_grid
 
+    @property
+    def regularization_mesh_grid(self):
+        if self._source_plane_mesh_grid is None:
+            return super().regularization_mesh_grid
+        return self._source_plane_mesh_grid
+
     @property
     def adapt_data(self):
         if self._adapt_data is None:

autoarray/inversion/regularization/abstract.py

@@ -44,6 +44,22 @@ def validate_coefficient(coefficient, name: str = "coefficient"):
     )
 
 
+def regularization_mesh_grid_from(linear_obj):
+    """
+    The (y,x) positions a kernel regularization measures pixel-to-pixel distances on.
+
+    Mappers expose `regularization_mesh_grid` (the regular lattice for the adaptive rectangular meshes, see
+    `InterpolatorRectangular.source_plane_mesh_grid`). Other linear objects, including duck-typed ones in
+    downstream libraries, only define `source_plane_mesh_grid`, which is used for them.
+    """
+    grid = getattr(linear_obj, "regularization_mesh_grid", None)
+
+    if grid is None:
+        return linear_obj.source_plane_mesh_grid
+
+    return grid
+
+
 class AbstractRegularization:
     is_split_regularization = False
     """

autoarray/inversion/regularization/exponential_kernel.py

@@ -7,6 +7,7 @@
 
 from autoarray.inversion.regularization.abstract import AbstractRegularization
 from autoarray.inversion.regularization.abstract import validate_coefficient
+from autoarray.inversion.regularization.abstract import regularization_mesh_grid_from
 
 
 def exp_cov_matrix_from(
@@ -161,7 +162,7 @@ def regularization_matrix_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray
 
         covariance_matrix = exp_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             jitter=self.jitter_value,
             jitter_relative=self.jitter_relative,
             xp=xp,
@@ -185,7 +186,7 @@ def log_det_regularization_matrix_term_from(
         """
         covariance_matrix = exp_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             jitter=self.jitter_value,
             jitter_relative=self.jitter_relative,
             xp=xp,
@@ -217,7 +218,7 @@ def regularization_term_from(
 
         covariance_matrix = exp_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             jitter=self.jitter_value,
             jitter_relative=self.jitter_relative,
             xp=xp,

autoarray/inversion/regularization/gaussian_kernel.py

@@ -7,6 +7,7 @@
 
 from autoarray.inversion.regularization.abstract import AbstractRegularization
 from autoarray.inversion.regularization.abstract import validate_coefficient
+from autoarray.inversion.regularization.abstract import regularization_mesh_grid_from
 
 
 def gauss_cov_matrix_from(
@@ -150,7 +151,7 @@ def regularization_matrix_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray
 
         covariance_matrix = gauss_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             jitter=self.jitter_value,
             jitter_relative=self.jitter_relative,
             xp=xp,
@@ -198,7 +199,7 @@ def log_det_regularization_matrix_term_from(
         """
         covariance_matrix = gauss_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             jitter=self.jitter_value,
             jitter_relative=self.jitter_relative,
             xp=xp,
@@ -236,7 +237,7 @@ def regularization_term_from(
 
         covariance_matrix = gauss_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             jitter=self.jitter_value,
             jitter_relative=self.jitter_relative,
             xp=xp,

autoarray/inversion/regularization/matern_adapt_kernel.py

@@ -15,6 +15,7 @@
 )
 from autoarray.inversion.regularization.adapt import adapt_regularization_weights_from
 from autoarray.inversion.regularization.abstract import validate_coefficient
+from autoarray.inversion.regularization.abstract import regularization_mesh_grid_from
 
 
 class MaternAdaptKernel(MaternKernel):
@@ -146,7 +147,7 @@ def regularization_matrix_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray
             linear_obj=linear_obj, xp=xp
         )
 
-        pixel_points = linear_obj.source_plane_mesh_grid.array
+        pixel_points = regularization_mesh_grid_from(linear_obj).array
 
         covariance_matrix = matern_cov_matrix_from(
             scale=self.scale,
@@ -178,7 +179,7 @@ def log_det_regularization_matrix_term_from(
 
         covariance_matrix = matern_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             nu=self.nu,
             weights=kernel_weights,
             jitter=self.jitter_value,
@@ -210,7 +211,7 @@ def regularization_term_from(
 
         covariance_matrix = matern_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             nu=self.nu,
             weights=kernel_weights,
             jitter=self.jitter_value,

autoarray/inversion/regularization/matern_kernel.py

@@ -8,6 +8,7 @@
 
 from autoarray.inversion.regularization.abstract import AbstractRegularization
 from autoarray.inversion.regularization.abstract import validate_coefficient
+from autoarray.inversion.regularization.abstract import regularization_mesh_grid_from
 
 
 def kv_xp(v, z, xp=np):
@@ -367,7 +368,7 @@ def regularization_matrix_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray
         """
         covariance_matrix = matern_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             nu=self.nu,
             jitter=self.jitter_value,
             jitter_relative=self.jitter_relative,
@@ -390,7 +391,7 @@ def log_det_regularization_matrix_term_from(
         """
         covariance_matrix = matern_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             nu=self.nu,
             jitter=self.jitter_value,
             jitter_relative=self.jitter_relative,
@@ -419,7 +420,7 @@ def regularization_term_from(
         """
         covariance_matrix = matern_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             nu=self.nu,
             jitter=self.jitter_value,
             jitter_relative=self.jitter_relative,