Merge pull request #42 from sibirrer/dsp_msd

Double source plane likelihood with the MST effect and sampling of gamma_pl

hierarc/Likelihood/LensLikelihood/base_lens_likelihood.py

@@ -1,5 +1,8 @@
 __author__ = "sibirrer"
 
+from hierarc.Likelihood.LensLikelihood.double_source_plane import (
+    beta_double_source_plane,
+)
 
 LIKELIHOOD_TYPES = [
     "DdtGaussian",
@@ -15,6 +18,7 @@
     "Mag",
     "TDMag",
     "TDMagMagnitude",
+    "DSPL",
 ]
 
 
@@ -26,6 +30,7 @@ def __init__(
         z_lens,
         z_source,
         likelihood_type,
+        z_source2=None,
         name="name",
         normalized=False,
         kwargs_lens_properties=None,
@@ -37,6 +42,7 @@ def __init__(
         :param z_source: source redshift
         :param name: string (optional) to name the specific lens
         :param likelihood_type: string to specify the likelihood type
+        :param z_source2: redshift of the second source for the double source plane lens type "DSP"
         :param normalized: bool, if True, returns the normalized likelihood, if False, separates the constant prefactor
          (in case of a Gaussian 1/(sigma sqrt(2 pi)) ) to compute the reduced chi2 statistics
         :param kwargs_lens_properties: keyword arguments of the lens properties
@@ -46,6 +52,7 @@ def __init__(
         self._name = name
         self.z_lens = z_lens
         self.z_source = z_source
+        self.z_source2 = z_source2
         self.likelihood_type = likelihood_type
         if kwargs_lens_properties is None:
             kwargs_lens_properties = {}
@@ -97,14 +104,16 @@ def __init__(
                 DdtHistLikelihood,
             )
 
-            self._lens_type = DdtHistLikelihood(z_lens, z_source, **kwargs_likelihood)
+            self._lens_type = DdtHistLikelihood(
+                z_lens, z_source, normalized=normalized, **kwargs_likelihood
+            )
         elif likelihood_type == "DdtHistKDE":
             from hierarc.Likelihood.LensLikelihood.ddt_hist_likelihood import (
                 DdtHistKDELikelihood,
             )
 
             self._lens_type = DdtHistKDELikelihood(
-                z_lens, z_source, **kwargs_likelihood
+                z_lens, z_source, normalized=normalized, **kwargs_likelihood
             )
         elif likelihood_type == "DdtHistKin":
             from hierarc.Likelihood.LensLikelihood.ddt_hist_kin_likelihood import (
@@ -140,6 +149,12 @@ def __init__(
             )
 
             self._lens_type = TDMagMagnitudeLikelihood(**kwargs_likelihood)
+        elif likelihood_type == "DSPL":
+            from hierarc.Likelihood.LensLikelihood.double_source_plane import (
+                DSPLikelihood,
+            )
+
+            self._lens_type = DSPLikelihood(normalized=normalized, **kwargs_likelihood)
         else:
             raise ValueError(
                 "likelihood_type %s not supported! Supported are %s."
@@ -154,17 +169,29 @@ def num_data(self):
         return self._lens_type.num_data
 
     def log_likelihood(
-        self, ddt, dd, kin_scaling=None, sigma_v_sys_error=None, mu_intrinsic=None
+        self,
+        ddt,
+        dd,
+        beta_dsp=None,
+        kin_scaling=None,
+        sigma_v_sys_error=None,
+        mu_intrinsic=None,
+        gamma_pl=None,
+        lambda_mst=None,
     ):
         """
 
         :param ddt: time-delay distance [physical Mpc]
         :param dd: angular diameter distance to the lens [physical Mpc]
+        :param beta_dsp: ratio of reduced deflection angles between first and second source redshift,
+         dds1 / ds1 * ds2 / dds2
         :param kin_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted
          dimensionless quantity J (proportional to sigma_v^2) of the anisotropy model in the sampling relative to the
          anisotropy model used to derive the prediction and covariance matrix in the init of this class.
         :param sigma_v_sys_error: unaccounted uncertainty in the velocity dispersion measurement
         :param mu_intrinsic: float, intrinsic source brightness (in magnitude)
+        :param gamma_pl: power-law density slope of main deflector (=2 being isothermal) (only used for DSP likelihood)
+        :param lambda_mst: mass-sheet transform at the main deflector (only used for DSP likelihood)
         :return: natural logarithm of the likelihood of the data given the model
         """
         if self.likelihood_type in [
@@ -187,6 +214,10 @@ def log_likelihood(
             return self._lens_type.log_likelihood(mu_intrinsic=mu_intrinsic)
         elif self.likelihood_type in ["TDMag", "TDMagMagnitude"]:
             return self._lens_type.log_likelihood(ddt=ddt, mu_intrinsic=mu_intrinsic)
+        elif self.likelihood_type in ["DSPL"]:
+            return self._lens_type.log_likelihood(
+                beta_dsp=beta_dsp, gamma_pl=gamma_pl, lambda_mst=lambda_mst
+            )
         else:
             raise ValueError(
                 "likelihood type %s not fully supported." % self.likelihood_type
@@ -231,3 +262,21 @@ def sigma_v_prediction(self, ddt, dd, kin_scaling=None):
         if self.likelihood_type in ["DdtHistKin", "IFUKinCov", "DdtGaussKin"]:
             return self._lens_type.sigma_v_prediction(ddt, dd, kin_scaling)
         return None, None
+
+    def beta_dsp(self, cosmo):
+        """Model prediction of ratio of Einstein radii theta_E_1 / theta_E_2 or scaled
+        deflection angles. Only computes it when likelihood is DSP.
+
+        :param cosmo: ~astropy.cosmology instance
+        :return: beta
+        """
+        if self.likelihood_type == "DSPL":
+            beta = beta_double_source_plane(
+                z_lens=self.z_lens,
+                z_source_1=self.z_source,
+                z_source_2=self.z_source2,
+                cosmo=cosmo,
+            )
+        else:
+            beta = None
+        return beta

hierarc/Likelihood/LensLikelihood/double_source_plane.py

@@ -6,48 +6,38 @@ class DSPLikelihood(object):
 
     def __init__(
         self,
-        z_lens,
-        z_source_1,
-        z_source_2,
         beta_dspl,
         sigma_beta_dspl,
         normalized=False,
     ):
         """
-        :param z_lens: lens redshift
-        :param z_source_1: redshift of first source
-        :param z_source_2: redshift of second source
         :param beta_dspl: measured ratio of Einstein rings theta_E_1 / theta_E_2
-        :param sigma_beta_dspl:
+        :param sigma_beta_dspl: 1-sigma uncertainty in the measurement of the Einstein radius ratio
         :param normalized: normalize the likelihood
         :type normalized: boolean
         """
-        self._z_lens = z_lens
-        self._z_source_1 = z_source_1
-        self._z_source_2 = z_source_2
         self._beta_dspl = beta_dspl
         self._sigma_beta_dspl = sigma_beta_dspl
         self._normalized = normalized
 
-    def lens_log_likelihood(
+    def log_likelihood(
         self,
-        cosmo,
-        kwargs_lens=None,
-        kwargs_kin=None,
-        kwargs_source=None,
-        kwargs_los=None,
+        beta_dsp,
+        gamma_pl=2,
+        lambda_mst=1,
     ):
         """Log likelihood of the data given a model.
 
-        :param cosmo: astropy.cosmology instance
-        :param kwargs_lens: keyword arguments of lens
-        :param kwargs_kin: keyword arguments of kinematics
-        :param kwargs_source: keyword argument of source
-        :param kwargs_los: keyword argument list of line of sight
+        :param beta_dsp: scaled deflection angles alpha_1 / alpha_2 as ratio between
+            z_source and z_source2 source planes
+        :param gamma_pl: power-law density slope of main deflector (=2 being isothermal)
+        :param lambda_mst: mass-sheet transform at the main deflector
         :return: log likelihood of data given model
         """
-        beta_model = self._beta_model(cosmo)
-        log_l = -0.5 * ((beta_model - self._beta_dspl) / self._sigma_beta_dspl) ** 2
+        theta_E_ratio = beta2theta_e_ratio(
+            beta_dsp, gamma_pl=gamma_pl, lambda_mst=lambda_mst
+        )
+        log_l = -0.5 * ((theta_E_ratio - self._beta_dspl) / self._sigma_beta_dspl) ** 2
         if self._normalized:
             log_l -= 1 / 2.0 * np.log(2 * np.pi * self._sigma_beta_dspl**2)
         return log_l
@@ -59,18 +49,6 @@ def num_data(self):
         """
         return 1
 
-    def _beta_model(self, cosmo):
-        """Model prediction of ratio of Einstein radii theta_E_1 / theta_E_2 or scaled
-        deflection angles.
-
-        :param cosmo: astropy.cosmology instance
-        :return: beta
-        """
-        beta = beta_double_source_plane(
-            self._z_lens, self._z_source_1, self._z_source_2, cosmo=cosmo
-        )
-        return beta
-
 
 def beta_double_source_plane(z_lens, z_source_1, z_source_2, cosmo):
     """Model prediction of ratio of Einstein radii theta_E_1 / theta_E_2 or scaled
@@ -79,7 +57,7 @@ def beta_double_source_plane(z_lens, z_source_1, z_source_2, cosmo):
     :param z_lens: lens redshift
     :param z_source_1: source_1 redshift
     :param z_source_2: source_2 redshift
-    :param cosmo: astropy.cosmology instance
+    :param cosmo: ~astropy.cosmology instance
     :return: beta
     """
     ds1 = cosmo.angular_diameter_distance(z=z_source_1).value
@@ -88,3 +66,16 @@ def beta_double_source_plane(z_lens, z_source_1, z_source_2, cosmo):
     dds2 = cosmo.angular_diameter_distance_z1z2(z1=z_lens, z2=z_source_2).value
     beta = dds1 / ds1 * ds2 / dds2
     return beta
+
+
+def beta2theta_e_ratio(beta_dsp, gamma_pl=2, lambda_mst=1):
+    """Calculates Einstein radii ratio for a power-law + MST profile with given
+    parameters.
+
+    :param beta_dsp: scaled deflection angles alpha_1 / alpha_2 as ratio between
+        z_source and z_source2 source planes
+    :param gamma_pl: power-law density slope of main deflector (=2 being isothermal)
+    :param lambda_mst: mass-sheet transform at the main deflector
+    :return: theta_E1 / theta_E2
+    """
+    return (beta_dsp - (1 - lambda_mst) * (1 - beta_dsp)) ** (1 / (gamma_pl - 1))

hierarc/Likelihood/LensLikelihood/mag_likelihood.py

@@ -19,9 +19,9 @@ def __init__(
 
         :param amp_measured: array, amplitudes of measured fluxes of image positions
         :param cov_amp_measured: 2d array, error covariance matrix of the measured amplitudes, in linear space
-        for given magnitude zero point
+         for given magnitude zero point
         :param magnitude_zero_point: magnitude zero point for which the image amplitudes and covariance matrix are
-        defined
+         defined
         :param magnification_model: mean magnification of the model prediction (array with number of images)
         :param cov_magnification_model: 2d array (image amplitudes); model lensing magnification covariances
         """

hierarc/Likelihood/cosmo_likelihood.py

@@ -33,10 +33,10 @@ def __init__(
         :param kwargs_model: model settings for ParamManager() class
         :type kwargs_model: dict
         :param kwargs_bounds: keyword arguments of the lower and upper bounds and parameters that are held fixed.
-        Includes:
-        'kwargs_lower_lens', 'kwargs_upper_lens', 'kwargs_fixed_lens',
-        'kwargs_lower_kin', 'kwargs_upper_kin', 'kwargs_fixed_kin'
-        'kwargs_lower_cosmo', 'kwargs_upper_cosmo', 'kwargs_fixed_cosmo'
+         Includes:
+         'kwargs_lower_lens', 'kwargs_upper_lens', 'kwargs_fixed_lens',
+         'kwargs_lower_kin', 'kwargs_upper_kin', 'kwargs_fixed_kin'
+         'kwargs_lower_cosmo', 'kwargs_upper_cosmo', 'kwargs_fixed_cosmo'
         :param KDE_likelihood_chain: (Likelihood.chain.Chain). Chain object to be evaluated with a kernel density
          estimator
         :param kwargs_kde_likelihood: keyword argument for the KDE likelihood, see KDELikelihood module for options
@@ -99,9 +99,9 @@ def __init__(
             if "z_source" in kwargs_lens:
                 if kwargs_lens["z_source"] > z_max:
                     z_max = kwargs_lens["z_source"]
-            if "z_source_2" in kwargs_lens:
-                if kwargs_lens["z_source_2"] > z_max:
-                    z_max = kwargs_lens["z_source_2"]
+            if "z_source2" in kwargs_lens:
+                if kwargs_lens["z_source2"] > z_max:
+                    z_max = kwargs_lens["z_source2"]
         self._z_max = z_max
 
     def likelihood(self, args, kwargs_cosmo_interp=None):
@@ -124,10 +124,10 @@ def likelihood(self, args, kwargs_cosmo_interp=None):
             # assert we are not in a crazy cosmological situation that prevents computing the angular distance integral
             h0, ok, om = kwargs_cosmo["h0"], kwargs_cosmo["ok"], kwargs_cosmo["om"]
             for lens in self._kwargs_lens_list:
-                if "z_source" in lens:
+                if "z_source2" in lens:
+                    z = lens["z_source2"]
+                elif "z_source" in lens:
                     z = lens["z_source"]
-                elif "z_source_2" in lens:
-                    z = lens["z_source_2"]
                 else:
                     z = 1100
                 cut = ok * (1.0 + z) ** 2 + om * (1.0 + z) ** 3 + (1.0 - om - ok)
@@ -174,11 +174,10 @@ def cosmo_instance(self, kwargs_cosmo):
         """
 
         :param kwargs_cosmo: cosmology parameter keyword argument list
-        :return: astropy.cosmology (or equivalent interpolation scheme class)
+        :return: ~astropy.cosmology (or equivalent interpolation scheme class)
         """
-        if hasattr(kwargs_cosmo, "ang_diameter_distances") and hasattr(
-            kwargs_cosmo, "redshifts"
-        ):
+        print(kwargs_cosmo, "test kwargs_cosmo")
+        if "ang_diameter_distances" in kwargs_cosmo and "redshifts" in kwargs_cosmo:
             # in that case we use directly the interpolation mode to approximate angular diameter distances
             cosmo = CosmoInterp(
                 ang_dist_list=kwargs_cosmo["ang_diameter_distances"],