merge complete

Author: Jammy2211

autolens/__init__.py

@@ -123,4 +123,4 @@
 
 conf.instance.register(__file__)
 
-__version__ = "2024.9.21.2"
+__version__ = "2024.11.13.2"

autolens/aggregator/tracer.py

@@ -41,13 +41,15 @@ def _tracer_from(
         galaxies = instance.galaxies
 
         if hasattr(instance, "extra_galaxies"):
-            galaxies = galaxies + fit.instance.extra_galaxies
+            if fit.instance.extra_galaxies is not None:
+                galaxies = galaxies + fit.instance.extra_galaxies
 
     else:
         galaxies = fit.instance.galaxies
 
         if hasattr(fit.instance, "extra_galaxies"):
-            galaxies = galaxies + fit.instance.extra_galaxies
+            if fit.instance.extra_galaxies is not None:
+                galaxies = galaxies + fit.instance.extra_galaxies
 
     try:
         cosmology = instance.cosmology
@@ -56,17 +58,18 @@ def _tracer_from(
 
     tracer = Tracer(galaxies=galaxies, cosmology=cosmology)
 
-    if len(fit.children) > 0:
-        logger.info(
-            """
-            Using database for a fit with multiple summed Analysis objects.
-
-            Tracer objects do not fully support this yet (e.g. model parameters which vary over analyses may be incorrect)
-            so proceed with caution!
-            """
-        )
-
-        return [tracer] * len(fit.children)
+    if fit.children is not None:
+        if len(fit.children) > 0:
+            logger.info(
+                """
+                Using database for a fit with multiple summed Analysis objects.
+    
+                Tracer objects do not fully support this yet (e.g. model parameters which vary over analyses may be incorrect)
+                so proceed with caution!
+                """
+            )
+
+            return [tracer] * len(fit.children)
 
     return [tracer]
 

autolens/analysis/result.py

@@ -1,3 +1,4 @@
+import logging
 import os
 import numpy as np
 from typing import Optional, Union
@@ -15,6 +16,8 @@
 from autolens.lens.tracer import Tracer
 from autolens.point.solver import PointSolver
 
+logger = logging.getLogger(__name__)
+
 
 class Result(AgResultDataset):
     @property
@@ -101,8 +104,76 @@ def image_plane_multiple_image_positions(self) -> aa.Grid2DIrregular:
             source_plane_coordinate=self.source_plane_centre.in_list[0],
         )
 
+        if multiple_images.shape[0] == 1:
+            return self.image_plane_multiple_image_positions_for_single_image_from()
+
         return aa.Grid2DIrregular(values=multiple_images)
 
+    def image_plane_multiple_image_positions_for_single_image_from(
+        self, increments: int = 20
+    ) -> aa.Grid2DIrregular:
+        """
+        If the standard point solver only locates one multiple image, finds one or more additional images, which are
+        not technically multiple image in the point source regime, but are close enough to it they can be used
+        in a position threshold likelihood.
+
+        This is performed by incrementally moving the source-plane centre's coordinates towards the centre of the
+        source-plane at (0.0", 0.0"). This ensures that the centre will eventually go inside the caustic, where
+        multiple images are formed.
+
+        To move the source-plane centre, the original source-plane centre is multiplied by a factor that decreases
+        from 1.0 to 0.0 in increments of 1/increments. For example, if the source-plane centre is (1.0", -0.5") and
+        the `factor` is 0.5, the input source-plane centre is (0.5", -0.25").
+
+        The multiple images are always computed for the same mass model, thus they will always be valid multiple images
+        for the model being fitted, but as the factor decrease the multiple images may move furhter from their observed
+        positions.
+
+        Parameters
+        ----------
+        increments
+            The number of increments the source-plane centre is moved to compute multiple images.
+        """
+
+        logger.info(
+            """
+        Could not find multiple images for maximum likelihood lens model.
+        
+        Incrementally moving source centre inwards towards centre of source-plane until caustic crossing occurs
+        and multiple images are formed.        
+        """
+        )
+
+        grid = self.analysis.dataset.mask.derive_grid.all_false
+
+        centre = self.source_plane_centre.in_list[0]
+
+        solver = PointSolver.for_grid(
+            grid=grid,
+            pixel_scale_precision=0.001,
+        )
+
+        for i in range(1, increments):
+            factor = 1.0 - (1.0 * (i / increments))
+
+            multiple_images = solver.solve(
+                tracer=self.max_log_likelihood_tracer,
+                source_plane_coordinate=(centre[0] * factor, centre[1] * factor),
+            )
+
+            if multiple_images.shape[0] > 1:
+                return aa.Grid2DIrregular(values=multiple_images)
+
+        logger.info(
+            """
+        Could not find multiple images for maximum likelihood lens model, even after incrementally moving the source
+        centre inwards to the centre of the source-plane.
+
+        Set the multiple image postiions to two images at (1.0", 1.0") so code continues to run.
+        """
+        )
+        return aa.Grid2DIrregular(values=[(1.0, 1.0), (1.0, 1.0)])
+
     def positions_threshold_from(
         self,
         factor=1.0,
@@ -168,7 +239,47 @@ def positions_likelihood_from(
         minimum_threshold=None,
         use_resample=False,
         positions: Optional[aa.Grid2DIrregular] = None,
+        mass_centre_radial_distance_min: float = None,
     ) -> Union[PositionsLHPenalty, PositionsLHResample]:
+        """
+        Returns a `PositionsLH` object from the result of a lens model-fit, where the maximum log likelihood mass
+        and source models are used to determine the multiple image positions in the image-plane and source-plane
+        and ray-trace them to the source-plane to determine the threshold.
+
+        In chained fits, for example the SLaM pipelines, this means that a simple initial fit (e.g. SIE mass model,
+        parametric source) can be used to determine the multiple image positions and threshold for a more complex
+        subsequent fit (e.g. power-law mass model, pixelized source).
+
+        The mass model central image is removed from the solution, as this is rarely physically observed and therefore
+        should not be included in the likelihood penalty or resampling. It is removed by setting a positive
+        magnification threshold in the `PointSolver`. For strange lens models the central image may still be
+        solved for, in which case the `mass_centre_radial_distance_min` parameter can be used to remove it.
+
+        Parameters
+        ----------
+        factor
+            The value the computed threshold is multiplied by to make the position threshold larger or smaller than the
+            maximum log likelihood model's threshold.
+        minimum_threshold
+            The output threshold is rounded up to this value if it is below it, to avoid extremely small threshold
+            values.
+        use_resample
+            If `False` the `PositionsLH` object is created using the `PositionsLHPenalty` class, which uses the
+            threshold to apply a penalty term to the likelihood. If `True` the `PositionsLH` object is created using
+            the `PositionsLHResample` class, which resamples the positions to the threshold.
+        positions
+            If input, these positions are used instead of the computed multiple image positions from the lens mass
+            model.
+        mass_centre_radial_distance_min
+            The minimum radial distance from the mass model centre that a multiple image position must be to be
+            included in the likelihood penalty or resampling. If `None` all positions are used. This is an additional
+            method to remove central images that may make it through the point solver's magnification threshold.
+
+        Returns
+        -------
+        The `PositionsLH` object used to apply a likelihood penalty or resample the positions.
+        """
+
         if os.environ.get("PYAUTOFIT_TEST_MODE") == "1":
             return None
 
@@ -177,6 +288,18 @@ def positions_likelihood_from(
             if positions is None
             else positions
         )
+
+        if mass_centre_radial_distance_min is not None:
+            mass_centre = self.max_log_likelihood_tracer.extract_attribute(
+                cls=ag.mp.MassProfile, attr_name="centre"
+            )
+
+            distances = positions.distances_to_coordinate_from(
+                coordinate=mass_centre[0]
+            )
+
+            positions = positions[distances > mass_centre_radial_distance_min]
+
         threshold = self.positions_threshold_from(
             factor=factor, minimum_threshold=minimum_threshold, positions=positions
         )

autolens/config/visualize/plots.yaml

@@ -7,6 +7,10 @@
     psf: false
   positions:                                 # Settings for plots with resampling image-positions on (e.g. the image).
     image_with_positions: true
+  point_dataset:                             # Settings for plots of point source datasets (e.g. PointDatasetPlotter).
+    subplot_dataset: true                    # Plot subplot containing all dataset quantities (e.g. the data, noise-map, etc.)?
+    positions: false                         # Plot the positions of th multiple images of the point source in the image plane?
+    fluxes: false                            # Plot the fluxes of the multiple images of the point source in the image plane?
   fit:                                       # Settings for plots of all fits (e.g. FitImagingPlotter, FitInterferometerPlotter).
     subplot_fit: true                        # Plot subplot of all fit quantities for any dataset (e.g. the model data, residual-map, etc.)?
     all_at_end_png: true                     # Plot all individual plots listed below as .png (even if False)?
@@ -25,6 +29,9 @@
     subtracted_images_of_planes: false       # Plot individual plots of each plane's subtracted image?
     plane_images_of_planes: false            # Plot individual plots of each plane's image (e.g. in the source plane)?
   fit_imaging: {}                            # Settings for plots of fits to imaging datasets (e.g. FitImagingPlotter).
+  fit_point_dataset:                         # Settings for plots of fits to point source datasets (e.g. FitPointDatasetPlotter).
+    positions: false                         # Plot the positions of th multiple images of the point source in the image plane?
+    fluxes: false                            # Plot the fluxes of the multiple images of the point source in the image plane?
   tracer:                               # Settings for plots of tracers (e.g. TracerPlotter).
     subplot_tracer: true                # Plot subplot of all quantities in each tracer (e.g. images, convergence)?
     all_at_end_png: true                     # Plot all individual plots listed below as .png (even if False)?
@@ -48,7 +55,7 @@
     all_at_end_png: true                     # Plot all individual plots listed below as .png (even if False)?
     all_at_end_fits: true                    # Plot all individual plots listed below as .fits (even if False)?
     all_at_end_pdf: false                    # Plot all individual plots listed below as publication-quality .pdf (even if False)?
-    errors: false
+    reconstruction_noise_map: false
     reconstructed_image: false
     reconstruction: false
     regularization_weights: false

autolens/imaging/plot/fit_imaging_plotters.py

@@ -148,6 +148,8 @@ def figures_2d_of_planes(
         subtracted_image: bool = False,
         model_image: bool = False,
         plane_image: bool = False,
+        plane_noise_map: bool = False,
+        plane_signal_to_noise_map: bool = False,
         use_source_vmax: bool = False,
         zoom_to_brightest: bool = True,
         interpolate_to_uniform: bool = False,
@@ -178,6 +180,14 @@ def figures_2d_of_planes(
             Whether to make a 2D plot (via `imshow`) of the image of a plane in its source-plane (e.g. unlensed).
             Depending on how the fit is performed, this could either be an image of light profiles of the reconstruction
             of an `Inversion`.
+        plane_noise_map
+            Whether to make a 2D plot of the noise-map of a plane in its source-plane, where the
+            noise map can only be computed when a pixelized source reconstruction is performed and they correspond to
+            the noise map in each reconstructed pixel as given by the inverse curvature matrix.
+        plane_signal_to_noise_map
+            Whether to make a 2D plot of the signal-to-noise map of a plane in its source-plane,
+            where the signal-to-noise map values can only be computed when a pixelized source reconstruction and they
+            are the ratio of reconstructed flux to error in each pixel.
         use_source_vmax
             If `True`, the maximum value of the lensed source (e.g. in the image-plane) is used to set the `vmax` of
             certain plots (e.g. the `data`) in order to ensure the lensed source is visible compared to the lens.
@@ -283,6 +293,8 @@ def figures_2d_of_planes(
 
                 elif self.tracer.planes[plane_index].has(cls=aa.Pixelization):
 
+                    pix = self.tracer.planes[plane_index].cls_list_from(cls=aa.Pixelization)[0]
+
                     inversion_plotter = self.inversion_plotter_of_plane(
                         plane_index=plane_index
                     )
@@ -294,8 +306,41 @@ def figures_2d_of_planes(
                         interpolate_to_uniform=interpolate_to_uniform
                     )
 
-        if use_source_vmax:
-            self.mat_plot_2d.cmap.kwargs.pop("vmax")
+            if use_source_vmax:
+                try:
+                    self.mat_plot_2d.cmap.kwargs.pop("vmax")
+                except KeyError:
+                    pass
+
+            if plane_noise_map:
+
+                if self.tracer.planes[plane_index].has(cls=aa.Pixelization):
+
+                    inversion_plotter = self.inversion_plotter_of_plane(
+                        plane_index=plane_index
+                    )
+
+                    inversion_plotter.figures_2d_of_pixelization(
+                        pixelization_index=0,
+                        reconstruction_noise_map=True,
+                        zoom_to_brightest=zoom_to_brightest,
+                        interpolate_to_uniform=interpolate_to_uniform
+                    )
+
+            if plane_signal_to_noise_map:
+
+                if self.tracer.planes[plane_index].has(cls=aa.Pixelization):
+
+                    inversion_plotter = self.inversion_plotter_of_plane(
+                        plane_index=plane_index
+                    )
+
+                    inversion_plotter.figures_2d_of_pixelization(
+                        pixelization_index=0,
+                        signal_to_noise_map=True,
+                        zoom_to_brightest=zoom_to_brightest,
+                        interpolate_to_uniform=interpolate_to_uniform
+                    )
 
     def subplot_of_planes(self, plane_index: Optional[int] = None):
         """
@@ -614,7 +659,12 @@ def subplot_mappings_of_plane(self, plane_index: Optional[int] = None, auto_file
                     "total_mappings_pixels"
                 ]
 
-                pix_indexes = inversion_plotter.inversion.brightest_pixel_list_from(
+                mapper = inversion_plotter.inversion.cls_list_from(cls=aa.AbstractMapper)[0]
+                mapper_valued = aa.MapperValued(
+                    values=inversion_plotter.inversion.reconstruction_dict[mapper],
+                    mapper=mapper,
+                )
+                pix_indexes = mapper_valued.max_pixel_list_from(
                     total_pixels=total_pixels, filter_neighbors=True
                 )