pyfcfc — Python bindings for the Fast Correlation Function Calculator

pyfcfc wraps the two pair-counting components of the C package FCFC (by Cheng Zhao) so that they can be driven entirely from Python with in-memory NumPy catalogues:

module	FCFC component	data	line of sight
pyfcfc.boxes	FCFC_2PT_BOX	periodic-box Cartesian coordinates	box z axis
pyfcfc.sky	FCFC_2PT	Cartesian comoving coordinates, or {RA, Dec, redshift} with on-the-fly cosmological conversion	pair midpoint

Pair counts, correlation functions, multipoles and projected correlation functions are computed with k-d tree or ball tree dual-tree traversal, OpenMP parallelisation and (optionally) SIMD vectorisation, and returned as plain NumPy arrays in a Python dictionary.

The package also ships utilities (pyfcfc.utils) to accumulate split-random measurements, integrate correlation functions, and convert results into pycorr states, so that pair counting can be done with (fast) FCFC while the analysis uses the pycorr API (estimators, rebinning, covariance tools, saving...). The optional pyfcfc.external module converts results into lsstypes containers (Count2, Count2Correlation, Count2CorrelationPoles, ...), giving access to replica covariances and serialisation.

Compared to pycorr (Corrfunc engine), pyfcfc pair counts and correlation functions agree to machine precision on identical inputs (see examples/); the wall-clock speedup is typically 1.5–2.5x on the counting step (and larger with SIMD enabled).

---

Table of contents

• Installation
• Quick start
• The py_compute_cf function
  • Configuration options
  • Soft histogramming (optional method)
  • Result dictionary
• Conventions
  • Pair counts and normalizations
  • Multipole integration: midpoint vs exact
• Working with pycorr
• Accumulating split-random measurements
• Tests
• Examples
• Gotchas and notes
• Changes with respect to the original port
• Acknowledgements

---

Installation

Requirements: a C compiler with OpenMP support (gcc/clang), Python ≥ 3.9, NumPy, SciPy; Cython and NumPy at build time (handled automatically by pip).

pip install .            # from the repository root
# or, for development:
pip install -e .
# or with the Makefile:
make

A one-shot script that installs the system build dependencies, creates a virtual environment with all (optional) extras and installs pyfcfc in editable mode is provided:

bash scripts/setup_dev_env.sh [VENV_DIR]

SIMD acceleration. FCFC's vectorised counting kernels are enabled with

PYFCFC_WITH_SIMD=1 pip install . --no-build-isolation

This implies -march=native: build on a machine whose instruction set (AVX2/AVX512) is available at run time. Compilation of the SIMD kernels is memory hungry (a couple of GB of RAM); PYFCFC_NO_LTO=1 and PYFCFC_EXTRA_CFLAGS="-O2 -mno-avx512f" can reduce the footprint (the latter restricts the build to AVX2). Besides the hard kernels, this also vectorises the soft histogramming kernels (leaf-pair range checks in the box component, libmvec exp in the windowed softmax on glibc Linux, to which -lmvec is then linked automatically); see the soft-histogramming section below for measured speedups and numerical guarantees.

Other build environment variables: PYFCFC_MARCH_NATIVE=1 (native code generation without SIMD), PYFCFC_NO_LTO=1 (disable link-time optimisation), PYFCFC_EXTRA_CFLAGS="..." (extra compiler flags).

Optional dependencies for the extras:

• pycorr (+ Corrfunc, its default engine) for the interoperability utilities: pip install git+https://github.com/cosmodesi/pycorr.git and pip install git+https://github.com/cosmodesi/Corrfunc@desi;
• matplotlib to run the plotting examples.

The number of OpenMP threads is read from OMP_NUM_THREADS at run time.

Quick start

import numpy as np
from pyfcfc.boxes import py_compute_cf

data = ...   # (N_d, 3) float array, periodic-box coordinates
rand = ...   # (N_r, 3) float array
w_d, w_rand = np.ones(len(data)), np.ones(len(rand))

s_edges = np.arange(0, 151, 5.)
results = py_compute_cf(
    [data, rand], [w_d, w_rand],
    s_edges,          # separation bin edges
    None,             # pi bin edges (only for bin=2)
    48,               # number of mu bins in [0, 1) (bin=1)
    label=['D', 'R'],                 # catalogue labels
    bin=1,                            # binning scheme: 0 iso, 1 (s,mu), 2 (s_perp,pi)
    pair=['DD', 'DR', 'RR'],          # pair counts to evaluate
    cf=['(DD - 2 * DR + RR) / RR'],   # correlation function estimator(s)
    multipole=[0, 2, 4],              # Legendre multipoles to integrate
    box=1000.,                        # periodic box side(s)
)
s, xi0, xi2, xi4 = results['s'], *results['multipoles'][0]

Survey-like data, converting {RA, Dec, z} internally:

from pyfcfc.sky import py_compute_cf

results = py_compute_cf(
    [rdd_data, rdd_rand], [w_data, w_rand],
    s_edges, None, 48,
    label=['D', 'R'], bin=1, pair=['DD', 'DR', 'RR'],
    cf=['(DD - 2 * DR + RR) / RR'], multipole=[0, 2, 4],
    convert=True, omega_m=0.31, omega_l=0.69, eos_w=-1.0,
)

Configuration files (conf=...) are still supported (catalogue/position entries are ignored, since the data come from Python); keyword arguments override the file, exactly as command-line options do in FCFC.

The py_compute_cf function

py_compute_cf(data_cats, data_wts, sedges, pedges=None, nmu=1, **kwargs)

• data_cats: sequence of catalogues, each convertible to an (N, 3) array. Non-contiguous arrays, lists and integer/float32 inputs are converted internally (a copy is made; your arrays are never modified). float32 is used natively only when all positions and weights are float32; otherwise everything is promoted to float64.
• data_wts: sequence of per-object weights, each of length N. Pass arrays of ones for unweighted catalogues (this automatically enables FCFC's faster, exact integer counting path).
• sedges: separation (or s_perp) bin edges; strictly increasing, non-negative, finite. For the box component, smax < box/2.
• pedges: pi bin edges, required for bin=2.
• nmu: number of mu bins in [0, 1), used for bin=1.
• **kwargs: FCFC configuration options (below). Python booleans are translated to 'T'/'F', sequences to FCFC array literals ([0, 2, 4]), and unknown options raise ValueError with suggestions instead of being silently ignored.

How configuration parsing works (and what is left of the parameter files)

All options reach FCFC through a single libcfg parameter table, whether they come from keyword arguments or from a configuration file: the keywords are translated to command-line-style --key=value strings (which take priority), and an optional conf= file is read afterwards with lower priority, exactly as in the FCFC CLI. Since the catalogues and the bin edges always come from the function arguments, the upstream file-based inputs were removed from the table (CATALOG, POSITION, ASCII_*, SELECTION, SEP_BIN_*, MU_BIN_NUM, PI_BIN_*), as were the CLI-only --help/--version/--template callbacks and the FITS/HDF5 catalogue readers (excluded from the build).

What remains of the file machinery, and stays supported:

• the conf= overlay itself (any option in the table below);
• the optional outputs pair_output, cf_output, mp_output, wp_output (with out_format and overwrite); when a PAIR_COUNT_FILE already exists, FCFC reads the pair counts from it instead of recomputing them (tested: write/read round-trip);
• cmvdst_file (sky component): a {redshift, comoving distance} table replacing the internal distance integration;
• the verbose=True configuration summary.

Configuration options

option (keyword)	boxes	sky	meaning
label	✔	✔	single-uppercase-letter catalogue labels (defaults to A, B, C, ...)
pair	✔	✔	pair counts to evaluate, e.g. ['DD','DR','RR']
bin	✔	✔	binning scheme: 0 = isotropic, 1 = (s, mu), 2 = (s_perp, pi)
soft_hist	✔	✔	soft histogramming: spread every pair over all separation bins with softmax Gaussian-kernel weights (experimental)
soft_sigma	✔	✔	Gaussian kernel width in the units of sedges (default: half the mean bin width); requires soft_hist
soft_sigma_log	✔	✔	interpret soft_sigma as a dimensionless fraction f with a scale-dependent kernel width sigma(s) = f * s, in the units of sedges (for logarithmic separation bins); requires soft_hist
soft_mu	✔	✔	soft histogramming of the mu dimension of the (s, mu) scheme: softmax Gaussian kernel over the mu bin centres (bin=1 only)
soft_mu_sigma	✔	✔	mu kernel width in mu units (default: half the mu bin width); requires soft_mu
soft_pi	✔	✔	soft histogramming of the pi dimension of the (s_perp, pi) scheme (bin=2), making w_p differentiable
soft_pi_sigma	✔	✔	pi kernel width in the units of pedges (default: half the mean pi bin width); requires soft_pi
soft_pad	✔	✔	edge padding of the soft counting range in kernel widths (default: 2; 0 disables): removes the truncated-kernel-tail bias of the outermost bins, see below
cf	✔	✔	correlation function estimator expression(s), e.g. ['(DD - 2*DR + RR)/RR']; '@@' denotes the analytic random-random counts (box)
multipole	✔	✔	multipole orders to integrate (requires bin=1 and cf)
wp	✔	✔	compute the projected correlation function (requires bin=2 and cf)
box	✔		periodic box side length(s)
convert		✔	convert {RA, Dec, z} to comoving Cartesian coordinates
omega_m, omega_l, eos_w		✔	fiducial cosmology for the conversion
cmvdst_err, cmvdst_file		✔	accuracy / lookup table for the distance integration
data_struct	✔	✔	0 = k-d tree (default), 1 = ball tree
weight	✔	✔	per-catalogue weight flags; set automatically ('1' when all weights are exactly 1, enabling the integer path)
pair_output, cf_output, mp_output, wp_output	✔	✔	optional files where FCFC writes its results
out_format, overwrite	✔	✔	output file options
verbose	✔	✔	print FCFC's configuration summary and progress (default: False)

Errors inside the C code (bad configuration, missing pairs, ...) raise Python exceptions (RuntimeError) and never crash the interpreter; all allocated memory is released on both success and failure paths.

Soft histogramming (optional method)

Besides the standard ("hard") histogram pair counting of FCFC, both modules support an optional soft histogramming method, enabled with the soft_hist=True keyword (or SOFT_HIST = T in a configuration file):

res = py_compute_cf([data, rand], [w_d, w_r], s_edges, None, 48,
                    soft_hist=True, soft_sigma=2.0,   # kernel width, units of sedges
                    label=['D', 'R'], bin=1, pair=['DD', 'DR', 'RR'], box=1000.)

Instead of adding the weight of every pair to exactly one separation bin, the weight is distributed over all separation bins with the softmax of a Gaussian kernel centred on the bins,

x_i = (s - c_i)^2 / (2 sigma^2),    W_i = softmax_i(-x_i),    cnt_i += w_pair * W_i,

where c_i are the bin centres and sigma = soft_sigma (soft_sigma=None / SOFT_SIGMA = 0 selects automatically half the mean separation bin width). Since the softmax weights sum to one, the total weight of all counted pairs is conserved exactly (over the padded counting range when soft_pad > 0; see below): the soft histogram is a kernel-smoothed version of the hard one, and the two are directly comparable (soft counts converge to the hard counts in the sigma -> 0 limit). With soft_sigma_log=True (SOFT_SIGMA_LOG = T) the width is instead scale dependent: soft_sigma is then interpreted as a dimensionless fraction f and the kernel width becomes sigma(s) = f * s (in the units of sedges, evaluated at each pair's separation). This is the natural choice for logarithmic separation bins, where a fixed sigma would wash out the small bins while barely smoothing the large ones: locally the bin width is ~ s · dln s, so the ratio sigma / bin width = f / dln s is the same at every scale. The default f (soft_sigma=None) is half the mean logarithmic bin width, i.e. sigma equals half the local bin width everywhere. The mu and pi kernels always take fixed absolute widths (soft_mu_sigma in mu units, soft_pi_sigma in the units of pedges): FCFC's mu bins are uniform on [0, 1) by construction, and no log scaling is offered for pi.

Edge leakage and the soft_pad padding. The Gaussian kernel has infinite support, but the tree traversal only counts pairs inside the binning range [sedges[0], sedges[-1]): the kernel tails reaching beyond the outermost edges are truncated. For the outermost bins the estimator then averages xi(s) over a half kernel -- the weight that leaks out of the last bin towards larger separations (where xi is smaller for a steep correlation function) is not compensated by weight leaking in from beyond sedges[-1], because those pairs are never counted. The result is a positive bias of the last separation bin (a smaller negative one at the first, where the s^2 volume factor suppresses the missing pairs), growing with the steepness of xi and with the kernel width: for a xi ~ s^-1.8 power law with the default scale-free kernel the last bin carries ~ +9% excess (relative to the kernel-smoothed value the interior bins measure), ~ +5% for gamma = 0.8 and ~ +12% for gamma = 2.6. soft_pad = z (the default, 2 kernel widths) removes it: the counting range is extended by z kernel widths beyond the outermost edges, realised with buffer bins continuing the separation grid (geometrically for log-spaced bins, arithmetically otherwise). The softmax then normalises over the real and buffer centres, pairs beyond the original edges contribute their decaying kernel tails to the outermost real bins, and the buffer-bin counts are discarded -- the returned results refer to the original bins, with their original edges. soft_pad=0 disables the padding and restores the unpadded estimator exactly. Details:

• the padding is quantised to whole buffer bins per side, ceil(z * sigma_edge / spacing): with the default kernel (f = 0.5 * dln s) any 0 < z <= 2 gives one buffer bin per side (z_eff = 2);
• measured wall-clock overhead of the default on the soft counts (pairs scale with the shell volume, ~ exp(3 z f) - 1 for log-spaced bins): +32% at 0.1 dex bins, +69% at 0.2 dex, +20% for linear bins with the default kernel, +29% for a survey-like catalogue; soft_pad=0 restores the previous speed exactly;
• in a periodic box the extended maximum edge is clamped below L/2 (FCFC's minimum-image limit), so setups binning close to half the box get reduced (possibly no) upper padding;
• the lower edge is never padded below zero;
• soft_pad cannot be combined with a configuration file (conf=), whose box size and kernel options the Python wrapper cannot see;
• with padding, the weight conservation and the sigma -> 0 limit hold on the extended grid (the real bins alone keep less than the total, the rest sits in the discarded buffer bins).

Properties and limitations:

• the mu (bin=1) and pi (bin=2) dimensions remain hard-binned unless soft_mu=True / soft_pi=True are set: soft binning them makes the 2D counts smooth in both arguments, hence the multipoles xi_ell(s) and the projected correlation function w_p (linear operators on xi(s, mu) / xi(s_perp, pi)) become differentiable as well;
• pair counts become fractional and are always accumulated in double precision (cf_setup forces the weighted counting path in soft mode);
• the tree traversal, periodic wrapping and per-pair range checks replicate the scalar hard kernels exactly, so both methods operate on the same set of pairs;
• the soft kernels are OpenMP-parallel with the same thread decomposition as the hard kernels; for linear separation bins the kernel weights are generated with an exact multiplicative recurrence (two exp() calls per pair) over an adaptive window, so the per-pair cost is O(min(ns, sigma / ds)) instead of O(ns);
• PYFCFC_WITH_SIMD=1 builds accelerate the soft method: in the box component the leaf-pair loops evaluate the coordinate differences, periodic wraps and range checks of 4 (AVX2) or 8 (AVX-512) pairs at a time, with vector arithmetic that matches the scalar expressions bit-for-bit (no FMA contraction), so exactly the same pairs are binned; in both components the windowed stable softmax (non-linear bins and soft_sigma_log kernels) evaluates its Gaussian terms with glibc's libmvec vector exponential on x86-64 Linux (scalar exp elsewhere). Measured single-thread gains of such a build (2.5 GHz Xeon, AVX-512): -20...-30% for box s / sigma(s)=f*s binning, -11% for box (s, mu), -12% for box (s_perp, pi), -8...-12% for the survey-like component; the survey-like leaf loops stay scalar on purpose (there the pairs that reach the leaves mostly pass the range checks, so a vector pre-filter does not pay off). Counts differ from non-SIMD builds only by floating-point summation order and the ulp-level accuracy of the vector exponential (~1e-15 relative). The adjoint (VJP) traversals remain scalar but share the vectorised softmax, so the backward pass still replays the forward weights exactly. Single-threaded the soft method is typically a few times slower than the hard kernels (dominated by the kernel evaluation), and it parallelises like them; MPI builds are not supported.

Validation against O(N^2) brute-force soft references, the conservation laws and the sigma -> 0 limit are covered by tests/test_soft_hist.py (unpadded) and tests/test_soft_pad.py (padding plans, padded counts vs brute force on the extended grid, adjoints, clamping, validation); examples/example_soft_hist.py plots the soft histograms against the hard ones (counts, ratios and xi(s)).

Adjoint (reverse-mode derivatives). Because the soft histogram is a smooth function of the pair separations, both modules also expose a VJP (vector-Jacobian product) entry point,

out = soft_hist_vjp(data_cats, data_wts, sedges, pedges, nmu,
                    pair='DD', cotangent=g, soft_sigma=2.0,
                    grad_wts=False, grad_sigma=False, box=1000.)
# out['grads']     : dict label -> (N, 3) coordinate cotangents
# out['grad_wts']  : dict label -> (N,) weight cotangents (or None)
# out['grad_sigma']: kernel-width cotangent (or None)

which propagates cotangents from the raw pair count bins back to the catalogue coordinates (and, optionally, weights and kernel width), i.e. exactly what a reverse-mode autodiff engine needs to backpropagate through the pair counting step. With W_i the softmax weights of a pair, cbar = sum_i W_i c_i and bin cotangents g_i,

dW_i/ds = W_i (c_i - cbar) / sigma^2,
s_bar   = (w / sigma^2) * sum_i g_i W_i (c_i - cbar),

chained to the coordinates through s (or s_perp for bin=2; for the survey-like component through the midpoint-LOS decomposition ds_perp/dr1 = (1 - q)(d - q m)/s_perp with q = d.m/|m|^2). Discrete choices (range checks, mu / pi bins, kernel window) are held constant, as in any autodiff engine. The adjoint replays the forward pair selection with its own index-preserving k-d tree (FCFC's counting trees store permuted coordinates without an index map), is OpenMP-parallel, and is validated against central finite differences in tests/test_soft_vjp.py (including the weight and sigma cotangents, and the invariant that a constant cotangent yields zero coordinate gradient). examples/example_soft_hist_jax.py wraps the forward and the adjoint into a jax.custom_vjp primitive and demonstrates it: gradient validation against finite differences, and a differentiable Alcock-Paczyński-like stretch recovery with jax.grad + Adam (figures/soft_hist_jax.png). examples/example_soft_hist_multipoles.py registers the soft 2D counts as a jax.custom_vjp primitive and differentiates the multipoles xi_0, xi_2 (Landy-Szalay + FCFC's Legendre integration, in JAX), recovering Alcock-Paczyński stretch parameters from them (figures/soft_hist_multipoles.png). examples/example_soft_hist_wp.py makes the projected correlation function differentiable (soft s and pi bins) and recovers a transverse stretch from w_p alone (figures/soft_hist_wp.png). examples/example_soft_hist_rppi.py treats the full 2D xi(r_p, pi) map as the differentiable object: it computes the adjoint Jacobian maps d xi(r_p, pi)/d theta of an anisotropic stretch model with one-hot cotangents and validates them bin-by-bin against finite differences (figures/soft_hist_rppi.png). examples/example_soft_hist_fit.py shows the same machinery without any autodiff framework: a 3-parameter rigid displacement of the galaxy catalogue (one shift per axis) is recovered by Adam on the loss 1/2 |DR(theta) - DR(theta_*)|^2, with dL/dtheta obtained by contracting the soft_hist_vjp coordinate cotangents (gradient check vs finite differences and the rigid-translation null direction included; figures/soft_hist_fit.png).

JAX integration. The optional pyfcfc.jax module (requires jax) packages the forward + adjoint into a ready-made autodiff primitive: soft_pair_counts(positions, weights, module=..., s_edges=..., ...) is a jax.custom_vjp function returning the raw pair counts as JAX arrays, differentiable with respect to the catalogue positions and weights (its forward calls py_compute_cf, its backward soft_hist_vjp). The primitive is created once per static configuration and cached, so repeated calls - also inside jax.jit-compiled outer functions - reuse the same callable and hence JAX's tracing/compilation caches; the C library runs on host arrays registered via jax.pure_callback, so the primitive works eagerly and under jax.jit, for both the forward pass and jax.grad (jit correctness and compilation caching are covered by tests/test_soft_vjp_jax.py). Importing pyfcfc.jax does not enable JAX's x64 mode: the C library accumulates in double precision internally, but the module follows the precision of the JAX inputs, returning float32 counts and float32 gradients for float32 catalogues (float64 inputs use the full double-precision path when x64 is enabled). The helpers pair_normalizations, xi_landy_szalay, multipoles, projected_wp and legpoly implement FCFC's conventions in JAX, so a whole xi / multipole / w_p pipeline becomes differentiable with plain jax.grad. examples/example_soft_hist_wp.py builds a differentiable w_p on realistic logarithmic r_p bins this way (scale-free kernel sigma(s) = f * s via soft_sigma_log, transverse stretch recovered with Adam; figures/soft_hist_wp.png, and a kernel-width comparison in figures/soft_hist_wp_sigmas.png).

JIT caching and recompilation. Under jax.jit, only the first call for a given input signature (array shapes/dtypes plus the static configuration) pays the tracing + XLA compilation cost; every subsequent call reuses the cached executable - there is no recompilation. examples/benchmark_jit_cache.py measures exactly this (DD counts, 20 log-spaced s bins 0.5-40 Mpc/h, scale-free kernel sigma(s) = f * s; recorded on a 2-vCPU Linux box, jax 0.10.2, 2 OpenMP threads - absolute times are machine-dependent, the pattern is not):

N	jitted call	first run	all subsequent runs	speedup
400	forward	32 ms	0.9 ms each	~35x
400	jax.grad	24 ms	1.0 ms each	~24x
4000	forward	~30-60 ms	~10 ms each	~3-6x

Per-call series of the jitted forward at N = 400 (ms): [32.2, 0.9, 0.9, 0.9, 0.9, 0.8, 0.9, 0.9, 0.9, 0.9] - only the first call compiles. At N = 400 a warm jitted call is also ~17x cheaper than an eager (non-jit) one, because jax.jit additionally amortises the Python-side argument packing of soft_pair_counts; at N = 4000 the C pair counting itself (~10 ms here) dominates every call. A new input shape (e.g. a different number of points) triggers exactly one extra compilation; repeated calls with it reuse that one too. The tests assert both the compilation-cache sizes (where the installed jax exposes fn._cache_size()) and a first-vs-subsequent timing check that works on any jax version (tests/test_soft_vjp_jax.py). The module adapts to the jax.pure_callback API across jax versions (vmap_method on jax >= 0.4.31, vectorized on older releases); run python examples/benchmark_jit_cache.py (optionally --json FILE) to reproduce the numbers on your machine.

Result dictionary

key	shape / content
labels	list of catalogue labels
number	dict label → number of objects
weighted_number	dict label → sum of weights (equals number for unit weights)
normalization	dict pair → pair-count normalization (see conventions)
pairs	dict pair → normalized pair counts, shape (ns,), (ns, nmu) or (ns, np); plus bin edge arrays smin, smax and (mumin, mumax) or (pimin, pimax)
s	(ns,) centers of the separation (s_perp) bins
cf	(ncf, ns[, nmu or np]) estimator values, if cf was given
multipoles	(ncf, nl, ns) multipole moments, if multipole was given
poles	list of the multipole orders computed
projected	(ncf, ns) projected correlation function w_p, if wp=True

The pair-count arrays are normalized counts; multiply by normalization[pair] to recover the raw (weighted) pair counts.

Conventions

Pair counts are always computed

Every pair listed in pair (or PAIR_COUNT in a configuration file) is evaluated and returned in results['pairs'], whether or not a CF estimator references it: the pair counts are the main product of the code, and estimators/multipoles/projections are optional extras computed from them (here, or externally with pyfcfc.utils, pycorr or lsstypes). Passing a dummy estimator such as cf=['DD'] is therefore not needed (for historical reasons the original port required an estimator to keep the requested binning scheme; this was fixed in 0.2.0 and is covered by regression tests). The only pairs not counted from scratch are those with an existing PAIR_COUNT_FILE output, which FCFC reads instead.

Pair counts and normalizations

• Auto pairs (e.g. DD) are counted twice per unordered pair (ordered-pair convention), and the normalization is (Σw)² − Σw² — the same convention as pycorr's wnorm.
• Cross pairs (e.g. DR) are counted once per (i, j) combination, with normalization Σw₁ × Σw₂.
• For bin=1 and bin=2, only mu ≥ 0 (pi ≥ 0) is binned, with mu (pi) the absolute value of the cosine (parallel separation).
• Bins are left-closed, right-open; pairs with exactly mu == 1 are dropped (as in FCFC compiled without WITH_MU_ONE).
• The box component uses the periodic minimum-image metric; the sky component defines the line of sight as the pair midpoint direction (equivalent to pycorr's los='midpoint').

Multipole integration: midpoint vs exact

FCFC integrates multipoles with the midpoint rule, ξ_ℓ(s) = (2ℓ+1) Σ_j ξ(s, μ_j) P_ℓ(μ_j) Δμ. pycorr instead uses the exact integral of the Legendre polynomial over each bin. The two agree up to O(Δμ²): for ℓ = 0 they coincide (the difference cancels), while for ℓ = 2, 4 the difference can reach ~10⁻² at small scales when ξ is large — this is a scheme difference, not a counting error (the raw pair counts of both codes agree to machine precision).

pyfcfc.utils.compute_multipoles supports both schemes:

from pyfcfc.utils import compute_multipoles
mp_mid   = compute_multipoles(xi_smu, [0, 2, 4])                    # FCFC's scheme
mp_exact = compute_multipoles(xi_smu, [0, 2, 4], method='exact')    # pycorr's scheme

so you can reproduce either code's numbers exactly from pyfcfc pair counts (both also accept ignore_nan=True to skip empty bins, as pycorr does).

Working with pycorr

from pyfcfc.utils import pairs_to_pycorr
from pycorr import TwoPointCorrelationFunction

state = pairs_to_pycorr(
    results,                                  # pyfcfc result dict
    'landyszalay',                            # pycorr estimator name
    dict(DD='D1D2', DR=('D1R2', 'R1D2'), RR='R1R2'),   # label mapping
    box_size=1000.,                           # periodic box (optional)
)
np.save('counts.pkl.npy', state)
tpcf = TwoPointCorrelationFunction.load('counts.pkl.npy')
s, xi0 = tpcf(ell=0, return_sep=True)
rebinned = tpcf[::2, ::4]                     # pycorr rebinning works

Notes:

• Reversible pair counts must be mapped to both pycorr names (as for DR above); pyfcfc's |mu| counts are mirrored to pycorr's full [-1, 1] range (exactly for auto pairs, statistically for cross pairs).
• For the natural estimator in a periodic box, the analytic random-random counts are generated automatically when box_size is given.
• The pycorr state records los_type='z' for box measurements and 'midpoint' for sky measurements, matching FCFC's conventions.
• In rppi mode, pycorr stores symmetric pi edges [-pimax, pimax]; the conversion mirrors pyfcfc's pi ≥ 0 counts accordingly.

Working with lsstypes (optional)

lsstypes provides container types for LSS measurements (replica algebra, covariances, serialisation). The optional module pyfcfc.external converts pyfcfc results into them:

pip install git+https://github.com/adematti/lsstypes   # optional extra

from pyfcfc.external import to_lsstypes, to_lsstypes_correlation, to_lsstypes_counts

counts = to_lsstypes_counts(results)              # {'DD': Count2, 'DR': ..., 'RD': ..., 'RR': ...}
corr = to_lsstypes_correlation(results, 'landyszalay')   # Count2Correlation
poles = to_lsstypes(results, ells=[0, 2, 4])      # Count2CorrelationPoles
wp = to_lsstypes(results_spi, project='wp')       # Count2CorrelationWp

• project accepts 'poles' (default), 'wedges', 'wp', 'binned', or None (returns the un-projected Count2Correlation);
• multipoles use lsstypes' exact per-bin Legendre integration, i.e. compute_multipoles(..., method='exact') on the pyfcfc side;
• for the 'natural' estimator in a periodic box without RR counts, the analytic random-random counts are synthesized (same convention as FCFC's @@);
• catalogue labels other than D/R are handled with pair_mapping, e.g. {'AA': 'DD', 'AB': ('DR', 'RD'), 'BB': 'RR'};
• replica tools work out of the box: lsstypes.mean/cov/sum over converted objects, and lsstypes.write/read (.pkl, .h5, .txt). Note that lsstypes.sum weights replicas by the DD normalization, while pyfcfc.utils.add_pair_counts accumulates each pair type with its own normalization (exact raw-count addition); the two coincide for unweighted, equal-size replicas.

See notebooks/04_lsstypes.ipynb for a complete walkthrough and tests/test_lsstypes.py (11 tests, skipped when lsstypes is absent).

Accumulating split-random measurements

from pyfcfc.utils import add_pair_counts

total = None
for chunk in random_chunks:
    res = py_compute_cf([data, chunk.pos], [w_data, chunk.w], ...)
    total = res if total is None else add_pair_counts(total, res)

add_pair_counts combines the raw pair counts exactly (it is a normalization-weighted average of the normalized counts, with summed normalizations), returns a new dictionary, and drops stale estimator-dependent keys (cf, multipoles, ...), which must be recomputed from the combined counts.

Tests

pip install pytest
python -m pytest tests -v          # 61 tests

The suite covers, among other things:

• exact agreement with brute-force O(N²) counting (periodic box and survey-like data; iso/(s, mu)/(s_perp, pi); linear, logarithmic and irregular bins; unit and random weights; k-d and ball trees);
• cross-validation against pycorr/Corrfunc on identical catalogues (pair counts, natural and Landy-Szalay estimators, multipoles, w_p, isotropic mode, state conversion + save/load + rebinning);
• soft histogramming (optional method): agreement with the brute-force O(N^2) soft references (box and sky, all schemes, both trees), exact conservation of the total pair weight with respect to the hard method, and the sigma -> 0 limit;
• API robustness: input validation and conversion, error propagation, immutability of caller arrays, determinism, dtype handling;
• runtime hygiene: GIL release during counting and absence of gross memory leaks on success and failure paths.

Tests that need pycorr/Corrfunc are skipped automatically when those are not installed.

Examples and notebooks

examples/README.md contains the guided tour of the examples, a feature-coverage matrix and the benchmarking guide. In short:

script	content	figure
examples/example_features_tour.py	start here: every feature in one commented script (result dict, binning schemes, weights, conf files, trees, dtypes, cross-pairs, integrations, sky, errors, verbose)	—
examples/example_validation_bruteforce.py	pyfcfc vs brute-force counting, all binning schemes	figures/validation_bruteforce.png
examples/example_box_vs_pycorr.py	periodic box: xi(s), xi_0/2/4, the 2D xi(s_perp, pi), w_p, natural & LS estimators, pyfcfc vs pycorr + timings	figures/box_vs_pycorr.png, figures/box_vs_pycorr_natural.png, figures/box_vs_pycorr_2d.png
examples/example_survey_vs_pycorr.py	survey-like data from {RA, Dec, z}: multipoles and w_p vs pycorr (midpoint LOS)	figures/survey_vs_pycorr.png
examples/example_pycorr_interop.py	split-random accumulation → pycorr state → save/load/rebin	figures/pycorr_interop.png
examples/example_benchmark.py	performance study vs pycorr: sizes, binning schemes, thread scaling, SIMD reporting	figures/benchmark.png, figures/benchmark_results.csv

Run them all with make examples (or individually with python examples/example_*.py; matplotlib and, for the pycorr ones, pycorr + Corrfunc are required).

Interactive versions are provided as executed Jupyter notebooks in notebooks/:

notebook	content
01_quickstart_box.ipynb	result dictionary, xi(s), multipoles, w_p, conventions
02_survey_and_pycorr.ipynb	survey-like data, midpoint LOS, pycorr cross-check, state save/load/rebin
03_benchmarking.ipynb	how to benchmark properly (sizes, threads, SIMD) with a small in-notebook run
04_lsstypes.ipynb	conversion to lsstypes containers: counts, estimators, poles, replica covariance, serialisation

Typical output of the box example (40k data / 80k randoms, 2 threads):

max |pyfcfc - pycorr|:
  xi(s)  : 6.5e-13
  multipoles, exact per-bin Legendre integration: xi_0 1.8e-13, xi_2 9.1e-13, xi_4 1.2e-12
  wp     : 1.6e-11

Gotchas and notes

• multipole (and wp) require a cf estimator: multipoles are integrated from the estimator's xi(s, mu), as in FCFC. If you only need pair counts but still want the correct binning, pass a trivial estimator such as cf=['DD'].
• For bin=1, pass nmu ≥ 2; pycorr comparisons need linear mu bins. nmu is ignored (clamped to 1) for bin=0 and bin=2, and pedges is ignored unless bin=2.
• The s-mode and multipole results are deterministic for unit weights (integer counting); with non-trivial weights, repeated runs can differ at the ~1e-15 relative level due to OpenMP reduction ordering.
• Empty (s, mu) cells produce NaN in cf when the random-random counts vanish there; use ignore_nan=True in the integration helpers or mask them.
• Catalogue labels must be unique uppercase letters (at most 26 catalogues).
• The vendored FCFC-main/etc/*.conf files are the upstream CLI templates: several of their entries (CATALOG, POSITION, ...) do not apply to library mode (the data come from Python), and their CATALOG_LABEL is empty, so pass label= (and pair=, bin=, ...) explicitly or copy the template and fill them in. Unknown entries in a configuration file are reported by libcfg and otherwise ignored.

Repository layout

path	content
pyfcfc/	the Python package: boxes and sky Cython extensions, utils, optional external
FCFC-main/	vendored FCFC C sources (upstream v1.0.1 with selected v1.1.0 fixes backported and the Python-binding modifications described below), plus the new pyfcfc_helpers.{c,h} accessor layer
tests/	pytest suite (brute-force and pycorr cross-validation, API robustness, lsstypes)
examples/	runnable scripts + figures/ outputs (see examples/README.md)
notebooks/	executed Jupyter notebooks
scripts/	setup_dev_env.sh environment bootstrap
setup.py, pyproject.toml, Makefile	build & packaging

Changes with respect to the original port

Version 0.2.0 is a substantial hardening of the original wrapper. The most important fixes:

Correctness bugs

1. Input arrays were mutated in place. cf_setup rescaled the separation/pi bin edges inside the caller's NumPy buffers (the C code aliased them), silently corrupting any array reused across calls. Bins are now deep-copied on entry.
2. Uninitialized variable in the auto-pair normalization. The hack that switched to pycorr's (Σw)² − Σw² normalization read an uninitialized accumulator (wt2sum), i.e. undefined behaviour. Replaced by the clean upstream FCFC v1.1.0 implementation (a w2 field computed when the trees are built).
3. weighted_number reported the number of objects, not the sum of weights, breaking normalizations and the pycorr state conversion for weighted catalogues.
4. w_p used rescaled pi bins (upstream FCFC bug fixed in v1.1.0, backported): dpi is now taken from the unrescaled bins.
5. The multipole orders (poles) and catalogue labels were read from freed memory (the configuration structure is destroyed before the results are extracted on the Python side). Both are deep-copied now.
6. Dangling pointers in the command-line argument array: the bytes objects backing argv were freed when the helper function returned, before the C parser ran. They are now kept alive for the duration of the call (and argv is freed afterwards).
7. Memory leaks / ownership hazards: the initial cf_init() result was leaked on every call; input catalogues leaked on configuration errors; the survey-like weight-buffer branch could leak/overwrite buffers. Ownership of the input data is now explicit and covered by tests.
8. Linear-bin detection compared only the first and last bin widths with == (and mutated the configuration structure); it is now a tolerant check over all widths, and nonlinear bins with coincidental first/last widths can no longer be misclassified.
9. Backported upstream v1.1.0 tree fix (gather_nodes dropping shallow leaves), which can lose pairs in some tree shapes.
10. SIMD padding: with SIMD enabled, FCFC's kernels read a few elements past the end of the coordinate/weight arrays; the Python-provided buffers were unpadded (potential out-of-bounds reads). They are now padded exactly like FCFC's own input routines.
11. pairs_to_pycorr: the rppi conversion used one-sided pi edges without the mirror+halve transformation (wrong normalization in pycorr), the isotropic mode halved counts that should not be halved, iso states crashed on 1-D bin arrays, and los_type was recorded as firstpoint although FCFC uses the midpoint convention. All fixed, and validated against native pycorr counts.

Robustness and usability

• The Cython layer no longer duplicates the C structure layouts (a silent-corruption hazard): all access goes through an accessor layer (pyfcfc_helpers.c) that treats CF/DATA as opaque.
• Full input validation with informative Python exceptions: shapes, dtypes, monotonic/finite/non-negative bin edges, label consistency, unknown configuration keywords (with typo suggestions).
• Python booleans, lists, tuples and NumPy arrays are accepted as keyword values; weight is set automatically so that unweighted catalogues use the exact integer counting path.
• The GIL is released during pair counting; the library is quiet by default (verbose=True restores FCFC's progress output).
• Result arrays are plain contiguous NumPy arrays with deterministic shapes (cf is no longer squeezed), and poles is returned.
• Modern packaging: pyproject.toml, no numpy.compat/distutils imports (works with NumPy 2 and Cython 3), no hardcoded system paths, optional SIMD/march/LTO flags, stale generated C sources removed.
• Dead code from the CLI era was pruned: read_bins (file-based bin edges), the usage/version/conf_template callbacks with their stale texts, the commented-out conf_print blocks, and the FITS/HDF5 catalogue readers (excluded from the build); read_ascii.c stays because it backs cmvdst_file.
• Pair counts read back from an existing PAIR_COUNT_FILE now get their catalogue indices resolved like computed ones; previously their labels were left unset (out-of-bounds label reads in the original port), so the results dictionary could not name them correctly.
• The number of mu bins is now forced to 1 outside the (s, mu) scheme (upstream enforced this in the configuration parser, which the Python binding bypasses): passing nmu > 1 with bin=0 used to overflow the analytic-RR buffer whenever an estimator with @@ was requested, corrupting memory; pi edges are likewise ignored unless bin=2.
• Pair counts requested via pair are always evaluated and returned, even when no CF estimator references them (the estimator only controls the optional cf/multipoles/projected products); a defensive check now raises if FCFC ever returns an unevaluated pair.
• Repository hygiene: the legacy test/ scripts and artefacts of the original port (NERSC-dependent scripts, stale pair-count/.pkl.npy/ .png outputs, command_line_help.txt, duplicated *_original.conf files) and the cython_debug/ and main.zip leftovers were removed; their functionality is superseded by tests/, examples/ and notebooks/.

Acknowledgements

FCFC is written by Cheng Zhao and distributed under the MIT licence (see LICENSE.txt); this wrapper inherits it. If you use FCFC in research, please cite the FCFC paper (Zhao et al.). pycorr is developed by the DESI collaboration (cosmodesi/pycorr).