In this notebook, we explore the functionality of the FFTPower
algorithm, which can compute the 1D power spectrum $P(k)$, 2D power spectrum $P(k,\mu)$, and multipoles $P_\ell(k)$. The algorithm is suitable for use on data sets in periodic simulation boxes, as the power spectrum is computed via a single FFT of the density mesh.
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
from nbodykit.lab import *
from nbodykit import setup_logging, style
import matplotlib.pyplot as plt
plt.style.use(style.notebook)
setup_logging()
We start by generating a mock catalog of biased objects ($b_1 = 2$ ) at a redshif $z=0.55$. We use the Planck 2015 cosmology and the Eisenstein-Hu linear power spectrum fitting formula. We generate the catalog in a box of side length $L = 1380 \ \mathrm{Mpc}/h$ with a constant number density $\bar{n} = 3 \times 10^{-3} \ h^{3} \mathrm{Mpc}^{-3}$.
redshift = 0.55
cosmo = cosmology.Planck15
Plin = cosmology.LinearPower(cosmo, redshift, transfer='EisensteinHu')
b1 = 2.0
cat = LogNormalCatalog(Plin=Plin, nbar=3e-4, BoxSize=1380., Nmesh=256, bias=b1, seed=42)
We update the Position
column to add redshift-space distortions along the z
axis of the box using the VelocityOffset
column.
# add RSD
line_of_sight = [0,0,1]
cat['RSDPosition'] = cat['Position'] + cat['VelocityOffset'] * line_of_sight
In this section, we compute and plot the 1D power spectrum $P(k)$ of the log-normal mock.
We must first convert our CatalogSource
object to a MeshSource
, by setting up the mesh and specifying which interpolation kernel we wish to use. Here, we use "TSC" interpolation, and specify via compensated=True
that we wish to correct for the effects of the interpolation window in Fourier space.
# convert to a MeshSource, using TSC interpolation on 256^3 mesh
mesh = cat.to_mesh(window='tsc', Nmesh=256, compensated=True, position='RSDPosition')
We compute the 1D power spectrum by specifying mode
as "1d". We also choose the desired linear k
binning by specifying the bin spacing via dk
and the minimum k
value via kmin
.
# compute the power, specifying desired linear k-binning
r = FFTPower(mesh, mode='1d', dk=0.005, kmin=0.01)
[ 000010.90 ] 0: 10-30 20:44 CatalogMesh INFO painted 787121 out of 787121 objects to mesh [ 000010.90 ] 0: 10-30 20:44 CatalogMesh INFO mean particles per cell is 0.0469161 [ 000010.91 ] 0: 10-30 20:44 CatalogMesh INFO sum is 787120 [ 000010.91 ] 0: 10-30 20:44 CatalogMesh INFO normalized the convention to 1 + delta [ 000011.25 ] 0: 10-30 20:44 CatalogMesh INFO field: (LogNormalCatalog(seed=42, bias=2) as CatalogMesh) painting done
The result is computed when initializing the FFTPower
class and the power spectrum results are stored as a BinnedStatistic
object as the power
attribute.
# the result is stored at "power" attribute
Pk = r.power
print(Pk)
<BinnedStatistic: dims: (k: 115), variables: ('k', 'power', 'modes')>
The coords
attribute of the BinnedStatistic
object specifies the coordinate grid for the binned result, which in this case, is just the center values of the k
bins.
By default, the power is computed up to the 1D Nyquist frequency, which is defined as $k_\mathrm{Nyq} = \pi N_\mathrm{mesh} / L_\mathrm{box}$, which in this case is equal to $k_\mathrm{Nyq} = \pi \cdot 256 / 1380 = 0.5825 \ h \ \mathrm{Mpc}^{-1}$.
print(Pk.coords)
{'k': array([ 0.0125, 0.0175, 0.0225, 0.0275, 0.0325, 0.0375, 0.0425, 0.0475, 0.0525, 0.0575, 0.0625, 0.0675, 0.0725, 0.0775, 0.0825, 0.0875, 0.0925, 0.0975, 0.1025, 0.1075, 0.1125, 0.1175, 0.1225, 0.1275, 0.1325, 0.1375, 0.1425, 0.1475, 0.1525, 0.1575, 0.1625, 0.1675, 0.1725, 0.1775, 0.1825, 0.1875, 0.1925, 0.1975, 0.2025, 0.2075, 0.2125, 0.2175, 0.2225, 0.2275, 0.2325, 0.2375, 0.2425, 0.2475, 0.2525, 0.2575, 0.2625, 0.2675, 0.2725, 0.2775, 0.2825, 0.2875, 0.2925, 0.2975, 0.3025, 0.3075, 0.3125, 0.3175, 0.3225, 0.3275, 0.3325, 0.3375, 0.3425, 0.3475, 0.3525, 0.3575, 0.3625, 0.3675, 0.3725, 0.3775, 0.3825, 0.3875, 0.3925, 0.3975, 0.4025, 0.4075, 0.4125, 0.4175, 0.4225, 0.4275, 0.4325, 0.4375, 0.4425, 0.4475, 0.4525, 0.4575, 0.4625, 0.4675, 0.4725, 0.4775, 0.4825, 0.4875, 0.4925, 0.4975, 0.5025, 0.5075, 0.5125, 0.5175, 0.5225, 0.5275, 0.5325, 0.5375, 0.5425, 0.5475, 0.5525, 0.5575, 0.5625, 0.5675, 0.5725, 0.5775, 0.5825])}
The input parameters to the algorithm, as well as the meta-data computed during the calculation, are stored in the attrs
dictionary attribute. A key attribute is the Poisson shot noise, stored as the "shotnoise" key.
# print out the meta-data
for k in Pk.attrs:
print("%s = %s" %(k, str(Pk.attrs[k])))
Nmesh = [256 256 256] BoxSize = [ 1380. 1380. 1380.] dk = 0.005 kmin = 0.01 Lx = 1380.0 Ly = 1380.0 Lz = 1380.0 volume = 2628072000.0 mode = 1d los = [0, 0, 1] Nmu = 1 poles = [] N1 = 787121 N2 = 787121 shotnoise = 3338.84116927
Now, we plot the 1D power, first subtracting out the shot noise. Note that the power is complex, so we only plot the real part.
# print the shot noise subtracted P(k)
plt.loglog(Pk['k'], Pk['power'].real - Pk.attrs['shotnoise'])
# format the axes
plt.xlabel(r"$k$ [$h \ \mathrm{Mpc}^{-1}$]")
plt.ylabel(r"$P(k)$ [$h^{-3}\mathrm{Mpc}^3$]")
plt.xlim(0.01, 0.6)
(0.01, 0.6)
In this section, we compute and plot the 2D power spectrum $P(k,\mu)$. For illustration, we compute results using a line-of-sight that is both parallel and perpendicular to the direction that we added redshift-space distortions.
Here, we compute $P(k,\mu)$ where $\mu$ is defined with respect to the z
axis (los=[0,0,1]
), using 5 $\mu$ bins ranging from $\mu=0$ to $\mu=1$.
# compute the 2D power
r = FFTPower(mesh, mode='2d', dk=0.005, kmin=0.01, Nmu=5, los=[0,0,1])
Pkmu = r.power
print(Pkmu)
[ 000014.28 ] 0: 10-30 20:44 CatalogMesh INFO painted 787121 out of 787121 objects to mesh [ 000014.29 ] 0: 10-30 20:44 CatalogMesh INFO mean particles per cell is 0.0469161 [ 000014.29 ] 0: 10-30 20:44 CatalogMesh INFO sum is 787120 [ 000014.29 ] 0: 10-30 20:44 CatalogMesh INFO normalized the convention to 1 + delta [ 000014.71 ] 0: 10-30 20:44 CatalogMesh INFO field: (LogNormalCatalog(seed=42, bias=2) as CatalogMesh) painting done
<BinnedStatistic: dims: (k: 115, mu: 5), variables: ('k', 'mu', 'power', 'modes')>
The coords
attribute of the result now gives the centers of both the k
and mu
bins.
print(Pkmu.coords)
{'k': array([ 0.0125, 0.0175, 0.0225, 0.0275, 0.0325, 0.0375, 0.0425, 0.0475, 0.0525, 0.0575, 0.0625, 0.0675, 0.0725, 0.0775, 0.0825, 0.0875, 0.0925, 0.0975, 0.1025, 0.1075, 0.1125, 0.1175, 0.1225, 0.1275, 0.1325, 0.1375, 0.1425, 0.1475, 0.1525, 0.1575, 0.1625, 0.1675, 0.1725, 0.1775, 0.1825, 0.1875, 0.1925, 0.1975, 0.2025, 0.2075, 0.2125, 0.2175, 0.2225, 0.2275, 0.2325, 0.2375, 0.2425, 0.2475, 0.2525, 0.2575, 0.2625, 0.2675, 0.2725, 0.2775, 0.2825, 0.2875, 0.2925, 0.2975, 0.3025, 0.3075, 0.3125, 0.3175, 0.3225, 0.3275, 0.3325, 0.3375, 0.3425, 0.3475, 0.3525, 0.3575, 0.3625, 0.3675, 0.3725, 0.3775, 0.3825, 0.3875, 0.3925, 0.3975, 0.4025, 0.4075, 0.4125, 0.4175, 0.4225, 0.4275, 0.4325, 0.4375, 0.4425, 0.4475, 0.4525, 0.4575, 0.4625, 0.4675, 0.4725, 0.4775, 0.4825, 0.4875, 0.4925, 0.4975, 0.5025, 0.5075, 0.5125, 0.5175, 0.5225, 0.5275, 0.5325, 0.5375, 0.5425, 0.5475, 0.5525, 0.5575, 0.5625, 0.5675, 0.5725, 0.5775, 0.5825]), 'mu': array([ 0.1, 0.3, 0.5, 0.7, 0.9])}
We plot each the power for each of the 5 $\mu$ bins, and we see the effects of redshift-space distortions as a function of $\mu$.
# plot each mu bin
for i in range(Pkmu.shape[1]):
Pk = Pkmu[:,i] # select the ith mu bin
label = r'$\mu$=%.1f' % (Pkmu.coords['mu'][i])
plt.loglog(Pk['k'], Pk['power'].real - Pk.attrs['shotnoise'], label=label)
# format the axes
plt.legend(loc=0, ncol=2)
plt.xlabel(r"$k$ [$h \ \mathrm{Mpc}^{-1}$]")
plt.ylabel(r"$P(k, \mu)$ [$h^{-3}\mathrm{Mpc}^3$]")
plt.xlim(0.01, 0.6)
(0.01, 0.6)
Now, we specify the line-of-sight as the x
axis and again compute the 2D power.
r = FFTPower(mesh, mode='2d', dk=0.005, kmin=0.01, Nmu=5, los=[1,0,0])
Pkmu = r.power
[ 000018.03 ] 0: 10-30 20:44 CatalogMesh INFO painted 787121 out of 787121 objects to mesh [ 000018.03 ] 0: 10-30 20:44 CatalogMesh INFO mean particles per cell is 0.0469161 [ 000018.03 ] 0: 10-30 20:44 CatalogMesh INFO sum is 787120 [ 000018.03 ] 0: 10-30 20:44 CatalogMesh INFO normalized the convention to 1 + delta [ 000018.45 ] 0: 10-30 20:44 CatalogMesh INFO field: (LogNormalCatalog(seed=42, bias=2) as CatalogMesh) painting done
Again we plot the 2D power for each of the $\mu$ bins, but now that $\mu$ is not defined along the same axis as the redshift-space distortions, we measure isotropic power as a function of $\mu$.
# plot each mu bin
for i in range(Pkmu.shape[1]):
Pk = Pkmu[:,i] # select the ith mu bin
label = r'$\mu$=%.1f' % (Pkmu.coords['mu'][i])
plt.loglog(Pk['k'], Pk['power'].real - Pk.attrs['shotnoise'], label=label)
# format the axes
plt.legend(loc=0, ncol=2)
plt.xlabel(r"$k$ [$h \ \mathrm{Mpc}^{-1}$]")
plt.ylabel(r"$P(k, \mu)$ [$h^{-3}\mathrm{Mpc}^3$]")
plt.xlim(0.01, 0.6)
(0.01, 0.6)
In this section, we also measure the power spectrum multipoles, $P_\ell$, which projects the 2D power on to a basis defined by Legendre weights. The desired multipole numbers $\ell$ should be specified as the poles
keyword.
# compute the 2D power AND ell=0,2,4 multipoles
r = FFTPower(mesh, mode='2d', dk=0.005, kmin=0.01, Nmu=5, los=[0,0,1], poles=[0,2,4])
[ 000021.70 ] 0: 10-30 20:44 CatalogMesh INFO painted 787121 out of 787121 objects to mesh [ 000021.70 ] 0: 10-30 20:44 CatalogMesh INFO mean particles per cell is 0.0469161 [ 000021.70 ] 0: 10-30 20:44 CatalogMesh INFO sum is 787120 [ 000021.70 ] 0: 10-30 20:44 CatalogMesh INFO normalized the convention to 1 + delta [ 000022.04 ] 0: 10-30 20:44 CatalogMesh INFO field: (LogNormalCatalog(seed=42, bias=2) as CatalogMesh) painting done
Now, there is an additional attribute poles
which stores the $P_\ell(k)$ result.
poles = r.poles
print(poles)
print("variables = ", poles.variables)
<BinnedStatistic: dims: (k: 115), variables: 5 total> variables = ['k', 'power_0', 'power_2', 'power_4', 'modes']
We plot each multipole, subtracting the shot noise only from the monopole $P_0$. The multipoles are stored using the variable names "power_0", "power_2", etc for $\ell=0,2,$ etc.
for ell in [0, 2, 4]:
label = r'$\ell=%d$' % (ell)
P = poles['power_%d' %ell].real
if ell == 0: P = P - poles.attrs['shotnoise']
plt.plot(poles['k'], poles['k'] * P, label=label)
# format the axes
plt.legend(loc=0)
plt.xlabel(r"$k$ [$h \ \mathrm{Mpc}^{-1}$]")
plt.ylabel(r"$k \ P_\ell$ [$h^{-2} \mathrm{Mpc}^2$]")
plt.xlim(0.01, 0.6)
(0.01, 0.6)