Merge pull request #19 from chinyitan/main

Allow 2-parameter lambda_sampling

hierarc/Likelihood/cosmo_likelihood.py

@@ -17,7 +17,8 @@ def __init__(self, kwargs_likelihood_list, cosmology, kwargs_bounds, sne_likelih
                  kwargs_sne_likelihood=None, KDE_likelihood_chain=None, kwargs_kde_likelihood=None,
                  ppn_sampling=False,
                  lambda_mst_sampling=False, lambda_mst_distribution='delta', anisotropy_sampling=False,
-                 kappa_ext_sampling=False, kappa_ext_distribution='NONE', alpha_lambda_sampling=False,
+                 kappa_ext_sampling=False, kappa_ext_distribution='NONE',
+                 alpha_lambda_sampling=False, beta_lambda_sampling=False,
                  lambda_ifu_sampling=False, lambda_ifu_distribution='NONE', sigma_v_systematics=False,
                  sne_apparent_m_sampling=False, sne_distribution='GAUSSIAN', z_apparent_m_anchor=0.1,
                  log_scatter=False,
@@ -44,7 +45,9 @@ def __init__(self, kwargs_likelihood_list, cosmology, kwargs_bounds, sne_likelih
         independently
         :param lambda_ifu_distribution: string, distribution function of the lambda_ifu parameter
         :param alpha_lambda_sampling: bool, if True samples a parameter alpha_lambda, which scales lambda_mst linearly
-         according to a predefined quantity of the lens
+         according to the lens posterior kwargs 'lambda_scaling_property'
+        :param beta_lambda_sampling: bool, if True samples a parameter beta_lambda, which scales lambda_mst linearly
+         according to the lens posterior kwargs 'lambda_scaling_property_beta'
         :param kappa_ext_sampling: bool, if True samples a global external convergence parameter
         :param kappa_ext_distribution: string, distribution function of the kappa_ext parameter
         :param anisotropy_sampling: bool, if True adds a global stellar anisotropy parameter that alters the single lens
@@ -75,6 +78,7 @@ def __init__(self, kwargs_likelihood_list, cosmology, kwargs_bounds, sne_likelih
                                   lambda_ifu_sampling=lambda_ifu_sampling,
                                   lambda_ifu_distribution=lambda_ifu_distribution,
                                   alpha_lambda_sampling=alpha_lambda_sampling,
+                                  beta_lambda_sampling=beta_lambda_sampling,
                                   sne_apparent_m_sampling=sne_apparent_m_sampling,
                                   sne_distribution=sne_distribution, z_apparent_m_anchor=z_apparent_m_anchor,
                                   sigma_v_systematics=sigma_v_systematics,

hierarc/Likelihood/hierarchy_likelihood.py

@@ -13,7 +13,8 @@ class LensLikelihood(TransformedCosmography, LensLikelihoodBase, AnisotropyScali
     def __init__(self, z_lens, z_source, name='name', likelihood_type='TDKin', anisotropy_model='NONE',
                  ani_param_array=None, ani_scaling_array_list=None, ani_scaling_array=None,
                  num_distribution_draws=50, kappa_ext_bias=False, kappa_pdf=None, kappa_bin_edges=None, mst_ifu=False,
-                 lambda_scaling_property=0, normalized=False, kwargs_lens_properties=None, **kwargs_likelihood):
+                 lambda_scaling_property=0,lambda_scaling_property_beta=0,
+                 normalized=False, kwargs_lens_properties=None, **kwargs_likelihood):
         """
 
         :param z_lens: lens redshift
@@ -37,6 +38,8 @@ def __init__(self, z_lens, z_source, name='name', likelihood_type='TDKin', aniso
          in sampling this lens.
         :param lambda_scaling_property: float (optional), scaling of
          lambda_mst = lambda_mst_global + alpha * lambda_scaling_property
+        :param lambda_scaling_property_beta: float (optional), scaling of
+         lambda_mst = lambda_mst_global + beta * lambda_scaling_property_beta
         :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
@@ -60,6 +63,7 @@ def __init__(self, z_lens, z_source, name='name', likelihood_type='TDKin', aniso
         else:
             self._draw_kappa = False
         self._lambda_scaling_property = lambda_scaling_property
+        self._lambda_scaling_property_beta = lambda_scaling_property_beta
 
     def lens_log_likelihood(self, cosmo, kwargs_lens=None, kwargs_kin=None, kwargs_source=None):
         """
@@ -199,7 +203,7 @@ def check_dist(self, kwargs_lens, kwargs_kin, kwargs_source):
         return False
 
     def draw_lens(self, lambda_mst=1, lambda_mst_sigma=0, kappa_ext=0, kappa_ext_sigma=0, gamma_ppn=1, lambda_ifu=1,
-                  lambda_ifu_sigma=0, alpha_lambda=0):
+                  lambda_ifu_sigma=0, alpha_lambda=0, beta_lambda=0):
         """
         draws a realization of a specific model from the hyper-parameter distribution
 
@@ -212,13 +216,17 @@ def draw_lens(self, lambda_mst=1, lambda_mst_sigma=0, kappa_ext=0, kappa_ext_sig
         :param lambda_ifu_sigma: secondary lambda_mst_sigma parameter for subset of lenses specified for
         :param alpha_lambda: float, linear slope of the lambda_int scaling relation with lens quantity
          self._lambda_scaling_property
+        :param beta_lambda: float, a second linear slope of the lambda_int scaling relation with lens quantity
+         self._lambda_scaling_property_beta
         :return: draw from the distributions
         """
         if self._mst_ifu is True:
-            lambda_lens = lambda_ifu + alpha_lambda * self._lambda_scaling_property
+            lambda_lens = lambda_ifu + alpha_lambda * self._lambda_scaling_property \
+                          + beta_lambda * self._lambda_scaling_property_beta
             lambda_mst_draw = np.random.normal(lambda_lens, lambda_ifu_sigma)
         else:
-            lambda_lens = lambda_mst + alpha_lambda * self._lambda_scaling_property
+            lambda_lens = lambda_mst + alpha_lambda * self._lambda_scaling_property \
+                          + beta_lambda * self._lambda_scaling_property_beta
             lambda_mst_draw = np.random.normal(lambda_lens, lambda_mst_sigma)
         if self._draw_kappa is True:
             kappa_ext_draw = self._kappa_dist.draw_one

hierarc/Sampling/ParamManager/lens_param.py

@@ -7,7 +7,7 @@ class LensParam(object):
     """
     def __init__(self, lambda_mst_sampling=False, lambda_mst_distribution='NONE', kappa_ext_sampling=False,
                  kappa_ext_distribution='NONE', lambda_ifu_sampling=False, lambda_ifu_distribution='NONE',
-                 alpha_lambda_sampling=False, kwargs_fixed=None, log_scatter=False):
+                 alpha_lambda_sampling=False, beta_lambda_sampling=False, kwargs_fixed=None, log_scatter=False):
         """
 
         :param lambda_mst_sampling: bool, if True adds a global mass-sheet transform parameter in the sampling
@@ -19,6 +19,8 @@ def __init__(self, lambda_mst_sampling=False, lambda_mst_distribution='NONE', ka
         :param lambda_ifu_distribution: string, distribution function of the lambda_ifu parameter
         :param alpha_lambda_sampling: bool, if True samples a parameter alpha_lambda, which scales lambda_mst linearly
          according to a predefined quantity of the lens
+        :param beta_lambda_sampling: bool, if True samples a parameter beta_lambda, which scales lambda_mst linearly
+         according to a predefined quantity of the lens
         :param log_scatter: boolean, if True, samples the Gaussian scatter amplitude in log space (and thus flat prior in log)
         :param kwargs_fixed: keyword arguments that are held fixed through the sampling
         """
@@ -29,6 +31,7 @@ def __init__(self, lambda_mst_sampling=False, lambda_mst_distribution='NONE', ka
         self._kappa_ext_sampling = kappa_ext_sampling
         self._kappa_ext_distribution = kappa_ext_distribution
         self._alpha_lambda_sampling = alpha_lambda_sampling
+        self._beta_lambda_sampling = beta_lambda_sampling
         self._log_scatter = log_scatter
         if kwargs_fixed is None:
             kwargs_fixed = {}
@@ -89,6 +92,12 @@ def param_list(self, latex_style=False):
                     list.append(r'$\alpha_{\lambda}$')
                 else:
                     list.append('alpha_lambda')
+        if self._beta_lambda_sampling is True:
+            if 'beta_lambda' not in self._kwargs_fixed:
+                if latex_style is True:
+                    list.append(r'$\beta_{\lambda}$')
+                else:
+                    list.append('beta_lambda')
         return list
 
     def args2kwargs(self, args, i=0):
@@ -146,6 +155,12 @@ def args2kwargs(self, args, i=0):
             else:
                 kwargs['alpha_lambda'] = args[i]
                 i += 1
+        if self._beta_lambda_sampling is True:
+            if 'beta_lambda' in self._kwargs_fixed:
+                kwargs['beta_lambda'] = self._kwargs_fixed['beta_lambda']
+            else:
+                kwargs['beta_lambda'] = args[i]
+                i += 1
         return kwargs, i
 
     def kwargs2args(self, kwargs):
@@ -182,4 +197,7 @@ def kwargs2args(self, kwargs):
         if self._alpha_lambda_sampling is True:
             if 'alpha_lambda' not in self._kwargs_fixed:
                 args.append(kwargs['alpha_lambda'])
+        if self._beta_lambda_sampling is True:
+            if 'beta_lambda' not in self._kwargs_fixed:
+                args.append(kwargs['beta_lambda'])
         return args

hierarc/Sampling/ParamManager/param_manager.py

@@ -11,7 +11,8 @@ class for managing the parameters involved
     def __init__(self, cosmology, ppn_sampling=False, lambda_mst_sampling=False, lambda_mst_distribution='NONE',
                  anisotropy_sampling=False, anisotropy_model='OM', anisotropy_distribution='NONE',
                  kappa_ext_sampling=False, kappa_ext_distribution='NONE', lambda_ifu_sampling=False,
-                 lambda_ifu_distribution='NONE', alpha_lambda_sampling=False, sigma_v_systematics=False,
+                 lambda_ifu_distribution='NONE', alpha_lambda_sampling=False, beta_lambda_sampling=False,
+                 sigma_v_systematics=False,
                  sne_apparent_m_sampling=False, sne_distribution='GAUSSIAN', z_apparent_m_anchor=0.1,
                  log_scatter=False,
                  kwargs_lower_cosmo=None, kwargs_upper_cosmo=None,
@@ -28,6 +29,8 @@ def __init__(self, cosmology, ppn_sampling=False, lambda_mst_sampling=False, lam
         independently
         :param alpha_lambda_sampling: bool, if True samples a parameter alpha_lambda, which scales lambda_mst linearly
          according to a predefined quantity of the lens
+        :param beta_lambda_sampling: bool, if True samples a parameter beta_lambda, which scales lambda_mst linearly
+         according to a predefined quantity of the lens
         :param lambda_ifu_distribution: string, distribution function of the lambda_ifu parameter
         :param anisotropy_sampling: bool, if True adds a global stellar anisotropy parameter that alters the single lens
         kinematic prediction
@@ -52,6 +55,7 @@ def __init__(self, cosmology, ppn_sampling=False, lambda_mst_sampling=False, lam
                                      kappa_ext_sampling=kappa_ext_sampling,
                                      kappa_ext_distribution=kappa_ext_distribution,
                                      alpha_lambda_sampling=alpha_lambda_sampling,
+                                     beta_lambda_sampling=beta_lambda_sampling,
                                      log_scatter=log_scatter,
                                      kwargs_fixed=kwargs_fixed_lens)
         self._source_param = SourceParam(sne_apparent_m_sampling=sne_apparent_m_sampling,

test/test_Sampling/test_ParamManager/test_lens_param.py

@@ -7,28 +7,28 @@ class TestLensParam(object):
     def setup(self):
         self._param = LensParam(lambda_mst_sampling=True, lambda_mst_distribution='GAUSSIAN', kappa_ext_sampling=True,
                  kappa_ext_distribution='GAUSSIAN', lambda_ifu_sampling=True, lambda_ifu_distribution='GAUSSIAN',
-                                alpha_lambda_sampling=True, kwargs_fixed={})
+                                alpha_lambda_sampling=True, beta_lambda_sampling=True,  kwargs_fixed={})
 
         kwargs_fixed = {'lambda_mst': 1, 'lambda_mst_sigma': 0.1, 'lambda_ifu': 1.1, 'lambda_ifu_sigma': 0.2,
-                        'kappa_ext': 0.01, 'kappa_ext_sigma': 0.03, 'alpha_lambda': 0}
+                        'kappa_ext': 0.01, 'kappa_ext_sigma': 0.03, 'alpha_lambda': 0, 'beta_lambda': 0 }
         self._param_fixed = LensParam(lambda_mst_sampling=True, lambda_mst_distribution='GAUSSIAN', kappa_ext_sampling=True,
                                 kappa_ext_distribution='GAUSSIAN', lambda_ifu_sampling=True,
-                                lambda_ifu_distribution='GAUSSIAN', alpha_lambda_sampling=True,
+                                lambda_ifu_distribution='GAUSSIAN', alpha_lambda_sampling=True, beta_lambda_sampling=True,
                                 kwargs_fixed=kwargs_fixed)
         self._param_log_scatter = LensParam(lambda_mst_sampling=True, lambda_mst_distribution='GAUSSIAN', kappa_ext_sampling=True,
                  kappa_ext_distribution='GAUSSIAN', lambda_ifu_sampling=True, lambda_ifu_distribution='GAUSSIAN',
-                                alpha_lambda_sampling=True, log_scatter=True, kwargs_fixed={})
+                                alpha_lambda_sampling=True, beta_lambda_sampling=True, log_scatter=True, kwargs_fixed={})
 
     def test_param_list(self):
         param_list = self._param.param_list(latex_style=False)
-        assert len(param_list) == 7
+        assert len(param_list) == 8
         param_list = self._param.param_list(latex_style=True)
-        assert len(param_list) == 7
+        assert len(param_list) == 8
 
         param_list = self._param_log_scatter.param_list(latex_style=False)
-        assert len(param_list) == 7
+        assert len(param_list) == 8
         param_list = self._param_log_scatter.param_list(latex_style=True)
-        assert len(param_list) == 7
+        assert len(param_list) == 8
 
         param_list = self._param_fixed.param_list(latex_style=False)
         assert len(param_list) == 0
@@ -37,7 +37,7 @@ def test_param_list(self):
 
     def test_args2kwargs(self):
         kwargs = {'lambda_mst': 1.1, 'lambda_mst_sigma': 0.1, 'lambda_ifu': 1.1, 'lambda_ifu_sigma': 0.2,
-                        'kappa_ext': 0.01, 'kappa_ext_sigma': 0.03, 'alpha_lambda': 0.1}
+                        'kappa_ext': 0.01, 'kappa_ext_sigma': 0.03, 'alpha_lambda': 0.1, 'beta_lambda': 0.1}
         args = self._param.kwargs2args(kwargs)
         kwargs_new, i = self._param.args2kwargs(args, i=0)
         args_new = self._param.kwargs2args(kwargs_new)