abacusnbody.analysis

power_spectrum module

The power spectrum module contains various useful tools for computing power spectra (wedges, multipoles and beyond) in the cubic box.

abacusnbody.analysis.power_spectrum.calc_power(pos, Lbox, kbins=None, mubins=None, k_max=None, logk=False, paste='TSC', nmesh=128, compensated=True, interlaced=True, w=None, pos2=None, w2=None, poles=None, squeeze_mu_axis=True, nthread=2, dtype=<class 'numpy.float32'>)

Compute the 3D power spectrum given particle positions by first painting them on a cubic mesh and then applying Fourier transforms and mode counting. Outputs (k,mu) wedges by default; can also output Legendre multipoles.

posarray_like

particle positions, shape (N,3)

Lboxfloat

box size of the simulation.

kbinsint, array_like, or None, optional

An int indicating the number of bins of k, which ranges from 0 to k_max if logk, and 2pi/L to k_max if not. Or an array-like, which will be used as-is. Default is None, which sets kbins to nmesh.

mubinsint, None, or array_like, optional

An int indicating the number of bins of mu. mu ranges from 0 to 1. Or an array-like of bin edges, which will be used as-is. Default of None sets mubins to 1.

k_maxfloat, optional

maximum k wavenumber. Default is None, which sets k_max to k_Nyquist of the mesh.

logkbool, optional

Logarithmic or linear k bins. Ignored if kbins is array-like. Default is False.

pastestr, optional

particle pasting approach (CIC or TSC). Default is ‘TSC’.

nmeshint, optional

size of the 3d array along x and y dimension. Default is 128.

compensatedbool, optional

want to apply first-order compensated filter? Default is True.

interlacedbool, optional

want to apply interlacing? Default is True.

warray_like, optional

weights for each particle.

pos2array_like, optional

second set of particle positions, shape (N,3)

polesNone or list of int, optional

Legendre multipoles of the power spectrum or correlation function. Default of None gives the monopole.

squeeze_mu_axisbool, optional

Remove the mu axis from the output arrays if it has length 1. Default: True

nthreadint, optional

Number of numba threads to use

dtypenp.dtype, optional

Data type of the field

Returns:

power – The power spectrum in an astropy Table of length nbins_k. The columns are:

• k_mid: arithmetic bin centers of the k wavenumbers, shape (nbins_k,)

• k_avg: mean wavenumber per (k, mu) wedge, shape (nbins_k,nbins_mu)

• mu_mid: arithmetic bin centers of the mu angles, shape (nbins_k,nbins_mu)

• power: mean power spectrum per (k, mu) wedge, shape (nbins_k,nbins_mu)

• N_mode: number of modes per (k, mu) wedge, shape (nbins_k,nbins_mu)

If multipoles are requested via poles, the table includes:

• poles: mean Legendre multipole coefficients, shape (nbins_k,len(poles))

• N_mode_poles: number of modes per pole, shape (nbins_k,len(poles))

The meta field of the table will have metadata about the power spectrum.

Return type:

astropy.Table

abacusnbody.analysis.power_spectrum.calc_pk_from_deltak(field_fft, Lbox, k_bin_edges, mu_bin_edges, field2_fft=None, poles=array([], dtype=int64), squeeze_mu_axis=True, nthread=2)

Calculate the power spectrum of a given Fourier field, with binning in (k,mu). Optionally computes Legendre multipoles in k bins.

Parameters:

• field_fft (array_like) – Fourier 3D field.

• Lbox (float) – box size of the simulation.

• k_bin_edges (array_like) – edges of the k wavenumbers.

• mu_bin_edges (array_like) – edges of the mu angles.

• field2_fft (array_like, optional) – second Fourier 3D field, used in cross-correlation.

• poles (np.ndarray, optional) – Legendre multipoles of the power spectrum or correlation function. Probably has to be a Numpy array, or Numba will complain.

• squeeze_mu_axis (bool, optional) – Remove the mu axis from the output arrays if it has length 1. Default: True

• nthread (int, optional) – Number of numba threads to use

Returns:

• power (array_like) – mean power spectrum per (k, mu) wedge.

• N_mode (array_like) – number of modes per (k, mu) wedge.

• binned_poles (array_like) – mean power spectrum per k for each Legendre multipole.

• N_mode_poles (array_like) – number of modes per k.

• k_avg (array_like) – mean wavenumber per (k, mu) wedge.

abacusnbody.analysis.power_spectrum.pk_to_xi(Pk, Lbox, r_bins, poles=[0, 2, 4])

Transform 3D power spectrum into correlation function multipoles.

Parameters:

• Pk (array_like) – 3D power spectrum.

• Lbox (float) – box size of the simulation.

• r_bins (array_like) – r separation bins.

• poles (array_like) – Legendre multipoles of the power spectrum or correlation function.

Returns:

• r_binc (array_like) – r separation bin centers.

• binned_poles (array_like) – correlation function multipoles.

• Npoles (array_like) – number of modes per r bin.

abacusnbody.analysis.power_spectrum.project_3d_to_poles(k_bin_edges, raw_p3d, Lbox, poles)

Project 3D power spectrum into multipoles of the power spectrum.

Parameters:

• k_bin_edges (array_like) – edges of the k wavenumbers.

• raw_p3d (array_like) – array containing the power spectrum modes.

• Lbox (float) – box size of the simulation.

• poles (array_like) – Legendre multipoles of the power spectrum or correlation function.

Returns:

• binned_poles (array_like) – mean power spectrum per k for each Legendre multipole.

• Npoles (array_like) – number of modes per k.

abacusnbody.analysis.power_spectrum.get_k_mu_edges(Lbox, k_max, kbins, mubins, logk)

Obtain bin edges of the k wavenumbers and mu angles.

Parameters:

• Lbox (float) – box size of the simulation.

• k_max (float) – maximum k wavenumber.

• kbins (int or array_like) – An int indicating the number of bins of k, which ranges from 0 to k_max if logk, and 2pi/L to k_max if not. Or an array-like, which will be returned unchanged.

• mubins (int or array_like) – An int indicating the number of bins of mu, which ranges from 0 to 1. Or an array-like, which will be returned unchanged.

• logk (bool) – Logarithmic or linear k bins

Returns:

• k_bin_edges – edges of the k wavenumbers.

• mu_bin_edges – edges of the mu angles.

tpcf_corrfunc module

abacusnbody.analysis.tpcf_corrfunc.tpcf_multipole(s_mu_tcpf_result, mu_bins, order=0)

Calculate the multipoles of the two point correlation function after first computing ~halotools.mock_observables.s_mu_tpcf. This is copied over from halotools. Original author was Duncan Campbell.

Parameters:

• s_mu_tcpf_result (np.ndarray) – 2-D array with the two point correlation function calculated in bins of \(s\) and \(\mu\). See ~halotools.mock_observables.s_mu_tpcf.

• mu_bins (array_like) – array of \(\mu = \cos(\theta_{\rm LOS})\) bins for which s_mu_tcpf_result has been calculated. Must be between [0,1].

• order (int, optional) – order of the multpole returned.

Returns:

xi_l – multipole of s_mu_tcpf_result of the indicated order.

Return type:

np.array

Examples

For demonstration purposes we create a randomly distributed set of points within a periodic cube of length 250 Mpc/h.:

>>> Npts = 100
>>> Lbox = 250.
>>> x = np.random.uniform(0, Lbox, Npts)
>>> y = np.random.uniform(0, Lbox, Npts)
>>> z = np.random.uniform(0, Lbox, Npts)

We transform our x, y, z points into the array shape used by the pair-counter by taking the transpose of the result of numpy.vstack. This boilerplate transformation is used throughout the ~halotools.mock_observables sub-package:

>>> sample1 = np.vstack((x,y,z)).T

First, we calculate the correlation function using ~halotools.mock_observables.s_mu_tpcf.:

>>> from halotools.mock_observables import s_mu_tpcf
>>> s_bins  = np.linspace(0.01, 25, 10)
>>> mu_bins = np.linspace(0, 1, 15)
>>> xi_s_mu = s_mu_tpcf(sample1, s_bins, mu_bins, period=Lbox)

Then, we can calculate the quadrapole of the correlation function:

>>> xi_2 = tpcf_multipole(xi_s_mu, mu_bins, order=2)

abacusnbody.analysis.tpcf_corrfunc.calc_xirppi_fast(x1, y1, z1, rpbins, pimax, pi_bin_size, lbox, Nthread, num_cells=20, x2=None, y2=None, z2=None)

abacusnbody.analysis.tpcf_corrfunc.calc_multipole_fast(x1, y1, z1, sbins, lbox, Nthread, nbins_mu=50, num_cells=20, x2=None, y2=None, z2=None, orders=[0, 2])

abacusnbody.analysis.tpcf_corrfunc.calc_wp_fast(x1, y1, z1, rpbins, pimax, lbox, Nthread, num_cells=30, x2=None, y2=None, z2=None)

tsc module

abacusnbody.analysis.tsc.tsc_parallel(pos, densgrid, box, weights=None, nthread=-1, wrap=True, npartition=None, sort=False, coord=0, verbose=False, offset=0.0)

A parallel implementation of TSC mass assignment using numba. The algorithm partitions the particles into stripes of sufficient width that their TSC clouds don’t overlap, and then does all the even-numbered stripes in parallel, followed by the odd-numbered ones.

This method parallelizes well and can exceed, e.g., 500 million particles per second with 32 cores on a modern processor for a cache-resident grid size. That’s 20 ms for 10^7 particles, which is number density 1e-3 in a 2 Gpc/h box.

The algorithm is expected to be bound by memory bandwidth, rather than CPU. Sometimes using, e.g., half of all CPUs will be faster than using all of them.

npartition is a tuning parameter. Generally it should be at least 2*nthread, so that all threads have work to do in both passes. Sometimes an even finer partitioning can produce a favorable ordering in memory of particles for TSC. Sorting the particles within each stripe produces an even more favorable ordering, but the current implementation of sorting is slow enough that it’s usually not worth it.

Parameters:

• pos (ndarray of shape (n,3)) – The particles, in domain [0,box)

• densgrid (ndarray, tuple, or int) – Either an ndarray in which to write the density, or a tuple/int indicating the shape of the array to allocate. Can be 2D or 3D; ints are interpreted as 3D cubic grids. Anisotropic grids are also supported (nx != ny != nz).

• box (float) – The domain size. Positions are expected in domain [0,box) (but may be wrapped; see wrap).

• weights (ndarray of shape (n,), optional) – Particle weights/masses. Default: None

• nthread (int, optional) – Number of threads, for both the parallel partition and the TSC. Values < 0 use numba.config.NUMBA_NUM_THREADS, which is usually all CPUs. Default: -1

• wrap (bool, optional) – Apply an in-place periodic wrap to any out-of-bounds particle positions, bringing them back to domain [0,box). This is on by default because it’s generally fast compared to TSC. Default: True

• npartition (int, optional) – Number of stripes in which to partition the positions. This is a tuning parameter (with certain constraints on the max value); a value of None will use the default of (no larger than) 2*nthread. Default: None

• sort (bool, optional) – Sort the particles along the coord coordinate within each partition stripe. This can affect performance. Default: False

• coord (int, optional) – The coordinate on which to partition. coord = 0 means x, coord = 1 means y, etc. Default: 0

• verbose (bool, optional) – Print some information about settings and timings. Default: False

Returns:

dens – The density grid, which may be newly allocated, or the same as the input densgrid argument if that argument was an ndarray.

Return type:

ndarray

abacusnbody.analysis.tsc.partition_parallel(pos, npartition, boxsize, weights=None, coord=0, nthread=-1, sort=False)

A parallel partition. Partitions a set of positions into npartition pieces, using the coord coordinate (coord=0 partitions on x, coord=1 partitions on y, etc.).

The particle copy stage is coded as a scatter rather than a gather.

Note that this function is expected to be bound by memory bandwidth rather than CPU.

Parameters:

• pos (ndarray of shape (n,3)) – The positions, in domain [0,boxsize)

• npartition (int) – The number of partitions

• boxsize (float) – The domain of the particles

• weights (ndarray of shape (n,), optional) – Particle weights. Default: None

• coord (int, optional) – The coordinate to partition on. 0 is x, 1 is y, etc. Default: 0 (x coordinate)

• nthread (int, optional) – Number of threads to parallelize over (using Numba threading). Values < 0 use numba.config.NUMBA_NUM_THREADS, which is usually all CPUs. Default: -1

• sort (bool, optional) – Sort the particles on the coord coordinate within each partition. Can speed up subsequent TSC, but generally slow and not worth it. Default: False

Returns:

• partitioned (ndarray like pos) – The particles, in partitioned order

• part_starts (ndarray, shape (npartition + 1,), dtype int64) – The index in partitioned where each partition starts

• wpart (ndarray or None) – The weights, in partitioned order; or None if weights not given.