Cosmology Object Oriented Package (COOP)



Author: Zhiqi Huang (email: huangzhq25@sysu.edu.cn)

Software Cosmology Advanced Topics Forthcoming
  • Fisher-MCMC blended forecast
  • Second-order Boltzmann code
  • Quick large scale structure simulator
  • Scalar-field lattice simulator

Download

[back to top]

I assume you are using LINUX or MacOS. In a terminal you run

linux terminal

git clone https://github.com/zqhuang/COOP

Install

[back to top]

Once you have cloned a copy from Github, enter the COOP directory:

linux terminal

cd COOP

The key file for compiling COOP is configure.in. Open this file with your favorite text editor. (Here I assume emacs).

linux terminal

emacs configure.in

You see that this file defines the compiler/compiler flags, external libraries etc. To make things easier, you can start with a version that is closest to the settings on your platform. For example, if you are using gfortran compiler on a Linux platform, you close the file configure.in and overwrite it with

linux terminal

cp configure_linux_gfortran.in configure.in

Now open it again

linux terminal

emacs configure.in

Edit the compiler options. COOP can be compiled without any external libraries. If you are not sure about how to set an option, you can have a look at configure_minimum.in, which gives the "minimal settings".


When you are done with configure.in, you run

linux terminal

./clear.sh

then

linux terminal

./compile.sh

The first shell script clear.sh removes all previously compiled objects, and the second script compile.sh compiles all the codes. You see that clear.sh is not necessary for the first-time compile.

If error occurs, you may need to upgrade your compiler to a more modern version. COOP has been tested for gcc 4.8.0 up and ifort 14.0.0 up. Early versions of compilers may not support the object-oriented features in COOP. If you think it is due to a bug of the code rather than compiler versions, please do send me an email (zhiqihuang AT gmail DOT com).


Update

[back to top]

If you do not want to save any of your changes to the code, you can run

linux terminal

./autoupdate.sh
This shell scripts synchronize your COOP copy to the latest version, except for the configure.in file. If you also want to get configure.in file updated, do
linux terminal

cp configure.in configure_mycopy.in
and
linux terminal

git checkout configure.in

Then you manually merge your local copy configure_mycopy.in and the updated configure.in.


If you do not want to overwrite your local changes, you do

linux terminal

git status

You will see a list of modified files. Save a copy of each modified file before you run autoupdate.sh.

Note:

To solve conflict of one single file, let's say xxx, you make a local copy of it
linux terminal

cp xxx xxx_copy

pull it from github

linux terminal

git checkout xxx

then manually merge the local copy and the remote copy.

Once all conflicts are resolved, you can use
linux terminal

git pull origin master

to update your copy.


Citation

[back to top]

I am too lazy to write a software paper. You can cite this paper where I first mentioned COOP and made it publicly available.


The metric and units conventions.

[back to top]

Metric: ds2 = a2(τ)[-(1+2Φ)dτ2 + (1 - 2 Ψ) dx2]

Units: time and lengths are measured with H0-1; densities are meaured with H02Mp2; wavenumbers are measured with H0.

Choosing Dark Energy Models.

[back to top]

Choose a Dark Energy model in configure.in:

DARK_ENERGY_MODEL = LAMBDA (LambdaCDM model)

DARK_ENERGY_MODEL = COUPLED_DE (Coupled CDM-DE model; see e.g. arXiv:1305.7457)

for DARK_ENERGY_MODEL = COUPLED_DE, the input contains a free function 1+w(a) (named as wp1 in the source code). It is defined as the effective background 1+w(a) in GR that gives the same expansion history. In this model Ωc is defined differently from the standard convention: it is defined as the ratio between CDM density extrapolated from early time with an a-3 law and the critical density.

DARK_ENERGY_MODEL = EFT (effective field theory Dark Energy description; based on aXiv: 1411.3712).

For DARK_ENERGY_MODEL = EFT, the input contains 6 free functions: αM(a); αT(a); αB(a); αH(a); αK(a); and 1+w(a) (named as wp1 in the source code). Here w(a) can have two different meanings. See examples in firstorder/params.ini. There is a switch "w_is_background".

w_is_background = .false.
means that I treat w as p_D/ rho_d, with p_D and rho_D defined in eq. (114) of this paper . This is the default option.

w_is_background = .true.
means I treat w as the effective w in GR which gives the same expansion history. In other words, if you keep w fixed and change alpha_M, the background evolution remains the same.

The α functions are defined in this paper.

For MCMC runs, you can add
w_is_background = T
to the ini file to enable this feature.

Solving Linear Perturbations.

[back to top]

Enter the firstorder folder:

linux terminal

cd firstorder/

Compile the linear order solver:

linux terminal

make SolvePert

Set the parameters and output options in params.ini and run it

linux terminal

./SolvePert params.ini

Computing CMB power spectrum.

[back to top]

Enter the firstorder folder:

linux terminal

cd firstorder/

Compile the Cl solver:

linux terminal

make CalcCls

Set the cosmological parameters and output options in params.ini and run it

linux terminal

./CalcCls params.ini

Computing matter power spectrum.

[back to top]

Enter the firstorder folder:

linux terminal

cd firstorder/

Compile the Matter Power solver:

linux terminal

make CalcMP

Set the cosmological parameters and output options in params.ini and run it

linux terminal

./CalcMP params.ini

MCMC and Planck 2015 likelihood.

[back to top]

To use the MCMC with Planck 2015 data ( optionally + BAO + SN + HST), you have to change configure.in and recompile the code.

################ changes in configure.in ###############

## set MPI
USE_MPI = YES
MPI_COMPILER_ALIAS = mpif90
#(in your platform it may be a different alias, e.g. mpif77)

#The code is safe for MPI/OPENMP hybrid, so you can keep
USE_OPENMP = YES

##For MCM runs with many parameters. It is strongly recommended that you use a LAPACK library to guarantee the speed and stability:
##on Mac OS you set LAPACK = APPLE, otherwise set
LAPACK = your_lapack_path

#link to Planck 2015 Likelihood (Planck 2013 likelihood has been deprecated and no longer works with the code)
CLIKPATH = XXXX/plc-2.0/
##replace XXXX with the path to your Planck 2015 likelihood installation

##################################################

That's all what you need to change in configure.in.

Now recompile the code

linux terminal

./clear.sh
and
linux terminal

./compile.sh

If this does not work out straightly, you may need to hack the compiler flags in configure.in. Most problems are caused by various conflicts between intel fortran compiler and different versions of openmpi. Things can be quite painful if the compilers on your platform are not up to date.


The MCMC code is in the folder forecast/. Try

linux terminal

cd forecast/
and
linux terminal

./MCMC myinis/std6.ini

I have a few sample ini files in forecast/myinis/. If you want to run EFT DE, look at the file forecast/myinis/eftde.ini.

For a parallel MCMC run on cluster (do not run on the laptop as it will burn your laptop), you change the ini file

############## changes in the ini file ###############

action = MCMC
##action = TEST is the same as action = 4 in cosmomc, and action = MCMC is the same as action = 0 in cosmomc
feedback = 1
overwrite = F
##for a run with checkpoint;

#Then pick up the data sets you want to use. In most cases I use the standard Planck TTTEEE + lowTEB + lensing + BAO + SN + HST:
use_CMB = T
use_SN = T
use_BAO = T
use_HST = T
use_lensing = T
cmb_dataset1 = XXXX/hi_l/plik/plik_dx11dr2_HM_v18_TTTEEE.clik
cmb_dataset2 = XXXX/low_l/bflike/lowl_SMW_70_dx11d_2014_10_03_v5c_Ap.clik
cmb_lensing_dataset = XXXX/lensing/smica_g30_ftl_full_pp.clik_lensing
#you need to replace XXXX with the path of Planck 2015 data directory (by default the data files are in the plc_2.0 folder in Planck 2015 installation).

##################################################

In the cluster I am using, each node contains two sockets, each socket contains 4 cores. So I use 4 nodes (4x2x4 = 32 cores) to run 8 parallel chains. Each chain uses 4 cores:

linux terminal

mpirun -np 8 --map-by node:PE:4 ./MCMC myinis/eftde.ini > eftde_run1.log

This can be very platform-dependent, you may need to ask the cluster system admin about how to optimally run an MPI/OPENMP hybrid code.

Fisher Forecast.

[back to top]

For fisher matrix calculation with many parameters, it is recommended that you link COOP to a LAPACK library. This is not mandatory, however, and may not even matter when you only have a few tens of parameters.

Enter the forecast directory:

linux terminal

cd forecast/

Run FISHER (if it does not exist, run "make FISHER" first)

linux terminal

./FISHER fisher_inis/fisher_lcdm.ini

The calculation typically takes a few minutes on a desktop or a fast laptop. The output files are ROOT_std.txt, ROOT_fisher.txt and ROOT_cov.txt, where ROOT is defined in the ini file that you passed to FISHER (in this example it is fisher_lcdm.ini). If ROOT contains a directory, make sure the directory exists.

************************ details about the ini file ********************

In the ini file, you specify the number of parameters and the names of parameters. Not all the parameters have to be used.

params_slow defines the parameters that affect the background or linear transfer functions. params_fast defines the parameters for primordial power spectra. All the other parameters are treated as nuisance parameters.

For each parameter xxx you must have a line:

param[xxx] = fiducial step1 step2 Gaussian_prior_width

fiducial (the fiducial value) is mandatory. step1, step2 and Gaussian_prior_width are optional. When step1 is missing, the parameter xxx is treated as a fixed parameter. Otherwise you must specify both step1 and step2. Step1 must be different from step2, and (fiducial + step1) and (fiducial + step2) must be both physical. For example, for tensor to scalar ratio r, if you have used a fiducial r = 0, step1 and step2 must be all positive. The order of magnitude of step1 and step2 should be reasonable (not too far away from the expected standard deviation), but the actual values of step1 and step2 do not matter for the fisher matrix calculation. You should find essentially no effect on the final results if you increase both step1 and step2 by a factor of 2. Finally, if Gaussian_prior_width is given, a Gaussian prior will be applied to this parameter.

In the ini file you also need to define the number of observations and an ini file for each observation, which defines the configurations of each observation. You can read these ini files for more details.

************************ details about the output files ********************

ROOT_std.txt gives the standard deviations of each parameter.

ROOT_fisher.txt contains: 1st line, the list of parameters (fixed parameters or unconstrained parameters are skipped); 2nd line to (n+1)-th line: the n x n fisher matrix of the parameters listed in the 1st line. Each line contains a row of the matrix.

ROOT_cov.txt contains: 1st line, the list of parameters (fixed parameters or unconstrained parameters are skipped); 2nd line, the fiducial values of parameters listed in the 1st line; 3rd line to (n+2)-th line: the n x n covariance matrix (inverse of fisher matrix) of the parameters listed in the 1st line. Each line contains a row of the matrix.

************************* defining your own model of primordial power spectrum ***********************************************

Let us assume that you have a model:

Ps(k) = As (k/kpivot)(ns-1) { 1 + μ cos[ln(k/kpivot) + φ]}

Pt(k) = r As

In COOP kpivot is fixed to be 0.05 Mpc-1 by default. You have five parameters: As, ns, μ, φ and r. In COOP the user-defined models are parameterized with up to 10 parameters: user_pp1, user_pp2, ..., user_pp10. You can use user_pp1, user_pp2, user_pp3, user_pp4 and user_pp5 for As, ns, μ, φ and r, respectively.

Open the file include/user_defined_primordial_power.h. The user-defined models are labeled as model 101, model 102, ....; You can edit model 101 by modifying the subroutine:

subroutine coop_user_defined_primordial_power_101(kbykpiv, ps, pt, cosmology, args)
COOP_REAL::kbykpiv, ps, pt
type(coop_cosmology_firstorder)::cosmology
type(coop_arguments)::args
!!these lines are entered by you
!!you define some macros to make the code more readable
#define SCALAR_AS USER_PARAM(1)
#define SCALAR_NS USER_PARAM(2)
#define MU USER_PARAM(3)
#define PHI USER_PARAM(4)
#define TENSOR_R USER_PARAM(5)
ps = SCALAR_AS * kbykpiv ** (SCALAR_NS - 1.) * (1. + MU * cos(log(kbykpiv)+PHI))
pt = TENSOR_R * SCALAR_AS
!!undefine the macros so that they don't mess up with other models
#undef SCALAR_AS
#undef SCALAR_NS
#undef MU
#undef PHI
#undef TENSOR_R
!!end of your modification
end subroutine coop_user_defined_primordial_power_101

Once this is done, recompile COOP by running clear.sh and compile.sh in the COOP directory.

The sample ini file fisher_inis/fisher_modified_pk.ini shows how to use your parametrization. First you need to specify that you want to use model 101:

pp_genre = 101

Then you put fiducial values, step1, step2 for each of the parameters user_pp1 ( = As), user_pp2 (= ns), user_pp3 (=μ), user_pp4 (= φ), and user_pp5 (=r).

param[user_pp1] = ...
param[user_pp2] = ...
param[user_pp3] = ...
param[user_pp4] = ...
param[user_pp5] = ...

For the other unused user_pp6, user_pp7 ... just fix them to be zero (or whatever values).

You are all set! You can run ./FISHER fisher_inis/fisher_modified_pk.ini and get your result.

***************************** modifying the window functions ****************

In most cases because there are very few modes at the survey scale (kmin), the smearing of survey window is not an important effect. For power spectrum with sharp features, however, it is not the case. You may want to define your own window function for a forecast. This is very simple. Open forecast/fisher_inis/params_mpk_with_windows.ini and read the comments there.

**************************** visualizing the constraints with 2D contours ***********

If you have asymptote or pyplot installed, you can visualize the fisher forecast with COOP program PLOTF. Say you ran FISHER and obtained fisher_out/lcdm_Planck_cov.txt, you want to plot the 2D constraint on Ωbh2 (named ombh2 in your ini file) and Ωch2 (named omch2 in your ini file), you do

linux terminal

./PLOTF -out myfigure.txt -xvar ombh2 -yvar omch2 -cov1 fisher_out/lcdm_Planck_cov.txt -color1 black -xlabel '$\Omega_b h^2$' -ylabel '$\Omega_c h^2$'

You should get a "COOP figure" in ascII text form: myfigure.txt. To convert it to a pdf figure, you need to run asymptote

linux terminal

../utils/fasy.sh myfigure.txt

, in which case you get a pdf figure myfigure.pdf; or run pyplot

linux terminal

python ../utils/pypl.py myfigure.txt myfigure.pdf

and get a pdf output as well.

PLOTF can plot multiple contours in one single figure and support many other features. Run

linux terminal

./PLOTF

without any argument. It will print an inline documentation.

COOP objects: list, dictionary, function, figure, file...

[back to top]

documentation TBD

Making a simple plot with COOP only.

[back to top]

documentation TBD

Making plots with COOP+Asymptote.

[back to top]

documentation TBD

Making plots with COOP+pyplot.

[back to top]

documentation TBD

Genereal-purpose MCMC.

[back to top]

documentation TBD

CMB map analysis.

[back to top]

Modify configure.in and compile COOP with cfitsiio and Healpix.

Enter the mapio folder:

linux terminal

cd mapio/

Suppose we want to stack a temperature map "map1.fits" on the peaks of "map2.fits".

The first step is to find the peaks in "map2.fits". The program "GetPeaks" does this job. Run

linux terminal

./GetPeaks -map map2.fits -out map2peaks.dat
Now the peaks of "map2.fits" are saved in file "map2peaks.dat".

The next step is to stack "map1.fits" at the peaks. Run

linux terminal

./Stack 0map map1.fits -peaks map2peaks.dat -field T -out stackedfig
If map1.fits is a IQU map or QU map and you want to do the usual Q_r-U_r stacking, you can do
linux terminal

./Stack 0map map1.fits -peaks map2peaks.dat -field QrUr -out stackedfig

The output "stackedfig.txt" is a COOP format figure. Read this documentation to understand how to convert a COOP figure to other figure formats.

You can run "./GetPeaks" and "./Stack" without arguments to get help information. Very often you want to process the maps before stacking them. The most powerful tool is "MSMAP". For example, if you want to convert a temperature map "map1.fits" to "T, Q_T, U_T" maps, run
linux terminal

./MSMAP
After the prompt "enter input file and press enter", enter
linux terminal

map1.fits
and press enter. Then enter
linux terminal

I2TQTUT
A map "map1_conv_TQTUT.fits" is produced.

CMB power spectrum analysis.

[back to top]

For CMB power spectrum analysis you need to compile COOP with cfitsio and Healpix libraries.

Enter the mapio folder:

linux terminal

cd mapio/
Compile XFASTER (if it does not exist yet):
linux terminal

make XFASTER
To run XFASTER you need data maps, fiducial model power spectrum, signal simulation maps, and noise simulation maps. The file spider.ini gives a good example about how to run XFASTER. For a full run, you first set action = DO_MASK and run
linux terminal

./XFASTER spider.ini
XFASTER generates masks. You then set action = DO_ALL and run XFASTER again. The output "COOP figures" can be converted to pdf files using asymptote or pyplot.

MCMC chains postprocessing

[back to top]

documentation TBD