CosmoMC Readme

July 2014. Check the web page for the latest version.


Contents

Introduction

CosmoMC is a Fortran 2003/2008 Markov-Chain Monte-Carlo (MCMC) engine for exploring cosmological parameter space, together with code for analysing Monte-Carlo samples and importance sampling (plus a suite of python scripts for building grids of runs and plotting and presenting results). The code does brute force (but accurate) theoretical matter power spectrum and Cl calculations with CAMB. See the paper for an introduction, descriptions, and typical results from some pre-WMAP data. 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.

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 fortran programs supplied cosmomc and getdist. The first 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 "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). 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

You will need a Fortran 2003 (or higher) compiler - Intel Fortran 14 or higher works (earlier versions probably will not). Others may, but gfortran does not support Fortran 2003 well enough at the moment (see gfortran bug list if you'd like to help; bugs are being worked on and hope to have gfortran support in a few months).

If not using Intel also need to link to LAPACK (for doing matrix diagonalization, etc) - you may need to edit the Makefile to specify where this on your system (Intel compilers use MPL).

Using Visual Fortran there's no need to use the Makefile, just open the project file in the source folder, and set params.ini as the program argument. This is not set up to compile with MPI, so mainly useful for development.

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). If necessary edit the job_script file as appropriate for your machine, and also if necessary the submitJob function in python/jobQueue.py.

    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 batch1/ directory, set up for Planck/BICEP analysis 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

    Analysing samples and plotting

    The getdist program analyses text files produced by the cosmomc program. These 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). The getdist program could be used completely independently of the cosmomc program.

    Run getdist distparams.ini to process the chains specified in the parameter input file distparams.ini. This should be fairly self-explanatory, in the same format as the cosmomc parameter input file.

    GetDist Parameters

    The .ini file comments should explain the other options.

    Output Text Files

    Plotting

    If GetDist produces scripts files to make simple 1D, 2D and 3D plots. These can be either python or matlab, set plot_ext=py or plot_ext=m for which you prefer. The script files produced are called

    To compare two different sets of chains set compare_num=1 in the .ini file, and compare1 to the root name of some chains you have previously run GetDist on.

    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 python and Matlab 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 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 determimes 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. There is also a FAQ in the CosmoloGUI readme.


    Antony Lewis.