CAMB Python example notebook

Run it online yourself in Binder.

In [1]:
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import sys, platform, os
import matplotlib
from matplotlib import pyplot as plt
import numpy as np

#Assume installed from github using "git clone --recursive https://github.com/cmbant/CAMB.git"
#This file is then in the docs folders
camb_path = os.path.realpath(os.path.join(os.getcwd(),'..'))
sys.path.insert(0,camb_path)
import camb
from camb import model, initialpower
print('Using CAMB %s installed at %s'%(camb.__version__,os.path.dirname(camb.__file__)))
Using CAMB 1.0.3 installed at c:\work\dist\git\camb\camb
In [2]:
#Set up a new set of parameters for CAMB
pars = camb.CAMBparams()
#This function sets up CosmoMC-like settings, with one massive neutrino and helium set using BBN consistency
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122, mnu=0.06, omk=0, tau=0.06)
pars.InitPower.set_params(As=2e-9, ns=0.965, r=0)
pars.set_for_lmax(2500, lens_potential_accuracy=0);
In [3]:
#calculate results for these parameters
results = camb.get_results(pars)
In [4]:
#get dictionary of CAMB power spectra
powers =results.get_cmb_power_spectra(pars, CMB_unit='muK')
for name in powers: print(name)
total
unlensed_scalar
unlensed_total
lensed_scalar
tensor
lens_potential
In [5]:
#plot the total lensed CMB power spectra versus unlensed, and fractional difference
totCL=powers['total']
unlensedCL=powers['unlensed_scalar']
print(totCL.shape)
#Python CL arrays are all zero based (starting at L=0), Note L=0,1 entries will be zero by default.
#The different CL are always in the order TT, EE, BB, TE (with BB=0 for unlensed scalar results).
ls = np.arange(totCL.shape[0])
fig, ax = plt.subplots(2,2, figsize = (12,12))
ax[0,0].plot(ls,totCL[:,0], color='k')
ax[0,0].plot(ls,unlensedCL[:,0], color='r')
ax[0,0].set_title('TT')
ax[0,1].plot(ls[2:], 1-unlensedCL[2:,0]/totCL[2:,0]);
ax[0,1].set_title(r'$\Delta TT$')
ax[1,0].plot(ls,totCL[:,1], color='k')
ax[1,0].plot(ls,unlensedCL[:,1], color='r')
ax[1,0].set_title(r'$EE$')
ax[1,1].plot(ls,totCL[:,3], color='k')
ax[1,1].plot(ls,unlensedCL[:,3], color='r')
ax[1,1].set_title(r'$TE$');
for ax in ax.reshape(-1): ax.set_xlim([2,2500]);
(2551, 4)
In [6]:
# The lensing B modes are non-linear, so need to be calculated carefully if you want them accurate (even at low ell)
# Need both high lmax, non-linear lensing and high k 
# lens_potential_accuracy=1 turns on the latter, and can be increased to check precision

pars.set_for_lmax(2500, lens_potential_accuracy=1)
results = camb.get_results(pars)
lmax2500CL = results.get_lensed_scalar_cls(CMB_unit='muK')
                  
pars.set_for_lmax(4000, lens_potential_accuracy=1)
results = camb.get_results(pars)
lmax4000CL = results.get_lensed_scalar_cls(CMB_unit='muK')

pars.set_for_lmax(4000, lens_potential_accuracy=2)
results = camb.get_results(pars)
accCL = results.get_lensed_scalar_cls(CMB_unit='muK')

pars.set_for_lmax(6000, lens_potential_accuracy=4)
results = camb.get_results(pars)
refCL = results.get_lensed_scalar_cls(CMB_unit='muK')

fig, ax = plt.subplots(1,2, figsize = (12,4))
ax[0].plot(ls,totCL[:len(ls),2], color='C0')
ax[0].plot(ls,lmax2500CL[:len(ls),2], color='C1')
ax[0].plot(ls,lmax4000CL[:len(ls),2], color='C2')
ax[0].plot(ls,accCL[:len(ls),2], color='C3')
ax[0].plot(ls,refCL[:len(ls),2], color='k')

ax[0].set_xlim([2,2500])
ax[0].set_xlabel(r'$\ell$',fontsize=13)
ax[0].set_ylabel(r'$\ell(\ell+1)C_\ell^{BB}/2\pi\,[\mu {\rm K}^2]$', fontsize=13)

ax[1].plot(ls[2:],totCL[2:len(ls),2]/refCL[2:len(ls),2]-1, color='C0')
ax[1].plot(ls[2:],lmax2500CL[2:len(ls),2]/refCL[2:len(ls),2]-1, color='C1')
ax[1].plot(ls[2:],lmax4000CL[2:len(ls),2]/refCL[2:len(ls),2]-1, color='C2')
ax[1].plot(ls[2:],accCL[2:len(ls),2]/refCL[2:len(ls),2]-1, color='C3')

ax[1].axhline(0,color='k')
ax[1].set_xlim([2,2500])
ax[1].set_xlabel(r'$\ell$',fontsize=13)
ax[1].set_ylabel('fractional error', fontsize=13);
ax[1].legend(['Default accuracy','lmax=2500, lens_potential_accuracy=1',
           'lmax=4000, lens_potential_accuracy=1','lmax=4000, lens_potential_accuracy=2']);
In [7]:
#You can calculate spectra for different primordial power spectra without recalculating everything
#for example, let's plot the BB spectra as a function of r
pars.set_for_lmax(4000, lens_potential_accuracy=1)
pars.WantTensors = True
results = camb.get_transfer_functions(pars)
lmax=2000
rs = np.linspace(0,0.2,6)
for r in rs:
    inflation_params = initialpower.InitialPowerLaw()
    inflation_params.set_params(ns=0.96, r=r)
    results.power_spectra_from_transfer(inflation_params) #warning OK here, not changing scalars
    cl = results.get_total_cls(lmax, CMB_unit='muK')
    plt.loglog(np.arange(lmax+1),cl[:,2])
plt.xlim([2,lmax])
plt.legend(["$r = %s$"%r for r in  rs], loc='lower right');
plt.ylabel(r'$\ell(\ell+1)C_\ell^{BB}/ (2\pi \mu{\rm K}^2)$')
plt.xlabel(r'$\ell$');
WARNING:root:power_spectra_from_transfer with non-linear lensing does not recalculate the non-linear correction
In [8]:
#Now get matter power spectra and sigma8 at redshift 0 and 0.8
pars = camb.CAMBparams()
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122)
pars.InitPower.set_params(ns=0.965)
#Note non-linear corrections couples to smaller scales than you want
pars.set_matter_power(redshifts=[0., 0.8], kmax=2.0)

#Linear spectra
pars.NonLinear = model.NonLinear_none
results = camb.get_results(pars)
kh, z, pk = results.get_matter_power_spectrum(minkh=1e-4, maxkh=1, npoints = 200)
s8 = np.array(results.get_sigma8())

#Non-Linear spectra (Halofit)
pars.NonLinear = model.NonLinear_both
results.calc_power_spectra(pars)
kh_nonlin, z_nonlin, pk_nonlin = results.get_matter_power_spectrum(minkh=1e-4, maxkh=1, npoints = 200)
Note: redshifts have been re-sorted (earliest first)
In [9]:
print(results.get_sigma8())
[0.53298468 0.80251083]
In [10]:
for i, (redshift, line) in enumerate(zip(z,['-','--'])):
    plt.loglog(kh, pk[i,:], color='k', ls = line)
    plt.loglog(kh_nonlin, pk_nonlin[i,:], color='r', ls = line)
plt.xlabel('k/h Mpc');
plt.legend(['linear','non-linear'], loc='lower left');
plt.title('Matter power at z=%s and z= %s'%tuple(z));
In [11]:
#Plot CMB lensing potential power for various values of w
pars = camb.CAMBparams()
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122)
pars.InitPower.set_params(As=2e-9, ns=0.965)
pars.set_for_lmax(2000, lens_potential_accuracy=1)

ws = np.linspace(-1.5, -0.6, 5)
for w in ws:
    pars.set_dark_energy(w=w, wa=0, dark_energy_model='fluid') 
    results = camb.get_results(pars)
    cl = results.get_lens_potential_cls(lmax=2000)
    plt.loglog(np.arange(2001), cl[:,0])

plt.legend(ws)
plt.ylabel('$[L(L+1)]^2C_L^{\phi\phi}/2\pi$')
plt.xlabel('$L$')
plt.xlim([2,2000]);

You can view the parameters (as used by fortran CAMB internals) by just printing the parameter object. See the docs for parameter and structure descriptions

In [12]:
# parameters can also be read from text .ini files, for example this sets up a best-fit 
# Planck 2018 LCDM cosmology (base_plikHM_TTTEEE_lowl_lowE_lensing). 
# [Use planck_2018_acc.ini if you need high-ell and/or accurate BB and CMB lensng spectra at beyond-Planck accuracy]

pars=camb.read_ini(os.path.join(camb_path,'inifiles','planck_2018.ini'))
print(pars)
class: <CAMBparams>
 WantCls = True
 WantTransfer = True
 WantScalars = True
 WantTensors = False
 WantVectors = False
 WantDerivedParameters = True
 Want_cl_2D_array = True
 Want_CMB = True
 Want_CMB_lensing = True
 DoLensing = True
 NonLinear = NonLinear_both
 Transfer: <TransferParams>
   high_precision = True
   accurate_massive_neutrinos = False
   kmax = 1.3464234
   k_per_logint = 0
   PK_num_redshifts = 1
   PK_redshifts = [0.0]
 want_zstar = False
 want_zdrag = False
 min_l = 2
 max_l = 2700
 max_l_tensor = 600
 max_eta_k = 18000.0
 max_eta_k_tensor = 1200.0
 ombh2 = 0.0223828
 omch2 = 0.1201075
 omk = 0.0
 omnuh2 = 0.0006451439
 H0 = 67.32117
 TCMB = 2.7255
 YHe = 0.2454006
 num_nu_massless = 2.046
 num_nu_massive = 1
 nu_mass_eigenstates = 1
 share_delta_neff = True
 nu_mass_degeneracies = [0.0]
 nu_mass_fractions = [1.0]
 nu_mass_numbers = [1]
 InitPower: <InitialPowerLaw>
   tensor_parameterization = tensor_param_indeptilt
   ns = 0.9660499
   nrun = 0.0
   nrunrun = 0.0
   nt = 0.0
   ntrun = 0.0
   r = 0.0
   pivot_scalar = 0.05
   pivot_tensor = 0.05
   As = 2.100549e-09
   At = 0.0
 Recomb: <Recfast>
   min_a_evolve_Tm = 0.0011098779505118728
   RECFAST_fudge = 1.125
   RECFAST_fudge_He = 0.86
   RECFAST_Heswitch = 6
   RECFAST_Hswitch = True
   AGauss1 = -0.14
   AGauss2 = 0.079
   zGauss1 = 7.28
   zGauss2 = 6.73
   wGauss1 = 0.18
   wGauss2 = 0.33
 Reion: <TanhReionization>
   Reionization = True
   use_optical_depth = True
   redshift = 10.0
   optical_depth = 0.05430842
   delta_redshift = 0.5
   fraction = -1.0
   include_helium_fullreion = True
   helium_redshift = 3.5
   helium_delta_redshift = 0.5
   helium_redshiftstart = 6.0
   tau_solve_accuracy_boost = 1.0
   timestep_boost = 1.0
   max_redshift = 50.0
 DarkEnergy: <DarkEnergyPPF>
   w = -1.0
   wa = 0.0
   cs2 = 1.0
   use_tabulated_w = False
 NonLinearModel: <Halofit>
   Min_kh_nonlinear = 0.005
   halofit_version = mead
 Accuracy: <AccuracyParams>
   AccuracyBoost = 1.0
   lSampleBoost = 1.0
   lAccuracyBoost = 1.0
   AccuratePolarization = True
   AccurateBB = False
   AccurateReionization = True
   TimeStepBoost = 1.0
   BackgroundTimeStepBoost = 1.0
   IntTolBoost = 1.0
   SourcekAccuracyBoost = 1.0
   IntkAccuracyBoost = 1.0
   TransferkBoost = 1.0
   NonFlatIntAccuracyBoost = 1.0
   BessIntBoost = 1.0
   LensingBoost = 1.0
   NonlinSourceBoost = 1.0
   BesselBoost = 1.0
   LimberBoost = 1.0
   SourceLimberBoost = 1.0
   KmaxBoost = 1.0
   neutrino_q_boost = 1.0
 SourceTerms: <SourceTermParams>
   limber_windows = True
   limber_phi_lmin = 100
   counts_density = True
   counts_redshift = True
   counts_lensing = False
   counts_velocity = True
   counts_radial = False
   counts_timedelay = True
   counts_ISW = True
   counts_potential = True
   counts_evolve = False
   line_phot_dipole = False
   line_phot_quadrupole = False
   line_basic = True
   line_distortions = True
   line_extra = False
   line_reionization = False
   use_21cm_mK = True
 z_outputs = []
 scalar_initial_condition = initial_adiabatic
 InitialConditionVector = []
 OutputNormalization = 1
 Alens = 1.0
 MassiveNuMethod = Nu_trunc
 DoLateRadTruncation = True
 Evolve_baryon_cs = False
 Evolve_delta_xe = False
 Evolve_delta_Ts = False
 Do21cm = False
 transfer_21cm_cl = False
 Log_lvalues = False
 use_cl_spline_template = True
 SourceWindows = []
 CustomSources: <CustomSources>
   num_custom_sources = 0
   c_source_func = None
   custom_source_ell_scales = []
 
In [13]:
#The dark energy model can be changed as in the previous example, or by assigning to pars.DarkEnergy.
#e.g. use the PPF model
from camb.dark_energy import DarkEnergyPPF, DarkEnergyFluid
pars.DarkEnergy = DarkEnergyPPF(w=-1.2, wa=0.2)
print('w, wa model parameters:\n\n', pars.DarkEnergy)
results = camb.get_background(pars)

#or can also use a w(a) numerical function 
#(note this will slow things down; make your own dark energy class in fortran for best performance)
a = np.logspace(-5, 0, 1000)
w = -1.2 + 0.2 * (1 - a)
pars.DarkEnergy= DarkEnergyPPF()
pars.DarkEnergy.set_w_a_table(a, w)
print('Table-interpolated parameters (w and wa are set to estimated values at 0):\n\n' 
      ,pars.DarkEnergy)
results2 = camb.get_background(pars)

rho, _ = results.get_dark_energy_rho_w(a)
rho2, _ = results2.get_dark_energy_rho_w(a)
plt.plot(a, rho, color='k')
plt.plot(a, rho2, color='r', ls='--')
plt.ylabel(r'$\rho/\rho_0$')
plt.xlabel('$a$')
plt.xlim(0,1)
plt.title('Dark enery density');
w, wa model parameters:

 class: <DarkEnergyPPF>
 w = -1.2
 wa = 0.2
 cs2 = 1.0
 use_tabulated_w = False
 
Table-interpolated parameters (w and wa are set to estimated values at 0):

 class: <DarkEnergyPPF>
 w = -1.2
 wa = 0.1993346429183863
 cs2 = 1.0
 use_tabulated_w = True
 
In [14]:
#Get various background functions and derived parameters
pars = camb.CAMBparams()
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122)
results = camb.get_background(pars)
print('Derived parameter dictionary: %s'%results.get_derived_params())
Derived parameter dictionary: {'age': 13.74073863845039, 'zstar': 1090.5659735246988, 'rstar': 144.20177032101333, 'thetastar': 1.0445808915883217, 'DAstar': 13.804749012950973, 'zdrag': 1059.2190440609993, 'rdrag': 146.98003046025178, 'kd': 0.1406940577023048, 'thetad': 0.16175048908224343, 'zeq': 3441.168556293452, 'keq': 0.010502750486921146, 'thetaeq': 0.8080627625764873, 'thetarseq': 0.4470154731275294}
In [15]:
z = np.linspace(0,4,100)
DA = results.angular_diameter_distance(z)
plt.plot(z, DA)
plt.xlabel('$z$')
plt.ylabel(r'$D_A /\rm{Mpc}$')
plt.title('Angular diameter distance')
plt.ylim([0,2000])
plt.xlim([0,4]);
In [16]:
print('CosmoMC theta_MC parameter: %s'%results.cosmomc_theta())
CosmoMC theta_MC parameter: 0.010443417642533299
In [17]:
#You can also directly access some lower level quantities, for example the CMB transfer functions:
pars = camb.CAMBparams()
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122)
data = camb.get_transfer_functions(pars)
transfer = data.get_cmb_transfer_data()
print('Number of sources (T, E, phi..): %s; number of ell: %s; number of k: %s '%tuple(transfer.delta_p_l_k.shape))
Number of sources (T, E, phi..): 3; number of ell: 87; number of k: 2737 
In [18]:
#Plot the transfer functions as a function of k for various ell
fig, axs = plt.subplots(2,2, figsize=(12,8), sharex = True)
for ix, ax in zip([3, 20, 40, 60],axs.reshape(-1)):
    ax.plot(transfer.q,transfer.delta_p_l_k[0,ix,:])
    ax.set_title(r'$\ell = %s$'%transfer.l[ix])
    if ix>1: ax.set_xlabel(r'$k \rm{Mpc}$')