CosmoLib Basic Tutorial

• Copy Right and Citations
• Compiling the package
• Producing the mock data
• Fisher matrix analysis
• MCMC
• Transfer functions

An advanced tutorial with more features and examples can be find at http://www.cita.utoronto.ca/~zqhuang/CosmoLib/advdoc.php.

The download page is http://www.cita.utoronto.ca/~zqhuang/CosmoLib/download.php.

CosmoLib is a free software under the GNU GENERAL PUBLIC LICENSE (GPL, http://www.gnu.org/copyleft/gpl.html); Feel free to share or modify it.

If you have used the software in your publication, please cite

Setting compiler and library options

I assume you are in a linux environment. Download and extract CosmoLib package:

Here x.x should be replaced with the version number.

Compiler and library options are defined in the file configure.in, which is self-explanatory. The minimal requirement is that you need to have a fortran compiler, a C compiler and a preprocessor. You can obtain the GNU Compiler Collection gcc from http://gcc.gnu.org. It contains everything you need.

The default compilers are gfortran for Fortran and gcc (with preprocessor cpp) for C. If you use other compilers, you need to change the lines "FORTRAN_COMPILER = ..." and "C_COMPILER = ...". The flags are automatically set for gfortran + gcc + cpp and ifort + icc + fpp. For other compilers you should define the compiler flags.

There are some linking-to-library options in configure.in. If you have no idea how to link to these libraries, just set all of them to be null. CosmoLib can be compiled without any library. However, including some library, for example LAPACK, may improve the performance of the code.

Setting program options

Program options are defined in the file /headfiles/options.h, which, I hope, is also self-explanatory. At this moment you don't have to change anything there.

Compiling the entire package

To compile everything, run the shell script:

If you have changed the program options, you need to recompile. You run

The shell script clear.sh removes all compiled object files and executable files, except for those in the directory utils/. If you have modified some source files in utils/, you need to run

Go to the folder forecast/.

By default the compiled executable program is "cosm", run it:

You see an error message:

The story of the inifile

If you have experience using CosmoMC, you would know what the "inifile" is. It is a text file defining the program parameters. Let's take an example: the file forecast/lcdm.ini. Use your favorite text file editor to open it. (I assume emacs.)

Briefly read the file lcdm.ini. A line starting with # is a comment line. Most of the content is just specifications of the CMB (Cosmic Microwave Background), LSS (large scale structure, here referring to the spectroscopic galaxy surveys), SN (type Ia supernova) experiments. At this moment let's not worry about the details. We focus on the "action = ..." entry and "fidfile = ..." entry.

Set the entry action = -1, which, as explained in the comments, produces the mock CMB, LSS, SN data. Make sure you have set simulate_CMB = T (here T stands for true, and F stands for "false"), simulate_LSS = T and simulate_SN = T. Now let's just start running...oh, no, wait a second, where are the fiducial cosmological parameters? OK! If you have this question in your mind, you know what you are doing! The fiducial cosmological parameters are defined in another file called "fidfile". Notice there is an entry fidfile = lcdm.fid in the inifile lcdm.ini.

The story of the fidfile

Open the fidfile lcdm.fid.

Again all comment lines start with #. You find the first non-comment line:
data/LCDM.
This defines the root of your mock data files. The root can be either a relative path w.r.t. where you run the program, or an absolute path. Here "data/LCDM" is a relative path pointing to the path "forecast/data/LCDM...", because we are running the program in the directory forecast.

After running the program, the following files will be produced:
root_CMB.dat (in our example it is forecast/data/LCDM_CMB.dat)
root_LSS.dat (in our example it is forecast/data/LCDM_LSS.dat)
root_SN.dat (in our example it is forecast/data/LCDM_SN.dat)
root_CMB.Cls (in our example it is forecast/data/LCDM_CMB.Cls. This file saves l(l+1)C_l/(2\pi)T_cmb^2 in unit of \mu K^2. Note that CosmoLib and CAMB use different conventions for the sign of C_l^{TE}.)

You can change the root in the fidfile of course.

The second non-comment line is a line with bunch of numbers. These numbers are the fiducial parameters.

What exactly are these parameters? If you have set USE_COSMOMC_THETA = YES in the file headfiles/options.h before you compiling the code (still remember?), these parameters are Ω_bh², Ω_ch², θ, τ_re, ...; The complete list of indices of parameters are defined in headfiles/options.h (where you define the number of parameters in each "block"), headfiles/define_MCMC_params.h (where you define the actual index of each parameter), and forecast/PassParams.f90 (the code that passes the parameters). But let's not worry about the details for now. The point here is to let you know where fiducial parameters are defined.

The 3rd, 4th, and 5th non-comment lines in the fidfile are the lower-bounds, upper-bounds and step sizes of the parameters. You can fix a parameter by setting the step size to be zero, or by setting the upper-bound equal to the lower-bound. Parameters with upper-bound > lower-bound and step size > 0 are the varying parameters. (There is a special case where you can use negative step sizes, but let's not worry about the details for now.)

Producing the mock data

Having set action = -1 in the inifile and the fiducial parameters in the fidfile, you can now run

This time it should work. You see a lot of feedback information on the terminal. After about 10 seconds (if your laptop is not much worse than mine), the mock data are produced. Go to the directory forecast/data/. You should see the mock data files LCDM_CMB.dat, LCDM_LSS.dat, LCDM_SN.dat, and LCDM_CMB.Cls.

Produce the fisher matrix

The quick way to do forecast is to use the Fisher matrix. I assume you know everything about it. My task is to explain to you how to compute it with CosmoLib.

Change the entry action = -1 to action = 2 in the inifile lcdm.ini. As explained in the comments, "action = 2" means calculating the fisher matrix. The entry
output = chains/test
in the inifile defines the root of output files (still remember that the root of mock data is defined in fidfile?).

After running the program you will get files:
root_fisher.info (in our example it is forecast/chains/test_fisher.info)
root_CMB.fisher (in our example it is forecast/chains/test_CMB.fisher
root_LSS.fisher (in our example it is forecast/chains/test_LSS.fisher
root_SN.fisher (in our example it is forecast/chains/test_SN.fisher

OK. Now let's just run it (make sure you have action = 2 in the inifile):

After about a few minutes, you should get the result.

Go to the directory forecast/chains and check out the outputs.

test_CMB.fisher is just the CMB fisher matrix of the "varying parameters" (remember? step-size > 0 in the fidfile). Similarly test_LSS.fisher and test_SN.fisher are the LSS and SN fisher matrices.

A more useful file is test_fisher.info containing the indices of the varying parameters, the fiducial values of the varying parameters, the standard deviation of the varying parameters, and the covariance matrix (the inverse of the sum of all fisher matrices).

Fisher matrix : CMB + LSS v.s. CMB + SN

Now a challenging question: how do you compare CMB + LSS with CMB + SN?

Open the inifile lcdm.ini. Keep the action = 2 entry unchanged. You modify the entry simulate_LSS = T to simulate_LSS = F. Modify the entry output = chains/test to output = chains/test_CMBSN. Save the inifile and run

You obtain the following files:
forecast/chains/test_CMBSN_CMB.fisher
forecast/chains/test_CMBSN_SN.fisher
forecast/chains/test_CMBSN_fisher.info

This time the LSS fisher matrix file test_CMBSN_LSS.fisher is not produced, because you have set simulate_LSS = F in the inifile. The test_CMBSN_fisher.info contains the standard deviations and covariance matrix with only CMB + SN.

Now you repeat the same thing. This time you switch off SN by setting simulate_SN = F in the inifile. Do not forget to switch on LSS by setting simulate_LSS = T. Change the root output = test_CMBLSS. Run the program

You obtain
chains/test_CMBLSS_CMB.fisher
chains/test_CMBLSS_LSS.fisher
chains/test_CMBLSS_fisher.info

Producing the chains

Set action = 0 in the inifile forecast/lcdm.ini. Set the root output = chains/test. And set simulate_CMB = T, simulate_LSS = T, simulate_SN = T (or set to F if not wanted).

Remember in configure.in you have decided whether to use MPI or not. If you have MPI, you can generate 4 chains by running

The out chains will be
forecast/chains/test_1.txt
forecast/chains/test_2.txt
forecast/chains/test_3.txt
forecast/chains/test_4.txt

It will take a while to have the chains converged. By "a while" I mean at least a few hours, sometimes can be a few days. CosmoLib is slower than CAMB, because it is written in Newtonian gauge. Also the C_l integration scheme is different (optimized for primordial power spectra with sharp glitches, but slower in normal cases).

It will be significantly faster if you can run the program on a computer cluster with MPI and OPENMP together. For example, you can launch the program on 32 processors to produce 4 chains. Each chain is produced by 8 processors.

For example, on CITA computer cluster you can use the perl script "forecast/runMPI.pl". You submit the job with command

That will use 32 processors but only producing 4 chains.

If you want to run MCMC on your poor little laptop without MPI installed, just run

Only one chain will be produced:
forecast/chains/test_1.txt

Warning: Running MCMC on your laptop is not recommended. It is like frying your CPUs for many hours. Do not blame me if your laptop gets burned.

Postprocessing the chains

Chains have the same format as that in CosmoMC. You can analyze them using "getdist" in CosmoMC. However, as I claimed, CosmoLib contains everything you need. You can even make figures with only fortran and c codes in CosmoLib.

Now go to the directory postprocess/.

Open the file pcparams.ini. It is a inifile for chain postprocessing.

You set the entry chain_name = ../forecast/chains/test. This is the root of the chains. Note that the relative path is different from the one you set in forecast/lcdm.ini. The relative path should be w.r.t. the directory where you run the postprocessing program, which is postprocess/.

Edit the inifile pcparams.ini. Define the symbols for parameters
param1 = \Omega_b h^2
param2 = \Omega_c h^2
... CosmoLib knows how to deal with simple latex expressions.

You also want to plot 2D contours, for example, for n_s (the 11th parameter) and r (the 15th parameter). So you set
index_1_2d = 11
index_2_2d = 15
Set the color of the figure by setting, for example,
twod_color = skyblue

Finally, you set the output root to be plot_dir = plots/. (Do not forget the last slash in the path.)

Now run the postprocessing program:

Go to the directory postprocess/plots/ to check out the outputs.

postprocess/plots/test_1D.eps
Figure: 1D marginalized posterior probabilities.

postprocess/plots/test.tex
Table: 1D marginalized posterior probabilities.

postprocess/plots/test_p11_p15_2D.eps
Figure: 2D contours for n_s and r.

The way to compute transfer functions is extremely easy. Open your inifile (here I assume lcdm.ini). Set
action = 3
and set the entries "transfer_kmin, transfer_kmax, transfer_numk, transfer_zmin, transfer_zmax, transfer_num_redshift". For example, you want to see the redshift evolution of δ_c, δ_b, v_c,v_b for a given wavenumber k=0.1 h Mpc^-1. You set

transfer_kmin = 0.1
transfer_kmax = 0.1
transfer_num_k = 1
transfer_zmin = 0
transfer_zmax = 1000.
transfer_num_redshift = 100

The output directory and root are defined in your inifile (output = ***). The filename is "***_transfer_fixed_k.dat". The meaning of each column are explained in the output file itself.

If you want the transfer functions at a given redshift (let's say redshift z=1), but for a range of wave numbers (let's say 0.001 Mpc^-1<k/h<1 Mpc^-1). You set

transfer_kmin = 0.001
transfer_kmax = 1
transfer_num_k = 200
transfer_zmin = 1.
transfer_zmax = 1.
transfer_num_redshift = 1

END of the Basic Tutorial