Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

marciab000 #30

Merged
merged 9 commits into from
Oct 7, 2023
Merged

marciab000 #30

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
e6a2d70
sampler distribution and burnin
Sandeepharidasu Sep 28, 2023
6b28658
GP likelihood with individual scatter
Sandeepharidasu Sep 28, 2023
a3e8c31
new GP config
Sandeepharidasu Sep 28, 2023
620e596
quasar bug
Sandeepharidasu Sep 28, 2023
07f57c0
GP likelihood with individual intrinsic scatter
Sandeepharidasu Sep 28, 2023
7f0ca26
QSO likelihood fixed
Sandeepharidasu Sep 29, 2023
fa47dd0
fix QSO data and sampler change
Sandeepharidasu Sep 29, 2023
d225078
pantheon pos
Sandeepharidasu Oct 4, 2023
0aa5538
if any
Sandeepharidasu Oct 5, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
36 changes: 18 additions & 18 deletions Data/Pantheon_old/lcparam_full.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1028,21 +1028,21 @@
1.800000000000000044e+00 1.799601000000000006e+00 2.620110000000000028e+01 1.893198849999999993e+02 6.227817900000000151e+01
2.259999999999999787e+00 2.259600999999999971e+00 2.687699999999999889e+01 1.891565860000000043e+02 6.230917399999999873e+01
1.913999999999999924e+00 1.914828999999999892e+00 2.616090000000000160e+01 3.444306600000000174e+01 -5.256643000000000399e+00
1.300000000000000044e+00 1.301201999999999970e+00 2.552015000000000100e+01 0.000000000000000000e+00 0.000000000000000000e+00
1.340000000000000080e+00 1.341202000000000005e+00 2.546519999999999939e+01 0.000000000000000000e+00 0.000000000000000000e+00
1.020000000000000018e+00 1.021201999999999943e+00 2.488234999999999886e+01 0.000000000000000000e+00 0.000000000000000000e+00
7.349999999999999867e-01 7.362020000000000231e-01 2.370605000000000118e+01 0.000000000000000000e+00 0.000000000000000000e+00
1.120000000000000107e+00 1.121202000000000032e+00 2.506390000000000029e+01 0.000000000000000000e+00 0.000000000000000000e+00
1.229999999999999982e+00 1.231201999999999908e+00 2.528614999999999924e+01 0.000000000000000000e+00 0.000000000000000000e+00
1.229999999999999982e+00 1.231201999999999908e+00 2.549849999999999994e+01 0.000000000000000000e+00 0.000000000000000000e+00
8.539999999999999813e-01 8.552020000000000177e-01 2.431279999999999930e+01 0.000000000000000000e+00 0.000000000000000000e+00
1.370000000000000107e+00 1.371202000000000032e+00 2.551589999999999847e+01 0.000000000000000000e+00 0.000000000000000000e+00
9.751999999999999558e-01 9.764019999999999921e-01 2.470609999999999928e+01 0.000000000000000000e+00 0.000000000000000000e+00
9.699999999999999734e-01 9.712020000000000097e-01 2.511544999999999916e+01 0.000000000000000000e+00 0.000000000000000000e+00
7.399999999999999911e-01 7.412020000000000275e-01 2.384459999999999980e+01 0.000000000000000000e+00 0.000000000000000000e+00
1.389999999999999902e+00 1.391202000000000050e+00 2.536305000000000121e+01 0.000000000000000000e+00 0.000000000000000000e+00
1.304999999999999938e+00 1.306202000000000085e+00 2.526419999999999888e+01 0.000000000000000000e+00 0.000000000000000000e+00
9.350000000000000533e-01 9.362019999999999786e-01 2.426660000000000039e+01 0.000000000000000000e+00 0.000000000000000000e+00
1.014000000000000012e+00 1.015201999999999938e+00 2.485444999999999993e+01 0.000000000000000000e+00 0.000000000000000000e+00
1.314999999999999947e+00 1.316202000000000094e+00 2.554234999999999900e+01 0.000000000000000000e+00 0.000000000000000000e+00
1.092000000000000082e+00 1.093202000000000007e+00 2.463834999999999908e+01 0.000000000000000000e+00 0.000000000000000000e+00
1.300000000000000044e+00 1.301201999999999970e+00 2.552015000000000100e+01 5.315630000000000166e+01 -2.777959999999999852e+01
1.340000000000000080e+00 1.341202000000000005e+00 2.546519999999999939e+01 1.892880000000000109e+02 6.219140000000000157e+01
1.020000000000000018e+00 1.021201999999999943e+00 2.488234999999999886e+01 1.893369999999999891e+02 6.222820000000000107e+01
7.349999999999999867e-01 7.362020000000000231e-01 2.370605000000000118e+01 5.309308300000000003e+01 -2.774080000000000013e+01
1.120000000000000107e+00 1.121202000000000032e+00 2.506390000000000029e+01 1.890579999999999927e+02 6.220210000000000150e+01
1.229999999999999982e+00 1.231201999999999908e+00 2.528614999999999924e+01 1.890954999999999870e+02 6.230644399999999905e+01
1.229999999999999982e+00 1.231201999999999908e+00 2.549849999999999994e+01 1.892359999999999900e+02 6.221479999999999677e+01
8.539999999999999813e-01 8.552020000000000177e-01 2.431279999999999930e+01 1.891450000000000102e+02 6.226359999999999673e+01
1.370000000000000107e+00 1.371202000000000032e+00 2.551589999999999847e+01 5.304180000000000206e+01 -2.783050000000000068e+01
9.751999999999999558e-01 9.764019999999999921e-01 2.470609999999999928e+01 5.310560000000000258e+01 -2.775080000000000169e+01
9.699999999999999734e-01 9.712020000000000097e-01 2.511544999999999916e+01 1.895374999999999943e+02 6.231312199999999990e+01
7.399999999999999911e-01 7.412020000000000275e-01 2.384459999999999980e+01 5.307562500000000227e+01 -2.773626399999999848e+01
1.389999999999999902e+00 1.391202000000000050e+00 2.536305000000000121e+01 1.892249999999999943e+02 6.213969999999999771e+01
1.304999999999999938e+00 1.306202000000000085e+00 2.526419999999999888e+01 5.310329999999999728e+01 -2.777159999999999940e+01
9.350000000000000533e-01 9.362019999999999786e-01 2.426660000000000039e+01 1.893710000000000093e+02 6.219109999999999872e+01
1.014000000000000012e+00 1.015201999999999938e+00 2.485444999999999993e+01 3.542530000000000001e+01 -3.364799999999999791e+00
1.314999999999999947e+00 1.316202000000000094e+00 2.554234999999999900e+01 3.544369999999999976e+01 -3.382299999999999862e+00
1.092000000000000082e+00 1.093202000000000007e+00 2.463834999999999908e+01 1.873568960000000061e+02 1.849053000000000058e+00
58 changes: 57 additions & 1 deletion Notebooks/Validation_cosmo.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,62 @@
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 38,
"id": "6767992e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Frodo\n",
"Koekemoer\n",
"Patuxent\n",
"Rakke\n",
"SCP06C0\n"
]
}
],
"source": [
"# To make the positions fo teh pantheon dataset \n",
"\n",
"\"\"\" More importantly I notice some of teh supernovae no longer exist in the Panhteon+ compilation I dont know why !!\n",
"\"\"\"\n",
"\n",
"pant_old_names = np.loadtxt('../Data/Pantheon_old/lcparam_full_long.txt', dtype='str', usecols=0, skiprows=1)\n",
"pant_old = np.loadtxt('../Data/Pantheon_old/lcparam_full.txt', )\n",
"pant_new_names = np.loadtxt('../Data/Pantheon+/Pantheon+SH0ES.dat', skiprows=1, dtype='str', usecols=0 )\n",
"pant_new = np.loadtxt('../Data/Pantheon+/Pantheon+SH0ES.dat', skiprows=1, usecols=(26, 27))\n",
"# taken from github of Pantheon+\n",
"others = np.array([[53.093083,-27.74080], [189.095500, 62.306444], [189.537500, 62.313122], [53.075625, -27.736264], [187.356896, 1.849053]]) \n",
"\n",
"j = 0\n",
"for i in range(len(pant_old_names[1030:])):\n",
" if pant_old_names[1030+i] in pant_new_names:\n",
" \n",
" index_new = np.where(pant_new_names == pant_old_names[1030+i])\n",
" pant_old[1030+i, 3] = pant_new[index_new, 0]\n",
" pant_old[1030+i, 4] = pant_new[index_new, 1]\n",
" else:\n",
" # print('Not found')\n",
" print(pant_old_names[1030+i])\n",
" pant_old[1030+i, 3] = others[j, 0]\n",
" pant_old[1030+i, 4] = others[j, 1]\n",
" j += 1\n",
"\n",
"# np.savetxt('../Data/Pantheon_old/lcparam_full_2.txt', pant_old)\n"
]
},
{
"cell_type": "code",
"execution_count": 36,
"id": "5cc1d83f",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 3,
Expand Down Expand Up @@ -482,7 +538,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.8"
"version": "3.8.11"
}
},
"nbformat": 4,
Expand Down
12 changes: 10 additions & 2 deletions marcia/GPconfig.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@
filename = 'GPconfig.ini'

# self scale implies all the tasks have the same length scale paramter.
config['GENERAL'] = {'n_Tasks': n_Tasks, 'self_scale': True}
config['GENERAL'] = {'n_Tasks': n_Tasks, 'self_scale': True, 'ind_scatter': False}

config['KERNEL'] = {'Interpolate': True, 'Method': 'cubic', 'n_points': 100 }

Expand All @@ -39,6 +39,8 @@
config.set(Task, 'l_max', '10.0')
config.set(Task, 'sigma_f_min', '0.001')
config.set(Task, 'sigma_f_max', '10.0')
config.set(Task, 'sig_int_min', '0.001')
config.set(Task, 'sig_int_max', '5.0')

# if necessary to include a intrinsic scatter or an offset to the covariance matrix
config['INTRINSIC_SCATTER'] = {'sigma_int': True, 'offset': True}
Expand Down Expand Up @@ -72,6 +74,9 @@ def __init__(self, filename):
# Set the self_scale parameter
self.self_scale = config.getboolean('GENERAL', 'self_scale')

# Set the ind_scatter parameter
self.ind_scatter = config.getboolean('GENERAL', 'ind_scatter')

# Set the interpolation parameters
self.Interpolate = config.getboolean('KERNEL', 'Interpolate')
self.Method = config.get('KERNEL', 'Method')
Expand Down Expand Up @@ -103,7 +108,10 @@ def __init__(self, filename):
l_max = config.getfloat(Task, 'l_max')
sigma_f_min = config.getfloat(Task, 'sigma_f_min')
sigma_f_max = config.getfloat(Task, 'sigma_f_max')
sig_int_min = config.getfloat(Task, 'sig_int_min')
sig_int_max = config.getfloat(Task, 'sig_int_max')

# Create the GPparams class
self.__dict__[Task] = {'model': model, 'nu': nu, 'l_min': l_min, 'l_max': l_max, 'sigma_f_min': sigma_f_min, 'sigma_f_max': sigma_f_max}
self.__dict__[Task] = {'model': model, 'nu': nu, 'l_min': l_min, 'l_max': l_max, 'sigma_f_min': sigma_f_min, 'sigma_f_max': sigma_f_max, 'sig_int_min': sig_int_min, 'sig_int_max': sig_int_max}


10 changes: 6 additions & 4 deletions marcia/database.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,8 +90,10 @@ def get_data(self,choose):
return self.get_GRB()
elif choose == 'SNE':
return self.get_SNE()
elif choose == 'QSA':
return self.get_QSA()
elif choose == 'QSO':
return self.get_QSO()
elif choose == 'QSO_full':
return self.get_QSO_full()
elif choose == 'Planck_TT':
return self.get_CMB_planck_TT()
elif choose == 'Planck_EE':
Expand Down Expand Up @@ -245,8 +247,8 @@ def get_CMB_planck_TE(self):

@load_data_once
def get_data_pantheon_old(self):
datapath = os.path.join(__datapath__, 'Pantheon_old','lcparam_select.txt')
covpath = os.path.join(__datapath__, 'Pantheon_old','cov_select.txt')
datapath = os.path.join(__datapath__, 'Pantheon_old','lcparam_full.txt')
covpath = os.path.join(__datapath__, 'Pantheon_old','cov_full.txt')
zcmb, zhel, mb, ra, dec = np.loadtxt(datapath).T
cov = np.loadtxt(covpath)
return zcmb, zhel, mb, cov, ra, dec
Expand Down
95 changes: 66 additions & 29 deletions marcia/likelihood.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ def chisq_BAO_alam(self,theta):

return np.dot(dist- distance_theory, np.dot(np.linalg.inv(cov), dist - distance_theory))

def chisq_Pantheon_plus(self,theta):
def chisq_Pantheon_plus(self,theta, residual = False):
cmb_z, mb, covariance = self.db.get_pantheon_plus()
helio_z = self.db.get_pantheon_plus(Zhel=True)
distance_theory = self.theory.distance_modulus(theta, cmb_z, helio_z)
Expand All @@ -66,9 +66,12 @@ def chisq_Pantheon_plus(self,theta):
icov = np.linalg.inv(covariance)
self.inv_covariance['Pantheon_plus'] = icov

return np.dot(delta, np.dot(icov, delta))
if residual:
return delta, np.sqrt(np.diag(covariance))
else:
return np.dot(delta, np.dot(icov, delta))

def chisq_Pantheon_old(self,theta):
def chisq_Pantheon_old(self,theta, residual = False):
cmb_z, mb, covariance = self.db.get_pantheon_old()
helio_z = self.db.get_pantheon_old(Zhel=True)
distance_theory = self.theory.distance_modulus(theta, cmb_z, helio_z)
Expand All @@ -78,17 +81,22 @@ def chisq_Pantheon_old(self,theta):
else:
icov = np.linalg.inv(covariance)
self.inv_covariance['Pantheon_old'] = icov

return np.dot(delta, np.dot(icov, delta))
if residual:
return delta, np.sqrt(np.diag(covariance))
else:
return np.dot(delta, np.dot(icov, delta))

def chisq_QSO_dm(self,theta):
def chisq_QSO(self,theta, residual = False):
p = self.params(theta)
#pdb.set_trace()
z,dm,cov = self.db.get_QSO()
distance_theory = self.theory.distance_modulus(theta, z,z)
distance_theory = self.theory.distance_modulus(theta, z,z) - p.M_b
delta = dm - distance_theory
var = np.diag(cov)
return np.sum(delta**2./(var + p.qso_sigma**2))
if residual:
return delta, np.sqrt(var)
else:
return np.sum(delta**2./(var + p.qso_sigma**2))

def chisq_QSO_full(self,theta):
p = self.params(theta)
Expand Down Expand Up @@ -147,13 +155,20 @@ def __init__(self, data, GP_filename=None, cosmo_filename=None):
self.nmodel = self.nTasks
self.self_scale = MTGP['self_scale']
self.scatter = MTGP['sigma_int']
self.offset = MTGP['offset']

# extra paramter per dataset to be added to the theta
self.ind_scatter = MTGP['ind_scatter']

# 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']])
if self.ind_scatter:
self.priors.append([MTGP[f'Task_{i}']['sig_int_min'], MTGP[f'Task_{i}']['sig_int_max']])


self.priors = np.array(self.priors)

Expand Down Expand Up @@ -184,44 +199,68 @@ def set_theta(self, theta):
# 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

if self.self_scale == False:
params_GP = []
if len(theta)>=2*self.nTasks:
for i in range(self.nTasks):
params_GP.append([theta[2*i], theta[2*i+1]])
else:
raise ValueError(f'The number of hyperparameters does not match the number of tasks, there must be atleast {2*self.nTasks} hyperparameters')

elif self.self_scale == True:
params_GP =[]
if len(theta)>=self.nTasks+1:
for i in range(self.nTasks):
params_GP.append([theta[i], theta[self.nTasks]])
else:
raise ValueError(f'The number of hyperparameters does not match the number of tasks, there must be atleast {self.nTasks+1} hyperparameters')


if self.ind_scatter == True and len(theta)>=self.nTasks+2:
params_scatter = theta[-self.nTasks:]
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]])
params_scatter = [self.scatter]*self.nTasks

return params
return params_GP, params_scatter

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)
def logLike(self,theta_GP, theta_scatter):
chi2 = self.chisq(theta_GP, theta_scatter)
return -0.5*chi2

def logProb(self, theta):

lp = self.logPrior(theta)
theta = self.set_theta(theta)
theta_GP, theta_scatter = self.set_theta(theta)
if not np.isfinite(lp):
return -np.inf
return lp + self.logLike(theta)
return lp + self.logLike(theta_GP, theta_scatter)

def chisq(self,theta):
def chisq(self, theta_GP, theta_scatter):
# We define the GP likelihood here
chi2 = 0

# To set the scatter covariance matrix which is a diagonal matrix of size len of ntasks
scatter_matrix = self.offset*np.eye(len(self.D_covmat))
scatter_diag = []
for i in range(self.nTasks):
# create arrays of length of the data for each task
scatter_diag = np.append(scatter_diag, len(np.array(self.d.x[f'x{i}']))*[theta_scatter[i]**2.])

int_scatter_covmat = scatter_matrix + np.diag(scatter_diag)

# 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))
GP_covmat = self.GP_kernels.Cov_Mat(theta_GP)

# To calculate the total covariance matrix
Total_covmat = GP_covmat + self.D_covmat + int_scatter_covmat

# To perform cholesky decomposition
L_inv = np.linalg.inv(Total_covmat)
Expand All @@ -235,6 +274,4 @@ def chisq(self,theta):
y_inv = np.dot(L_inv, y)
chi2 = np.dot(y, y_inv.T) + log_det - offset

return chi2


return chi2
Loading