Merge pull request #7 from sibirrer/main

Update hierarc fork

.pre-commit-config.yaml

@@ -10,7 +10,7 @@ ci:
 
 repos:
   - repo: https://github.com/psf/black
-    rev: 24.3.0
+    rev: 24.8.0
     hooks:
       - id: black
         # It is recommended to specify the latest version of Python
@@ -19,7 +19,7 @@ repos:
         # https://pre-commit.com/#top_level-default_language_version
         language_version: python3
   - repo: https://github.com/psf/black
-    rev: 24.3.0
+    rev: 24.8.0
     hooks:
       - id: black-jupyter
         language_version: python3

HISTORY.rst

@@ -28,3 +28,8 @@ History
 * double source plane likelihood
 * Pantheon+ likelihood
 * improved API
+
+1.1.3 (2024-06-27)
+------------------
+
+* composite model likelihood and fitting

codecov.yml

@@ -1,3 +1,5 @@
+codecov:
+  token: 399bc344-9039-4f2b-9ade-fee2aac10b78
 comment:                  # this is a top-level key
   layout: " diff, flags, files"
   behavior: default

hierarc/Diagnostics/goodness_of_fit.py

@@ -28,6 +28,7 @@ def plot_ddt_fit(
         cosmo,
         kwargs_lens,
         kwargs_kin,
+        kwargs_los,
         color_measurement=None,
         color_prediction=None,
         redshift_trend=False,
@@ -64,7 +65,9 @@ def plot_ddt_fit(
                     ddt_model_sigma,
                     dd_model_mean,
                     dd_model_sigma,
-                ) = likelihood.ddt_dd_model_prediction(cosmo, kwargs_lens=kwargs_lens)
+                ) = likelihood.ddt_dd_model_prediction(
+                    cosmo, kwargs_lens=kwargs_lens, kwargs_los=kwargs_los
+                )
 
                 ddt_name_list.append(name)
                 ddt_model_mean_list.append(ddt_model_mean)
@@ -134,13 +137,14 @@ def plot_ddt_fit(
         ax.legend()
         return f, ax
 
-    def kin_fit(self, cosmo, kwargs_lens, kwargs_kin):
+    def kin_fit(self, cosmo, kwargs_lens, kwargs_kin, kwargs_los):
         """Plots the prediction and the uncorrelated error bars on the individual lenses
         currently works for likelihood classes 'TDKinGaussian', 'KinGaussian'.
 
         :param cosmo: astropy.cosmology instance
         :param kwargs_lens: lens model parameter keyword arguments
         :param kwargs_kin: kinematics model keyword arguments
+        :param kwargs_los: line of sight list of dictionaries
         :return: list of name, measurement, measurement errors, model prediction, model
             prediction error
         """
@@ -160,7 +164,10 @@ def kin_fit(self, cosmo, kwargs_lens, kwargs_kin):
                 sigma_v_predict_mean,
                 cov_error_predict,
             ) = likelihood.sigma_v_measured_vs_predict(
-                cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin
+                cosmo,
+                kwargs_lens=kwargs_lens,
+                kwargs_kin=kwargs_kin,
+                kwargs_los=kwargs_los,
             )
 
             if sigma_v_measurement is not None:
@@ -188,6 +195,7 @@ def plot_kin_fit(
         cosmo,
         kwargs_lens,
         kwargs_kin,
+        kwargs_los,
         color_measurement=None,
         color_prediction=None,
     ):
@@ -197,21 +205,24 @@ def plot_kin_fit(
         :param cosmo: astropy.cosmology instance
         :param kwargs_lens: lens model parameter keyword arguments
         :param kwargs_kin: kinematics model keyword arguments
+        :param kwargs_los: line of sight list of dictionaries
         :param color_measurement: color of measurement
         :param color_prediction: color of model prediction
         :return: fig, axes of matplotlib instance
         """
-        logL = self._sample_likelihood.log_likelihood(cosmo, kwargs_lens, kwargs_kin)
+        logL = self._sample_likelihood.log_likelihood(
+            cosmo, kwargs_lens, kwargs_kin, kwargs_los=kwargs_los
+        )
         print(logL, "log likelihood")
         (
             sigma_v_name_list,
             sigma_v_measurement_list,
             sigma_v_measurement_error_list,
             sigma_v_model_list,
             sigma_v_model_error_list,
-        ) = self.kin_fit(cosmo, kwargs_lens, kwargs_kin)
+        ) = 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,
@@ -256,6 +267,7 @@ def plot_ifu_fit(
         cosmo,
         kwargs_lens,
         kwargs_kin,
+        kwargs_los,
         lens_index,
         bin_edges,
         show_legend=True,
@@ -268,6 +280,7 @@ def plot_ifu_fit(
         :param cosmo: astropy.cosmology instance
         :param kwargs_lens: lens model parameter keyword arguments
         :param kwargs_kin: kinematics model keyword arguments
+        :param kwargs_los: line of sight list of dictionaries
         :param lens_index: int, index in kwargs_lens to be plotted (needs to be of type
             'IFUKinCov')
         :param bin_edges: radial bin edges in arc seconds. If number, then uniform
@@ -293,7 +306,7 @@ def plot_ifu_fit(
             sigma_v_predict_mean,
             cov_error_predict,
         ) = likelihood.sigma_v_measured_vs_predict(
-            cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin
+            cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin, kwargs_los=kwargs_los
         )
 
         if len(np.atleast_1d(bin_edges)) < 2:

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,9 @@ def __init__(
         num_kin_sampling=1000,
         multi_observations=False,
         cosmo_fiducial=None,
+        gamma_in_scaling=None,
+        log_m2l_scaling=None,
+        gamma_pl_scaling=None,
     ):
         """
 
@@ -62,6 +65,9 @@ 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)
+        :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
 
@@ -105,4 +111,12 @@ 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,
+            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):
@@ -124,7 +127,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,
-            "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,
         }
+
+        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/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/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