pypowspec

Python (Cython) bindings for powspec by Cheng Zhao - a fast C/FFTW code for auto- and cross-power-spectrum multipoles of galaxy catalogs and simulation boxes.

pypowspec feeds the C engine in-memory NumPy arrays instead of ASCII files and returns the results as dictionaries/arrays, with all settings available as Python keyword arguments:

import numpy as np
import pypowspec

pk = pypowspec.compute_auto_box(x, y, z,                # coordinates in [0, L)
                                box_size=1000.0,        # Mpc/h
                                grid_size=256,
                                poles=(0, 2, 4),
                                assign="TSC",
                                interlace=True)

pk.k                     # bin centres        pk.multipoles   # (nbin, nl)
pk.shot_noise            # estimator metadata pk.nmodes       # modes per bin

It also extends the original code with power spectra of simulation boxes with random catalogs (e.g. shifted randoms after BAO reconstruction) and of pre-gridded NumPy meshes.

Version 0.2 is a substantial reliability overhaul of the original wrapper: fixed shot-noise and log-binning bugs, an OpenMP data race, memory leaks and crash paths; a modern build; a keyword-based API that no longer touches configuration files at all (settings travel as command-line options to the original C parser); and a 103-test pytest suite validating the pipeline against an independent NumPy reference implementation. See REFACTORING.md and CHANGELOG.md.

Contents

• Installation
• Estimators
• API
  • Simulation boxes
  • Boxes with randoms (reconstruction)
  • Surveys / light cones
  • Pre-gridded meshes
  • Cross spectra
• Configuration options
• Results
• Units and conventions
• Legacy configuration-file interface
• Tests and examples
• Threading and performance
• Credits and license

Installation

Requirements:

• a C compiler with OpenMP support (optional but recommended),
• FFTW3 (double precision), e.g.
  • Debian/Ubuntu: sudo apt install libfftw3-dev
  • RHEL/Fedora: sudo dnf install fftw-devel
  • conda: conda install -c conda-forge fftw
  • macOS (Homebrew): brew install fftw libomp
• Python >= 3.9 with NumPy (Cython and a NumPy-2 build environment are pulled automatically by pip).

Then, from a clone of this repository:

pip install .            # regular install
pip install -e .         # editable/development install
# or, for the classic in-tree build:
make inplace             # python setup.py build_ext --inplace

Environment variables understood by setup.py:

variable	effect
POWSPEC_FFTW_DIR (or FFTW_DIR)	prefix of a non-system FFTW (include/, lib/)
POWSPEC_DISABLE_OMP=1	build single-threaded (no OpenMP)
POWSPEC_DISABLE_FFTW_OMP=1	OpenMP in the wrapper, serial FFTW
CC	C compiler used for the feature probes

The OpenMP flavour (-DOMP -fopenmp + libfftw3_omp/libfftw3_threads) is detected with real test compilations and disabled automatically if unavailable. An RPATH is recorded for non-system FFTW prefixes.

Verify the installation:

python -m pytest tests    # 103 tests, ~10 s, no network or data files needed

Estimators

All estimators are FFT-based (mesh assignment -> FFT -> mode counting):

• Periodic boxes (compute_auto_box, compute_cross_box): the classic simulation-box estimator with shot noise V / sum(w) subtracted and multipoles measured with respect to a fixed line of sight (default [0, 0, 1]).
• Boxes with randoms (compute_auto_box_rand, compute_cross_box_rand): the field data - alpha * randoms with alpha = sum(w_data)/sum(w_rand), and the shot noise of both samples, P_shot = V (sum(w_d^2) + alpha^2 sum(w_r^2)) / (sum w_d)^2 (= V/N_d + V/N_r for unit weights). This is the estimator used for reconstructed boxes with shifted randoms.
• Surveys / light cones (compute_auto_lc, compute_cross_lc): FKP-style estimator with completeness and FKP weights and n(z) normalisation, optional internal conversion of [RA, Dec, z] to comoving coordinates (Legendre-Gauss integration of the fiducial cosmology, or a user-provided z -> chi table), and multipoles with a local line of sight (Hand et al. 2017 decomposition).
• Meshes (compute_auto_box_mesh): the box estimator applied to a pre-gridded number-counts mesh.

API

Every function accepts the coordinates/weights as 1-D array-likes (float32 or float64; non-contiguous views such as data[:, 0] are fine), validates the inputs up front, and returns a PowerSpectrumResult.

Simulation boxes

pypowspec.compute_auto_box(x, y, z, w=None, *,
                           powspec_conf_file=None, output_file=None, **conf)

• x, y, z: Cartesian coordinates in [0, box_size) (mandatory conf setting: box_size);
• w: weights (default: ones);
• output_file: also write the powspec ASCII output;
• **conf: any configuration option, e.g. box_size=1000, grid_size=256, bin_size=0.01, poles=(0, 2, 4), assign="TSC", interlace=True, los=(0, 0, 1), verbose=False.

Boxes with randoms (reconstruction)

pypowspec.compute_auto_box_rand(x, y, z, w, rand_x, rand_y, rand_z, rand_w,
                                *, powspec_conf_file=None, output_file=None,
                                **conf)

Surveys / light cones

pypowspec.compute_auto_lc(x, y, z, wcomp, wfkp, nz,
                          rand_x, rand_y, rand_z, rand_wcomp, rand_wfkp,
                          rand_nz, *, powspec_conf_file=None,
                          output_file=None, convert=True, **conf)

• coordinates are [RA, Dec, z] in degrees when convert=True (the default), otherwise pre-converted comoving x, y, z;
• convert=True requires a fiducial cosmology: omega_m (mandatory), optionally omega_l (default 1 - omega_m), eos_w (default -1), cmvdst_err, or a z_cmvdst_file table;
• wcomp/wfkp may be None (ones); nz (radial number density, may be constant) is required for the normalisation;
• box_size is optional - the box is sized automatically from the catalog extent plus box_pad (default 2%). bin_size is mandatory unless box_size is given.

Pre-gridded meshes

pypowspec.compute_auto_box_mesh(mesh, *, powspec_conf_file=None,
                                output_file=None, **conf)

• mesh: (Ng, Ng, Ng) array of number counts (not overdensities), C-contiguity/dtype handled automatically;
• assign must match the scheme used to grid the mesh so that the aliasing correction deconvolves the right window;
• box_size is mandatory; grid_size defaults to Ng; interlacing is disabled (an explicit interlace=True raises).

Cross spectra

pypowspec.compute_cross_box(x1, y1, z1, w1, x2, y2, z2, w2, *,
                            powspec_conf_file=None, output_auto=None,
                            output_cross=None, **conf)
pypowspec.compute_cross_box_rand(..., powspec_conf_file=None,
                                 output_auto=None, output_cross=None, **conf)
pypowspec.compute_cross_lc(..., powspec_conf_file=None, output_auto=None,
                           output_cross=None, convert=True, **conf)

Both auto spectra and the cross spectrum are computed in one call (shared meshes and FFTs). output_auto takes two file names (pass ["", ""] to skip computing the auto spectra - the corresponding entries of the result are then NaN); output_cross takes one.

Configuration options

Keyword arguments are validated by pypowspec.conf.PowspecConf and passed to the original C parser as command-line options (PowspecConf.to_argv)

• no configuration file is created at any point, yet the semantics are exactly those of upstream powspec:

keyword	default	meaning
box_size	-	box side length (scalar or 3-vector). Mandatory for boxes/meshes
box_pad	C default 0.02	padding for automatic box sizing (surveys)
grid_size	256	mesh cells per side (1 < Ng <= 65536)
assign	"TSC"	"NGP", "CIC", "TSC", "PCS" (or 0-3)
interlace	False	interlaced grids (aliasing reduction)
poles	(0, 2, 4)	multipoles, 0 <= ell <= 6 (sorted, deduplicated)
bin_size	2*pi/box_size	k-bin width (or log10 ratio if log_scale)
kmin, kmax	0, Nyquist	bin range
log_scale	False	logarithmic k bins (kmin > 0 required)
los	[0, 0, 1]	line-of-sight unit vector (boxes; normalised for you)
omega_m, omega_l, eos_w	-, 1-omega_m, -1	cosmology for the coordinate conversion
cmvdst_err	C default 1e-8	conversion integration tolerance
z_cmvdst_file	-	z -> chi table replacing the cosmology
output_header, overwrite, verbose	True, 1, False	output-file options

Results

Functions return a PowerSpectrumResult: a dict (backwards compatible) whose entries are also available as attributes.

key	shape/type	meaning
k	(nbin,)	bin centres
kmin, kmax	(nbin,)	bin edges
kavg	(nbin,)	mean wavenumber of the modes in the bin (0 for empty bins)
nmodes	(nbin,) int	modes per bin (half-space doubling as in the C code)
poles	list	multipole orders, matching axis 1 of the multipole arrays
multipoles	(nbin, nl)	auto spectra (shot-noise subtracted)
auto_multipoles	(2, nbin, nl)	auto spectra of both catalogs (cross functions)
cross_multipoles	(nbin, nl)	cross spectrum (no shot-noise subtraction)
n_data_objects, w_data_objects	int/float (or lists)	data catalog size / total weight
n_rand_objects, w_rand_objects, alpha	int/float (or lists)	randoms metadata (when used)
shot_noise, normalisation	float (or lists)	estimator metadata
dk, log_scale	float, bool	binning info

Units and conventions

• Lengths/wavenumbers are in your input units, typically Mpc/h and h/Mpc; nz in (Mpc/h)^-3.
• Coordinates of simulation boxes must lie in [0, box_size).
• The DC mode (k = 0) is counted in nmodes but excluded from power sums (upstream convention); with kmin = 0 the first bin therefore usually reports 0.
• Odd multipoles of periodic boxes vanish by construction (the estimator folds +/-k).
• Survey multipoles use the local line of sight; box multipoles the global los vector.
• Repeated identical calls are bitwise reproducible except for the OpenMP-atomic mesh deposition of weighted catalogs (~1e-15 relative).

Legacy configuration-file interface

The historical interface is preserved:

from pypowspec.pypowspec import compute_auto_box      # deprecated shim
from pypowspec import _pypowspec                      # low-level module

pk = _pypowspec.compute_auto_box(x, y, z, w, "my_powspec.conf",
                                 output_file="out.powspec")

The functions take the same positional arguments as before (except that powspec_conf_file is now optional), return the same dictionary keys, and additionally accept pathlib.Path objects and a conf_options list of --key=value strings that are forwarded to the C parser (taking priority over any file, and sufficient on their own - powspec_conf_file=None is fully supported). When a configuration file is used, its OUTPUT_* and OVERWRITE entries are overridden by the function arguments (command-line priority), so a file prepared for the stand-alone program keeps working. Catalog-related keys (DATA_CATALOG, formatters, selections, ...) are ignored - the arrays you pass in replace them. Reserved options (--conf, --auto, --cross, --help, --template) are rejected in conf_options; the last two would otherwise terminate the process from inside the C code.

Tests and examples

python -m pytest tests            # 103 tests, ~10 s

The suite includes an independent NumPy reference implementation of the box estimator (tests/reference.py: NGP/CIC assignment incl. the C boundary semantics, interlacing, alias corrections, half-space mode counting, multipoles, shot noise), analytic checks (sinusoidal field with exact amplitude, white noise, data == randoms invariants), a scipy-quadrature check of the coordinate conversion, output-file format checks, regression tests for every fixed bug, and robustness tests (memory-leak, thread-safety, GIL release, invalid inputs).

Runnable examples with plots live in examples/:

1. 01_auto_power_spectrum_box.ipynb - white-noise sanity check, clustered mock, assignment schemes and interlacing;
2. 02_cross_spectra_and_randoms.ipynb - cross spectra of biased tracers, reconstruction-style random subtraction;
3. 03_survey_lightcone.ipynb - RA/Dec/z conversion, FKP-style survey estimator;
4. 04_power_spectrum_from_mesh.ipynb - the mesh interface and window consistency.

Threading and performance

• The C engine parallelises mesh assignment, FFTs and mode counting with OpenMP; honour OMP_NUM_THREADS to control it.
• The GIL is released during the computation, and concurrent calls are serialised by an internal lock (FFTW planning is not thread safe), so pypowspec can be used from thread pools.
• Memory: roughly (16-40) * Ng^3 bytes per call depending on options (fields, FFT plans and auxiliary arrays); the auxiliary arrays that the original fork allocated unnecessarily for boxes with randoms (~1.7 GB at Ng = 512) are no longer allocated.

Credits and license

• The C code is powspec by Cheng Zhao (MIT license; some files GPLv3 - see LICENSE.txt and LICENSE_MIT.txt), with the wrapper-specific modifications documented in REFACTORING.md.
• The Python wrapper was originally written by Daniel Forero-Sanchez (https://github.com/dforero0896/pypowspec).