Example1:hi_class paper fig2

In this entry we will show how to obtain the figure 2 of hi_class paper (arXiv:1605.06102).

With hi_class (./class)

Compute the models

First of all, we need to create the corresponding ini-files. The figure's caption tells us that they used "lcdm" as expansion model and "propto_omega" as gravity model. In addition, they chose x_k = 1, x_t = x_m = 0, M*^2_ini = 1. This configuration can be summarized in the following lines:

Omega_Lambda = 0
Omega_fld = 0
Omega_smg = -1
expansion_model = lcdm
gravity_model = propto_omega
                #x_k, x_b, x_m, x_t, M*^2_ini
parameters_smg = 1., 0.625, 0., 0., 1.
output = tCl,pCl,lCl,mPk
root = output/hi_class-paper-fig2/b0.625_  # Choose the path you want
write parameters = yeap  # We suggest writting them to be sure what file goes with what parameters

To compute the others models, we only have to change the parameters_smg line to one of:

parameters_smg = 1., 0., 0., 0., 1.
parameters_smg = 1., 0.625, 0., 0., 1.
parameters_smg = 1., 1.25, 0., 0., 1.
parameters_smg = 1., 1.875, 0., 0., 1.
parameters_smg = 1., 2.5, 0., 0., 1.

Before computing the model, we must create the output folder. In this case:

mkdir output/hi_class-paper/

Then, from the hi_class_public directory, run ./class /path/to/my_ini_file.ini This will create the requested files, i.e. b0.625_cl.dat, b0.625_pk.dat, b0.625_parameters.ini and b0.625_unused_parameters, in the output/hi_class-paper-fig2.

After doing the same for the other parameters, we only have to compute the LCDM results to compare with. The minimal ini-file would be

output = tCl,pCl,lCl,mPk
root = output/hi_class-paper-fig2/lcdm_
write parameters = yeap

And, again:

./class /path/to/my_lcdm_ini_file.ini

Now we can plot the results and compare the modified gravity models with LCDM

Plot the results (with python)

Let's create the following python script (or jupyter notebook) in the output directory (i.e. output/hi_class-paper-fig2):

# In[1]:
import numpy as np
from matplotlib import pyplot as plt

# In[2]:

#This way of calculating the relative deviation works because the the x-axis arrays are equal.
def reldev(y1, y2):
    return 100*(y2 - y1)/y1

# In[3]:

# Load files:
lcdmCl = np.loadtxt('./lcdm_cl.dat', unpack=True)
lcdmPk = np.loadtxt('./lcdm_pk.dat', unpack=True)
cl = {}
pk = {}
for alpha_b in [0, 0.625, 1.25, 1.875, 2.5]:
    cl[alpha_b] = np.loadtxt('./b{}_cl.dat'.format(alpha_b), unpack=True)  # {alpha_b: [l, (l(l+1)/2pi)cl^TT, ...]}
    pk[alpha_b] = np.loadtxt('./b{}_pk.dat'.format(alpha_b), unpack=True)  # {alpha_b: [k, P(k)]}

# In[4]:

f, ax = plt.subplots(2, 2, figsize=(15, 10))

# First plot LCDM results
clList = lcdmCl
pkList = lcdmPk
label = r'$\Lambda$'

ax[0, 0].semilogx(clList[0], clList[1], '--', label=label)
ax[0, 0].set_ylabel(r'$[l(l+1)/2\pi]c^{TT}_l$')

ax[0, 1].loglog(pkList[0], pkList[1], '--', label=label)
ax[0, 1].set_ylabel(r'$P(k)$')

ax[1, 0].semilogx(clList[0], reldev(clList[1], clList[1]), '--', label=label)
ax[1, 0].set_xlabel('l')
ax[1, 0].set_ylabel(r'rev.dev [%]')

ax[1, 1].semilogx(pkList[0], reldev(pkList[1], pkList[1]), '--', label=label)
ax[1, 1].set_xlabel('k')
ax[1, 1].set_ylabel(r'rev.dev [%]')

# Then plot the modified gravity results
for alpha_b in cl:
    clList = cl[alpha_b]
    pkList = pk[alpha_b]
    label = r'$\alpha_B={}$'.format(alpha_b)
    
    ax[0, 0].semilogx(clList[0], clList[1], label=label)    
    ax[0, 1].loglog(pkList[0], pkList[1] , label=label)
    ax[1, 0].semilogx(clList[0], reldev(lcdmCl[1], clList[1]), label=label)
    ax[1, 1].semilogx(pkList[0], reldev(lcdmPk[1], pkList[1]), label=label)

plt.legend(loc=0)
plt.show()
plt.close()

Now, from the output directory run:

python my_script.py

With classy

Recall that when we work with classy we do not use ini-files nor external files to write to. Nevertheless, we have to fill class with our parameter selection. This is done through cosmo.set(params). But first things first.

Since we are working with python, what is the langague of classy, we must start our program importing the needed modules:

# In[1]:

import numpy as np
from matplotlib import pyplot as plt
from classy import Class

Then, we would define the function which is going to calculate the relative deviation between outputs:

# In[2]:

#This way of calculating the relative deviation works because the the x-axis arrays are equal.
def reldev(y1, y2):
    return 100*(y2 - y1)/y1

With all the preliminaries done, we can start working with classy. First, we must initialize the Class class:

# In[3]:

cosmo = Class()

As before, we compute the LCDM apart from the other quantities. As we said, we need first state what parameters are needed. In this case, we only have to tell classy we want it to compute perturbations:

# # LCDM part

# In[4]:

paramsLCDM = {
    'output': 'tCl,pCl,lCl,mPk'
}


# In[5]:

cosmo.set(paramsLCDM)


# In[6]:

cosmo.compute()

We can access now the computed results stored in memory. We will use cosmo.raw_cl to obatin the raw cl's and cosmo.pk to obtain the power spectrum:

# In[7]:

lcdmCl = cosmo.raw_cl()


# In[8]:

kList = np.logspace(10**(-5), 10**(-1), 3000)
lcdmPk = np.array([cosmo.pk(k, 0) for k in klist])

Finally, we must clean the structures in order to save memory.

# In[9]:

cosmo.struct_cleanup()

We can start now with the alpha's part. In this case, there will be more parameters to input, as there were before. In this case, as we want to compute for different paramters_smg we loop through its values and do all computation (similar to the LCDM part) in it.

# # Alphas' part

# In[10]:

# We copy the parameters of one of the parameters.ini
params = {
    'Omega_Lambda': 0,
    'Omega_fld': 0,
    'Omega_smg': -1,
    'expansion_model': 'lcdm',
    'gravity_model': 'propto_omega',
                    #x_k, x_b, x_m, x_t, M*^2_ini
    #parameters_smg = 1., 0.625, 0., 0., 1.,
    'output': 'tCl,pCl,lCl,mPk'
}


# In[11]:

x_k = 1
x_b = 0
x_m = 0
x_t = 0
M = 1


# In[12]:

cl = {}
pk = {}

for x_b in [0, 0.625, 1.25, 1.875, 2.5]:
    params['parameters_smg'] = "{}, {}, {}, {}, {}".format(x_k, x_b, x_m, x_t, M)
    cosmo.set(params)
    cosmo.compute()
    cl[x_b] = cosmo.raw_cl()
    pk[x_b] = np.array([cosmo.pk(k,0) for k in kList])
    cosmo.struct_cleanup()

With everything computed, we can plot the results:

# In[13]:

f, ax = plt.subplots(2, 2, figsize=(15, 10))

# First plot LCDM results
clList = [np.arange(len(lcdmCl['tt'][2:]))+2]
clList.append(lcdmCl['tt'][2:]*(clList[0]*(clList[0]+1)/(2*np.pi)))
pkList = lcdmPk
label = r'$\Lambda$'

ax[0, 0].semilogx(clList[0], clList[1], '--', label=label)
ax[0, 0].set_ylabel(r'$[l(l+1)/2\pi]c^{TT}_l$')

ax[0, 1].loglog(pkList[0], pkList[1], '--', label=label)
ax[0, 1].set_ylabel(r'$P(k)$')

ax[1, 0].semilogx(clList[0], reldev(clList[1], clList[1]), '--', label=label)
ax[1, 0].set_xlabel('l')
ax[1, 0].set_ylabel(r'rev.dev [%]')

ax[1, 1].semilogx(pkList[0], reldev(pkList[1], pkList[1]), '--', label=label)
ax[1, 1].set_xlabel('k')
ax[1, 1].set_ylabel(r'rev.dev [%]')

# Then plot the modified gravity results
lcdmClList = list(clList)  # Make a copy of the array, otherwise it would be referenced to the same memory address.
for alpha_b in cl:
    clList[1] = cl[alpha_b]['tt'][2:]*(clList[0]*(clList[0]+1)/(2*np.pi))
    pkList = pk[alpha_b]
    label = r'$\alpha_B={}$'.format(alpha_b)
    
    ax[0, 0].semilogx(clList[0], clList[1], label=label)    
    ax[0, 1].loglog(kList, pkList , label=label)
    ax[1, 0].semilogx(clList[0], reldev(lcdmClList[1], clList[1]), label=label)
    ax[1, 1].semilogx(kList, reldev(lcdmPk, pkList), label=label)

plt.legend(loc=0)
plt.show()
plt.close()