CosmoMC Readme

June 2021. Check the web page for the latest version.


Contents

Specific documentation

Introduction

CosmoMC is a Fortran 2008 Markov-Chain Monte-Carlo (MCMC) engine for exploring cosmological parameter space, together with Fortran and python code for analysing Monte-Carlo samples and importance sampling (plus a suite of scripts for building grids of runs, plotting and presenting results). The code does brute force (but accurate) theoretical matter power spectrum and Cl calculations with CAMB. See the original paper for an introduction and descriptions, and up-to-date sampling algorithm details. It can also be compiled as a generic sampler without using any cosmology codes.

On a multi-processor machine you can start to get good results in a few of hours. On single processors you'll need to set aside rather longer. You can also run on a cluster.

You may also want to consider the new Cobaya sampler, which is based on Python and also allows use of other samplers or cosmology codes.

By default CosmoMC uses a simple Metropolis algorithm or an optimized fast-slow sampling method (which works for likelihood with many fast nuisance parameters like Planck). The program takes as inputs estimates of central values and posterior uncertainties of the various parameters. The proposal density can use information about parameter correlations from a supplied covariance matrix: use one if possible as it will significantly improve performance. Covariance matrices are supplied for common sets of default base parameters. If you compile and run with MPI (to run across nodes in a cluster), there is an option to dynamically learn the proposal matrix from the covariance of the post-burn-in samples so far. The MPI option also allows you to terminate the computation automatically when a particular convergence criterion is matched. MPI is strongly recommended.

There are two programs supplied cosmomc and getdist. The first is a compiled Fortran code does the actual Monte-Carlo and produces sets of .txt chain files and (optionally) .data output files (the binary .data files include the theoretical CMB power spectra etc.). The "cosmomc" program can also do post processing on .data files, for example doing importance sampling with new data. The "getdist" program analyses the .txt files calculating statistics and outputs files for the requested 1D, 2D and 3D plots (and could be used independently of the main cosmomc program), and there are Fortran and python versions. The Python version is is the most up to date, and output is not the same as the fortran version (which is mainly provided for backwards checking and reproducing Planck results). GetDist is included with CosmoMC, but can be also run quite separately from the standalone python package. The "cosmomc" program also does post processing on .data files, for example doing importance sampling with new data.

Please e-mail details of any bugs or enhancements to Antony Lewis. If you have any questions please ask in the CosmoCoffee computers and software forum. You can also read answers to other people's questions there.

Downloading and Compiling

To run compile and run CosmoMC you will need a Fortran 2008 compiler. Various options include If you just want to analyse chains and run plots, you only need python 2.7+ or 3.5+ (see plotting below).

You then need to download and install the program:

Note that the latest version of CosmoMC is not the same as used for the Planck 2018 final analysis. See the github planck2018 branch for the exact CAMB and other settings used for that.

Some people have reported problems using Intel MPI 2019, earlier versions should be OK.

By default the build uses MPI. Using MPI simplifies running several chains and proposal optimization. MPI can be used with OpenMP for parallelization of each chain. Generally you want to use OpenMP to use all the shared-memory processors on each node of a cluster, and MPI to run multiple chains on different nodes (the program can also just be run on a single CPU).

For code editing some project files are available for standard development environments:

See BibTex file for relevant citations.

Running and input parameters

See the supplied params.ini file for a fairly self-explanatory list of input parameters and options. The file_root entry gives the root name for all files produced. Running using MPI on a cluster is recommended if possible as you can automatically handle convergence testing and stopping.

  • MPI runs on a cluster

    When you run with MPI (compiled with -DMPI, the default), a number of chains will be run according to the number of processes assigned by MPI when you run the program. Output file names will have "_1","_2" appended to the root name automatically. You may be able to run over 4 nodes using e.g.
    mpirun -np 4 ./cosmomc params.ini
    There is also a supplied script python/runMPI.py that you may be able to adapt for submitting jobs to a job queue, e.g.
    python python/runMPI.py params
    to run using the params.ini parameters file. Run runMPI.py without arguments to see a list of options (number of chains, cores, wall time, etc). You will probably need to set up a job script template file for your computer cluster (a couple of samples are supplied). See Configuring job settings for your machine. If you want to run batches of chains with different parameter or data combinations, see the help file on setting up parameter grids.

    In the parameters file you can set MPI_Converge_Stop to a convergence criterion (see Convergence Diagnostics. Small numbers are better convergence; generally need R-1 < 0.1, but may want much smaller especially for importance sampling or accurate confidence limits. You can also directly impose an error limit on confidence values - see the parameter input file for details). If you set MPI_LearnPropose = T, the proposal density will be updated using the covariance matrix of the last half of all samples (across chains) so far.

    Using MPI_LearnPropose will significantly improve performance if you are adding new parameters for which you don't have a pre-computed covariance matrix. However adjusting the proposal density is not strictly Markovian, though asymptotically it is as the covariance will converge. The burn in period should therefore be taken to be larger when learning the proposal density to allow for the time it takes the covariance matrix to converge sufficiently (though in practice it probably doesn't matter much in most cases). Note also that as the proposal density changes the number of steps between independent samples is not constant (i.e. the correlation length should decrease along the chain until the covariance has converged). The stopping criterion uses the last half of the samples in each chain, minus a (shortish) initial burn in period. If you don't have a starting covariance matrix a fairly safe thing to do is set ignore_rows=0.3 in the GetDist parameter file to skip the first chunk of each chain.

    You can also set the MPI_R_StopProposeUpdate parameter to stop updating the proposal density after it has converged to R-1 smaller than a specified number. After this point the chain will be strictly Markovian. The number of burn in rows can be calculated by looking at the params.log file for the number of output lines at which R-1 has reached the required value.

    If things go wrong check the .log and any error files in your cosmomc/scripts directory.

    Input Parameters

    The input .ini files contain sets of name-value pairs, in the form key = value.

    Keys setting input values for specific parameters can be set using numbered key values, param[omegabh2]=xxx yy... for omegabh2. The names of the various parameters are set in .paramnames files: for cosmology parameters see params_CMB.paramnames. Likelihoods may also have nuisance parameters that need to be varied: see the .paramnames files in the ./data/ directory.

    .ini files can include other ini files using DEFAULT(xx.ini) or INCLUDE(yy.ini), so there's no need to repeat lots of settings when you do different similar runs. The sample test.ini file includes .ini files in the batch2/ directory, set up for Planck runs. The action=4 setting just calculates test likelihoods for one model (the central value of your parameter ranges) and writes out the result. Set action=0 to run a test chain. Using DEFAULT in a .ini files means you can include default settings, but override ones of interest for a particular run.
    Output files
    Configuring job settings for your machine

    The runMPI.py and runbatch.py submission scripts use a job template file to submit jobs to the queue on your cluster. The job_script files is supplied as sample, and you may need to make your own template specific to the configuration of your machine, depending on what queuing system you have (e.g. PBS, SLURM, etc.) Recommended steps are

    If necessary you can also edit the submitJob function in python/paramgrid/jobqueue.py, but most things can be configured via the template and .bashrc settings. You can change the number of chains, nodes, walltime etc when submitting a job by using optional parameters to runMPI.py and runbatch.py: see -h help for available options.

    Analysing samples and plotting

    There are now several ways you can analyse chain results: The GetDist programs analyses text chain files produced by the cosmomc program. These are a group of plain text files of the form
    xxx_1.txt
    xxx_2.txt
    ...
    xxx.paramnames
    xxx.ranges
    
    where "xxx" is some root filename, the .txt files are separate chain files, xxx.paramnames gives the names (and optionally latex labels) for each parameter, and xxx.ranges optionally specified bounds for each (named) parameter. The chain .txt files are in the format

        weight like param1 param2 param3 ...

    The weight gives the number of samples (or importance weight) with these parameters. like gives -log(likelihood), and param1, param2... are the values of the parameters at each point.

    The .paramnames file lists the names of the parameters, one per line, optionally followed by a latex label. Names cannot include spaces, and if they end in "*" they are interpreted as derived (rather than MCMC) parameters, e.g.

    x1   x_1
    y1   y_1
    x2   x_2
    xy*  x_1+y_1
    
    The .ranges file gives hard bounds for the parameters, e.g.
    x1  -5 5
    x2   0 N
    
    Note that not all parameters need to be specified, and "N" can be used to denote that a particular upper or lower limit is unbounded. The ranges are used to determine densities and plot bounds if there are samples near the boundary; if there are no samples anywhere near the boundary the ranges have no affect on plot bounds, which are chosen appropriately for the range of the samples.

    The GetDist program, scripts and gui could be used completely independently of the CosmoMC program if you produce chains (or samples) in the same format (see the pip-installable python package). If there is just one file of samples, you can also use xxx.txt without the _n number tag that's used to denote separate chains.

    Running the GetDist program

    The GetDist program processes the chains and produces a set of output files giving marginalized constraints, convergence diagnostics, etc, and optionally some scripts for producing plots. This is convenient for bulk production of parameter constraints and first-look plotting; use the python modules in your own scripts if you need more fine-grained control.

    Running using the python version

    python python/GetDist.py [distparams.ini] [chains/chainroot]
    or using compiled Fortran:
    ./getdist distparams.ini [chains/chainroot]
    Here the parameter input file distparams.ini determines analysis parameters; in the python version this is optional, and if omitted defaults will be used. The chain root name to analyse can be given in the .ini file, or you can specify it as an optional parameter on the command line to quickly analyse different chains with the same settings. The parameter input file should be fairly self-explanatory, and is in the same format as the cosmomc parameter input file. To make a new .ini file you can edit to customize settings use
    python python/GetDist.py --make_param_file myparams.ini

    Note that the outputs from the python and fortran versions of GetDist are not identical due to slightly different algorithms, but should usually be consistent to better than the sampling noise for parameter means and limits. Both can produce pre-computed densities for plotting (in the plot_data_dir directory), but the python plotting scripts do not require this - output will be faster and make far fewer files if plot_data_dir is not set.

    GetDist Parameters

    The .ini file comments should explain the other options.

    Output Text Files
    These are all output by the GetDist program, and most can also be viewed on demand by using the menu options in the
    GetDist GUI.

    Plotting

    GetDist can be used to produce basic simple scripts for various types of plots. You can also just write plotting scripts, or use the GetDist GUI, in which case there is no need to run GetDist first.

    If no_plots=F, GetDist produces scripts files to make simple 1D, 2D and 3D plots. The script files produced are called

    Parameter labels are set in distparams.ini - if any are blank the parameter is ignored. You can also specify which parameters to plot, or if parameters are not specified for the 2D plots or the colour of the 3D plots getdist automatically works out the most correlated variables and uses them. The data files used by auto-generated plotting scripts are output to the plot_data directory.

    Convergence diagnostics

    The getdist program will output convergence diagnostics, both short summary information when getdist is run, and also more detailed information in the file_root.converge file. When running with MPI the first two of the parameters below can also be calculated when running the chains for comparison with a stopping criterion (see the .ini input file).
    Differences between GetDist and MPI run-time statistics

    GetDist will cut out ignore_rows from the beginning of each chain, then compute the R statistic using all of the remaining samples. The MPI run-time statistic uses the last half of all of the samples. In addition, GetDist will use all the parameters, including derived parameters. If a derived parameter has poor convergence this may show up when running GetDist but not when running the chain (however the eigenvalues of covariance of means is computed using only base parameters). The run-time values also use thinned samples (by default every one in ten), whereas GetDist will use all of them. GetDist will allow you to use only subsets of the chains.

    Parameterizations

    Performance of the MCMC can be improved by using parameters which have a close to Gaussian posterior distribution. The default parameters (which get implicit flat priors) are

    Parameters like H0 and omegal (ΩΛ) are derived from the above. Using theta rather than H0 is more efficient as it is much less correlated with other parameters. There is an implicit prior 40 < H0 < 100 (which can be changed). The .txt chain files list derived parameters after the base parameters. The list of parameter names and labels used in the default parameterization is listed in the supplied params_CMB.paramnames file.

    Since the program uses a covariance matrix for the parameters, it knows about (or will learn about) linear combination degeneracies. In particular ln[10^10 A_s] - 2*tau is well constrained, since exp(-2tau)A_s determines the overall amplitude of the observed CMB anisotropy (thus the above parameterization explores the tau-A degeneracy efficiently). The supplied covariance matrix will do this even if you add new parameters.

    Changing parameters does in principle change the results as each base parameter has a flat prior. However for well constrained parameters this effect is very small. In particular using theta rather than H_0 has a small effect on marginalized results.

    The above parameterization does make use of some knowledge about the physics, in particular the (approximate) formula for the sound horizon. To change parameterization make a new .paramnames file, then change sources/params_CMB.f90 to change the mapping of physical parameters to MCMC array indices, and also to read in your new .paramnames file.

    Likelihoods that you use may also have their own nuisance parameters.

    Using CosmoMC as a generic sampler

    You can use CosmoMC to sample any custom function that you want to provide. To do this You can use named or un-named parameters. Making your own .paramnames file to name your parameters usually makes things easier (except possibly if you have a very large number of parameters where number labels make sense).

    Programming

    The code makes use of Fortran 2003 features and structures. There is automated code documentation, which includes modules, classes and subroutines. File and class documentation gives inheritance trees to help understand the structure, for example GeneralTypes.f90 contains the key base classes that other things are inherited from.

    You are encouraged to examine what the code is doing and consider carefully changes you may wish to make. For example, the results can depend on the parameterization. You may also want to use different CAMB modules, e.g. slow-roll parameters for inflation, or use a fast approximator. The main source code files you may want to modify are

    Add-ons and extra datasets

    Various extended data sets require some code modifications. After testing these may be included as part of the public program, but meanwhile can be installed manually or obtained from the cosmomc git repository (access available on request).

    Version History

    Reference links

    See the BibTex file of CosmoMC references you should cite (including data likelihoods), along with some references of potential interest. These are the two main CosmoMC papers and some general sampling references:

    FAQ

    1. What are the dotted lines on the plots?
      Dotted lines are mean likelihoods of samples, solid lines are marginalized probabilities. For Gaussian distributions they should be the same. For skew distributions, or if chains are poorly converged, they will not be. Sometime there is a much better fit (giving high mean likelihood) in a small region of parameter space which has very low marginalized probability. There is more discussion in the original paper.
       
    2. What's in the .likestats file produced by getdist?
      These are values for the best-fit sample, and projections of the n-Dimensional confidence region. Usually people only look at the best-fit values. The n-D limits give you some idea of the range of the posterior, and are much more conservative than the marginalized limits.
       
    3. I'm writing a paper "constraints on X". How do I add X to cosmomc?
      To add new theory parameters
      • add the definition of the physical parameter to CMBParams in CosmologyTypes.f90
      • Change CosmologyParameterizations.f90 to add a new parameterization or modify an existing one. In particular the InitWithSetNames (_Init) procedure has
        call SetTheoryParameterNumbers(num_slow_params,num_semi_slow_params)
        
        which sets the number of hard and semi-hard parameters (the latter being things like initial power spectrum parameters). Make a new .paramnames file with the new parameter names in and change and then
        call this%Initialize(Ini,Names, 'params_CMB.paramnames', Config)
        
        to set your new parameter names file.
      • Modify Calculator_CAMB.f90 (the CAMBCalc_CMBToCAMB procedure) to pass your new CMBParams parameters to CAMB to actually change the result
      • Edit your .ini files to add range and width settings for the new parameter names

      If instead your parameter is a likelihood parameter (nuisance parameter), the new parameter can be handled entirely inside your likelihood function code (see supernovae_JLA.f90 for a simple example that has two nuisance parameters): you just set the .paramnames file (which determines the number of new parameters), and then use them when passed in to your likelihood function. See Likelihood_Cosmology.f90 for likelihood classes that your likelihood can extend.
       
    4. Why do some chains sometimes appear to get stuck?
      Usually this is because the starting position for the chain is a long way from the best fit region. Since the marginal distributions of e.q. A_s are rather narrow, it can take a while for chains to move from into an acceptable region of A_s exp(-2τ). The cure is to check your starting parameter values and start widths (especially make sure the widths are not too wide-- too small is better than too wide when running with MPI), or to use a sampling method that is more robust (e.g. use sampling_method = 7). If you are patient, stuck chains should eventually find a sensible region of parameter space anyway. Occasionally the staring position may be in a corner of parameter space so that prior ranges prevent any reasonable proposed moves. In this case check your starting values and ranges, or just try re-starting the chain (a different random starting position will possibly be OK).
       
    5. How to I simulate futuristic CMB data? See CosmoCoffee.

    Feel free to ask questions (and read answers to other people's) on the CosmoCoffee software forum.


    Antony Lewis.