Mesh Interpolation Windows¶

In this notebook we explore the different window options for interpolating discrete objects on to a mesh. We compute the 1D power spectrum \(P(k)\) of a density field interpolated from a log-normal mock using different windows. The windows we test are:

Cloud in Cell: cic (default in nbodykit)
Triangular Shaped Cloud: tsc
Lanczos Kernel with \(a=2\) and \(a=3\): lanczos2 and lanczos3
Symmetric Daubechies wavelet with sizes 6, 12, and 20: sym6, sym12, and sym20

We also include timing tests when using each of these windows. The CIC window (default) is the fastest and can be considerably faster than the other kernels, especially the wavelet windows with large sizes.

When computing the power spectrum \(P(k)\), we de-convolve the effects of the interpolation on the measured power (compensated=True) for the CIC and TSC windows, but do not apply any corrections for the other windows.

For more information on using the Daubechies wavelets for power spectrum measurements, see Cui et al. 2008.

In [1]:

%matplotlib inline
%config InlineBackend.figure_format = 'retina'

In [2]:

from nbodykit.lab import *
from nbodykit import style
from scipy.interpolate import InterpolatedUnivariateSpline as spline

import matplotlib.pyplot as plt
plt.style.use(style.notebook)

Initalizing a Log-normal Mock¶

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

In [3]:

redshift = 0.55
cosmo = cosmology.Planck15
Plin = cosmology.EHPower(cosmo, redshift)

cat = LogNormalCatalog(Plin=Plin, nbar=3e-4, BoxSize=1380., Nmesh=256, bias=2.0, seed=42)

Generating the “Truth”¶

We generate the “truth” power spectrum by using a higher resolution mesh with Nmesh=512, using the TSC window interpolation method. With a higher resolution mesh, the Nyquist frequency is larger, decreasing the effects of aliasing at a given \(k\) value.

In [4]:

mesh = cat.to_mesh(window='tsc', Nmesh=512, compensated=True)
r = FFTPower(mesh, mode='1d') # hi-resolution mesh
truth = r.power

truth = spline(truth['k'], truth['power'].real - truth.attrs['shotnoise'])

Comparing Different Windows¶

We compare the results when using different windows to measure \(P(k)\) as compared to the high-resolution “truth” measured in the previous section. We compute results using a lower-resultion mesh, Nmesh=256, and the effects of aliasing are clear, causing the measured power to differ from the “truth”. The goal is to identify which windows perform better/worse at minimizing this differences.

In [5]:

for window in ['cic', 'tsc', 'sym6', 'sym12', 'sym20', 'lanczos2', 'lanczos3']:

    print("computing power for window '%s'\n" %window + "-"*32)
    compensated = True if window in ['cic', 'tsc'] else False
    mesh = cat.to_mesh(Nmesh=256, window=window, compensated=compensated, interlaced=False)

    # compute the power
    %time r = FFTPower(mesh, mode='1d')
    Pk = r.power

    plt.plot(Pk['k'], (Pk['power'].real - Pk.attrs['shotnoise']) / truth(Pk['k']), label=window)

plt.legend(loc=0, ncol=2, fontsize=16)
plt.xlabel(r"$k$ [$h \ \mathrm{Mpc}^{-1}$]")
plt.ylabel(r"$P(k) / P(k)^\mathrm{truth}$")
plt.xlim(0.01, 0.6)
plt.ylim(0.9, 1.1)

computing power for window 'cic'
--------------------------------
CPU times: user 1.38 s, sys: 96.5 ms, total: 1.48 s
Wall time: 1.48 s
computing power for window 'tsc'
--------------------------------
CPU times: user 1.41 s, sys: 95.2 ms, total: 1.5 s
Wall time: 1.5 s
computing power for window 'sym6'
--------------------------------
CPU times: user 4.57 s, sys: 102 ms, total: 4.67 s
Wall time: 4.68 s
computing power for window 'sym12'
--------------------------------
CPU times: user 11 s, sys: 121 ms, total: 11.1 s
Wall time: 11.1 s
computing power for window 'sym20'
--------------------------------
CPU times: user 19.1 s, sys: 274 ms, total: 19.4 s
Wall time: 19.8 s
computing power for window 'lanczos2'
--------------------------------
CPU times: user 1.92 s, sys: 97.8 ms, total: 2.02 s
Wall time: 2.03 s
computing power for window 'lanczos3'
--------------------------------
CPU times: user 3.56 s, sys: 124 ms, total: 3.69 s
Wall time: 3.78 s

Out[5]:

(0.9, 1.1)

../_images/cookbook_interpolation-windows_8_2.png

In this figure, we see that each window kernel has a specific \(k_\mathrm{max}\) where we can reasonably trust the measured results. The wavelet kernels sym12 and sym20 perform the best but taking over an order magnitude longer to compute than the default kernel (CIC). The Triangular Shaped Cloud window performs the next best, followed by the sym6 kernel and the CIC kernel. The Lanczos kernels do not perform well and should not be used by the user without additional compensation corrections (currently not implemented by default in nbodykit).

Recommendations

Given the increased speed costs of the wavelet windows, we recommend users use either the CIC or TSC windows. For optimized solutions providing the most precise results in the fastest amount of time, users should test the effects of setting interlaced=True (see this tutorial) while also trying various mesh sizes by changing the Nmesh parameter. The most robust results are obtained by running a series of convergence tests by computing results with a high-resolution mesh, which reduces the contributions of aliasing at fixed \(k\). This allows the user to robustly determine the best configuration for their desired accuracy.