Repository: sibirrer/hierArc (language: python)

Issue:
Sklearn KDE not working as attempted
https://github.com/sibirrer/hierArc/blob/0a2148460002548992a54f3d9d0276a74e0d9607/hierarc/Likelihood/LensLikelihood/ddt_hist_likelihood.py#L39

This routine originates from https://github.com/shsuyu/H0LiCOW-public/blob/master/H0_inference_code/lensutils.py but did not work with the current sklearn version. @aymgal @martin-millon  do you know whether the notebooks for the sampling of the Wong et al. paper are a jour with the latest versions? In the meantime I replaced it with the scipy gaussian kde which may not have the accuracy we require.


Discussion / hints:
With what I have in hand I cannot answer but @martin-millon sure will.
It was working for me with lenstronomy 1.3.0. I did not update it to the latest version yet. Let me run it again with the latest version and see if I have the same results. 
@martin-millon I am not pointing on a lenstronomy dependent routine. It is the sklearn KDE with weights attached, effectively line 185 or 218 in the file on the H0LiCOW repo of the definition likelihood_hist()
OK, sorry, I am running the script with sklearn v0.21.3. The latest version of sklearn is now 0.22.2 and I did not update the script to the latest version. 
Thank you! I could fix it. The requirement was `scikit-learn` and I installed `sklearn`
OK, I'll change the requirements in the documentation then. 

Fix patch:
diff --git a/hierarc/LensPosterior/kin_constraints_composite.py b/hierarc/LensPosterior/kin_constraints_composite.py
index 9d00ab4..f43b62b 100644
--- a/hierarc/LensPosterior/kin_constraints_composite.py
+++ b/hierarc/LensPosterior/kin_constraints_composite.py
@@ -7,6 +7,7 @@
 from lenstronomy.Util import constants as const
 from lenstronomy.Analysis.light_profile import LightProfileAnalysis
 from lenstronomy.LightModel.light_model import LightModel
+from lenstronomy.LensModel.Profiles.gnfw import GNFW
 
 
 class KinConstraintsComposite(KinConstraints):
@@ -16,7 +17,7 @@ def __init__(
         z_source,
         gamma_in_array,
         log_m2l_array,
-        kappa_s_array,
+        alpha_Rs_array,
         r_s_angle_array,
         theta_E,
         theta_E_error,
@@ -42,6 +43,7 @@ def __init__(
         num_kin_sampling=1000,
         multi_observations=False,
         rho0_array=None,
+        kappa_s_array=None,
         r_s_array=None,
         is_m2l_population_level=True,
     ):
@@ -53,7 +55,7 @@ def __init__(
         :param log_m2l_array: array of log10(mass-to-light ratios) of the stellar component,
             needs to be in the unit/scaling such that m2l / sigma_crit * amp in the
             kwargs_lens_light provides the convergence amplitude of the stars
-        :param kappa_s_array: array of generalized NFW profile's convergence normalization at the scale radius
+        :param alpha_Rs_array: array of the deflection (angular units) at projected Rs
         :param r_s_angle_array: array of halo scale radii in arcsecond
         :param theta_E: Einstein radius (in arc seconds)
         :param theta_E_error: 1-sigma error on Einstein radius
@@ -88,6 +90,7 @@ def __init__(
         :param multi_observations: bool, if True, interprets kwargs_aperture and
             kwargs_seeing as lists of multiple observations
         :param rho0_array: array of halo mass normalizations in M_sun / Mpc^3
+        :param kappa_s_array: array of generalized NFW profile's convergence normalization at the scale radius
         :param r_s_array: array of halo scale radii in Mpc
         """
 
@@ -152,16 +155,23 @@ def __init__(
             log_m2l_scaling=log_m2l_scaling,
         )
 
-        if self._check_arrays(kappa_s_array, r_s_angle_array):
-            self._kappa_s_array = kappa_s_array
+        if self._check_arrays(alpha_Rs_array, r_s_angle_array):
+            self._halo_normalization_array = alpha_Rs_array
+            self._is_normalization_alpha_Rs = True
+            self._r_scale_angle_array = r_s_angle_array
+        elif self._check_arrays(kappa_s_array, r_s_angle_array):
+            self._halo_normalization_array = kappa_s_array
+            self._is_normalization_alpha_Rs = False
             self._r_scale_angle_array = r_s_angle_array
         elif self._check_arrays(rho0_array, r_s_array):
-            self._kappa_s_array, self._r_scale_angle_array = self.get_kappa_s_r_s_angle(
+            kappa_s_array, self._r_scale_angle_array = self.get_kappa_s_r_s_angle(
                 rho0_array, r_s_array
             )
+            self._halo_normalization_array = kappa_s_array
+            self._is_normalization_alpha_Rs = False
         else:
             raise ValueError(
-                "Both kappa_s_array and r_s_angle_array, or rho0_array and r_s_array must be arrays of the same length!"
+                "Both alpha_Rs_array and r_s_angle_array, kappa_s_array and r_s_angle_array, or rho0_array and r_s_array must be arrays of the same length!"
             )
 
         self.gamma_in_array = gamma_in_array
@@ -172,10 +182,10 @@ def __init__(
         self._gamma_in_prior_std = gamma_in_prior_std
 
         if not is_m2l_population_level and not self._check_arrays(
-            self._kappa_s_array, log_m2l_array
+            self._halo_normalization_array, log_m2l_array
         ):
             raise ValueError(
-                "log_m2l_array must have the same length as rho0_array or kappa_s_array!"
+                "log_m2l_array must have the same length as alpha_Rs_array, kappa_s_array, or rho0_array!"
             )
 
     @staticmethod
@@ -210,22 +220,24 @@ def draw_lens(self, no_error=False):
         if no_error is True:
             if self._is_m2l_population_level:
                 return (
-                    np.mean(self._kappa_s_array),
+                    np.mean(self._halo_normalization_array),
                     np.mean(self._r_scale_angle_array),
                     self._r_eff,
                     1,
                 )
             else:
                 return (
-                    np.mean(self._kappa_s_array),
+                    np.mean(self._halo_normalization_array),
                     np.mean(self._r_scale_angle_array),
                     np.mean(self.log_m2l_array),
                     self._r_eff,
                     1,
                 )
 
-        random_index = np.random.randint(low=0, high=len(self._kappa_s_array))
-        kappa_s_draw = self._kappa_s_array[random_index]
+        random_index = np.random.randint(
+            low=0, high=len(self._halo_normalization_array)
+        )
+        halo_normalization_draw = self._halo_normalization_array[random_index]
         r_scale_angle_draw = self._r_scale_angle_array[random_index]
 
         # we make sure no negative r_eff are being sampled
@@ -235,11 +247,11 @@ def draw_lens(self, no_error=False):
         r_eff_draw = delta_r_eff * self._r_eff
 
         if self._is_m2l_population_level:
-            return kappa_s_draw, r_scale_angle_draw, r_eff_draw, delta_r_eff
+            return halo_normalization_draw, r_scale_angle_draw, r_eff_draw, delta_r_eff
         else:
             log_m2l_draw = self.log_m2l_array[random_index]
             return (
-                kappa_s_draw,
+                halo_normalization_draw,
                 r_scale_angle_draw,
                 log_m2l_draw,
                 r_eff_draw,
@@ -290,8 +302,8 @@ def j_kin_draw_composite(
             the mean values instead
         :return: dimensionless kinematic component J() Birrer et al. 2016, 2019
         """
-        kappa_s_draw, r_scale_angle_draw, r_eff_draw, delta_r_eff = self.draw_lens(
-            no_error=no_error
+        halo_normalization_draw, r_scale_angle_draw, r_eff_draw, delta_r_eff = (
+            self.draw_lens(no_error=no_error)
         )
 
         kwargs_lens_stars = copy.deepcopy(self._kwargs_lens_light[0])
@@ -304,11 +316,20 @@ def j_kin_draw_composite(
         for kwargs in kwargs_light:
             kwargs["sigma"] *= delta_r_eff
 
+        # Input is alpha_Rs
+        if self._is_normalization_alpha_Rs:
+            alpha_Rs_draw = halo_normalization_draw
+        # Input is kappa_s
+        else:
+            alpha_Rs_draw = GNFW().kappa_s_to_alpha_Rs(
+                halo_normalization_draw, r_scale_angle_draw, gamma_in
+            )
+
         kwargs_lens = [
             {
                 "Rs": r_scale_angle_draw,
                 "gamma_in": gamma_in,
-                "kappa_s": kappa_s_draw,
+                "alpha_Rs": alpha_Rs_draw,
                 "center_x": 0,
                 "center_y": 0,
             },
@@ -336,7 +357,7 @@ def j_kin_draw_composite_m2l(self, kwargs_anisotropy, gamma_in, no_error=False):
         :return: dimensionless kinematic component J() Birrer et al. 2016, 2019
         """
         (
-            kappa_s_draw,
+            halo_normalization_draw,
             r_scale_angle_draw,
             log_m2l_draw,
             r_eff_draw,
@@ -353,11 +374,20 @@ def j_kin_draw_composite_m2l(self, kwargs_anisotropy, gamma_in, no_error=False):
         for kwargs in kwargs_light:
             kwargs["sigma"] *= delta_r_eff
 
+        # Input is alpha_Rs
+        if self._is_normalization_alpha_Rs:
+            alpha_Rs_draw = halo_normalization_draw
+        # Input is kappa_s
+        else:
+            alpha_Rs_draw = GNFW().kappa_s_to_alpha_Rs(
+                halo_normalization_draw, r_scale_angle_draw, gamma_in
+            )
+
         kwargs_lens = [
             {
                 "Rs": r_scale_angle_draw,
                 "gamma_in": gamma_in,
-                "kappa_s": kappa_s_draw,
+                "alpha_Rs": alpha_Rs_draw,
                 "center_x": 0,
                 "center_y": 0,
             },


Test patch:
diff --git a/test/test_LensPosterior/test_kin_constraints_composite.py b/test/test_LensPosterior/test_kin_constraints_composite.py
index 5b7b779..ef39153 100644
--- a/test/test_LensPosterior/test_kin_constraints_composite.py
+++ b/test/test_LensPosterior/test_kin_constraints_composite.py
@@ -81,6 +81,7 @@ def test_likelihoodconfiguration_om(self):
             z_source=z_source,
             gamma_in_array=gamma_in_array,
             log_m2l_array=log_m2l_array,
+            alpha_Rs_array=[],
             kappa_s_array=rho0_array,
             r_s_angle_array=r_s_array,
             theta_E=theta_E,
@@ -120,6 +121,7 @@ def test_likelihoodconfiguration_om(self):
             z_source=z_source,
             gamma_in_array=gamma_in_array,
             log_m2l_array=log_m2l_array,
+            alpha_Rs_array=[],
             kappa_s_array=rho0_array,
             r_s_angle_array=r_s_array,
             theta_E=theta_E,
@@ -147,6 +149,7 @@ def test_likelihoodconfiguration_om(self):
             z_source=z_source,
             gamma_in_array=gamma_in_array,
             log_m2l_array=log_m2l_array,
+            alpha_Rs_array=[],
             kappa_s_array=[],
             r_s_angle_array=[],
             theta_E=theta_E,
@@ -237,6 +240,7 @@ def test_likelihoodconfiguration_gom(self):
             z_source=z_source,
             gamma_in_array=gamma_in_array,
             log_m2l_array=log_m2l_array,
+            alpha_Rs_array=[],
             kappa_s_array=rho0_array,
             r_s_angle_array=r_s_array,
             theta_E=theta_E,
@@ -342,6 +346,7 @@ def test_likelihoodconfiguration_om(self):
             z_source=z_source,
             gamma_in_array=gamma_in_array,
             log_m2l_array=log_m2l_array,
+            alpha_Rs_array=[],
             kappa_s_array=rho0_array,
             r_s_angle_array=r_s_array,
             theta_E=theta_E,
@@ -445,6 +450,7 @@ def test_likelihoodconfiguration_gom(self):
             z_source=z_source,
             gamma_in_array=gamma_in_array,
             log_m2l_array=log_m2l_array,
+            alpha_Rs_array=[],
             kappa_s_array=rho0_array,
             r_s_angle_array=r_s_array,
             theta_E=theta_E,
@@ -474,6 +480,221 @@ def test_likelihoodconfiguration_gom(self):
         )
 
 
+class TestKinConstraintsCompositeAlphaRs(object):
+    def setup_method(self):
+        pass
+
+    def test_likelihoodconfiguration_om(self):
+        anisotropy_model = "OM"
+        kwargs_aperture = {
+            "aperture_type": "shell",
+            "r_in": 0,
+            "r_out": 3 / 2.0,
+            "center_ra": 0.0,
+            "center_dec": 0,
+        }
+        kwargs_seeing = {"psf_type": "GAUSSIAN", "fwhm": 1.4}
+
+        # numerical settings (not needed if power-law profiles with Hernquist light distribution is computed)
+        kwargs_numerics_galkin = {
+            "interpol_grid_num": 50,  # numerical interpolation, should converge -> infinity
+            "log_integration": True,
+            # log or linear interpolation of surface brightness and mass models
+            "max_integrate": 100,
+            "min_integrate": 0.001,
+        }  # lower/upper bound of numerical integrals
+
+        # redshift
+        z_lens = 0.5
+        z_source = 1.5
+
+        cosmo = FlatLambdaCDM(H0=70, Om0=0.3)
+        theta_E = 1.0
+        r_eff = 1
+        gamma = 2.1
+
+        kwargs_mge_light = {
+            "grid_spacing": 0.1,
+            "grid_num": 10,
+            "n_comp": 20,
+            "center_x": 0,
+            "center_y": 0,
+        }
+
+        kwargs_kin_api_settings = {
+            "multi_observations": False,
+            "kwargs_numerics_galkin": kwargs_numerics_galkin,
+            "kwargs_mge_light": kwargs_mge_light,
+            "sampling_number": 50,
+            "num_kin_sampling": 50,
+            "num_psf_sampling": 50,
+        }
+
+        kwargs_lens_light = [
+            {
+                "R_sersic": 2,
+                "amp": 1,
+                "n_sersic": 2,
+                "center_x": 0,
+                "center_y": 0,
+            }
+        ]
+        lens_light_model_list = ["SERSIC"]
+
+        gamma_in_array = np.linspace(0.1, 2.9, 5)
+        log_m2l_array = np.linspace(0.1, 1, 5)
+
+        rho0_array = np.random.normal(5, 0, 100)
+        r_s_array = np.random.normal(0.1, 0, 100)
+
+        # compute likelihood
+        kin_constraints = KinConstraintsComposite(
+            z_lens=z_lens,
+            z_source=z_source,
+            gamma_in_array=gamma_in_array,
+            log_m2l_array=log_m2l_array,
+            alpha_Rs_array=rho0_array,
+            kappa_s_array=[],
+            r_s_angle_array=r_s_array,
+            theta_E=theta_E,
+            theta_E_error=0.01,
+            gamma=gamma,
+            gamma_error=0.02,
+            r_eff=r_eff,
+            r_eff_error=0.05,
+            sigma_v_measured=[200],
+            sigma_v_error_independent=[10],
+            sigma_v_error_covariant=0,
+            kwargs_aperture=kwargs_aperture,
+            kwargs_seeing=kwargs_seeing,
+            anisotropy_model=anisotropy_model,
+            kwargs_lens_light=kwargs_lens_light,
+            lens_light_model_list=lens_light_model_list,
+            **kwargs_kin_api_settings
+        )
+
+        kin_constraints.draw_lens(no_error=True)
+
+        kwargs_likelihood = kin_constraints.hierarchy_configuration(num_sample_model=5)
+        kwargs_likelihood["normalized"] = False
+        ln_class = LensLikelihood(
+            gamma_in_sampling=True, log_m2l_sampling=True, **kwargs_likelihood
+        )
+        kwargs_kin = {"a_ani": 1}
+        ln_class.lens_log_likelihood(
+            cosmo, kwargs_lens={"gamma_in": 2, "log_m2l": 1}, kwargs_kin=kwargs_kin
+        )
+
+
+class TestKinConstraintsCompositeAlphaRsM2l(object):
+    def setup_method(self):
+        pass
+
+    def test_likelihoodconfiguration_om(self):
+        anisotropy_model = "OM"
+        kwargs_aperture = {
+            "aperture_type": "shell",
+            "r_in": 0,
+            "r_out": 3 / 2.0,
+            "center_ra": 0.0,
+            "center_dec": 0,
+        }
+        kwargs_seeing = {"psf_type": "GAUSSIAN", "fwhm": 1.4}
+
+        # numerical settings (not needed if power-law profiles with Hernquist light distribution is computed)
+        kwargs_numerics_galkin = {
+            "interpol_grid_num": 50,  # numerical interpolation, should converge -> infinity
+            "log_integration": True,
+            # log or linear interpolation of surface brightness and mass models
+            "max_integrate": 100,
+            "min_integrate": 0.001,
+        }  # lower/upper bound of numerical integrals
+
+        # redshift
+        z_lens = 0.5
+        z_source = 1.5
+
+        cosmo = FlatLambdaCDM(H0=70, Om0=0.3)
+        theta_E = 1.0
+        r_eff = 1
+        gamma = 2.1
+
+        kwargs_mge_light = {
+            "grid_spacing": 0.1,
+            "grid_num": 10,
+            "n_comp": 20,
+            "center_x": 0,
+            "center_y": 0,
+        }
+
+        kwargs_kin_api_settings = {
+            "multi_observations": False,
+            "kwargs_numerics_galkin": kwargs_numerics_galkin,
+            "kwargs_mge_light": kwargs_mge_light,
+            "sampling_number": 50,
+            "num_kin_sampling": 50,
+            "num_psf_sampling": 50,
+        }
+
+        kwargs_lens_light = [
+            {
+                "R_sersic": 2,
+                "amp": 1,
+                "n_sersic": 2,
+                "center_x": 0,
+                "center_y": 0,
+            }
+        ]
+        lens_light_model_list = ["SERSIC"]
+
+        gamma_in_array = np.linspace(0.1, 2.9, 5)
+
+        log_m2l_array = np.random.uniform(0.1, 1, 5)
+        rho0_array = np.random.normal(5, 0, 5)
+        r_s_array = np.random.normal(0.1, 0, 5)
+
+        # compute likelihood
+        kin_constraints = KinConstraintsComposite(
+            z_lens=z_lens,
+            z_source=z_source,
+            gamma_in_array=gamma_in_array,
+            log_m2l_array=log_m2l_array,
+            alpha_Rs_array=rho0_array,
+            kappa_s_array=[],
+            r_s_angle_array=r_s_array,
+            theta_E=theta_E,
+            theta_E_error=0.01,
+            gamma=gamma,
+            gamma_error=0.02,
+            r_eff=r_eff,
+            r_eff_error=0.05,
+            sigma_v_measured=[200],
+            sigma_v_error_independent=[10],
+            sigma_v_error_covariant=0,
+            kwargs_aperture=kwargs_aperture,
+            kwargs_seeing=kwargs_seeing,
+            anisotropy_model=anisotropy_model,
+            kwargs_lens_light=kwargs_lens_light,
+            lens_light_model_list=lens_light_model_list,
+            is_m2l_population_level=False,
+            gamma_in_prior_mean=2,
+            gamma_in_prior_std=0.5,
+            **kwargs_kin_api_settings
+        )
+
+        kin_constraints.draw_lens(no_error=True)
+
+        kwargs_likelihood = kin_constraints.hierarchy_configuration(num_sample_model=5)
+        kwargs_likelihood["normalized"] = False
+        ln_class = LensLikelihood(
+            gamma_in_sampling=True, log_m2l_sampling=False, **kwargs_likelihood
+        )
+        kwargs_kin = {"a_ani": 1}
+        ln_class.lens_log_likelihood(
+            cosmo, kwargs_lens={"gamma_in": 2, "log_m2l": 0}, kwargs_kin=kwargs_kin
+        )
+
+
 class TestRaise(unittest.TestCase):
     def test_raise(self):
         with self.assertRaises(ValueError):
@@ -532,6 +753,7 @@ def test_raise(self):
                 z_source=z_source,
                 gamma_in_array=gamma_in_array,
                 log_m2l_array=log_m2l_array,
+                alpha_Rs_array=[],
                 kappa_s_array=rho0_array,
                 r_s_angle_array=r_s_array,
                 theta_E=theta_E,
@@ -609,6 +831,7 @@ def test_raise(self):
                 z_source=z_source,
                 gamma_in_array=gamma_in_array,
                 log_m2l_array=log_m2l_array,
+                alpha_Rs_array=[],
                 kappa_s_array=rho0_array,
                 r_s_angle_array=r_s_array,
                 theta_E=theta_E,
@@ -683,6 +906,7 @@ def test_raise(self):
                 z_source=z_source,
                 gamma_in_array=gamma_in_array,
                 log_m2l_array=log_m2l_array,
+                alpha_Rs_array=[],
                 kappa_s_array=[],
                 r_s_angle_array=[],
                 theta_E=theta_E,
@@ -760,6 +984,7 @@ def test_raise_m2l(self):
                 z_source=z_source,
                 gamma_in_array=gamma_in_array,
                 log_m2l_array=log_m2l_array,
+                alpha_Rs_array=[],
                 kappa_s_array=rho0_array,
                 r_s_angle_array=r_s_array,
                 theta_E=theta_E,
@@ -838,6 +1063,7 @@ def test_raise_m2l(self):
                 z_source=z_source,
                 gamma_in_array=gamma_in_array,
                 log_m2l_array=log_m2l_array,
+                alpha_Rs_array=[],
                 kappa_s_array=rho0_array,
                 r_s_angle_array=r_s_array,
                 theta_E=theta_E,
@@ -913,6 +1139,7 @@ def test_raise_m2l(self):
                 z_source=z_source,
                 gamma_in_array=gamma_in_array,
                 log_m2l_array=log_m2l_array,
+                alpha_Rs_array=[],
                 kappa_s_array=rho0_array,
                 r_s_angle_array=r_s_array,
                 theta_E=theta_E,