I would like to find an estimation of omega_BD parameter and constrain it in Brans-Dicke context, one suggests me to make a "grid" of likelihood. For this, I am using the data from Planck.
Below the script implementing this method with HI-CLASS since it can generate [math] of Brans-Dicke :
Code: Select all
import numpy as np
import matplotlib.pyplot as plt
from hi_classy import Class # Importing Hi-CLASS
import clik # For Planck's likelihood
# Initializing Planck's likelihood
path_to_planck_likelihood = './baseline/plc_3.0/hi_l/plik_lite/plik_lite_v22_TT.clik'
planck_likelihood = clik.clik(path_to_planck_likelihood)
lmax_tuple = planck_likelihood.get_lmax()
lmax = lmax_tuple[0]
n_cls = lmax - 1
def read_ini(filename):
with open(filename, 'r') as file:
lines = file.readlines()
parameters = {}
for line in lines:
line = line.strip()
if not line.startswith("#") and '=' in line:
key, value = line.split('=')
parameters[key.strip()] = value.strip()
return parameters
# Read the brans_dicke.ini file only once
params = read_ini('./brans_dicke.ini')
def get_spectrum(omega_BD):
# Use the previously initialized params dictionary and update only the value of omega_BD
omega_BD = float(omega_BD[0])
params['parameters_smg'] = '0.,' + str(omega_BD) + ', 1., 1e-3'
print('omega_BD = ', omega_BD)
print('params = ', params)
# Set up and run CLASS
cosmology = Class()
cosmology.set(params)
cosmology.compute()
# Obtain the relevant spectra
cl = cosmology.lensed_cl(2500)
spectrum = cl['tt']
spectrum = np.append(spectrum, [0]*9)
cosmology.struct_cleanup()
cosmology.empty()
return spectrum
# Continue with the Hi-CLASS initialization as in your previous code...
def compute_likelihood(omega_BD):
spectrum = get_spectrum([omega_BD])
# Add zeros first and then cut to the required length
spectrum_extended = np.append(spectrum, [0, 0, 0])
cls = spectrum_extended[:2510]
cls = cls.reshape(1, -1)
# Compute chi-squared using Planck's likelihood for the given theoretical Cl's
chi2 = planck_likelihood(cls)
return -0.5 * chi2
# Set up a range of omega_BD values for which you want to compute the likelihood
omega_BD_values = np.linspace(100, 1e5, 100) # Adjust the range and number of points as needed
# Compute the likelihood for each omega_BD value
likelihoods = [compute_likelihood(val) for val in omega_BD_values]
print(likelihoods[:10])
# Plot the likelihood as a function of omega_BD
plt.plot(omega_BD_values, likelihoods)
plt.xlabel('omega_BD')
plt.ylabel('Likelihood')
plt.title('Likelihood as a function of omega_BD')
plt.grid(True)
plt.show()
Code: Select all
Parameter file for Brans Dicke theory
-> See Avilez+13 (1303.4330) or Bellini+17 (1709.09135)
Omega_Lambda = 0
Omega_fld = 0
Omega_smg = -1
1) Brans-Dicke theory: G_2 -> (omega_BD/phi)X -V0,
G_4 -> phi M_p^2/2
initial field value sets M_*^2 (effective Planck mass)
gravity_model = brans dicke
V0, omega_BD, phi_ini, phi'_ini
parameters_smg = 0.7, 800, 1., 1e-3
2) Options for tuning of effective M_pl using phi_ini
Note that (b,c) are very close for large omega_BD (viable models)
a) No specifications uses the input_value
#comment M_pl_today_smg line below
b) Specifying M_pl_today fixes phi_0
(overwrites phi_ini)
M_pl_today_smg = 1.0
c) Specifying M_pl_today AND normalize_G_NR
sets G_eff -> 1 for non-relativistic matter
(overwrites phi_ini AND M_pl_today_smg)
#normalize_G_NR = yes
3) Brans Dicke runs into some numerical problems at very early times
a_min_stability_test_smg = 1e-6
------ other parameters --------
root = output/brans_dicke_
output = tCl, pCl, lCl, mPk
lensing = yes
#output = tCl,mPk
output_background_smg = 10
write parameters = yeap
write background = yup
write thermodynamics = yeahh
#input_verbose = 4
#background_verbose = 4
#output_verbose = 1
#thermodynamics_verbose = 1
#perturbations_verbose = 2
#spectra_verbose = 2
equal to nan :
For example, here is a typical output during the run :
Code: Select all
omega_BD = 88900.0
params = {'Omega_Lambda': '0', 'Omega_fld': '0', 'Omega_smg': '-1', 'gravity_model': 'brans dicke', 'parameters_smg': '0.,88900.0, 1., 1e-3', 'M_pl_today_smg': '1.0', 'a_min_stability_test_smg': '1e-6', 'root': 'output/brans_dicke_', 'output': 'tCl, pCl, lCl, mPk', 'lensing': 'yes', 'output_background_smg': '10', 'write parameters': 'yeap', 'write background': 'yup', 'write thermodynamics': 'yeahh'}
M_pl = 1.000092e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.000092e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.001092e+00, want 1.000000e+00, param=1.001000e+00
M_pl = 9.999997e-01, want 1.000000e+00, param=9.999093e-01
M_pl = 9.999997e-01, want 1.000000e+00, param=9.999093e-01
M_pl = 9.989996e-01, want 1.000000e+00, param=9.989093e-01
omega_BD = 89909.09090909091
params = {'Omega_Lambda': '0', 'Omega_fld': '0', 'Omega_smg': '-1', 'gravity_model': 'brans dicke', 'parameters_smg': '0.,89909.09090909091, 1., 1e-3', 'M_pl_today_smg': '1.0', 'a_min_stability_test_smg': '1e-6', 'root': 'output/brans_dicke_', 'output': 'tCl, pCl, lCl, mPk', 'lensing': 'yes', 'output_background_smg': '10', 'write parameters': 'yeap', 'write background': 'yup', 'write thermodynamics': 'yeahh'}
M_pl = 1.000091e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.000091e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.001091e+00, want 1.000000e+00, param=1.001000e+00
M_pl = 9.999997e-01, want 1.000000e+00, param=9.999103e-01
M_pl = 9.999996e-01, want 1.000000e+00, param=9.999103e-01
M_pl = 9.989996e-01, want 1.000000e+00, param=9.989103e-01
omega_BD = 90918.18181818182
params = {'Omega_Lambda': '0', 'Omega_fld': '0', 'Omega_smg': '-1', 'gravity_model': 'brans dicke', 'parameters_smg': '0.,90918.18181818182, 1., 1e-3', 'M_pl_today_smg': '1.0', 'a_min_stability_test_smg': '1e-6', 'root': 'output/brans_dicke_', 'output': 'tCl, pCl, lCl, mPk', 'lensing': 'yes', 'output_background_smg': '10', 'write parameters': 'yeap', 'write background': 'yup', 'write thermodynamics': 'yeahh'}
M_pl = 1.000090e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.000090e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.001090e+00, want 1.000000e+00, param=1.001000e+00
M_pl = 9.999997e-01, want 1.000000e+00, param=9.999113e-01
M_pl = 9.999997e-01, want 1.000000e+00, param=9.999113e-01
M_pl = 9.989996e-01, want 1.000000e+00, param=9.989113e-01
omega_BD = 91927.27272727274
params = {'Omega_Lambda': '0', 'Omega_fld': '0', 'Omega_smg': '-1', 'gravity_model': 'brans dicke', 'parameters_smg': '0.,91927.27272727274, 1., 1e-3', 'M_pl_today_smg': '1.0', 'a_min_stability_test_smg': '1e-6', 'root': 'output/brans_dicke_', 'output': 'tCl, pCl, lCl, mPk', 'lensing': 'yes', 'output_background_smg': '10', 'write parameters': 'yeap', 'write background': 'yup', 'write thermodynamics': 'yeahh'}
M_pl = 1.000089e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.000089e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.001089e+00, want 1.000000e+00, param=1.001000e+00
M_pl = 9.999997e-01, want 1.000000e+00, param=9.999123e-01
M_pl = 9.999997e-01, want 1.000000e+00, param=9.999123e-01
M_pl = 9.989996e-01, want 1.000000e+00, param=9.989123e-01
omega_BD = 92936.36363636363
params = {'Omega_Lambda': '0', 'Omega_fld': '0', 'Omega_smg': '-1', 'gravity_model': 'brans dicke', 'parameters_smg': '0.,92936.36363636363, 1., 1e-3', 'M_pl_today_smg': '1.0', 'a_min_stability_test_smg': '1e-6', 'root': 'output/brans_dicke_', 'output': 'tCl, pCl, lCl, mPk', 'lensing': 'yes', 'output_background_smg': '10', 'write parameters': 'yeap', 'write background': 'yup', 'write thermodynamics': 'yeahh'}
M_pl = 1.000088e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.000088e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.001088e+00, want 1.000000e+00, param=1.001000e+00
M_pl = 9.999996e-01, want 1.000000e+00, param=9.999132e-01
M_pl = 9.999996e-01, want 1.000000e+00, param=9.999132e-01
M_pl = 9.989996e-01, want 1.000000e+00, param=9.989132e-01
omega_BD = 93945.45454545454
params = {'Omega_Lambda': '0', 'Omega_fld': '0', 'Omega_smg': '-1', 'gravity_model': 'brans dicke', 'parameters_smg': '0.,93945.45454545454, 1., 1e-3', 'M_pl_today_smg': '1.0', 'a_min_stability_test_smg': '1e-6', 'root': 'output/brans_dicke_', 'output': 'tCl, pCl, lCl, mPk', 'lensing': 'yes', 'output_background_smg': '10', 'write parameters': 'yeap', 'write background': 'yup', 'write thermodynamics': 'yeahh'}
M_pl = 1.000087e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.000087e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.001087e+00, want 1.000000e+00, param=1.001000e+00
M_pl = 9.999997e-01, want 1.000000e+00, param=9.999142e-01
M_pl = 9.999997e-01, want 1.000000e+00, param=9.999142e-01
M_pl = 9.989996e-01, want 1.000000e+00, param=9.989142e-01
omega_BD = 94954.54545454546
params = {'Omega_Lambda': '0', 'Omega_fld': '0', 'Omega_smg': '-1', 'gravity_model': 'brans dicke', 'parameters_smg': '0.,94954.54545454546, 1., 1e-3', 'M_pl_today_smg': '1.0', 'a_min_stability_test_smg': '1e-6', 'root': 'output/brans_dicke_', 'output': 'tCl, pCl, lCl, mPk', 'lensing': 'yes', 'output_background_smg': '10', 'write parameters': 'yeap', 'write background': 'yup', 'write thermodynamics': 'yeahh'}
M_pl = 1.000086e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.000086e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.001086e+00, want 1.000000e+00, param=1.001000e+00
M_pl = 9.999997e-01, want 1.000000e+00, param=9.999151e-01
M_pl = 9.999997e-01, want 1.000000e+00, param=9.999151e-01
M_pl = 9.989996e-01, want 1.000000e+00, param=9.989151e-01
omega_BD = 95963.63636363637
params = {'Omega_Lambda': '0', 'Omega_fld': '0', 'Omega_smg': '-1', 'gravity_model': 'brans dicke', 'parameters_smg': '0.,95963.63636363637, 1., 1e-3', 'M_pl_today_smg': '1.0', 'a_min_stability_test_smg': '1e-6', 'root': 'output/brans_dicke_', 'output': 'tCl, pCl, lCl, mPk', 'lensing': 'yes', 'output_background_smg': '10', 'write parameters': 'yeap', 'write background': 'yup', 'write thermodynamics': 'yeahh'}
M_pl = 1.000085e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.000085e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.001085e+00, want 1.000000e+00, param=1.001000e+00
M_pl = 9.999997e-01, want 1.000000e+00, param=9.999160e-01
M_pl = 9.999997e-01, want 1.000000e+00, param=9.999160e-01
M_pl = 9.989996e-01, want 1.000000e+00, param=9.989160e-01
omega_BD = 96972.72727272728
params = {'Omega_Lambda': '0', 'Omega_fld': '0', 'Omega_smg': '-1', 'gravity_model': 'brans dicke', 'parameters_smg': '0.,96972.72727272728, 1., 1e-3', 'M_pl_today_smg': '1.0', 'a_min_stability_test_smg': '1e-6', 'root': 'output/brans_dicke_', 'output': 'tCl, pCl, lCl, mPk', 'lensing': 'yes', 'output_background_smg': '10', 'write parameters': 'yeap', 'write background': 'yup', 'write thermodynamics': 'yeahh'}
M_pl = 1.000084e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.000084e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.001084e+00, want 1.000000e+00, param=1.001000e+00
M_pl = 9.999997e-01, want 1.000000e+00, param=9.999169e-01
M_pl = 9.999997e-01, want 1.000000e+00, param=9.999169e-01
M_pl = 9.989997e-01, want 1.000000e+00, param=9.989169e-01
omega_BD = 97981.81818181819
params = {'Omega_Lambda': '0', 'Omega_fld': '0', 'Omega_smg': '-1', 'gravity_model': 'brans dicke', 'parameters_smg': '0.,97981.81818181819, 1., 1e-3', 'M_pl_today_smg': '1.0', 'a_min_stability_test_smg': '1e-6', 'root': 'output/brans_dicke_', 'output': 'tCl, pCl, lCl, mPk', 'lensing': 'yes', 'output_background_smg': '10', 'write parameters': 'yeap', 'write background': 'yup', 'write thermodynamics': 'yeahh'}
M_pl = 1.000083e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.000083e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.001084e+00, want 1.000000e+00, param=1.001000e+00
M_pl = 9.999997e-01, want 1.000000e+00, param=9.999177e-01
M_pl = 9.999997e-01, want 1.000000e+00, param=9.999177e-01
M_pl = 9.989996e-01, want 1.000000e+00, param=9.989177e-01
omega_BD = 98990.90909090909
params = {'Omega_Lambda': '0', 'Omega_fld': '0', 'Omega_smg': '-1', 'gravity_model': 'brans dicke', 'parameters_smg': '0.,98990.90909090909, 1., 1e-3', 'M_pl_today_smg': '1.0', 'a_min_stability_test_smg': '1e-6', 'root': 'output/brans_dicke_', 'output': 'tCl, pCl, lCl, mPk', 'lensing': 'yes', 'output_background_smg': '10', 'write parameters': 'yeap', 'write background': 'yup', 'write thermodynamics': 'yeahh'}
M_pl = 1.000083e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.000083e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.001083e+00, want 1.000000e+00, param=1.001000e+00
M_pl = 9.999998e-01, want 1.000000e+00, param=9.999186e-01
M_pl = 9.999998e-01, want 1.000000e+00, param=9.999186e-01
M_pl = 9.989997e-01, want 1.000000e+00, param=9.989186e-01
omega_BD = 100000.0
params = {'Omega_Lambda': '0', 'Omega_fld': '0', 'Omega_smg': '-1', 'gravity_model': 'brans dicke', 'parameters_smg': '0.,100000.0, 1., 1e-3', 'M_pl_today_smg': '1.0', 'a_min_stability_test_smg': '1e-6', 'root': 'output/brans_dicke_', 'output': 'tCl, pCl, lCl, mPk', 'lensing': 'yes', 'output_background_smg': '10', 'write parameters': 'yeap', 'write background': 'yup', 'write thermodynamics': 'yeahh'}
M_pl = 1.000082e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.000082e+00, want 1.000000e+00, param=1.000000e+00
M_pl = 1.001082e+00, want 1.000000e+00, param=1.001000e+00
M_pl = 9.999997e-01, want 1.000000e+00, param=9.999194e-01
M_pl = 9.999997e-01, want 1.000000e+00, param=9.999194e-01
M_pl = 9.989997e-01, want 1.000000e+00, param=9.989194e-01
[array([nan]), array([nan]), array([nan]), array([nan]), array([nan]), array([nan]), array([nan]), array([nan]), array([nan]), array([nan])]
I think that I need a right input file brans_dicke.ini
Any help is welcome,
Regards
