From 07a956589e4296de9dffc44b5c435dc8b4cb0296 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 14 May 2025 20:48:55 +0100 Subject: [PATCH 01/11] extended PositionsLHPenalty to use plane index dict --- autolens/analysis/positions.py | 84 ++++++++++++++++-------- autolens/point/max_separation.py | 8 ++- setup.py | 2 +- test_autolens/analysis/test_positions.py | 3 +- 4 files changed, 67 insertions(+), 30 deletions(-) diff --git a/autolens/analysis/positions.py b/autolens/analysis/positions.py index 90a7cfe84..5c78ff16f 100644 --- a/autolens/analysis/positions.py +++ b/autolens/analysis/positions.py @@ -2,6 +2,7 @@ from typing import Optional, Union from os import path import os +from typing import Dict import autoarray as aa import autofit as af @@ -21,7 +22,12 @@ class AbstractPositionsLH: - def __init__(self, positions: aa.Grid2DIrregular, threshold: float): + def __init__( + self, + positions: aa.Grid2DIrregular, + threshold: float, + plane_index_positions_dict: Optional[Dict[int, aa.Grid2DIrregular]] = None, + ): """ The `PositionsLH` objects add a penalty term to the likelihood of the **PyAutoLens** `log_likelihood_function` defined in the `Analysis` classes. @@ -54,6 +60,9 @@ def __init__(self, positions: aa.Grid2DIrregular, threshold: float): self.positions = positions self.threshold = threshold + if plane_index_positions_dict is None: + self.plane_index_positions_dict = {-1: positions} + def log_likelihood_function_positions_overwrite( self, instance: af.ModelInstance, analysis: AnalysisDataset ) -> Optional[float]: @@ -77,21 +86,26 @@ def output_positions_info(self, output_path: str, tracer: Tracer): ------- """ - positions_fit = SourceMaxSeparation( - data=self.positions, noise_map=None, tracer=tracer - ) + with open_(path.join(output_path, "positions.info"), "w+") as f: - distances = positions_fit.data.distances_to_coordinate_from( - coordinate=(0.0, 0.0) - ) + for plane_index, positions in self.plane_index_positions_dict.items(): - with open_(path.join(output_path, "positions.info"), "w+") as f: - f.write(f"Positions: \n {self.positions} \n\n") - f.write(f"Radial Distance from (0.0, 0.0): \n {distances} \n\n") - f.write(f"Threshold = {self.threshold} \n") - f.write( - f"Max Source Plane Separation of Maximum Likelihood Model = {positions_fit.max_separation_of_source_plane_positions}" - ) + positions_fit = SourceMaxSeparation( + data=self.positions, noise_map=None, tracer=tracer + ) + + distances = positions_fit.data.distances_to_coordinate_from( + coordinate=(0.0, 0.0) + ) + + f.write(f"Plane Index: {plane_index} \n") + f.write(f"Positions: \n {self.positions} \n\n") + f.write(f"Radial Distance from (0.0, 0.0): \n {distances} \n\n") + f.write(f"Threshold = {self.threshold} \n") + f.write( + f"Max Source Plane Separation of Maximum Likelihood Model = {positions_fit.max_separation_of_source_plane_positions}" + ) + f.write("") class PositionsLHResample(AbstractPositionsLH): @@ -141,15 +155,20 @@ def log_likelihood_function_positions_overwrite( if not tracer.has(cls=ag.mp.MassProfile) or len(tracer.planes) == 1: return - positions_fit = SourceMaxSeparation( - data=self.positions, noise_map=None, tracer=tracer - ) + for plane_index, positions in self.plane_index_positions_dict.items(): + + positions_fit = SourceMaxSeparation( + data=self.positions, + noise_map=None, + tracer=tracer, + plane_index=plane_index, + ) - if not positions_fit.max_separation_within_threshold(self.threshold): - if os.environ.get("PYAUTOFIT_TEST_MODE") == "1": - return + if not positions_fit.max_separation_within_threshold(self.threshold): + if os.environ.get("PYAUTOFIT_TEST_MODE") == "1": + return - raise exc.RayTracingException + raise exc.RayTracingException class PositionsLHPenalty(AbstractPositionsLH): @@ -264,15 +283,26 @@ def log_likelihood_penalty_from(self, tracer: Tracer) -> Optional[float]: if not tracer.has(cls=ag.mp.MassProfile) or len(tracer.planes) == 1: return - positions_fit = SourceMaxSeparation( - data=self.positions, noise_map=None, tracer=tracer - ) + log_likelihood_penalty = 0.0 + + for plane_index, positions in self.plane_index_positions_dict.items(): - if not positions_fit.max_separation_within_threshold(self.threshold): - return self.log_likelihood_penalty_factor * ( - positions_fit.max_separation_of_source_plane_positions - self.threshold + positions_fit = SourceMaxSeparation( + data=self.positions, + noise_map=None, + tracer=tracer, + plane_index=plane_index, ) + if not positions_fit.max_separation_within_threshold(self.threshold): + + log_likelihood_penalty += self.log_likelihood_penalty_factor * ( + positions_fit.max_separation_of_source_plane_positions + - self.threshold + ) + + return log_likelihood_penalty + def log_likelihood_function_positions_overwrite( self, instance: af.ModelInstance, analysis: AnalysisDataset ) -> Optional[float]: diff --git a/autolens/point/max_separation.py b/autolens/point/max_separation.py index 2b60f0968..f6fcfa591 100644 --- a/autolens/point/max_separation.py +++ b/autolens/point/max_separation.py @@ -11,6 +11,7 @@ def __init__( data: aa.Grid2DIrregular, noise_map: Optional[aa.ArrayIrregular], tracer: Tracer, + plane_index: int = -1, ): """ Given a positions dataset, which is a list of positions with names that associated them to model source @@ -28,11 +29,16 @@ def __init__( The object that defines the ray-tracing of the strong lens system of galaxies. noise_value The noise-value assumed when computing the log likelihood. + plane_index + The index of the plane in the `Tracer` that the source-plane positions are computed from. This is typically + the last plane in the `Tracer`, which is the source-plane. """ self.data = data self.noise_map = noise_map - self.source_plane_positions = tracer.traced_grid_2d_list_from(grid=data)[-1] + self.source_plane_positions = tracer.traced_grid_2d_list_from(grid=data)[ + plane_index + ] @property def furthest_separations_of_source_plane_positions(self) -> aa.ArrayIrregular: diff --git a/setup.py b/setup.py index 4da4962ed..4c77d449f 100644 --- a/setup.py +++ b/setup.py @@ -5,4 +5,4 @@ setup( version=version, -) \ No newline at end of file +) diff --git a/test_autolens/analysis/test_positions.py b/test_autolens/analysis/test_positions.py index bf579252b..fa393e4d5 100644 --- a/test_autolens/analysis/test_positions.py +++ b/test_autolens/analysis/test_positions.py @@ -42,7 +42,8 @@ def test__output_positions_info(): with open_(positions_file, "r") as f: output = f.readlines() - assert "Positions" in output[0] + assert "Plane Index: -1" in output[0] + assert "Positions" in output[1] os.remove(positions_file) From f8109168b4dcf88e0a8b95efdf0bd82979cd7151 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 14 May 2025 20:58:26 +0100 Subject: [PATCH 02/11] double position input works --- autolens/analysis/positions.py | 36 ++++++++++++------- .../imaging/model/test_analysis_imaging.py | 27 ++++++++++++++ 2 files changed, 50 insertions(+), 13 deletions(-) diff --git a/autolens/analysis/positions.py b/autolens/analysis/positions.py index 5c78ff16f..ce69e97cc 100644 --- a/autolens/analysis/positions.py +++ b/autolens/analysis/positions.py @@ -24,8 +24,8 @@ class AbstractPositionsLH: def __init__( self, - positions: aa.Grid2DIrregular, threshold: float, + positions: Optional[aa.Grid2DIrregular] = None, plane_index_positions_dict: Optional[Dict[int, aa.Grid2DIrregular]] = None, ): """ @@ -50,19 +50,24 @@ def __init__( If the maximum separation of any two source plane coordinates is above the threshold the penalty term is applied. """ - if len(positions) == 1: - raise exc.PositionsException( - f"The positions input into the Positions have length one " - f"(e.g. it is only one (y,x) coordinate and therefore cannot be compared with other images).\n\n" - "Please input more positions into the Positions." - ) self.positions = positions self.threshold = threshold + self.plane_index_positions_dict = plane_index_positions_dict + if plane_index_positions_dict is None: self.plane_index_positions_dict = {-1: positions} + for positions in self.plane_index_positions_dict.values(): + + if len(positions) == 1: + raise exc.PositionsException( + f"The positions input into the PositionsLikelihood object have length one " + f"(e.g. it is only one (y,x) coordinate and therefore cannot be compared with other images).\n\n" + "Please input more positions into the Positions." + ) + def log_likelihood_function_positions_overwrite( self, instance: af.ModelInstance, analysis: AnalysisDataset ) -> Optional[float]: @@ -91,7 +96,7 @@ def output_positions_info(self, output_path: str, tracer: Tracer): for plane_index, positions in self.plane_index_positions_dict.items(): positions_fit = SourceMaxSeparation( - data=self.positions, noise_map=None, tracer=tracer + data=positions, noise_map=None, tracer=tracer ) distances = positions_fit.data.distances_to_coordinate_from( @@ -99,7 +104,7 @@ def output_positions_info(self, output_path: str, tracer: Tracer): ) f.write(f"Plane Index: {plane_index} \n") - f.write(f"Positions: \n {self.positions} \n\n") + f.write(f"Positions: \n {positions} \n\n") f.write(f"Radial Distance from (0.0, 0.0): \n {distances} \n\n") f.write(f"Threshold = {self.threshold} \n") f.write( @@ -158,7 +163,7 @@ def log_likelihood_function_positions_overwrite( for plane_index, positions in self.plane_index_positions_dict.items(): positions_fit = SourceMaxSeparation( - data=self.positions, + data=positions, noise_map=None, tracer=tracer, plane_index=plane_index, @@ -174,9 +179,10 @@ def log_likelihood_function_positions_overwrite( class PositionsLHPenalty(AbstractPositionsLH): def __init__( self, - positions: aa.Grid2DIrregular, threshold: float, log_likelihood_penalty_factor: float = 1e8, + positions: Optional[aa.Grid2DIrregular] = None, + plane_index_positions_dict: Optional[Dict[int, aa.Grid2DIrregular]] = None, ): """ The `PositionsLH` objects add a penalty term to the likelihood of the **PyAutoLens** `log_likelihood_function` @@ -207,7 +213,11 @@ def __init__( A factor which multiplies how far source pixels do not trace within the threshold of one another, with a larger factor producing a larger penalty making the non-linear parameter space gradient steeper. """ - super().__init__(positions=positions, threshold=threshold) + super().__init__( + positions=positions, + threshold=threshold, + plane_index_positions_dict=plane_index_positions_dict, + ) self.log_likelihood_penalty_factor = log_likelihood_penalty_factor @@ -288,7 +298,7 @@ def log_likelihood_penalty_from(self, tracer: Tracer) -> Optional[float]: for plane_index, positions in self.plane_index_positions_dict.items(): positions_fit = SourceMaxSeparation( - data=self.positions, + data=positions, noise_map=None, tracer=tracer, plane_index=plane_index, diff --git a/test_autolens/imaging/model/test_analysis_imaging.py b/test_autolens/imaging/model/test_analysis_imaging.py index 170ab566c..c5a15f065 100644 --- a/test_autolens/imaging/model/test_analysis_imaging.py +++ b/test_autolens/imaging/model/test_analysis_imaging.py @@ -106,6 +106,33 @@ def test__positions__likelihood_overwrites__changes_likelihood(masked_imaging_7x assert analysis_log_likelihood == pytest.approx(-22048700558.9052, 1.0e-4) +def test__positions__likelihood_overwrites__changes_likelihood__double_source_plane_example(masked_imaging_7x7): + + lens = al.Galaxy(redshift=0.5, mass=al.mp.IsothermalSph(centre=(0.05, 0.05))) + source_0 = al.Galaxy(redshift=1.0, light=al.lp.SersicSph(centre=(0.05, 0.05))) + source_1 = al.Galaxy(redshift=2.0, light=al.lp.SersicSph(centre=(0.05, 0.05))) + + model = af.Collection(galaxies=af.Collection(lens=lens, source_0=source_0, source_1=source_1)) + + instance = model.instance_from_unit_vector([]) + + plane_index_positions_dict = { + 1: al.Grid2DIrregular([(1.0, 100.0), (200.0, 2.0)]), + 2: al.Grid2DIrregular([(1.0, 100.0), (200.0, 2.0)]) + } + + positions_likelihood = al.PositionsLHPenalty( + plane_index_positions_dict=plane_index_positions_dict, threshold=0.01 + ) + + analysis = al.AnalysisImaging( + dataset=masked_imaging_7x7, positions_likelihood=positions_likelihood + ) + analysis_log_likelihood = analysis.log_likelihood_function(instance=instance) + + assert analysis_log_likelihood == pytest.approx(-44140499647.28964, 1.0e-4) + + def test__profile_log_likelihood_function(masked_imaging_7x7): pixelization = al.Pixelization( mesh=al.mesh.Rectangular(shape=(3, 3)), From 2278a463b98474eb3c284e67b72b7f598aea2a04 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 21 May 2025 13:05:16 +0100 Subject: [PATCH 03/11] added plane_index_via_redshift_from --- autolens/analysis/result.py | 24 ++++++++++++++++-------- autolens/lens/tracer.py | 25 +++++++++++++++++++++++++ test_autolens/lens/test_tracer.py | 12 ++++++++++++ 3 files changed, 53 insertions(+), 8 deletions(-) diff --git a/autolens/analysis/result.py b/autolens/analysis/result.py index 45a026101..849c969b3 100644 --- a/autolens/analysis/result.py +++ b/autolens/analysis/result.py @@ -1,7 +1,7 @@ import logging import os import numpy as np -from typing import Optional, Union +from typing import List, Optional, Union import autoarray as aa import autogalaxy as ag @@ -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. @@ -179,6 +179,7 @@ def positions_threshold_from( factor=1.0, minimum_threshold=None, positions: Optional[aa.Grid2DIrregular] = None, + plane_index_list : Optional[List[int]] = None, ) -> float: """ Compute a new position threshold from these results corresponding to the image-plane multiple image positions @@ -211,11 +212,18 @@ def positions_threshold_from( by `factor` and rounded up to the `threshold`. """ - positions = ( - self.image_plane_multiple_image_positions - if positions is None - else positions - ) + if plane_index_list is None: + plane_index_list = [-1] + + plane_index_positions_dict = {} + + for plane_index in plane_index_list: + + plane_index_positions_dict[plane_index] = ( + self.image_plane_multiple_image_positions + if positions is None + else positions + ) tracer = Tracer(galaxies=self.max_log_likelihood_galaxies) @@ -281,7 +289,7 @@ def positions_likelihood_from( """ if os.environ.get("PYAUTOFIT_TEST_MODE") == "1": - return None + return positions = ( self.image_plane_multiple_image_positions diff --git a/autolens/lens/tracer.py b/autolens/lens/tracer.py index 9ac72dcae..d9cfe4425 100644 --- a/autolens/lens/tracer.py +++ b/autolens/lens/tracer.py @@ -135,6 +135,31 @@ def planes(self) -> List[ag.Galaxies]: plane_redshifts=self.plane_redshifts, ) + def plane_index_via_redshift_from(self, redshift : float) -> Optional[int]: + """ + Returns the index of a plane at a given redshift. + + This is used to determine the index of a plane in the tracer, which is useful for multi-plane ray-tracing + calculations. The index of the plane may, for example, be used to extract grid of a specific plane in the + tracer after multi-plane ray-tracing calculations. + + A tolerance of 1e-8 is used to determine if the input redshift is close to the redshift of a plane. If + no matching plane is found, None is returned. + + Parameters + ---------- + redshift + The redshift of the plane to find the index of. + + Returns + ------- + The index of the plane that matches the input redshift. + """ + + for plane_index, plane_redshift in enumerate(self.plane_redshifts): + if np.isclose(redshift, plane_redshift, atol=1e-8): + return plane_index + @classmethod def sliced_tracer_from( cls, diff --git a/test_autolens/lens/test_tracer.py b/test_autolens/lens/test_tracer.py index 23e5bda5f..fc3866d54 100644 --- a/test_autolens/lens/test_tracer.py +++ b/test_autolens/lens/test_tracer.py @@ -78,6 +78,18 @@ def test__planes(): assert tracer.planes == [[g1, g1], [g2, g2], [g3]] +def test__plane_index_via_redshift_from(): + + g1 = al.Galaxy(redshift=1) + g2 = al.Galaxy(redshift=2) + g3 = al.Galaxy(redshift=3) + + tracer = al.Tracer(galaxies=[g1, g2, g3]) + + assert tracer.plane_index_via_redshift_from(redshift=2.0) == 1 + assert tracer.plane_index_via_redshift_from(redshift=3.0) == 2 + assert tracer.plane_index_via_redshift_from(redshift=3.001) == None + def test__upper_plane_index_with_light_profile(): g0 = al.Galaxy(redshift=0.5) g1 = al.Galaxy(redshift=1.0) From 75f965ee7e8f24aaa2d40874e32774724be353bf Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 21 May 2025 13:19:02 +0100 Subject: [PATCH 04/11] switch to plane redshift in source max --- autolens/analysis/positions.py | 46 ++++++++++++----- autolens/lens/tracer.py | 50 +++++++++---------- autolens/point/max_separation.py | 14 ++++-- .../imaging/model/test_analysis_imaging.py | 8 +-- 4 files changed, 72 insertions(+), 46 deletions(-) diff --git a/autolens/analysis/positions.py b/autolens/analysis/positions.py index ce69e97cc..f3e5e408c 100644 --- a/autolens/analysis/positions.py +++ b/autolens/analysis/positions.py @@ -26,7 +26,7 @@ def __init__( self, threshold: float, positions: Optional[aa.Grid2DIrregular] = None, - plane_index_positions_dict: Optional[Dict[int, aa.Grid2DIrregular]] = None, + plane_redshift_positions_dict: Optional[Dict[float, aa.Grid2DIrregular]] = None, ): """ The `PositionsLH` objects add a penalty term to the likelihood of the **PyAutoLens** `log_likelihood_function` @@ -41,6 +41,13 @@ def __init__( via `positions=aa.Grid2DIrregular([(1.0, 0.0), (-1.0, 0.0)]` and they do not trace within `threshold=0.3` of one another, the mass model will be rejected and a new model sampled. + The default behaviour assumes a single lens plane and single source plane, meaning the input `positions` + are the image-plane coordinates of one source galaxy at a specifc redshift. + + For multiple source planes, the `plane_redshift_positions_dict` dictionary can be used to input multiple + sets of image positions, where the key is the redshift of the source plane and the value is the + corresponding image-plane coordinates. + Parameters ---------- positions @@ -49,17 +56,21 @@ def __init__( threshold If the maximum separation of any two source plane coordinates is above the threshold the penalty term is applied. + plane_redshift_positions_dict + A dictionary of plane redshifts and the corresponding image-plane coordinates of the lensed source + multiple images. The key is the redshift of the source plane and the value is the corresponding + image-plane coordinates. """ self.positions = positions self.threshold = threshold - self.plane_index_positions_dict = plane_index_positions_dict + self.plane_redshift_positions_dict = plane_redshift_positions_dict - if plane_index_positions_dict is None: - self.plane_index_positions_dict = {-1: positions} + if plane_redshift_positions_dict is None: + self.plane_redshift_positions_dict = {None: positions} - for positions in self.plane_index_positions_dict.values(): + for positions in self.plane_redshift_positions_dict.values(): if len(positions) == 1: raise exc.PositionsException( @@ -93,10 +104,12 @@ def output_positions_info(self, output_path: str, tracer: Tracer): """ with open_(path.join(output_path, "positions.info"), "w+") as f: - for plane_index, positions in self.plane_index_positions_dict.items(): + plane_index = 0 + + for plane_redshift, positions in self.plane_redshift_positions_dict.items(): positions_fit = SourceMaxSeparation( - data=positions, noise_map=None, tracer=tracer + data=positions, noise_map=None, tracer=tracer, plane_redshift=plane_redshift ) distances = positions_fit.data.distances_to_coordinate_from( @@ -104,6 +117,7 @@ def output_positions_info(self, output_path: str, tracer: Tracer): ) f.write(f"Plane Index: {plane_index} \n") + f.write(f"Plane Redshift: {plane_redshift} \n") f.write(f"Positions: \n {positions} \n\n") f.write(f"Radial Distance from (0.0, 0.0): \n {distances} \n\n") f.write(f"Threshold = {self.threshold} \n") @@ -112,6 +126,8 @@ def output_positions_info(self, output_path: str, tracer: Tracer): ) f.write("") + plane_index += 1 + class PositionsLHResample(AbstractPositionsLH): """ @@ -160,13 +176,13 @@ def log_likelihood_function_positions_overwrite( if not tracer.has(cls=ag.mp.MassProfile) or len(tracer.planes) == 1: return - for plane_index, positions in self.plane_index_positions_dict.items(): + for plane_redshift, positions in self.plane_redshift_positions_dict.items(): positions_fit = SourceMaxSeparation( data=positions, noise_map=None, tracer=tracer, - plane_index=plane_index, + plane_redshift=plane_redshift, ) if not positions_fit.max_separation_within_threshold(self.threshold): @@ -182,7 +198,7 @@ def __init__( threshold: float, log_likelihood_penalty_factor: float = 1e8, positions: Optional[aa.Grid2DIrregular] = None, - plane_index_positions_dict: Optional[Dict[int, aa.Grid2DIrregular]] = None, + plane_redshift_positions_dict: Optional[Dict[int, aa.Grid2DIrregular]] = None, ): """ The `PositionsLH` objects add a penalty term to the likelihood of the **PyAutoLens** `log_likelihood_function` @@ -212,11 +228,15 @@ def __init__( log_likelihood_penalty_factor A factor which multiplies how far source pixels do not trace within the threshold of one another, with a larger factor producing a larger penalty making the non-linear parameter space gradient steeper. + plane_redshift_positions_dict + A dictionary of plane redshifts and the corresponding image-plane coordinates of the lensed source + multiple images. The key is the redshift of the source plane and the value is the corresponding + image-plane coordinates. """ super().__init__( positions=positions, threshold=threshold, - plane_index_positions_dict=plane_index_positions_dict, + plane_redshift_positions_dict=plane_redshift_positions_dict, ) self.log_likelihood_penalty_factor = log_likelihood_penalty_factor @@ -295,13 +315,13 @@ def log_likelihood_penalty_from(self, tracer: Tracer) -> Optional[float]: log_likelihood_penalty = 0.0 - for plane_index, positions in self.plane_index_positions_dict.items(): + for plane_redshift, positions in self.plane_redshift_positions_dict.items(): positions_fit = SourceMaxSeparation( data=positions, noise_map=None, tracer=tracer, - plane_index=plane_index, + plane_redshift=plane_redshift, ) if not positions_fit.max_separation_within_threshold(self.threshold): diff --git a/autolens/lens/tracer.py b/autolens/lens/tracer.py index d9cfe4425..07e480c08 100644 --- a/autolens/lens/tracer.py +++ b/autolens/lens/tracer.py @@ -135,31 +135,6 @@ def planes(self) -> List[ag.Galaxies]: plane_redshifts=self.plane_redshifts, ) - def plane_index_via_redshift_from(self, redshift : float) -> Optional[int]: - """ - Returns the index of a plane at a given redshift. - - This is used to determine the index of a plane in the tracer, which is useful for multi-plane ray-tracing - calculations. The index of the plane may, for example, be used to extract grid of a specific plane in the - tracer after multi-plane ray-tracing calculations. - - A tolerance of 1e-8 is used to determine if the input redshift is close to the redshift of a plane. If - no matching plane is found, None is returned. - - Parameters - ---------- - redshift - The redshift of the plane to find the index of. - - Returns - ------- - The index of the plane that matches the input redshift. - """ - - for plane_index, plane_redshift in enumerate(self.plane_redshifts): - if np.isclose(redshift, plane_redshift, atol=1e-8): - return plane_index - @classmethod def sliced_tracer_from( cls, @@ -826,6 +801,31 @@ def cls_list_from(self, cls: Type) -> List: return cls_list + def plane_index_via_redshift_from(self, redshift : float) -> Optional[int]: + """ + Returns the index of a plane at a given redshift. + + This is used to determine the index of a plane in the tracer, which is useful for multi-plane ray-tracing + calculations. The index of the plane may, for example, be used to extract grid of a specific plane in the + tracer after multi-plane ray-tracing calculations. + + A tolerance of 1e-8 is used to determine if the input redshift is close to the redshift of a plane. If + no matching plane is found, None is returned. + + Parameters + ---------- + redshift + The redshift of the plane to find the index of. + + Returns + ------- + The index of the plane that matches the input redshift. + """ + + for plane_index, plane_redshift in enumerate(self.plane_redshifts): + if np.isclose(redshift, plane_redshift, atol=1e-8): + return plane_index + @property def plane_indexes_with_pixelizations(self) -> List[int]: """ diff --git a/autolens/point/max_separation.py b/autolens/point/max_separation.py index f6fcfa591..957c91153 100644 --- a/autolens/point/max_separation.py +++ b/autolens/point/max_separation.py @@ -11,7 +11,7 @@ def __init__( data: aa.Grid2DIrregular, noise_map: Optional[aa.ArrayIrregular], tracer: Tracer, - plane_index: int = -1, + plane_redshift: float = Optional[None], ): """ Given a positions dataset, which is a list of positions with names that associated them to model source @@ -29,13 +29,19 @@ def __init__( The object that defines the ray-tracing of the strong lens system of galaxies. noise_value The noise-value assumed when computing the log likelihood. - plane_index - The index of the plane in the `Tracer` that the source-plane positions are computed from. This is typically - the last plane in the `Tracer`, which is the source-plane. + plane_redshift + The redshift of the plane in the `Tracer` that the source-plane positions are computed from. This is + often the last plane in the `Tracer`, which is the source-plane. """ self.data = data self.noise_map = noise_map + + if plane_redshift is None: + plane_index = -1 + else: + plane_index = tracer.plane_index_via_redshift_from(redshift=plane_redshift) + self.source_plane_positions = tracer.traced_grid_2d_list_from(grid=data)[ plane_index ] diff --git a/test_autolens/imaging/model/test_analysis_imaging.py b/test_autolens/imaging/model/test_analysis_imaging.py index c5a15f065..8e4be1ec8 100644 --- a/test_autolens/imaging/model/test_analysis_imaging.py +++ b/test_autolens/imaging/model/test_analysis_imaging.py @@ -116,13 +116,13 @@ def test__positions__likelihood_overwrites__changes_likelihood__double_source_pl instance = model.instance_from_unit_vector([]) - plane_index_positions_dict = { - 1: al.Grid2DIrregular([(1.0, 100.0), (200.0, 2.0)]), - 2: al.Grid2DIrregular([(1.0, 100.0), (200.0, 2.0)]) + plane_redshift_positions_dict = { + 1.0: al.Grid2DIrregular([(1.0, 100.0), (200.0, 2.0)]), + 2.0: al.Grid2DIrregular([(1.0, 100.0), (200.0, 2.0)]) } positions_likelihood = al.PositionsLHPenalty( - plane_index_positions_dict=plane_index_positions_dict, threshold=0.01 + plane_redshift_positions_dict=plane_redshift_positions_dict, threshold=0.01 ) analysis = al.AnalysisImaging( From cdc8be5a00ca87d625dc52784795dcda3c8a7701 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 21 May 2025 13:28:01 +0100 Subject: [PATCH 05/11] refactor point naming --- autolens/point/fit/abstract.py | 6 +-- .../point/fit/positions/image/abstract.py | 2 +- .../point/fit/positions/source/separations.py | 2 +- autolens/point/max_separation.py | 16 +++---- autolens/point/mock/mock_solver.py | 2 +- autolens/point/solver/point_solver.py | 11 +++-- autolens/point/solver/shape_solver.py | 48 +++++++++---------- autolens/point/solver/step.py | 2 +- 8 files changed, 45 insertions(+), 44 deletions(-) diff --git a/autolens/point/fit/abstract.py b/autolens/point/fit/abstract.py index ade14c149..2b5ee9a4d 100644 --- a/autolens/point/fit/abstract.py +++ b/autolens/point/fit/abstract.py @@ -145,7 +145,7 @@ def source_plane_coordinate(self) -> Tuple[float, float]: return self.profile.centre @property - def source_plane_index(self) -> int: + def plane_index(self) -> int: """ Returns the integer plane index containing the point source galaxy, which is used when computing the deflection angles of image-plane positions from the tracer. @@ -160,7 +160,7 @@ def source_plane_index(self) -> int: return self.tracer.extract_plane_index_of_profile(profile_name=self.name) @property - def source_plane_redshift(self) -> float: + def plane_redshift(self) -> float: """ Returns the redshift of the plane containing the point source galaxy, which is used when computing the deflection angles of image-plane positions from the tracer. @@ -173,4 +173,4 @@ def source_plane_redshift(self) -> float: ------- The redshift of the plane containing the point-source galaxy. """ - return self.tracer.planes[self.source_plane_index].redshift + return self.tracer.planes[self.plane_index].redshift diff --git a/autolens/point/fit/positions/image/abstract.py b/autolens/point/fit/positions/image/abstract.py index 5ab480dc0..9a0597d62 100644 --- a/autolens/point/fit/positions/image/abstract.py +++ b/autolens/point/fit/positions/image/abstract.py @@ -101,5 +101,5 @@ def model_data(self) -> aa.Grid2DIrregular: return self.solver.solve( tracer=self.tracer, source_plane_coordinate=self.source_plane_coordinate, - source_plane_redshift=self.source_plane_redshift, + plane_redshift=self.plane_redshift, ) diff --git a/autolens/point/fit/positions/source/separations.py b/autolens/point/fit/positions/source/separations.py index 2c35269b8..95ec8c952 100644 --- a/autolens/point/fit/positions/source/separations.py +++ b/autolens/point/fit/positions/source/separations.py @@ -95,7 +95,7 @@ def model_data(self) -> aa.Grid2DIrregular: deflections = self.tracer.deflections_yx_2d_from(grid=self.data) else: deflections = self.tracer.deflections_between_planes_from( - grid=self.data, plane_i=0, plane_j=self.source_plane_index + grid=self.data, plane_i=0, plane_j=self.plane_index ) return self.data.grid_2d_via_deflection_grid_from(deflection_grid=deflections) diff --git a/autolens/point/max_separation.py b/autolens/point/max_separation.py index 957c91153..2cb19a2ce 100644 --- a/autolens/point/max_separation.py +++ b/autolens/point/max_separation.py @@ -42,23 +42,23 @@ def __init__( else: plane_index = tracer.plane_index_via_redshift_from(redshift=plane_redshift) - self.source_plane_positions = tracer.traced_grid_2d_list_from(grid=data)[ + self.plane_positions = tracer.traced_grid_2d_list_from(grid=data)[ plane_index ] @property - def furthest_separations_of_source_plane_positions(self) -> aa.ArrayIrregular: + def furthest_separations_of_plane_positions(self) -> aa.ArrayIrregular: """ Returns the furthest distance of every source-plane (y,x) coordinate to the other source-plane (y,x) coordinates. For example, for the following source-plane positions: - source_plane_positions = [[(0.0, 0.0), (0.0, 1.0), (0.0, 3.0)] + plane_positions = [[(0.0, 0.0), (0.0, 1.0), (0.0, 3.0)] The returned furthest distances are: - source_plane_positions = [3.0, 2.0, 3.0] + plane_positions = [3.0, 2.0, 3.0] Returns ------- @@ -66,11 +66,11 @@ def furthest_separations_of_source_plane_positions(self) -> aa.ArrayIrregular: The further distances of every set of grouped source-plane coordinates the other source-plane coordinates that it is grouped with. """ - return self.source_plane_positions.furthest_distances_to_other_coordinates + return self.plane_positions.furthest_distances_to_other_coordinates @property - def max_separation_of_source_plane_positions(self) -> float: - return max(self.furthest_separations_of_source_plane_positions) + def max_separation_of_plane_positions(self) -> float: + return max(self.furthest_separations_of_plane_positions) def max_separation_within_threshold(self, threshold) -> bool: - return self.max_separation_of_source_plane_positions <= threshold + return self.max_separation_of_plane_positions <= threshold diff --git a/autolens/point/mock/mock_solver.py b/autolens/point/mock/mock_solver.py index 4e1e5d062..16366fa01 100644 --- a/autolens/point/mock/mock_solver.py +++ b/autolens/point/mock/mock_solver.py @@ -9,6 +9,6 @@ def solve( self, tracer, source_plane_coordinate, - source_plane_redshift: Optional[float] = None, + plane_redshift: Optional[float] = None, ): return self.model_positions diff --git a/autolens/point/solver/point_solver.py b/autolens/point/solver/point_solver.py index 2d160d583..f98d5a129 100644 --- a/autolens/point/solver/point_solver.py +++ b/autolens/point/solver/point_solver.py @@ -22,7 +22,7 @@ def solve( self, tracer: OperateDeflections, source_plane_coordinate: Tuple[float, float], - source_plane_redshift: Optional[float] = None, + plane_redshift: Optional[float] = None, ) -> aa.Grid2DIrregular: """ Solve for the image plane coordinates that are traced to the source plane coordinate. @@ -37,11 +37,12 @@ def solve( Parameters ---------- source_plane_coordinate - The source plane coordinate to trace to the image plane. + The plane coordinate to trace to the image plane, which by default in the source-plane coordinate + but could be a coordinate in another plane is `plane_redshift` is input. tracer The tracer that traces the image plane coordinates to the source plane - source_plane_redshift - The redshift of the source plane coordinate. + plane_redshift + The redshift of the plane coordinate, which for multi-plane systems may not be the source-plane. Returns ------- @@ -50,7 +51,7 @@ def solve( kept_triangles = super().solve_triangles( tracer=tracer, shape=Point(*source_plane_coordinate), - source_plane_redshift=source_plane_redshift, + plane_redshift=plane_redshift, ) filtered_means = self._filter_low_magnification( diff --git a/autolens/point/solver/shape_solver.py b/autolens/point/solver/shape_solver.py index f398f41cd..64547a9aa 100644 --- a/autolens/point/solver/shape_solver.py +++ b/autolens/point/solver/shape_solver.py @@ -186,10 +186,10 @@ def n_steps(self) -> int: return math.ceil(math.log2(self.scale / self.pixel_scale_precision)) @staticmethod - def _source_plane_grid( + def _plane_grid( tracer: OperateDeflections, grid: aa.type.Grid2DLike, - source_plane_redshift: Optional[float] = None, + plane_redshift: Optional[float] = None, ) -> aa.type.Grid2DLike: """ Calculate the source plane grid from the image plane grid. @@ -204,16 +204,16 @@ def _source_plane_grid( The source plane grid computed by applying the deflections to the image plane grid. """ - source_plane_index = -1 + plane_index = -1 - if source_plane_redshift is not None: + if plane_redshift is not None: for redshift in tracer.plane_redshifts: - source_plane_index += 1 - if redshift == source_plane_redshift: + plane_index += 1 + if redshift == plane_redshift: break deflections = tracer.deflections_between_planes_from( - grid=grid, plane_i=0, plane_j=source_plane_index + grid=grid, plane_i=0, plane_j=plane_index ) # noinspection PyTypeChecker return grid.grid_2d_via_deflection_grid_from(deflection_grid=deflections) @@ -223,7 +223,7 @@ def solve_triangles( self, tracer: OperateDeflections, shape: Shape, - source_plane_redshift: Optional[float] = None, + plane_redshift: Optional[float] = None, ) -> AbstractTriangles: """ Solve for the image plane coordinates that are traced to the source plane coordinate. @@ -241,7 +241,7 @@ def solve_triangles( The tracer to use to trace the image plane coordinates to the source plane. shape The shape in the source plane for which we want to identify the image plane coordinates. - source_plane_redshift + plane_redshift The redshift of the source plane. Returns @@ -257,7 +257,7 @@ def solve_triangles( self.steps( tracer=tracer, shape=shape, - source_plane_redshift=source_plane_redshift, + plane_redshift=plane_redshift, ) ) final_step = steps[-1] @@ -286,27 +286,27 @@ def _filter_low_magnification( mask = np.abs(magnifications.array) > self.magnification_threshold return np.where(mask[:, None], points, np.nan) - def _source_triangles( + def _plane_triangles( self, tracer: OperateDeflections, triangles: aa.AbstractTriangles, - source_plane_redshift, + plane_redshift, ): """ Filter the triangles to keep only those that meet the solver condition """ - source_plane_grid = self._source_plane_grid( + plane_grid = self._plane_grid( tracer=tracer, grid=aa.Grid2DIrregular(triangles.vertices), - source_plane_redshift=source_plane_redshift, + plane_redshift=plane_redshift, ) - return triangles.with_vertices(source_plane_grid.array) + return triangles.with_vertices(plane_grid.array) def steps( self, tracer: OperateDeflections, shape: Shape, - source_plane_redshift: Optional[float] = None, + plane_redshift: Optional[float] = None, ) -> Iterator[Step]: """ Iterate over the steps of the triangle solver algorithm. @@ -315,7 +315,7 @@ def steps( ---------- tracer The tracer to use to trace the image plane coordinates to the source plane. - source_plane_redshift + plane_redshift The redshift of the source plane. shape The shape in the source plane for which we want to identify the image plane coordinates. @@ -326,13 +326,13 @@ def steps( """ initial_triangles = self.initial_triangles for number in range(self.n_steps): - source_triangles = self._source_triangles( + plane_triangles = self._plane_triangles( tracer=tracer, triangles=initial_triangles, - source_plane_redshift=source_plane_redshift, + plane_redshift=plane_redshift, ) - indexes = source_triangles.containing_indices(shape=shape) + indexes = plane_triangles.containing_indices(shape=shape) kept_triangles = initial_triangles.for_indexes(indexes=indexes) neighbourhood = kept_triangles @@ -347,7 +347,7 @@ def steps( filtered_triangles=kept_triangles, neighbourhood=neighbourhood, up_sampled=up_sampled, - source_triangles=source_triangles, + plane_triangles=plane_triangles, ) initial_triangles = up_sampled @@ -376,7 +376,7 @@ def find_magnification( self, tracer: OperateDeflections, shape: Shape, - source_plane_redshift: Optional[float] = None, + plane_redshift: Optional[float] = None, ) -> float: """ Find the magnification of the shape in the source plane. @@ -387,7 +387,7 @@ def find_magnification( A tracer that traces the image plane to the source plane. shape The shape of an image plane pixel. - source_plane_redshift + plane_redshift The redshift of the source plane. Returns @@ -397,6 +397,6 @@ def find_magnification( kept_triangles = super().solve_triangles( tracer=tracer, shape=shape, - source_plane_redshift=source_plane_redshift, + plane_redshift=plane_redshift, ) return kept_triangles.area / shape.area diff --git a/autolens/point/solver/step.py b/autolens/point/solver/step.py index c25909677..27a8ae779 100644 --- a/autolens/point/solver/step.py +++ b/autolens/point/solver/step.py @@ -38,7 +38,7 @@ class Step: filtered_triangles: aa.AbstractTriangles neighbourhood: aa.AbstractTriangles up_sampled: aa.AbstractTriangles - source_triangles: aa.AbstractTriangles + plane_triangles: aa.AbstractTriangles def tree_flatten(self): return ( From 0d878df9476f3797431dc03b7753653fd06df95e Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 21 May 2025 14:03:48 +0100 Subject: [PATCH 06/11] mostly implemented but now making an API change --- autolens/analysis/positions.py | 7 ++- autolens/analysis/result.py | 85 +++++++++++++++------------ autolens/point/max_separation.py | 8 ++- autolens/point/solver/shape_solver.py | 11 ++-- test_autolens/analysis/test_result.py | 2 +- 5 files changed, 64 insertions(+), 49 deletions(-) diff --git a/autolens/analysis/positions.py b/autolens/analysis/positions.py index f3e5e408c..8f4cfa43b 100644 --- a/autolens/analysis/positions.py +++ b/autolens/analysis/positions.py @@ -70,6 +70,9 @@ def __init__( if plane_redshift_positions_dict is None: self.plane_redshift_positions_dict = {None: positions} + if len(self.plane_redshift_positions_dict) == 1: + self.positions = list(self.plane_redshift_positions_dict.values())[0] + for positions in self.plane_redshift_positions_dict.values(): if len(positions) == 1: @@ -122,7 +125,7 @@ def output_positions_info(self, output_path: str, tracer: Tracer): f.write(f"Radial Distance from (0.0, 0.0): \n {distances} \n\n") f.write(f"Threshold = {self.threshold} \n") f.write( - f"Max Source Plane Separation of Maximum Likelihood Model = {positions_fit.max_separation_of_source_plane_positions}" + f"Max Source Plane Separation of Maximum Likelihood Model = {positions_fit.max_separation_of_plane_positions}" ) f.write("") @@ -327,7 +330,7 @@ def log_likelihood_penalty_from(self, tracer: Tracer) -> Optional[float]: if not positions_fit.max_separation_within_threshold(self.threshold): log_likelihood_penalty += self.log_likelihood_penalty_factor * ( - positions_fit.max_separation_of_source_plane_positions + positions_fit.max_separation_of_plane_positions - self.threshold ) diff --git a/autolens/analysis/result.py b/autolens/analysis/result.py index 849c969b3..cfb0b64d7 100644 --- a/autolens/analysis/result.py +++ b/autolens/analysis/result.py @@ -53,7 +53,7 @@ def max_log_likelihood_positions_threshold(self) -> float: tracer=self.max_log_likelihood_tracer, ) - return positions_fits.max_separation_of_source_plane_positions + return positions_fits.max_separation_of_plane_positions @property def source_plane_light_profile_centre(self) -> aa.Grid2DIrregular: @@ -83,8 +83,7 @@ def source_plane_centre(self) -> aa.Grid2DIrregular: """ return self.source_plane_light_profile_centre - @property - def image_plane_multiple_image_positions(self) -> aa.Grid2DIrregular: + def image_plane_multiple_image_positions(self, plane_redshift : Optional[float] = None) -> aa.Grid2DIrregular: """ Backwards ray-trace the source-plane centres (see above) to the image-plane via the mass model, to determine the multiple image position of the source(s) in the image-plane. @@ -102,6 +101,7 @@ def image_plane_multiple_image_positions(self) -> aa.Grid2DIrregular: multiple_images = solver.solve( tracer=self.max_log_likelihood_tracer, source_plane_coordinate=self.source_plane_centre.in_list[0], + plane_redshift=plane_redshift, ) if multiple_images.shape[0] == 1: @@ -110,7 +110,7 @@ def image_plane_multiple_image_positions(self) -> aa.Grid2DIrregular: return aa.Grid2DIrregular(values=multiple_images) def image_plane_multiple_image_positions_for_single_image_from( - self, increments: int = 20 + self, plane_redshift : Optional[float] = None, increments: int = 20 ) -> aa.Grid2DIrregular: """ If the standard point solver only locates one multiple image, finds one or more additional images, which are @@ -159,6 +159,7 @@ def image_plane_multiple_image_positions_for_single_image_from( multiple_images = solver.solve( tracer=self.max_log_likelihood_tracer, source_plane_coordinate=(centre[0] * factor, centre[1] * factor), + plane_redshift=plane_redshift, ) if multiple_images.shape[0] > 1: @@ -179,7 +180,7 @@ def positions_threshold_from( factor=1.0, minimum_threshold=None, positions: Optional[aa.Grid2DIrregular] = None, - plane_index_list : Optional[List[int]] = None, + plane_redshift : Optional[float] = None, ) -> float: """ Compute a new position threshold from these results corresponding to the image-plane multiple image positions @@ -192,6 +193,10 @@ def positions_threshold_from( This is used for non-linear search chaining, specifically updating the position threshold of a new model-fit using the maximum likelihood model of a previous search. + + The above behaviour assumes a single lens and single source plane, if there are multiple source planes (e.g. + a double source plane lens) the input `plane_redshift_list` can be used to specify which source planea + the position threshold is computed for. Parameters ---------- @@ -211,28 +216,22 @@ def positions_threshold_from( The maximum source plane separation of this results maximum likelihood `Tracer` multiple images multiplied by `factor` and rounded up to the `threshold`. """ - - if plane_index_list is None: - plane_index_list = [-1] - - plane_index_positions_dict = {} - - for plane_index in plane_index_list: - - plane_index_positions_dict[plane_index] = ( - self.image_plane_multiple_image_positions - if positions is None - else positions - ) + plane_redshift_positions_dict = {} tracer = Tracer(galaxies=self.max_log_likelihood_galaxies) + positions = ( + self.image_plane_multiple_image_positions(plane_redshift=plane_redshift) + if positions is None + else positions + ) + positions_fits = SourceMaxSeparation( - data=positions, noise_map=None, tracer=tracer + data=positions, noise_map=None, tracer=tracer, plane_redshift=plane_redshift ) threshold = factor * np.max( - positions_fits.max_separation_of_source_plane_positions + positions_fits.max_separation_of_plane_positions ) if minimum_threshold is not None: @@ -248,6 +247,7 @@ def positions_likelihood_from( use_resample=False, positions: Optional[aa.Grid2DIrregular] = None, mass_centre_radial_distance_min: float = None, + plane_redshift_list: Optional[List[int]] = None, ) -> Union[PositionsLHPenalty, PositionsLHResample]: """ Returns a `PositionsLH` object from the result of a lens model-fit, where the maximum log likelihood mass @@ -291,30 +291,43 @@ def positions_likelihood_from( if os.environ.get("PYAUTOFIT_TEST_MODE") == "1": return - positions = ( - self.image_plane_multiple_image_positions - if positions is None - else positions - ) + plane_redshift_list = [None] if plane_redshift_list is None else plane_redshift_list + + plane_redshift_positions_dict = {} + threshold_list = [] - 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" + for plane_redshift in plane_redshift_list: + + positions = ( + self.image_plane_multiple_image_positions(plane_redshift=plane_redshift) + if positions is None + else positions ) - distances = positions.distances_to_coordinate_from( - coordinate=mass_centre[0] + 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, plane_redshift=plane_redshift ) - positions = positions[distances > mass_centre_radial_distance_min] + threshold_list.append(threshold) - threshold = self.positions_threshold_from( - factor=factor, minimum_threshold=minimum_threshold, positions=positions - ) + plane_redshift_positions_dict[plane_redshift] = positions + + threshold = np.max(threshold_list) if not use_resample: - return PositionsLHPenalty(positions=positions, threshold=threshold) - return PositionsLHResample(positions=positions, threshold=threshold) + return PositionsLHPenalty(plane_redshift_positions_dict=plane_redshift_positions_dict, threshold=threshold) + return PositionsLHResample(plane_redshift_positions_dict=plane_redshift_positions_dict, threshold=threshold) class ResultDataset(Result): diff --git a/autolens/point/max_separation.py b/autolens/point/max_separation.py index 2cb19a2ce..20c8f6ae7 100644 --- a/autolens/point/max_separation.py +++ b/autolens/point/max_separation.py @@ -37,10 +37,12 @@ def __init__( self.data = data self.noise_map = noise_map - if plane_redshift is None: - plane_index = -1 - else: + print(plane_redshift) + + try: plane_index = tracer.plane_index_via_redshift_from(redshift=plane_redshift) + except TypeError: + plane_index = -1 self.plane_positions = tracer.traced_grid_2d_list_from(grid=data)[ plane_index diff --git a/autolens/point/solver/shape_solver.py b/autolens/point/solver/shape_solver.py index 64547a9aa..e3fb06e44 100644 --- a/autolens/point/solver/shape_solver.py +++ b/autolens/point/solver/shape_solver.py @@ -204,13 +204,10 @@ def _plane_grid( The source plane grid computed by applying the deflections to the image plane grid. """ - plane_index = -1 - - if plane_redshift is not None: - for redshift in tracer.plane_redshifts: - plane_index += 1 - if redshift == plane_redshift: - break + if plane_redshift is None: + plane_index = -1 + else: + plane_index = tracer.plane_index_via_redshift_from(redshift=plane_redshift) deflections = tracer.deflections_between_planes_from( grid=grid, plane_i=0, plane_j=plane_index diff --git a/test_autolens/analysis/test_result.py b/test_autolens/analysis/test_result.py index 34917ccfc..32416c150 100644 --- a/test_autolens/analysis/test_result.py +++ b/test_autolens/analysis/test_result.py @@ -225,7 +225,7 @@ def test__image_plane_multiple_image_positions(analysis_imaging_7x7): samples_summary=samples_summary, analysis=analysis_imaging_7x7 ) - multiple_images = result.image_plane_multiple_image_positions + multiple_images = result.image_plane_multiple_image_positions() assert pytest.approx((0.968719, 0.366210), 1.0e-2) in multiple_images.in_list From ca5023cd609ee353f1e64170a25a769e0c7f41a2 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 21 May 2025 14:17:32 +0100 Subject: [PATCH 07/11] update to positions_likelihood_list --- autolens/analysis/analysis/dataset.py | 18 +-- autolens/analysis/analysis/lens.py | 19 +-- autolens/analysis/positions.py | 129 +++++++----------- .../imaging/model/test_analysis_imaging.py | 18 ++- 4 files changed, 81 insertions(+), 103 deletions(-) diff --git a/autolens/analysis/analysis/dataset.py b/autolens/analysis/analysis/dataset.py index 68e4b9535..654c9aaf4 100644 --- a/autolens/analysis/analysis/dataset.py +++ b/autolens/analysis/analysis/dataset.py @@ -1,9 +1,9 @@ import os import logging -from typing import Optional, Union +from typing import List, Optional, Union from autoconf import conf -from autoconf.dictable import to_dict, output_to_json +from autoconf.dictable import output_to_json import autofit as af import autoarray as aa @@ -27,8 +27,8 @@ class AnalysisDataset(AgAnalysisDataset, AnalysisLens): def __init__( self, dataset, - positions_likelihood: Optional[ - Union[PositionsLHResample, PositionsLHPenalty] + positions_likelihood_list: Optional[ + List[Union[PositionsLHResample, PositionsLHPenalty]] ] = None, adapt_image_maker: Optional[ag.AdaptImageMaker] = None, cosmology: ag.cosmo.LensingCosmology = ag.cosmo.Planck15(), @@ -49,10 +49,10 @@ def __init__( ---------- dataset The imaging, interferometer or other dataset that the model if fitted too. - positions_likelihood - An object which alters the likelihood function to include a term which accounts for whether - image-pixel coordinates in arc-seconds corresponding to the multiple images of the lensed source galaxy - trace close to one another in the source-plane. + positions_likelihood_list + Alters the likelihood function to include a term which accounts for whether image-pixel coordinates in + arc-seconds corresponding to the multiple images of each lensed source galaxy trace close to one another in + their source-plane. adapt_images Contains the adapt-images which are used to make a pixelization's mesh and regularization adapt to the reconstructed galaxy's morphology. @@ -78,7 +78,7 @@ def __init__( AnalysisLens.__init__( self=self, - positions_likelihood=positions_likelihood, + positions_likelihood_list=positions_likelihood_list, cosmology=cosmology, ) diff --git a/autolens/analysis/analysis/lens.py b/autolens/analysis/analysis/lens.py index 41ac474e3..78fd223bc 100644 --- a/autolens/analysis/analysis/lens.py +++ b/autolens/analysis/analysis/lens.py @@ -1,6 +1,6 @@ import logging import numpy as np -from typing import Dict, Optional, Union +from typing import Dict, List, Optional, Union import autofit as af import autoarray as aa @@ -22,8 +22,8 @@ class AnalysisLens: def __init__( self, - positions_likelihood: Optional[ - Union[PositionsLHResample, PositionsLHPenalty] + positions_likelihood_list: Optional[ + List[Union[PositionsLHResample, PositionsLHPenalty]] ] = None, cosmology: ag.cosmo.LensingCosmology = ag.cosmo.Planck15(), ): @@ -43,7 +43,7 @@ def __init__( The Cosmology assumed for this analysis. """ self.cosmology = cosmology - self.positions_likelihood = positions_likelihood + self.positions_likelihood_list = positions_likelihood_list def tracer_via_instance_from( self, @@ -123,10 +123,13 @@ def log_likelihood_positions_overwrite_from( The penalty value of the positions log likelihood, if the positions do not trace close in the source plane, else a None is returned to indicate there is no penalty. """ - if self.positions_likelihood is not None: + if self.positions_likelihood_list is not None: try: - return self.positions_likelihood.log_likelihood_function_positions_overwrite( - instance=instance, analysis=self - ) + log_likelihood = 0.0 + for positions_likelihood in self.positions_likelihood_list: + log_likelihood += positions_likelihood.log_likelihood_function_positions_overwrite( + instance=instance, analysis=self + ) + return log_likelihood except (ValueError, np.linalg.LinAlgError) as e: raise exc.FitException from e diff --git a/autolens/analysis/positions.py b/autolens/analysis/positions.py index 8f4cfa43b..0b2b5b8b5 100644 --- a/autolens/analysis/positions.py +++ b/autolens/analysis/positions.py @@ -26,7 +26,7 @@ def __init__( self, threshold: float, positions: Optional[aa.Grid2DIrregular] = None, - plane_redshift_positions_dict: Optional[Dict[float, aa.Grid2DIrregular]] = None, + plane_redshift: Optional[float] = None, ): """ The `PositionsLH` objects add a penalty term to the likelihood of the **PyAutoLens** `log_likelihood_function` @@ -44,9 +44,8 @@ def __init__( The default behaviour assumes a single lens plane and single source plane, meaning the input `positions` are the image-plane coordinates of one source galaxy at a specifc redshift. - For multiple source planes, the `plane_redshift_positions_dict` dictionary can be used to input multiple - sets of image positions, where the key is the redshift of the source plane and the value is the - corresponding image-plane coordinates. + For multiple source planes, the `plane_redshift` can be used to pair the image-plane positions with the + redshift of the source plane they trace too. Parameters ---------- @@ -56,31 +55,21 @@ def __init__( threshold If the maximum separation of any two source plane coordinates is above the threshold the penalty term is applied. - plane_redshift_positions_dict - A dictionary of plane redshifts and the corresponding image-plane coordinates of the lensed source - multiple images. The key is the redshift of the source plane and the value is the corresponding - image-plane coordinates. + plane_redshift + The plane redshift of the lensed source multiple images, which is only required if position threshold + for a double source plane lens system is being used where the specific plane is required. """ self.positions = positions self.threshold = threshold + self.plane_redshift = plane_redshift - self.plane_redshift_positions_dict = plane_redshift_positions_dict - - if plane_redshift_positions_dict is None: - self.plane_redshift_positions_dict = {None: positions} - - if len(self.plane_redshift_positions_dict) == 1: - self.positions = list(self.plane_redshift_positions_dict.values())[0] - - for positions in self.plane_redshift_positions_dict.values(): - - if len(positions) == 1: - raise exc.PositionsException( - f"The positions input into the PositionsLikelihood object have length one " - f"(e.g. it is only one (y,x) coordinate and therefore cannot be compared with other images).\n\n" - "Please input more positions into the Positions." - ) + if len(positions) == 1: + raise exc.PositionsException( + f"The positions input into the PositionsLikelihood object have length one " + f"(e.g. it is only one (y,x) coordinate and therefore cannot be compared with other images).\n\n" + "Please input more positions into the Positions." + ) def log_likelihood_function_positions_overwrite( self, instance: af.ModelInstance, analysis: AnalysisDataset @@ -107,29 +96,22 @@ def output_positions_info(self, output_path: str, tracer: Tracer): """ with open_(path.join(output_path, "positions.info"), "w+") as f: - plane_index = 0 - - for plane_redshift, positions in self.plane_redshift_positions_dict.items(): - - positions_fit = SourceMaxSeparation( - data=positions, noise_map=None, tracer=tracer, plane_redshift=plane_redshift - ) - - distances = positions_fit.data.distances_to_coordinate_from( - coordinate=(0.0, 0.0) - ) + positions_fit = SourceMaxSeparation( + data=self.positions, noise_map=None, tracer=tracer, plane_redshift=self.plane_redshift + ) - f.write(f"Plane Index: {plane_index} \n") - f.write(f"Plane Redshift: {plane_redshift} \n") - f.write(f"Positions: \n {positions} \n\n") - f.write(f"Radial Distance from (0.0, 0.0): \n {distances} \n\n") - f.write(f"Threshold = {self.threshold} \n") - f.write( - f"Max Source Plane Separation of Maximum Likelihood Model = {positions_fit.max_separation_of_plane_positions}" - ) - f.write("") + distances = positions_fit.data.distances_to_coordinate_from( + coordinate=(0.0, 0.0) + ) - plane_index += 1 + f.write(f"Plane Redshift: {self.plane_redshift} \n") + f.write(f"Positions: \n {self.positions} \n\n") + f.write(f"Radial Distance from (0.0, 0.0): \n {distances} \n\n") + f.write(f"Threshold = {self.threshold} \n") + f.write( + f"Max Source Plane Separation of Maximum Likelihood Model = {positions_fit.max_separation_of_plane_positions}" + ) + f.write("") class PositionsLHResample(AbstractPositionsLH): @@ -179,20 +161,18 @@ def log_likelihood_function_positions_overwrite( if not tracer.has(cls=ag.mp.MassProfile) or len(tracer.planes) == 1: return - for plane_redshift, positions in self.plane_redshift_positions_dict.items(): - - positions_fit = SourceMaxSeparation( - data=positions, - noise_map=None, - tracer=tracer, - plane_redshift=plane_redshift, - ) + positions_fit = SourceMaxSeparation( + data=self.positions, + noise_map=None, + tracer=tracer, + plane_redshift=self.plane_redshift, + ) - if not positions_fit.max_separation_within_threshold(self.threshold): - if os.environ.get("PYAUTOFIT_TEST_MODE") == "1": - return + if not positions_fit.max_separation_within_threshold(self.threshold): + if os.environ.get("PYAUTOFIT_TEST_MODE") == "1": + return - raise exc.RayTracingException + raise exc.RayTracingException class PositionsLHPenalty(AbstractPositionsLH): @@ -201,7 +181,7 @@ def __init__( threshold: float, log_likelihood_penalty_factor: float = 1e8, positions: Optional[aa.Grid2DIrregular] = None, - plane_redshift_positions_dict: Optional[Dict[int, aa.Grid2DIrregular]] = None, + plane_redshift: Optional[float] = None, ): """ The `PositionsLH` objects add a penalty term to the likelihood of the **PyAutoLens** `log_likelihood_function` @@ -231,15 +211,14 @@ def __init__( log_likelihood_penalty_factor A factor which multiplies how far source pixels do not trace within the threshold of one another, with a larger factor producing a larger penalty making the non-linear parameter space gradient steeper. - plane_redshift_positions_dict - A dictionary of plane redshifts and the corresponding image-plane coordinates of the lensed source - multiple images. The key is the redshift of the source plane and the value is the corresponding - image-plane coordinates. + plane_redshift + The plane redshift of the lensed source multiple images, which is only required if position threshold + for a double source plane lens system is being used where the specific plane is required. """ super().__init__( positions=positions, threshold=threshold, - plane_redshift_positions_dict=plane_redshift_positions_dict, + plane_redshift=plane_redshift ) self.log_likelihood_penalty_factor = log_likelihood_penalty_factor @@ -318,21 +297,19 @@ def log_likelihood_penalty_from(self, tracer: Tracer) -> Optional[float]: log_likelihood_penalty = 0.0 - for plane_redshift, positions in self.plane_redshift_positions_dict.items(): - - positions_fit = SourceMaxSeparation( - data=positions, - noise_map=None, - tracer=tracer, - plane_redshift=plane_redshift, - ) + positions_fit = SourceMaxSeparation( + data=self.positions, + noise_map=None, + tracer=tracer, + plane_redshift=self.plane_redshift, + ) - if not positions_fit.max_separation_within_threshold(self.threshold): + if not positions_fit.max_separation_within_threshold(self.threshold): - log_likelihood_penalty += self.log_likelihood_penalty_factor * ( - positions_fit.max_separation_of_plane_positions - - self.threshold - ) + log_likelihood_penalty += self.log_likelihood_penalty_factor * ( + positions_fit.max_separation_of_plane_positions + - self.threshold + ) return log_likelihood_penalty diff --git a/test_autolens/imaging/model/test_analysis_imaging.py b/test_autolens/imaging/model/test_analysis_imaging.py index 8e4be1ec8..108d4ddea 100644 --- a/test_autolens/imaging/model/test_analysis_imaging.py +++ b/test_autolens/imaging/model/test_analysis_imaging.py @@ -57,7 +57,7 @@ def test__positions__resample__raises_exception(masked_imaging_7x7): ) analysis = al.AnalysisImaging( - dataset=masked_imaging_7x7, positions_likelihood=positions_likelihood + dataset=masked_imaging_7x7, positions_likelihood_list=[positions_likelihood] ) instance = model.instance_from_unit_vector([]) @@ -89,7 +89,7 @@ def test__positions__likelihood_overwrites__changes_likelihood(masked_imaging_7x ) analysis = al.AnalysisImaging( - dataset=masked_imaging_7x7, positions_likelihood=positions_likelihood + dataset=masked_imaging_7x7, positions_likelihood_list=[positions_likelihood] ) analysis_log_likelihood = analysis.log_likelihood_function(instance=instance) @@ -116,17 +116,15 @@ def test__positions__likelihood_overwrites__changes_likelihood__double_source_pl instance = model.instance_from_unit_vector([]) - plane_redshift_positions_dict = { - 1.0: al.Grid2DIrregular([(1.0, 100.0), (200.0, 2.0)]), - 2.0: al.Grid2DIrregular([(1.0, 100.0), (200.0, 2.0)]) - } - - positions_likelihood = al.PositionsLHPenalty( - plane_redshift_positions_dict=plane_redshift_positions_dict, threshold=0.01 + positions_likelihood_0 = al.PositionsLHPenalty( + plane_redshift=1.0, positions=al.Grid2DIrregular([(1.0, 100.0), (200.0, 2.0)]), threshold=0.01 + ) + positions_likelihood_1 = al.PositionsLHPenalty( + plane_redshift=2.0, positions=al.Grid2DIrregular([(1.0, 100.0), (200.0, 2.0)]), threshold=0.01 ) analysis = al.AnalysisImaging( - dataset=masked_imaging_7x7, positions_likelihood=positions_likelihood + dataset=masked_imaging_7x7, positions_likelihood_list=[positions_likelihood_0, positions_likelihood_1] ) analysis_log_likelihood = analysis.log_likelihood_function(instance=instance) From 341aba980d195669df507324a15c37fad1f3c25f Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 21 May 2025 14:52:19 +0100 Subject: [PATCH 08/11] position lieklihood list now used in analysis --- autolens/__init__.py | 3 +- autolens/analysis/analysis/dataset.py | 14 +- autolens/analysis/analysis/lens.py | 5 +- autolens/analysis/positions.py | 142 +++--------------- autolens/analysis/result.py | 119 +++++---------- autolens/imaging/model/analysis.py | 8 - autolens/imaging/model/visualizer.py | 26 +++- autolens/interferometer/model/analysis.py | 29 ++-- autolens/interferometer/model/visualizer.py | 26 +++- autolens/lens/mock/mock_tracer.py | 3 + autolens/point/max_separation.py | 2 - docs/general/demagnified_solutions.rst | 6 +- .../analysis/test_analysis_dataset.py | 4 +- test_autolens/analysis/test_analysis.py | 4 +- test_autolens/analysis/test_positions.py | 14 +- test_autolens/analysis/test_result.py | 67 +-------- .../imaging/model/test_analysis_imaging.py | 28 +--- .../model/test_analysis_interferometer.py | 26 +--- 18 files changed, 146 insertions(+), 380 deletions(-) diff --git a/autolens/__init__.py b/autolens/__init__.py index b4c26fcfc..520c102a4 100644 --- a/autolens/__init__.py +++ b/autolens/__init__.py @@ -89,8 +89,7 @@ from .lens.tracer import Tracer from .lens.sensitivity import SubhaloSensitivityResult from .lens.to_inversion import TracerToInversion -from .analysis.positions import PositionsLHResample -from .analysis.positions import PositionsLHPenalty +from .analysis.positions import PositionsLH from .imaging.simulator import SimulatorImaging from .imaging.fit_imaging import FitImaging from .imaging.model.analysis import AnalysisImaging diff --git a/autolens/analysis/analysis/dataset.py b/autolens/analysis/analysis/dataset.py index 654c9aaf4..3e7b03c8a 100644 --- a/autolens/analysis/analysis/dataset.py +++ b/autolens/analysis/analysis/dataset.py @@ -13,8 +13,7 @@ from autolens.analysis.analysis.lens import AnalysisLens from autolens.analysis.result import ResultDataset -from autolens.analysis.positions import PositionsLHResample -from autolens.analysis.positions import PositionsLHPenalty +from autolens.analysis.positions import PositionsLH from autolens import exc @@ -28,7 +27,7 @@ def __init__( self, dataset, positions_likelihood_list: Optional[ - List[Union[PositionsLHResample, PositionsLHPenalty]] + List[PositionsLH] ] = None, adapt_image_maker: Optional[ag.AdaptImageMaker] = None, cosmology: ag.cosmo.LensingCosmology = ag.cosmo.Planck15(), @@ -52,7 +51,8 @@ def __init__( positions_likelihood_list Alters the likelihood function to include a term which accounts for whether image-pixel coordinates in arc-seconds corresponding to the multiple images of each lensed source galaxy trace close to one another in - their source-plane. + their source-plane. This is a list, as it may support multiple planes, where a positions likelihood object + is input for each plane (e.g. double source plane lensing). adapt_images Contains the adapt-images which are used to make a pixelization's mesh and regularization adapt to the reconstructed galaxy's morphology. @@ -62,7 +62,7 @@ def __init__( Settings controlling how an inversion is fitted during the model-fit, for example which linear algebra formalism is used. raise_inversion_positions_likelihood_exception - If an inversion is used without the `positions_likelihood` it is likely a systematic solution will + If an inversion is used without the `positions_likelihood_list` it is likely a systematic solution will be inferred, in which case an Exception is raised before the model-fit begins to inform the user of this. This exception is not raised if this input is False, allowing the user to perform the model-fit anyway. @@ -124,7 +124,7 @@ def raise_exceptions(self, model): if has_pix: if ( - self.positions_likelihood is None + self.positions_likelihood_list is None and self.raise_inversion_positions_likelihood_exception and not conf.instance["general"]["test"][ "disable_positions_lh_inversion_check" @@ -133,7 +133,7 @@ def raise_exceptions(self, model): raise exc.AnalysisException( """ You have begun a model-fit which reconstructs the source using a pixelization. - However, you have not input a `positions_likelihood` object. + However, you have not input a `positions_likelihood_list` object. It is likely your model-fit will infer an inaccurate solution. Please read the following readthedocs page for a description of why this is, and how to set up diff --git a/autolens/analysis/analysis/lens.py b/autolens/analysis/analysis/lens.py index 78fd223bc..9675c3171 100644 --- a/autolens/analysis/analysis/lens.py +++ b/autolens/analysis/analysis/lens.py @@ -6,8 +6,7 @@ import autoarray as aa import autogalaxy as ag -from autolens.analysis.positions import PositionsLHResample -from autolens.analysis.positions import PositionsLHPenalty +from autolens.analysis.positions import PositionsLH from autolens.lens.tracer import Tracer from autolens.lens import tracer_util @@ -23,7 +22,7 @@ class AnalysisLens: def __init__( self, positions_likelihood_list: Optional[ - List[Union[PositionsLHResample, PositionsLHPenalty]] + List[PositionsLH] ] = None, cosmology: ag.cosmo.LensingCosmology = ag.cosmo.Planck15(), ): diff --git a/autolens/analysis/positions.py b/autolens/analysis/positions.py index 0b2b5b8b5..74d92a8f5 100644 --- a/autolens/analysis/positions.py +++ b/autolens/analysis/positions.py @@ -21,11 +21,12 @@ from autolens import exc -class AbstractPositionsLH: +class PositionsLH: def __init__( self, + positions: aa.Grid2DIrregular, threshold: float, - positions: Optional[aa.Grid2DIrregular] = None, + log_likelihood_penalty_factor: float = 1e8, plane_redshift: Optional[float] = None, ): """ @@ -37,9 +38,17 @@ def __init__( on the basis that this indicates an unphysical or inaccurate mass model. If they trace within the threshold the penalty term is not applied. + For the `PositionsLH` object, if the multiple image coordinates do not trace within the source-plane + threshold of one another a penalty to the likelihood is applied: + + `log_Likelihood_penalty_base - log_likelihood_penalty_factor * (max_source_plane_separation - threshold)` + + The penalty term reduces as the source-plane coordinates trace closer to one another, meaning that the + initial stages of the non-linear search can sample mass models that reduce the threshold. + For example, for one penalty term, if the multiple image coordinates are defined via `positions=aa.Grid2DIrregular([(1.0, 0.0), (-1.0, 0.0)]` and they do not trace within `threshold=0.3` of - one another, the mass model will be rejected and a new model sampled. + one another, the mass model will receive a large likelihood penalty. The default behaviour assumes a single lens plane and single source plane, meaning the input `positions` are the image-plane coordinates of one source galaxy at a specifc redshift. @@ -55,11 +64,13 @@ def __init__( threshold If the maximum separation of any two source plane coordinates is above the threshold the penalty term is applied. + log_likelihood_penalty_factor + A factor which multiplies how far source pixels do not trace within the threshold of one another, with a + larger factor producing a larger penalty making the non-linear parameter space gradient steeper. plane_redshift The plane redshift of the lensed source multiple images, which is only required if position threshold for a double source plane lens system is being used where the specific plane is required. """ - self.positions = positions self.threshold = threshold self.plane_redshift = plane_redshift @@ -71,10 +82,7 @@ def __init__( "Please input more positions into the Positions." ) - def log_likelihood_function_positions_overwrite( - self, instance: af.ModelInstance, analysis: AnalysisDataset - ) -> Optional[float]: - raise NotImplementedError + self.log_likelihood_penalty_factor = log_likelihood_penalty_factor def output_positions_info(self, output_path: str, tracer: Tracer): """ @@ -94,7 +102,7 @@ def output_positions_info(self, output_path: str, tracer: Tracer): ------- """ - with open_(path.join(output_path, "positions.info"), "w+") as f: + with open_(path.join(output_path, "positions.info"), "a+") as f: positions_fit = SourceMaxSeparation( data=self.positions, noise_map=None, tracer=tracer, plane_redshift=self.plane_redshift @@ -104,7 +112,11 @@ def output_positions_info(self, output_path: str, tracer: Tracer): coordinate=(0.0, 0.0) ) - f.write(f"Plane Redshift: {self.plane_redshift} \n") + if self.plane_redshift is None: + f.write(f"Plane Index: -1 \n") + else: + f.write(f"Plane Redshift: {self.plane_redshift} \n") + f.write(f"Positions: \n {self.positions} \n\n") f.write(f"Radial Distance from (0.0, 0.0): \n {distances} \n\n") f.write(f"Threshold = {self.threshold} \n") @@ -113,116 +125,6 @@ def output_positions_info(self, output_path: str, tracer: Tracer): ) f.write("") - -class PositionsLHResample(AbstractPositionsLH): - """ - The `PositionsLH` objects add a penalty term to the likelihood of the **PyAutoLens** `log_likelihood_function` - defined in the `Analysis` classes. - - The penalty term inspects the distance that the locations of the multiple images of the lensed source galaxy - trace within one another in the source-plane and penalizes solutions where they trace far from one another, - on the basis that this indicates an unphysical or inaccurate mass model. If they trace within the - threshold the penalty term is not applied. - - For the `PositionsLHResample` object, if the multiple image coordinates do not trace within the source-plane - threshold of one another the mass model is rejected and a new model is sampled. - - The penalty term rejects any model where the source-plane coordinates do not trace within the threshold, meaning - that the initial stages of the non-linear search may need to sample many mass models randomly in order to sample - an initial set that that trace within the threshold. - - Parameters - ---------- - positions - The arcsecond coordinates of the lensed source multiple images which are used to compute the likelihood - penalty. - threshold - If the maximum separation of any two source plane coordinates is above the threshold the penalty term - is applied. - """ - - def log_likelihood_function_positions_overwrite( - self, instance: af.ModelInstance, analysis: AnalysisDataset - ) -> Optional[float]: - """ - This is called in the `log_likelihood_function` of certain `Analysis` classes to add the penalty term of - this class, which rejects and resamples mass models which do not trace within the threshold of one another - in the source-plane. - - Parameters - ---------- - instance - The instance of the lens model that is being fitted for this iteration of the non-linear search. - analysis - The analysis class from which the log likliehood function is called. - """ - tracer = analysis.tracer_via_instance_from(instance=instance) - - if not tracer.has(cls=ag.mp.MassProfile) or len(tracer.planes) == 1: - return - - positions_fit = SourceMaxSeparation( - data=self.positions, - noise_map=None, - tracer=tracer, - plane_redshift=self.plane_redshift, - ) - - if not positions_fit.max_separation_within_threshold(self.threshold): - if os.environ.get("PYAUTOFIT_TEST_MODE") == "1": - return - - raise exc.RayTracingException - - -class PositionsLHPenalty(AbstractPositionsLH): - def __init__( - self, - threshold: float, - log_likelihood_penalty_factor: float = 1e8, - positions: Optional[aa.Grid2DIrregular] = None, - plane_redshift: Optional[float] = None, - ): - """ - The `PositionsLH` objects add a penalty term to the likelihood of the **PyAutoLens** `log_likelihood_function` - defined in the `Analysis` classes. - - The penalty term inspects the distance that the locations of the multiple images of the lensed source galaxy - trace within one another in the source-plane and penalizes solutions where they trace far from one another, - on the basis that this indicates an unphysical or inaccurate mass model. If they trace within the - threshold the penalty term is not applied. - - For the `PositionsLHPenalty` object, if the multiple image coordinates do not trace within the source-plane - threshold of one another a penalty to the likelihood is applied: - - `log_Likelihood_penalty_base - log_likelihood_penalty_factor * (max_source_plane_separation - threshold)` - - The penalty term reduces as the source-plane coordinates trace closer to one another, meaning that the - initial stages of the non-linear search can sample mass models that reduce the threshold. - - Parameters - ---------- - positions - The arcsecond coordinates of the lensed source multiple images which are used to compute the likelihood - penalty. - threshold - If the maximum separation of any two source plane coordinates is above the threshold the penalty term - is applied. - log_likelihood_penalty_factor - A factor which multiplies how far source pixels do not trace within the threshold of one another, with a - larger factor producing a larger penalty making the non-linear parameter space gradient steeper. - plane_redshift - The plane redshift of the lensed source multiple images, which is only required if position threshold - for a double source plane lens system is being used where the specific plane is required. - """ - super().__init__( - positions=positions, - threshold=threshold, - plane_redshift=plane_redshift - ) - - self.log_likelihood_penalty_factor = log_likelihood_penalty_factor - def log_likelihood_penalty_base_from( self, dataset: Union[aa.Imaging, aa.Interferometer] ) -> float: diff --git a/autolens/analysis/result.py b/autolens/analysis/result.py index cfb0b64d7..015abe0ad 100644 --- a/autolens/analysis/result.py +++ b/autolens/analysis/result.py @@ -8,8 +8,7 @@ from autogalaxy.analysis.result import ResultDataset as AgResultDataset -from autolens.analysis.positions import PositionsLHResample -from autolens.analysis.positions import PositionsLHPenalty +from autolens.analysis.positions import PositionsLH from autolens.point.max_separation import ( SourceMaxSeparation, ) @@ -27,34 +26,6 @@ def max_log_likelihood_tracer(self) -> Tracer: """ return self.analysis.tracer_via_instance_from(instance=self.instance) - @property - def max_log_likelihood_positions_threshold(self) -> float: - """ - If the `Analysis` has a `PositionsLH` object this add a penalty term to the likelihood of the - `log_likelihood_function`. - - This term is computed by ray-tracing a set of multiple image positions (E.g. corresponding to the lensed - sources brightest pixels) to the source-plane and computing their maximum separation. This separation is - then combined with a threshold to compute the likelihood term. - - This property returns the maximum separation of this `Analysis` object's multiple image positions for the - maximum log likelihood tracer. - - It therefore provides information on how closely the lens model was able to ray trace the multiple images to - one another in the source plane, and can be used for setting up the `PositionsLH` object in subsequent fits. - - Returns - ------- - - """ - positions_fits = SourceMaxSeparation( - data=self.analysis.positions_likelihood.positions, - noise_map=None, - tracer=self.max_log_likelihood_tracer, - ) - - return positions_fits.max_separation_of_plane_positions - @property def source_plane_light_profile_centre(self) -> aa.Grid2DIrregular: """ @@ -89,7 +60,15 @@ def image_plane_multiple_image_positions(self, plane_redshift : Optional[float] the multiple image position of the source(s) in the image-plane. These image-plane positions are used by the next search in a pipeline if automatic position updating is turned - on.""" + on. + + Parameters + ---------- + plane_redshift + The plane redshift of the lensed source, which is only required when a double (or triple) source plane lens + system is being analysed where the specific plane via its redshift is required to define whihch source + galaxy is used to compute the multiple images. + """ grid = self.analysis.dataset.mask.derive_grid.all_false @@ -131,6 +110,10 @@ def image_plane_multiple_image_positions_for_single_image_from( Parameters ---------- + plane_redshift + The plane redshift of the lensed source, which is only required when a double (or triple) source plane lens + system is being analysed where the specific plane via its redshift is required to define whihch source + galaxy is used to compute the multiple images. increments The number of increments the source-plane centre is moved to compute multiple images. """ @@ -209,6 +192,10 @@ def positions_threshold_from( positions If input, these positions are used instead of the computed multiple image positions from the lens mass model. + plane_redshift + The plane redshift of the lensed source, which is only required when a double (or triple) source plane lens + system is being analysed where the specific plane via its redshift is required to define whihch source + galaxy is used to compute the multiple images. Returns ------- @@ -216,7 +203,6 @@ def positions_threshold_from( The maximum source plane separation of this results maximum likelihood `Tracer` multiple images multiplied by `factor` and rounded up to the `threshold`. """ - plane_redshift_positions_dict = {} tracer = Tracer(galaxies=self.max_log_likelihood_galaxies) @@ -244,11 +230,10 @@ def positions_likelihood_from( self, factor=1.0, minimum_threshold=None, - use_resample=False, positions: Optional[aa.Grid2DIrregular] = None, mass_centre_radial_distance_min: float = None, - plane_redshift_list: Optional[List[int]] = None, - ) -> Union[PositionsLHPenalty, PositionsLHResample]: + plane_redshift: Optional[float] = None, + ) -> PositionsLH: """ 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 @@ -271,10 +256,6 @@ def positions_likelihood_from( 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. @@ -282,6 +263,10 @@ def positions_likelihood_from( 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. + plane_redshift + The plane redshift of the lensed source, which is only required when a double (or triple) source plane lens + system is being analysed where the specific plane via its redshift is required to define whihch source + galaxy is used to compute the multiple images. Returns ------- @@ -291,43 +276,28 @@ def positions_likelihood_from( if os.environ.get("PYAUTOFIT_TEST_MODE") == "1": return - plane_redshift_list = [None] if plane_redshift_list is None else plane_redshift_list - - plane_redshift_positions_dict = {} - threshold_list = [] - - for plane_redshift in plane_redshift_list: + positions = ( + self.image_plane_multiple_image_positions(plane_redshift=plane_redshift) + if positions is None + else positions + ) - positions = ( - self.image_plane_multiple_image_positions(plane_redshift=plane_redshift) - 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" ) - 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, plane_redshift=plane_redshift + distances = positions.distances_to_coordinate_from( + coordinate=mass_centre[0] ) - threshold_list.append(threshold) - - plane_redshift_positions_dict[plane_redshift] = positions - - threshold = np.max(threshold_list) + positions = positions[distances > mass_centre_radial_distance_min] - if not use_resample: - return PositionsLHPenalty(plane_redshift_positions_dict=plane_redshift_positions_dict, threshold=threshold) - return PositionsLHResample(plane_redshift_positions_dict=plane_redshift_positions_dict, threshold=threshold) + threshold = self.positions_threshold_from( + factor=factor, minimum_threshold=minimum_threshold, positions=positions, plane_redshift=plane_redshift + ) + + return PositionsLH(positions=positions, threshold=threshold, plane_redshift=plane_redshift) class ResultDataset(Result): @@ -344,15 +314,6 @@ def max_log_likelihood_tracer(self) -> Tracer: return self.analysis.tracer_via_instance_from(instance=instance) - @property - def positions(self): - """ - The (y,x) arc-second coordinates of the lensed sources brightest pixels, which are used for discarding mass - models which do not trace within a threshold in the source-plane of one another. - """ - if self.analysis.positions_likelihood is not None: - return self.analysis.positions_likelihood.positions - @property def source_plane_centre(self) -> aa.Grid2DIrregular: """ diff --git a/autolens/imaging/model/analysis.py b/autolens/imaging/model/analysis.py index 24a2ff8a1..1c2c727ad 100644 --- a/autolens/imaging/model/analysis.py +++ b/autolens/imaging/model/analysis.py @@ -205,14 +205,6 @@ def save_attributes(self, paths: af.DirectoryPaths): analysis.save_attributes(paths=paths) - if self.positions_likelihood is not None: - - paths.save_json( - name="positions", - object_dict=to_dict(self.positions_likelihood.positions), - prefix="dataset", - ) - def profile_log_likelihood_function( self, instance: af.ModelInstance, paths: Optional[af.DirectoryPaths] = None ) -> Tuple[Dict, Dict]: diff --git a/autolens/imaging/model/visualizer.py b/autolens/imaging/model/visualizer.py index b758ea5ff..de1a70993 100644 --- a/autolens/imaging/model/visualizer.py +++ b/autolens/imaging/model/visualizer.py @@ -35,10 +35,21 @@ def visualize_before_fit( plotter_interface.imaging(dataset=analysis.dataset) - if analysis.positions_likelihood is not None: + if analysis.positions_likelihood_list is not None: + + positions_list = [] + + for positions_likelihood in analysis.positions_likelihood_list: + + positions_list.append( + positions_likelihood.positions + ) + + positions = ag.Grid2DIrregular(positions_list) + plotter_interface.image_with_positions( image=analysis.dataset.data, - positions=analysis.positions_likelihood.positions, + positions=positions, ) if analysis.adapt_images is not None: @@ -80,10 +91,13 @@ def visualize( fit = analysis.fit_from(instance=instance) - if analysis.positions_likelihood is not None: - analysis.positions_likelihood.output_positions_info( - output_path=paths.output_path, tracer=fit.tracer - ) + if analysis.positions_likelihood_list is not None: + + for positions_likelihood in analysis.positions_likelihood_list: + + positions_likelihood.output_positions_info( + output_path=paths.output_path, tracer=fit.tracer + ) if fit.inversion is not None: try: diff --git a/autolens/interferometer/model/analysis.py b/autolens/interferometer/model/analysis.py index 3dec88777..a8cdeb214 100644 --- a/autolens/interferometer/model/analysis.py +++ b/autolens/interferometer/model/analysis.py @@ -2,8 +2,6 @@ import numpy as np from typing import Dict, Optional, Tuple, Union -from autoconf.dictable import to_dict - import autofit as af import autoarray as aa import autogalaxy as ag @@ -11,8 +9,7 @@ from autoarray.exc import PixelizationException from autolens.analysis.analysis.dataset import AnalysisDataset -from autolens.analysis.positions import PositionsLHResample -from autolens.analysis.positions import PositionsLHPenalty +from autolens.analysis.positions import PositionsLH from autolens.interferometer.model.result import ResultInterferometer from autolens.interferometer.model.visualizer import VisualizerInterferometer from autolens.interferometer.fit_interferometer import FitInterferometer @@ -31,8 +28,8 @@ class AnalysisInterferometer(AnalysisDataset): def __init__( self, dataset, - positions_likelihood: Optional[ - Union[PositionsLHResample, PositionsLHPenalty] + positions_likelihood_list: Optional[ + PositionsLH ] = None, adapt_image_maker: Optional[ag.AdaptImageMaker] = None, cosmology: ag.cosmo.LensingCosmology = ag.cosmo.Planck15(), @@ -60,10 +57,11 @@ def __init__( ---------- dataset The interferometer dataset that the model is fitted too. - positions_likelihood - An object which alters the likelihood function to include a term which accounts for whether - image-pixel coordinates in arc-seconds corresponding to the multiple images of the lensed source galaxy - trace close to one another in the source-plane. + positions_likelihood_list + Alters the likelihood function to include a term which accounts for whether image-pixel coordinates in + arc-seconds corresponding to the multiple images of each lensed source galaxy trace close to one another in + their source-plane. This is a list, as it may support multiple planes, where a positions likelihood object + is input for each plane (e.g. double source plane lensing). adapt_images Contains the adapt-images which are used to make a pixelization's mesh and regularization adapt to the reconstructed galaxy's morphology. @@ -72,7 +70,7 @@ def __init__( settings_inversion Settings controlling how an inversion is fitted, for example which linear algebra formalism is used. raise_inversion_positions_likelihood_exception - If an inversion is used without the `positions_likelihood` it is likely a systematic solution will + If an inversion is used without the `positions_likelihood_list` it is likely a systematic solution will be inferred, in which case an Exception is raised before the model-fit begins to inform the user of this. This exception is not raised if this input is False, allowing the user to perform the model-fit anyway. @@ -82,7 +80,7 @@ def __init__( """ super().__init__( dataset=dataset, - positions_likelihood=positions_likelihood, + positions_likelihood_list=positions_likelihood_list, adapt_image_maker=adapt_image_maker, cosmology=cosmology, settings_inversion=settings_inversion, @@ -261,13 +259,6 @@ def save_attributes(self, paths: af.DirectoryPaths): analysis.save_attributes(paths=paths) - if self.positions_likelihood is not None: - paths.save_json( - name="positions", - object_dict=to_dict(self.positions_likelihood.positions), - prefix="dataset", - ) - def profile_log_likelihood_function( self, instance: af.ModelInstance, paths: Optional[af.DirectoryPaths] = None ) -> Tuple[Dict, Dict]: diff --git a/autolens/interferometer/model/visualizer.py b/autolens/interferometer/model/visualizer.py index 74cb75eb7..2eade0b88 100644 --- a/autolens/interferometer/model/visualizer.py +++ b/autolens/interferometer/model/visualizer.py @@ -1,4 +1,5 @@ import autofit as af +import autogalaxy as ag from autolens.interferometer.model.plotter_interface import ( PlotterInterfaceInterferometer, @@ -34,10 +35,20 @@ def visualize_before_fit( plotter_interface.interferometer(dataset=analysis.interferometer) - if analysis.positions_likelihood is not None: + if analysis.positions_likelihood_list is not None: + + positions_list = [] + + for positions_likelihood in analysis.positions_likelihood_list: + positions_list.append( + positions_likelihood.positions + ) + + positions = ag.Grid2DIrregular(positions_list) + plotter_interface.image_with_positions( image=analysis.dataset.dirty_image, - positions=analysis.positions_likelihood.positions, + positions=positions ) if analysis.adapt_images is not None: @@ -81,10 +92,13 @@ def visualize( """ fit = analysis.fit_from(instance=instance) - if analysis.positions_likelihood is not None: - analysis.positions_likelihood.output_positions_info( - output_path=paths.output_path, tracer=fit.tracer - ) + if analysis.positions_likelihood_list is not None: + + for positions_likelihood in analysis.positions_likelihood_list: + + positions_likelihood.output_positions_info( + output_path=paths.output_path, tracer=fit.tracer + ) if fit.inversion is not None: try: diff --git a/autolens/lens/mock/mock_tracer.py b/autolens/lens/mock/mock_tracer.py index 70004a9e0..b4808cf19 100644 --- a/autolens/lens/mock/mock_tracer.py +++ b/autolens/lens/mock/mock_tracer.py @@ -8,6 +8,9 @@ def __init__( def traced_grid_2d_list_from(self, grid): return self._traced_grid_2d_list_from + def plane_index_via_redshift_from(self, redshift): + raise TypeError + class MockTracerPoint(MockTracer): def __init__( diff --git a/autolens/point/max_separation.py b/autolens/point/max_separation.py index 20c8f6ae7..7a88b4075 100644 --- a/autolens/point/max_separation.py +++ b/autolens/point/max_separation.py @@ -37,8 +37,6 @@ def __init__( self.data = data self.noise_map = noise_map - print(plane_redshift) - try: plane_index = tracer.plane_index_via_redshift_from(redshift=plane_redshift) except TypeError: diff --git a/docs/general/demagnified_solutions.rst b/docs/general/demagnified_solutions.rst index 6a2309d27..51297e87b 100644 --- a/docs/general/demagnified_solutions.rst +++ b/docs/general/demagnified_solutions.rst @@ -70,7 +70,7 @@ stars: The ``autolens_workspace`` also includes a Graphical User Interface for drawing lensed source positions via mouse click (https://github.com/Jammy2211/autolens_workspace/blob/release/scripts/imaging/preprocess/gui/positions.py). -Next, we create ``PositionsLHPenalty`` object, which has an input ``threshold``. +Next, we create ``PositionsLH`` object, which has an input ``threshold``. This requires that a mass model traces the multiple image ``positions`` specified above within the ``threshold`` value (e.g. 0.5") of one another in the source-plane. If this criteria is not met, a large penalty term is @@ -88,7 +88,7 @@ The penalty term is created and passed to an ``Analysis`` object as follows: .. code-block:: python - positions_likelihood = al.PositionsLHPenalty(positions=positions, threshold=0.3) + positions_likelihood = al.PositionsLH(positions=positions, threshold=0.3) analysis = al.AnalysisImaging( dataset=dataset, positions_likelihood=positions_likelihood @@ -155,7 +155,7 @@ These inputs are useful when using function to set the ``threshold`` in a new fi it allows us to make sure we do set too small a threshold that we remove genuinely physically mass models. For writing search chaining pipelines, a convenience method is available in the ``result`` which returns directly a -``PositionsLHPenalty`` object: +``PositionsLH`` object: .. code-block:: python diff --git a/test_autolens/analysis/analysis/test_analysis_dataset.py b/test_autolens/analysis/analysis/test_analysis_dataset.py index 4f93db1c4..df21c174c 100644 --- a/test_autolens/analysis/analysis/test_analysis_dataset.py +++ b/test_autolens/analysis/analysis/test_analysis_dataset.py @@ -30,12 +30,12 @@ def test__modify_before_fit__inversion_no_positions_likelihood__raises_exception with pytest.raises(exc.AnalysisException): analysis.modify_before_fit(paths=af.DirectoryPaths(), model=model) - positions_likelihood = al.PositionsLHPenalty( + positions_likelihood = al.PositionsLH( positions=al.Grid2DIrregular([(1.0, 100.0), (200.0, 2.0)]), threshold=0.01 ) analysis = al.AnalysisImaging( - dataset=masked_imaging_7x7, positions_likelihood=positions_likelihood + dataset=masked_imaging_7x7, positions_likelihood_list=[positions_likelihood] ) analysis.modify_before_fit(paths=af.DirectoryPaths(), model=model) diff --git a/test_autolens/analysis/test_analysis.py b/test_autolens/analysis/test_analysis.py index 93d4f1b2c..38fb73bd0 100644 --- a/test_autolens/analysis/test_analysis.py +++ b/test_autolens/analysis/test_analysis.py @@ -140,12 +140,12 @@ def test__modify_before_fit__inversion_no_positions_likelihood__raises_exception with pytest.raises(exc.AnalysisException): analysis.modify_before_fit(paths=af.DirectoryPaths(), model=model) - positions_likelihood = al.PositionsLHPenalty( + positions_likelihood = al.PositionsLH( positions=al.Grid2DIrregular([(1.0, 100.0), (200.0, 2.0)]), threshold=0.01 ) analysis = al.AnalysisImaging( - dataset=masked_imaging_7x7, positions_likelihood=positions_likelihood + dataset=masked_imaging_7x7, positions_likelihood_list=[positions_likelihood] ) analysis.modify_before_fit(paths=af.DirectoryPaths(), model=model) diff --git a/test_autolens/analysis/test_positions.py b/test_autolens/analysis/test_positions.py index fa393e4d5..176f4696b 100644 --- a/test_autolens/analysis/test_positions.py +++ b/test_autolens/analysis/test_positions.py @@ -10,14 +10,14 @@ def test__check_positions_on_instantiation(): - al.PositionsLHResample( + al.PositionsLH( positions=al.Grid2DIrregular([(1.0, 2.0), (3.0, 4.0)]), threshold=0.1 ) # Positions input with threshold but positions are length 1. with pytest.raises(exc.PositionsException): - al.PositionsLHResample( + al.PositionsLH( positions=al.Grid2DIrregular([(1.0, 2.0)]), threshold=0.1 ) @@ -27,7 +27,7 @@ def test__output_positions_info(): "{}".format(os.path.dirname(os.path.realpath(__file__))), "files" ) - positions_likelihood = al.PositionsLHResample( + positions_likelihood = al.PositionsLH( positions=al.Grid2DIrregular([(1.0, 2.0), (3.0, 4.0)]), threshold=0.1 ) @@ -52,7 +52,7 @@ def test__output_positions_info(): def make_settings_dict(): return { "type": "instance", - "class_path": "autolens.analysis.positions.PositionsLHPenalty", + "class_path": "autolens.analysis.positions.PositionsLH", "arguments": { "positions": { "type": "ndarray", @@ -66,20 +66,20 @@ def make_settings_dict(): def test_settings_from_dict(settings_dict): - assert isinstance(from_dict(settings_dict), al.PositionsLHPenalty) + assert isinstance(from_dict(settings_dict), al.PositionsLH) def test_file(): filename = "/tmp/temp.json" output_to_json( - al.PositionsLHPenalty( + al.PositionsLH( positions=al.Grid2DIrregular([(1.0, 2.0), (3.0, 4.0)]), threshold=0.1 ), filename, ) try: - assert isinstance(from_json(filename), al.PositionsLHPenalty) + assert isinstance(from_json(filename), al.PositionsLH) finally: os.remove(filename) diff --git a/test_autolens/analysis/test_result.py b/test_autolens/analysis/test_result.py index 32416c150..1c082058c 100644 --- a/test_autolens/analysis/test_result.py +++ b/test_autolens/analysis/test_result.py @@ -31,36 +31,6 @@ def test__max_log_likelihood_tracer( assert isinstance(result.max_log_likelihood_tracer.galaxies[0], al.Galaxy) -def test__max_log_likelihood_positions_threshold(masked_imaging_7x7): - positions_likelihood = al.PositionsLHResample( - positions=al.Grid2DIrregular(values=[(1.0, 1.0), [-1.0, -1.0]]), threshold=100.0 - ) - - analysis = al.AnalysisImaging( - dataset=masked_imaging_7x7, positions_likelihood=positions_likelihood - ) - - 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) - - assert result.max_log_likelihood_positions_threshold == pytest.approx( - 0.8309561230, 1.0e-4 - ) - - def test__source_plane_light_profile_centre(analysis_imaging_7x7): lens = al.Galaxy(redshift=0.5, light=al.lp.SersicSph(intensity=1.0)) @@ -280,17 +250,9 @@ def test__positions_likelihood_from(analysis_imaging_7x7): factor=0.1, minimum_threshold=0.2 ) - assert isinstance(positions_likelihood, al.PositionsLHPenalty) + assert isinstance(positions_likelihood, al.PositionsLH) assert positions_likelihood.threshold == pytest.approx(0.2, 1.0e-4) - positions_likelihood = result.positions_likelihood_from( - factor=0.1, minimum_threshold=0.2, use_resample=True - ) - - assert isinstance(positions_likelihood, al.PositionsLHResample) - assert positions_likelihood.threshold == pytest.approx(0.2, 1.0e-4) - - def test__positions_likelihood_from__mass_centre_radial_distance_min( analysis_imaging_7x7, ): @@ -314,7 +276,7 @@ def test__positions_likelihood_from__mass_centre_radial_distance_min( factor=0.1, minimum_threshold=0.2, mass_centre_radial_distance_min=0.1 ) - assert isinstance(positions_likelihood, al.PositionsLHPenalty) + assert isinstance(positions_likelihood, al.PositionsLH) assert len(positions_likelihood.positions) == 2 assert positions_likelihood.positions[0] == pytest.approx( (-1.00097656e00, 5.63818622e-04), 1.0e-4 @@ -335,31 +297,6 @@ def test__results_include_mask__available_as_property( assert (result.mask == masked_imaging_7x7.mask).all() -def test__results_include_positions__available_as_property( - analysis_imaging_7x7, masked_imaging_7x7, samples_summary_with_result -): - result = res.ResultDataset( - samples_summary=samples_summary_with_result, - analysis=analysis_imaging_7x7, - ) - - assert result.positions == None - - positions_likelihood = al.PositionsLHResample( - positions=al.Grid2DIrregular([(1.0, 100.0), (200.0, 2.0)]), threshold=1.0 - ) - - analysis = al.AnalysisImaging( - dataset=masked_imaging_7x7, positions_likelihood=positions_likelihood - ) - - result = res.ResultDataset( - samples_summary=samples_summary_with_result, - analysis=analysis, - ) - - assert (result.positions[0] == np.array([1.0, 100.0])).all() - def test___image_dict(analysis_imaging_7x7): galaxies = af.ModelInstance() diff --git a/test_autolens/imaging/model/test_analysis_imaging.py b/test_autolens/imaging/model/test_analysis_imaging.py index 108d4ddea..7453d6f59 100644 --- a/test_autolens/imaging/model/test_analysis_imaging.py +++ b/test_autolens/imaging/model/test_analysis_imaging.py @@ -44,28 +44,6 @@ def test__figure_of_merit__matches_correct_fit_given_galaxy_profiles( assert fit.log_likelihood == analysis_log_likelihood -def test__positions__resample__raises_exception(masked_imaging_7x7): - model = af.Collection( - galaxies=af.Collection( - lens=al.Galaxy(redshift=0.5, mass=al.mp.IsothermalSph()), - source=al.Galaxy(redshift=1.0), - ) - ) - - positions_likelihood = al.PositionsLHResample( - positions=al.Grid2DIrregular([(1.0, 100.0), (200.0, 2.0)]), threshold=0.01 - ) - - analysis = al.AnalysisImaging( - dataset=masked_imaging_7x7, positions_likelihood_list=[positions_likelihood] - ) - - instance = model.instance_from_unit_vector([]) - - with pytest.raises(exc.RayTracingException): - analysis.log_likelihood_function(instance=instance) - - def test__positions__likelihood_overwrites__changes_likelihood(masked_imaging_7x7): lens = al.Galaxy(redshift=0.5, mass=al.mp.IsothermalSph(centre=(0.05, 0.05))) source = al.Galaxy(redshift=1.0, light=al.lp.SersicSph(centre=(0.05, 0.05))) @@ -84,7 +62,7 @@ def test__positions__likelihood_overwrites__changes_likelihood(masked_imaging_7x assert fit.log_likelihood == pytest.approx(analysis_log_likelihood, 1.0e-4) assert analysis_log_likelihood == pytest.approx(-14.79034680979, 1.0e-4) - positions_likelihood = al.PositionsLHPenalty( + positions_likelihood = al.PositionsLH( positions=al.Grid2DIrregular([(1.0, 100.0), (200.0, 2.0)]), threshold=0.01 ) @@ -116,10 +94,10 @@ def test__positions__likelihood_overwrites__changes_likelihood__double_source_pl instance = model.instance_from_unit_vector([]) - positions_likelihood_0 = al.PositionsLHPenalty( + positions_likelihood_0 = al.PositionsLH( plane_redshift=1.0, positions=al.Grid2DIrregular([(1.0, 100.0), (200.0, 2.0)]), threshold=0.01 ) - positions_likelihood_1 = al.PositionsLHPenalty( + positions_likelihood_1 = al.PositionsLH( plane_redshift=2.0, positions=al.Grid2DIrregular([(1.0, 100.0), (200.0, 2.0)]), threshold=0.01 ) diff --git a/test_autolens/interferometer/model/test_analysis_interferometer.py b/test_autolens/interferometer/model/test_analysis_interferometer.py index 3f62c28b1..017d3d28b 100644 --- a/test_autolens/interferometer/model/test_analysis_interferometer.py +++ b/test_autolens/interferometer/model/test_analysis_interferometer.py @@ -39,28 +39,6 @@ def test__figure_of_merit__matches_correct_fit_given_galaxy_profiles(interferome assert fit.log_likelihood == analysis_log_likelihood -def test__positions__resample__raises_exception(interferometer_7, mask_2d_7x7): - model = af.Collection( - galaxies=af.Collection( - lens=al.Galaxy(redshift=0.5, mass=al.mp.IsothermalSph()), - source=al.Galaxy(redshift=1.0), - ) - ) - - positions_likelihood = al.PositionsLHResample( - positions=al.Grid2DIrregular([(1.0, 100.0), (200.0, 2.0)]), threshold=0.01 - ) - - analysis = al.AnalysisInterferometer( - dataset=interferometer_7, positions_likelihood=positions_likelihood - ) - - instance = model.instance_from_unit_vector([]) - - with pytest.raises(exc.RayTracingException): - analysis.log_likelihood_function(instance=instance) - - def test__positions__likelihood_overwrite__changes_likelihood( interferometer_7, mask_2d_7x7 ): @@ -81,12 +59,12 @@ def test__positions__likelihood_overwrite__changes_likelihood( assert fit.log_likelihood == analysis_log_likelihood assert analysis_log_likelihood == pytest.approx(-127914.36273, 1.0e-4) - positions_likelihood = al.PositionsLHPenalty( + positions_likelihood = al.PositionsLH( positions=al.Grid2DIrregular([(1.0, 100.0), (200.0, 2.0)]), threshold=0.01 ) analysis = al.AnalysisInterferometer( - dataset=interferometer_7, positions_likelihood=positions_likelihood + dataset=interferometer_7, positions_likelihood_list=[positions_likelihood] ) analysis_log_likelihood = analysis.log_likelihood_function(instance=instance) From c2f0cfd23ab009223cf32c7b3beadba2f96dd84a Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 21 May 2025 14:52:53 +0100 Subject: [PATCH 09/11] black --- autolens/analysis/analysis/dataset.py | 4 +-- autolens/analysis/analysis/lens.py | 4 +-- autolens/analysis/positions.py | 8 +++--- autolens/analysis/result.py | 29 ++++++++++++--------- autolens/interferometer/model/analysis.py | 4 +-- autolens/interferometer/model/visualizer.py | 7 ++--- autolens/lens/tracer.py | 2 +- autolens/point/max_separation.py | 4 +-- test_autolens/analysis/test_positions.py | 4 +-- test_autolens/analysis/test_result.py | 2 +- test_autolens/lens/test_tracer.py | 1 + 11 files changed, 32 insertions(+), 37 deletions(-) diff --git a/autolens/analysis/analysis/dataset.py b/autolens/analysis/analysis/dataset.py index 3e7b03c8a..c4e12d949 100644 --- a/autolens/analysis/analysis/dataset.py +++ b/autolens/analysis/analysis/dataset.py @@ -26,9 +26,7 @@ class AnalysisDataset(AgAnalysisDataset, AnalysisLens): def __init__( self, dataset, - positions_likelihood_list: Optional[ - List[PositionsLH] - ] = None, + positions_likelihood_list: Optional[List[PositionsLH]] = None, adapt_image_maker: Optional[ag.AdaptImageMaker] = None, cosmology: ag.cosmo.LensingCosmology = ag.cosmo.Planck15(), settings_inversion: aa.SettingsInversion = None, diff --git a/autolens/analysis/analysis/lens.py b/autolens/analysis/analysis/lens.py index 9675c3171..1641756b5 100644 --- a/autolens/analysis/analysis/lens.py +++ b/autolens/analysis/analysis/lens.py @@ -21,9 +21,7 @@ class AnalysisLens: def __init__( self, - positions_likelihood_list: Optional[ - List[PositionsLH] - ] = None, + positions_likelihood_list: Optional[List[PositionsLH]] = None, cosmology: ag.cosmo.LensingCosmology = ag.cosmo.Planck15(), ): """ diff --git a/autolens/analysis/positions.py b/autolens/analysis/positions.py index 74d92a8f5..c57899319 100644 --- a/autolens/analysis/positions.py +++ b/autolens/analysis/positions.py @@ -105,7 +105,10 @@ def output_positions_info(self, output_path: str, tracer: Tracer): with open_(path.join(output_path, "positions.info"), "a+") as f: positions_fit = SourceMaxSeparation( - data=self.positions, noise_map=None, tracer=tracer, plane_redshift=self.plane_redshift + data=self.positions, + noise_map=None, + tracer=tracer, + plane_redshift=self.plane_redshift, ) distances = positions_fit.data.distances_to_coordinate_from( @@ -209,8 +212,7 @@ def log_likelihood_penalty_from(self, tracer: Tracer) -> Optional[float]: if not positions_fit.max_separation_within_threshold(self.threshold): log_likelihood_penalty += self.log_likelihood_penalty_factor * ( - positions_fit.max_separation_of_plane_positions - - self.threshold + positions_fit.max_separation_of_plane_positions - self.threshold ) return log_likelihood_penalty diff --git a/autolens/analysis/result.py b/autolens/analysis/result.py index 015abe0ad..ab9a9875b 100644 --- a/autolens/analysis/result.py +++ b/autolens/analysis/result.py @@ -54,7 +54,9 @@ def source_plane_centre(self) -> aa.Grid2DIrregular: """ return self.source_plane_light_profile_centre - def image_plane_multiple_image_positions(self, plane_redshift : Optional[float] = None) -> aa.Grid2DIrregular: + def image_plane_multiple_image_positions( + self, plane_redshift: Optional[float] = None + ) -> aa.Grid2DIrregular: """ Backwards ray-trace the source-plane centres (see above) to the image-plane via the mass model, to determine the multiple image position of the source(s) in the image-plane. @@ -89,7 +91,7 @@ def image_plane_multiple_image_positions(self, plane_redshift : Optional[float] return aa.Grid2DIrregular(values=multiple_images) def image_plane_multiple_image_positions_for_single_image_from( - self, plane_redshift : Optional[float] = None, increments: int = 20 + self, plane_redshift: Optional[float] = None, increments: int = 20 ) -> aa.Grid2DIrregular: """ If the standard point solver only locates one multiple image, finds one or more additional images, which are @@ -149,7 +151,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. @@ -163,7 +165,7 @@ def positions_threshold_from( factor=1.0, minimum_threshold=None, positions: Optional[aa.Grid2DIrregular] = None, - plane_redshift : Optional[float] = None, + plane_redshift: Optional[float] = None, ) -> float: """ Compute a new position threshold from these results corresponding to the image-plane multiple image positions @@ -176,10 +178,10 @@ def positions_threshold_from( This is used for non-linear search chaining, specifically updating the position threshold of a new model-fit using the maximum likelihood model of a previous search. - + The above behaviour assumes a single lens and single source plane, if there are multiple source planes (e.g. a double source plane lens) the input `plane_redshift_list` can be used to specify which source planea - the position threshold is computed for. + the position threshold is computed for. Parameters ---------- @@ -216,9 +218,7 @@ def positions_threshold_from( data=positions, noise_map=None, tracer=tracer, plane_redshift=plane_redshift ) - threshold = factor * np.max( - positions_fits.max_separation_of_plane_positions - ) + threshold = factor * np.max(positions_fits.max_separation_of_plane_positions) if minimum_threshold is not None: if threshold < minimum_threshold: @@ -294,10 +294,15 @@ def positions_likelihood_from( positions = positions[distances > mass_centre_radial_distance_min] threshold = self.positions_threshold_from( - factor=factor, minimum_threshold=minimum_threshold, positions=positions, plane_redshift=plane_redshift + factor=factor, + minimum_threshold=minimum_threshold, + positions=positions, + plane_redshift=plane_redshift, + ) + + return PositionsLH( + positions=positions, threshold=threshold, plane_redshift=plane_redshift ) - - return PositionsLH(positions=positions, threshold=threshold, plane_redshift=plane_redshift) class ResultDataset(Result): diff --git a/autolens/interferometer/model/analysis.py b/autolens/interferometer/model/analysis.py index a8cdeb214..3fbced456 100644 --- a/autolens/interferometer/model/analysis.py +++ b/autolens/interferometer/model/analysis.py @@ -28,9 +28,7 @@ class AnalysisInterferometer(AnalysisDataset): def __init__( self, dataset, - positions_likelihood_list: Optional[ - PositionsLH - ] = None, + positions_likelihood_list: Optional[PositionsLH] = None, adapt_image_maker: Optional[ag.AdaptImageMaker] = None, cosmology: ag.cosmo.LensingCosmology = ag.cosmo.Planck15(), settings_inversion: aa.SettingsInversion = None, diff --git a/autolens/interferometer/model/visualizer.py b/autolens/interferometer/model/visualizer.py index 2eade0b88..fc5415bff 100644 --- a/autolens/interferometer/model/visualizer.py +++ b/autolens/interferometer/model/visualizer.py @@ -40,15 +40,12 @@ def visualize_before_fit( positions_list = [] for positions_likelihood in analysis.positions_likelihood_list: - positions_list.append( - positions_likelihood.positions - ) + positions_list.append(positions_likelihood.positions) positions = ag.Grid2DIrregular(positions_list) plotter_interface.image_with_positions( - image=analysis.dataset.dirty_image, - positions=positions + image=analysis.dataset.dirty_image, positions=positions ) if analysis.adapt_images is not None: diff --git a/autolens/lens/tracer.py b/autolens/lens/tracer.py index 07e480c08..037a05560 100644 --- a/autolens/lens/tracer.py +++ b/autolens/lens/tracer.py @@ -801,7 +801,7 @@ def cls_list_from(self, cls: Type) -> List: return cls_list - def plane_index_via_redshift_from(self, redshift : float) -> Optional[int]: + def plane_index_via_redshift_from(self, redshift: float) -> Optional[int]: """ Returns the index of a plane at a given redshift. diff --git a/autolens/point/max_separation.py b/autolens/point/max_separation.py index 7a88b4075..d7c9d50c8 100644 --- a/autolens/point/max_separation.py +++ b/autolens/point/max_separation.py @@ -42,9 +42,7 @@ def __init__( except TypeError: plane_index = -1 - self.plane_positions = tracer.traced_grid_2d_list_from(grid=data)[ - plane_index - ] + self.plane_positions = tracer.traced_grid_2d_list_from(grid=data)[plane_index] @property def furthest_separations_of_plane_positions(self) -> aa.ArrayIrregular: diff --git a/test_autolens/analysis/test_positions.py b/test_autolens/analysis/test_positions.py index 176f4696b..33175bdf1 100644 --- a/test_autolens/analysis/test_positions.py +++ b/test_autolens/analysis/test_positions.py @@ -17,9 +17,7 @@ def test__check_positions_on_instantiation(): # Positions input with threshold but positions are length 1. with pytest.raises(exc.PositionsException): - al.PositionsLH( - positions=al.Grid2DIrregular([(1.0, 2.0)]), threshold=0.1 - ) + al.PositionsLH(positions=al.Grid2DIrregular([(1.0, 2.0)]), threshold=0.1) def test__output_positions_info(): diff --git a/test_autolens/analysis/test_result.py b/test_autolens/analysis/test_result.py index 1c082058c..143961604 100644 --- a/test_autolens/analysis/test_result.py +++ b/test_autolens/analysis/test_result.py @@ -253,6 +253,7 @@ def test__positions_likelihood_from(analysis_imaging_7x7): assert isinstance(positions_likelihood, al.PositionsLH) assert positions_likelihood.threshold == pytest.approx(0.2, 1.0e-4) + def test__positions_likelihood_from__mass_centre_radial_distance_min( analysis_imaging_7x7, ): @@ -297,7 +298,6 @@ def test__results_include_mask__available_as_property( assert (result.mask == masked_imaging_7x7.mask).all() - def test___image_dict(analysis_imaging_7x7): galaxies = af.ModelInstance() galaxies.lens = al.Galaxy(redshift=0.5) diff --git a/test_autolens/lens/test_tracer.py b/test_autolens/lens/test_tracer.py index fc3866d54..c4f7d8bfc 100644 --- a/test_autolens/lens/test_tracer.py +++ b/test_autolens/lens/test_tracer.py @@ -90,6 +90,7 @@ def test__plane_index_via_redshift_from(): assert tracer.plane_index_via_redshift_from(redshift=3.0) == 2 assert tracer.plane_index_via_redshift_from(redshift=3.001) == None + def test__upper_plane_index_with_light_profile(): g0 = al.Galaxy(redshift=0.5) g1 = al.Galaxy(redshift=1.0) From b3eb06f6b0babb83d33d1498a9652a1773b251fb Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 21 May 2025 14:58:26 +0100 Subject: [PATCH 10/11] docs --- docs/general/demagnified_solutions.rst | 45 +++++++++++++++----------- 1 file changed, 27 insertions(+), 18 deletions(-) diff --git a/docs/general/demagnified_solutions.rst b/docs/general/demagnified_solutions.rst index 51297e87b..9c4816b94 100644 --- a/docs/general/demagnified_solutions.rst +++ b/docs/general/demagnified_solutions.rst @@ -91,30 +91,13 @@ The penalty term is created and passed to an ``Analysis`` object as follows: positions_likelihood = al.PositionsLH(positions=positions, threshold=0.3) analysis = al.AnalysisImaging( - dataset=dataset, positions_likelihood=positions_likelihood + dataset=dataset, positions_likelihood_list=[positions_likelihood] ) The threshold of 0.5" is large. For an accurate lens model we would anticipate the positions trace within < 0.01" of one another. However, we only want the threshold to aid the non-linear with the choice of mass model in the initial fit and remove demagnified solutions. -Resampling ----------- - -An alternative penalty term is available via the ``PositionsLHResample`` object, which rejects and resamples a lens -model if the ``positions``do not trace within the ``threshold`` of one another in the source plane. - -This is not the recommended option, as it is slower and can often lead to prolonged periods of the non-linear search -guessing and rejecting mass models. - -.. code-block:: python - - positions_likelihood = al.PositionsLHResample(positions=positions, threshold=0.3) - - analysis = al.AnalysisImaging( - dataset=dataset, positions_likelihood=positions_likelihood - ) - Auto Position Updates --------------------- @@ -172,4 +155,30 @@ This is often used to set up new ``Analysis`` objects with a positions penalty c positions_likelihood=result_1.positions_likelihood_from( factor=3.0, minimum_threshold=0.2 ), + ) + +Multiple Source Plane Systems +----------------------------- + +A double source plane system is a lens system where there are mutiple source-planes at different redshifts, meaning that +incuding the image-plane there are at least 3 planes. + +The ``PositionsLH`` class can have a `plane_redshift` input, which specifies the redshift of the source-plane +the positions are ray-traced to. + +Multiple ``PositionsLH`` objects can be passed to the ``Analysis`` object, which then applies the penalty term to +both source-planes independently such that a double source-plane system can be fitted with the penalty based likelihood +functionality. + +.. code-block:: python + + positions_likelihood_source_plane_0 = al.PositionsLH(positions=positions, threshold=0.3, plane_redshift=1.0) + positions_likelihood_source_plane_1 = al.PositionsLH(positions=positions, threshold=0.3, plane_redshift=2.0) + + analysis = al.AnalysisImaging( + dataset=dataset, positions_likelihood_list= + [ + positions_likelihood_source_plane_0, + positions_likelihood_source_plane_1 + ] ) \ No newline at end of file From 26b9b01511981e6ac38ce6a63483298e188f6fc1 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 21 May 2025 16:28:48 +0100 Subject: [PATCH 11/11] source also updates based on plane redshift --- autolens/analysis/plotter_interface.py | 6 +- autolens/analysis/result.py | 73 ++++++++++++++----- autolens/exc.py | 5 +- autolens/imaging/plot/fit_imaging_plotters.py | 10 +-- docs/general/demagnified_solutions.rst | 19 ++++- test_autolens/analysis/test_result.py | 16 ++-- 6 files changed, 91 insertions(+), 38 deletions(-) diff --git a/autolens/analysis/plotter_interface.py b/autolens/analysis/plotter_interface.py index 7d790358d..3654401b6 100644 --- a/autolens/analysis/plotter_interface.py +++ b/autolens/analysis/plotter_interface.py @@ -139,7 +139,7 @@ def should_plot(name): def image_with_positions(self, image: aa.Array2D, positions: aa.Grid2DIrregular): """ - Visualizes the positions of a model-fit, where these positions are used to resample lens models where + Visualizes the positions of a model-fit, where these positions are used to penalize lens models where the positions to do trace within an input threshold of one another in the source-plane. Images are output to the `image` folder of the `image_path`. When used with a non-linear search the `image_path` @@ -154,8 +154,8 @@ def image_with_positions(self, image: aa.Array2D, positions: aa.Grid2DIrregular) ---------- imaging The imaging dataset whose image the positions are overlaid. - position - The 2D (y,x) arc-second positions used to resample inaccurate mass models. + positions + The 2D (y,x) arc-second positions used to penalize inaccurate mass models. """ def should_plot(name): diff --git a/autolens/analysis/result.py b/autolens/analysis/result.py index ab9a9875b..9f971f6dd 100644 --- a/autolens/analysis/result.py +++ b/autolens/analysis/result.py @@ -26,8 +26,9 @@ def max_log_likelihood_tracer(self) -> Tracer: """ return self.analysis.tracer_via_instance_from(instance=self.instance) - @property - def source_plane_light_profile_centre(self) -> aa.Grid2DIrregular: + def source_plane_light_profile_centre_from( + self, plane_redshift: Optional[float] = None + ) -> aa.Grid2DIrregular: """ Return a light profile centre of one of the a galaxies in the maximum log likelihood `Tracer`'s source-plane. If there are multiple light profiles, the first light profile's centre is returned. @@ -35,14 +36,22 @@ def source_plane_light_profile_centre(self) -> aa.Grid2DIrregular: These centres are used by automatic position updating to determine the best-fit lens model's image-plane multiple-image positions. """ - centre = self.max_log_likelihood_tracer.planes[-1].extract_attribute( + try: + plane_index = self.max_log_likelihood_tracer.plane_index_via_redshift_from( + redshift=plane_redshift + ) + except TypeError: + plane_index = -1 + + centre = self.max_log_likelihood_tracer.planes[plane_index].extract_attribute( cls=ag.LightProfile, attr_name="centre" ) if centre is not None: return aa.Grid2DIrregular(values=[np.asarray(centre[0])]) - @property - def source_plane_centre(self) -> aa.Grid2DIrregular: + def source_plane_centre_from( + self, plane_redshift: Optional[float] = None + ) -> aa.Grid2DIrregular: """ Return the centre of a source-plane galaxy via the following criteria: @@ -52,7 +61,9 @@ def source_plane_centre(self) -> aa.Grid2DIrregular: These centres are used by automatic position updating to determine the multiple-images of a best-fit lens model (and thus tracer) by back-tracing the centres to the image plane via the mass model. """ - return self.source_plane_light_profile_centre + return self.source_plane_light_profile_centre_from( + plane_redshift=plane_redshift + ) def image_plane_multiple_image_positions( self, plane_redshift: Optional[float] = None @@ -79,9 +90,13 @@ def image_plane_multiple_image_positions( pixel_scale_precision=0.001, ) + source_plane_centre = self.source_plane_centre_from( + plane_redshift=plane_redshift + ) + multiple_images = solver.solve( tracer=self.max_log_likelihood_tracer, - source_plane_coordinate=self.source_plane_centre.in_list[0], + source_plane_coordinate=source_plane_centre.in_list[0], plane_redshift=plane_redshift, ) @@ -270,7 +285,7 @@ def positions_likelihood_from( Returns ------- - The `PositionsLH` object used to apply a likelihood penalty or resample the positions. + The `PositionsLH` object used to apply a likelihood penalty using the positions. """ if os.environ.get("PYAUTOFIT_TEST_MODE") == "1": @@ -319,8 +334,9 @@ def max_log_likelihood_tracer(self) -> Tracer: return self.analysis.tracer_via_instance_from(instance=instance) - @property - def source_plane_centre(self) -> aa.Grid2DIrregular: + def source_plane_centre_from( + self, plane_redshift: Optional[float] = None + ) -> aa.Grid2DIrregular: """ Return the centre of a source-plane galaxy via the following criteria: @@ -330,23 +346,46 @@ def source_plane_centre(self) -> aa.Grid2DIrregular: These centres are used by automatic position updating to determine the multiple-images of a best-fit lens model (and thus tracer) by back-tracing the centres to the image plane via the mass model. """ - if self.source_plane_inversion_centre is not None: - return self.source_plane_inversion_centre - elif self.source_plane_light_profile_centre is not None: - return self.source_plane_light_profile_centre + source_plane_inversion_centre = self.source_plane_inversion_centre_from( + plane_redshift=plane_redshift + ) - @property - def source_plane_inversion_centre(self) -> aa.Grid2DIrregular: + if source_plane_inversion_centre is not None: + return source_plane_inversion_centre + + return self.source_plane_light_profile_centre_from( + plane_redshift=plane_redshift + ) + + def source_plane_inversion_centre_from( + self, plane_redshift: Optional[float] = None + ) -> aa.Grid2DIrregular: """ Returns the centre of the brightest source pixel(s) of an `LEq`. These centres are used by automatic position updating to determine the best-fit lens model's image-plane multiple-image positions. """ + if self.max_log_likelihood_fit.inversion is not None: if self.max_log_likelihood_fit.inversion.has(cls=aa.AbstractMapper): + + try: + plane_index = ( + self.max_log_likelihood_tracer.plane_index_via_redshift_from( + redshift=plane_redshift + ) + ) + mapper_index = ( + self.max_log_likelihood_tracer.plane_indexes_with_pixelizations[ + plane_index + ] + ) + except TypeError: + mapper_index = 0 + inversion = self.max_log_likelihood_fit.inversion - mapper = inversion.cls_list_from(cls=aa.AbstractMapper)[0] + mapper = inversion.cls_list_from(cls=aa.AbstractMapper)[mapper_index] mapper_valued = aa.MapperValued( values=inversion.reconstruction_dict[mapper], diff --git a/autolens/exc.py b/autolens/exc.py index 93a18b918..83a73449a 100644 --- a/autolens/exc.py +++ b/autolens/exc.py @@ -8,10 +8,7 @@ class RayTracingException(af.exc.FitException): """ Raises exceptions associated with the `lens/tracer.py` module and `Tracer` class. - For example if the multiple image positions do not trace without a threshold of one another, in order to - resample inaccurate mass models during a model-fit. - - This exception inehrits from a `FitException`. This means that if this exception is raised during a model-fit in + This exception inherits from a `FitException`. This means that if this exception is raised during a model-fit in the analysis class's `log_likelihood_function` that model is resampled and does not terminate the code. """ diff --git a/autolens/imaging/plot/fit_imaging_plotters.py b/autolens/imaging/plot/fit_imaging_plotters.py index bcbf4e3b0..ac958cacf 100644 --- a/autolens/imaging/plot/fit_imaging_plotters.py +++ b/autolens/imaging/plot/fit_imaging_plotters.py @@ -213,11 +213,11 @@ def figures_2d_of_planes( plane_indexes = self.plane_indexes_from(plane_index=plane_index) - if use_source_vmax: - self.mat_plot_2d.cmap.kwargs["vmax"] = np.max(self.fit.model_images_of_planes_list[-1]) - for plane_index in plane_indexes: + if use_source_vmax: + self.mat_plot_2d.cmap.kwargs["vmax"] = np.max(self.fit.model_images_of_planes_list[plane_index]) + if subtracted_image: title = f"Subtracted Image of Plane {plane_index}" @@ -765,7 +765,7 @@ def figures_2d( if data: if use_source_vmax: - self.mat_plot_2d.cmap.kwargs["vmax"] = np.max(self.fit.model_images_of_planes_list[-1]) + self.mat_plot_2d.cmap.kwargs["vmax"] = np.max(self.fit.model_images_of_planes_list[1:]) self.mat_plot_2d.plot_array( array=self.fit.data, @@ -799,7 +799,7 @@ def figures_2d( if model_image: if use_source_vmax: - self.mat_plot_2d.cmap.kwargs["vmax"] = np.max(self.fit.model_images_of_planes_list[-1]) + self.mat_plot_2d.cmap.kwargs["vmax"] = np.max(self.fit.model_images_of_planes_list[1:]) self.mat_plot_2d.plot_array( array=self.fit.model_data, diff --git a/docs/general/demagnified_solutions.rst b/docs/general/demagnified_solutions.rst index 9c4816b94..4d8c7936a 100644 --- a/docs/general/demagnified_solutions.rst +++ b/docs/general/demagnified_solutions.rst @@ -150,7 +150,7 @@ This is often used to set up new ``Analysis`` objects with a positions penalty c .. code-block:: python - analysis_2 = al.AnalysisImaging( + analysis = al.AnalysisImaging( dataset=dataset, positions_likelihood=result_1.positions_likelihood_from( factor=3.0, minimum_threshold=0.2 @@ -181,4 +181,21 @@ functionality. positions_likelihood_source_plane_0, positions_likelihood_source_plane_1 ] + ) + +To set up an ``Analysis`` object with multiple ``PositionsLH`` objects from a result, each positions likelihood is +computed from the result for each ``plane_redshift``: + +.. code-block:: + + analysis = al.AnalysisImaging( + dataset=dataset, + positions_likelihood_list=[ + source_lp_result.positions_likelihood_from( + factor=3.0, minimum_threshold=0.2, plane_redshift=source_lp_result.instance.galaxies.source_1.redshift, + ), + source_lp_result.positions_likelihood_from( + factor=3.0, minimum_threshold=0.2, plane_redshift=source_lp_result.instance.galaxies.source_2.redshift, + ) + ], ) \ No newline at end of file diff --git a/test_autolens/analysis/test_result.py b/test_autolens/analysis/test_result.py index 143961604..9c136b1db 100644 --- a/test_autolens/analysis/test_result.py +++ b/test_autolens/analysis/test_result.py @@ -44,7 +44,7 @@ def test__source_plane_light_profile_centre(analysis_imaging_7x7): result = res.Result(samples_summary=samples_summary, analysis=analysis_imaging_7x7) - assert result.source_plane_light_profile_centre.in_list == [(1.0, 2.0)] + assert result.source_plane_light_profile_centre_from().in_list == [(1.0, 2.0)] source_0 = al.Galaxy( redshift=1.0, @@ -62,7 +62,7 @@ def test__source_plane_light_profile_centre(analysis_imaging_7x7): result = res.Result(samples_summary=samples_summary, analysis=analysis_imaging_7x7) - assert result.source_plane_light_profile_centre.in_list == [(1.0, 2.0)] + assert result.source_plane_light_profile_centre_from().in_list == [(1.0, 2.0)] source_0 = al.Galaxy( redshift=1.0, light=al.lp.SersicSph(centre=(1.0, 2.0), intensity=2.0) @@ -78,7 +78,7 @@ def test__source_plane_light_profile_centre(analysis_imaging_7x7): result = res.Result(samples_summary=samples_summary, analysis=analysis_imaging_7x7) - assert result.source_plane_light_profile_centre.in_list == [(5.0, 6.0)] + assert result.source_plane_light_profile_centre_from().in_list == [(5.0, 6.0)] tracer = al.Tracer(galaxies=[al.Galaxy(redshift=0.5)]) @@ -86,7 +86,7 @@ def test__source_plane_light_profile_centre(analysis_imaging_7x7): result = res.Result(samples_summary=samples_summary, analysis=analysis_imaging_7x7) - assert result.source_plane_light_profile_centre == None + assert result.source_plane_light_profile_centre_from() == None def test__source_plane_inversion_centre(analysis_imaging_7x7): @@ -116,7 +116,7 @@ def test__source_plane_inversion_centre(analysis_imaging_7x7): ) assert ( - result.source_plane_inversion_centre.in_list[0] + result.source_plane_inversion_centre_from().in_list[0] == mapper_valued.max_pixel_centre.in_list[0] ) @@ -131,7 +131,7 @@ def test__source_plane_inversion_centre(analysis_imaging_7x7): samples_summary=samples_summary, analysis=analysis_imaging_7x7 ) - assert result.source_plane_inversion_centre == None + assert result.source_plane_inversion_centre_from() == None lens = al.Galaxy(redshift=0.5, light=al.lp_linear.Sersic()) source = al.Galaxy(redshift=1.0) @@ -144,7 +144,7 @@ def test__source_plane_inversion_centre(analysis_imaging_7x7): samples_summary=samples_summary, analysis=analysis_imaging_7x7 ) - assert result.source_plane_inversion_centre == None + assert result.source_plane_inversion_centre_from() == None def test__source_plane_centre(analysis_imaging_7x7): @@ -169,7 +169,7 @@ def test__source_plane_centre(analysis_imaging_7x7): samples_summary=samples_summary, analysis=analysis_imaging_7x7 ) - assert result.source_plane_centre.in_list[0] == pytest.approx( + assert result.source_plane_centre_from().in_list[0] == pytest.approx( (-0.916666, -0.916666), 1.0e-4 )