feat: Mapper.mesh_magnification — source-plane magnification per mesh pixel

autoarray/inversion/mappers/abstract.py

@@ -581,7 +581,44 @@ def data_weight_total_for_pix_from(self) -> np.ndarray:
         )
 
     def data_pixel_area_for_pix_from(self) -> np.ndarray:
-        pass
+        """
+        Returns the image-plane area, in arc-seconds squared, that maps to every pixelization pixel.
+
+        Each column of the `mapping_matrix` gives the fraction of every image pixel that maps to that pixelization
+        pixel (over-sampled sub-pixels contribute their sub-fraction, spread over the pixels they interpolate to), so
+        its column sum is the number of image pixels mapping to it. Multiplying by the area of one image pixel gives
+        the image-plane area.
+
+        Summed over every pixelization pixel this is the total area of the mask.
+        """
+        pixel_area = self.mask.pixel_scales[0] * self.mask.pixel_scales[1]
+
+        return np.asarray(self.mapping_matrix).sum(axis=0) * pixel_area
+
+    @property
+    def mesh_magnification(self) -> np.ndarray:
+        """
+        The magnification of every pixelization pixel: the image-plane area that maps to it divided by its
+        source-plane area.
+
+        For a lensed source this is how much the lens magnifies the region of the source plane each mesh pixel
+        covers, summed over all the images of that region. The source-plane areas are the mesh's
+        `areas_for_magnification` (exact cell areas for a rectangular mesh, barycentric dual areas for a Delaunay mesh,
+        integrated basis weights for a k-nearest-neighbour mesh), so the image area attributed to each pixel follows
+        the same interpolation as the reconstruction and `sum(mesh_magnification * areas_for_magnification)` is the
+        total area of the mask.
+
+        Pixels with zero source-plane area, or with no image area mapping to them, have a magnification of zero.
+        """
+        image_areas = self.data_pixel_area_for_pix_from()
+        source_areas = np.asarray(self.mesh_geometry.areas_for_magnification)
+
+        return np.divide(
+            image_areas,
+            source_areas,
+            out=np.zeros_like(image_areas),
+            where=source_areas > 0.0,
+        )
 
     def mapped_to_source_from(self, array: Array2D) -> np.ndarray:
         """

autoarray/plot/inversion.py

@@ -34,6 +34,7 @@ def plot_inversion_reconstruction(
     vmin: Optional[float] = None,
     vmax: Optional[float] = None,
     use_log10: bool = False,
+    cb_unit: Optional[str] = None,
     zoom_to_brightest: bool = True,
     zoom_extent_scale: float = 1.0,
     arcsec_after_decimal: Optional[bool] = None,
@@ -74,6 +75,8 @@ def plot_inversion_reconstruction(
         Explicit colour scale limits.
     use_log10
         Apply ``LogNorm``.
+    cb_unit
+        Unit appended to the middle colour-bar label. ``None`` reads it from config; ``""`` shows no unit.
     zoom_to_brightest
         Pass through to ``mapper.extent_from``.
     arcsec_after_decimal
@@ -141,13 +144,27 @@ def plot_inversion_reconstruction(
         mapper.interpolator, (InterpolatorRectangular, InterpolatorRectangularUniform)
     ):
         raster_extent = _plot_rectangular(
-            ax, pixel_values, mapper, norm, colormap, extent, is_subplot=is_subplot
+            ax,
+            pixel_values,
+            mapper,
+            norm,
+            colormap,
+            extent,
+            is_subplot=is_subplot,
+            cb_unit=cb_unit,
         )
     elif isinstance(
         mapper.interpolator, (InterpolatorDelaunay, InterpolatorKNearestNeighbor)
     ):
         _plot_delaunay(
-            ax, pixel_values, mapper, norm, colormap, extent, is_subplot=is_subplot
+            ax,
+            pixel_values,
+            mapper,
+            norm,
+            colormap,
+            extent,
+            is_subplot=is_subplot,
+            cb_unit=cb_unit,
         )
 
     # --- overlays --------------------------------------------------------------
@@ -229,7 +246,7 @@ def _for_origin(yx):
 
 
 def _plot_rectangular(
-    ax, pixel_values, mapper, norm, colormap, extent, is_subplot=False
+    ax, pixel_values, mapper, norm, colormap, extent, is_subplot=False, cb_unit=None
 ):
     """Render a rectangular pixelization reconstruction onto *ax*.
 
@@ -303,7 +320,7 @@ def _plot_rectangular(
             aspect="auto",
             origin=_conf_imshow_origin(),
         )
-        _apply_colorbar(im, ax, is_subplot=is_subplot)
+        _apply_colorbar(im, ax, cb_unit=cb_unit, is_subplot=is_subplot)
         return pix_array.geometry.extent
     else:
         y_edges, x_edges = mapper.mesh_geometry.edges_transformed.T
@@ -316,11 +333,13 @@ def _plot_rectangular(
             norm=norm,
             cmap=colormap,
         )
-        _apply_colorbar(im, ax, is_subplot=is_subplot)
+        _apply_colorbar(im, ax, cb_unit=cb_unit, is_subplot=is_subplot)
         return None
 
 
-def _plot_delaunay(ax, pixel_values, mapper, norm, colormap, extent, is_subplot=False):
+def _plot_delaunay(
+    ax, pixel_values, mapper, norm, colormap, extent, is_subplot=False, cb_unit=None
+):
     """Render a Delaunay or KNN pixelization reconstruction onto *ax*.
 
     Uses ``ax.tripcolor`` with its default *flat* shading: each triangle is
@@ -400,4 +419,4 @@ def _plot_delaunay(ax, pixel_values, mapper, norm, colormap, extent, is_subplot=
         tc = ax.tripcolor(x, y, vals, cmap=colormap, norm=norm)
     from autoarray.plot.utils import _apply_colorbar
 
-    _apply_colorbar(tc, ax, is_subplot=is_subplot)
+    _apply_colorbar(tc, ax, cb_unit=cb_unit, is_subplot=is_subplot)

test_autoarray/inversion/pixelization/mappers/test_mesh_magnification.py

@@ -0,0 +1,84 @@
+"""
+`Mapper.data_pixel_area_for_pix_from` and `Mapper.mesh_magnification` (PyAutoArray#673).
+
+The tests reuse the identity lens of `mesh_geometry/test_areas_for_magnification.py`: the source-plane data grid is
+the image grid itself, so the true magnification of every source-plane region is exactly 1.
+
+The checks are:
+
+(a) the image-plane area mapped to the pixels sums to the mask area, and `sum(mesh_magnification * areas)` gives it
+    back, for every mesh;
+(b) away from the mesh edge each pixel's magnification is close to 1, within a per-mesh tolerance;
+(c) a pixel with zero source-plane area has a magnification of zero.
+
+Per-pixel tolerances (median |mu - 1| over pixels whose centre is inside |y|, |x| < 1.5 on this module's coarse 24x24
+grid, 2x2 sub-grid and 12x12-equivalent mesh) are set ~2x the measured value, quoted in brackets. The image area a
+pixel receives follows the mapper's interpolation weights, while its source-plane area follows
+`areas_for_magnification`, so the per-pixel error is set by how well those match per node:
+
+- Tight: RectangularUniform [9.1e-4], RectangularBilinearAdaptDensity [5.5e-4], RectangularRTUAdaptDensity
+  [1.3e-3], Delaunay [2.9e-3], KNearestNeighbor [1.3e-3].
+- Looser: RectangularRTUAdaptImage [7.9e-3], RectangularBilinearAdaptImage [1.6e-2] (a steep adapt-image CDF on this
+  coarse grid); KNNBarycentric [2.5e-2] (discontinuous clipped weights).
+- Documented approximation: DelaunayNN [9.2e-2] uses the Delaunay dual areas, exact only for the linear interpolant,
+  not for its Sibson weights. The error averages out over a source (phase-4 total magnification 4.3e-3) but not
+  per pixel.
+"""
+
+import numpy as np
+import pytest
+
+from test_autoarray.inversion.pixelization.mesh_geometry.test_areas_for_magnification import (
+    PIXEL_SCALE,
+    SHAPE,
+    _setup,
+)
+
+MEDIAN_TOLERANCE = {
+    "RectangularUniform": 2.0e-3,
+    "RectangularBilinearAdaptDensity": 1.5e-3,
+    "RectangularBilinearAdaptImage": 3.5e-2,
+    "RectangularRTUAdaptDensity": 3.0e-3,
+    "RectangularRTUAdaptImage": 1.6e-2,
+    "Delaunay": 6.0e-3,
+    "DelaunayNN": 0.2,
+    "KNearestNeighbor": 3.0e-3,
+    "KNNBarycentric": 5.0e-2,
+}
+
+MASK_AREA = (SHAPE * PIXEL_SCALE) ** 2
+
+
+@pytest.mark.parametrize("name", list(MEDIAN_TOLERANCE))
+def test__image_areas_and_magnification_integrate_to_the_mask_area(name):
+    _, mapper, areas, _, _ = _setup(name)
+
+    image_areas = mapper.data_pixel_area_for_pix_from()
+
+    assert image_areas.shape == (mapper.params,)
+    assert image_areas.sum() == pytest.approx(MASK_AREA, rel=1.0e-10)
+    assert np.sum(mapper.mesh_magnification * areas) == pytest.approx(
+        MASK_AREA, rel=1.0e-6
+    )
+
+
+@pytest.mark.parametrize("name", list(MEDIAN_TOLERANCE))
+def test__identity_lens__interior_pixels_have_unit_magnification(name):
+    _, mapper, areas, _, _ = _setup(name)
+
+    centres = np.asarray(mapper.source_plane_mesh_grid)
+    interior = (np.abs(centres) < 1.5).all(axis=1) & (areas > 0.0)
+
+    mu = mapper.mesh_magnification[interior]
+
+    assert np.all(np.isfinite(mapper.mesh_magnification))
+    assert np.median(np.abs(mu - 1.0)) < MEDIAN_TOLERANCE[name]
+
+
+def test__zero_area_pixels_have_zero_magnification():
+    _, mapper, _, _, _ = _setup("RectangularBilinearAdaptDensity")
+
+    areas = np.asarray(mapper.mesh_geometry.areas_for_magnification)
+
+    assert np.any(areas == 0.0)
+    assert np.all(mapper.mesh_magnification[areas == 0.0] == 0.0)