Simulation Box Power Spectrum/Multipoles (FFTPower
)¶
The FFTPower
class computes the 1d power spectrum \(P(k)\), 2d
power spectrum \(P(k,\mu)\), and/or multipoles \(P_\ell(k)\) for data
in a simulation box, using a Fast Fourier Transform (FFT). Here, we provide
a brief overview of the algorithm itself as well as the key things to know for
the user to get up and running quickly.
Note
To jump right into the FFTPower
algorithm, see this
cookbook recipe for a detailed
walk-through of the FFTPower
algorithm.
The Algorithm¶
The steps involved in computing the power spectrum via FFTPower
are as follows:
Generate data on a mesh
Data must be painted on to a discrete mesh to compute the power spectrum. There are several ways to generate data on a mesh (see Creating a Mesh), but the most common is painting a discrete catalog of objects on to a mesh (see Converting a CatalogSource to a Mesh and Painting Catalogs to a Mesh). The
FFTPower
class accepts input data in either the form of aMeshSource
or aCatalogSource
. In the latter case, the catalog is automatically converted to a mesh using the default parameters of theto_mesh()
function.When converting from a catalog to a mesh, users can customize the painting procedure via the options of the
to_mesh()
function. These options have important effects on the resulting power spectrum of the field in Fourier space. See Converting a CatalogSource to a Mesh for more details.FFT the mesh to Fourier space
Once the density field is painted to the mesh, the Fourier transform of the field \(\delta(\vx)\) is performed in parallel to obtain the complex modes of the overdensity field, \(\delta(\vk)\). The field is stored using the
ComplexField
object.Generate the 3D power spectrum on the mesh
The 3D power spectrum field is computed on the mesh, using
\[P(\mathbf{k}) = \delta(\mathbf{k}) \cdot \delta^\star(\mathbf{k}),\]where \(\delta^\star (\mathbf{k})\) is the complex conjugate of \(\delta(\mathbf{k})\).
Perform the binning in the specified basis
Finally, the 3D power defined on the mesh \(P(\mathbf{k})\) is binned using the basis specified by the user. The available options for binning are:
1D binning as a function of wavenumber \(k\)
2D binning as a function of wavenumber \(k\) and cosine of the angle to the line-of-sight \(\mu\)
Multipole binning as a function of \(k\) and multipole number \(\ell\)
The Functionality¶
Users can compute various quantities using the FFTPower
. We’ll discuss
the available functionality briefly in the sub-sections below.
Auto power spectra and cross power spectra¶
Both auto and cross spectra are supported. Users can compute cross power spectra
by passing a second mesh object to the FFTPower
class using
the second
keyword. The first mesh object should always be specified as
the first
argument.
1D Power Spectrum, \(P(k)\)¶
The 1D power spectrum \(P(k)\) can be computed by specifying the
mode
argument as “1d”. The wavenumber binning will be linear, and can be
customized by specifying the dk
and kmin
attributes. By default,
the edge of the last wavenumber bin is the
Nyquist frequency, given
by \(k_\mathrm{Nyq} = \pi N_\mathrm{mesh} / L_\mathrm{box}\). If dk
is not specified, then the fundamental mode of the box is used:
\(2\pi/L_\mathrm{box}\).
2D Power Spectrum, \(P(k,\mu)\)¶
The 2D power spectrum \(P(k,\mu)\) can be computed by specifying the
mode
argument as “2d”. The number of \(\mu\) bins is specified via
the Nmu
keyword. The bins range from \(\mu=0\) to \(\mu=1\).
Multipoles of \(P(k,\mu)\)¶
The FFTPower
class can also compute the multipoles of the 2D power
spectrum, defined as
where \(\mathcal{L}_\ell\) is the Legendre polynomial of order
\(\ell\). Users can specify which multipoles they wish to compute
by passing a list of the desired \(\ell\) values as the poles
keyword to the FFTPower
class.
For example, we can compute both \(P(k,\mu)\) and \(P_\ell(k)\) for a uniform catalog of objects using:
[2]:
from nbodykit.lab import UniformCatalog, FFTPower
cat = UniformCatalog(nbar=100, BoxSize=1.0, seed=42)
r = FFTPower(cat, mode='2d', Nmesh=32, Nmu=5, poles=[0,2,4])
/home/yfeng1/anaconda3/install/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
from ._conv import register_converters as _register_converters
The Results¶
The power spectrum results are stored in two attributes of the
initialized FFTPower
object:
power
and poles
. These attributes are
BinnedStatistic
objects, which
behave like structured numpy arrays and store
the measured results on a coordinate grid defined by the bins.
See Analyzing your Results for a full tutorial on using
the BinnedStatistic
class.
The power
attribute stores the following variables:
- k :
the mean value for each
k
bin
- muif
mode=2d
the mean value for each
mu
bin
- muif
- power :
complex array storing the real and imaginary components of the power
- modes :
the number of Fourier modes averaged together in each bin
The poles
attribute stores the following variables:
- k :
the mean value for each
k
bin
- power_L :
complex array storing the real and imaginary components for the \(\ell=L\) multipole
- modes :
the number of Fourier modes averaged together in each bin
Note that measured power results for bins where modes
is zero (no data points
to average over) are set to NaN
.
In our example, the power
and poles
attributes are:
[3]:
# the 2D power spectrum results
print("power = ", r.power)
print("variables = ", r.power.variables)
for name in r.power.variables:
var = r.power[name]
print("'%s' has shape %s and dtype %s" %(name, var.shape, var.dtype))
power = <BinnedStatistic: dims: (k: 16, mu: 5), variables: ('k', 'mu', 'power', 'modes')>
variables = ['k', 'mu', 'power', 'modes']
'k' has shape (16, 5) and dtype float64
'mu' has shape (16, 5) and dtype float64
'power' has shape (16, 5) and dtype complex128
'modes' has shape (16, 5) and dtype int64
[4]:
# the multipole results
print("poles = ", r.poles)
print("variables = ", r.poles.variables)
for name in r.poles.variables:
var = r.poles[name]
print("'%s' has shape %s and dtype %s" %(name, var.shape, var.dtype))
poles = <BinnedStatistic: dims: (k: 16), variables: 5 total>
variables = ['k', 'power_0', 'power_2', 'power_4', 'modes']
'k' has shape (16,) and dtype float64
'power_0' has shape (16,) and dtype complex128
'power_2' has shape (16,) and dtype complex128
'power_4' has shape (16,) and dtype complex128
'modes' has shape (16,) and dtype int64
These attributes also store meta-data computed during the power calculation
in the attrs
dictionary. Most importantly, the shotnoise
key
gives the Poisson shot noise, \(P_\mathrm{shot} = V / N\), where V
is the volume of the simulation box and N is the number of objects. The keys
N1
and N2
store the number of objects.
In our example, the meta-data is:
[5]:
for k in r.power.attrs:
print("%s = %s" %(k, str(r.power.attrs[k])))
Nmesh = [32 32 32]
BoxSize = [1. 1. 1.]
Lx = 1.0
Ly = 1.0
Lz = 1.0
volume = 1.0
mode = 2d
los = [0, 0, 1]
Nmu = 5
poles = [0, 2, 4]
dk = 6.283185307179586
kmin = 0.0
N1 = 96
N2 = 96
shotnoise = 0.010416666666666666
Note
The shot noise is not subtracted from any measured results. Users can
access the Poisson shot noise value in the meta-data attrs
dictionary.
Saving and Loading¶
Results can easily be saved and loaded from disk in a reproducible manner
using the FFTPower.save()
and FFTPower.load()
functions.
The save
function stores the state of the algorithm,
including the meta-data in the FFTPower.attrs
dictionary, in a
JSON plaintext format.
[6]:
# save to file
r.save("fftpower-example.json")
# load from file
r2 = FFTPower.load("fftpower-example.json")
print("power = ", r2.power)
print("poles = ", r2.poles)
print("attrs = ", r2.attrs)
power = <BinnedStatistic: dims: (k: 16, mu: 5), variables: ('k', 'mu', 'power', 'modes')>
poles = <BinnedStatistic: dims: (k: 16), variables: 5 total>
attrs = {'Nmesh': array([32, 32, 32]), 'BoxSize': array([1., 1., 1.]), 'Lx': 1.0, 'Ly': 1.0, 'Lz': 1.0, 'volume': 1.0, 'mode': '2d', 'los': [0, 0, 1], 'Nmu': 5, 'poles': [0, 2, 4], 'dk': 6.283185307179586, 'kmin': 0.0, 'N1': 96, 'N2': 96, 'shotnoise': 0.010416666666666666}
Common Pitfalls¶
The default configuration of nbodykit should lead to reasonable results
when using the FFTPower
algorithm. When performing custom, more complex
analyses, some of the more common pitfalls are:
When the results of
FFTPower
do not seem to make sense, the most common culprit is usually the configuration of the mesh, and whether or not the mesh is “compensated”. In the language of nbodykit, “compensated” refers to whether the effects of the interpolation window used to paint the density field have been de-convolved in Fourier space. See the Converting a CatalogSource to a Mesh section for detailed notes on this procedure.Be wary of normalization issues when painting weighted density fields. See Painting Catalogs to a Mesh for further details on painting meshes and Applying Functions to the Mesh for notes on applying arbitrary functions to the mesh while painting. See this cookbook recipe for examples of more advanced painting uses.