Software | Cosmology | Advanced Topics | Forthcoming |
---|---|---|---|
|
I assume you are using LINUX or MacOS. In a terminal you run
git clone https://github.com/zqhuang/COOP |
Once you have cloned a copy from Github, enter the COOP directory:
cd COOP |
The key file for compiling COOP is configure.in. Open this file with your favorite text editor. (Here I assume emacs).
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
cp configure_linux_gfortran.in configure.in |
Now open it again
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
./clear.sh |
then
./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).
If you do not want to save any of your changes to the code, you can run
./autoupdate.sh |
cp configure.in configure_mycopy.in |
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
git status |
You will see a list of modified files. Save a copy of each modified file before you run autoupdate.sh.
To solve conflict of one single file, let's say xxx, you make a local copy of it
cp xxx xxx_copy |
pull it from github
git checkout xxx |
then manually merge the local copy and the remote copy.
Once all conflicts are resolved, you can use
git pull origin master |
to update your copy.
I am too lazy to write a software paper. You can cite this paper where I first mentioned COOP and made it publicly available.
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.
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.
Enter the firstorder folder:
cd firstorder/ |
Compile the linear order solver:
make SolvePert |
Set the parameters and output options in params.ini and run it
./SolvePert params.ini |
Enter the firstorder folder:
cd firstorder/ |
Compile the Cl solver:
make CalcCls |
Set the cosmological parameters and output options in params.ini and run it
./CalcCls params.ini |
Enter the firstorder folder:
cd firstorder/ |
Compile the Matter Power solver:
make CalcMP |
Set the cosmological parameters and output options in params.ini and run it
./CalcMP params.ini |
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
./clear.sh |
./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
cd forecast/ |
./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:
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.
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:
cd forecast/ |
Run FISHER (if it does not exist, run "make FISHER" first)
./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
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
./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
../utils/fasy.sh myfigure.txt |
, in which case you get a pdf figure myfigure.pdf; or run pyplot
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
./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:
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
./GetPeaks -map map2.fits -out map2peaks.dat |
The next step is to stack "map1.fits" at the peaks. Run
./Stack -map map1.fits -peaks map2peaks.dat -field T -out stackedfig |
./Stack -map 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./MSMAP |
map1.fits |
I2TQTUT |
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:
cd mapio/ |
make XFASTER |
./XFASTER spider.ini |
MCMC chains postprocessing
[back to top]
documentation TBD