Merge pull request #326 from Jammy2211/feature/positions_lh_mass_centre

feature/positions_lh_mass_centre

autolens/analysis/result.py

@@ -165,7 +165,7 @@ def image_plane_multiple_image_positions_for_single_image_from(
                 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.
 
@@ -239,6 +239,7 @@ 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
@@ -249,6 +250,11 @@ def positions_likelihood_from(
         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
@@ -264,6 +270,10 @@ def positions_likelihood_from(
         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
         -------
@@ -278,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
         )

test_autolens/analysis/test_result.py

@@ -291,6 +291,39 @@ def test__positions_likelihood_from(analysis_imaging_7x7):
     assert positions_likelihood.threshold == pytest.approx(0.2, 1.0e-4)
 
 
+def test__positions_likelihood_from__mass_centre_radial_distance_min(
+    analysis_imaging_7x7,
+):
+    tracer = al.Tracer(
+        galaxies=[
+            al.Galaxy(
+                redshift=0.5,
+                mass=al.mp.Isothermal(
+                    centre=(0.1, 0.0), einstein_radius=1.0, ell_comps=(0.0, 0.0)
+                ),
+            ),
+            al.Galaxy(redshift=1.0, bulge=al.lp.SersicSph(centre=(0.0, 0.0))),
+        ]
+    )
+
+    samples_summary = al.m.MockSamplesSummary(max_log_likelihood_instance=tracer)
+
+    result = res.Result(samples_summary=samples_summary, analysis=analysis_imaging_7x7)
+
+    positions_likelihood = result.positions_likelihood_from(
+        factor=0.1, minimum_threshold=0.2, mass_centre_radial_distance_min=0.1
+    )
+
+    assert isinstance(positions_likelihood, al.PositionsLHPenalty)
+    assert len(positions_likelihood.positions) == 2
+    assert positions_likelihood.positions[0] == pytest.approx(
+        (-1.00097656e00, 5.63818622e-04), 1.0e-4
+    )
+    assert positions_likelihood.positions[1] == pytest.approx(
+        (1.00097656e00, -5.63818622e-04), 1.0e-4
+    )
+
+
 def test__results_include_mask__available_as_property(
     analysis_imaging_7x7, masked_imaging_7x7, samples_summary_with_result
 ):