diff --git a/Notebooks/Validation_cosmo.ipynb b/Notebooks/Validation_cosmo.ipynb index 0d016dd..84c1b18 100644 --- a/Notebooks/Validation_cosmo.ipynb +++ b/Notebooks/Validation_cosmo.ipynb @@ -24,19 +24,10 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 3, "id": "a53c96d1-5fb0-42ed-b5f0-4f0e1d88f856", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "The autoreload extension is already loaded. To reload it, use:\n", - " %reload_ext autoreload\n" - ] - } - ], + "outputs": [], "source": [ "%load_ext autoreload\n", "%autoreload 2\n", @@ -58,7 +49,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 6, "id": "b53aa1cb-b55c-4ef2-b814-ce68fd25098f", "metadata": {}, "outputs": [], @@ -355,10 +346,123 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 67, "id": "14d103a7", "metadata": {}, "outputs": [], + "source": [ + "lcdm = Cosmology('LCDM',['H0','Omega_m'])\n", + "wcdm = Cosmology('wCDM',['H0','Omega_m','w0'])" + ] + }, + { + "cell_type": "code", + "execution_count": 73, + "id": "44681cac", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0. , 3563.00252552, 5480.86846219, 6652.2955333 ,\n", + " 7453.46122345, 8043.88164985, 8501.73372259, 8870.0167495 ,\n", + " 9174.50291195, 9431.67380384])" + ] + }, + "execution_count": 73, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "lcdm._transverse_distance_([70, 0.3], zlist)" + ] + }, + { + "cell_type": "code", + "execution_count": 87, + "id": "53990820", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0. , 3756.12021805, 5743.63627523, 6933.36832763,\n", + " 7740.79000998, 8333.81846036, 8792.9217005 , 9161.87050455,\n", + " 9466.74024632, 9724.14483018])" + ] + }, + "execution_count": 87, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "wcdm._transverse_distance_( [70, 0.3, -1.3], zlist)" + ] + }, + { + "cell_type": "code", + "execution_count": 78, + "id": "b4612f32", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0. , 3756.12021805, 5743.63627523, 6933.36832763,\n", + " 7740.79000998, 8333.81846036, 8792.9217005 , 9161.87050455,\n", + " 9466.74024632, 9724.14483018])" + ] + }, + "execution_count": 78, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "wcdm._transverse_distance_( [70, 0.3, -1.3], zlist)" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "id": "b0ea5a83", + "metadata": {}, + "outputs": [], + "source": [ + "klcdm = Cosmology('kLCDM',['H0','Omega_m','Omega_k'])" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "id": "407eb913", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([4.28157123e+00, 3.49909018e+03, 5.42973331e+03, 6.66479343e+03,\n", + " 7.54168971e+03, 8.20704907e+03, 8.73514182e+03, 9.16806810e+03,\n", + " 9.53172140e+03, 9.84303874e+03])" + ] + }, + "execution_count": 36, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "klcdm.transverse_distance( [70, 0.3, 0.1], zlist)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b4d2effa", + "metadata": {}, + "outputs": [], "source": [] } ], diff --git a/marcia/__init__.py b/marcia/__init__.py index 3082767..ffbeafa 100644 --- a/marcia/__init__.py +++ b/marcia/__init__.py @@ -6,4 +6,4 @@ from .likelihood import Likelihood from .sampler import Sampler from .kernel import * -from .likelihood_GP import Likelihood_GP \ No newline at end of file +from .likelihood import Likelihood_GP \ No newline at end of file diff --git a/marcia/backend/cosmology.py b/marcia/backend/cosmology.py index d1657a8..7a3a0ea 100644 --- a/marcia/backend/cosmology.py +++ b/marcia/backend/cosmology.py @@ -5,6 +5,7 @@ @njit([f8[:](f8,f8,f8,f8,f8[:],f8[:])]) def hubble_rate(H0,Omega_m,Omega_b,Omega_k,de,z): Omega_r = 4.18343*10**-5./(H0/100.)**2. + #print('backend_H_r',(1. - Omega_m - Omega_k - Omega_b - Omega_r)*de) E2 = Omega_r*((1. + z)**4) + Omega_m*((1. + z)**3) + Omega_b*((1. + z)**3) + Omega_k*((1. + z)**2) + (1. - Omega_m - Omega_k - Omega_b - Omega_r)*de Hofz = H0*np.sqrt(E2) return Hofz @@ -44,14 +45,16 @@ def z_inp(z): def interpolate(z_inp,z,func): return np.interp(z,z_inp,func) -@njit -def Ly(y,t,H0,Omega_m,Omega_b,Omega_k,de): - return inv_hubble_rate(H0,Omega_m,Omega_b,Omega_k,de,np.array([t]))[0] -def transverse_distance(H0,Omega_m,Omega_b,Omega_k,de,z,clight): - y=odeint(Ly,0.0,z,args=(H0,Omega_m,Omega_b,Omega_k,de)) - return clight* y[:,0] +@njit +def transverse_distance(H0,clight,Omega_k,y): + if Omega_k > 0.0: + return clight/(np.sqrt(Omega_k)*H0)*np.sinh(np.sqrt(Omega_k) * H0*y) + elif Omega_k < 0.0: + return clight/(np.sqrt(np.abs(Omega_k))*H0)*np.sin((np.sqrt(np.abs(Omega_k)))* H0* y) + elif Omega_k == 0.0: + return clight* y @njit def distance_modulus(Mb,z2,d): diff --git a/marcia/cosmology.py b/marcia/cosmology.py index 4acbc00..82dd931 100644 --- a/marcia/cosmology.py +++ b/marcia/cosmology.py @@ -193,6 +193,7 @@ def __check_mandatory_parameters__(self): def dark_energy_f(self, parameters, z): p = self.param(parameters) + #print('dark_energy_f',p.w0) return cbackend.dark_energy_f_wCDM(p.w0, p.wa, z) def dark_energy_w(self,parameters, z): @@ -201,19 +202,14 @@ def dark_energy_w(self,parameters, z): return p.w0 + p.wa*(1-a) def _transverse_distance_(self,parameters, z): - # Use this print, in case needed for Error_handling. - # print '' + p = self.param(parameters) - # def Ly(y,t): - # return self.inv_hubble_rate(parameters,np.array([t])) + def Ly(y,t): + return self.inv_hubble_rate(parameters,np.array([t])) - # y=it.odeint(Ly,0.0,z) - # tint = np.array(y[:,0]) - # return self.clight* tint + y=it.odeint(Ly,0.0,z) + return cbackend.transverse_distance(p.H0,self.clight,p.Omega_k,y[:,0]) - p = self.param(parameters) - de = self.dark_energy_f(parameters, z) - return cbackend.transverse_distance(p.H0,p.Omega_m,p.Omega_b,p.Omega_k,de, z,self.clight) class LCDM(wCDM): def __init__(self, parameters,prior_file=None): @@ -326,22 +322,6 @@ def __check_mandatory_parameters__(self): assert len(self.param.parameters) == n, 'kwCDM: parameters are not correct' - def transverse_distance(self, parameters, z): - p = self.param(parameters) - def Ly(y,t): - return 1./self.hubble_rate(parameters,t) - - y=it.odeint(Ly,0.0,z) - tint = np.array(y[:,0]) - if p.Omega_k > 0.0: - return np.nan_to_num( self.clight/(np.sqrt(p.Omega_k)*p.H0)*np.sinh(np.sqrt(p.Omega_k) * p.H0*tint)) - - elif p.Omega_k < 0.0: - return np.nan_to_num( self.clight/(np.sqrt(abs(p.Omega_k))*p.H0)*np.sin((np.sqrt(abs(p.Omega_k)))* p.H0* tint)) - - elif p.Omega_k == 0.0: - return np.nan_to_num( self.clight* tint) - class kLCDM(kwCDM): def __init__(self, parameters,prior_file=None): @@ -394,21 +374,6 @@ def __check_mandatory_parameters__(self): n+=1 assert len(self.param.parameters) == n, 'kCPL3: parameters are not correct' - def transverse_distance(self, parameters, z): - p = self.param(parameters) - def Ly(y,t): - return 1./self.hubble_rate(parameters,t) - - y=it.odeint(Ly,0.0,z) - tint = np.array(y[:,0]) - if p.Omega_k > 0.0: - return np.nan_to_num( self.clight/(np.sqrt(p.Omega_k)*p.H0)*np.sinh(np.sqrt(p.Omega_k) * p.H0*tint)) - - elif p.Omega_k < 0.0: - return np.nan_to_num( self.clight/(np.sqrt(abs(p.Omega_k))*p.H0)*np.sin((np.sqrt(abs(p.Omega_k)))* p.H0* tint)) - - elif p.Omega_k == 0.0: - return np.nan_to_num( self.clight* tint) class kXCDM(XCDM): @@ -428,21 +393,5 @@ def __check_mandatory_parameters__(self): if self.Mbsample: n+=1 assert len(self.param.parameters) == n, 'kXCDM: parameters are not correct' - - def transverse_distance(self, parameters, z): - """transverse distance in Mpc/h""" - p = self.param(parameters) - def Ly(y,t): - return 1./self.hubble_rate(parameters,t) - y=it.odeint(Ly,0.0,z) - tint = np.array(y[:,0]) - if p.Omega_k > 0.0: - return np.nan_to_num( self.clight/(np.sqrt(p.Omega_k)*p.H0)*np.sinh(np.sqrt(p.Omega_k) * p.H0*tint)) - - elif p.Omega_k < 0.0: - return np.nan_to_num( self.clight/(np.sqrt(abs(p.Omega_k))*p.H0)*np.sin((np.sqrt(abs(p.Omega_k)))* p.H0* tint)) - - elif p.Omega_k == 0.0: - return np.nan_to_num( self.clight* tint) diff --git a/marcia/likelihood.py b/marcia/likelihood.py index 03a3d43..70021be 100644 --- a/marcia/likelihood.py +++ b/marcia/likelihood.py @@ -2,6 +2,11 @@ from marcia import Data from marcia import Cosmology as cosmo +import sys +from marcia import GPconfig +from marcia import CosmoConfig +from marcia.kernel import Kernel as kern + class Likelihood(object): # This contains the likelihoods for all the datasets implemented in the code as yet @@ -87,3 +92,109 @@ def logProb(self, theta): return -np.inf return lp + self.logLike(theta) + + +class Likelihood_GP(object): + # This contains the generalized likelihood for the GP and MTGP + def __init__(self, data, GP_filename=None, cosmo_filename=None): + + # Read form the GP config file + MTGP = GPconfig.GPConfig(GP_filename if GP_filename else 'GPconfig.ini').__dict__ + + # read the models and the number of tasks + self.models = MTGP['models'] + self.nTasks = MTGP['n_Tasks'] + self.nmodel = self.nTasks + self.self_scale = MTGP['self_scale'] + self.scatter = MTGP['sigma_int'] + + # To define the priors for the GP hyperparameters using the the GPconfig file + self.priors = [] + for i in range(self.nTasks): + i = i+1 # To start from 1 + self.priors.append([MTGP[f'Task_{i}']['sigma_f_min'], MTGP[f'Task_{i}']['sigma_f_max']]) + self.priors.append([MTGP[f'Task_{i}']['l_min'], MTGP[f'Task_{i}']['l_max']]) + + self.priors = np.array(self.priors) + + # validate the number of models and data + if len(self.models) != len(data): + raise ValueError('The number of models and the number of data are not equal') + + # read the data and set the data covarince matrix + self.d = Data(data) + self.x, self.y, self.cov = self.d() + self.D_covmat = self.cov + + # To set the mean function + self.mean = [0.0]*len(self.y) + # if MTGP['cosmo_model'] == None: + # self.mean = [0.0]*len(self.y) + # else: + # # To read the cosmology config file + # raise NotImplementedError('Cosmology mean is not implemented yet') + + # intialize the kernals from Kernel class + self.GP_kernels = kern(list(self.d.x.values()), filename=GP_filename) + + def set_theta(self, theta): + # To reorganise the theta values for the GP hyperparameters in accordanace with the GPconfig file + # theta contains no self-scale: [sigma_f_1, l_1, sigma_f_2, l_2, ...] + # or [sigma_f_1, sigma_f_2, ..., l] if self_scale is True + # and possibly with intrinsic scatter at the end [..., sigma_int] + # Final form is [[sigma_f_1, l_1], [sigma_f_2, l_2], ...] + + if self.self_scale == False and len(theta) != 2*self.nTasks and self.scatter == False: + raise ValueError('The number of hyperparameters does not match the number of tasks') + elif self.self_scale == False and len(theta) != 2*self.nTasks+1 and self.scatter == True: + pass + else: + # rearrange the theta values into a list of lists + params = [] + for i in range(self.nTasks): + params.append([theta[2*i], theta[2*i+1]]) + + return params + + def logPrior(self, theta): + # To set the priors for the GP hyperparameters + for (lower, upper), value in zip(self.priors, theta): + if not lower < value < upper: + # print(value, lower, upper) + # sys.exit() + return -np.inf + return 0.0 + + def logLike(self,theta): + chi2 = self.chisq(theta) + return -0.5*chi2 + + def logProb(self, theta): + lp = self.logPrior(theta) + theta = self.set_theta(theta) + if not np.isfinite(lp): + return -np.inf + return lp + self.logLike(theta) + + def chisq(self,theta): + # We define the GP likelihood here + chi2 = 0 + # To call the GP class to get the covariance matrix + GP_covmat = self.GP_kernels.Cov_Mat(theta) + Total_covmat = GP_covmat + self.D_covmat + 1e-6*np.eye(len(self.D_covmat)) + + # To perform cholesky decomposition + L_inv = np.linalg.inv(Total_covmat) + y = self.y - self.mean + + # To calculate the log determinant + log_det = np.linalg.slogdet(Total_covmat)[1] + offset = len(y)*np.log(2*np.pi) + + # To calculate the chi2 + y_inv = np.dot(L_inv, y) + chi2 = np.dot(y, y_inv.T) + log_det - offset + + return chi2 + + diff --git a/marcia/likelihood_GP.py b/marcia/likelihood_GP.py deleted file mode 100644 index 48d5cb4..0000000 --- a/marcia/likelihood_GP.py +++ /dev/null @@ -1,117 +0,0 @@ -import sys -import numpy as np -from marcia import Data -from marcia import GPconfig -from marcia import CosmoConfig -from marcia import Cosmology as cosmo -from marcia import Kernel as kern - - -class Likelihood_GP(object): - # This contains the generalized likelihood for the GP and MTGP - def __init__(self, data, GP_filename=None, cosmo_filename=None): - - # Read form the GP config file - MTGP = GPconfig.GPConfig(GP_filename if GP_filename else 'GPconfig.ini').__dict__ - - # read the models and the number of tasks - self.models = MTGP['models'] - self.nTasks = MTGP['n_Tasks'] - self.nmodel = self.nTasks - self.self_scale = MTGP['self_scale'] - self.scatter = MTGP['sigma_int'] - - # To define the priors for the GP hyperparameters using the the GPconfig file - self.priors = [] - for i in range(self.nTasks): - i = i+1 # To start from 1 - self.priors.append([MTGP[f'Task_{i}']['sigma_f_min'], MTGP[f'Task_{i}']['sigma_f_max']]) - self.priors.append([MTGP[f'Task_{i}']['l_min'], MTGP[f'Task_{i}']['l_max']]) - - self.priors = np.array(self.priors) - - # validate the number of models and data - if len(self.models) != len(data): - raise ValueError('The number of models and the number of data are not equal') - - # read the data and set the data covarince matrix - self.d = Data(data) - self.x, self.y, self.cov = self.d() - self.D_covmat = self.cov - - # To set the mean function - self.mean = [0.0]*len(self.y) - # if MTGP['cosmo_model'] == None: - # self.mean = [0.0]*len(self.y) - # else: - # # To read the cosmology config file - # raise NotImplementedError('Cosmology mean is not implemented yet') - - # intialize the kernals from Kernel class - self.GP_kernels = kern(list(self.d.x.values()), filename=GP_filename) - - def set_theta(self, theta): - # To reorganise the theta values for the GP hyperparameters in accordanace with the GPconfig file - # theta contains no self-scale: [sigma_f_1, l_1, sigma_f_2, l_2, ...] - # or [sigma_f_1, sigma_f_2, ..., l] if self_scale is True - # and possibly with intrinsic scatter at the end [..., sigma_int] - # Final form is [[sigma_f_1, l_1], [sigma_f_2, l_2], ...] - - if self.self_scale == False and len(theta) != 2*self.nTasks and self.scatter == False: - raise ValueError('The number of hyperparameters does not match the number of tasks') - elif self.self_scale == False and len(theta) != 2*self.nTasks+1 and self.scatter == True: - pass - else: - # rearrange the theta values into a list of lists - params = [] - for i in range(self.nTasks): - params.append([theta[2*i], theta[2*i+1]]) - - return params - - def logPrior(self, theta): - # To set the priors for the GP hyperparameters - for (lower, upper), value in zip(self.priors, theta): - if not lower < value < upper: - # print(value, lower, upper) - # sys.exit() - return -np.inf - return 0.0 - - def logLike(self,theta): - chi2 = self.chisq(theta) - return -0.5*chi2 - - def logProb(self, theta): - lp = self.logPrior(theta) - theta = self.set_theta(theta) - if not np.isfinite(lp): - return -np.inf - return lp + self.logLike(theta) - - def chisq(self,theta): - # We define the GP likelihood here - chi2 = 0 - # To call the GP class to get the covariance matrix - GP_covmat = self.GP_kernels.Cov_Mat(theta) - Total_covmat = GP_covmat + self.D_covmat + 1e-6*np.eye(len(self.D_covmat)) - - # To perform cholesky decomposition - L_inv = np.linalg.inv(Total_covmat) - y = self.y - self.mean - - # To calculate the log determinant - log_det = np.linalg.slogdet(Total_covmat)[1] - offset = len(y)*np.log(2*np.pi) - - # To calculate the chi2 - y_inv = np.dot(L_inv, y) - chi2 = np.dot(y, y_inv.T) + log_det - offset - - return chi2 - - - - - -