Gridding/interpolation

Contents

  • grid_points: Returns a map representing a sum of delta functions (or a “galaxies - randoms” difference map).

  • interpolate_points: Interpolates real-space map at a specified set of points.

  • apply_kernel_compensation: Modifies Fourier-space map ‘arr’ in-place, to debias interpolation/gridding.


kszx.grid_points(box, points, weights=None, rpoints=None, rweights=None, kernel=None, fft=False, spin=0, periodic=False, compensate=False, wscal=1.0)

Returns a map representing a sum of delta functions (or a “galaxies - randoms” difference map).

Function args:

  • box (kszx.Box): defines pixel size, bounding box size, and location of observer. See Box for more info.

  • points (2-d array): Sequence of points where the delta functions are located (usually galaxy locations). Array shape should be (n,d), where n is the number of interpolation points, and d is the box dimension (usually 3).

  • weights (either scalar, None, or 1-d array).

    • if weights is an array, then it should have length n, where n is the number of interpolation points.

    • if weights is a scalar, then all delta functions have equal weight.

    • if weights is None, then all delta functions have weight 1.

  • rpoints (2-d array or None): Optional sequence of points where delta functions with negative coefficients are located (representing a “random” catalog).

  • rweights (either scalar, None, or 1-d array). Weights for the rpoints array (same semantics as weights).

    NOTE: an additional multiplicative normalization is applied to the rweights so that sum(rweights) = -sum(weights). That is, if rweights are specified, then the (galaxies - randoms) maps returned by this function always has mean zero.

  • kernel (string): either 'cic' or 'cubic' (more options will be defined later).

  • fft (boolean): if True, then fft_r2c(.., spin) will be applied to the output map after gridding.

  • spin (integer): passed as ‘spin’ argument to fft_r2c(). (Only used if fft=True.)

  • compensate (boolean): If true, then kszx.apply_kernel_compensation() is called after fft_r2c(). (Only used if fft=True.)

  • periodic (boolean): if True, then the box has periodic boundary conditions.

  • wscal (float, default 1.0): overall scaling applied to weights.

Return value:

  • A numpy array representing a real-space (fft=False) or Fourier-space (fft=True) map.

    The real-space and Fourier-space array shapes are given by box.real_space_shape and box.fourier_space_shape, and are related as follows:

    \[\begin{split}\begin{align} (\mbox{real-space shape}) &= (n_0, n_1, \cdots, n_{d-1}) \\ (\mbox{Fourier-space shape}) &= (n_0, n_1, \cdots, \lfloor n_{d-1}/2 \rfloor + 1) \end{align}\end{split}\]

Notes

  • The points array should be specified in “observer coordinates”, not “grid coordinates”.

    (Reminder: in observer coordinates, the observer is at the origin, coordinates have units of distance, and the corners of the box are at box.{lpos,rpos}. See Box for more info.)

  • The normalization of the output map includes a factor (pixel volume)\(^{-1}\). (That is, the sum of the output array is \((\sum_j w_j)/V_{\rm pix}\), not \((\sum_j w_j)\).)

    This normalization best represents a weighted sum of delta functions \(f(x) = \sum_j w_j \delta^3(x-x_j)\). For example, if we FFT the output array with kszx.fft_r2c(), the result is a weighted sum of plane waves:

    \[f(k) = \sum_j w_j \exp(-i{\bf k}\cdot {\bf x_j})\]

    with no factor of box or pixel volume.

  • After calling grid_points(), you may want to call apply_kernel_compensation() to mitigate high-\(k\) biases. See apply_kernel_compensation() docstring for more info.


kszx.interpolate_points(box, arr, points, kernel, fft=False, spin=0, periodic=False)

Interpolates real-space map at a specified set of points.

Function args:

  • box (kszx.Box): defines pixel size, bounding box size, and location of observer. See Box for more info.

  • arr (numpy array): represents the map to be interppolated, either in real space (if fft=False) or Fourier space (if fft=True).

    The real-space and Fourier-space array shapes are given by box.real_space_shape and box.fourier_space_shape, and are related as follows:

\[\begin{split}\begin{align} (\mbox{real-space shape}) &= (n_0, n_1, \cdots, n_{d-1}) \\ (\mbox{Fourier-space shape}) &= (n_0, n_1, \cdots, \lfloor n_{d-1}/2 \rfloor + 1) \end{align}\end{split}\]

  • points (numpy array): Sequence of points where the map is to be interpolated. Array shape should be (n,d), where n is the number of interpolation points, and d is the box dimension (usually 3).

  • kernel (string): either 'cic' or 'cubic' (more options will be defined later).

  • fft (boolean): if True, then arr is a Fourier-space map, and fft_c2r(arr, spin) will be called before interpolating.

  • spin (integer): passed as ‘spin’ argument to fft_c2r(). (Only used if fft=True.)

  • periodic (boolean): if True, then the box has periodic boundary conditions.

Return value:

  • 1-d numpy array with length npoints, containing interpolated values.

Notes

  • The points array should be specified in “observer coordinates”, not “grid coordinates”.

    (Reminder: in observer coordinates, the observer is at the origin, coordinates have units of distance, and the corners of the box are at box.{lpos,rpos}. See Box for more info.)

  • Before calling interpolate_points(), you may want to call apply_kernel_compensation() to mitigate high-\(k\) biases. See apply_kernel_compensation() docstring for more info.


kszx.apply_kernel_compensation(box, arr, kernel, exponent=-0.5)

Modifies Fourier-space map ‘arr’ in-place, to debias interpolation/gridding.

Context: gridding kernels (see grid_points()) multiplicatively bias power spectrum estimation,

\[<P(k)>_{\rm estimated} = C(k) \, P(k)_{true}\]

Here, \(C(k)\) is a “compensation factor” satisfying \(0 \le C(k) \le 1\) which depends on both the magnitude and orientation of \(k\).

There is a similar bias which pertains to interpolation kernels, rather than gridding kernels (see interpolate_points()). Suppose we start with a Fourier-space map \(f(k)\), then Fourier transform and interpolate at random locations. One would expect that an interpolated value \(f_{\rm interp}\) has variance

\[\langle f_{\rm interp}^2 \rangle = \int \frac{d^3k}{(2\pi)^3} \, f(k)^2\]

However, the interpolation kernel produces a bias: the actual variance is

\[\langle f_{\rm interp}^2 \rangle = \int \frac{d^3k}{(2\pi)^3} \, C(k) f(k)^2\]

The function apply_kernel_compensation multiplies Fourier-space map arr in-place by \(C(k)^p\), where \(p\) is the exponent argument. Here are two common applications:

  1. Before calling estimate_power_spectrum() on one or more Fourier-space maps, you should call apply_kernel_compensation() on each map, to multiply by \(C(k)^{-1/2}\). This will mitigate the power spectrum estimation bias noted above.

  2. Before calling interpolate_points() on a map, you should call apply_kernel_compensation() on the map, to multiply by \(C(k)^{-1/2}\). This will mitigate the interpolation bias noted above. (This assumes that you start with the map in Fourier space, and FFT before interpolating.)

Function args:

  • box (kszx.Box): defines pixel size, bounding box size, and location of observer. See Box for more info.

  • arr: numpy array representing a Fourier-space map. The array shape should be given by box.fourier_space_shape and the dtype should be complex, see note below.

  • kernel (string): either 'cic' or 'cubic' (more options will be defined later).

  • exponent (float): array will be multiplied by C(k)**exponent. (The default value is exponent = -0.5, since this value arises in both applications above.)

Return value: None (the arr argument is modified in-place, by multiplying by C(k)**exponent).

Reminder: real-space and Fourier-space array shapes are given by box.real_space_shape and box.fourier_space_shape, and are related as follows:

\[\begin{split}\begin{align} (\mbox{real-space shape}) &= (n_0, n_1, \cdots, n_{d-1}) \\ (\mbox{Fourier-space shape}) &= (n_0, n_1, \cdots, \lfloor n_{d-1}/2 \rfloor + 1) \end{align}\end{split}\]