Authors: Alkistis Pourtsidou, Dida Markovic
Using part of http://camb.readthedocs.io/en/latest/CAMBdemo.html
To run this Jupyter notebook you need to have CAMB and the CAMB python package installed. In order to install the CAMB python package on your computer follow the instructions in http://camb.readthedocs.io/en/latest/
In [ ]:
# Set up science stuff
from __future__ import division
import numpy as np
import scipy, scipy.interpolate
# Set up all kinds of libraries
import os
# Import custom modules in this folder
import misc
# Importing CAMB
import camb
In [ ]:
# Options for running this notebook
### Paths
run_name = '' # Choose a name if you want, set to '' if don't want
outpath = './outputs/' # Where to save outputs
TS = misc.get_timestamp() # Arbitrary time stamp so linked output files are obvious
code_name = 'GC-Fish-'
### Set up the precision of the interpolations
sig_figs = None # If want to round the P(k) and kh
n_k_bin = 800 # Settings for CAMB P(k)
minkh = 1e-4 # Settings for CAMB P(k)
maxkh = 10.0 # Settings for CAMB P(k)
interp_type = 'linear' # Type of interpolation
k_per_logint = 200 # Setting for CAMB points per k decade
### Fiducial cosmological parameters
hubble_fid = 0.67
omegab_fid = 0.022445/hubble_fid**2
omegac_fid = 0.121203/hubble_fid**2
om0_fid = omegac_fid + omegab_fid
H00_fid = 100*hubble_fid
Ass_fid = 2.1265e-9
nss_fid = 0.96
# Dark Energy parameters
w0_fid = -1.0
wa_fid = 0.0
gamma_fid = 0.545
# Speed of light
c = 3.0e5
### Set up the survey, Euclid Definition Study Report (RedBook)
### https://arxiv.org/abs/1110.3193
Area = 15000.0 #deg^2
omegatot = Area*pow(np.pi/180,2)
Dzbin = 0.1
sig_p = 0.001
### Set up the redshift binned functions
zlist = np.arange(0.7,2.1,Dzbin)
Nzbins = len(zlist)
biaslist = [1.083, 1.125, 1.104, 1.126, 1.208, 1.243, 1.282, 1.292, 1.363, 1.497, 1.486, \
1.491, 1.573, 1.568]
dn3 = [2434.28, 4364.812, 4728.559, 4825.798, 4728.797, 4507.625, 4269.851, 3720.657, 3104.309, \
2308.975, 1541.831, 1474.707, 893.716, 497.613]
### Set up the derivative steps
z_dep_step_default = 1e-5
shape_step = 1e-3
### Set up the Fisher Matrix integration (not power spectrum!)
kmin = 0.001
kmax = 0.2
### Set up the Fisher matrix calculation
params = {
0: 'lnH',
1: 'lnDA',
2: 'lnfsig8',
3: 'lnbsig8',
4: 'Ps',
5: 'ns',
6: 'wb',
7: 'wm',
8: 'h'}
shape_params = [5,6,7,8]
In [ ]:
# Geometry functions for a spatially flat Universe
# Functions for DE
def w_integrand(z,w0=w0_fid,wa=wa_fid):
return (1 + w0+wa*z/(z+1)) / (1+z)
def DE_evo(zc,w0=w0_fid,wa=wa_fid):
return np.exp(3*scipy.integrate.romberg(lambda z: w_integrand(z,w0,wa), 0, zc))
# Define E(z) = H(z)/H0
def Ez(zc,w0=w0_fid,wa=wa_fid):
return np.sqrt((1-om0_fid)*DE_evo(zc,w0,wa) + om0_fid*pow(1+zc,3))
def Hz(zc,w0=w0_fid,wa=wa_fid):
return Ez(zc,w0,wa)*H00_fid
# Define the cosmological distances
def drdz(zp,w0=w0_fid,wa=wa_fid):
return (c/H00_fid)/Ez(zp,w0,wa)
def rcom(zc,w0=w0_fid,wa=wa_fid):
return scipy.integrate.romberg(lambda z: drdz(z,w0,wa),0,zc)
def DA(zc,w0=w0_fid,wa=wa_fid):
return rcom(zc,w0,wa)/(1+zc)
In [ ]:
# LCDM growth rate and growth factor
def fg(zz,w0=w0_fid,wa=wa_fid,gamma=gamma_fid):
omz=om0_fid*pow(1+zz,3)/pow(Ez(zz,w0,wa),2)
return pow(omz,gamma)
def Dg_dz(zz,w0=w0_fid,wa=wa_fid):
return -fg(zz,w0,wa)/(1+zz)
def Dgz(zc,w0=w0_fid,wa=wa_fid):
start_z = 0.0
ans = scipy.integrate.romberg(lambda z: Dg_dz(z,w0,wa), start_z, zc)
return np.exp(ans)
In [ ]:
# Function to get the sigma_8 and P(k) at z = 0.0 from CAMB
def get_fiducial(h=hubble_fid,
wb=omegab_fid*pow(hubble_fid,2),
wm=(omegac_fid+omegab_fid)*pow(hubble_fid,2),
ns=nss_fid):
# Reset CAMB
pars = camb.CAMBparams()
pars.set_cosmology(H0=100.0*h, ombh2=wb, omch2=wm-wb,omk=0,mnu=0)
pars.set_dark_energy() #LCDM (default)
pars.InitPower.set_params(ns=ns, r=0, As=Ass_fid)
# Set linear matter power spectrum at z=0: P(k,z=0)
pars.set_matter_power(redshifts=[0.],kmax=maxkh/hubble_fid * 1.1,k_per_logint=k_per_logint)
pars.NonLinear = camb.model.NonLinear_none
# Calculate the intermediate results for these parameters
results = camb.get_results(pars)
results.calc_power_spectra(pars)
# Calculate the CAMB power spectra to be interpolated
kh, _z, pk = results.get_linear_matter_power_spectrum(have_power_spectra=True,nonlinear=False)
if sig_figs is not None:
kh = misc.roundsf(kh,sig_figs)
pk = misc.roundsf(pk,sig_figs)
# Only one element both times, because always only at z=0
return results.get_sigma8()[0], scipy.interpolate.interp1d(kh,pk[0],kind=interp_type)
In [ ]:
# The power spectra depend on the fiducial cosmology
# Construct the observed power spectrum, P_gg(k,μ,z)
# See e.g. Seo & Eisenstein 2003 (astro-ph/0307460), Equation (14)
def Pgg(kk,mu,zc=None,bgs8=None,fgs8=None):
if bgs8 is None: bgs8 = bgs8_fid[round(zc,2)]
if fgs8 is None: fgs8 = fgs8_fid[round(zc,2)]
# The ratio P(k)/σ8**2 is always fiducial, because fgs8 and bgs8 are taken as
# independent parameters, so they should not be changing the σ8.
return pow(bgs8+fgs8*mu**2,2) * (pk0_fid(kk)/sig8_fid**2)
In [ ]:
# Calculate the survey properties
# Redshift errors
def photoz(kk,mu,zc):
sigz = sig_p*(1+zc)
return np.exp(-pow(kk*mu,2)*pow(c*sigz,2)/pow(Hz(zc),2))
# Survey (bin) volume [Mpc^3]
def dVsurdz(zz):
return omegatot*c*pow(rcom(zz),2)/(H00_fid*Ez(zz))
def Vsur(zc):
return scipy.integrate.romberg(dVsurdz,zc-Dzbin/2,zc+Dzbin/2)
def Pshot(zc):
return Vsur(zc)/Ngal
# Effective volume (the hubble**3 factors are for units consistency)
def Veff(kk,mu,zc):
return hubble_fid**3*Vsur(zc)* \
(Pgg(kk,mu,zc)*photoz(kk,mu,zc)/(Pgg(kk,mu,zc)*photoz(kk,mu,zc) \
+ hubble_fid**3*Pshot(zc)))**2
In [ ]:
# Get the power spectrum quantities in the fiducial cosmology
sig8_fid, pk0_fid = get_fiducial()
# These are the fiducial values for the redshift-dependent (DE "model independent") parameters
fgs8_fid = {round(zc,2): fg(zc)*sig8_fid*Dgz(zc) for zc in zlist}
bgs8_fid = {round(zc,2): bg*sig8_fid*Dgz(zc) for bg,zc in zip(biaslist,zlist)}
In [ ]:
#See, for example, Wang et al 2010, https://arxiv.org/abs/1006.3517
# Dark Energy model-independent (aka redshift-dependent) parameters first
# Shot noise derivative analytically since it is trivial
def dlnP_dPs(kk,mu,zc):
return 1.0/(Pgg(kk,mu,zc)*photoz(kk,mu,zc))
# NUMERICAL DERIVATIVES
# Growth function parameter, f*σ8(z)
def dlnP_dlnfsig8(kk,mu,zc,step_fs8=z_dep_step_default):
fs8_fid = fgs8_fid[round(zc,2)]
fs8_p = fs8_fid*(1+step_fs8)
fs8_m = fs8_fid*(1-step_fs8)
Pgg_fid = Pgg(kk,mu,zc)
Pgg_p = Pgg(kk,mu,zc,fgs8=fs8_p)
Pgg_m = Pgg(kk,mu,zc,fgs8=fs8_m)
return (fs8_fid/Pgg_fid)*(Pgg_p-Pgg_m)/(2*step_fs8*fs8_fid)
# Linear galaxy bias parameter, b*σ8(z)
def dlnP_dlnbsig8(kk,mu,zc,step_bs8=z_dep_step_default):
bs8_fid = bgs8_fid[round(zc,2)]
bs8_p = bs8_fid*(1+step_bs8)
bs8_m = bs8_fid*(1-step_bs8)
Pgg_fid = Pgg(kk,mu,zc)
Pgg_p = Pgg(kk,mu,zc,bgs8=bs8_p)
Pgg_m = Pgg(kk,mu,zc,bgs8=bs8_m)
return (bs8_fid/Pgg_fid)*(Pgg_p-Pgg_m)/(2*step_bs8*bs8_fid)
# Angular diameter distance parameter, D_A(z)
def dlnP_dlnDA(kk,mu,zc,step_DA=z_dep_step_default):
DA_fid = DA(zc)
DA_p = DA_fid*(1+step_DA)
DA_m = DA_fid*(1-step_DA)
Pgg_fid = Pgg(kk,mu,zc)
k_p = kk*(DA_fid/DA_p)*np.sqrt(1+mu**2*((DA_p/DA_fid)**2 - 1))
mu_p = mu*(DA_p/DA_fid)/np.sqrt(1+mu**2*((DA_p/DA_fid)**2 - 1))
Pgg_p = (DA_fid/DA_p)**2 * Pgg(k_p,mu_p,zc)
k_m = kk*(DA_fid/DA_m)*np.sqrt(1+mu**2*((DA_m/DA_fid)**2 - 1))
mu_m = mu*(DA_m/DA_fid)/np.sqrt(1+mu**2*((DA_m/DA_fid)**2 - 1))
Pgg_m = (DA_fid/DA_m)**2 * Pgg(k_m,mu_m,zc)
return (DA_fid/Pgg_fid)*(Pgg_p-Pgg_m)/(2*step_DA*DA_fid)
# Hubble function parameter, H(z)
def dlnP_dlnH(kk,mu,zc,step_H=z_dep_step_default):
H_fid = Hz(zc)
H_p = H_fid*(1+step_H)
H_m = H_fid*(1-step_H)
Pgg_fid = Pgg(kk,mu,zc)
k_p = kk*np.sqrt(1+mu**2*((H_p/H_fid)**2 - 1))
mu_p = mu*(H_p/H_fid)/np.sqrt(1+mu**2*((H_p/H_fid)**2 - 1))
Pgg_p = (H_p/H_fid) * Pgg(k_p,mu_p,zc)
k_m = kk*np.sqrt(1+mu**2*((H_m/H_fid)**2 - 1))
mu_m = mu*(H_m/H_fid)/np.sqrt(1+mu**2*((H_m/H_fid)**2 - 1))
Pgg_m = (H_m/H_fid) * Pgg(k_m,mu_m,zc)
return (H_fid/Pgg_fid)*(Pgg_p-Pgg_m)/(2*step_H*H_fid)
In [ ]:
# Parameters that change the broad-band shape of the power spectrum as well
# aka "shape parameters", which are redshift-independent
# Spectral index, n_s
step_ns = shape_step*nss_fid
ns_p = nss_fid + step_ns # DERIV PLUS
ns_m = nss_fid - step_ns # DERIV MINUS
sig8_ns_p, pk0_ns_p = get_fiducial(ns=ns_p);
sig8_ns_m, pk0_ns_m = get_fiducial(ns=ns_m);
def dlnP_dns(kk,zc):
return (pk0_ns_p(kk)/sig8_ns_p**2-pk0_ns_m(kk)/sig8_ns_m**2)/(2*step_ns)/(pk0_fid(kk)/sig8_fid**2)
# Baryon density, omega_b
wb_fid = omegab_fid*pow(hubble_fid,2)
step_wb = shape_step*wb_fid
wb_p = wb_fid + step_wb # DERIV PLUS
wb_m = wb_fid - step_wb # DERIV MINUS
sig8_wb_p, pk0_wb_p = get_fiducial(wb=wb_p);
sig8_wb_m, pk0_wb_m = get_fiducial(wb=wb_m);
def dlnP_dwb(kk,zc):
return (pk0_wb_p(kk)/sig8_wb_p**2-pk0_wb_m(kk)/sig8_wb_m**2)/(2*step_wb)/(pk0_fid(kk)/sig8_fid**2)
# Total matter density, omega_m
wm_fid = (omegab_fid + omegac_fid)*pow(hubble_fid,2)
step_wm = shape_step*wm_fid
wm_p = wm_fid + step_wm # DERIV PLUS
wm_m = wm_fid - step_wm # DERIV MINUS
sig8_wm_p, pk0_wm_p = get_fiducial(wm=wm_p);
sig8_wm_m, pk0_wm_m = get_fiducial(wm=wm_m);
def dlnP_dwm(kk,zc):
return (pk0_wm_p(kk)/sig8_wm_p**2-pk0_wm_m(kk)/sig8_wm_m**2)/(2*step_wm)/(pk0_fid(kk)/sig8_fid**2)
# Hubble parameter, h
step_h = shape_step*hubble_fid
h_p = hubble_fid + step_h # DERIV PLUS
h_m = hubble_fid - step_h # DERIV MINUS
sig8_h_p, pk0_h_p = get_fiducial(h=h_p);
sig8_h_m, pk0_h_m = get_fiducial(h=h_m);
def dlnP_dh(kk,zc):
return (pk0_h_p(kk)/sig8_h_p**2-pk0_h_m(kk)/sig8_h_m**2)/(2*step_h)/(pk0_fid(kk)/sig8_fid**2)
In [ ]:
# Auxiliary bits needed for Fisher matrix calculation
def dF(kk,mu,zc):
return (1./(8*np.pi*np.pi))*pow(kk,2)*deriv_i(kk,mu,zc)*deriv_j(kk,mu,zc)*Veff(kk,mu,zc)
# 2D integration function
def integrate2D(dfun, kgrid, mugrid):
muint = [scipy.integrate.romb(dfun.T[i], dx=np.diff(mugrid)[0]) for i in range(kgrid.size)]
return scipy.integrate.romb(muint, dx=np.diff(kgrid)[0])
In [ ]:
%%time
# Create array of zeros
Npar = len(params)
Npar_shape = len(shape_params)
s = (Npar-Npar_shape)*Nzbins + Npar_shape
Fishermat = np.zeros([s,s])
# Loop over redshift, k and mu and get all the Fisher matrix entries by calling the derivatives
kgrid = np.linspace(kmin, kmax, 513)
mugrid = np.linspace(-1., 1., 513)
K, MU = np.meshgrid(kgrid, mugrid)
for zi in range(0,Nzbins):
zc = zlist[zi]
bg = biaslist[zi]
zmin = zc-Dzbin/2
zmax = zc+Dzbin/2
Ngal = dn3[zi]*Area*Dzbin
for i in range(0,Npar):
if i not in shape_params:
k=zi*(Npar-Npar_shape) + i
else:
k = s + (i - Npar)
def deriv_i(kk,mu,z):
if i==0: return dlnP_dlnH(kk,mu,z)
elif i==1: return dlnP_dlnDA(kk,mu,z)
elif i==2: return dlnP_dlnfsig8(kk,mu,z)
elif i==3: return dlnP_dlnbsig8(kk,mu,z)
elif i==4: return dlnP_dPs(kk,mu,z)
elif i==5: return dlnP_dns(kk,z)
elif i==6: return dlnP_dwb(kk,z)
elif i==7: return dlnP_dwm(kk,z)
elif i==8: return dlnP_dh(kk,z)
else: print "Error: index out of range"
for j in range(0,Npar):
if j not in shape_params:
l=zi*(Npar-Npar_shape) + j
else:
l = s + (j - Npar)
if j>=i:
def deriv_j(kk,mu,z):
if j==0: return dlnP_dlnH(kk,mu,z)
elif j==1: return dlnP_dlnDA(kk,mu,z)
elif j==2: return dlnP_dlnfsig8(kk,mu,z)
elif j==3: return dlnP_dlnbsig8(kk,mu,z)
elif j==4: return dlnP_dPs(kk,mu,z)
elif j==5: return dlnP_dns(kk,z)
elif j==6: return dlnP_dwb(kk,z)
elif j==7: return dlnP_dwm(kk,z)
elif j==8: return dlnP_dh(kk,z)
else: print "Error: index out of range"
if l>=k:
Fishermat[k][l] += integrate2D(dF(K,MU,zc),kgrid,mugrid)
else:
Fishermat[k,l] = Fishermat[l,k]
In [ ]:
# Save the resulting Fisher matrix into a file
header = 'lnH_0.7 lnDa_0.7 lnfs8_0.7 lnbs8_0.7 Ps_0.7 ' +\
'lnH_0.8 lnDa_0.8 lnfs8_0.8 lnbs8_0.8 Ps_0.8 ' +\
'lnH_0.9 lnDa_0.9 lnfs8_0.9 lnbs8_0.9 Ps_0.9 ' +\
'lnH_1.0 lnDa_1.0 lnfs8_1.0 lnbs8_1.0 Ps_1.0 ' +\
'lnH_1.1 lnDa_1.1 lnfs8_1.1 lnbs8_1.1 Ps_1.1 ' +\
'lnH_1.2 lnDa_1.2 lnfs8_1.2 lnbs8_1.2 Ps_1.2 ' +\
'lnH_1.3 lnDa_1.3 lnfs8_1.3 lnbs8_1.3 Ps_1.3 ' +\
'lnH_1.4 lnDa_1.4 lnfs8_1.4 lnbs8_1.4 Ps_1.4 ' +\
'lnH_1.5 lnDa_1.5 lnfs8_1.5 lnbs8_1.5 Ps_1.5 ' +\
'lnH_1.6 lnDa_1.6 lnfs8_1.6 lnbs8_1.6 Ps_1.6 ' +\
'lnH_1.7 lnDa_1.7 lnfs8_1.7 lnbs8_1.7 Ps_1.7 ' +\
'lnH_1.8 lnDa_1.8 lnfs8_1.8 lnbs8_1.8 Ps_1.8 ' +\
'lnH_1.9 lnDa_1.9 lnfs8_1.9 lnbs8_1.9 Ps_1.9 ' +\
'lnH_2.0 lnDa_2.0 lnfs8_2.0 lnbs8_2.0 Ps_2.0 ' +\
'ns wb wm h'
# Generate an output file name
newfile = os.path.abspath(os.path.join(outpath, code_name + '_' + interp_type))
# Save file
git = misc.GitEnv()
np.savetxt(newfile+TS+'_'+git.get_hash(7)+'_'+run_name+'.txt', Fishermat, header=header)
# Print for testing purposes
print np.sqrt(np.diag(scipy.linalg.inv(Fishermat)))
Alkistis & Dida, 2017