Tracer model with partition to different surface brightness models

docs/lenstronomy.Analysis.rst

@@ -4,51 +4,66 @@ lenstronomy.Analysis package
 Submodules
 ----------
 
+lenstronomy.Analysis.image\_reconstruction module
+-------------------------------------------------
+
+.. automodule:: lenstronomy.Analysis.image_reconstruction
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 lenstronomy.Analysis.kinematics\_api module
 -------------------------------------------
 
 .. automodule:: lenstronomy.Analysis.kinematics_api
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
 
 lenstronomy.Analysis.lens\_profile module
 -----------------------------------------
 
 .. automodule:: lenstronomy.Analysis.lens_profile
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
 
 lenstronomy.Analysis.light2mass module
 --------------------------------------
 
 .. automodule:: lenstronomy.Analysis.light2mass
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
 
 lenstronomy.Analysis.light\_profile module
 ------------------------------------------
 
 .. automodule:: lenstronomy.Analysis.light_profile
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+lenstronomy.Analysis.multi\_patch\_reconstruction module
+--------------------------------------------------------
+
+.. automodule:: lenstronomy.Analysis.multi_patch_reconstruction
+   :members:
+   :undoc-members:
+   :show-inheritance:
 
 lenstronomy.Analysis.td\_cosmography module
 -------------------------------------------
 
 .. automodule:: lenstronomy.Analysis.td_cosmography
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
+   :members:
+   :undoc-members:
+   :show-inheritance:
 
 Module contents
 ---------------
 
 .. automodule:: lenstronomy.Analysis
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:

lenstronomy/Analysis/kinematics_api.py

@@ -310,6 +310,12 @@ def kinematic_light_profile(self, kwargs_lens_light, r_eff=None, MGE_fit=False,
             return None, {'r_eff': r_eff}
         light_profile_list = []
         kwargs_light = []
+        if Hernquist_approx is True:
+            if r_eff is None:
+                raise ValueError('r_eff needs to be pre-computed and specified when using the Hernquist approximation')
+            light_profile_list = ['HERNQUIST']
+            kwargs_light = [{'Rs': r_eff * 0.551, 'amp': 1.}]
+            return light_profile_list, kwargs_light
         if model_kinematics_bool is None:
             model_kinematics_bool = [True] * len(kwargs_lens_light)
         for i, light_model in enumerate(self._lens_light_model_list):
@@ -321,19 +327,14 @@ def kinematic_light_profile(self, kwargs_lens_light, r_eff=None, MGE_fit=False,
                     kwargs_lens_light_i['e1'] = 0
                     kwargs_lens_light_i['e2'] = 0
                 kwargs_light.append(kwargs_lens_light_i)
-        if Hernquist_approx is True:
-            if r_eff is None:
-                raise ValueError('r_eff needs to be pre-computed and specified when using the Hernquist approximation')
-            light_profile_list = ['HERNQUIST']
-            kwargs_light = [{'Rs': r_eff * 0.551, 'amp': 1.}]
-        else:
-            if MGE_fit is True:
-                if kwargs_mge is None:
-                    raise ValueError('kwargs_mge must be provided to compute the MGE')
-                amps, sigmas, center_x, center_y = self._lensLightProfile.multi_gaussian_decomposition(
-                    kwargs_lens_light, model_bool_list=model_kinematics_bool, r_h=r_eff, **kwargs_mge)
-                light_profile_list = ['MULTI_GAUSSIAN']
-                kwargs_light = [{'amp': amps, 'sigma': sigmas}]
+
+        if MGE_fit is True:
+            if kwargs_mge is None:
+                raise ValueError('kwargs_mge must be provided to compute the MGE')
+            amps, sigmas, center_x, center_y = self._lensLightProfile.multi_gaussian_decomposition(
+                kwargs_lens_light, model_bool_list=model_kinematics_bool, r_h=r_eff, **kwargs_mge)
+            light_profile_list = ['MULTI_GAUSSIAN']
+            kwargs_light = [{'amp': amps, 'sigma': sigmas}]
         return light_profile_list, kwargs_light
 
     def kinematics_modeling_settings(self, anisotropy_model, kwargs_numerics_galkin, analytic_kinematics=False,

lenstronomy/Cosmo/lens_cosmo.py

@@ -27,15 +27,6 @@ def __init__(self, z_lens, z_source, cosmo=None):
         self.background = Background(cosmo=cosmo)
         self.nfw_param = NFWParam(cosmo=cosmo)
 
-    def a_z(self, z):
-        """
-        convert redshift into scale factor
-
-        :param z: redshift
-        :return: scale factor
-        """
-        return 1. / (1. + z)
-
     @property
     def h(self):
         return self.background.cosmo.H(0).value / 100.

lenstronomy/Cosmo/nfw_param.py

@@ -31,25 +31,29 @@ def rhoc_z(self, z):
         :return: critical density of the universe at redshift z in physical units [h^2 M_sun Mpc^-3]
         """
         return self.rhoc * (self.cosmo.efunc(z)) ** 2
-        #return self.rhoc*(1+z)**3
+        # return self.rhoc*(1+z)**3
 
-    def M200(self, rs, rho0, c):
+    @staticmethod
+    def M200(rs, rho0, c):
         """
-        M(R_200) calculation for NFW profile
+        Calculation of the mass enclosed r_200 for NFW profile defined as
+
+        .. math::
+            M_{200} = 4 \\pi \\rho_0^{3} * \\left(\\log(1+c) - c / (1 + c)  \\right))
 
         :param rs: scale radius
         :type rs: float
-        :param rho0: density normalization (characteristic density)
+        :param rho0: density normalization (characteristic density) in units mass/[distance unit of rs]^3
         :type rho0: float
         :param c: concentration
         :type c: float [4,40]
-        :return: M(R_200) density
+        :return: M(R_200) mass in units of rho0 * rs^3
         """
         return 4 * np.pi * rho0 * rs ** 3 * (np.log(1. + c) - c / (1. + c))
 
     def r200_M(self, M, z):
         """
-        computes the radius R_200 crit of a halo of mass M in physical distances M/h
+        computes the radius R_200 crit of a halo of mass M in physical mass M/h
 
         :param M: halo mass in M_sun/h
         :type M: float or numpy array
@@ -87,13 +91,15 @@ def c_rho0(self, rho0, z):
         :return: concentration parameter c
         """
         if not hasattr(self, '_c_rho0_interp'):
-            c_array = np.linspace(0.1, 10, 100)
+            c_array = np.linspace(0.1, 30, 100)
             rho0_array = self.rho0_c(c_array, z)
             from scipy import interpolate
-            self._c_rho0_interp = interpolate.InterpolatedUnivariateSpline(rho0_array, c_array, w=None, bbox=[None, None], k=3)
+            self._c_rho0_interp = interpolate.InterpolatedUnivariateSpline(rho0_array, c_array, w=None,
+                                                                           bbox=[None, None], k=3)
         return self._c_rho0_interp(rho0)
 
-    def c_M_z(self, M, z):
+    @staticmethod
+    def c_M_z(M, z):
         """
         fitting function of http://moriond.in2p3.fr/J08/proceedings/duffy.pdf for the mass and redshift dependence of
         the concentration parameter
@@ -118,6 +124,9 @@ def nfw_Mz(self, M, z):
         rho_s in  h^2/Mpc^3 (physical)
         Rs in Mpc/h physical
         c unit less
+
+        :param M: Mass in physical M_sun/h
+        :param z: redshift
         """
         c = self.c_M_z(M, z)
         r200 = self.r200_M(M, z)

lenstronomy/ImSim/tracer_model.py

@@ -14,7 +14,7 @@ class TracerModelSource(ImageModel):
     def __init__(self, data_class, psf_class=None, lens_model_class=None, source_model_class=None,
                  lens_light_model_class=None, point_source_class=None, extinction_class=None,
                  tracer_source_class=None, kwargs_numerics=None, likelihood_mask=None,
-                 psf_error_map_bool_list=None, kwargs_pixelbased=None):
+                 psf_error_map_bool_list=None, kwargs_pixelbased=None, tracer_partition=None):
         """
 
         :param data_class: ImageData() instance
@@ -31,11 +31,17 @@ def __init__(self, data_class, psf_class=None, lens_model_class=None, source_mod
          Indicates whether PSF error map is used for the point source model stated as the index.
         :param kwargs_pixelbased: keyword arguments with various settings related to the pixel-based solver
          (see SLITronomy documentation) being applied to the point sources.
+        :param tracer_partition: in case of tracer models for specific sub-parts of the surface brightness model
+         [[list of light profiles, list of tracer profiles], [list of light profiles, list of tracer profiles], [...], ...]
+        :type tracer_partition: None or list
         """
         if likelihood_mask is None:
             likelihood_mask = np.ones_like(data_class.data)
         self.likelihood_mask = np.array(likelihood_mask, dtype=bool)
         self._mask1d = util.image2array(self.likelihood_mask)
+        if tracer_partition is None:
+            tracer_partition = [[None, None]]
+        self._tracer_partition = tracer_partition
         super(TracerModelSource, self).__init__(data_class, psf_class=psf_class, lens_model_class=lens_model_class,
                                                 source_model_class=source_model_class,
                                                 lens_light_model_class=lens_light_model_class,
@@ -60,16 +66,21 @@ def tracer_model(self, kwargs_tracer_source, kwargs_lens, kwargs_source, kwargs_
         :param kwargs_tracer_source:
         :param kwargs_lens:
         :param kwargs_source:
-        :return:
+        :return: model predicted observed tracer component
         """
-        source_light = self._source_surface_brightness_analytical_numerics(kwargs_source, kwargs_lens,
-                                                                           kwargs_extinction,
-                                                                           kwargs_special=kwargs_special,
-                                                                           de_lensed=de_lensed)
-        source_light_conv = self.ImageNumerics.re_size_convolve(source_light, unconvolved=False)
-        source_light_conv[source_light_conv < 10 ** (-20)] = 10 ** (-20)
-        tracer = self._tracer_model_source(kwargs_tracer_source, kwargs_lens, de_lensed=de_lensed)
-        tracer_brightness_conv = self.ImageNumerics.re_size_convolve(tracer * source_light, unconvolved=False)
+        tracer_brightness_conv = np.zeros_like(self.Data.data)
+        source_light_conv = np.zeros_like(self.Data.data)
+        for [k_light, k_tracer] in self._tracer_partition:
+            source_light_k = self._source_surface_brightness_analytical_numerics(kwargs_source, kwargs_lens,
+                                                                               kwargs_extinction,
+                                                                               kwargs_special=kwargs_special,
+                                                                               de_lensed=de_lensed, k=k_light)
+            source_light_conv_k = self.ImageNumerics.re_size_convolve(source_light_k, unconvolved=False)
+            source_light_conv_k[source_light_conv_k < 10 ** (-20)] = 10 ** (-20)
+            tracer_k = self._tracer_model_source(kwargs_tracer_source, kwargs_lens, de_lensed=de_lensed, k=k_tracer)
+            tracer_brightness_conv_k = self.ImageNumerics.re_size_convolve(tracer_k * source_light_k, unconvolved=False)
+            tracer_brightness_conv += tracer_brightness_conv_k
+            source_light_conv += source_light_conv_k
         return tracer_brightness_conv / source_light_conv
 
     def _tracer_model_source(self, kwargs_tracer_source, kwargs_lens, de_lensed=False, k=None):

lenstronomy/LensModel/Profiles/cored_density.py

@@ -11,7 +11,20 @@ class CoredDensity(LensProfileBase):
     """
     class for a uniform cored density dropping steep in the outskirts
     This profile is e.g. featured in Blum et al. 2020 https://arxiv.org/abs/2001.07182v1
-    3d rho(r) = 2/pi * Sigma_crit R_c**3 * (R_c**2 + r**2)**(-2)
+
+    ..math::
+        \\rho_c(r) = \\frac{2}{\\pi} \\Sigma_{c} R_c^3 \\left(R_c^2 + r^2 \\right)^{-2}
+
+    with the convergence profile as
+
+    ..math::
+        \\kappa_c(\\theta) = \\left(1 + \\frac{\\theta^2}{\\theta_c^2} \\right)^{-3/2}.
+
+    An approximate mass-sheet degeneracy can then be written as
+
+    ..math::
+        \\kappa_{\\lambda_c}(\\theta) = \\lambda_c \\kappa(\\theta) + (1-\\lambda_c) \\kappa_c(\\theta).
+
 
     """
     _s = 0.000001  # numerical limit for minimal radius

lenstronomy/LensModel/Profiles/sersic_utils.py

@@ -196,7 +196,7 @@ def _R_stable(self, R):
         """
         return np.maximum(self._smoothing, R)
 
-    def _r_sersic(self, R, R_sersic, n_sersic, max_R_frac=100.0, alpha=1.0, R_break=0.0):
+    def _r_sersic(self, R, R_sersic, n_sersic, max_R_frac=1000.0, alpha=1.0, R_break=0.0):
         """
 
         :param R: radius (array or float)

lenstronomy/LensModel/profile_list_base.py

@@ -52,7 +52,9 @@ def _load_model_instances(self, lens_model_list, custom_class=None, lens_redshif
         func_list = []
         imported_classes = {}
         for i, lens_type in enumerate(lens_model_list):
-            # those models require a new instance per profile as certain pre-computations are relevant per individual profile
+            # those models require a new instance per profile as some pre-computations are different when parameters or
+            # other settings are changed. For example, the 'INTERPOL' model needs to know the specific map to be
+            # interpolated.
             if lens_type in ['NFW_MC', 'CHAMELEON', 'DOUBLE_CHAMELEON', 'TRIPLE_CHAMELEON', 'NFW_ELLIPSE_GAUSS_DEC',
                              'CTNFW_GAUSS_DEC', 'INTERPOL', 'INTERPOL_SCALED', 'NIE', 'NIE_SIMPLE']:
                 lensmodel_class = self._import_class(lens_type, custom_class, z_lens=lens_redshift_list[i],

lenstronomy/LightModel/Profiles/sersic.py

@@ -27,7 +27,7 @@ class Sersic(SersicUtil):
     lower_limit_default = {'amp': 0, 'R_sersic': 0, 'n_sersic': 0.5, 'center_x': -100, 'center_y': -100}
     upper_limit_default = {'amp': 100, 'R_sersic': 100, 'n_sersic': 8, 'center_x': 100, 'center_y': 100}
 
-    def function(self, x, y, amp, R_sersic, n_sersic, center_x=0, center_y=0, max_R_frac=100.0):
+    def function(self, x, y, amp, R_sersic, n_sersic, center_x=0, center_y=0, max_R_frac=1000.0):
         """
 
         :param x:
@@ -37,7 +37,7 @@ def function(self, x, y, amp, R_sersic, n_sersic, center_x=0, center_y=0, max_R_
         :param n_sersic: Sersic index
         :param center_x: center in x-coordinate
         :param center_y: center in y-coordinate
-        :param max_R_frac: maximum window outside of which the mass is zeroed, in units of R_sersic (float)
+        :param max_R_frac: maximum window outside which the mass is zeroed, in units of R_sersic (float)
         :return: Sersic profile value at (x, y)
         """
         R = self.get_distance_from_center(x, y, e1=0, e2=0, center_x=center_x, center_y=center_y)
@@ -66,19 +66,19 @@ class SersicElliptic(SersicUtil):
     upper_limit_default = {'amp': 100, 'R_sersic': 100, 'n_sersic': 8, 'e1': 0.5, 'e2': 0.5, 'center_x': 100,
                            'center_y': 100}
 
-    def function(self, x, y, amp, R_sersic, n_sersic, e1, e2, center_x=0, center_y=0, max_R_frac=100.0):
+    def function(self, x, y, amp, R_sersic, n_sersic, e1, e2, center_x=0, center_y=0, max_R_frac=1000.0):
         """
 
         :param x:
         :param y:
         :param amp: surface brightness/amplitude value at the half light radius
         :param R_sersic: half light radius (either semi-major axis or product average of semi-major and semi-minor axis)
         :param n_sersic: Sersic index
-        :param e1: eccentricity parameter
-        :param e2: eccentricity parameter
+        :param e1: eccentricity parameter e1
+        :param e2: eccentricity parameter e2
         :param center_x: center in x-coordinate
         :param center_y: center in y-coordinate
-        :param max_R_frac: maximum window outside of which the mass is zeroed, in units of R_sersic (float)
+        :param max_R_frac: maximum window outside which the mass is zeroed, in units of R_sersic (float)
         :return: Sersic profile value at (x, y)
         """
 
@@ -91,7 +91,7 @@ def function(self, x, y, amp, R_sersic, n_sersic, e1, e2, center_x=0, center_y=0
 @export
 class CoreSersic(SersicUtil):
     """
-    this class contains the Core-Sersic function introduced by e.g Trujillo et al. 2004
+    this class contains the Core-Sersic function introduced by e.g. Trujillo et al. 2004
 
     .. math::
 
@@ -113,7 +113,7 @@ class CoreSersic(SersicUtil):
                            'center_y': 100}
 
     def function(self, x, y, amp, R_sersic, Rb, n_sersic, gamma, e1, e2, center_x=0, center_y=0, alpha=3.0,
-                 max_R_frac=100.0):
+                 max_R_frac=1000.0):
         """
         :param x:
         :param y:
@@ -122,18 +122,18 @@ def function(self, x, y, amp, R_sersic, Rb, n_sersic, gamma, e1, e2, center_x=0,
         :param Rb: "break" core radius
         :param n_sersic: Sersic index
         :param gamma: inner power-law exponent
-        :param e1: eccentricity parameter
-        :param e2: eccentricity parameter
+        :param e1: eccentricity parameter e1
+        :param e2: eccentricity parameter e2
         :param center_x: center in x-coordinate
         :param center_y: center in y-coordinate
         :param alpha: sharpness of the transition between the cusp and the outer Sersic profile (float)
-        :param max_R_frac: maximum window outside of which the mass is zeroed, in units of R_sersic (float)
+        :param max_R_frac: maximum window outside which the mass is zeroed, in units of R_sersic (float)
         :return: Cored Sersic profile value at (x, y)
         """
         # TODO: max_R_frac not implemented
         R_ = self.get_distance_from_center(x, y, e1, e2, center_x, center_y)
         R = self._R_stable(R_)
         bn = self.b_n(n_sersic)
         result = amp * (1 + (Rb / R) ** alpha) ** (gamma / alpha) * \
-                 np.exp(-bn * (((R ** alpha + Rb ** alpha) / R_sersic ** alpha) ** (1. / (alpha * n_sersic)) - 1.))
+                  np.exp(-bn * (((R ** alpha + Rb ** alpha) / R_sersic ** alpha) ** (1. / (alpha * n_sersic)) - 1.))
         return np.nan_to_num(result)