Merge pull request #38 from sibirrer/kin_scaling_restructuring

major restructuring of parameter distribution and design of the inter…

hierarc/Diagnostics/goodness_of_fit.py

@@ -222,7 +222,7 @@ def plot_kin_fit(
             sigma_v_model_error_list,
         ) = self.kin_fit(cosmo, kwargs_lens, kwargs_kin, kwargs_los)
 
-        f, ax = plt.subplots(1, 1, figsize=(int(len(sigma_v_name_list) / 2), 4))
+        f, ax = plt.subplots(1, 1, figsize=(max(int(len(sigma_v_name_list) / 2), 1), 4))
         ax.errorbar(
             np.arange(len(sigma_v_name_list)),
             sigma_v_measurement_list,

hierarc/LensPosterior/base_config.py

@@ -1,9 +1,9 @@
 from lenstronomy.Analysis.td_cosmography import TDCosmography
 from hierarc.LensPosterior.imaging_constraints import ImageModelPosterior
-from hierarc.LensPosterior.anisotropy_config import AnisotropyConfig
+from hierarc.LensPosterior.kin_scaling_config import KinScalingConfig
 
 
-class BaseLensConfig(TDCosmography, ImageModelPosterior, AnisotropyConfig):
+class BaseLensConfig(TDCosmography, ImageModelPosterior, KinScalingConfig):
     """This class contains and manages the base configurations of the lens posteriors
     and makes sure that they are universally applied consistently through the different
     likelihood definitions."""
@@ -33,6 +33,8 @@ def __init__(
         num_kin_sampling=1000,
         multi_observations=False,
         cosmo_fiducial=None,
+        gamma_in_scaling=None,
+        log_m2l_scaling=None,
     ):
         """
 
@@ -62,6 +64,8 @@ def __init__(
             light profile
         :param cosmo_fiducial: astropy.cosmology instance, if None,
             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)
         """
         self._z_lens, self._z_source = z_lens, z_source
 
@@ -105,4 +109,10 @@ def __init__(
         ImageModelPosterior.__init__(
             self, theta_E, theta_E_error, gamma, gamma_error, r_eff, r_eff_error
         )
-        AnisotropyConfig.__init__(self, anisotropy_model, r_eff)
+        KinScalingConfig.__init__(
+            self,
+            anisotropy_model,
+            r_eff,
+            gamma_in_scaling=gamma_in_scaling,
+            log_m2l_scaling=log_m2l_scaling,
+        )

hierarc/LensPosterior/ddt_kin_constraints.py

@@ -124,7 +124,8 @@ def hierarchy_configuration(self, num_sample_model=20):
             "j_model": j_model_list,
             "error_cov_measurement": error_cov_measurement,
             "error_cov_j_sqrt": error_cov_j_sqrt,
-            "ani_param_array": self.ani_param_array,
-            "ani_scaling_array_list": ani_scaling_array_list,
+            "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_array_list,
         }
         return kwargs_likelihood

hierarc/LensPosterior/ddt_kin_gauss_constraints.py

@@ -123,7 +123,8 @@ def hierarchy_configuration(self, num_sample_model=20):
             "j_model": j_model_list,
             "error_cov_measurement": error_cov_measurement,
             "error_cov_j_sqrt": error_cov_j_sqrt,
-            "ani_param_array": self.ani_param_array,
-            "ani_scaling_array_list": ani_scaling_array_list,
+            "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_array_list,
         }
         return kwargs_likelihood

hierarc/LensPosterior/kin_constraints.py

@@ -36,6 +36,8 @@ def __init__(
         num_kin_sampling=1000,
         multi_observations=False,
         cosmo_fiducial=None,
+        gamma_in_scaling=None,
+        log_m2l_scaling=None,
     ):
         """
 
@@ -77,6 +79,8 @@ def __init__(
             kwargs_seeing as lists of multiple observations
         :param cosmo_fiducial: astropy.cosmology instance, if None,
             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)
         """
         self._sigma_v_measured = np.array(sigma_v_measured)
         self._sigma_v_error_independent = np.array(sigma_v_error_independent)
@@ -111,6 +115,8 @@ def __init__(
             num_kin_sampling=num_kin_sampling,
             multi_observations=multi_observations,
             cosmo_fiducial=cosmo_fiducial,
+            gamma_in_scaling=gamma_in_scaling,
+            log_m2l_scaling=log_m2l_scaling,
         )
 
     def j_kin_draw(self, kwargs_anisotropy, no_error=False):
@@ -171,8 +177,9 @@ def hierarchy_configuration(self, num_sample_model=20):
             "j_model": j_model_list,
             "error_cov_measurement": error_cov_measurement,
             "error_cov_j_sqrt": error_cov_j_sqrt,
-            "ani_param_array": self.ani_param_array,
-            "ani_scaling_array_list": ani_scaling_array_list,
+            "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_array_list,
         }
         return kwargs_likelihood
 
@@ -245,11 +252,16 @@ def _anisotropy_scaling_relative(self, j_ani_0):
 
         if self._anisotropy_model == "GOM":
             ani_scaling_array_list = [
-                np.zeros((len(self.ani_param_array[0]), len(self.ani_param_array[1])))
+                np.zeros(
+                    (
+                        len(self.kin_scaling_param_array[0]),
+                        len(self.kin_scaling_param_array[1]),
+                    )
+                )
                 for _ in range(num_data)
             ]
-            for i, a_ani in enumerate(self.ani_param_array[0]):
-                for j, beta_inf in enumerate(self.ani_param_array[1]):
+            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
                     )
@@ -260,7 +272,7 @@ def _anisotropy_scaling_relative(self, j_ani_0):
                         )  # 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.ani_param_array:
+            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):

hierarc/LensPosterior/kin_constraints_composite.py

@@ -115,6 +115,12 @@ def __init__(
 
         lens_model_list = ["GNFW", "MULTI_GAUSSIAN_KAPPA"]
 
+        # log_m2l is interpolated when sampled on the population level, otherwise marginalized
+        if is_m2l_population_level:
+            log_m2l_scaling = log_m2l_array
+        else:
+            log_m2l_scaling = None
+
         super(KinConstraintsComposite, self).__init__(
             z_lens,
             z_source,
@@ -142,6 +148,8 @@ def __init__(
             num_psf_sampling=num_psf_sampling,
             num_kin_sampling=num_kin_sampling,
             multi_observations=multi_observations,
+            gamma_in_scaling=gamma_in_array,
+            log_m2l_scaling=log_m2l_scaling,
         )
 
         if self._check_arrays(kappa_s_array, r_s_angle_array):
@@ -392,16 +400,19 @@ def hierarchy_configuration(self, num_sample_model=20):
             "j_model": j_model_list,
             "error_cov_measurement": error_cov_measurement,
             "error_cov_j_sqrt": error_cov_j_sqrt,
-            "ani_param_array": self.ani_param_array,
-            "gamma_in_array": self.gamma_in_array,
-            "log_m2l_array": self.log_m2l_array,
-            "param_scaling_grid_list": ani_scaling_grid_list,
+            # "ani_param_array": self.kin_scaling_param_array,
+            # "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,
         }
 
-        if not self._is_m2l_population_level:
-            kwargs_likelihood["log_m2l_array"] = None
+        # if not self._is_m2l_population_level:
+        #    kwargs_likelihood["log_m2l_array"] = None
         return kwargs_likelihood
 
     def anisotropy_scaling(self):
@@ -438,16 +449,16 @@ def _anisotropy_scaling_relative(self, j_ani_0):
             ani_scaling_grid_list = [
                 np.zeros(
                     (
-                        len(self.ani_param_array[0]),
-                        len(self.ani_param_array[1]),
+                        len(self.kin_scaling_param_array[0]),
+                        len(self.kin_scaling_param_array[1]),
                         len(self.gamma_in_array),
                         len(self.log_m2l_array),
                     )
                 )
                 for _ in range(num_data)
             ]
-            for i, a_ani in enumerate(self.ani_param_array[0]):
-                for j, beta_inf in enumerate(self.ani_param_array[1]):
+            for i, a_ani in enumerate(self.kin_scaling_param_array[0]):
+                for j, beta_inf in enumerate(self.kin_scaling_param_array[1]):
                     for k, g_in in enumerate(self.gamma_in_array):
                         for l, log_m2l in enumerate(self.log_m2l_array):
                             kwargs_anisotropy = self.anisotropy_kwargs(
@@ -466,14 +477,14 @@ def _anisotropy_scaling_relative(self, j_ani_0):
             ani_scaling_grid_list = [
                 np.zeros(
                     (
-                        len(self.ani_param_array),
+                        len(self.kin_scaling_param_array[0]),
                         len(self.gamma_in_array),
                         len(self.log_m2l_array),
                     )
                 )
                 for _ in range(num_data)
             ]
-            for i, a_ani in enumerate(self.ani_param_array):
+            for i, a_ani in enumerate(self.kin_scaling_param_array[0]):
                 for k, g_in in enumerate(self.gamma_in_array):
                     for l, log_m2l in enumerate(self.log_m2l_array):
                         kwargs_anisotropy = self.anisotropy_kwargs(a_ani)
@@ -500,15 +511,15 @@ def _anisotropy_scaling_relative_m2l(self, j_ani_0):
             ani_scaling_grid_list = [
                 np.zeros(
                     (
-                        len(self.ani_param_array[0]),
-                        len(self.ani_param_array[1]),
+                        len(self.kin_scaling_param_array[0]),
+                        len(self.kin_scaling_param_array[1]),
                         len(self.gamma_in_array),
                     )
                 )
                 for _ in range(num_data)
             ]
-            for i, a_ani in enumerate(self.ani_param_array[0]):
-                for j, beta_inf in enumerate(self.ani_param_array[1]):
+            for i, a_ani in enumerate(self.kin_scaling_param_array[0]):
+                for j, beta_inf in enumerate(self.kin_scaling_param_array[1]):
                     for k, g_in in enumerate(self.gamma_in_array):
                         kwargs_anisotropy = self.anisotropy_kwargs(
                             a_ani=a_ani, beta_inf=beta_inf
@@ -524,13 +535,13 @@ def _anisotropy_scaling_relative_m2l(self, j_ani_0):
             ani_scaling_grid_list = [
                 np.zeros(
                     (
-                        len(self.ani_param_array),
+                        len(self.kin_scaling_param_array[0]),
                         len(self.gamma_in_array),
                     )
                 )
                 for _ in range(num_data)
             ]
-            for i, a_ani in enumerate(self.ani_param_array):
+            for i, a_ani in enumerate(self.kin_scaling_param_array[0]):
                 for k, g_in in enumerate(self.gamma_in_array):
                     kwargs_anisotropy = self.anisotropy_kwargs(a_ani)
                     j_kin_ani = self.j_kin_draw_composite_m2l(

hierarc/LensPosterior/anisotropy_config.py → hierarc/LensPosterior/kin_scaling_config.py

@@ -1,36 +1,51 @@
 import numpy as np
 
 
-class AnisotropyConfig(object):
+class KinScalingConfig(object):
     """Class to manage the anisotropy model and parameters for the Posterior
     processing."""
 
-    def __init__(self, anisotropy_model, r_eff):
+    def __init__(
+        self, anisotropy_model, r_eff, gamma_in_scaling=None, log_m2l_scaling=None
+    ):
         """
 
         :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)
         """
         self._r_eff = r_eff
         self._anisotropy_model = anisotropy_model
+        self._param_name_list = []
 
         if self._anisotropy_model == "OM":
-            self._ani_param_array = np.array(
-                [0.1, 0.2, 0.5, 1, 2, 5]
-            )  # used for r_ani OsipkovMerritt anisotropy description
+            self._ani_param_array = [np.array([0.1, 0.2, 0.5, 1, 2, 5])]
+            # used for r_ani OsipkovMerritt anisotropy description
+            self._param_name_list = ["a_ani"]
         elif self._anisotropy_model == "GOM":
             self._ani_param_array = [
                 np.array([0.1, 0.2, 0.5, 1, 2, 5]),
                 np.array([0, 0.5, 0.8, 1]),
             ]
+            self._param_name_list = ["a_ani", "beta_inf"]
         elif self._anisotropy_model == "const":
-            self._ani_param_array = np.linspace(
-                -0.49, 1, 7
-            )  # used for constant anisotropy description
+            self._ani_param_array = [
+                np.linspace(-0.49, 1, 7)
+            ]  # used for constant anisotropy description
+            self._param_name_list = ["a_ani"]
+        elif self._anisotropy_model == "NONE":
+            self._param_name_list = []
         else:
             raise ValueError(
                 "anisotropy model %s not supported." % self._anisotropy_model
             )
+        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))
 
     @property
     def kwargs_anisotropy_base(self):
@@ -57,14 +72,22 @@ def kwargs_anisotropy_base(self):
         return kwargs_anisotropy_0
 
     @property
-    def ani_param_array(self):
+    def kin_scaling_param_array(self):
         """
 
         :return: numpy array of anisotropy parameter values to be explored
         """
         return self._ani_param_array
 
-    def anisotropy_kwargs(self, a_ani, beta_inf=None):
+    @property
+    def param_name_list(self):
+        """List of parameters in same order as interpolated.
+
+        :return:
+        """
+        return self._param_name_list
+
+    def anisotropy_kwargs(self, a_ani=None, beta_inf=None):
         """
 
         :param a_ani: anisotropy parameter