From 8394294ca3094a53f6ef977946f38cd299d5984d Mon Sep 17 00:00:00 2001 From: martin-millon Date: Fri, 30 Jul 2021 16:18:52 +0200 Subject: [PATCH 01/28] conf change --- lenstronomy/Conf/conf_default.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lenstronomy/Conf/conf_default.yaml b/lenstronomy/Conf/conf_default.yaml index c5b2526c0..763f9414f 100644 --- a/lenstronomy/Conf/conf_default.yaml +++ b/lenstronomy/Conf/conf_default.yaml @@ -12,5 +12,5 @@ numba: conventions: - sersic_major_axis: False # if True, defines the half-light radius of the Sersic light profile along the semi-major axis (which is the Galfit convention) + sersic_major_axis: True # if True, defines the half-light radius of the Sersic light profile along the semi-major axis (which is the Galfit convention) # if False, uses the product average of semi-major and semi-minor axis as the convention (default definition for all light profiles in lenstronomy other than the Sersic profile) \ No newline at end of file From cc58883888f69fe3308f9e4ba1f46b3a2ae494d9 Mon Sep 17 00:00:00 2001 From: martin-millon Date: Fri, 7 Jan 2022 16:59:59 +0100 Subject: [PATCH 02/28] adding Thin disk model --- lenstronomy/LightModel/Profiles/thin_disk.py | 45 ++++++++++++++++++++ lenstronomy/LightModel/light_model_base.py | 5 ++- lenstronomy/LightModel/linear_basis.py | 4 +- 3 files changed, 51 insertions(+), 3 deletions(-) create mode 100644 lenstronomy/LightModel/Profiles/thin_disk.py diff --git a/lenstronomy/LightModel/Profiles/thin_disk.py b/lenstronomy/LightModel/Profiles/thin_disk.py new file mode 100644 index 000000000..4433c2a86 --- /dev/null +++ b/lenstronomy/LightModel/Profiles/thin_disk.py @@ -0,0 +1,45 @@ +__author__ = 'martin-millon' + +# this file contains a class to make a thin disk profile (Shakura & Sunyaev 1973) +import numpy as np + +__all__ = ['ThinDisk'] + +class ThinDisk(object): + """ + this class contains functions to evaluate an thin disk model + + .. math:: + I(R) = I_0 / \\left[\\exp(\\xi) - 1\\right] + + with :math:`I_0 = amp` + and + with :math:`\\xi (R) = \\left[ R / R_0 \\right] ^{3/4} \\left[ 1 - \\sqrt{R_in/R}\\right] ^{-1/4}` + + """ + + param_names = ['amp', 'R_0', 'R_in', 'center_x', 'center_y'] + lower_limit_default = {'amp': 0, 'R_0': 0.0, 'R_in': 0.0, 'center_x': -100, 'center_y': -100} + upper_limit_default = {'amp': 100000, 'R_0': 100, 'R_in': 10, 'center_x': 100, 'center_y': 100} + + def function(self, x, y, amp, R_0, R_in, center_x=0, center_y=0): + """ + + :param x: + :param y: + :param amp: surface brightness/amplitude value at the half light radius + :param R_0: scale radius + :param R_in: inner radius of the profile + :param center_x: center in x-coordinate + :param center_y: center in y-coordinate + :return: Thin disk profile value at (x, y) + """ + R = np.sqrt((x - center_x) ** 2 + (y - center_y)**2) + profile = np.zeros_like(R) * 1e-10 + if R_in <= 0: + profile = amp / (np.exp((R / R_0) ** (3. / 4.)) - 1) + else: + profile = np.where(R <= R_in, profile, + amp / (np.exp((R / R_0) ** (0.75) * (1 - np.sqrt(R_in / R)) ** (-0.25)) - 1)) + + return profile diff --git a/lenstronomy/LightModel/light_model_base.py b/lenstronomy/LightModel/light_model_base.py index cb35bfaf6..6c3fa9de6 100644 --- a/lenstronomy/LightModel/light_model_base.py +++ b/lenstronomy/LightModel/light_model_base.py @@ -15,7 +15,7 @@ 'SERSIC', 'SERSIC_ELLIPSE', 'CORE_SERSIC', 'SHAPELETS', 'SHAPELETS_POLAR', 'SHAPELETS_POLAR_EXP', 'HERNQUIST', 'HERNQUIST_ELLIPSE', 'PJAFFE', 'PJAFFE_ELLIPSE', 'UNIFORM', 'POWER_LAW', 'NIE', 'CHAMELEON', 'DOUBLE_CHAMELEON', 'TRIPLE_CHAMELEON', 'INTERPOL', 'SLIT_STARLETS', - 'SLIT_STARLETS_GEN2'] + 'SLIT_STARLETS_GEN2', 'THINDISK'] class LightModelBase(object): @@ -108,6 +108,9 @@ def __init__(self, light_model_list, smoothing=0.001, sersic_major_axis=None): elif profile_type == 'SLIT_STARLETS_GEN2': from lenstronomy.LightModel.Profiles.starlets import SLIT_Starlets self.func_list.append(SLIT_Starlets(second_gen=True)) + elif profile_type == 'THINDISK': + from lenstronomy.LightModel.Profiles.thin_disk import ThinDisk + self.func_list.append(ThinDisk()) else: raise ValueError('No light model of type %s found! Supported are the following models: %s' % (profile_type, _MODELS_SUPPORTED)) self._num_func = len(self.func_list) diff --git a/lenstronomy/LightModel/linear_basis.py b/lenstronomy/LightModel/linear_basis.py index 3737e30e7..98bd63a37 100644 --- a/lenstronomy/LightModel/linear_basis.py +++ b/lenstronomy/LightModel/linear_basis.py @@ -47,7 +47,7 @@ def functions_split(self, x, y, kwargs_list, k=None): if bool_list[i] is True: if model in ['SERSIC', 'SERSIC_ELLIPSE', 'CORE_SERSIC', 'HERNQUIST', 'HERNQUIST_ELLIPSE', 'PJAFFE', 'PJAFFE_ELLIPSE', 'GAUSSIAN', 'GAUSSIAN_ELLIPSE', 'POWER_LAW', 'NIE', 'CHAMELEON', - 'DOUBLE_CHAMELEON', 'TRIPLE_CHAMELEON', 'UNIFORM', 'INTERPOL', 'ELLIPSOID']: + 'DOUBLE_CHAMELEON', 'TRIPLE_CHAMELEON', 'UNIFORM', 'INTERPOL', 'ELLIPSOID', 'THINDISK']: kwargs_new = kwargs_list[i].copy() new = {'amp': 1} kwargs_new.update(new) @@ -134,7 +134,7 @@ def update_linear(self, param, i, kwargs_list): for k, model in enumerate(self.profile_type_list): if model in ['SERSIC', 'SERSIC_ELLIPSE', 'CORE_SERSIC', 'HERNQUIST', 'PJAFFE', 'PJAFFE_ELLIPSE', 'HERNQUIST_ELLIPSE', 'GAUSSIAN', 'GAUSSIAN_ELLIPSE', 'POWER_LAW', 'NIE', 'CHAMELEON', - 'DOUBLE_CHAMELEON', 'TRIPLE_CHAMELEON', 'UNIFORM', 'INTERPOL', 'ELLIPSOID']: + 'DOUBLE_CHAMELEON', 'TRIPLE_CHAMELEON', 'UNIFORM', 'INTERPOL', 'ELLIPSOID', 'THINDISK']: kwargs_list[k]['amp'] = param[i] i += 1 elif model in ['MULTI_GAUSSIAN', 'MULTI_GAUSSIAN_ELLIPSE']: From 7a5bf6e61f8e121b4c8ebab730c3251ba5575c9f Mon Sep 17 00:00:00 2001 From: martin-millon Date: Thu, 20 Jan 2022 14:18:44 +0100 Subject: [PATCH 03/28] adding Thin disk model --- lenstronomy/LightModel/Profiles/thin_disk.py | 4 ++-- lenstronomy/Plots/plot_util.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/lenstronomy/LightModel/Profiles/thin_disk.py b/lenstronomy/LightModel/Profiles/thin_disk.py index 4433c2a86..e08d3792e 100644 --- a/lenstronomy/LightModel/Profiles/thin_disk.py +++ b/lenstronomy/LightModel/Profiles/thin_disk.py @@ -35,9 +35,9 @@ def function(self, x, y, amp, R_0, R_in, center_x=0, center_y=0): :return: Thin disk profile value at (x, y) """ R = np.sqrt((x - center_x) ** 2 + (y - center_y)**2) - profile = np.zeros_like(R) * 1e-10 + profile = np.zeros_like(R) if R_in <= 0: - profile = amp / (np.exp((R / R_0) ** (3. / 4.)) - 1) + profile = amp / (np.exp((R / R_0) ** (0.75)) - 1) else: profile = np.where(R <= R_in, profile, amp / (np.exp((R / R_0) ** (0.75) * (1 - np.sqrt(R_in / R)) ** (-0.25)) - 1)) diff --git a/lenstronomy/Plots/plot_util.py b/lenstronomy/Plots/plot_util.py index aea26980d..dc9ee1134 100644 --- a/lenstronomy/Plots/plot_util.py +++ b/lenstronomy/Plots/plot_util.py @@ -122,10 +122,10 @@ def plot_line_set(ax, coords, line_set_list_x, line_set_list_y, origin=None, fli if isinstance(line_set_list_x, list): for i in range(len(line_set_list_x)): x_c, y_c = coords.map_coord2pix(line_set_list_x[i], line_set_list_y[i]) - ax.plot((x_c + 0.5) * pixel_width_x + origin[0], (y_c + 0.5) * pixel_width + origin[1], *args, **kwargs) # ',', color=color) + ax.plot((x_c) * pixel_width_x + origin[0], (y_c) * pixel_width + origin[1], *args, **kwargs) # ',', color=color) else: x_c, y_c = coords.map_coord2pix(line_set_list_x, line_set_list_y) - ax.plot((x_c + 0.5) * pixel_width_x + origin[0], (y_c + 0.5) * pixel_width + origin[1], *args, **kwargs) # ',', color=color) + ax.plot((x_c) * pixel_width_x + origin[0], (y_c) * pixel_width + origin[1], *args, **kwargs) # ',', color=color) return ax From 8c7b69c7467bbd5688dcfb8c14d518da7a1da09f Mon Sep 17 00:00:00 2001 From: martin-millon Date: Mon, 7 Feb 2022 11:19:03 +0100 Subject: [PATCH 04/28] change lower limit for lens light in fitting sequence --- lenstronomy/Workflow/fitting_sequence.py | 6 ++++-- lenstronomy/Workflow/update_manager.py | 7 ++++++- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/lenstronomy/Workflow/fitting_sequence.py b/lenstronomy/Workflow/fitting_sequence.py index 4784250c6..47ca1173f 100644 --- a/lenstronomy/Workflow/fitting_sequence.py +++ b/lenstronomy/Workflow/fitting_sequence.py @@ -428,7 +428,8 @@ def update_settings(self, kwargs_model={}, kwargs_constraints={}, kwargs_likelih lens_remove_fixed=[], source_remove_fixed=[], lens_light_remove_fixed=[], ps_remove_fixed=[], cosmo_remove_fixed=[], change_source_lower_limit=None, change_source_upper_limit=None, - change_lens_lower_limit=None, change_lens_upper_limit=None): + change_lens_lower_limit=None, change_lens_upper_limit=None, + change_lens_light_lower_limit=None, change_lens_light_upper_limit=None): """ updates lenstronomy settings "on the fly" @@ -446,6 +447,7 @@ def update_settings(self, kwargs_model={}, kwargs_constraints={}, kwargs_likelih :param ps_remove_fixed: [[i_model, ['param1', 'param2',...], [...]] :param cosmo_remove_fixed: ['param1', 'param2',...] :param change_lens_lower_limit: [[i_model, ['param_name', ...], [value1, value2, ...]]] + :param change_lens_light_lower_limit: [[i_model, ['param_name', ...], [value1, value2, ...]]] :return: 0, the settings are overwritten for the next fitting step to come """ self._updateManager.update_options(kwargs_model, kwargs_constraints, kwargs_likelihood) @@ -453,7 +455,7 @@ def update_settings(self, kwargs_model={}, kwargs_constraints={}, kwargs_likelih ps_add_fixed, cosmo_add_fixed, lens_remove_fixed, source_remove_fixed, lens_light_remove_fixed, ps_remove_fixed, cosmo_remove_fixed) self._updateManager.update_limits(change_source_lower_limit, change_source_upper_limit, change_lens_lower_limit, - change_lens_upper_limit) + change_lens_upper_limit, change_lens_light_lower_limit, change_lens_light_upper_limit) return 0 def set_param_value(self, **kwargs): diff --git a/lenstronomy/Workflow/update_manager.py b/lenstronomy/Workflow/update_manager.py index 959a5790d..2f3fd8564 100644 --- a/lenstronomy/Workflow/update_manager.py +++ b/lenstronomy/Workflow/update_manager.py @@ -217,7 +217,8 @@ def update_options(self, kwargs_model, kwargs_constraints, kwargs_likelihood): return kwargs_model_updated, kwargs_constraints_updated, kwargs_likelihood_updated def update_limits(self, change_source_lower_limit=None, change_source_upper_limit=None, - change_lens_lower_limit=None, change_lens_upper_limit=None,): + change_lens_lower_limit=None, change_lens_upper_limit=None, + change_lens_light_lower_limit=None, change_lens_light_upper_limit=None): """ updates the limits (lower and upper) of the update manager instance @@ -235,6 +236,10 @@ def update_limits(self, change_source_lower_limit=None, change_source_upper_limi self._lens_lower = self._update_limit(change_lens_lower_limit, self._lens_lower) if not change_lens_upper_limit is None: self._lens_upper = self._update_limit(change_lens_upper_limit, self._lens_upper) + if not change_lens_light_lower_limit is None: + self._lens_light_lower = self._update_limit(change_lens_light_lower_limit, self._lens_light_lower) + if not change_lens_light_upper_limit is None: + self._lens_light_upper = self._update_limit(change_lens_light_upper_limit, self._lens_light_upper) @staticmethod def _update_limit(change_limit, kwargs_limit_previous): From 069bbb71f206404943075555b046ca47f10e7fd4 Mon Sep 17 00:00:00 2001 From: Martin Millon Date: Fri, 28 Oct 2022 16:47:29 -0700 Subject: [PATCH 05/28] Thin disk updates + plots --- lenstronomy/LightModel/Profiles/thin_disk.py | 100 ++++++++++++++++++- lenstronomy/LightModel/light_model_base.py | 8 +- lenstronomy/LightModel/linear_basis.py | 11 +- lenstronomy/Plots/lens_plot.py | 5 +- lenstronomy/Plots/model_band_plot.py | 1 + 5 files changed, 117 insertions(+), 8 deletions(-) diff --git a/lenstronomy/LightModel/Profiles/thin_disk.py b/lenstronomy/LightModel/Profiles/thin_disk.py index e08d3792e..ee6963c50 100644 --- a/lenstronomy/LightModel/Profiles/thin_disk.py +++ b/lenstronomy/LightModel/Profiles/thin_disk.py @@ -2,8 +2,9 @@ # this file contains a class to make a thin disk profile (Shakura & Sunyaev 1973) import numpy as np +from lenstronomy.Util import param_util -__all__ = ['ThinDisk'] +__all__ = ['ThinDisk', 'ThinDiskEllipse'] class ThinDisk(object): """ @@ -43,3 +44,100 @@ def function(self, x, y, amp, R_0, R_in, center_x=0, center_y=0): amp / (np.exp((R / R_0) ** (0.75) * (1 - np.sqrt(R_in / R)) ** (-0.25)) - 1)) return profile + +class ThinDiskEllipse(object): + """ + this class contains functions to evaluate an elliptical thin disk model + + """ + + param_names = ['amp', 'R_0', 'R_in', 'e1','e2','center_x', 'center_y'] + lower_limit_default = {'amp': 0, 'R_0': 0.0, 'R_in': 0.0, 'e1': -0.5, 'e2': -0.5,'center_x': -100, 'center_y': -100} + upper_limit_default = {'amp': 100000, 'R_0': 100, 'R_in': 10, 'e1': 0.5, 'e2':-0.5,'center_x': 100, 'center_y': 100} + + def function(self, x, y, amp, R_0, R_in, e1, e2, center_x=0, center_y=0): + """ + + :param x: + :param y: + :param amp: surface brightness/amplitude value at the half light radius + :param R_0: scale radius + :param R_in: inner radius of the profile + :param center_x: center in x-coordinate + :param center_y: center in y-coordinate + :param e1: eccentricity parameter + :param e2: eccentricity parameter + :return: Thin disk profile value at (x, y) + """ + + x_, y_ = param_util.transform_e1e2_product_average(x, y, e1, e2, center_x, center_y) + R_circular = np.sqrt((x - center_x) ** 2 + (y - center_y) ** 2) + R = np.sqrt((x_) ** 2 + (y_)**2) + profile = np.zeros_like(R) + if R_in <= 0: + profile = amp / (np.exp((R / R_0) ** (0.75)) - 1) + else: + profile = np.where(R <= R_in, profile, + amp / (np.exp((R / R_0) ** (0.75) * (1 - np.sqrt(R_in / R)) ** (-0.25)) - 1)) + + return profile + + +class ThinDiskEccentric(object): + """ + this class contains functions to evaluate an eccentric thin disk model (as in Eracleous et al. 1995) + + """ + + param_names = ['amp', 'R_0', 'R_in','e', 'phi_0','center_x', 'center_y'] + lower_limit_default = {'amp': 0, 'R_0': 0.0, 'R_in': 0.0, 'e': 0., 'phi_0': 0.,'center_x': -100, 'center_y': -100} + upper_limit_default = {'amp': 100000, 'R_0': 100, 'R_in': 10, 'e': 1.0, 'phi_0':2*np.pi,'center_x': 100, 'center_y': 100} + + def function(self, x, y, amp, R_0, R_in, e, phi_0, center_x=0, center_y=0): + """ + + :param x: + :param y: + :param amp: surface brightness/amplitude value at the half light radius + :param R_0: scale radius + :param R_in: inner radius of the profile + :param center_x: center in x-coordinate + :param center_y: center in y-coordinate + :param phi_0: angle : 0 corresponds to rthe apocenter along the x direction + :param e: eccentricity of the orbit + :return: Thin disk profile value at (x, y) + """ + + x_, y_ = self.eccentric_coordinate_system(x, y, e, phi_0, center_x, center_y) + R = np.sqrt((x_) ** 2 + (y_)**2) + #rescaling of the caracteristique scales + R_0 = R_0 / (1+e) + R_in = R_in / (1+e) + + profile = np.zeros_like(R) + if R_in <= 0: + profile = amp / (np.exp((R / R_0) ** (0.75)) - 1) + else: + profile = np.where(R <= R_in, profile, + amp / (np.exp((R / R_0) ** (0.75) * (1 - np.sqrt(R_in / R)) ** (-0.25)) - 1)) + + return profile + + def eccentric_coordinate_system(self, x, y, e, phi_0, center_x, center_y): + """ + maps the coordinates x, y into an eccentric coordinate system + + :param x: x-coordinate + :param y: y-coordinate + :param e: eccentricity + :param center_x: center of distortion + :param center_y: center of distortion + :return: distorted coordinates x', y' + """ + R = np.sqrt((x- center_x)**2 + (y-center_y)**2) + phi = np.arctan2((y-center_y) , (x-center_x)) + a = (R / (1 + e)) * (1 - e * np.cos(phi - phi_0)) + + xt1 = a * np.cos(phi) + xt2 = a * np.sin(phi) + return xt1, xt2 \ No newline at end of file diff --git a/lenstronomy/LightModel/light_model_base.py b/lenstronomy/LightModel/light_model_base.py index 15816b1a2..c5d131c1b 100644 --- a/lenstronomy/LightModel/light_model_base.py +++ b/lenstronomy/LightModel/light_model_base.py @@ -15,7 +15,7 @@ 'SERSIC', 'SERSIC_ELLIPSE', 'CORE_SERSIC', 'SHAPELETS', 'SHAPELETS_POLAR', 'SHAPELETS_POLAR_EXP', 'HERNQUIST', 'HERNQUIST_ELLIPSE', 'PJAFFE', 'PJAFFE_ELLIPSE', 'UNIFORM', 'POWER_LAW', 'NIE', 'CHAMELEON', 'DOUBLE_CHAMELEON', 'TRIPLE_CHAMELEON', 'INTERPOL', 'SLIT_STARLETS', - 'SLIT_STARLETS_GEN2', 'THINDISK'] + 'SLIT_STARLETS_GEN2', 'THINDISK', 'THINDISK_ELLIPSE','THINDISK_ECCENTRIC'] class LightModelBase(object): @@ -111,6 +111,12 @@ def __init__(self, light_model_list, smoothing=0.001, sersic_major_axis=None): elif profile_type == 'THINDISK': from lenstronomy.LightModel.Profiles.thin_disk import ThinDisk self.func_list.append(ThinDisk()) + elif profile_type == 'THINDISK_ELLIPSE': + from lenstronomy.LightModel.Profiles.thin_disk import ThinDiskEllipse + self.func_list.append(ThinDiskEllipse()) + elif profile_type == 'THINDISK_ECCENTRIC': + from lenstronomy.LightModel.Profiles.thin_disk import ThinDiskEccentric + self.func_list.append(ThinDiskEccentric()) else: raise ValueError('No light model of type %s found! Supported are the following models: %s' % (profile_type, _MODELS_SUPPORTED)) diff --git a/lenstronomy/LightModel/linear_basis.py b/lenstronomy/LightModel/linear_basis.py index 8783861e8..951023fd8 100644 --- a/lenstronomy/LightModel/linear_basis.py +++ b/lenstronomy/LightModel/linear_basis.py @@ -48,7 +48,8 @@ def functions_split(self, x, y, kwargs_list, k=None): if bool_list[i] is True: if model in ['SERSIC', 'SERSIC_ELLIPSE', 'CORE_SERSIC', 'HERNQUIST', 'HERNQUIST_ELLIPSE', 'PJAFFE', 'PJAFFE_ELLIPSE', 'GAUSSIAN', 'GAUSSIAN_ELLIPSE', 'POWER_LAW', 'NIE', 'CHAMELEON', - 'DOUBLE_CHAMELEON', 'TRIPLE_CHAMELEON', 'UNIFORM', 'INTERPOL', 'ELLIPSOID', 'THINDISK']: + 'DOUBLE_CHAMELEON', 'TRIPLE_CHAMELEON', 'UNIFORM', 'INTERPOL', 'ELLIPSOID', 'THINDISK', + 'THINDISK_ELLIPSE','THINDISK_ECCENTRIC']: kwargs_new = kwargs_list[i].copy() new = {'amp': 1} kwargs_new.update(new) @@ -102,7 +103,8 @@ def num_param_linear_list(self, kwargs_list): for i, model in enumerate(self.profile_type_list): if model in ['SERSIC', 'SERSIC_ELLIPSE', 'CORE_SERSIC', 'HERNQUIST', 'HERNQUIST_ELLIPSE', 'PJAFFE', 'PJAFFE_ELLIPSE', 'GAUSSIAN', 'GAUSSIAN_ELLIPSE', 'POWER_LAW', 'NIE', 'CHAMELEON', - 'DOUBLE_CHAMELEON', 'TRIPLE_CHAMELEON', 'UNIFORM', 'INTERPOL', 'ELLIPSOID']: + 'DOUBLE_CHAMELEON', 'TRIPLE_CHAMELEON', 'UNIFORM', 'INTERPOL', 'ELLIPSOID', 'THINDISK', + 'THINDISK_ELLIPSE','THINDISK_ECCENTRIC']: n_list += [1] elif model in ['MULTI_GAUSSIAN', 'MULTI_GAUSSIAN_ELLIPSE']: num = len(kwargs_list[i]['sigma']) @@ -135,7 +137,8 @@ def update_linear(self, param, i, kwargs_list): for k, model in enumerate(self.profile_type_list): if model in ['SERSIC', 'SERSIC_ELLIPSE', 'CORE_SERSIC', 'HERNQUIST', 'PJAFFE', 'PJAFFE_ELLIPSE', 'HERNQUIST_ELLIPSE', 'GAUSSIAN', 'GAUSSIAN_ELLIPSE', 'POWER_LAW', 'NIE', 'CHAMELEON', - 'DOUBLE_CHAMELEON', 'TRIPLE_CHAMELEON', 'UNIFORM', 'INTERPOL', 'ELLIPSOID', 'THINDISK']: + 'DOUBLE_CHAMELEON', 'TRIPLE_CHAMELEON', 'UNIFORM', 'INTERPOL', 'ELLIPSOID', 'THINDISK', + 'THINDISK_ELLIPSE','THINDISK_ECCENTRIC']: kwargs_list[k]['amp'] = param[i] i += 1 elif model in ['MULTI_GAUSSIAN', 'MULTI_GAUSSIAN_ELLIPSE']: @@ -187,7 +190,7 @@ def check_positive_flux_profile(self, kwargs_list): if 'amp' in kwargs_list[k]: if model in ['SERSIC', 'SERSIC_ELLIPSE', 'CORE_SERSIC', 'HERNQUIST', 'PJAFFE', 'PJAFFE_ELLIPSE', 'HERNQUIST_ELLIPSE', 'GAUSSIAN', 'GAUSSIAN_ELLIPSE', 'POWER_LAW', 'NIE', 'CHAMELEON', - 'DOUBLE_CHAMELEON']: + 'DOUBLE_CHAMELEON', 'THINDISK', 'THINDISK_ELLIPSE','THINDISK_ECCENTRIC']: if kwargs_list[k]['amp'] < 0: pos_bool = False break diff --git a/lenstronomy/Plots/lens_plot.py b/lenstronomy/Plots/lens_plot.py index 1b875edfa..78cbb5623 100644 --- a/lenstronomy/Plots/lens_plot.py +++ b/lenstronomy/Plots/lens_plot.py @@ -102,7 +102,7 @@ def convergence_plot(ax, pixel_grid, lens_model, kwargs_lens, extent=None, vmin= def caustics_plot(ax, pixel_grid, lens_model, kwargs_lens, fast_caustic=True, coord_inverse=False, color_crit='r', - color_caustic='g', pixel_offset=False, *args, **kwargs): + color_caustic='g', pixel_offset=False, plot_critical=True, *args, **kwargs): """ :param ax: matplotlib axis instance @@ -145,7 +145,8 @@ def caustics_plot(ax, pixel_grid, lens_model, kwargs_lens, fast_caustic=True, co plot_util.plot_line_set(ax, pixel_grid, ra_caustic_list, dec_caustic_list, color=color_caustic, origin=origin, flipped_x=coord_inverse, points_only=points_only, pixel_offset=pixel_offset, *args, **kwargs) - plot_util.plot_line_set(ax, pixel_grid, ra_crit_list, dec_crit_list, color=color_crit, origin=origin, + if plot_critical: + plot_util.plot_line_set(ax, pixel_grid, ra_crit_list, dec_crit_list, color=color_crit, origin=origin, flipped_x=coord_inverse, points_only=points_only, pixel_offset=pixel_offset, *args, **kwargs) return ax diff --git a/lenstronomy/Plots/model_band_plot.py b/lenstronomy/Plots/model_band_plot.py index faae1aed9..7787d3c5c 100644 --- a/lenstronomy/Plots/model_band_plot.py +++ b/lenstronomy/Plots/model_band_plot.py @@ -328,6 +328,7 @@ def source_plot(self, ax, numPix, deltaPix_source, center=None, v_min=None, cax = divider.append_axes("right", size="5%", pad=0.05) cb = plt.colorbar(im, cax=cax) cb.set_label(colorbar_label, fontsize=font_size) + plot_util.scale_bar(ax, d_s, dist=0.1, text='0.1"', color='w', flipped=False, font_size=font_size) if with_caustics is True: ra_caustic_list, dec_caustic_list = self._caustics() From 54066123f6799e88b9d715c7ec4853787f5e9212 Mon Sep 17 00:00:00 2001 From: Martin Millon Date: Wed, 1 Mar 2023 10:17:20 -0800 Subject: [PATCH 06/28] adding JWST config files --- .../SimulationAPI/ObservationConfig/JWST.py | 95 +++++++++++++++++++ 1 file changed, 95 insertions(+) create mode 100644 lenstronomy/SimulationAPI/ObservationConfig/JWST.py diff --git a/lenstronomy/SimulationAPI/ObservationConfig/JWST.py b/lenstronomy/SimulationAPI/ObservationConfig/JWST.py new file mode 100644 index 000000000..7bbc26ec7 --- /dev/null +++ b/lenstronomy/SimulationAPI/ObservationConfig/JWST.py @@ -0,0 +1,95 @@ +"""Provisional JWST instrument and observational settings. +ZP can be found here : https://jwst-docs.stsci.edu/files/182256933/182256934/1/1669487685625/NRC_ZPs_0995pmap.txt + +Sky Brightness needs to be derived from the ETC +""" +import lenstronomy.Util.util as util + +__all__ = ['JWST'] + + +NIRCAM_F200W_band_obs = {'exposure_time': 3600., + 'sky_brightness': 29.52, #this is derived using the ETC + 'magnitude_zero_point': 28.00, + #'detector': 'NRCA1', + 'num_exposures': 1, + 'seeing': None, + 'psf_type': 'PIXEL', + } + + +NIRCAM_F356W_band_obs = {'exposure_time': 3600., + 'sky_brightness': 28.39, #this is derived using the ETC + 'magnitude_zero_point': 26.47, + #'detector': 'NRCALONG', + 'num_exposures': 1, + 'seeing': None, + 'psf_type': 'PIXEL' # note kernel_point_source (the PSF map) must be provided separately + } + +""" +:keyword exposure_time: exposure time per image (in seconds) +:keyword sky_brightness: sky brightness (in magnitude per square arcseconds in units of electrons) +:keyword magnitude_zero_point: magnitude in which 1 count (e-) per second per arcsecond square is registered +:keyword num_exposures: number of exposures that are combined (depends on coadd_years) +:keyword seeing: Full-Width-at-Half-Maximum (FWHM) of PSF +:keyword psf_type: string, type of PSF ('GAUSSIAN' and 'PIXEL' supported) +""" + + +class JWST(object): + """ + class contains JWST instrument and observation configurations + """ + + def __init__(self, band='F200W', psf_type='PIXEL', coadd_years=None): + """ + + :param band: string, 'F200W' or 'F356W' supported. Determines obs dictionary. + :param psf_type: string, type of PSF ('GAUSSIAN', 'PIXEL' supported). + :param coadd_years: int, number of years corresponding to num_exposures in obs dict. Currently supported: None. + """ + + if band == 'F200W': + self.obs = NIRCAM_F200W_band_obs + self.arm = 'short' + elif band == 'F356W': + self.obs = NIRCAM_F356W_band_obs + self.arm = 'long' + + else: + raise ValueError("band %s not supported!"% band) + + if psf_type == 'GAUSSIAN': + self.obs['psf_type'] = 'GAUSSIAN' + elif psf_type != 'PIXEL': + raise ValueError("psf_type %s not supported!" % psf_type) + + if coadd_years is not None: + raise ValueError(" %s coadd_years not supported! " + "You may manually adjust num_exposures in obs dict if required." % coadd_years) + + # NIRCAM camera settings + if self.arm == 'short': + self.camera = {'read_noise': 15.77, + 'pixel_scale': 0.031, + 'ccd_gain': 2.05, + } + elif self.arm == 'long': + self.camera = {'read_noise': 13.25, + 'pixel_scale': 0.063, + 'ccd_gain': 1.82, + } + """ + :keyword read_noise: std of noise generated by read-out (in units of electrons) + :keyword pixel_scale: scale (in arcseconds) of pixels + :keyword ccd_gain: electrons/ADU (analog-to-digital unit). + """ + + def kwargs_single_band(self): + """ + + :return: merged kwargs from camera and obs dicts + """ + kwargs = util.merge_dicts(self.camera, self.obs) + return kwargs From b50e1c57c4630f0916d78409fc3fd20b2cf80cea Mon Sep 17 00:00:00 2001 From: Martin Millon Date: Wed, 1 Mar 2023 10:31:43 -0800 Subject: [PATCH 07/28] adding tests --- .../test_ObservationConfig/test_JWST.py | 67 +++++++++++++++++++ 1 file changed, 67 insertions(+) create mode 100644 test/test_SimulationAPI/test_ObservationConfig/test_JWST.py diff --git a/test/test_SimulationAPI/test_ObservationConfig/test_JWST.py b/test/test_SimulationAPI/test_ObservationConfig/test_JWST.py new file mode 100644 index 000000000..249557809 --- /dev/null +++ b/test/test_SimulationAPI/test_ObservationConfig/test_JWST.py @@ -0,0 +1,67 @@ +import unittest +from lenstronomy.SimulationAPI.ObservationConfig.JWST import JWST +from lenstronomy.SimulationAPI.observation_api import Instrument, SingleBand +import lenstronomy.Util.util as util + + +class TestJWST(unittest.TestCase): + + def setUp(self): + self.F200W = JWST() # default is F200W + self.F356W = JWST(band='F356W') + self.F356W2 = JWST(band='F356W', psf_type='GAUSSIAN') + + kwargs_F200W = self.F200W.kwargs_single_band() + kwargs_F356W = self.F356W.kwargs_single_band() + kwargs_F356W2 = self.F356W2.kwargs_single_band() + + self.F200W_band = SingleBand(**kwargs_F200W) + self.F356W_band = SingleBand(**kwargs_F356W) + self.F356W2_band = SingleBand(**kwargs_F356W2) + + # dictionaries mapping JWST kwargs to SingleBand kwargs + self.camera_settings = {'read_noise': '_read_noise', + 'pixel_scale': 'pixel_scale', + 'ccd_gain': 'ccd_gain'} + self.obs_settings = {'exposure_time': '_exposure_time', + 'sky_brightness': '_sky_brightness_', + 'magnitude_zero_point': '_magnitude_zero_point', + 'num_exposures': '_num_exposures', + 'seeing': '_seeing', + 'psf_type': '_psf_type'} + + self.instrument = Instrument(**self.F200W.camera) + + def test_JWST_class(self): + default = self.F200W + explicit_F200W = JWST(band='F200W', psf_type='PIXEL') + self.assertEqual(explicit_F200W.camera, default.camera) + self.assertEqual(explicit_F200W.obs, default.obs) + + with self.assertRaises(ValueError): + bad_band = JWST(band='g') + + with self.assertRaises(ValueError): + bad_psf = JWST(psf_type='blah') + + with self.assertRaises(ValueError): + bad_coadd_years = JWST(coadd_years=100) + + def test_JWST_camera(self): + # comparing camera settings in JWST instance with those in Instrument instance + for config, setting in self.camera_settings.items(): + self.assertEqual(self.F200W.camera[config], getattr(self.instrument, setting), msg=f"{config} did not match") + + def test_JWST_obs(self): + # comparing obs settings in JWST instance with those in SingleBand instance + for config, setting in self.obs_settings.items(): + self.assertEqual(self.F200W.obs[config], getattr(self.F200W_band, setting), msg=f"{config} did not match") + self.assertEqual(self.F356W.obs[config], getattr(self.F356W_band, setting), msg=f"{config} did not match") + self.assertEqual(self.F356W2.obs[config], getattr(self.F356W2_band, setting), msg=f"{config} did not match") + + def test_kwargs_single_band(self): + kwargs_F200W = util.merge_dicts(self.F200W.camera, self.F200W.obs) + self.assertEqual(self.F200W.kwargs_single_band(), kwargs_F200W) + +if __name__ == '__main__': + unittest.main() From 489b6087f91e03ad38d59bebbc5f8e56d4beb1c2 Mon Sep 17 00:00:00 2001 From: martin-millon Date: Tue, 9 Jan 2024 15:48:41 +0100 Subject: [PATCH 08/28] changing search windows --- lenstronomy/ImSim/image_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lenstronomy/ImSim/image_model.py b/lenstronomy/ImSim/image_model.py index c18751b4b..233dc4128 100644 --- a/lenstronomy/ImSim/image_model.py +++ b/lenstronomy/ImSim/image_model.py @@ -68,7 +68,7 @@ def __init__( search_window=np.max(self.Data.width), x_center=x_center, y_center=y_center, - min_distance=self.Data.pixel_width, + min_distance=self.Data.pixel_width/4., only_from_unspecified=True, ) self._psf_error_map = self.PSF.psf_error_map_bool From c4865a1ed3c7496d803110c73dc81ba804798361 Mon Sep 17 00:00:00 2001 From: martin-millon Date: Sat, 18 May 2024 14:31:58 -0700 Subject: [PATCH 09/28] typo fix --- lenstronomy/LightModel/Profiles/thin_disk.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lenstronomy/LightModel/Profiles/thin_disk.py b/lenstronomy/LightModel/Profiles/thin_disk.py index ee6963c50..7528c218b 100644 --- a/lenstronomy/LightModel/Profiles/thin_disk.py +++ b/lenstronomy/LightModel/Profiles/thin_disk.py @@ -4,11 +4,11 @@ import numpy as np from lenstronomy.Util import param_util -__all__ = ['ThinDisk', 'ThinDiskEllipse'] +__all__ = ['ThinDisk', 'ThinDiskEllipse', 'ThinDiskEccentric'] class ThinDisk(object): """ - this class contains functions to evaluate an thin disk model + This class contains functions to evaluate a thin disk model .. math:: I(R) = I_0 / \\left[\\exp(\\xi) - 1\\right] @@ -135,7 +135,7 @@ def eccentric_coordinate_system(self, x, y, e, phi_0, center_x, center_y): :return: distorted coordinates x', y' """ R = np.sqrt((x- center_x)**2 + (y-center_y)**2) - phi = np.arctan2((y-center_y) , (x-center_x)) + phi = np.arctan2((y-center_y), (x-center_x)) a = (R / (1 + e)) * (1 - e * np.cos(phi - phi_0)) xt1 = a * np.cos(phi) From f21bfc72cbf7903f567d7b2fa96f350aa9d35355 Mon Sep 17 00:00:00 2001 From: martin-millon Date: Wed, 29 May 2024 11:33:13 -0700 Subject: [PATCH 10/28] first STARRED implementation --- lenstronomy/Workflow/psf_fitting.py | 230 ++++++++++++++++++------- test/test_Workflow/test_psf_fitting.py | 37 ++++ test_requirements.txt | 1 + 3 files changed, 201 insertions(+), 67 deletions(-) diff --git a/lenstronomy/Workflow/psf_fitting.py b/lenstronomy/Workflow/psf_fitting.py index 2c2de1388..29058b1c2 100644 --- a/lenstronomy/Workflow/psf_fitting.py +++ b/lenstronomy/Workflow/psf_fitting.py @@ -1,3 +1,4 @@ +import matplotlib.pyplot as plt from lenstronomy.Data.psf import PSF import lenstronomy.Util.util as util import lenstronomy.Util.image_util as image_util @@ -7,6 +8,7 @@ import numpy as np import copy from scipy import ndimage +import warnings __all__ = ["PsfFitting"] @@ -75,14 +77,14 @@ def calc_cornermask(kernelsize, psf_symmetry): return mask def update_iterative( - self, - kwargs_psf, - kwargs_params, - num_iter=10, - keep_psf_error_map=True, - no_break=True, - verbose=True, - **kwargs_psf_update + self, + kwargs_psf, + kwargs_params, + num_iter=10, + keep_psf_error_map=True, + no_break=True, + verbose=True, + **kwargs_psf_update ): """ @@ -131,6 +133,7 @@ def update_iterative( kwargs_psf_new, _kwargs_params, corner_mask, **kwargs_psf_update ) + print(logL_best, logL_after) if logL_after > logL_best: kwargs_psf_final = copy.deepcopy(kwargs_psf_new) error_map_final = copy.deepcopy(error_map) @@ -158,24 +161,27 @@ def update_iterative( return kwargs_psf_final def update_psf( - self, - kwargs_psf, - kwargs_params, - corner_mask=None, - stacking_method="median", - psf_symmetry=1, - psf_iter_factor=0.2, - block_center_neighbour=0, - error_map_radius=None, - block_center_neighbour_error_map=None, - new_procedure=True, - corner_symmetry=None, + self, + kwargs_psf, + kwargs_params, + corner_mask=None, + stacking_method="median", + psf_symmetry=1, + psf_iter_factor=0.2, + block_center_neighbour=0, + error_map_radius=None, + block_center_neighbour_error_map=None, + new_procedure=True, + corner_symmetry=None, + use_starred=False, + kwargs_starred=None, ): """ :param kwargs_psf: keyword arguments to construct the PSF() class :param kwargs_params: keyword arguments of the parameters of the model components (e.g. 'kwargs_lens' etc) :param stacking_method: 'median', 'mean'; the different estimates of the PSF are stacked and combined together. + Not used if use_starred is True. The choices are: 'mean': mean of pixel values as the estimator (not robust to outliers) 'median': median of pixel values as the estimator (outlier rejection robust but needs >2 point sources in the @@ -203,6 +209,10 @@ def update_psf( 2) creates a psf with no symmetry 3) adds the corner_symmetry psf (which has information at the corners) to the odd symmetry PSF, in the regions where the odd-symmetry PSF does not have information + :param use_starred: boolean, if True, uses the STARRED method to estimate the PSF (https://ui.adsabs.harvard.edu/abs/2023JOSS....8.5340M/abstract) + :param kwargs_starred: dictionary, keyword arguments for the starred.procedures.psf_routines.update_PSF() method. + Example: kwargs_starred = {'lambda_scales':2., 'lambda_hf':2., 'lambda_positivity':0.} + :return: kwargs_psf_new, logL_after, error_map """ if block_center_neighbour_error_map is None: @@ -243,7 +253,8 @@ def update_psf( kernel_old = psf_class.kernel_point_source kernel_size = len(kernel_old) - if not new_procedure: + # TODO: refractor this if new_procedure is adopted + if not new_procedure and not use_starred: image_single_point_source_list = self.image_single_point_source( self._image_model_class, kwargs_params ) @@ -297,7 +308,7 @@ def update_psf( n_kernel = len(kwargs_psf["kernel_point_source"]) kernel_new = kernel_util.cut_psf(kernel_new, psf_size=n_kernel) - else: + elif not use_starred: kernel_old_high_res = psf_class.kernel_point_source_supersampled( supersampling_factor=point_source_supersampling_factor ) @@ -339,14 +350,99 @@ def update_psf( point_source_supersampling_factor, error_map_radius=error_map_radius, ) + else: + # use STARRED empirical PSF reconstruction routines + try: + from starred.procedures.psf_routines import run_multi_steps_PSF_reconstruction + from starred.psf.psf import PSF as StarredPSF + from starred.psf.parameters import ParametersPSF + except ImportError: + raise ImportError("STARRED is not installed. Please install STARRED to use this feature. " + "STARRED available by typing 'pip install starred-astro'.") + + image_single_point_source_list = self.image_single_point_source( + self._image_model_class, kwargs_params + ) + + # get the non-aligned cutouts + psf_kernel_list = self.point_like_source_cutouts( + x_pos=x_, + y_pos=y_, + image_list=image_single_point_source_list, + cutout_size=kernel_size, + ) + # get the noise maps + sigma2_maps_list = self.point_like_source_cutouts( + x_pos=x_, + y_pos=y_, + image_list=[self._image_model_class.Data.C_D for _ in range(len(x_))], + cutout_size=kernel_size, + ) + + #STARRED does not like 0 or negative value in the noise maps... so we replace them with the median + sigma2_maps_list = np.asarray(sigma2_maps_list) + psf_kernel_list = np.asarray(psf_kernel_list) + #normalise the data + norm = psf_kernel_list[0].max() + psf_kernel_list /= norm + sigma2_maps_list /= norm ** 2 + + indneg = np.where(sigma2_maps_list <= 0) + if len(indneg[0]) > 0: + warnings.warn("Negative or zero values in the noise maps. Replacing these pixels with the median value.") + sigma2_maps_list[indneg] = np.median(sigma2_maps_list) + + # ideally the narrow PSF kernel should be saved for the next iteration + # todo: masks are not used here, but could be implemented in the future + N, image_size, _ = np.shape(psf_kernel_list) + + #setup the STARRED model + model = StarredPSF(image_size=image_size, number_of_sources=N, + upsampling_factor=point_source_supersampling_factor, elliptical_moffat=True) + + kwargs_init, kwargs_fixed, kwargs_up, kwargs_down = model.smart_guess(psf_kernel_list, fixed_background=False) + parameters = ParametersPSF(kwargs_init, kwargs_fixed, kwargs_up=kwargs_up, kwargs_down=kwargs_down) + + model, parameters, loss, kwargs_partial_list, LogL_list, loss_history_list = run_multi_steps_PSF_reconstruction( + psf_kernel_list, model, + parameters, sigma2_maps_list, masks=None, **kwargs_starred) + + # kernel_new, narrow_psf_new, psf_kernel_list, _ = starred.procedures.psf_routines.update_PSF(kernel_old, + # psf_kernel_list, + # sigma2_maps_list, + # masks=None, + # subsampling_factor=point_source_supersampling_factor, + # **kwargs_starred) + + psf_kernel_list = model.model(**kwargs_partial_list[-1], high_res=False) #return model for each point source at the resolution of the data + psf_kernel_list = [psf_k / np.sum(psf_k) for psf_k in psf_kernel_list] + kernel_new = model.get_full_psf(**kwargs_partial_list[-1], norm=True, high_res=True) #subsampled PSF + kernel_new = psf_iter_factor * kernel_new + (1.0 - psf_iter_factor) * kernel_old + + error_map = self.error_map_estimate_new( + kernel_new, + psf_kernel_list, + ra_image, + dec_image, + point_amp, + point_source_supersampling_factor, + error_map_radius=error_map_radius, + ) + + # fig, ax = plt.subplots(1,3, figsize=(15,12)) + # ax[0].imshow(np.log10(kernel_new)) + # ax[1].imshow(np.log10(kernel_old)) + # ax[2].imshow(kernel_new- kernel_old) + # plt.show() + kwargs_psf_new["kernel_point_source"] = kernel_new - # if 'psf_error_map' in kwargs_psf_new: - # kwargs_psf_new['psf_error_map'] *= 10 self._image_model_class.update_psf(PSF(**kwargs_psf_new)) logL_after, _ = self._image_model_class.likelihood_data_given_model( **kwargs_params ) + if use_starred: + logL_after =0. return kwargs_psf_new, logL_after, error_map def image_single_point_source(self, image_model_class, kwargs_params): @@ -397,15 +493,15 @@ def _point_sources_list(image_model_class, kwargs_ps, kwargs_lens, k=None): return point_list def cutout_psf( - self, - ra_image, - dec_image, - x, - y, - image_list, - kernel_size, - kernel_init, - block_center_neighbour=0, + self, + ra_image, + dec_image, + x, + y, + image_list, + kernel_size, + kernel_init, + block_center_neighbour=0, ): """ @@ -440,15 +536,15 @@ def cutout_psf( return kernel_list def psf_estimate_individual( - self, - ra_image, - dec_image, - point_amp, - residuals, - cutout_size, - kernel_guess, - supersampling_factor, - block_center_neighbour, + self, + ra_image, + dec_image, + point_amp, + residuals, + cutout_size, + kernel_guess, + supersampling_factor, + block_center_neighbour, ): """ @@ -597,13 +693,13 @@ def cutout_psf_single(x, y, image, mask, kernel_size, kernel_init): @staticmethod ####Correct this based on Maverick Oh's note about lack of flux conservation. def combine_psf( - kernel_list_new, - kernel_old, - factor=1.0, - stacking_option="median", - symmetry=1, - corner_symmetry=None, - corner_mask=None, + kernel_list_new, + kernel_old, + factor=1.0, + stacking_option="median", + symmetry=1, + corner_symmetry=None, + corner_mask=None, ): ## TODO: Account for image edges/masked pixels. """Updates psf estimate based on old kernel and several new estimates. @@ -694,14 +790,14 @@ def combine_psf( return kernel_return def error_map_estimate_new( - self, - psf_kernel, - psf_kernel_list, - ra_image, - dec_image, - point_amp, - supersampling_factor, - error_map_radius=None, + self, + psf_kernel, + psf_kernel_list, + ra_image, + dec_image, + point_amp, + supersampling_factor, + error_map_radius=None, ): """Relative uncertainty in the psf model (in quadrature) per pixel based on residuals achieved in the image. @@ -741,7 +837,7 @@ def error_map_estimate_new( residuals_i = np.abs(kernel_low - kernel_low_i) residuals_i -= np.sqrt(C_D_cutout) / amp_i residuals_i[residuals_i < 0] = 0 - error_map_list[i, :, :] = residuals_i**2 + error_map_list[i, :, :] = residuals_i ** 2 error_map = np.median(error_map_list, axis=0) error_map[kernel_low > 0] /= kernel_low[kernel_low > 0] ** 2 @@ -759,14 +855,14 @@ def error_map_estimate_new( return error_map def error_map_estimate( - self, - kernel, - star_cutout_list, - amp, - x_pos, - y_pos, - error_map_radius=None, - block_center_neighbour=0, + self, + kernel, + star_cutout_list, + amp, + x_pos, + y_pos, + error_map_radius=None, + block_center_neighbour=0, ): """Provides a psf_error_map based on the goodness of fit of the given PSF kernel on the point source cutouts, their estimated amplitudes and positions. @@ -822,7 +918,7 @@ def error_map_estimate( residual[residual < 0] = 0 # estimate relative error per star residual /= amp_i - error_map_list[i, :, :] = residual**2 * mask_i + error_map_list[i, :, :] = residual ** 2 * mask_i mask_list[i, :, :] = mask_i # take median absolute error for each pixel # TODO: only for pixels that are not masked diff --git a/test/test_Workflow/test_psf_fitting.py b/test/test_Workflow/test_psf_fitting.py index 0d98f0256..03aeb2428 100644 --- a/test/test_Workflow/test_psf_fitting.py +++ b/test/test_Workflow/test_psf_fitting.py @@ -1,5 +1,6 @@ __author__ = "sibirrer" +import matplotlib.pyplot as plt import pytest import numpy as np import copy @@ -178,6 +179,42 @@ def test_update_psf(self): diff_new = np.sum((kernel_new - kernel_true) ** 2) assert diff_old > diff_new + #test STARRED + kwargs_psf_iter_starred = { + "stacking_method": "median", + "error_map_radius": 0.5, + "psf_iter_factor": 0.2, + "new_procedure": False, + "use_starred": True, + "kwargs_starred": {'verbose':False, 'lambda_scales':3, 'lambda_hf':3}, + } + + kwargs_psf_return_starred, improved_bool_starred, error_map_starred = self.psf_fitting.update_psf( + kwargs_psf, self.kwargs_params, **kwargs_psf_iter_starred + ) + assert improved_bool_starred + # fig, ax = plt.subplots(2,3) + # ax[0,0].imshow(kwargs_psf_return_starred["kernel_point_source"]) + # ax[0,1].imshow(kernel_old) + # ax[0,2].imshow(kwargs_psf_return_starred["kernel_point_source"]- kernel_old) + # + # ax[1,0].imshow(kwargs_psf_return_starred["kernel_point_source"]) + # ax[1,1].imshow(kernel_true) + # ax[1,2].imshow(kwargs_psf_return_starred["kernel_point_source"] - kernel_true) + # plt.show() + # + # fig, ax = plt.subplots(1,3) + # ax[0].imshow(kwargs_psf_return_starred["kernel_point_source"]) + # ax[1].imshow(error_map_starred) + # ax[2].imshow(error_map) + # plt.show() + + kernel_new_starred = kwargs_psf_return_starred["kernel_point_source"] + diff_new_starred = np.sum((kernel_new_starred- kernel_true) ** 2) + + # print(diff_new_starred, diff_new, diff_old) + assert diff_old > diff_new_starred + def test_calc_corner_mask(self): kernel_old = np.ones((101, 101)) nsymmetry = 4 diff --git a/test_requirements.txt b/test_requirements.txt index 7ca73df3f..08e3c232d 100644 --- a/test_requirements.txt +++ b/test_requirements.txt @@ -2,3 +2,4 @@ colossus git+https://github.com/mattgomer/SKiNN@main#egg=SKiNN git+https://github.com/aymgal/coolest@main#egg=coolest +git+https://gitlab.com/cosmograil/starred.git@main#egg=starred From 2ea497a65a312bc4203333102f1fe67ba0ae7464 Mon Sep 17 00:00:00 2001 From: martin-millon Date: Wed, 29 May 2024 12:15:29 -0700 Subject: [PATCH 11/28] remove specifics from my fork --- lenstronomy/Conf/conf_default.yaml | 2 +- lenstronomy/LightModel/Profiles/thin_disk.py | 143 ------------------ lenstronomy/Plots/model_band_plot.py | 1 - .../SimulationAPI/ObservationConfig/JWST.py | 122 +++++++-------- .../test_ObservationConfig/test_JWST.py | 62 +++++--- 5 files changed, 104 insertions(+), 226 deletions(-) delete mode 100644 lenstronomy/LightModel/Profiles/thin_disk.py diff --git a/lenstronomy/Conf/conf_default.yaml b/lenstronomy/Conf/conf_default.yaml index 6131b766d..7359015d9 100644 --- a/lenstronomy/Conf/conf_default.yaml +++ b/lenstronomy/Conf/conf_default.yaml @@ -12,6 +12,6 @@ numba: conventions: - sersic_major_axis: True # if True, defines the half-light radius of the Sersic light profile along the semi-major axis (which is the Galfit convention) + sersic_major_axis: False # if True, defines the half-light radius of the Sersic light profile along the semi-major axis (which is the Galfit convention) # if False, uses the product average of semi-major and semi-minor axis as the convention (default definition for all light profiles in lenstronomy other than the Sersic profile) diff --git a/lenstronomy/LightModel/Profiles/thin_disk.py b/lenstronomy/LightModel/Profiles/thin_disk.py deleted file mode 100644 index 7528c218b..000000000 --- a/lenstronomy/LightModel/Profiles/thin_disk.py +++ /dev/null @@ -1,143 +0,0 @@ -__author__ = 'martin-millon' - -# this file contains a class to make a thin disk profile (Shakura & Sunyaev 1973) -import numpy as np -from lenstronomy.Util import param_util - -__all__ = ['ThinDisk', 'ThinDiskEllipse', 'ThinDiskEccentric'] - -class ThinDisk(object): - """ - This class contains functions to evaluate a thin disk model - - .. math:: - I(R) = I_0 / \\left[\\exp(\\xi) - 1\\right] - - with :math:`I_0 = amp` - and - with :math:`\\xi (R) = \\left[ R / R_0 \\right] ^{3/4} \\left[ 1 - \\sqrt{R_in/R}\\right] ^{-1/4}` - - """ - - param_names = ['amp', 'R_0', 'R_in', 'center_x', 'center_y'] - lower_limit_default = {'amp': 0, 'R_0': 0.0, 'R_in': 0.0, 'center_x': -100, 'center_y': -100} - upper_limit_default = {'amp': 100000, 'R_0': 100, 'R_in': 10, 'center_x': 100, 'center_y': 100} - - def function(self, x, y, amp, R_0, R_in, center_x=0, center_y=0): - """ - - :param x: - :param y: - :param amp: surface brightness/amplitude value at the half light radius - :param R_0: scale radius - :param R_in: inner radius of the profile - :param center_x: center in x-coordinate - :param center_y: center in y-coordinate - :return: Thin disk profile value at (x, y) - """ - R = np.sqrt((x - center_x) ** 2 + (y - center_y)**2) - profile = np.zeros_like(R) - if R_in <= 0: - profile = amp / (np.exp((R / R_0) ** (0.75)) - 1) - else: - profile = np.where(R <= R_in, profile, - amp / (np.exp((R / R_0) ** (0.75) * (1 - np.sqrt(R_in / R)) ** (-0.25)) - 1)) - - return profile - -class ThinDiskEllipse(object): - """ - this class contains functions to evaluate an elliptical thin disk model - - """ - - param_names = ['amp', 'R_0', 'R_in', 'e1','e2','center_x', 'center_y'] - lower_limit_default = {'amp': 0, 'R_0': 0.0, 'R_in': 0.0, 'e1': -0.5, 'e2': -0.5,'center_x': -100, 'center_y': -100} - upper_limit_default = {'amp': 100000, 'R_0': 100, 'R_in': 10, 'e1': 0.5, 'e2':-0.5,'center_x': 100, 'center_y': 100} - - def function(self, x, y, amp, R_0, R_in, e1, e2, center_x=0, center_y=0): - """ - - :param x: - :param y: - :param amp: surface brightness/amplitude value at the half light radius - :param R_0: scale radius - :param R_in: inner radius of the profile - :param center_x: center in x-coordinate - :param center_y: center in y-coordinate - :param e1: eccentricity parameter - :param e2: eccentricity parameter - :return: Thin disk profile value at (x, y) - """ - - x_, y_ = param_util.transform_e1e2_product_average(x, y, e1, e2, center_x, center_y) - R_circular = np.sqrt((x - center_x) ** 2 + (y - center_y) ** 2) - R = np.sqrt((x_) ** 2 + (y_)**2) - profile = np.zeros_like(R) - if R_in <= 0: - profile = amp / (np.exp((R / R_0) ** (0.75)) - 1) - else: - profile = np.where(R <= R_in, profile, - amp / (np.exp((R / R_0) ** (0.75) * (1 - np.sqrt(R_in / R)) ** (-0.25)) - 1)) - - return profile - - -class ThinDiskEccentric(object): - """ - this class contains functions to evaluate an eccentric thin disk model (as in Eracleous et al. 1995) - - """ - - param_names = ['amp', 'R_0', 'R_in','e', 'phi_0','center_x', 'center_y'] - lower_limit_default = {'amp': 0, 'R_0': 0.0, 'R_in': 0.0, 'e': 0., 'phi_0': 0.,'center_x': -100, 'center_y': -100} - upper_limit_default = {'amp': 100000, 'R_0': 100, 'R_in': 10, 'e': 1.0, 'phi_0':2*np.pi,'center_x': 100, 'center_y': 100} - - def function(self, x, y, amp, R_0, R_in, e, phi_0, center_x=0, center_y=0): - """ - - :param x: - :param y: - :param amp: surface brightness/amplitude value at the half light radius - :param R_0: scale radius - :param R_in: inner radius of the profile - :param center_x: center in x-coordinate - :param center_y: center in y-coordinate - :param phi_0: angle : 0 corresponds to rthe apocenter along the x direction - :param e: eccentricity of the orbit - :return: Thin disk profile value at (x, y) - """ - - x_, y_ = self.eccentric_coordinate_system(x, y, e, phi_0, center_x, center_y) - R = np.sqrt((x_) ** 2 + (y_)**2) - #rescaling of the caracteristique scales - R_0 = R_0 / (1+e) - R_in = R_in / (1+e) - - profile = np.zeros_like(R) - if R_in <= 0: - profile = amp / (np.exp((R / R_0) ** (0.75)) - 1) - else: - profile = np.where(R <= R_in, profile, - amp / (np.exp((R / R_0) ** (0.75) * (1 - np.sqrt(R_in / R)) ** (-0.25)) - 1)) - - return profile - - def eccentric_coordinate_system(self, x, y, e, phi_0, center_x, center_y): - """ - maps the coordinates x, y into an eccentric coordinate system - - :param x: x-coordinate - :param y: y-coordinate - :param e: eccentricity - :param center_x: center of distortion - :param center_y: center of distortion - :return: distorted coordinates x', y' - """ - R = np.sqrt((x- center_x)**2 + (y-center_y)**2) - phi = np.arctan2((y-center_y), (x-center_x)) - a = (R / (1 + e)) * (1 - e * np.cos(phi - phi_0)) - - xt1 = a * np.cos(phi) - xt2 = a * np.sin(phi) - return xt1, xt2 \ No newline at end of file diff --git a/lenstronomy/Plots/model_band_plot.py b/lenstronomy/Plots/model_band_plot.py index 4f69284c8..33111334b 100644 --- a/lenstronomy/Plots/model_band_plot.py +++ b/lenstronomy/Plots/model_band_plot.py @@ -711,7 +711,6 @@ def source_plot( cax = divider.append_axes("right", size="5%", pad=0.05) cb = plt.colorbar(im, cax=cax) cb.set_label(colorbar_label, fontsize=font_size) - plot_util.scale_bar(ax, d_s, dist=0.1, text='0.1"', color='w', flipped=False, font_size=font_size) if with_caustics is True: ra_caustic_list, dec_caustic_list = self._caustics() diff --git a/lenstronomy/SimulationAPI/ObservationConfig/JWST.py b/lenstronomy/SimulationAPI/ObservationConfig/JWST.py index a4a2fb8a0..e4631e6cb 100644 --- a/lenstronomy/SimulationAPI/ObservationConfig/JWST.py +++ b/lenstronomy/SimulationAPI/ObservationConfig/JWST.py @@ -6,44 +6,42 @@ import lenstronomy.Util.util as util -__all__ = ['JWST'] - - -NIRCAM_F200W_band_obs = {'exposure_time': 3600., - 'sky_brightness': 29.52, #this is derived using the ETC - 'magnitude_zero_point': 28.00, - #'detector': 'NRCA1', - 'num_exposures': 1, - 'seeing': None, - 'psf_type': 'PIXEL', - } - - -NIRCAM_F356W_band_obs = {'exposure_time': 3600., - 'sky_brightness': 28.39, #this is derived using the ETC - 'magnitude_zero_point': 26.47, - #'detector': 'NRCALONG', - 'num_exposures': 1, - 'seeing': None, - 'psf_type': 'PIXEL' # note kernel_point_source (the PSF map) must be provided separately - } - -""" -:keyword exposure_time: exposure time per image (in seconds) -:keyword sky_brightness: sky brightness (in magnitude per square arcseconds in units of electrons) -:keyword magnitude_zero_point: magnitude in which 1 count (e-) per second per arcsecond square is registered -:keyword num_exposures: number of exposures that are combined (depends on coadd_years) -:keyword seeing: Full-Width-at-Half-Maximum (FWHM) of PSF -:keyword psf_type: string, type of PSF ('GAUSSIAN' and 'PIXEL' supported) -""" +__all__ = ["JWST"] + + +NIRCAM_F200W_band_obs = { + "exposure_time": 3600.0, + "sky_brightness": 29.52, # this is derived using the ETC + "magnitude_zero_point": 28.00, + #'detector': 'NRCA1', + "num_exposures": 1, + "seeing": None, + "psf_type": "PIXEL", +} + + +NIRCAM_F356W_band_obs = { + "exposure_time": 3600.0, + "sky_brightness": 28.39, # this is derived using the ETC + "magnitude_zero_point": 26.47, + #'detector': 'NRCALONG', + "num_exposures": 1, + "seeing": None, + "psf_type": "PIXEL", # note kernel_point_source (the PSF map) must be provided separately +} + +# - keyword exposure_time: exposure time per image (in seconds) +# - keyword sky_brightness: sky brightness (in magnitude per square arcseconds in units of electrons) +# - keyword magnitude_zero_point: magnitude in which 1 count (e-) per second per arcsecond square is registered +# - keyword num_exposures: number of exposures that are combined (depends on coadd_years) +# - keyword seeing: Full-Width-at-Half-Maximum (FWHM) of PSF +# - keyword psf_type: string, type of PSF ('GAUSSIAN' and 'PIXEL' supported) class JWST(object): - """ - class contains JWST instrument and observation configurations - """ + """Class contains JWST instrument and observation configurations.""" - def __init__(self, band='F200W', psf_type='PIXEL', coadd_years=None): + def __init__(self, band="F200W", psf_type="PIXEL", coadd_years=None): """ :param band: string, 'F200W' or 'F356W' supported. Determines obs dictionary. @@ -51,41 +49,45 @@ def __init__(self, band='F200W', psf_type='PIXEL', coadd_years=None): :param coadd_years: int, number of years corresponding to num_exposures in obs dict. Currently supported: None. """ - if band == 'F200W': + if band == "F200W": self.obs = NIRCAM_F200W_band_obs - self.arm = 'short' - elif band == 'F356W': + self.arm = "short" + elif band == "F356W": self.obs = NIRCAM_F356W_band_obs - self.arm = 'long' + self.arm = "long" else: - raise ValueError("band %s not supported!"% band) + raise ValueError("band %s not supported!" % band) - if psf_type == 'GAUSSIAN': - self.obs['psf_type'] = 'GAUSSIAN' - elif psf_type != 'PIXEL': + if psf_type == "GAUSSIAN": + self.obs["psf_type"] = "GAUSSIAN" + elif psf_type != "PIXEL": raise ValueError("psf_type %s not supported!" % psf_type) if coadd_years is not None: - raise ValueError(" %s coadd_years not supported! " - "You may manually adjust num_exposures in obs dict if required." % coadd_years) + raise ValueError( + " %s coadd_years not supported! " + "You may manually adjust num_exposures in obs dict if required." + % coadd_years + ) # NIRCAM camera settings - if self.arm == 'short': - self.camera = {'read_noise': 15.77, - 'pixel_scale': 0.031, - 'ccd_gain': 2.05, - } - elif self.arm == 'long': - self.camera = {'read_noise': 13.25, - 'pixel_scale': 0.063, - 'ccd_gain': 1.82, - } - """ - :keyword read_noise: std of noise generated by read-out (in units of electrons) - :keyword pixel_scale: scale (in arcseconds) of pixels - :keyword ccd_gain: electrons/ADU (analog-to-digital unit). - """ + if self.arm == "short": + self.camera = { + "read_noise": 15.77, + "pixel_scale": 0.031, + "ccd_gain": 2.05, + } + elif self.arm == "long": + self.camera = { + "read_noise": 13.25, + "pixel_scale": 0.063, + "ccd_gain": 1.82, + } + + # - keyword read_noise: std of noise generated by read-out (in units of electrons) + # - keyword pixel_scale: scale (in arcseconds) of pixels + # - keyword ccd_gain: electrons/ADU (analog-to-digital unit). def kwargs_single_band(self): """ @@ -93,4 +95,4 @@ def kwargs_single_band(self): :return: merged kwargs from camera and obs dicts """ kwargs = util.merge_dicts(self.camera, self.obs) - return kwargs + return kwargs \ No newline at end of file diff --git a/test/test_SimulationAPI/test_ObservationConfig/test_JWST.py b/test/test_SimulationAPI/test_ObservationConfig/test_JWST.py index 249557809..810a095e6 100644 --- a/test/test_SimulationAPI/test_ObservationConfig/test_JWST.py +++ b/test/test_SimulationAPI/test_ObservationConfig/test_JWST.py @@ -5,11 +5,10 @@ class TestJWST(unittest.TestCase): - def setUp(self): self.F200W = JWST() # default is F200W - self.F356W = JWST(band='F356W') - self.F356W2 = JWST(band='F356W', psf_type='GAUSSIAN') + self.F356W = JWST(band="F356W") + self.F356W2 = JWST(band="F356W", psf_type="GAUSSIAN") kwargs_F200W = self.F200W.kwargs_single_band() kwargs_F356W = self.F356W.kwargs_single_band() @@ -20,29 +19,33 @@ def setUp(self): self.F356W2_band = SingleBand(**kwargs_F356W2) # dictionaries mapping JWST kwargs to SingleBand kwargs - self.camera_settings = {'read_noise': '_read_noise', - 'pixel_scale': 'pixel_scale', - 'ccd_gain': 'ccd_gain'} - self.obs_settings = {'exposure_time': '_exposure_time', - 'sky_brightness': '_sky_brightness_', - 'magnitude_zero_point': '_magnitude_zero_point', - 'num_exposures': '_num_exposures', - 'seeing': '_seeing', - 'psf_type': '_psf_type'} + self.camera_settings = { + "read_noise": "_read_noise", + "pixel_scale": "pixel_scale", + "ccd_gain": "ccd_gain", + } + self.obs_settings = { + "exposure_time": "_exposure_time", + "sky_brightness": "_sky_brightness_", + "magnitude_zero_point": "_magnitude_zero_point", + "num_exposures": "_num_exposures", + "seeing": "_seeing", + "psf_type": "_psf_type", + } self.instrument = Instrument(**self.F200W.camera) def test_JWST_class(self): default = self.F200W - explicit_F200W = JWST(band='F200W', psf_type='PIXEL') + explicit_F200W = JWST(band="F200W", psf_type="PIXEL") self.assertEqual(explicit_F200W.camera, default.camera) self.assertEqual(explicit_F200W.obs, default.obs) with self.assertRaises(ValueError): - bad_band = JWST(band='g') + bad_band = JWST(band="g") with self.assertRaises(ValueError): - bad_psf = JWST(psf_type='blah') + bad_psf = JWST(psf_type="blah") with self.assertRaises(ValueError): bad_coadd_years = JWST(coadd_years=100) @@ -50,18 +53,35 @@ def test_JWST_class(self): def test_JWST_camera(self): # comparing camera settings in JWST instance with those in Instrument instance for config, setting in self.camera_settings.items(): - self.assertEqual(self.F200W.camera[config], getattr(self.instrument, setting), msg=f"{config} did not match") + self.assertEqual( + self.F200W.camera[config], + getattr(self.instrument, setting), + msg=f"{config} did not match", + ) def test_JWST_obs(self): # comparing obs settings in JWST instance with those in SingleBand instance for config, setting in self.obs_settings.items(): - self.assertEqual(self.F200W.obs[config], getattr(self.F200W_band, setting), msg=f"{config} did not match") - self.assertEqual(self.F356W.obs[config], getattr(self.F356W_band, setting), msg=f"{config} did not match") - self.assertEqual(self.F356W2.obs[config], getattr(self.F356W2_band, setting), msg=f"{config} did not match") + self.assertEqual( + self.F200W.obs[config], + getattr(self.F200W_band, setting), + msg=f"{config} did not match", + ) + self.assertEqual( + self.F356W.obs[config], + getattr(self.F356W_band, setting), + msg=f"{config} did not match", + ) + self.assertEqual( + self.F356W2.obs[config], + getattr(self.F356W2_band, setting), + msg=f"{config} did not match", + ) def test_kwargs_single_band(self): kwargs_F200W = util.merge_dicts(self.F200W.camera, self.F200W.obs) self.assertEqual(self.F200W.kwargs_single_band(), kwargs_F200W) -if __name__ == '__main__': - unittest.main() + +if __name__ == "__main__": + unittest.main() \ No newline at end of file From 99ea7554b129c715059b7a073132162b78676995 Mon Sep 17 00:00:00 2001 From: martin-millon Date: Wed, 29 May 2024 12:18:00 -0700 Subject: [PATCH 12/28] remove specifics from my fork --- lenstronomy/Conf/conf_default.yaml | 1 - 1 file changed, 1 deletion(-) diff --git a/lenstronomy/Conf/conf_default.yaml b/lenstronomy/Conf/conf_default.yaml index 7359015d9..64b481d37 100644 --- a/lenstronomy/Conf/conf_default.yaml +++ b/lenstronomy/Conf/conf_default.yaml @@ -14,4 +14,3 @@ conventions: sersic_major_axis: False # if True, defines the half-light radius of the Sersic light profile along the semi-major axis (which is the Galfit convention) # if False, uses the product average of semi-major and semi-minor axis as the convention (default definition for all light profiles in lenstronomy other than the Sersic profile) - From c55f5d45c68313c3a49454c4eea913ba43ca401d Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 29 May 2024 19:36:46 +0000 Subject: [PATCH 13/28] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../SimulationAPI/ObservationConfig/JWST.py | 2 +- lenstronomy/Workflow/psf_fitting.py | 201 ++++++++++-------- .../test_ObservationConfig/test_JWST.py | 2 +- test/test_Workflow/test_psf_fitting.py | 12 +- 4 files changed, 125 insertions(+), 92 deletions(-) diff --git a/lenstronomy/SimulationAPI/ObservationConfig/JWST.py b/lenstronomy/SimulationAPI/ObservationConfig/JWST.py index e4631e6cb..64f6e8a00 100644 --- a/lenstronomy/SimulationAPI/ObservationConfig/JWST.py +++ b/lenstronomy/SimulationAPI/ObservationConfig/JWST.py @@ -95,4 +95,4 @@ def kwargs_single_band(self): :return: merged kwargs from camera and obs dicts """ kwargs = util.merge_dicts(self.camera, self.obs) - return kwargs \ No newline at end of file + return kwargs diff --git a/lenstronomy/Workflow/psf_fitting.py b/lenstronomy/Workflow/psf_fitting.py index 29058b1c2..b77335f20 100644 --- a/lenstronomy/Workflow/psf_fitting.py +++ b/lenstronomy/Workflow/psf_fitting.py @@ -77,14 +77,14 @@ def calc_cornermask(kernelsize, psf_symmetry): return mask def update_iterative( - self, - kwargs_psf, - kwargs_params, - num_iter=10, - keep_psf_error_map=True, - no_break=True, - verbose=True, - **kwargs_psf_update + self, + kwargs_psf, + kwargs_params, + num_iter=10, + keep_psf_error_map=True, + no_break=True, + verbose=True, + **kwargs_psf_update ): """ @@ -161,20 +161,20 @@ def update_iterative( return kwargs_psf_final def update_psf( - self, - kwargs_psf, - kwargs_params, - corner_mask=None, - stacking_method="median", - psf_symmetry=1, - psf_iter_factor=0.2, - block_center_neighbour=0, - error_map_radius=None, - block_center_neighbour_error_map=None, - new_procedure=True, - corner_symmetry=None, - use_starred=False, - kwargs_starred=None, + self, + kwargs_psf, + kwargs_params, + corner_mask=None, + stacking_method="median", + psf_symmetry=1, + psf_iter_factor=0.2, + block_center_neighbour=0, + error_map_radius=None, + block_center_neighbour_error_map=None, + new_procedure=True, + corner_symmetry=None, + use_starred=False, + kwargs_starred=None, ): """ @@ -353,12 +353,16 @@ def update_psf( else: # use STARRED empirical PSF reconstruction routines try: - from starred.procedures.psf_routines import run_multi_steps_PSF_reconstruction + from starred.procedures.psf_routines import ( + run_multi_steps_PSF_reconstruction, + ) from starred.psf.psf import PSF as StarredPSF from starred.psf.parameters import ParametersPSF except ImportError: - raise ImportError("STARRED is not installed. Please install STARRED to use this feature. " - "STARRED available by typing 'pip install starred-astro'.") + raise ImportError( + "STARRED is not installed. Please install STARRED to use this feature. " + "STARRED available by typing 'pip install starred-astro'." + ) image_single_point_source_list = self.image_single_point_source( self._image_model_class, kwargs_params @@ -379,33 +383,55 @@ def update_psf( cutout_size=kernel_size, ) - #STARRED does not like 0 or negative value in the noise maps... so we replace them with the median + # STARRED does not like 0 or negative value in the noise maps... so we replace them with the median sigma2_maps_list = np.asarray(sigma2_maps_list) psf_kernel_list = np.asarray(psf_kernel_list) - #normalise the data + # normalise the data norm = psf_kernel_list[0].max() psf_kernel_list /= norm - sigma2_maps_list /= norm ** 2 + sigma2_maps_list /= norm**2 indneg = np.where(sigma2_maps_list <= 0) if len(indneg[0]) > 0: - warnings.warn("Negative or zero values in the noise maps. Replacing these pixels with the median value.") + warnings.warn( + "Negative or zero values in the noise maps. Replacing these pixels with the median value." + ) sigma2_maps_list[indneg] = np.median(sigma2_maps_list) # ideally the narrow PSF kernel should be saved for the next iteration # todo: masks are not used here, but could be implemented in the future N, image_size, _ = np.shape(psf_kernel_list) - #setup the STARRED model - model = StarredPSF(image_size=image_size, number_of_sources=N, - upsampling_factor=point_source_supersampling_factor, elliptical_moffat=True) + # setup the STARRED model + model = StarredPSF( + image_size=image_size, + number_of_sources=N, + upsampling_factor=point_source_supersampling_factor, + elliptical_moffat=True, + ) - kwargs_init, kwargs_fixed, kwargs_up, kwargs_down = model.smart_guess(psf_kernel_list, fixed_background=False) - parameters = ParametersPSF(kwargs_init, kwargs_fixed, kwargs_up=kwargs_up, kwargs_down=kwargs_down) + kwargs_init, kwargs_fixed, kwargs_up, kwargs_down = model.smart_guess( + psf_kernel_list, fixed_background=False + ) + parameters = ParametersPSF( + kwargs_init, kwargs_fixed, kwargs_up=kwargs_up, kwargs_down=kwargs_down + ) - model, parameters, loss, kwargs_partial_list, LogL_list, loss_history_list = run_multi_steps_PSF_reconstruction( - psf_kernel_list, model, - parameters, sigma2_maps_list, masks=None, **kwargs_starred) + ( + model, + parameters, + loss, + kwargs_partial_list, + LogL_list, + loss_history_list, + ) = run_multi_steps_PSF_reconstruction( + psf_kernel_list, + model, + parameters, + sigma2_maps_list, + masks=None, + **kwargs_starred + ) # kernel_new, narrow_psf_new, psf_kernel_list, _ = starred.procedures.psf_routines.update_PSF(kernel_old, # psf_kernel_list, @@ -414,10 +440,16 @@ def update_psf( # subsampling_factor=point_source_supersampling_factor, # **kwargs_starred) - psf_kernel_list = model.model(**kwargs_partial_list[-1], high_res=False) #return model for each point source at the resolution of the data + psf_kernel_list = model.model( + **kwargs_partial_list[-1], high_res=False + ) # return model for each point source at the resolution of the data psf_kernel_list = [psf_k / np.sum(psf_k) for psf_k in psf_kernel_list] - kernel_new = model.get_full_psf(**kwargs_partial_list[-1], norm=True, high_res=True) #subsampled PSF - kernel_new = psf_iter_factor * kernel_new + (1.0 - psf_iter_factor) * kernel_old + kernel_new = model.get_full_psf( + **kwargs_partial_list[-1], norm=True, high_res=True + ) # subsampled PSF + kernel_new = ( + psf_iter_factor * kernel_new + (1.0 - psf_iter_factor) * kernel_old + ) error_map = self.error_map_estimate_new( kernel_new, @@ -435,14 +467,13 @@ def update_psf( # ax[2].imshow(kernel_new- kernel_old) # plt.show() - kwargs_psf_new["kernel_point_source"] = kernel_new self._image_model_class.update_psf(PSF(**kwargs_psf_new)) logL_after, _ = self._image_model_class.likelihood_data_given_model( **kwargs_params ) if use_starred: - logL_after =0. + logL_after = 0.0 return kwargs_psf_new, logL_after, error_map def image_single_point_source(self, image_model_class, kwargs_params): @@ -493,15 +524,15 @@ def _point_sources_list(image_model_class, kwargs_ps, kwargs_lens, k=None): return point_list def cutout_psf( - self, - ra_image, - dec_image, - x, - y, - image_list, - kernel_size, - kernel_init, - block_center_neighbour=0, + self, + ra_image, + dec_image, + x, + y, + image_list, + kernel_size, + kernel_init, + block_center_neighbour=0, ): """ @@ -536,15 +567,15 @@ def cutout_psf( return kernel_list def psf_estimate_individual( - self, - ra_image, - dec_image, - point_amp, - residuals, - cutout_size, - kernel_guess, - supersampling_factor, - block_center_neighbour, + self, + ra_image, + dec_image, + point_amp, + residuals, + cutout_size, + kernel_guess, + supersampling_factor, + block_center_neighbour, ): """ @@ -693,13 +724,13 @@ def cutout_psf_single(x, y, image, mask, kernel_size, kernel_init): @staticmethod ####Correct this based on Maverick Oh's note about lack of flux conservation. def combine_psf( - kernel_list_new, - kernel_old, - factor=1.0, - stacking_option="median", - symmetry=1, - corner_symmetry=None, - corner_mask=None, + kernel_list_new, + kernel_old, + factor=1.0, + stacking_option="median", + symmetry=1, + corner_symmetry=None, + corner_mask=None, ): ## TODO: Account for image edges/masked pixels. """Updates psf estimate based on old kernel and several new estimates. @@ -790,14 +821,14 @@ def combine_psf( return kernel_return def error_map_estimate_new( - self, - psf_kernel, - psf_kernel_list, - ra_image, - dec_image, - point_amp, - supersampling_factor, - error_map_radius=None, + self, + psf_kernel, + psf_kernel_list, + ra_image, + dec_image, + point_amp, + supersampling_factor, + error_map_radius=None, ): """Relative uncertainty in the psf model (in quadrature) per pixel based on residuals achieved in the image. @@ -837,7 +868,7 @@ def error_map_estimate_new( residuals_i = np.abs(kernel_low - kernel_low_i) residuals_i -= np.sqrt(C_D_cutout) / amp_i residuals_i[residuals_i < 0] = 0 - error_map_list[i, :, :] = residuals_i ** 2 + error_map_list[i, :, :] = residuals_i**2 error_map = np.median(error_map_list, axis=0) error_map[kernel_low > 0] /= kernel_low[kernel_low > 0] ** 2 @@ -855,14 +886,14 @@ def error_map_estimate_new( return error_map def error_map_estimate( - self, - kernel, - star_cutout_list, - amp, - x_pos, - y_pos, - error_map_radius=None, - block_center_neighbour=0, + self, + kernel, + star_cutout_list, + amp, + x_pos, + y_pos, + error_map_radius=None, + block_center_neighbour=0, ): """Provides a psf_error_map based on the goodness of fit of the given PSF kernel on the point source cutouts, their estimated amplitudes and positions. @@ -918,7 +949,7 @@ def error_map_estimate( residual[residual < 0] = 0 # estimate relative error per star residual /= amp_i - error_map_list[i, :, :] = residual ** 2 * mask_i + error_map_list[i, :, :] = residual**2 * mask_i mask_list[i, :, :] = mask_i # take median absolute error for each pixel # TODO: only for pixels that are not masked diff --git a/test/test_SimulationAPI/test_ObservationConfig/test_JWST.py b/test/test_SimulationAPI/test_ObservationConfig/test_JWST.py index 810a095e6..2ed0dff0c 100644 --- a/test/test_SimulationAPI/test_ObservationConfig/test_JWST.py +++ b/test/test_SimulationAPI/test_ObservationConfig/test_JWST.py @@ -84,4 +84,4 @@ def test_kwargs_single_band(self): if __name__ == "__main__": - unittest.main() \ No newline at end of file + unittest.main() diff --git a/test/test_Workflow/test_psf_fitting.py b/test/test_Workflow/test_psf_fitting.py index 03aeb2428..ca2dc11d8 100644 --- a/test/test_Workflow/test_psf_fitting.py +++ b/test/test_Workflow/test_psf_fitting.py @@ -179,18 +179,20 @@ def test_update_psf(self): diff_new = np.sum((kernel_new - kernel_true) ** 2) assert diff_old > diff_new - #test STARRED + # test STARRED kwargs_psf_iter_starred = { "stacking_method": "median", "error_map_radius": 0.5, "psf_iter_factor": 0.2, "new_procedure": False, "use_starred": True, - "kwargs_starred": {'verbose':False, 'lambda_scales':3, 'lambda_hf':3}, + "kwargs_starred": {"verbose": False, "lambda_scales": 3, "lambda_hf": 3}, } - kwargs_psf_return_starred, improved_bool_starred, error_map_starred = self.psf_fitting.update_psf( - kwargs_psf, self.kwargs_params, **kwargs_psf_iter_starred + kwargs_psf_return_starred, improved_bool_starred, error_map_starred = ( + self.psf_fitting.update_psf( + kwargs_psf, self.kwargs_params, **kwargs_psf_iter_starred + ) ) assert improved_bool_starred # fig, ax = plt.subplots(2,3) @@ -210,7 +212,7 @@ def test_update_psf(self): # plt.show() kernel_new_starred = kwargs_psf_return_starred["kernel_point_source"] - diff_new_starred = np.sum((kernel_new_starred- kernel_true) ** 2) + diff_new_starred = np.sum((kernel_new_starred - kernel_true) ** 2) # print(diff_new_starred, diff_new, diff_old) assert diff_old > diff_new_starred From 0cdc4f28437eb53fd48c7a31dd3202cb801a4f75 Mon Sep 17 00:00:00 2001 From: martin-millon Date: Fri, 31 May 2024 11:00:37 -0700 Subject: [PATCH 14/28] Final Starred implementation --- lenstronomy/Util/kernel_util.py | 2 +- lenstronomy/Workflow/psf_fitting.py | 37 +++++++++----------------- test/test_Workflow/test_psf_fitting.py | 16 ----------- 3 files changed, 13 insertions(+), 42 deletions(-) diff --git a/lenstronomy/Util/kernel_util.py b/lenstronomy/Util/kernel_util.py index 57221e80d..c2c55bccc 100644 --- a/lenstronomy/Util/kernel_util.py +++ b/lenstronomy/Util/kernel_util.py @@ -403,7 +403,7 @@ def degrade_kernel(kernel_super, degrading_factor): def averaging_odd_kernel(kernel_super, degrading_factor): """""" n_kernel = len(kernel_super) - numPix = int(round(n_kernel / degrading_factor + 0.5)) + numPix = int(round(n_kernel / degrading_factor)) if numPix % 2 == 0: numPix += 1 n_high = numPix * degrading_factor diff --git a/lenstronomy/Workflow/psf_fitting.py b/lenstronomy/Workflow/psf_fitting.py index b77335f20..fdfeb51e1 100644 --- a/lenstronomy/Workflow/psf_fitting.py +++ b/lenstronomy/Workflow/psf_fitting.py @@ -133,7 +133,6 @@ def update_iterative( kwargs_psf_new, _kwargs_params, corner_mask, **kwargs_psf_update ) - print(logL_best, logL_after) if logL_after > logL_best: kwargs_psf_final = copy.deepcopy(kwargs_psf_new) error_map_final = copy.deepcopy(error_map) @@ -363,11 +362,15 @@ def update_psf( "STARRED is not installed. Please install STARRED to use this feature. " "STARRED available by typing 'pip install starred-astro'." ) + if point_source_supersampling_factor != 3: + warnings.warn("Point source subsampling factor of 3 is highly recommended when using Starred PSF iteration routine") + kernel_old_high_res = psf_class.kernel_point_source_supersampled( + supersampling_factor=point_source_supersampling_factor + ) image_single_point_source_list = self.image_single_point_source( self._image_model_class, kwargs_params ) - # get the non-aligned cutouts psf_kernel_list = self.point_like_source_cutouts( x_pos=x_, @@ -416,7 +419,6 @@ def update_psf( parameters = ParametersPSF( kwargs_init, kwargs_fixed, kwargs_up=kwargs_up, kwargs_down=kwargs_down ) - ( model, parameters, @@ -433,26 +435,19 @@ def update_psf( **kwargs_starred ) - # kernel_new, narrow_psf_new, psf_kernel_list, _ = starred.procedures.psf_routines.update_PSF(kernel_old, - # psf_kernel_list, - # sigma2_maps_list, - # masks=None, - # subsampling_factor=point_source_supersampling_factor, - # **kwargs_starred) - psf_kernel_list = model.model( - **kwargs_partial_list[-1], high_res=False - ) # return model for each point source at the resolution of the data - psf_kernel_list = [psf_k / np.sum(psf_k) for psf_k in psf_kernel_list] - kernel_new = model.get_full_psf( + **kwargs_partial_list[-1], high_res=True, + ) # return model for each point source at the super-sampled resolution + psf_kernel_list = [psf_k / np.sum(psf_k) for psf_k in psf_kernel_list] #normalise the models + kernel_new_high_res = model.get_full_psf( **kwargs_partial_list[-1], norm=True, high_res=True ) # subsampled PSF kernel_new = ( - psf_iter_factor * kernel_new + (1.0 - psf_iter_factor) * kernel_old + psf_iter_factor * kernel_new_high_res + (1.0 - psf_iter_factor) * kernel_old_high_res ) error_map = self.error_map_estimate_new( - kernel_new, + kernel_new_high_res, psf_kernel_list, ra_image, dec_image, @@ -461,19 +456,11 @@ def update_psf( error_map_radius=error_map_radius, ) - # fig, ax = plt.subplots(1,3, figsize=(15,12)) - # ax[0].imshow(np.log10(kernel_new)) - # ax[1].imshow(np.log10(kernel_old)) - # ax[2].imshow(kernel_new- kernel_old) - # plt.show() - kwargs_psf_new["kernel_point_source"] = kernel_new self._image_model_class.update_psf(PSF(**kwargs_psf_new)) logL_after, _ = self._image_model_class.likelihood_data_given_model( **kwargs_params ) - if use_starred: - logL_after = 0.0 return kwargs_psf_new, logL_after, error_map def image_single_point_source(self, image_model_class, kwargs_params): @@ -834,7 +821,7 @@ def error_map_estimate_new( residuals achieved in the image. :param psf_kernel: PSF kernel (super-sampled) - :param psf_kernel_list: list of individual best PSF kernel estimates + :param psf_kernel_list: list of individual best PSF kernel estimates (super-sampled) :param ra_image: image positions in angles :param dec_image: image positions in angles :param point_amp: image amplitude diff --git a/test/test_Workflow/test_psf_fitting.py b/test/test_Workflow/test_psf_fitting.py index ca2dc11d8..a0fa3b5dc 100644 --- a/test/test_Workflow/test_psf_fitting.py +++ b/test/test_Workflow/test_psf_fitting.py @@ -195,22 +195,6 @@ def test_update_psf(self): ) ) assert improved_bool_starred - # fig, ax = plt.subplots(2,3) - # ax[0,0].imshow(kwargs_psf_return_starred["kernel_point_source"]) - # ax[0,1].imshow(kernel_old) - # ax[0,2].imshow(kwargs_psf_return_starred["kernel_point_source"]- kernel_old) - # - # ax[1,0].imshow(kwargs_psf_return_starred["kernel_point_source"]) - # ax[1,1].imshow(kernel_true) - # ax[1,2].imshow(kwargs_psf_return_starred["kernel_point_source"] - kernel_true) - # plt.show() - # - # fig, ax = plt.subplots(1,3) - # ax[0].imshow(kwargs_psf_return_starred["kernel_point_source"]) - # ax[1].imshow(error_map_starred) - # ax[2].imshow(error_map) - # plt.show() - kernel_new_starred = kwargs_psf_return_starred["kernel_point_source"] diff_new_starred = np.sum((kernel_new_starred - kernel_true) ** 2) From e64707667ef9dfe986bed4f510a08297418395ef Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 31 May 2024 18:01:12 +0000 Subject: [PATCH 15/28] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- lenstronomy/Workflow/psf_fitting.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/lenstronomy/Workflow/psf_fitting.py b/lenstronomy/Workflow/psf_fitting.py index fdfeb51e1..8cf1d4f15 100644 --- a/lenstronomy/Workflow/psf_fitting.py +++ b/lenstronomy/Workflow/psf_fitting.py @@ -363,7 +363,9 @@ def update_psf( "STARRED available by typing 'pip install starred-astro'." ) if point_source_supersampling_factor != 3: - warnings.warn("Point source subsampling factor of 3 is highly recommended when using Starred PSF iteration routine") + warnings.warn( + "Point source subsampling factor of 3 is highly recommended when using Starred PSF iteration routine" + ) kernel_old_high_res = psf_class.kernel_point_source_supersampled( supersampling_factor=point_source_supersampling_factor @@ -432,18 +434,22 @@ def update_psf( parameters, sigma2_maps_list, masks=None, - **kwargs_starred + **kwargs_starred, ) psf_kernel_list = model.model( - **kwargs_partial_list[-1], high_res=True, + **kwargs_partial_list[-1], + high_res=True, ) # return model for each point source at the super-sampled resolution - psf_kernel_list = [psf_k / np.sum(psf_k) for psf_k in psf_kernel_list] #normalise the models + psf_kernel_list = [ + psf_k / np.sum(psf_k) for psf_k in psf_kernel_list + ] # normalise the models kernel_new_high_res = model.get_full_psf( **kwargs_partial_list[-1], norm=True, high_res=True ) # subsampled PSF kernel_new = ( - psf_iter_factor * kernel_new_high_res + (1.0 - psf_iter_factor) * kernel_old_high_res + psf_iter_factor * kernel_new_high_res + + (1.0 - psf_iter_factor) * kernel_old_high_res ) error_map = self.error_map_estimate_new( @@ -821,7 +827,8 @@ def error_map_estimate_new( residuals achieved in the image. :param psf_kernel: PSF kernel (super-sampled) - :param psf_kernel_list: list of individual best PSF kernel estimates (super-sampled) + :param psf_kernel_list: list of individual best PSF kernel estimates (super- + sampled) :param ra_image: image positions in angles :param dec_image: image positions in angles :param point_amp: image amplitude From 09a2ac841bc6b93059d615de36036a0b39c555f5 Mon Sep 17 00:00:00 2001 From: martin-millon Date: Fri, 31 May 2024 11:29:42 -0700 Subject: [PATCH 16/28] updating test-requirements for automatic tests --- test/test_Workflow/test_psf_fitting.py | 1 - test_requirements.txt | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/test/test_Workflow/test_psf_fitting.py b/test/test_Workflow/test_psf_fitting.py index a0fa3b5dc..558538ecd 100644 --- a/test/test_Workflow/test_psf_fitting.py +++ b/test/test_Workflow/test_psf_fitting.py @@ -1,6 +1,5 @@ __author__ = "sibirrer" -import matplotlib.pyplot as plt import pytest import numpy as np import copy diff --git a/test_requirements.txt b/test_requirements.txt index 08e3c232d..0090a924d 100644 --- a/test_requirements.txt +++ b/test_requirements.txt @@ -2,4 +2,4 @@ colossus git+https://github.com/mattgomer/SKiNN@main#egg=SKiNN git+https://github.com/aymgal/coolest@main#egg=coolest -git+https://gitlab.com/cosmograil/starred.git@main#egg=starred +starred-astro>=1.4.2 From 1376b22f9028b93727fa49ffda87ea351891f9d5 Mon Sep 17 00:00:00 2001 From: martin-millon Date: Fri, 31 May 2024 11:35:02 -0700 Subject: [PATCH 17/28] updating test-requirements for automatic tests --- setup.py | 1 + 1 file changed, 1 insertion(+) diff --git a/setup.py b/setup.py index b28f85846..c5cfd9f43 100644 --- a/setup.py +++ b/setup.py @@ -68,6 +68,7 @@ def run_tests(self): "zeus-mcmc>=2.4.0", "nautilus-sampler>=0.2.1", "coolest", + "starred-astro>=1.4.2", ] PACKAGE_PATH = os.path.abspath(os.path.join(__file__, os.pardir)) From ef4020ec62f989110837490cb13875166863c086 Mon Sep 17 00:00:00 2001 From: martin-millon Date: Fri, 31 May 2024 16:46:15 -0700 Subject: [PATCH 18/28] adding PSF masks warning messages --- lenstronomy/Workflow/psf_fitting.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/lenstronomy/Workflow/psf_fitting.py b/lenstronomy/Workflow/psf_fitting.py index 8cf1d4f15..e3b505973 100644 --- a/lenstronomy/Workflow/psf_fitting.py +++ b/lenstronomy/Workflow/psf_fitting.py @@ -252,7 +252,7 @@ def update_psf( kernel_old = psf_class.kernel_point_source kernel_size = len(kernel_old) - # TODO: refractor this if new_procedure is adopted + # TODO: refractor this if new_procedure is definitively adopted if not new_procedure and not use_starred: image_single_point_source_list = self.image_single_point_source( self._image_model_class, kwargs_params @@ -366,6 +366,8 @@ def update_psf( warnings.warn( "Point source subsampling factor of 3 is highly recommended when using Starred PSF iteration routine" ) + if psf_symmetry != 1: + warnings.warn("Starred PSF fitting routine does not assume any PSF symmetry. Setting psf_symmetry=1.") kernel_old_high_res = psf_class.kernel_point_source_supersampled( supersampling_factor=point_source_supersampling_factor @@ -403,10 +405,17 @@ def update_psf( ) sigma2_maps_list[indneg] = np.median(sigma2_maps_list) - # ideally the narrow PSF kernel should be saved for the next iteration - # todo: masks are not used here, but could be implemented in the future N, image_size, _ = np.shape(psf_kernel_list) + # todo: Lenstronomy is not supporting custom PSF mask for the moment and propagating the corner mask makes no sense for STARRED, which is not using psf_symetry. + if corner_mask is not None: + warnings.warn("Corner mask is not used in Starred PSF fitting routine. Corner_symmetry is ignored.") + #possible implementation of corner_mask in STARRED could look like if custom masks can be propagated one day. + #starred requires a mask at the original resolution and has opposite convention (0 is masked) + # corner_mask_starred = kernel_util.degrade_kernel(np.invert(corner_mask), point_source_supersampling_factor) + # corner_mask_starred = np.repeat(corner_mask_starred[np.newaxis, :, :], N, axis=0) + + # setup the STARRED model model = StarredPSF( image_size=image_size, From 5d8661c58ba4747b9da17cd6785c89d402c264a0 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 31 May 2024 23:47:33 +0000 Subject: [PATCH 19/28] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- lenstronomy/Workflow/psf_fitting.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/lenstronomy/Workflow/psf_fitting.py b/lenstronomy/Workflow/psf_fitting.py index e3b505973..b767325d8 100644 --- a/lenstronomy/Workflow/psf_fitting.py +++ b/lenstronomy/Workflow/psf_fitting.py @@ -367,7 +367,9 @@ def update_psf( "Point source subsampling factor of 3 is highly recommended when using Starred PSF iteration routine" ) if psf_symmetry != 1: - warnings.warn("Starred PSF fitting routine does not assume any PSF symmetry. Setting psf_symmetry=1.") + warnings.warn( + "Starred PSF fitting routine does not assume any PSF symmetry. Setting psf_symmetry=1." + ) kernel_old_high_res = psf_class.kernel_point_source_supersampled( supersampling_factor=point_source_supersampling_factor @@ -409,13 +411,14 @@ def update_psf( # todo: Lenstronomy is not supporting custom PSF mask for the moment and propagating the corner mask makes no sense for STARRED, which is not using psf_symetry. if corner_mask is not None: - warnings.warn("Corner mask is not used in Starred PSF fitting routine. Corner_symmetry is ignored.") - #possible implementation of corner_mask in STARRED could look like if custom masks can be propagated one day. - #starred requires a mask at the original resolution and has opposite convention (0 is masked) + warnings.warn( + "Corner mask is not used in Starred PSF fitting routine. Corner_symmetry is ignored." + ) + # possible implementation of corner_mask in STARRED could look like if custom masks can be propagated one day. + # starred requires a mask at the original resolution and has opposite convention (0 is masked) # corner_mask_starred = kernel_util.degrade_kernel(np.invert(corner_mask), point_source_supersampling_factor) # corner_mask_starred = np.repeat(corner_mask_starred[np.newaxis, :, :], N, axis=0) - # setup the STARRED model model = StarredPSF( image_size=image_size, From 8bc19f1d25de48f4b7e7e0bc8758f83f335dee1e Mon Sep 17 00:00:00 2001 From: martin-millon Date: Tue, 4 Jun 2024 12:25:36 -0700 Subject: [PATCH 20/28] Re-running the CI pipeline --- lenstronomy/Workflow/psf_fitting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lenstronomy/Workflow/psf_fitting.py b/lenstronomy/Workflow/psf_fitting.py index e3b505973..ff926c06d 100644 --- a/lenstronomy/Workflow/psf_fitting.py +++ b/lenstronomy/Workflow/psf_fitting.py @@ -393,7 +393,7 @@ def update_psf( # STARRED does not like 0 or negative value in the noise maps... so we replace them with the median sigma2_maps_list = np.asarray(sigma2_maps_list) psf_kernel_list = np.asarray(psf_kernel_list) - # normalise the data + # Starred works best with normalised the data norm = psf_kernel_list[0].max() psf_kernel_list /= norm sigma2_maps_list /= norm**2 From aee2e5b20a80657f023f529da1d908d070384a86 Mon Sep 17 00:00:00 2001 From: martin-millon Date: Tue, 4 Jun 2024 14:56:15 -0700 Subject: [PATCH 21/28] starred installation from repo --- test_requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test_requirements.txt b/test_requirements.txt index 0090a924d..697194d76 100644 --- a/test_requirements.txt +++ b/test_requirements.txt @@ -2,4 +2,4 @@ colossus git+https://github.com/mattgomer/SKiNN@main#egg=SKiNN git+https://github.com/aymgal/coolest@main#egg=coolest -starred-astro>=1.4.2 +git+https://gitlab.com/cosmograil/starred@main#egg=starred-astro From 2d62d4cb8a6f7b89737de2129c2431aec098b1bb Mon Sep 17 00:00:00 2001 From: martin-millon Date: Tue, 4 Jun 2024 16:43:59 -0700 Subject: [PATCH 22/28] starred installation from repo --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index c5cfd9f43..9e0b375c6 100644 --- a/setup.py +++ b/setup.py @@ -68,7 +68,7 @@ def run_tests(self): "zeus-mcmc>=2.4.0", "nautilus-sampler>=0.2.1", "coolest", - "starred-astro>=1.4.2", + "starred-astro>=1.4.3", ] PACKAGE_PATH = os.path.abspath(os.path.join(__file__, os.pardir)) From c18babd637bfca706aabfa9bf5f7c434d8887f2a Mon Sep 17 00:00:00 2001 From: martin-millon Date: Tue, 4 Jun 2024 17:04:40 -0700 Subject: [PATCH 23/28] changed starred requirements --- test/test_Workflow/test_psf_fitting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_Workflow/test_psf_fitting.py b/test/test_Workflow/test_psf_fitting.py index 558538ecd..eb300654b 100644 --- a/test/test_Workflow/test_psf_fitting.py +++ b/test/test_Workflow/test_psf_fitting.py @@ -182,7 +182,7 @@ def test_update_psf(self): kwargs_psf_iter_starred = { "stacking_method": "median", "error_map_radius": 0.5, - "psf_iter_factor": 0.2, + "psf_iter_factor": 1., "new_procedure": False, "use_starred": True, "kwargs_starred": {"verbose": False, "lambda_scales": 3, "lambda_hf": 3}, From 87edb6e0bd483d03ef74219b68987fe1210cecc5 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 5 Jun 2024 00:05:12 +0000 Subject: [PATCH 24/28] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- test/test_Workflow/test_psf_fitting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_Workflow/test_psf_fitting.py b/test/test_Workflow/test_psf_fitting.py index eb300654b..741b94473 100644 --- a/test/test_Workflow/test_psf_fitting.py +++ b/test/test_Workflow/test_psf_fitting.py @@ -182,7 +182,7 @@ def test_update_psf(self): kwargs_psf_iter_starred = { "stacking_method": "median", "error_map_radius": 0.5, - "psf_iter_factor": 1., + "psf_iter_factor": 1.0, "new_procedure": False, "use_starred": True, "kwargs_starred": {"verbose": False, "lambda_scales": 3, "lambda_hf": 3}, From fec6dda2ecb60c02e2cfd8b46191f88abdc094e4 Mon Sep 17 00:00:00 2001 From: martin-millon Date: Wed, 5 Jun 2024 10:55:37 -0700 Subject: [PATCH 25/28] testing from another fork of SKiNN --- test_requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test_requirements.txt b/test_requirements.txt index 697194d76..08d72ea33 100644 --- a/test_requirements.txt +++ b/test_requirements.txt @@ -1,5 +1,5 @@ # requirements for tests colossus -git+https://github.com/mattgomer/SKiNN@main#egg=SKiNN +git+https://github.com/martin-millon/SKiNN@main#egg=SKiNN git+https://github.com/aymgal/coolest@main#egg=coolest git+https://gitlab.com/cosmograil/starred@main#egg=starred-astro From 0b43e24a59ec964de9e25120180f1e706a9dfa2d Mon Sep 17 00:00:00 2001 From: martin-millon Date: Wed, 5 Jun 2024 13:23:50 -0700 Subject: [PATCH 26/28] improving coverage --- test/test_Workflow/test_psf_fitting.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/test_Workflow/test_psf_fitting.py b/test/test_Workflow/test_psf_fitting.py index 741b94473..969e91a5b 100644 --- a/test/test_Workflow/test_psf_fitting.py +++ b/test/test_Workflow/test_psf_fitting.py @@ -179,10 +179,13 @@ def test_update_psf(self): assert diff_old > diff_new # test STARRED + self.psf_fitting._image_model_class.Data._C_D[42,70] = -1 #introducing one negative value in one of the noise maps cutouts to test warning message kwargs_psf_iter_starred = { "stacking_method": "median", "error_map_radius": 0.5, "psf_iter_factor": 1.0, + "psf_symmetry": 2, #to test warning message + "corner_symmetry": 2, #to test warning message "new_procedure": False, "use_starred": True, "kwargs_starred": {"verbose": False, "lambda_scales": 3, "lambda_hf": 3}, From eb90df07c0db8e35450f3ed525c48da7fc596b00 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 5 Jun 2024 20:24:23 +0000 Subject: [PATCH 27/28] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- test/test_Workflow/test_psf_fitting.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/test/test_Workflow/test_psf_fitting.py b/test/test_Workflow/test_psf_fitting.py index 969e91a5b..0a0f9dc1b 100644 --- a/test/test_Workflow/test_psf_fitting.py +++ b/test/test_Workflow/test_psf_fitting.py @@ -179,13 +179,15 @@ def test_update_psf(self): assert diff_old > diff_new # test STARRED - self.psf_fitting._image_model_class.Data._C_D[42,70] = -1 #introducing one negative value in one of the noise maps cutouts to test warning message + self.psf_fitting._image_model_class.Data._C_D[42, 70] = ( + -1 + ) # introducing one negative value in one of the noise maps cutouts to test warning message kwargs_psf_iter_starred = { "stacking_method": "median", "error_map_radius": 0.5, "psf_iter_factor": 1.0, - "psf_symmetry": 2, #to test warning message - "corner_symmetry": 2, #to test warning message + "psf_symmetry": 2, # to test warning message + "corner_symmetry": 2, # to test warning message "new_procedure": False, "use_starred": True, "kwargs_starred": {"verbose": False, "lambda_scales": 3, "lambda_hf": 3}, From 6e07853983f5e3531ade8e09794bc64b4dee9fac Mon Sep 17 00:00:00 2001 From: martin-millon Date: Thu, 6 Jun 2024 14:18:36 -0700 Subject: [PATCH 28/28] switching back to SKiNN main branch --- test_requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test_requirements.txt b/test_requirements.txt index 08d72ea33..697194d76 100644 --- a/test_requirements.txt +++ b/test_requirements.txt @@ -1,5 +1,5 @@ # requirements for tests colossus -git+https://github.com/martin-millon/SKiNN@main#egg=SKiNN +git+https://github.com/mattgomer/SKiNN@main#egg=SKiNN git+https://github.com/aymgal/coolest@main#egg=coolest git+https://gitlab.com/cosmograil/starred@main#egg=starred-astro