quasar likelihood best

Author: antolonappan

Notebooks/likelihood.ipynb

@@ -133,6 +133,11 @@ def distance_modulus(self,parameters,z1, z2):
             mu = cbackend.distance_modulus(Mb,z2,d)
         return mu
 
+    def luminosity_distance(self,parameters,z):
+        p = self.param(parameters)
+        d = self.transverse_distance(parameters, z)
+        return d * (1+z)
+    
     # This is to provide the Mb alone
     def Abs_M(self,parameters):
         p = self.param(parameters)
@@ -184,12 +189,12 @@ def __check_mandatory_parameters__(self):
         assert 'H0' in self.param.parameters, 'wCDM: H0 is not defined in the parameters'
         assert 'w0' in self.param.parameters, 'wCDM: w0 is not defined in the parameters'
         assert 'Omega_m' in self.param.parameters, 'wCDM: Omega_m is not defined in the parameters'
-        n=3
-        if self.rdsample or self.Obsample:
-            n+=1
-        if self.Mbsample:
-            n+=1
-        assert len(self.param.parameters) == n, 'wCDM: parameters are not correct'
+        # n=3
+        # if self.rdsample or self.Obsample:
+        #     n+=1
+        # if self.Mbsample:
+        #     n+=1
+        # assert len(self.param.parameters) == n, 'wCDM: parameters are not correct'
 
     def dark_energy_f(self, parameters, z):
         p = self.param(parameters)
@@ -218,12 +223,12 @@ def __init__(self, parameters,prior_file=None):
     def __check_mandatory_parameters__(self):
         assert 'H0' in self.param.parameters, 'LCDM: H0 is not defined in the parameters'
         assert 'Omega_m' in self.param.parameters, 'LCDM: Omega_m is not defined in the parameters'
-        n=2
-        if self.rdsample or self.Obsample:
-            n+=1
-        if self.Mbsample:
-            n+=1
-        assert len(self.param.parameters) == n, 'LCDM: parameters are not correct'
+        # n=2
+        # if self.rdsample or self.Obsample:
+        #     n+=1
+        # if self.Mbsample:
+        #     n+=1
+        # assert len(self.param.parameters) == n, 'LCDM: parameters are not correct'
 
 class CPL(wCDM):
     def __init__(self, parameters,prior_file=None):

marcia/cosmology.py

@@ -6,14 +6,14 @@
 
 
 __path__ = os.path.dirname(os.path.realpath(__file__))
-__datapath__ = os.path.join(__path__, 'Data')
+__datapath__ = os.path.join(__path__, '..','Data')
 
 
 
 class Data:
 
     def __init__(self,data,file_fs8=0,Lambda=1,b=1,sigma_sys=0.7571,H0=1,):
-        datalist = ['CC','BAO-alam','BAO-zhao','GR','Lya','GRB','SNE','QSA', 'Planck_TT', 'Planck_EE', 'Planck_TE','Pantheon_plus']
+        datalist = ['CC','BAO-alam','BAO-zhao','GR','Lya','GRB','SNE','QSO','QSO_full', 'Planck_TT', 'Planck_EE', 'Planck_TE','Pantheon_plus']
         if type(data) is str:
             assert data in datalist, f'{data} is not in {datalist}'
             self.data = [data]
@@ -186,12 +186,23 @@ def get_SNE(self):
         sigma = data[:,2]
         covar = sigma * corr * sigma.T
         return x,y,covar
+    
+    
 
-    def get_QSA(self):
-        data = loadtxt(os.path.join(__datapath__, 'Quasars','DL_all_short.txt'))
-        x = data[:,0]
-        y = data[:,1]
-        sigma = data[:,2]
+    
+    @load_data_once
+    def get_QSO_data(self):
+        data = loadtxt(os.path.join(__datapath__, 'Quasars','table3.dat'), usecols=(3,4,5,6,7,9,10, 11, 12))
+        z, lnFUV, lnFUV_err, lnFX, lnFX_err = data[:,0], data[:,1], data[:,2], data[:,3], data[:,4]
+        DM, dDM = data[:,7], data[:,8]
+        return z, lnFUV, lnFUV_err, lnFX, lnFX_err, DM, dDM
+    
+    def get_QSO_full(self):
+        z, lnFUV, lnFUV_err, lnFX, lnFX_err, _, _ = self.get_QSO_data()
+        return z, lnFUV, lnFUV_err, lnFX, lnFX_err
+
+    def get_QSO(self):
+        x,_,_,_,_,y,sigma = self.get_QSO_data()
         covar = np.diag(sigma**2)
         return x,y,covar
 

marcia/database.py

@@ -6,6 +6,9 @@
 from marcia import GPconfig
 from marcia import CosmoConfig
 from marcia.kernel import Kernel as kern
+from marcia import Params
+
+import pdb
 
 
 class Likelihood(object):
@@ -16,6 +19,7 @@ def __init__(self, model,parameters,data,prior_file=None):
         self.theory = cosmo(model,parameters,prior_file) # This needs to be done here to make the code more autonomous
         # This list sets the available data to work with, has to be updated when ever a new dataset is added
         self.priors = self.theory.priors
+        self.params = Params(parameters)
         self.data = data
         self.db = Data(data)
         self.inv_covariance = {}
@@ -64,6 +68,31 @@ def chisq_Pantheon_plus(self,theta):
             self.inv_covariance['Pantheon_plus'] = icov
 
         return np.dot(delta, np.dot(icov, delta))
+    
+    def chisq_QSO_dm(self,theta):
+        p = self.params(theta)
+        #pdb.set_trace()
+        z,dm,cov = self.db.get_QSO()
+        distance_theory = self.theory.distance_modulus(theta, z,z)
+        delta = dm - distance_theory
+        var = np.diag(cov)
+        return np.sum(delta**2./(var + p.qso_sigma**2))
+    
+    def chisq_QSO_full(self,theta):
+        p = self.params(theta)
+        z, lnFUV, lnFUV_err, lnFX, lnFX_err = self.db.get_QSO_full()
+        yi = lnFX
+        dyi = lnFX_err
+        xi = lnFUV
+        dxi = lnFUV_err
+        DL = self.theory.luminosity_distance(theta,z)*3.086e24
+
+        psi = p.qso_beta + p.qso_gamma*(xi) + 2*(p.qso_gamma -1)*(np.log10(DL)) + (p.qso_gamma-1)*np.log10(4*np.pi)
+
+        si2 = p.qso_gamma**2 * dxi**2 + dyi**2 + np.exp( 2*np.log(p.qso_sigma) )
+
+        return np.sum((psi - yi)**2/si2 + np.log(2*np.pi*si2))
+
 
 
     def chisq(self,theta):

marcia/likelihood.py

@@ -13,11 +13,14 @@ wc = '$w_c$'
 fa = '$f_a$'
 fb = '$f_b$'
 fc = '$f_c$'
+qso_beta = '$\beta_{QSO}$'
+qso_sigma = '$\sigma_{QSO}$'
+qso_gamma = '$\gamma_{QSO}$'
 
 
 [Priors]
 r_d = [100.0,180.0]
-M_b = [-20.0,-18.0]
+M_b = [-19.5,-18.7]
 Omega_b = [0.01,0.1]
 H0 = [61.0,81.0]
 Omega_m = [0.1,0.5]
@@ -30,6 +33,9 @@ wc = [-150.0,80.0]
 fa = [-20.0,20.0]
 fb = [-50.0,50.0]
 fc = [-50.0,50.0]
+qso_beta = [0,15]
+qso_sigma = [0,1]
+qso_gamma = [0,2]
 
 [Defaults]
 r_d = 0
@@ -46,4 +52,7 @@ wc = 0
 fa = 1
 fb = 0
 fc = 0
+qso_beta = 8
+qso_sigma = 0.2
+qso_gamma = 0.6