Demonstrations of RESPRESSO package¶
In this notebook, I present how to use the module to compute the nonlinear power spectrum.
%pylab inline
import respresso
Some parameters for plotting:
plt.rcParams['font.family'] = 'sans-serif'
plt.rcParams['font.size'] = 18
plt.rcParams['axes.linewidth'] = 1.5
plt.rcParams['xtick.major.size'] = 5
plt.rcParams['ytick.major.size'] = 5
plt.rcParams['xtick.minor.size'] = 3
plt.rcParams['ytick.minor.size'] = 3
plt.rcParams['xtick.minor.visible'] = True
plt.rcParams['ytick.minor.visible'] = True
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'
Below is the main RESPRESSO object for the reconstruction. All the relevant precomputed data are loaded when instantiated.
respresso_obj = respresso.respresso_core()
Demonstration 1: from tabulated linear power spectrum to its nonlinear counterpart¶
In the first example, I start from a linear power spectrum of your target cosmology normalized at the target redshift prepared as a simple 2-columns data table. After the table is read, I create its spline interpolater and give it to the RESPRESSO object. The reconstruction is done after determining the path for the reconstruction. I use InterpolatedUniveriateSpline module in scipy to interpolate the linear power spectrum table.
from scipy.interpolate import InterpolatedUnivariateSpline as ius
First, prepare a linear power spectrum interpolater for the target cosmology from a table. The format is $k,[h\,\mathrm{Mpc}^{-1}]$ vs $P(k)\,[h^{-3}\mathrm{Mpc}^3]$. Our $2\pi$ convention in Fourier Transofrm is: $$\delta(\mathbf{x}) = \int\frac{\mathrm{d}^3k}{(2\pi)^3}e^{i\mathbf{k}\cdot\mathbf{x}} \delta_\mathbf{k}$$
plin_target = np.loadtxt('test_data/linear_power/plin_wmap5_z1.dat').T
plin_target_spl = ius(plin_target[0,],plin_target[1,])
By the following command, I am telling the RESPRESSO object the target power spectrum.
respresso_obj.set_target(plin_target_spl)
Now, let the object find the "path" in the linear power spectrum space based on the distance between the two models. It prints the number of intermediate steps. You can force it to take a given number of non-negative number of intermediate steps by giving it an argument, if you want.
respresso_obj.find_path()
Check the estimated maximum wavenumber below which the accuracy is within $\sim1\%$. This does not always guarantee you $1\%$ especially when the target cosmology is very far away from the fiducial. The value is internally computed and stored in the "find_path" process above.
kmax = respresso_obj.get_kmax()
print 'Rough estimate of the maximum wavenumber: k_max =', kmax, ' h/Mpc'
Everything is ready now. The "reconstruct" routine gives you the reconstructed nonlinear power spectrum at internal wavenumber samples $k_\mathrm{sample}\,[h\,\mathrm{Mpc}^{-1}]$. You can retrieve $k_\mathrm{sample}$ by the "get_kinternal" routine. The internal wavenumbers are always the same, so you can reuse the "kwave" array in what follows.
kwave = respresso_obj.get_kinternal()
pnl_rec = respresso_obj.reconstruct()
I show below the linear and nonlinear power spectra for the target cosmology. The vertical allow marks the location of $k_\mathrm{max}$.
plt.figure(figsize=(10,6))
plt.loglog(kwave,pnl_rec,label='nonlinear')
plt.loglog(kwave,plin_target_spl(kwave),':',label='linear')
plt.annotate("", xy=(kmax, 50), xytext=(kmax, 20),arrowprops=dict(arrowstyle="->",lw=5))
plt.xlabel('$k\,[h\,\mathrm{Mpc}^{-1}]$')
plt.ylabel('$P(k)\,[h^{-3}\mathrm{Mpc}^{3}]$')
plt.legend()
Demonstration 2: from tabulated linear power spectrum at z=0 to its nonlinear counterpart¶
Very often, one has a data table of the linear power spectrum normalized at z=0 and would like to compute the nonlinear power spectrum at a different redshift. RESPRESSO internally has the relevant routine to rescale the linear power spectrum following the linear growth factor. The function, "linearGrowth", simply takes the scale factor, Omega_m0, and dark energy w parameter and gives the linear growth factor (for flat wCDM cosmology). Let's try z=0.5 this time for the same WMAP5 cosmology.
z = 0.5
ascale = 1/(1+z)
Om, w = 0.279, -1
D2 = respresso.linearGrowth(ascale,Om,w)**2
plin_target = np.loadtxt('test_data/linear_power/plin_wmap5_z0.dat').T
plin_target_spl = ius(plin_target[0,],D2*plin_target[1,])
Then, the rest is exactly the same as before:
respresso_obj.set_target(plin_target_spl)
respresso_obj.find_path()
kmax = respresso_obj.get_kmax()
print 'Rough estimate of the maximum wavenumber: k_max =', kmax, ' h/Mpc'
kwave = respresso_obj.get_kinternal()
pnl_rec = respresso_obj.reconstruct()
plt.figure(figsize=(10,6))
plt.loglog(kwave,pnl_rec,label='nonlinear')
plt.loglog(kwave,plin_target_spl(kwave),':',label='linear')
plt.annotate("", xy=(kmax, 70), xytext=(kmax, 30),arrowprops=dict(arrowstyle="->",lw=5))
plt.xlabel('$k\,[h\,\mathrm{Mpc}^{-1}]$')
plt.ylabel('$P(k)\,[h^{-3}\mathrm{Mpc}^{3}]$')
plt.legend()
Demonstration 3: from CAMB transfer function to nonlinear power spectrum¶
Let's try the same thing from a table of T(k) obtained by the CAMB routine (you have to arrange the format to $k\,[h\mathrm{Mpc}^{-1}]$ vs $T_{matter}(k)$ before passing it). I provide an auxiliary class called "Cosmology", which accepts cosmological parameters and a linear matter transfer function. It computes linear matter power spectrum at a desired redshift, which you can pass to the main RESPRESSO routine to do the same as in the previous demonstration.
c1 = respresso.Cosmology(0.234,0.734,2.37301e-09,0.961,0.002,-1.,'test_data/transfer/tkwmap3.dat') # Om0, hubble, As, ns, k0[1/Mpc], w and transfer table
Then compute the linear power spectrum and create its spline interpolator:
z = 0.5
ks = np.logspace(-3,1,400)
plin_data = c1.get_plin(ks,z)
plin_spl = ius(ks,plin_data)
The remaining procedures are exactly the same as previously:
respresso_obj.set_target(plin_spl)
respresso_obj.find_path()
kmax = respresso_obj.get_kmax()
print 'Rough estimate of the maximum wavenumber: k_max =', kmax, ' h/Mpc'
kwave = respresso_obj.get_kinternal()
pnl_rec = respresso_obj.reconstruct()
plt.figure(figsize=(10,6))
plt.loglog(kwave,pnl_rec,label='nonlinear')
plt.loglog(kwave,plin_spl(kwave),':',label='linear')
plt.annotate("", xy=(kmax, 60), xytext=(kmax, 25),arrowprops=dict(arrowstyle="->",lw=5))
plt.xlabel('$k\,[h\,\mathrm{Mpc}^{-1}]$')
plt.ylabel('$P(k)\,[h^{-3}\mathrm{Mpc}^{3}]$')
plt.legend()
Demonstration 4: from cosmological paramters to nonlinear power spectrum¶
Let's try the same thing from your cosmological paramters. Below I use the public Boltzman code CLASS to compute the linear power spectrum. Please refer to their webpage http://class-code.net/ and/or https://github.com/lesgourg/class_public/wiki/Python-wrapper to learn how to install its python wrapper "classy".
from classy import Class
Prepare a Class instance and tell it a list of cosmological parameters and then execute it by "compute":
cosmo = Class()
params = {'output': 'mPk','T_cmb': '2.7255','h': '0.7','Omega_b': '0.03','Omega_cdm': '0.12','tau_reio': '0.079','n_s': '0.96','ln10^{10}A_s': '3.','Omega_k': '0.0','k_pivot': '0.05','P_k_max_h/Mpc': '10.0','z_max_pk': '5'}
cosmo.set(params)
cosmo.compute()
Compute the linear power spectrum at the desired redshift:
z = 1.
ks = np.logspace(-3,1,400)
plin_data = [cosmo.pk(ks[i]*cosmo.h(), z)*cosmo.h()**3 for i in range(400)]
plin_spl = ius(ks,plin_data)
The below are exactly the same as before again.
respresso_obj.set_target(plin_spl)
respresso_obj.find_path()
kmax = respresso_obj.get_kmax()
print 'Rough estimate of the maximum wavenumber: k_max =', kmax, ' h/Mpc'
kwave = respresso_obj.get_kinternal()
pnl_rec = respresso_obj.reconstruct()
plt.figure(figsize=(10,6))
plt.loglog(kwave,pnl_rec,label='nonlinear')
plt.loglog(kwave,plin_spl(kwave),':',label='linear')
plt.annotate("", xy=(kmax, 20), xytext=(kmax, 6),arrowprops=dict(arrowstyle="->",lw=5))
plt.xlabel('$k\,[h\,\mathrm{Mpc}^{-1}]$')
plt.ylabel('$P(k)\,[h^{-3}\mathrm{Mpc}^{3}]$')
plt.legend()
