Merge pull request #43 from sibirrer/dsp_msd

verbose feature added in likelihood

hierarc/Likelihood/cosmo_likelihood.py

@@ -104,17 +104,24 @@ def __init__(
                     z_max = kwargs_lens["z_source2"]
         self._z_max = z_max
 
-    def likelihood(self, args, kwargs_cosmo_interp=None):
+    def likelihood(self, args, kwargs_cosmo_interp=None, verbose=False):
         """
 
         :param args: list of sampled parameters
         :param kwargs_cosmo_interp: interpolated angular diameter distances with
          'ang_diameter_distances' and 'redshifts', and optionally 'ok' and 'K' in none-flat scenarios
         :type kwargs_cosmo_interp: dict
+        :param verbose: If true, prints intermediate outputs of likelihood calculation
+        :type verbose: bool
         :return: log likelihood of the combined lenses
         """
         for i in range(0, len(args)):
             if args[i] < self._lower_limit[i] or args[i] > self._upper_limit[i]:
+                if verbose:
+                    print(
+                        "Parameter %i with value %s out of range [%s, %s]"
+                        % (i, args[i], self._lower_limit[i], self._upper_limit[i])
+                    )
                 return -np.inf
 
         kwargs_cosmo, kwargs_lens, kwargs_kin, kwargs_source, kwargs_los = (
@@ -135,6 +142,8 @@ def likelihood(self, args, kwargs_cosmo_interp=None):
                     return -np.inf
             # make sure that Omega_DE is not negative...
             if 1.0 - om - ok <= 0:
+                if verbose:
+                    print("curvature unphysical with 1-Om - Ok <= 0")
                 return -np.inf
         if kwargs_cosmo_interp is not None:
             kwargs_cosmo = {**kwargs_cosmo, **kwargs_cosmo_interp}
@@ -145,6 +154,7 @@ def likelihood(self, args, kwargs_cosmo_interp=None):
             kwargs_kin=kwargs_kin,
             kwargs_source=kwargs_source,
             kwargs_los=kwargs_los,
+            verbose=verbose,
         )
 
         if self._sne_evaluate is True:
@@ -176,7 +186,6 @@ def cosmo_instance(self, kwargs_cosmo):
         :param kwargs_cosmo: cosmology parameter keyword argument list
         :return: ~astropy.cosmology (or equivalent interpolation scheme class)
         """
-        print(kwargs_cosmo, "test kwargs_cosmo")
         if "ang_diameter_distances" in kwargs_cosmo and "redshifts" in kwargs_cosmo:
             # in that case we use directly the interpolation mode to approximate angular diameter distances
             cosmo = CosmoInterp(

hierarc/Likelihood/hierarchy_likelihood.py

@@ -175,16 +175,19 @@ def lens_log_likelihood(
         kwargs_kin=None,
         kwargs_source=None,
         kwargs_los=None,
+        verbose=False,
     ):
         """Log likelihood of the data of a lens given a model (defined with hyper-
         parameters) and cosmology.
 
         :param cosmo: astropy.cosmology instance
-        :param kwargs_lens: keywords of the hyper parameters of the lens model
-        :param kwargs_kin: keyword arguments of the kinematic model hyper parameters
+        :param kwargs_lens: keywords of the hyperparameters of the lens model
+        :param kwargs_kin: keyword arguments of the kinematic model hyperparameters
         :param kwargs_source: keyword argument of the source model (such as SNe)
         :param kwargs_los: list of keyword arguments of global line of sight
             distributions
+        :param verbose: If true, prints intermediate outputs of likelihood calculation
+        :type verbose: bool
         :return: log likelihood of the data given the model
         """
 
@@ -209,6 +212,8 @@ def lens_log_likelihood(
             kwargs_los=kwargs_los,
             cosmo=cosmo,
         )
+        if verbose:
+            print("log likelihood of lens %s = %s" % (self._name, a))
         return np.nan_to_num(a)
 
     def hyper_param_likelihood(

hierarc/Likelihood/lens_sample_likelihood.py

@@ -53,6 +53,7 @@ def log_likelihood(
         kwargs_kin=None,
         kwargs_source=None,
         kwargs_los=None,
+        verbose=False,
     ):
         """
 
@@ -61,6 +62,8 @@ def log_likelihood(
         :param kwargs_kin: keyword arguments of the kinematic model
         :param kwargs_source: keyword argument of the source model (such as SNe)
         :param kwargs_los: line of sight keyword argument list
+        :param verbose: If true, prints intermediate outputs of likelihood calculation
+        :type verbose: bool
         :return: log likelihood of the combined lenses
         """
         log_likelihood = 0
@@ -71,6 +74,7 @@ def log_likelihood(
                 kwargs_kin=kwargs_kin,
                 kwargs_source=kwargs_source,
                 kwargs_los=kwargs_los,
+                verbose=verbose,
             )
         return log_likelihood
 

test/test_Likelihood/test_cosmo_likelihood.py

@@ -99,29 +99,29 @@ def custom_prior(
 
         kwargs_cosmo = {"h0": self.H0_true, "om": self.omega_m_true, "ok": 0}
         args = cosmoL.param.kwargs2args(kwargs_cosmo=kwargs_cosmo)
-        logl = cosmoL.likelihood(args=args)
-        logl_prior = cosmoL_prior.likelihood(args=args)
+        logl = cosmoL.likelihood(args=args, verbose=True)
+        logl_prior = cosmoL_prior.likelihood(args=args, verbose=True)
         npt.assert_almost_equal(logl - logl_prior, 1, decimal=8)
 
         kwargs = {"h0": self.H0_true * 0.99, "om": self.omega_m_true, "ok": 0}
         args = cosmoL.param.kwargs2args(kwargs_cosmo=kwargs)
-        logl_sigma = cosmoL.likelihood(args=args)
+        logl_sigma = cosmoL.likelihood(args=args, verbose=True)
         npt.assert_almost_equal(logl - logl_sigma, 0.5, decimal=2)
 
         kwargs = {"h0": 100, "om": 1.0, "ok": 0.1}
         args = cosmoL.param.kwargs2args(kwargs)
-        logl = cosmoL.likelihood(args=args)
+        logl = cosmoL.likelihood(args=args, verbose=True)
         assert logl == -np.inf
 
         kwargs = {"h0": 100, "om": 0.1, "ok": -0.6}
         args = cosmoL.param.kwargs2args(kwargs)
-        logl = cosmoL.likelihood(args=args)
+        logl = cosmoL.likelihood(args=args, verbose=True)
         assert logl == -np.inf
 
         # outside the prior limit
         kwargs = {"h0": 1000, "om": 0.3, "ok": -0.1}
         args = cosmoL.param.kwargs2args(kwargs)
-        logl = cosmoL.likelihood(args=args)
+        logl = cosmoL.likelihood(args=args, verbose=True)
         assert logl == -np.inf
 
     def test_cosmo_instance(self):
@@ -232,7 +232,7 @@ def test_sne_likelihood_integration(self):
         )
         kwargs_cosmo = {"h0": self.H0_true, "om": self.omega_m_true, "ok": 0}
         args = cosmoL.param.kwargs2args(kwargs_cosmo=kwargs_cosmo)
-        logl = cosmoL.likelihood(args=args)
+        logl = cosmoL.likelihood(args=args, verbose=True)
         assert logl < 0
 
     def test_kde_likelihood_integration(self):
@@ -253,7 +253,7 @@ def test_kde_likelihood_integration(self):
         )
         kwargs_cosmo = {"h0": self.H0_true, "om": self.omega_m_true, "ok": 0}
         args = cosmoL.param.kwargs2args(kwargs_cosmo=kwargs_cosmo)
-        logl = cosmoL.likelihood(args=args)
+        logl = cosmoL.likelihood(args=args, verbose=True)
         assert logl < 0