Brans-Dicke - Estimation of a single parameter omega_BD with a grid of Likelihood

Use of Cobaya. camb, CLASS, cosmomc, compilers, etc.
Post Reply
Fabien Dournac
Posts: 20
Joined: March 27 2010
Affiliation: CERFACS
Contact:

Brans-Dicke - Estimation of a single parameter omega_BD with a grid of Likelihood

Post by Fabien Dournac » November 02 2023

Hello,

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()
And the brans_dicke.ini input file :

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
Unfortunately, during the execution, I get systematically param values lower than M_pl : this implies at the end that likelihood is always
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])]
How to circumvent this issue of param values that are always lower than M_pl and want parameters ?

I think that I need a right input file brans_dicke.ini

Any help is welcome,

Regards

Fabien Dournac
Posts: 20
Joined: March 27 2010
Affiliation: CERFACS
Contact:

Re: Brans-Dicke - Estimation of a single parameter omega_BD with a grid of Likelihood

Post by Fabien Dournac » November 06 2023

Hi,

after multiple investigations, I found out how to produce the likelihood of omega_BD in Brans-Dicke context.

You can see the following plot :

Image

After higher value for omega_BD, I get a plateau with small fluctuations (like a Chi2). You can notice that -2*log(L) on y coordinates are pretty high ( ~ +1e8) : is it normal ?

More generally, does Planck likelihood from clik return log(likelihood) or simple likelihood ?

Regards

Post Reply