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:

  1. 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 a MeshSource or a CatalogSource. In the latter case, the catalog is automatically converted to a mesh using the default parameters of the to_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.

  2. 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.

  3. 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})\).

  4. 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

\[P_\ell(k) = (2\ell + 1) \int_0^1 d\mu P(k,\mu) \mathcal{L}_\ell(\mu),\]

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

  • 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.