------------------
==== ModeCode ====
------------------

ReadMe for modecode_mn -- MultiNest enabled version.

June 2012.

M.J. Mortonson, H.V. Peiris, and R. Easther, arXiv:1007.4205
R. Easther and H. Peiris, arXiv:1112.0326
J. Norena, C. Wagner, L. Verde, H. Peiris, and R. Easther, arXiv:1202.0304

In this ReadMe "ModeCode" refers to the Multinest enabled version of 
ModeCode.

ModeCode computes the primordial scalar and tensor power spectra for 
models of inflation. It can be run as a standalone code, or combined 
with CAMB to compute CMB angular power spectra and CosmoMC+Multinest to 
perform likelihood analysis and constrain parameters of inflationary 
models.

ModeCode is integrated with the nested sampling package MultiNest, and 
permits both parameter estimation and the computation of Bayesian 
evidence.

The code can be downloaded at:
http://zuserver2.star.ucl.ac.uk/~hiranya/ModeCode/ModeCode/

Multinest can be obtained here:
http://ccpforge.cse.rl.ac.uk/gf/project/multinest/

CosmoMC can be obtained here:
http://cosmologist.info/cosmomc/


Version history:
----------------
June 2012 - added Hubble slow roll reconstruction to the choice of
  potentials, see arXiv:1202.0304; updated to Jan2012 version of 
  cosmomc/CAMB.
December 4, 2011 - updated to interface with MultiNest, adds priors
  on post-inflationary dynamics, computation of Bayesian Evidence.
February 17, 2011 - fixed bug with output of derived parameters when 
  sampling from priors, added params_modpk_F.paramnames for running code
  with use_modpk=F. (available as modecode.tar.gz)
February 10, 2011 - updated to January 2011 version of CAMB.
July 21, 2010 - initial release.


The currently available files are:

 (1) modecode_mn.tar.gz, including the following directories:

   - cosmomc/ - sample files for running CosmoMC with ModeCode + Multinest

   - cosmomc/camb/ - sample driver for standalone code 
     (driver_modpk.f90), new files used by ModeCode, and modified CAMB files

   - cosmomc/source/ - modified CosmoMC code

   - cosmomc/multinest/ - multinest source

   - Models/*/ - files for running CosmoMC with the Multinest sampler 
     for the 7 models included in ModeCode (note that additional models 
     can be easily added as described below), where * = m2phi2 
     (quadratic potential), lphi4 (quartic), neq1 (n=1), neq2ov3 
     (n=2/3), natural, hilltop, and reconstruction (Hubble slow-roll 
     reconstruction). In addition there is a  distparams_modpk.paramnames 
     file corresponding to each model, which must be included via the 
     parameter_names variable when running getdist.

As usual, *.tar.gz files can be unzipped and extracted using the command 
'tar -xzvf *.tar.gz'.

The current version of the code is based on the January 2012 version of 
CAMB and CosmoMC and the April 2011 Multinest.

Only new or modified CAMB/CosmoMC files are provided in the downloads, so 
any applications beyond using the standalone version of ModeCode will 
require adding the downloaded files to an existing copy of the basic 
CAMB and CosmoMC codes.

In each of the modified files, changes to the original code are marked 
by 'MODIFIED P(K) ... END MODIFIED P(K)' so that they can be easily 
located.

==== Setting up and compiling ModeCode ====

After extracting files from modecode_mn.tar.gz, the ModeCode files must 
be copied into existing versions of CAMB/CosmoMC for use with those 
codes. One easy way to do this is to start in the cosmomc/ directory for 
the ModeCode files, and use 'cp -r . cosmomc_target/' where 
cosmomc_target is the base directory for the unmodified code.

*** Note that this will overwrite original versions of some of the 
CAMB/CosmoMC files; if you do not want this to happen, rename the files 
first before copying (they will also need to be renamed in the 
Makefiles). ***

Once the new and modified code files have been copied into the 
appropriate directories, the ModeCode versions of CAMB and CosmoMC can 
be compiled with the usual Makefiles, after adjusting the compiler options 
in the camb/Makefile and source/Makefile supplied with ModeCode to work
on your system. You will also need to supply any likelihood codes (e.g. 
WMAP) that do not ship with the standard CAMB/CosmoMC.

Executing the make command in the source directory will produce the 
executables cosmomc_multinest AND cosmomc.  The cosmomc_multinest must 
be used to compute evidence with the multinest sampler.

To compile the sample driver program for the standalone version of 
ModeCode without compiling CAMB or CosmoMC, run 'make powspec' in the 
camb/ directory.

==== General description of ModeCode parameters and options ====

Inflationary models are described in ModeCode using an array of 
parameters for the potential V(phi), vparams(i) for i = 1, 2, etc., and 
the parameter N_pivot (not used for HSR models, see below), which 
accounts for the uncertainty in the post-inflationary epoch between the 
inflation and reheating. N_pivot is the number of inflationary e-folds 
remaining when the mode of a given scale infl_pivot_k exits the horizon. 
The value of infl_pivot_k can be specified in units of Mpc^{-1}.



-----------------
ModeCode options:
-----------------

- use_modpk (logical): a general flag to turn ModeCode on or off. 
Setting use_modpk=F recovers the original functionality of CAMB/CosmoMC 
while retaining the Multinest sampler. When running the code with 
use_modpk=F, ParamNamesFile in params_modpk.ini should be changed to 
params_modpk_F.paramnames.

- modpk_physical_priors (logical) Excludes models where the scalar
  spectral amplitude is vastly different from the observed
  value. As implemented, this cut excludes models who which
  10^-11<A_s<10^-7

- modpk_rho_reheat (real): Reheating **density** in GeV^4 Note that
  the temperature / energy ~ modpk_rho_reheat^(1/4)

- modpk_w_primordial_lower (real): Lower limit on effective equation
  of state parameter during between the end of inflation and
  rho=rho_reheat

- modpk_w_primordial_upper (real): Upper limit on effective equation
  of state parameter during between the end of inflation and
  rho=rho_reheat

- potential_choice (integer): specifies which inflationary model to use.
  The translation between potential_choice and V(phi) is set in
  camb/modpk_potential.f90. The models corresponding to different values
  of potential_choice by default are listed below.

- vparams_num (integer): number of parameters needed to describe V(phi)
  (for general reheating models, the total number of MCMC parameters is
  vparams_num+1). The maximum number of parameters is set to 9 by 
  default, but can be increased in camb/modpk_modules.f90 by changing 
  max_vparams.

- phi_init (real): starting guess for the initial value of phi; for all 
  models provided with the code, this setting is superceded by the 
  'initialphi' function in camb/modpk_potential.f90.

- infl_pivot_k (real): ``pivot'' wavenumber for evaluating N_pivot, in Mpc^{-1}.
  Note that this parameter is named k_pivot, when ModeCode is called with CAMB
  instead of CosmoMC/Multinest.

- N_pivot (real): number of e-folds of inflation remaining after the
  mode with wavenumber infl_pivot_k leaves the horizon (note that 
  N_pivot is treated as a parameter to be varied rather than a fixed 
  setting).

- instreheat (logical): whether or not to assume instant reheating, 
  which fixes the value of N_pivot; if instreheat=T, the values chosen for 
  infl_pivot_k and N_pivot are ignored.

- slowroll_infl_end (logical): whether or not to determine when 
  inflation ends by the breakdown of slow roll conditions, i.e. 
  epsilon_H=1. If slowroll_infl_end=F, the end of inflation is assumed 
  to occur when phi=phi_infl_end.

- phi_infl_end (real): final value of phi during inflation for models 
  that do not end via slow roll violation (slowroll_infl_end=F). If 
  slowroll_infl_end=T, the value of this parameter is ignored.

- vnderivs (logical): whether to use numerical derivatives of the 
  potential (vnderivs=T) or analytic forms supplied in 
  camb/modpk_potential.f90 (vnderivs=F). The latter option is STRONGLY 
  recommended for all models for which the derivatives of V(phi) can be 
  computed and expressed analytically, as the use of numerical derivatives 
  may lead to inaccurate results for certain models.

- action (integer) action = 3 to use multinest, action=0 to use the
  standard MCMC sampler (the use of the MCMC sampler with this version
  of ModeCode has not been tested, and is currently unsupported)

---------------------------
Default models in ModeCode:
---------------------------
- potential_choice = 1: quadratic
      V(phi) = m^2 phi^2 / 2
      vparams(1) = log_10(m^2)
- potential_choice = 2: natural
      V(phi) = Lambda^4 [1+cos(phi/f)]
      vparams(1) = log_10(Lambda), vparams(2) = log_10(f)
- potential_choice = 3: quartic
      V(phi) = lambda phi^4 / 4
      vparams(1) = log_10(lambda)
- potential_choice = 4: linear
      V(phi) = lambda phi
      vparams(1) = log_10(lambda)
- potential_choice = 5: exponent n=2/3
      V(phi) = (3/2) lambda phi^{2/3}
      vparams(1) = log_10(lambda)
- potential_choice = 6: hilltop
      V(phi) = Lambda^4 - lambda phi^4 / 4
      vparams(1) = log_10(Lambda), vparams(2) = log_10(lambda)

--------------------------------------
Using ModeCode for HSR reconstruction:
--------------------------------------
J. Norena, C. Wagner, L. Verde, H. Peiris, R. Easther, arXiv:1202.0304

In the Hubble slow-roll (HSR) reconstruction the Hubble parameter is 
modeled by a finite polynomial:
H(phi) = H_star (1 + A_1 phi + A_2 phi^2 + ... + A_N phi^N).

The coefficients A_i are related by a one-to-one correspondence to the 
Hubble slow-roll parameters (epsilon, eta, xi, ...), and H_star sets the 
overall energy scale.

The corresponding potential is then given by:
V(phi) = M_pl^2 H_star^2 [3(1 + A_1 phi + ... + A_N phi^N)^2
         - 2(A_1 + ... + N A_N phi^(N-1))]

At the moment, the HSR potential is implemented in ModeCode for the 
first three HSR parameters, i.e. epsilon, eta and xi. However, the 
extension to higher orders is straightforward. To allow for uniform and 
log priors on epsilon, one can choose between:

- potential_choice = 7: HSR with eps, eta, xi
      vparams(1) = epsilon
      vparams(2) = eta
      vparams(3) = xi
      vparams(4) = log(10^10 A_SR)

 and

- potential_choice = 8: HSR with eps, eta, xi
      vparams(1) = log_10(epsilon)
      vparams(2) = eta
      vparams(3) = xi
      vparams(4) = log(10^10 A_SR)

where the HSR parameters are given at phi=0. A_SR is the amplitude of 
the curvature power spectrum at the pivot scale (i.e., the mode which 
exits the horizon at phi=0) computed at second order in the HSR 
parameters. This relation between the HSR parameters and the amplitude 
is then used to set the value of H_star.

Note that N_pivot is not used. Instead the parameter 
reconstruction_Nefold_limit specifies the minimum number of inflationary 
e-folds after the scale given by infl_min_k exits the horizon. In 
addition, the smallest scale for which one requires that it exits the 
horizon still during inflation is given by infl_max_k.

In summary, the parameters needed for the HSR reconstruction are:

- infl_pivot_k (real): specifies the mode which exits the horizon when 
      phi=0, in Mpc^{-1}. In addition, derived parameters like the spectral
      tilt or the tensor-to-scalar ratio etc. are evaluated at this scale.
      Note that this parameter is named k_pivot, when ModeCode is called with 
      CAMB instead of CosmoMC/Multinest.

- infl_min_k (real):  largest scale for which one requires that it exits 
      the horizon during inflation, in Mpc^{-1}.

- infl_max_k (real):  smallest scale for which one requires that it exits 
      the horizon during inflation, in Mpc^{-1}.

- reconstruction_Nefold_limit (real): minimum number of inflationary
      e-folds counted from infl_min_k.


==== Computing power spectra (standalone code) ====

The sample driver included with the code (camb/driver_modpk.f90) can be 
run using the command 'powspec' in the camb/ directory. It is set up to 
compute the scalar and tensor power spectra at 500 k values between 
5x10^{-4} Mpc^{-1} and 5 Mpc^{-1} for a natural inflation model. The 
code outputs (k, P_s(k), P_t(k)) and also computes the scalar and tensor 
amplitudes and spectral tilts, as well as the running of the scalar 
spectral index, at k_pivot = 0.05 Mpc^{-1}. This program can be easily 
modified to compute the spectra for different models of inflation and/or 
different values of k.

Warnings that phi_init is inconsistent and the value of phi_init is 
being rescaled are a normal byproduct of the algorithm that searches for 
self-consistent initial conditions for inflation. If there are an 
excessive number of these warnings for each model evaluated, it may help 
to change the user-supplied phi_init value or the function used in 
'initialphi' in camb/modpk_potential.f90 for the inflationary model in 
question.

Other warnings or errors produced by ModeCode typically indicate that 
the chosen model parameters do not have a physically acceptable 
inflationary solution.

==== Using ModeCode with CAMB ====

ModeCode can be used within CAMB by running 'camb params_modpk.ini' in 
the camb/ directory. The ModeCode options and parameters are set in the 
section of params_modpk.ini marked 'MODIFIED P(K)', which follows the 
entries for the initial power spectrum parameters from the unmodified 
version of CAMB. Note that these original spectral parameter values 
(e.g. scalar_amp(1), scalar_spectral_index(1), etc.) are ignored if 
use_modpk=T since the initial power spectra are entirely specified by 
the values of N_pivot and the vparams array.


==== Using ModeCode with CosmoMC+Multinest ====

The use of the ModeCode version of CosmoMC+Multinest is largely 
unchanged from the original version of the code.

ModeCode-specific options have been added to the params_modpk.ini files, 
and which these are are marked by 'MODIFIED P(K)'.

When CAMB+ModeCode is run as a standalone code, N_pivot is a free 
parameter specified by the user. Within CosmoMC+Multinest N_pivot is an 
independent parameter.  If instant reheating is assumed (instreheat=T), 
then N_pivot should be fixed by setting 'param[Npivot] = 50 50 50 0 0' 
in params_modpk.ini (the choice of 50 here is unimportant since the code 
ignores this value for instant reheating models).

N_pivot and the vparams array (named vpar1, vpar2, etc. in CosmoMC) are 
added on the end of the list of the usual CosmoMC parameters, following 
A_SZ. Elements in the vparams array that are not used by a particular 
inflationary model should be fixed in the MCMC analysis. If max_vparams 
in camb/modpk_modules.f90 is changed to allow for models with more than 
9 V(phi) parameters, the value of num_vpar in source/settings.f90 should 
also be changed so that num_vpar=max_vparams+1.

The chain files output by the ModeCode version of CosmoMC have several 
additional derived parameters which are listed at the end of 
params_modpk.paramnames:

  - modpk_Npivot; if instreheat=F, this should always be equal to the 
    chain parameter N_pivot, but if instreheat=T it will be the value of 
    N_pivot computed by the code to satisfy the matching equation for 
    instant reheating models. Note that for HSR reconstruction this value
    is meaningless.

  - modpk_ns, modpk_nt (scalar and tensor spectral tilt)

  - modpk_nrun (scalar spectral running dn_s/dlnk)

  - modpk_logA (ln(10^{10}A_s), i.e. the usual CosmoMC scalar amplitude 
    parameter)

  - modpk_r (tensor-to-scalar ratio)

  - modpk_w Effective primordial equation of state parameter (see 
    arXiv:1112.0326 for details)

Each of these spectral parameters is computed directly from the 
primordial power spectra computed by ModeCode at the pivot scale.

==== Using ModeCode with GetDist ====

Output files can be processed with GetDist, but the appropriate 
distparams.ini file must have the parameter_names variable set to 
distparams_modpk.paramnames from the corresponding subdirectory within 
/Models. 

The Multinest-enabled version of ModeCode only writes out parameters 
which are varied by the sampler, followed in sequence by the derived 
parameters in subsequent columns. The distparams_modpk.paramnames file 
must therefore contain only this set of parameters (i.e. leaving out 
parameters which were not varied by the sampler) in order for the MatLab 
plots produced by GetDist to have the correct axis labels. Only a single 
set of examples with distparams_modpk.paramnames is given for each model 
within ModeCode - the user might need to modify the 
distparams_modpk.paramnames file if varying a different subset of 
parameters than this default set.

In summary:

- take the disparams_modpk.paramnames file provided in Models/POTENTIAL_TYPE/

- delete any lines which are in there corresponding to parameters which 
you have not varied - leave the derived parameters (flagged by *) untouched.

- point distparams.ini at the modified disparams_modpk.paramnames file 
using the parameter_names flag.

==== Adding new inflationary models ====

Single field inflationary models beyond those provided can be simply 
added to ModeCode by adding the following functions to 
camb/modpk_potential.f90:

  - V(phi) in the function pot(phi)

  - the first derivative of V in the function dVdphi(phi)

  - the second derivative of V in the function d2Vdphi2(phi)

  - an approximate expression for phi(N_pivot) (e.g., derived using slow 
    roll relations) in the function initialphi; this is used to compute a 
    reasonable first guess for the initial conditions given the shape of 
    the potential

Although only specification of V(phi) is absolutely necessary (with 
vnderivs=T the code will attempt to compute numerical derivatives of the 
potential, and the function initialphi defaults to the set value of 
phi_init if no phi(N_pivot) relation is given), the results of the code 
are generally much more reliable if all four of these functions can be 
provided.

The 8 potentials choices provided in the code can be used as templates for 
each of these functions. For example, a new form of the potential corresponding 
to potential_choice=9 would require adding a 'case(9)' statement to each 
of the four functions described above with the V(phi) parameters and 
functions specified following the examples provided.
 
==== Sampling from priors ====

In addition to computing the inflationary power spectrum we have 
modified CosmoMC to generate a uniform sample from the prior, without 
including constraints from data. This is useful, for example, for 
studying the effects of flat priors on free parameters (such as N_pivot 
and the vparams array) on the distributions of derived parameters (such 
as n_s and r). Isolating the effect of the priors helps determine which 
constraints are coming from the data and which are primarily a 
consequence of the chosen priors.

The directories Models/natural/ and Models/hilltop/ each contain an 
example file params_priors.ini with the settings needed to sample 
uniformly from the priors for these inflationary models.  This option is 
largely untested within the Multinest enabled version of CosmoMC but 
should still work using the cosmomc executable.

The main changes to the usual settings are the following:

- sampling_method = 7 (performs uniform sampling within the priors 
instead of using Metropolis-Hastings sampling)

- use_x = F for all data sets x

- all parameters that do not affect the conversion between independent 
parameters and derived parameters should be fixed; for example, if the 
derived parameters are the power spectrum parameters A_s, n_s, r, etc., 
all parameters other than the vparams array and N_pivot (for general 
reheating) should be held fixed at some value.

- additional options that are not necessary for prior sampling can be 
turned off, e.g. bbn_consistency=F, compute_tensors=F (ModeCode will 
still compute the derived parameters r and n_t in this case, as these 
are purely functions of the primordial spectra and not the CMB angular 
power spectra)

- The prior can be sampled with CosmoMC+Multinest itself. In this case 
the number of live points should be increased to at least 10,000 and 
then the code should be terminated before it starts sampling the 
likelihood.