Merge 8f7a925 into bc48fa7

hierarc/Data/SNe/README.rst

@@ -0,0 +1,3 @@
+# Pantheon data set
+
+Source: https://github.com/dscolnic/Pantheon/blob/master/Binned_data/lcparam_DS17f.txt

hierarc/Data/SNe/pantheon_binned_lcparam_DS17f.txt

@@ -0,0 +1,41 @@
+#name zcmb zhel dz mb dmb x1 dx1 color dcolor 3rdvar d3rdvar cov_m_s cov_m_c cov_s_c set ra dec biascor
+0 0.014 0.014 0.0 14.57001926 0.0311 0 0 0 0 0 0 0 0 0 0 0 0
+1 0.0194 0.0194 0.0 15.2279603507 0.02645 0 0 0 0 0 0 0 0 0 0 0 0
+2 0.0264 0.0264 0.0 15.934102317 0.0231 0 0 0 0 0 0 0 0 0 0 0 0
+3 0.0329 0.0329 0.0 16.4217829558 0.0235 0 0 0 0 0 0 0 0 0 0 0 0
+4 0.0396 0.0396 0.0 16.836418956 0.02845 0 0 0 0 0 0 0 0 0 0 0 0
+5 0.0475 0.0475 0.0 17.2346439906 0.0334 0 0 0 0 0 0 0 0 0 0 0 0
+6 0.056 0.056 0.0 17.5758459622 0.03215 0 0 0 0 0 0 0 0 0 0 0 0
+7 0.064 0.064 0.0 17.9104543501 0.0439 0 0 0 0 0 0 0 0 0 0 0 0
+8 0.0721 0.0721 0.0 18.1780074589 0.03865 0 0 0 0 0 0 0 0 0 0 0 0
+9 0.0811 0.0811 0.0 18.5543855764 0.04225 0 0 0 0 0 0 0 0 0 0 0 0
+10 0.0889 0.0889 0.0 18.7006554037 0.03635 0 0 0 0 0 0 0 0 0 0 0 0
+11 0.1001 0.1001 0.0 19.0388242428 0.03515 0 0 0 0 0 0 0 0 0 0 0 0
+12 0.1071 0.1071 0.0 19.1720185353 0.03355 0 0 0 0 0 0 0 0 0 0 0 0
+13 0.1195 0.1195 0.0 19.3711295152 0.026 0 0 0 0 0 0 0 0 0 0 0 0
+14 0.1278 0.1278 0.0 19.5555026238 0.02635 0 0 0 0 0 0 0 0 0 0 0 0
+15 0.1396 0.1396 0.0 19.8185930766 0.02335 0 0 0 0 0 0 0 0 0 0 0 0
+16 0.1519 0.1519 0.0 19.9445025196 0.0237 0 0 0 0 0 0 0 0 0 0 0 0
+17 0.1635 0.1635 0.0 20.0642537602 0.02585 0 0 0 0 0 0 0 0 0 0 0 0
+18 0.1778 0.1778 0.0 20.3389008923 0.01965 0 0 0 0 0 0 0 0 0 0 0 0
+19 0.1906 0.1906 0.0 20.5047456137 0.02205 0 0 0 0 0 0 0 0 0 0 0 0
+20 0.2067 0.2067 0.0 20.6854007995 0.0214 0 0 0 0 0 0 0 0 0 0 0 0
+21 0.2216 0.2216 0.0 20.8610322507 0.0233 0 0 0 0 0 0 0 0 0 0 0 0
+22 0.2405 0.2405 0.0 21.0518440051 0.0222 0 0 0 0 0 0 0 0 0 0 0 0
+23 0.2558 0.2558 0.0 21.2037930719 0.02025 0 0 0 0 0 0 0 0 0 0 0 0
+24 0.2762 0.2762 0.0 21.3636137766 0.0221 0 0 0 0 0 0 0 0 0 0 0 0
+25 0.2972 0.2972 0.0 21.5790014334 0.0216 0 0 0 0 0 0 0 0 0 0 0 0
+26 0.3215 0.3215 0.0 21.813277279 0.0214 0 0 0 0 0 0 0 0 0 0 0 0
+27 0.3453 0.3453 0.0 21.9665751024 0.0234 0 0 0 0 0 0 0 0 0 0 0 0
+28 0.3708 0.3708 0.0 22.1428940056 0.02145 0 0 0 0 0 0 0 0 0 0 0 0
+29 0.4049 0.4049 0.0 22.3799192408 0.03225 0 0 0 0 0 0 0 0 0 0 0 0
+30 0.4355 0.4355 0.0 22.5579272526 0.0254 0 0 0 0 0 0 0 0 0 0 0 0
+31 0.4738 0.4738 0.0 22.797376642 0.02935 0 0 0 0 0 0 0 0 0 0 0 0
+32 0.5174 0.5174 0.0 23.0011594448 0.02685 0 0 0 0 0 0 0 0 0 0 0 0
+33 0.5742 0.5742 0.0 23.3004345026 0.0245 0 0 0 0 0 0 0 0 0 0 0 0
+34 0.6299 0.6299 0.0 23.5036280324 0.031 0 0 0 0 0 0 0 0 0 0 0 0
+35 0.724 0.724 0.0 23.8666293039 0.027 0 0 0 0 0 0 0 0 0 0 0 0
+36 0.821 0.821 0.0 24.2445519795 0.0248 0 0 0 0 0 0 0 0 0 0 0 0
+37 0.9511 0.9511 0.0 24.6411323707 0.0276 0 0 0 0 0 0 0 0 0 0 0 0
+38 1.2336 1.2336 0.0 25.3039940191 0.05635 0 0 0 0 0 0 0 0 0 0 0 0
+39 1.6123 1.6123 0.0 25.9259729107 0.0735 0 0 0 0 0 0 0 0 0 0 0 0

hierarc/Likelihood/SneLikelihood/sne_likelihood.py

@@ -0,0 +1,150 @@
+"""
+This is a lightweight version of the COSMOMC/Cobaya sampler: https://github.com/CobayaSampler/cobaya/blob/71b87842d12c6a04eec182c39b6bef1cd9a987af/cobaya/likelihoods/_base_classes/_sn_prototype.py#L287
+It uses the binned Pantheon data: https://github.com/dscolnic/Pantheon/blob/master/Binned_data/lcparam_DS17f.txt
+And computes the cosmographic likelihood.
+The main difference is that this class is compatible with the hierArc cosmology module for evaluating likelihoods.
+This likelihood does NOT include systematics!
+
+.. |br| raw:: html
+   <br />
+.. note::
+   - If you use ``sn.pantheon``, please cite:|br|
+     Scolnic, D. M. et al,
+     `The Complete Light-curve Sample of Spectroscopically
+     Confirmed Type Ia Supernovae from Pan-STARRS1 and
+     Cosmological Constraints from The Combined Pantheon Sample`
+     `(arXiv:1710.00845) <https://arxiv.org/abs/1710.00845>`_
+
+
+:Synopsis: Supernovae likelihood, from CosmoMC's JLA module, for Pantheon and JLA samples.
+:Author: Alex Conley, Marc Betoule, Antony Lewis (see source for more specific authorship)
+
+"""
+__author__ = 'sibirrer'
+
+# Global
+import numpy as np
+import os
+
+# Local
+import hierarc
+
+_twopi = 2 * np.pi
+_SAMPLE_NAME_SUPPORTED = ['Pantheon_binned']
+
+
+class SneLikelihood(object):
+    """
+    Base likelihood class for evaluating Sne likelihoods
+    """
+    def __init__(self, sample_name='Pantheon_binned', pec_z=0.001):
+        """
+
+        :param sample_name: string, name of data sample
+        :param pec_z: float, peculiar velocity in units of redshift (ignored when binned data products are used)
+        """
+        if sample_name not in _SAMPLE_NAME_SUPPORTED:
+            raise ValueError('Sample name %s not supported. Please chose the Sne sample name among %s.'
+                             % (sample_name, _SAMPLE_NAME_SUPPORTED))
+        if sample_name == 'Pantheon_binned':
+            self._data_file = os.path.join(os.path.dirname(hierarc.__file__), 'Data', 'SNe',
+                                           'pantheon_binned_lcparam_DS17f.txt')
+            pec_z = 0
+
+        self._pecz = pec_z
+        cols = None
+
+        self.names = []
+        ix = 0
+        with open(self._data_file, 'r') as f:
+            lines = f.readlines()
+            for line in lines:
+                if '#' in line:
+                    cols = line[1:].split()
+                    for rename, new in zip(
+                            ['mb', 'color', 'x1', '3rdvar', 'd3rdvar',
+                             'cov_m_s', 'cov_m_c', 'cov_s_c'],
+                            ['mag', 'colour', 'stretch', 'third_var',
+                             'dthird_var', 'cov_mag_stretch',
+                             'cov_mag_colour', 'cov_stretch_colour']):
+                        if rename in cols:
+                            cols[cols.index(rename)] = new
+
+                    zeros = np.zeros(len(lines) - 1)
+                    self.set = zeros.copy()
+                    for col in cols:
+                        setattr(self, col, zeros.copy())
+                elif line.strip():
+                    if cols is None: raise ImportError('Data file must have comment header')
+                    vals = line.split()
+                    for i, (col, val) in enumerate(zip(cols, vals)):
+                        if col == 'name':
+                            self.names.append(val)
+                        else:
+                            getattr(self, col)[ix] = np.float64(val)
+                    ix += 1
+        # Check whether required instances are read in
+        assert hasattr(self, 'dz')
+        # TODO: make read-in such that the arguments required are explicitly matched ('zcmb', 'zhel', 'dz', 'mag', 'dmb')
+        # spectroscopic redshift error. ATTENTION! This value =0 for binned data. In this code the value is not used.
+        # Cobaya also does not use it!
+        self.z_var = self.dz ** 2
+        self.mag_var = self.dmb ** 2
+
+        self.nsn = ix
+
+        # jla_prep
+        zfacsq = 25.0 / np.log(10.0) ** 2
+        self.pre_vars = self.mag_var + zfacsq * self._pecz ** 2 * (
+                (1.0 + self.zcmb) / (self.zcmb * (1 + 0.5 * self.zcmb))) ** 2
+
+        self._inv_cov = self._inverse_covariance_matrix()
+
+    def _inverse_covariance_matrix(self):
+        """
+        inverse error covariance matrix. Combines redshift uncertainties (to first order) and magnitude uncertainties
+
+        :return: inverse covariance matrix (2d numpy array)
+        """
+        # here is the option for adding an additional covariance matrix term of the calibration and/or systematic
+        # errors in the evolution of the Sne population
+        invcovmat = np.zeros((len(self.pre_vars), len(self.pre_vars)))
+        invcovmat_diag = invcovmat.diagonal()  # if invcovmat is a matrix, then this is invcovmat.diagonal()
+
+        delta = self.pre_vars
+        np.fill_diagonal(invcovmat, invcovmat_diag + delta)
+        invcov = np.linalg.inv(invcovmat)
+        return invcov
+
+    def logp_lum_dist(self, lum_dists):
+        """
+
+        :param lum_dists: numpy array of luminosity distances to the measured supernovae bins
+         (units do not matter since normalization is subtracted off for the likelihood)
+        :return: log likelihood of the data given the luminosity distances
+        """
+
+        invvars = 1.0 / self.pre_vars
+        wtval = np.sum(invvars)
+        # uncertainty weighted estimated normalization of magnitude (maximum likelihood value)
+        estimated_scriptm = np.sum((self.mag - lum_dists) * invvars) / wtval
+        diffmag = self.mag - lum_dists - estimated_scriptm
+        invcovmat = self._inv_cov
+
+        invvars = invcovmat.dot(diffmag)
+        amarg_A = invvars.dot(diffmag)
+
+        amarg_B = np.sum(invvars)
+        amarg_E = np.sum(invcovmat)
+        chi2 = amarg_A + np.log(amarg_E / _twopi) - amarg_B ** 2 / amarg_E
+        return - chi2 / 2
+
+    def log_likelihood(self, cosmo):
+        """
+
+        :param cosmo: instance of a class to compute angular diameter distances on arrays
+        :return: log likelihood of the data given the specified cosmology
+        """
+        angular_diameter_distances = cosmo.angular_diameter_distance(self.zcmb).value
+        lum_dists = (5 * np.log10((1 + self.zhel) * (1 + self.zcmb) * angular_diameter_distances))
+        return self.logp_lum_dist(lum_dists)

hierarc/Likelihood/cosmo_likelihood.py

@@ -1,6 +1,7 @@
 from hierarc.Likelihood.lens_sample_likelihood import LensSampleLikelihood
 from hierarc.Sampling.ParamManager.param_manager import ParamManager
 from hierarc.Likelihood.cosmo_interp import CosmoInterp
+from hierarc.Likelihood.SneLikelihood.sne_likelihood import SneLikelihood
 import numpy as np
 
 
@@ -9,7 +10,8 @@ class CosmoLikelihood(object):
     this class contains the likelihood function of the Strong lensing analysis
     """
 
-    def __init__(self, kwargs_likelihood_list, cosmology, kwargs_bounds, ppn_sampling=False,
+    def __init__(self, kwargs_likelihood_list, cosmology, kwargs_bounds, sne_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,
                  lambda_ifu_sampling=False, lambda_ifu_distribution='NONE', sigma_v_systematics=False,
@@ -26,6 +28,8 @@ def __init__(self, kwargs_likelihood_list, cosmology, kwargs_bounds, ppn_samplin
         'kwargs_lower_lens', 'kwargs_upper_lens', 'kwargs_fixed_lens',
         'kwargs_lower_kin', 'kwargs_upper_kin', 'kwargs_fixed_kin'
         'kwargs_lower_cosmo', 'kwargs_upper_cosmo', 'kwargs_fixed_cosmo'
+        :param sne_likelihood: (string), optional. Sampling supernovae relative expansion history likelihood, see
+         SneLikelihood module for options
         :param ppn_sampling:post-newtonian parameter sampling
         :param lambda_mst_sampling: bool, if True adds a global mass-sheet transform parameter in the sampling
         :param lambda_mst_distribution: string, defines the distribution function of lambda_mst
@@ -78,6 +82,13 @@ def __init__(self, kwargs_likelihood_list, cosmology, kwargs_bounds, ppn_samplin
         self._num_redshift_interp = num_redshift_interp
         self._cosmo_fixed = cosmo_fixed
         z_max = 0
+        if sne_likelihood is not None:
+            self._sne_likelihood = SneLikelihood(sample_name=sne_likelihood)
+            z_max = np.max(self._sne_likelihood.zcmb)
+            self._sne_evaluate = True
+        else:
+            self._sne_evaluate = False
+
         for kwargs_lens in kwargs_likelihood_list:
             if kwargs_lens['z_source'] > z_max:
                 z_max = kwargs_lens['z_source']
@@ -108,6 +119,8 @@ def likelihood(self, args):
         logL = self._likelihoodLensSample.log_likelihood(cosmo=cosmo, kwargs_lens=kwargs_lens, kwargs_kin=kwargs_kin,
                                                          kwargs_source=kwargs_source)
 
+        if self._sne_evaluate is True:
+            logL += self._sne_likelihood.log_likelihood(cosmo=cosmo)
         if self._prior_add is True:
             logL += self._custom_prior(kwargs_cosmo, kwargs_lens, kwargs_kin)
         return logL

test/test_Likelihood/test_LensLikelihood/test_kin_likelihood.py

@@ -90,6 +90,17 @@ def test_log_likelihood_marg(self):
         logl_average = np.log(l_sum/num_sample)
         npt.assert_almost_equal(logl, logl_average, decimal=1)
 
+    def test_log_likelihood_invert_error(self):
+        """
+        test that in the case of a none-invertable covariance matrix -np.inf is returned by the likelihood statement
+
+        """
+        ifu_likelihood = KinLikelihood(z_lens=0.5, z_source=2., sigma_v_measurement=[100],
+                                       j_model=[1], error_cov_measurement=[0], error_cov_j_sqrt=[0], normalized=True,
+                                       sigma_sys_error_include=True)
+        logl = ifu_likelihood.log_likelihood(ddt=1, dd=1, aniso_scaling=1)
+        assert logl == -np.inf
+
 
 class TestRaise(unittest.TestCase):
 

test/test_Likelihood/test_SneLikelihood/test_sne_likelihood.py

@@ -0,0 +1,39 @@
+from hierarc.Likelihood.SneLikelihood.sne_likelihood import SneLikelihood
+import os
+import unittest
+import pytest
+import numpy.testing as npt
+
+
+class TestSnePantheon(object):
+
+    def setup(self):
+        self.pantheon_likelihood = SneLikelihood(sample_name='Pantheon_binned')
+
+    def test_import_pantheon(self):
+
+        assert os.path.exists(self.pantheon_likelihood._data_file)
+        assert len(self.pantheon_likelihood.zcmb) == 40
+        assert len(self.pantheon_likelihood.zhel) == 40
+        assert len(self.pantheon_likelihood.mag) == 40
+        assert len(self.pantheon_likelihood.mag_var) == 40
+        assert len(self.pantheon_likelihood.z_var) == 40
+
+    def test_log_likelihood(self):
+
+        # cosmo instance
+        from astropy.cosmology import WMAP9 as cosmo
+        logL = self.pantheon_likelihood.log_likelihood(cosmo=cosmo)
+        npt.assert_almost_equal(logL, -29.447965959089437, decimal=4)
+
+
+class TestRaise(unittest.TestCase):
+
+    def test_raise(self):
+
+        with self.assertRaises(ValueError):
+            base = SneLikelihood(sample_name='BAD')
+
+
+if __name__ == '__main__':
+    pytest.main()

test/test_Likelihood/test_cosmo_interp.py

@@ -37,6 +37,22 @@ def test_angular_diameter_distance(self):
         npt.assert_almost_equal(da_interp / da, 1, decimal=3)
         assert da.unit == da_interp.unit
 
+    def test_angular_diameter_distance_array(self):
+        # test for array input
+        z1 = 1.
+        z2 = 2.
+        da_z1 = self.cosmo.angular_diameter_distance(z=[z1])
+        da_z2 = self.cosmo.angular_diameter_distance(z=[z2])
+        da_interp = self.cosmo_interp.angular_diameter_distance(z=[z1, z2])
+        npt.assert_almost_equal(da_interp[0] / da_z1, 1, decimal=3)
+        npt.assert_almost_equal(da_interp[1] / da_z2, 1, decimal=3)
+        assert da_z1.unit == da_interp.unit
+        assert len(da_interp) == 2
+
+        da_z12 = self.cosmo.angular_diameter_distance(z=[z1, z2])
+        npt.assert_almost_equal(da_z12[0] / da_z1, 1, decimal=3)
+        npt.assert_almost_equal(da_z12[1] / da_z2, 1, decimal=3)
+
     def test_angular_diameter_distance_z1z2(self):
         z1 = .3
         z2 = 2.

test/test_Likelihood/test_cosmo_likelihood.py

@@ -115,6 +115,14 @@ def test_cosmo_instance(self):
         npt.assert_almost_equal(dd_astropy, dd_fixed, decimal=1)
         npt.assert_almost_equal(dd_astropy, dd_fixed_interp, decimal=1)
 
+    def test_sne_likelihood_integration(self):
+        cosmoL = CosmoLikelihood([], self.cosmology, self.kwargs_bounds, sne_likelihood='Pantheon_binned',
+                                 interpolate_cosmo=True, num_redshift_interp=100, cosmo_fixed=None)
+        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)
+        assert logl < 0
+
 
 if __name__ == '__main__':
     pytest.main()