Merge pull request #40 from sibirrer/gamma_pl_sampling

Processing of kinematics with explicit degree of freedom in the power-law slope

hierarc/LensPosterior/base_config.py

@@ -35,6 +35,7 @@ def __init__(
         cosmo_fiducial=None,
         gamma_in_scaling=None,
         log_m2l_scaling=None,
+        gamma_pl_scaling=None,
     ):
         """
 
@@ -66,6 +67,7 @@ def __init__(
             uses astropy's default cosmology
         :param gamma_in_scaling: array of gamma_in parameter to be interpolated (optional, otherwise None)
         :param log_m2l_scaling: array of log_m2l parameter to be interpolated (optional, otherwise None)
+        :param gamma_pl_scaling: array of power-law density profile slopes to be interpolated (optional, otherwise None)
         """
         self._z_lens, self._z_source = z_lens, z_source
 
@@ -115,4 +117,6 @@ def __init__(
             r_eff,
             gamma_in_scaling=gamma_in_scaling,
             log_m2l_scaling=log_m2l_scaling,
+            gamma_pl_scaling=gamma_pl_scaling,
+            gamma_pl_mean=gamma,
         )

hierarc/LensPosterior/ddt_kin_constraints.py

@@ -36,6 +36,7 @@ def __init__(
         num_psf_sampling=100,
         num_kin_sampling=1000,
         multi_observations=False,
+        gamma_pl_scaling=None,
     ):
         """
 
@@ -67,6 +68,7 @@ def __init__(
         :param kappa_ext_sigma: 1-sigma distribution uncertainty from which the ddt constraints are coming from
         :param multi_observations: bool, if True, interprets kwargs_aperture and kwargs_seeing as lists of multiple
          observations
+        :param gamma_pl_scaling: array of mass density profile power-law slope values (optional, otherwise None)
         """
         self._ddt_sample, self._ddt_weights = ddt_samples, ddt_weights
         self._kappa_ext_mean, self._kappa_ext_sigma = kappa_ext, kappa_ext_sigma
@@ -96,6 +98,7 @@ def __init__(
             num_psf_sampling=num_psf_sampling,
             num_kin_sampling=num_kin_sampling,
             multi_observations=multi_observations,
+            gamma_pl_scaling=gamma_pl_scaling,
         )
 
     def hierarchy_configuration(self, num_sample_model=20):
@@ -128,4 +131,15 @@ def hierarchy_configuration(self, num_sample_model=20):
             "j_kin_scaling_param_axes": self.kin_scaling_param_array,
             "j_kin_scaling_grid_list": ani_scaling_array_list,
         }
+
+        prior_list = []
+        if "gamma_pl" in self._param_name_list:
+            prior_list.append(["gamma_pl", self._gamma, self._gamma_error])
+        # TODO: make sure to add other priors if needed or available
+        # if "gamma_in" in self._param_name_list:
+        #    prior_list.append(["gamma_in"])
+        kwargs_likelihood["prior_list"] = prior_list
+        # if "gamma_pl" in self._param_name_list:
+        #    kwargs_likelihood["gamma_pl_sampling"] = True
+
         return kwargs_likelihood

hierarc/LensPosterior/imaging_constraints.py

@@ -18,23 +18,28 @@ def __init__(self, theta_E, theta_E_error, gamma, gamma_error, r_eff, r_eff_erro
         self._gamma, self._gamma_error = gamma, gamma_error
         self._r_eff, self._r_eff_error = r_eff, r_eff_error
 
-    def draw_lens(self, no_error=False):
+    def draw_lens(self, gamma_pl=None, no_error=False):
         """
 
         :param no_error: bool, if True, does not render from the uncertainty but uses the mean values instead
+        :param gamma_pl: power law slope, if None, draws from measurement uncertainty, otherwise takes at fixed value
+        :type gamma_pl: float or None
         :return: theta_E, gamma, r_eff, delta_r_eff
         """
         if no_error is True:
             return self._theta_E, self._gamma, self._r_eff, 1
         theta_E_draw = np.maximum(
             np.random.normal(loc=self._theta_E, scale=self._theta_E_error), 0
         )
-        gamma_draw = np.random.normal(loc=self._gamma, scale=self._gamma_error)
-        # distributions are drawn in the range [1, 3)
-        # the power-law slope gamma=3 is divergent in mass in the center and values close close to =3 may be unstable
-        # to compute the kinematics for.
-        gamma_draw = np.maximum(gamma_draw, 1.0)
-        gamma_draw = np.minimum(gamma_draw, 2.999)
+        if gamma_pl is None:
+            gamma_draw = np.random.normal(loc=self._gamma, scale=self._gamma_error)
+            # distributions are drawn in the range [1, 3)
+            # the power-law slope gamma=3 is divergent in mass in the center and values close to =3 may be unstable
+            # to compute the kinematics for.
+            gamma_draw = np.maximum(gamma_draw, 1.0)
+            gamma_draw = np.minimum(gamma_draw, 2.999)
+        else:
+            gamma_draw = gamma_pl
         # we make sure no negative r_eff are being sampled
         delta_r_eff = np.maximum(
             np.random.normal(loc=1, scale=self._r_eff_error / self._r_eff), 0.001

hierarc/LensPosterior/kin_constraints.py

@@ -38,6 +38,7 @@ def __init__(
         cosmo_fiducial=None,
         gamma_in_scaling=None,
         log_m2l_scaling=None,
+        gamma_pl_scaling=None,
     ):
         """
 
@@ -81,6 +82,7 @@ def __init__(
             uses astropy's default
         :param gamma_in_scaling: array of gamma_in parameter to be interpolated (optional, otherwise None)
         :param log_m2l_scaling: array of log_m2l parameter to be interpolated (optional, otherwise None)
+        :param gamma_pl_scaling: array of mass density profile power-law slope values (optional, otherwise None)
         """
         self._sigma_v_measured = np.array(sigma_v_measured)
         self._sigma_v_error_independent = np.array(sigma_v_error_independent)
@@ -117,18 +119,22 @@ def __init__(
             cosmo_fiducial=cosmo_fiducial,
             gamma_in_scaling=gamma_in_scaling,
             log_m2l_scaling=log_m2l_scaling,
+            gamma_pl_scaling=gamma_pl_scaling,
         )
 
-    def j_kin_draw(self, kwargs_anisotropy, no_error=False):
+    def j_kin_draw(self, kwargs_anisotropy, gamma_pl=None, no_error=False):
         """One simple sampling realization of the dimensionless kinematics of the model.
 
         :param kwargs_anisotropy: keyword argument of anisotropy setting
+        :param gamma_pl: power law slope, if None, draws from measurement uncertainty,
+            otherwise takes at fixed value
+        :type gamma_pl: float or None
         :param no_error: bool, if True, does not render from the uncertainty but uses
             the mean values instead
         :return: dimensionless kinematic component J() Birrer et al. 2016, 2019
         """
         theta_E_draw, gamma_draw, r_eff_draw, delta_r_eff = self.draw_lens(
-            no_error=no_error
+            gamma_pl=gamma_pl, no_error=no_error
         )
         kwargs_lens = [
             {"theta_E": theta_E_draw, "gamma": gamma_draw, "center_x": 0, "center_y": 0}
@@ -181,6 +187,15 @@ def hierarchy_configuration(self, num_sample_model=20):
             "j_kin_scaling_param_axes": self.kin_scaling_param_array,
             "j_kin_scaling_grid_list": ani_scaling_array_list,
         }
+        prior_list = []
+        if "gamma_pl" in self.param_name_list:
+            prior_list.append(["gamma_pl", self._gamma, self._gamma_error])
+        # TODO: make sure to add other priors if needed or available
+        # if "gamma_in" in self._param_name_list:
+        #    prior_list.append(["gamma_in"])
+        kwargs_likelihood["prior_list"] = prior_list
+        # if "gamma_pl" in self.param_name_list:
+        #    kwargs_likelihood["gamma_pl_sampling"] = True
         return kwargs_likelihood
 
     def model_marginalization(self, num_sample_model=20):
@@ -196,7 +211,9 @@ def model_marginalization(self, num_sample_model=20):
             (num_sample_model, num_data)
         )  # matrix that contains the sampled J() distribution
         for i in range(num_sample_model):
-            j_kin = self.j_kin_draw(self.kwargs_anisotropy_base, no_error=False)
+            j_kin = self.j_kin_draw(
+                self.kwargs_anisotropy_base, no_error=False, **self.kwargs_lens_base
+            )
             j_kin_matrix[i, :] = j_kin
 
         error_cov_j_sqrt = np.cov(np.sqrt(j_kin_matrix.T))
@@ -238,7 +255,9 @@ def anisotropy_scaling(self):
 
         :return: anisotropy scaling grid along the axes defined by ani_param_array
         """
-        j_ani_0 = self.j_kin_draw(self.kwargs_anisotropy_base, no_error=True)
+        j_ani_0 = self.j_kin_draw(
+            self.kwargs_anisotropy_base, no_error=True, **self.kwargs_lens_base
+        )
         return self._anisotropy_scaling_relative(j_ani_0)
 
     def _anisotropy_scaling_relative(self, j_ani_0):
@@ -249,34 +268,52 @@ def _anisotropy_scaling_relative(self, j_ani_0):
             scaling
         """
         num_data = len(self._sigma_v_measured)
-
-        if self._anisotropy_model == "GOM":
-            ani_scaling_array_list = [
-                np.zeros(
-                    (
-                        len(self.kin_scaling_param_array[0]),
-                        len(self.kin_scaling_param_array[1]),
-                    )
+        len_list = [len(a) for a in self.kin_scaling_param_array]
+        ani_scaling_array_list = [np.zeros(len_list) for _ in range(num_data)]
+        num = self.num_scaling_dim
+        if num == 1:
+            for i, param in enumerate(self.kin_scaling_param_array[0]):
+                param_array = [param]
+                kwargs_ani, kwargs_lens = self.param_array2kwargs(
+                    param_array=param_array
                 )
-                for _ in range(num_data)
-            ]
-            for i, a_ani in enumerate(self.kin_scaling_param_array[0]):
-                for j, beta_inf in enumerate(self.kin_scaling_param_array[1]):
-                    kwargs_anisotropy = self.anisotropy_kwargs(
-                        a_ani=a_ani, beta_inf=beta_inf
+                kwargs_anisotropy = self.anisotropy_kwargs(**kwargs_ani)
+                j_kin_ani = self.j_kin_draw(
+                    kwargs_anisotropy, no_error=True, **kwargs_lens
+                )
+                for s, j_kin in enumerate(j_kin_ani):
+                    ani_scaling_array_list[s][i] = j_kin / j_ani_0[s]
+        elif num == 2:
+            for i, param_i in enumerate(self.kin_scaling_param_array[0]):
+                for j, param_j in enumerate(self.kin_scaling_param_array[1]):
+                    param_array = [param_i, param_j]
+                    kwargs_ani, kwargs_lens = self.param_array2kwargs(
+                        param_array=param_array
+                    )
+                    kwargs_anisotropy = self.anisotropy_kwargs(**kwargs_ani)
+                    j_kin_ani = self.j_kin_draw(
+                        kwargs_anisotropy, no_error=True, **kwargs_lens
                     )
-                    j_kin_ani = self.j_kin_draw(kwargs_anisotropy, no_error=True)
-                    for k, j_kin in enumerate(j_kin_ani):
-                        ani_scaling_array_list[k][i, j] = (
-                            j_kin / j_ani_0[k]
-                        )  # perhaps change the order
-        elif self._anisotropy_model in ["OM", "const"]:
-            ani_scaling_array_list = [[] for _ in range(num_data)]
-            for a_ani in self.kin_scaling_param_array[0]:
-                kwargs_anisotropy = self.anisotropy_kwargs(a_ani)
-                j_kin_ani = self.j_kin_draw(kwargs_anisotropy, no_error=True)
-                for i, j_kin in enumerate(j_kin_ani):
-                    ani_scaling_array_list[i].append(j_kin / j_ani_0[i])
+                    for s, j_kin in enumerate(j_kin_ani):
+                        ani_scaling_array_list[s][i, j] = j_kin / j_ani_0[s]
+        elif num == 3:
+            for i, param_i in enumerate(self.kin_scaling_param_array[0]):
+                for j, param_j in enumerate(self.kin_scaling_param_array[1]):
+                    for k, param_k in enumerate(self.kin_scaling_param_array[2]):
+                        param_array = [param_i, param_j, param_k]
+                        kwargs_ani, kwargs_lens = self.param_array2kwargs(
+                            param_array=param_array
+                        )
+                        kwargs_anisotropy = self.anisotropy_kwargs(**kwargs_ani)
+                        j_kin_ani = self.j_kin_draw(
+                            kwargs_anisotropy, no_error=True, **kwargs_lens
+                        )
+                        for s, j_kin in enumerate(j_kin_ani):
+                            ani_scaling_array_list[s][i, j, k] = j_kin / j_ani_0[s]
         else:
-            raise ValueError("anisotropy model %s not valid." % self._anisotropy_model)
+            ValueError(
+                "Kin scaling with parameter dimension %s not supported, chose between 1-3."
+                % num
+            )
+
         return ani_scaling_array_list

hierarc/LensPosterior/kin_constraints_composite.py

@@ -404,13 +404,20 @@ def hierarchy_configuration(self, num_sample_model=20):
             # "gamma_in_array": self.gamma_in_array,
             # "log_m2l_array": self.log_m2l_array,
             # "param_scaling_grid_list": ani_scaling_grid_list,
-            "gamma_in_prior_mean": self._gamma_in_prior_mean,
-            "gamma_in_prior_std": self._gamma_in_prior_std,
             "kin_scaling_param_list": self.param_name_list,
             "j_kin_scaling_param_axes": self.kin_scaling_param_array,
             "j_kin_scaling_grid_list": ani_scaling_grid_list,
         }
-
+        prior_list = None
+        if (
+            self.gamma_in_array is not None
+            and self._gamma_in_prior_mean is not None
+            and self._gamma_in_prior_std is not None
+        ):
+            prior_list = [
+                ["gamma_in", self._gamma_in_prior_mean, self._gamma_in_prior_std]
+            ]
+        kwargs_likelihood["prior_list"] = prior_list
         # if not self._is_m2l_population_level:
         #    kwargs_likelihood["log_m2l_array"] = None
         return kwargs_likelihood

hierarc/LensPosterior/kin_scaling_config.py

@@ -1,19 +1,29 @@
 import numpy as np
+from hierarc.Likelihood.kin_scaling import KinScalingParamManager
 
 
-class KinScalingConfig(object):
+class KinScalingConfig(KinScalingParamManager):
     """Class to manage the anisotropy model and parameters for the Posterior
     processing."""
 
     def __init__(
-        self, anisotropy_model, r_eff, gamma_in_scaling=None, log_m2l_scaling=None
+        self,
+        anisotropy_model,
+        r_eff,
+        gamma_in_scaling=None,
+        log_m2l_scaling=None,
+        gamma_pl_scaling=None,
+        gamma_pl_mean=None,
     ):
         """
 
-        :param anisotropy_model: type of stellar anisotropy model. Supported are 'OM' and 'GOM' or 'const', see details in lenstronomy.Galkin module
+        :param anisotropy_model: type of stellar anisotropy model. Supported are 'OM' and 'GOM' or 'const',
+         see details in lenstronomy.Galkin module
         :param r_eff: half-light radius of the deflector galaxy
         :param gamma_in_scaling: array of gamma_in parameter to be interpolated (optional, otherwise None)
         :param log_m2l_scaling: array of log_m2l parameter to be interpolated (optional, otherwise None)
+        :param gamma_pl_scaling: array of power-law density profile slopes to be interpolated (optional, otherwise None)
+        :param gamma_pl_mean: mean gamma_pl upon which the covariances are calculated
         """
         self._r_eff = r_eff
         self._anisotropy_model = anisotropy_model
@@ -40,12 +50,23 @@ def __init__(
             raise ValueError(
                 "anisotropy model %s not supported." % self._anisotropy_model
             )
+        self._gamma_in_scaling = gamma_in_scaling
+        self._log_m2l_scaling = log_m2l_scaling
+        self._gamma_pl_scaling = gamma_pl_scaling
+        self._gamma_pl_mean = gamma_pl_mean
+
         if gamma_in_scaling is not None:
             self._param_name_list.append("gamma_in")
             self._ani_param_array.append(np.array(gamma_in_scaling))
         if log_m2l_scaling is not None:
             self._param_name_list.append("log_m2l")
             self._ani_param_array.append(np.array(log_m2l_scaling))
+        if gamma_pl_scaling is not None:
+            self._param_name_list.append("gamma_pl")
+            self._ani_param_array.append(np.array(gamma_pl_scaling))
+        KinScalingParamManager.__init__(
+            self, j_kin_scaling_param_name_list=self._param_name_list
+        )
 
     @property
     def kwargs_anisotropy_base(self):
@@ -71,11 +92,26 @@ def kwargs_anisotropy_base(self):
             )
         return kwargs_anisotropy_0
 
+    @property
+    def kwargs_lens_base(self):
+        """
+
+        :return: keyword arguments of lens model parameters that are getting interpolated
+        """
+        kwargs_base = {}
+        if "gamma_in" in self._param_name_list:
+            kwargs_base["gamma_in"] = np.mean(self._gamma_in_scaling)
+        if "log_m2l" in self._param_name_list:
+            kwargs_base["log_m2l"] = np.mean(self._log_m2l_scaling)
+        if "gamma_pl" in self._param_name_list:
+            kwargs_base["gamma_pl"] = self._gamma_pl_mean
+        return kwargs_base
+
     @property
     def kin_scaling_param_array(self):
         """
 
-        :return: numpy array of anisotropy parameter values to be explored
+        :return: numpy array of kinematic scaling parameter values to be explored, list of 1D arrays
         """
         return self._ani_param_array
 
@@ -92,7 +128,7 @@ def anisotropy_kwargs(self, a_ani=None, beta_inf=None):
 
         :param a_ani: anisotropy parameter
         :param beta_inf: anisotropy at infinity (only used for 'GOM' model)
-        :return: list of anisotropy keyword arguments, value of anisotropy parameter list
+        :return: list of anisotropy keyword arguments for GalKin module
         """
 
         if self._anisotropy_model == "OM":

hierarc/Likelihood/cosmo_likelihood.py

@@ -60,7 +60,10 @@ def __init__(
             normalized=normalized,
             kwargs_global_model=kwargs_model,
         )
-        self.param = ParamManager(cosmology, **kwargs_model, **kwargs_bounds)
+        gamma_pl_num = self._likelihoodLensSample.gamma_pl_num
+        self.param = ParamManager(
+            cosmology, gamma_pl_num=gamma_pl_num, **kwargs_model, **kwargs_bounds
+        )
         self._lower_limit, self._upper_limit = self.param.param_bounds
         self._prior_add = False
         if custom_prior is not None: