Merge pull request #16 from sibirrer/sne_custom_likelihood

Sne custom likelihood

hierarc/Likelihood/SneLikelihood/sne_likelihood.py

@@ -1,175 +1,35 @@
-"""
-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', 'Pantheon']
-_PATH_2_DATA = os.path.join(os.path.dirname(hierarc.__file__), 'Data', 'SNe')
+from hierarc.Likelihood.SneLikelihood.sne_likelihood_from_file import SneLikelihoodFromFile
+from hierarc.Likelihood.SneLikelihood.sne_likelihood_custom import CustomSneLikelihood
 
 
 class SneLikelihood(object):
     """
-    Base likelihood class for evaluating Sne likelihoods
+    Supernovae likelihood
+    This class supports custom likelihoods as well as likelihoods from the Pantheon sample from file
     """
-    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(_PATH_2_DATA, 'pantheon_binned_lcparam_DS17f.txt')
-            self._cov_file = None
-            pec_z = 0
-
-        if sample_name == 'Pantheon':
-            self._data_file = os.path.join(_PATH_2_DATA, 'pantheon_lcparam_full_long_zhel.txt')
-            self._cov_file = os.path.join(_PATH_2_DATA, 'pantheon_sys_full_long.txt')
-        self._pec_z = 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
-        # variance in the bolometric magnitude distribution of the same for the redshift and type of the SNe
-        self.mag_var = self.dmb ** 2
-
-        self.nsn = ix
-        self._cov = self._read_covmat(self._cov_file)
-
-        # jla_prep
-        zfacsq = 25.0 / np.log(10.0) ** 2
-        # adding peculiar redshift uncertainties to be added to the diagonal variance elements of the covariance matrix
-        self.pre_vars = self.mag_var + zfacsq * self._pec_z ** 2 * (
-                (1.0 + self.zcmb) / (self.zcmb * (1 + 0.5 * self.zcmb))) ** 2
-
-        self._inv_cov = self._inverse_covariance_matrix()
-
-    def _read_covmat(self, filename):
-        """
-        reads in covariance matrix file and returns it as a numpy matrix
-
-        :param filename: string, absolute path of covariance matrix file
-        :return: nxn covariance matrix
-        """
-        if filename is None:
-            return np.zeros((self.nsn, self.nsn))
-        cov = np.loadtxt(filename)
-        if np.isscalar(cov[0]) and cov[0] ** 2 + 1 == len(cov):
-            cov = cov[1:]
-        return cov.reshape((self.nsn, self.nsn))
-
-    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)
+    def __init__(self, sample_name='CUSTOM', **kwargs_sne_likelihood):
         """
-        # 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 = self._cov
-        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, estimated_scriptm=None):
+        :param sample_name: string, either 'CUSTOM' or a specific name supported by SneLikelihoodFromFile() class
+        :param kwargs_sne_likelihood: keyword arguments to initiate likelihood class
         """
+        if sample_name == 'CUSTOM':
+            self._likelihood = CustomSneLikelihood(**kwargs_sne_likelihood)
+        else:
+            self._likelihood = SneLikelihoodFromFile(sample_name=sample_name, **kwargs_sne_likelihood)
+        self.zhel = self._likelihood.zhel
+        self.zcmb = self._likelihood.zcmb
 
-        :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)
-        :param estimated_scriptm: mean magnitude at lum_dist=0 (optional)
-        :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)
-
-        if estimated_scriptm is None:
-            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, apparent_m_z=None, z_anchor=0.1):
+    def log_likelihood(self, cosmo, apparent_m_z=None, sigma_m_z=None, z_anchor=0.1):
         """
 
         :param cosmo: instance of a class to compute angular diameter distances on arrays
         :param apparent_m_z: mean apparent magnitude of SN Ia at z=z_anchor (optional)
         :param z_anchor: redshift where definition of apparent_m_z is set (only applicable when apparent_m_z != None)
+        :param sigma_m_z: 1-sigma scatter in magnitude in the intrinsic SNe brightness distribution not accounted-for
+         by the covariance matrix
         :return: log likelihood of the data given the specified cosmology
         """
         angular_diameter_distances = cosmo.angular_diameter_distance(self.zcmb).value
@@ -178,4 +38,4 @@ def log_likelihood(self, cosmo, apparent_m_z=None, z_anchor=0.1):
         ang_dist_anchor = cosmo.angular_diameter_distance(z_anchor).value
         lum_dist_anchor = (5 * np.log10((1 + z_anchor) * (1 + z_anchor) * ang_dist_anchor))
 
-        return self.logp_lum_dist(lum_dists - lum_dist_anchor, apparent_m_z)
+        return self._likelihood.log_likelihood_lum_dist(lum_dists - lum_dist_anchor, apparent_m_z, sigma_m_z)

hierarc/Likelihood/SneLikelihood/sne_likelihood_custom.py

@@ -0,0 +1,70 @@
+import numpy as np
+
+_twopi = 2 * np.pi
+
+
+class CustomSneLikelihood(object):
+    """
+    class method for an arbitrary apparent magnitude likelihood of a Sne sample where the error and systematic
+    covariance matrix is described in astronomical magnitude space
+    """
+
+    def __init__(self, mag_mean, cov_mag, zhel, zcmb):
+        """
+
+        :param mag_mean: array of mean astronomical magnitudes of the sample of Sne
+        :param cov_mag: error covariance matrix of the magnitudes to result in relative distance moduli
+         (including measurement and systematic uncertainties
+        :param zhel: array, redshift of the exploding shell
+        :param zcmb: array, CMB-corrected redshift of the Sne
+        """
+        self.zhel = zhel
+        self.zcmb = zcmb
+        self.mag = mag_mean
+        self._cov_mag = cov_mag
+        self._inv_cov_mag_input = np.linalg.inv(cov_mag)
+        self.num_sne = len(mag_mean)
+
+    def log_likelihood_lum_dist(self, lum_dists, estimated_scriptm=None, sigma_m_z=None):
+        """
+
+        :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)
+        :param estimated_scriptm: mean magnitude at lum_dist=0 (optional)
+        :param sigma_m_z: 1-sigma scatter in magnitude in the intrinsic SNe brightness distribution not accounted-for
+         by the covariance matrix
+        :return: log likelihood of the data given the luminosity distances
+        """
+        cov_mag, inv_cov = self._inverse_covariance_matrix(sigma_m_z)
+        pre_vars = cov_mag.diagonal()
+        invvars = 1.0 / pre_vars
+        wtval = np.sum(invvars)
+
+        # uncertainty weighted estimated normalization of magnitude (maximum likelihood value)
+        if estimated_scriptm is None:
+            estimated_scriptm = np.sum((self.mag - lum_dists) * invvars) / wtval
+        diffmag = self.mag - lum_dists - estimated_scriptm
+
+        lnlikelihood = -diffmag.dot(inv_cov.dot(diffmag)) / 2.
+        sign_det, lndet = np.linalg.slogdet(cov_mag)
+        lnlikelihood -= 1 / 2. * (self.num_sne * np.log(2 * np.pi) + lndet)
+        return lnlikelihood
+
+    def _inverse_covariance_matrix(self, sigma_m_z=None):
+        """
+        inverse error covariance matrix. Combines redshift uncertainties (to first order) and magnitude uncertainties
+        as well as intrinsic scatter uncertainties
+
+        :param sigma_m_z: float, 1-sigma additional intrinsic magnitude uncertainty of the distribution, not
+        accounted-for in the original covariance matrix
+        :return: covariance matrix, 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
+        if sigma_m_z is None:
+            return self._cov_mag, self._inv_cov_mag_input
+        # cov_mag_diag = self._cov_mag.diagonal()
+        cov_mag = self._cov_mag + np.diag(np.ones(self.num_sne) * sigma_m_z**2)
+        # np.fill_diagonal(self._cov_mag, cov_mag_diag + np.ones(self.num_sne) * sigma_m_z**2)
+        invcov = np.linalg.inv(cov_mag)
+        return cov_mag, invcov