#Creating a modified CAMB module

In this tutorial we will walk through creating a CosmoSIS module for a modified version of CAMB, the Code For Anisotropies in the Microwave Background.

Modifying CAMB to implement new physics or new data is a fairly common task in cosmology, but is fairly involved, and so worth covering in detail. We won't discuss here the content of the actual modifications you want to make (e.g. modifying gravity, computing new observables, etc) - instead this document is about what changes you want to make to let your modified CAMB talk to CosmoSIS.

The theory

For the purposes of this tutorial we'll pretend that you are making a new version of CAMB to implement a modified cosmology based on the ~~crackpot~~ independent research Space Mirror theory, with two new input parameters:

the mirror radius
the mirror reflectivity

One new configuration option:

the reflection count limit

And two new outputs:

the theory sanity level at the last scattering surface
the effective light attenuation as a function of redshift

Steps we will cover:

Make sure you really want to modify CAMB.
Make space for the module.
Copy existing CAMB code into the module.
Modify the CAMB files.
Add any new configuration options
Load any new inputs into the module.
Save any new outputs from the module.
Create an ini file to test the modified CAMB.
Run your CAMB!

1. Make sure you really want to modify CAMB.

If what you want to calculate can be derived just from the outputs of CAMB then it is best implemented as a new module rather than a modification to CAMB. For example, lots of observables are integrals over the matter power spectrum with some window function(s). Code to do this should be a new module since it can be decoupled from the actual calculation of that power spectrum.

If you want to implement new physics, though such as modified gravity or a new primordial power spectrum, then you probably need to modify CAMB.

2. Make space for the module.

CosmoSIS comes with a tool to create a code repository to contain new modules. We need a name for the project, and a name for the module(s) within the project:

#!bash 
    cosmosis-new-module space_mirror mod_camb

This will create a directory, modules/space_mirror/mod_camb, and put a Makefile in it. It will also initialize a git repository to store the code in, and show instructions on how to connect this repository to the BitBucket website so it will be stored safely and easy to share.

3. Copy existing CAMB code into the module.

If you look in the directory cosmosis-standard-library/boltzmann/camb you'll see these files:

#!bash 
    Makefile                 camb_interface.F90       camb_module.F90          patches
    camb_Jan15               module.yaml

We will need to copy all of these over to our new directory:

#!bash 
    cp -r cosmosis-standard-library/boltzmann/camb/* modules/space_mirror/mod_camb/

From now on you should work on this copied directory.

Other CAMB base versions

If you already have a modified camb based on a different version of CAMB that isn't in CosmoSIS (we currently have Nov13 as well as Jan15) then you will need to do a little more work. In the patches directory there is a collection of patch files that modify the vanilla CAMB files to make then work with CosmoSIS (it changes things like error handling and maximum number of redshift slices). In that case you will need to apply the patches in those files using the patch command:

#!bash 
    cd camb_Oct13
    patch < reionization_errorhandling.patch
    # ... and similar for the other patch files

4. Modify the CAMB files.

This is largely your problem! And it won't be fun - CAMB is painfully complicated.

If you already have your modifications in a separate CAMB source code then you will need to merge your changes into the CosmoSIS ones - this should not be very hard, and you can probably use the patch program as described above to apply the CosmoSIS changes to your files.

You should be modifying the files in the copy modules/space_mirror/mod_camb/Camb_Jan15 of course.

5. Add any new configuration options

We said that we wanted to add a new configuration option - let's call it nreflect_max and assume it is an integer. By configuration option we mean that it will be set at the start of the run and will not vary. We therefore want to read it from the parameter ini file.

Open the file camb_interface.F90, which provides functions and parameters connecting CAMB to CosmoSIS, and scroll down to the function camb_initial_setup:

#!fortran

    function camb_initial_setup(block, mode, fixed_mode) result(status)

This function is called just once, when the pipeline is first created. It reads any fixed parameters, like lmax, the feedback level, and what calculations should be done, from the main ini file.

The CosmoSIS camb interface can be called in four modes, depending on what calculations are needed:

background (fastest; background evolution only)
thermal (background + thermal history)
cmb (background + thermal history + CMB)
all (slowest; background + thermal history + CMB + matter power)

You may have parameters that only need to be read in certain modes, just as, for example, the lmax parameter is not needed in background or thermal modes. First we will add these parameters to the interface file so they can be recorded and set when needed.

We will assume that there is a global variable somewhere in your code that needs to be set - this may mean adding use statement at the top of this function if your new variable is in a module somewhere. If your new variable is not global and is instead in the CAMBParams CP structure then you should create new variables at the top of the interface file, set those here, and then set your CP variable from it when you read in the cosmological parameters.

Scroll down to the line status = 0. This is where we will read in the fixed parameters - let's assume that our parameter always needs to be read - if it is only needed for perturbation calculations we could encase it in the if (mode .ne. CAMB_MODE_BG) test:

#!fortran

    status = 0

    !Our new parameter
    status = status + datablock_get_int(block, option_section, "reflection_count_limit", reflection_count_limit)


    !End of the new part - this is the original from now on ...
    !We do not use the CMB lmax if only using the background mode
    if (mode .ne. CAMB_MODE_BG) then
        status = status + datablock_get_int_default(block, option_section, "lmax", default_lmax, standard_lmax)
        status = status + datablock_get_int_default(block, option_section, "k_eta_max_scalar", 2*standard_lmax, k_eta_max_scalar)
    endif

If our parameter was a double instead of a float, we would have used datablock_get_double, and if there had been a sensible default we would have used datablock_get_int_default and added a default value argument second from the end.

Our fixed parameters will now be read at the start of the chain, and if it is not found or is invalid (an integer value passed to a real variable, for example) then an error message will be produced.

6. Load any new inputs into the module.

You probably also have new, varied cosmological or other physical parameters in your new CAMB that you want to sample over in your MCMC. In this section we'll show how to get those parameters from the datablock that cosmosis will pass to your module.

Scroll down to the function:

#!fortran
    function camb_interface_set_params(block, params, background_only) result(status)

This function is called every time a new set of cosmological parameters is generated (i.e. at each step in the chain). It fills in params which is of type CambParams, and will later be used to set the camb global object called CP. The parameter background_only is True or False and we can read different parameters depending what is required. We could also read different values depending on what parameters were set in the setup function above.

Let's imagine that one of our new parameters is only required if perturbations are to be used, and the other is always needed. And further we'll see what to do if one of those parameters can have a default value where if it is not set then the default should be assumed. In that case we would modify the lines in this file like this:

#!fortran

    status = status + datablock_get_double(block, cosmo, "omega_nu", params%omegan)
    status = status + datablock_get_double(block, cosmo, "omega_k", params%omegak)
    status = status + datablock_get_double(block, cosmo, "hubble", params%H0)
    ! These lines above are the regular cosmological parameters from the Vanilla CAMB.
    !Here are is our first new parameter, thre radius of the space mirror:
    status = status + datablock_get_double(block, "space_mirror", "radius", params%mirror_radius)

    !This "if" statement is always here.  We'll add a new parameter to be used only if perturbations
    !are needed:
    if (perturbations) then

    !Our new parameter, the space mirror reflectivity.  Let's assume it has a default of 1.0.
    !We also need to put it in a named section, which in this case we just call "space_mirror".
    !This will be the name of the section ini the values.ini file later on.  There are a number
    !of pre-defined section names, like cosmological_parameters_section, matter_power_lin_section,
    !and plenty more, and higher up in the code we have set cosmo = cosmological_parameters_section
    status = status + datablock_get_double_default(block, "space_mirror", "reflectivity", 1.0D0, params%mirror_reflectivity)    

    ! And now back to the original code ...
    status = status + datablock_get_double(block, cosmo, "n_s",     params%initpower%an(1))
    status = status + datablock_get_double_default(block, cosmo, "k_s", default_pivot_scalar, params%initpower%k_0_scalar)

Now the code will read variables called reflectivity and radius from a datablock section called space_mirror. We'll see later how to make sure that the CosmoSIS samplers generate these values in this section in the first place.

7. Save any new outputs from the module.

Your new modified CAMB might just change how the CMB spectra, cosmic expansion, matter power, etc. are calculated. But it might also generate some new scalar or vector quantities that you want to save, either for debugging/plotting or for some later module to use. In this section we will look at how to save these quantities.

There are several functions in the file for saving different CAMB outputs - one for the transfer functions/matter power, one for sigma_8, one for the CMB, and one for distance measures. It's easiest to modify one of these functions rather than writing your own, and here we will show an example modifying the distance measures function.

Our two new pretend quantities to save are a scalar number called sanity and a vector attenuation. We assume that the former has been calculated already and saved in the CambParams CP global type, say an integer CP%sanity, and the latter has a function to calculate it, maybe called get_attenuation.

Scroll down to the function camb_interface_save_da. After the code saves the various distance measures and derived parameters, we will add our save code:

#!fortran
    status = status + datablock_put_double_array_1d(block, dist, "H", distance)
    status = status + datablock_put_metadata(block, dist, "H", "unit", "Mpc/c")
    !^^^ Original code

    !At this point in the code the vector "z" contains a series of values of z
    !that the distance measures D_A, D_M, etc. are sampled at.  We will use this 
    !to save our vector too, and will steal the vector "distance" to do it in.
    !We make sure first, though, that the contents of "distance" will not be needed.

    !First make the vector we want to save.  We could have made it earlier, of course.
    do i=1,nz
        distance(i) = get_attenuation(z(i))
    enddo
    !Now save both the redshift and attenuation in our own section
    status = status + datablock_put_double_array_1d(block, "space_mirror", "z", z)
    status = status + datablock_put_double_array_1d(block, "space_mirror", "attenuation", distance)

    !Now we've done that we also want to save our scalar value, the theory sanity
    status = status + datablock_put_int(block, "space_mirror", "sanity", CP%sanity)

We've now saved all the results we need back to the datablock. They could be used by future modules in the pipeline, or if we run the test sampler, saved to file.

8. Setting up an ini file to test the pipeline

Now all our code is ready we can compile and test it. Our module copied over the Makefile that CAMB used, so if you have any new Fortran code files to add then you must add them in the subdirectory Makefile. The upper directory Makefile just uses the camb library. Once you're ready you can compile with make, fix the inevitable errors, and then iterate until it compiles.

Once it is compiled we are ready to test with an ini file. We will base it off demo number 1, which runs the Vanilla CAMB code.

Copy the demo 1 files to work on them:

#!bash

    cp demos/demo1.ini modcamb.ini
    cp demos/values1.ini modcamb_values.ini

Now we make a few changes to our new modcamb.ini file.

Change where the resulting cosmological outputs are saved to:

#!ini

[test]
save_dir=modcamb_output_1

Change our pipeline to use the modified camb

And also get rid of Halofit for now, and change where cosmological parameters are read from:

#!ini
[pipeline]
modules = consistency mod_camb
values = modcamb_values.ini

We leave in the handy consistency module which translates between different ways of expressing the parameters (e.g. it knows that omega_m = omega_c+omega_b, etc.).

Tell CosmoSIS where to find our new modified camb

This section should be added anywhere in the ini file. It also contains our new configuration parameter, reflection_count_limit:

#!ini

[mod_camb]
file = modules/space_mirror/mod_camb/camb.so
mode=all
lmax=2500
feedback=2
reflection_count_limit = 476

The file entry points to the library that the Makefile will have made.

Set our new parameters in the values file.

At the bottom of modcamb_values.ini we would add a new section with our new parameters in:

#!ini
;  ... the original file...
w = -1.0
wa = 0.0

; and now our new part.  Because we chose a section name ``space_mirror`` when loading our parameters in the code,
; that's what we should do here:

[space_mirror]
; For now we just put fixed values for these parameters because we are just testing.
; Later we can set ranges to vary over.
radius = 1000000000.0
reflectivity = 0.95

9. Running!

We are now ready to run our modified camb:

#!bash

    cosmosis modcamb.ini

The resulting cosmological data will be saved in a directory modcamb_output_1, and you can make plots from it with the postprocess command:

#!bash

    postprocess modcamb.ini -o plots -p modcamb

Wiki

cosmosis / camb_mod_tutorial