feat: image↔source mapping regions, clumps, subplot_mappings (P1)

autoarray/__init__.py

@@ -37,6 +37,8 @@
 from .inversion.mesh.border_relocator import BorderRelocator
 from .inversion.pixelization import Pixelization
 from .inversion.mappers.abstract import Mapper
+from .inversion.mappings.mapping import ImageRegion
+from .inversion.mappings.mapping import Mapping
 from .inversion.mesh.image_mesh.abstract import AbstractImageMesh
 from .inversion.mesh.mesh.abstract import AbstractMesh
 from .inversion.mesh.interpolator.rectangular import InterpolatorRectangular

autoarray/config/visualize/general.yaml

@@ -8,7 +8,9 @@ general:
   output_format: show                   # Default output format: "show" displays the figure interactively via plt.show(), "png"/"pdf"/etc. saves to file.
 inversion:
   reconstruction_vmax_factor: 0.5
-  total_mappings_pixels: 8              # The number of source pixels used when plotting the subplot_mappings of a pixelization.
+  total_mappings: 5                     # The maximum number of source clumps drawn by subplot_mappings.
+  mappings_threshold: 0.5               # A source pixel joins a clump if its reconstructed value exceeds this fraction of the reconstruction's maximum.
+  mappings_min_pixels: 3                # Connected groups of source pixels smaller than this are not drawn as a clump.
 zoom:
   plane_percent: 0.01
   inversion_percent: 0.01               # Plots of an Inversion's reconstruction use the reconstructed data's bright value multiplied by this factor.

autoarray/inversion/inversion/abstract.py

@@ -1311,3 +1311,121 @@ def max_pixel_centre(self, mapper_index: int = 0) -> Grid2DIrregular:
         )
 
         return max_pixel_centre
+
+    def source_clumps_from(
+        self,
+        mapper_index: int = 0,
+        threshold: float = 0.5,
+        min_pixels: int = 3,
+        total_clumps: Optional[int] = None,
+        pix_indexes: Optional[List] = None,
+    ) -> List[np.ndarray]:
+        """
+        Returns the clumps of the reconstruction, each a group of connected mesh pixels which are
+        bright relative to the reconstruction's maximum value.
+
+        Mesh pixels whose reconstructed value exceeds `threshold * max(reconstruction)` are kept,
+        split into connected components over the mesh neighbour graph, filtered by `min_pixels` and
+        ordered brightest first.
+
+        In gravitational lensing a clump is a distinct bright structure of the source galaxy. The
+        `threshold` is therefore the knob which decides what counts as a distinct structure: ~0.5
+        isolates one smooth source, ~0.2 merges two nearby galaxies into one clump and ~0.8 splits
+        a source into its individual star-forming knots.
+
+        Parameters
+        ----------
+        mapper_index
+            The index of the mapper in the inversion whose reconstruction is clumped, where there
+            may be multiple mappers in the inversion.
+        threshold
+            Mesh pixels are in a clump if their reconstructed value exceeds this fraction of the
+            reconstruction's maximum value.
+        min_pixels
+            Connected groups with fewer than this many mesh pixels are discarded.
+        total_clumps
+            The maximum number of clumps returned, keeping the brightest. `None` returns all of them.
+        pix_indexes
+            An explicit list of mesh pixel index groups which bypasses the clump finding entirely,
+            returning them unchanged.
+
+        Returns
+        -------
+        A list of 1D integer arrays, one per clump, ordered by decreasing peak reconstructed value.
+        """
+        from autoarray.inversion.mappings.mapping import connected_components_from
+        from autoarray.inversion.mappings.mapping import pix_index_groups_from
+
+        if pix_indexes is not None:
+            return pix_index_groups_from(pix_indexes=pix_indexes)
+
+        mapper = self.cls_list_from(cls=Mapper)[mapper_index]
+
+        reconstruction = np.asarray(self.reconstruction_dict[mapper])
+
+        if reconstruction.size == 0:
+            return []
+
+        max_value = float(np.max(reconstruction))
+
+        if max_value <= 0.0:
+            return []
+
+        indexes = np.where(reconstruction > threshold * max_value)[0]
+
+        if indexes.size == 0:
+            return []
+
+        clumps = connected_components_from(indexes=indexes, neighbors=mapper.neighbors)
+
+        clumps = [clump for clump in clumps if clump.shape[0] >= min_pixels]
+
+        clumps.sort(key=lambda clump: -float(np.max(reconstruction[clump])))
+
+        if total_clumps is not None:
+            clumps = clumps[:total_clumps]
+
+        return clumps
+
+    def mappings_from(
+        self,
+        mapper_index: int = 0,
+        weight_threshold: float = 0.0,
+        **clump_kwargs,
+    ) -> List["Mapping"]:
+        """
+        Returns the `Mapping` of every clump of the reconstruction, pairing each source-plane clump
+        with the image-plane regions (the multiple images) it maps to.
+
+        Parameters
+        ----------
+        mapper_index
+            The index of the mapper in the inversion whose reconstruction is mapped, where there may
+            be multiple mappers in the inversion.
+        weight_threshold
+            A data pixel is in an image-plane region if its summed mapping weight to the clump
+            exceeds this value.
+        clump_kwargs
+            The `threshold`, `min_pixels`, `total_clumps` and `pix_indexes` inputs
+            of `source_clumps_from`.
+
+        Returns
+        -------
+        One `Mapping` per source-plane clump, ordered by decreasing peak reconstructed value.
+        """
+        from autoarray.inversion.mappings import mapping as mapping_module
+
+        mapper = self.cls_list_from(cls=Mapper)[mapper_index]
+
+        reconstruction = np.asarray(self.reconstruction_dict[mapper])
+
+        clumps = self.source_clumps_from(mapper_index=mapper_index, **clump_kwargs)
+
+        peak_values = [float(np.max(reconstruction[clump])) for clump in clumps]
+
+        return mapping_module.mappings_from(
+            mapper=mapper,
+            pix_indexes=clumps,
+            weight_threshold=weight_threshold,
+            peak_values=peak_values,
+        )

autoarray/inversion/mappers/abstract.py

@@ -619,7 +619,9 @@ def extent_from(
                 x_min, x_max, y_min, y_max = extent
                 x_centre = 0.5 * (x_min + x_max)
                 y_centre = 0.5 * (y_min + y_max)
-                target_half = 0.5 * max(x_max - x_min, y_max - y_min) * zoom_extent_scale
+                target_half = (
+                    0.5 * max(x_max - x_min, y_max - y_min) * zoom_extent_scale
+                )
 
                 bound = geometry_util.extent_symmetric_from(
                     extent=self.source_plane_mesh_grid.geometry.extent
@@ -630,8 +632,8 @@ def extent_from(
                     y_centre - bound[2],
                     bound[3] - y_centre,
                 )
-                bound_cap_half = 0.7 * 0.5 * min(
-                    bound[1] - bound[0], bound[3] - bound[2]
+                bound_cap_half = (
+                    0.7 * 0.5 * min(bound[1] - bound[0], bound[3] - bound[2])
                 )
                 final_half = min(target_half, max_allowable_half, bound_cap_half)
 
@@ -668,3 +670,41 @@ def mesh_pixels_per_image_pixels(self):
             values=mesh_pixels_per_image_pixels,
             mask=self.mask,
         )
+
+    def mappings_from(
+        self,
+        pix_indexes,
+        weight_threshold: float = 0.0,
+        min_pixels: int = 1,
+    ) -> List["Mapping"]:
+        """
+        Returns the `Mapping` of one or more groups of mesh pixels, pairing each source-plane group
+        with the image-plane regions (the multiple images) it maps to.
+
+        This is the bare-`Mapper` route into the mapping objects, used when there is no inversion to
+        find clumps in and the mesh pixels of interest are chosen by hand (e.g. in a tutorial). Each
+        returned `Mapping` therefore has a `peak_value` of `None`.
+
+        Parameters
+        ----------
+        pix_indexes
+            The mesh pixel indexes, either a flat sequence (one group) or a nested sequence (one
+            group per entry).
+        weight_threshold
+            A data pixel is in an image-plane region if its summed mapping weight to the group
+            exceeds this value.
+        min_pixels
+            Connected image-plane regions with fewer than this many pixels are discarded.
+
+        Returns
+        -------
+        One `Mapping` per group of mesh pixels.
+        """
+        from autoarray.inversion.mappings import mapping as mapping_module
+
+        return mapping_module.mappings_from(
+            mapper=self,
+            pix_indexes=pix_indexes,
+            weight_threshold=weight_threshold,
+            min_pixels=min_pixels,
+        )

autoarray/inversion/mappings/__init__.py

@@ -0,0 +1,7 @@
+from autoarray.inversion.mappings.mapping import ImageRegion
+from autoarray.inversion.mappings.mapping import Mapping
+from autoarray.inversion.mappings.mapping import connected_components_from
+from autoarray.inversion.mappings.mapping import contours_from_bool_native
+from autoarray.inversion.mappings.mapping import image_regions_from
+from autoarray.inversion.mappings.mapping import image_regions_from_slim_mask
+from autoarray.inversion.mappings.mapping import source_contours_from