diff --git a/autoarray/inversion/mappers/abstract.py b/autoarray/inversion/mappers/abstract.py index 83a37d883..154b2cfa4 100644 --- a/autoarray/inversion/mappers/abstract.py +++ b/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: """ diff --git a/autoarray/plot/inversion.py b/autoarray/plot/inversion.py index 7b8da7e2c..a109736ea 100644 --- a/autoarray/plot/inversion.py +++ b/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) diff --git a/test_autoarray/inversion/pixelization/mappers/test_mesh_magnification.py b/test_autoarray/inversion/pixelization/mappers/test_mesh_magnification.py new file mode 100644 index 000000000..dd1307c74 --- /dev/null +++ b/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)