Two-point correlation functions

Stable estimators

Angular two-point correlation function w(θ).

Provides: - treecorr-based measurement via w_theta() - JAX-native Landy-Szalay and Davis-Peebles kernels operating on pair counts

sum_stat.twopcf.angular.w_theta(gal, rand, theta_bins, estimator='landy-szalay', sep_units='arcmin', n_threads=4)[source]

Angular two-point correlation function w(θ).

Pair counting is performed by treecorr; the estimator arithmetic is done in JAX so it is differentiable with respect to the pair counts.

Parameters:

  • gal (GalaxyCatalogue) – Galaxy catalogue (ra, dec required).

  • rand (GalaxyCatalogue) – Random catalogue (ra, dec required).

  • theta_bins (ndarray, shape (n_bins+1,)) – Angular separation bin edges in units of sep_units.

  • estimator (str) – "landy-szalay" (default) or "davis-peebles".

  • sep_units (str) – Angular units for theta_bins: "arcmin", "deg", or "radians".

  • n_threads (int) – Number of OpenMP threads for treecorr.

Returns:

  • theta_centres (ndarray, shape (n_bins,)) – Geometric mean separation per bin [same units as theta_bins].

  • w (ndarray, shape (n_bins,)) – Angular correlation function (dimensionless).

  • var_w (ndarray, shape (n_bins,)) – Variance estimate from treecorr (Poisson / shot noise).

  • Performance

  • ———–

  • ~4 s/call (N_gal=5000, N_rand=50000, n_bins=20, treecorr, 4 threads, CPU x86-64)

Parameters:
Return type:

tuple[ndarray, ndarray, ndarray]

sum_stat.twopcf.angular.w_theta_from_pair_counts(dd, dr, rr, n_gal, n_rand, estimator='landy-szalay')[source]

Compute w(θ) from pre-computed pair counts using the JAX estimator.

Useful for jackknife/bootstrap resampling where pair counts per region are accumulated and then combined.

Parameters:

  • dd, dr, rr (ndarray, shape (n_bins,)) – Pair count arrays.

  • n_gal, n_rand (int) – Number of galaxies and randoms.

  • estimator (str) – "landy-szalay" or "davis-peebles".

Returns:

w (ndarray, shape (n_bins,))

Parameters:
Return type:

ndarray

Projected two-point correlation function w_p(r_p).

Provides: - corrfunc-based measurement via wp() (Landy-Szalay, pair counting)

sum_stat.twopcf.projected.wp(gal, rand, rp_bins, pi_max, cosmo, n_threads=4)[source]

Projected correlation function w_p(r_p) via Landy-Szalay estimator.

Uses Corrfunc DDrppi_mocks with pre-computed comoving distances to count pairs in (r_p, π) bins, then integrates over π up to pi_max.

Notes

The previous treecorr-based implementation using FisherRperp with min_rpar=0 produced all-negative wp due to a signed-rpar bug: treecorr uses r_par = r2 - r1 (by catalog index, not distance), so min_rpar=0 discards the majority of auto-correlation (DD, RR) pairs while only halving cross-correlation (DR) pairs. Corrfunc with is_comoving_dist=True is the recommended approach for wp on mock catalogues.

Parameters:

  • gal (GalaxyCatalogue) – Galaxy catalogue (ra, dec, redshift required).

  • rand (GalaxyCatalogue) – Random catalogue (ra, dec, redshift required).

  • rp_bins (ndarray, shape (n_bins+1,)) – Projected separation bin edges [Mpc].

  • pi_max (float) – Maximum line-of-sight integration limit [Mpc]. Rounded to the nearest integer because corrfunc uses dπ = 1 Mpc bins.

  • cosmo (FlatLambdaCDM) – Cosmology for redshift-to-comoving-distance conversion.

  • n_threads (int) – Number of OpenMP threads for corrfunc.

Returns:

  • rp_centres (ndarray, shape (n_bins,)) – Geometric mean projected separation per bin [Mpc].

  • wp_vals (ndarray, shape (n_bins,)) – Projected correlation function [Mpc].

  • var_wp (ndarray, shape (n_bins,)) – Poisson variance estimate from DD pair counts [Mpc^2]. Bins with zero DD pairs are assigned np.inf.

  • Performance

  • ———–

  • ~3 s/call (N_gal=5000, N_rand=50000, n_bins=20, pi_max=100, 4 threads, CPU x86-64)

Parameters:
Return type:

tuple[ndarray, ndarray, ndarray]

Development estimators

Warning

The estimators below (3D ξ(r) and multipoles ξ:sub:`ℓ`(s)) are in active development. Their results have not been validated against reference data and the APIs may change without notice. Do not use for scientific analysis.

3D two-point correlation function ξ(r).

Provides: - treecorr-based measurement via xi_r()

sum_stat.twopcf.spatial.xi_r(gal, rand, r_bins, cosmo, n_threads=4)[source]

3D two-point correlation function ξ(r) via Landy-Szalay estimator.

Comoving distances are computed from redshifts using the provided cosmology. treecorr operates in 3D Euclidean metric.

Parameters:

  • gal (GalaxyCatalogue) – Galaxy catalogue (ra, dec, redshift required).

  • rand (GalaxyCatalogue) – Random catalogue (ra, dec, redshift required).

  • r_bins (ndarray, shape (n_bins+1,)) – Comoving separation bin edges [Mpc].

  • cosmo (FlatLambdaCDM) – Cosmology for redshift-to-distance conversion.

  • n_threads (int) – Number of OpenMP threads for treecorr.

Returns:

  • r_centres (ndarray, shape (n_bins,)) – Mean comoving separation per bin [Mpc].

  • xi (ndarray, shape (n_bins,)) – 3D correlation function (dimensionless).

  • var_xi (ndarray, shape (n_bins,)) – Variance estimate (Poisson / shot noise).

  • Performance

  • ———–

  • ~8 s/call (N_gal=5000, N_rand=50000, n_bins=20, treecorr, 4 threads, CPU x86-64)

Parameters:
Return type:

tuple[ndarray, ndarray, ndarray]

Redshift-space multipoles of the two-point correlation function.

Computes ξ_ℓ(s) for ℓ = 0 (monopole), 2 (quadrupole), 4 (hexadecapole).

Provides: - treecorr-based measurement via xi_multipoles() - JAX Legendre decomposition kernel via _legendre_decompose_jax()

sum_stat.twopcf.multipoles.xi_multipoles(gal, rand, s_bins, cosmo, ells=(0, 2, 4), n_mu_bins=100, n_pi_bins=50, n_threads=4)[source]

Redshift-space multipoles ξ_ℓ(s).

Pair counts are accumulated on a 2D (rp, π) grid using treecorr’s FisherRperp metric (one call per π slice), then rebinned to (s, μ) and decomposed into Legendre multipoles by _legendre_decompose_jax().

Notes

  1. Loop over n_pi_bins linear π slices from 0 to s_bins[-1].

  2. For each slice run a separate NNCorrelation with FisherRperp to obtain ξ(rp, π) in the requested rp bins.

  3. Compute s = √(rp² + π²) and μ = π/s for each cell centre and redistribute ξ into the (n_s, n_mu_bins) (s, μ) grid by weighted averaging.

  4. Apply Legendre decomposition.

Parameters:

  • gal (GalaxyCatalogue) – Galaxy catalogue (ra, dec, redshift required).

  • rand (GalaxyCatalogue) – Random catalogue (ra, dec, redshift required).

  • s_bins (ndarray, shape (n_bins+1,)) – Redshift-space separation bin edges [Mpc].

  • cosmo (FlatLambdaCDM) – Cosmology for distance conversion.

  • ells (tuple of int) – Multipole orders, default (0, 2, 4).

  • n_mu_bins (int) – Number of μ bins in the intermediate (s, μ) grid (default 100).

  • n_pi_bins (int) – Number of linear π slices from 0 to s_max (default 50).

  • n_threads (int) – Number of OpenMP threads for treecorr.

Returns:

  • s_centres (ndarray, shape (n_bins,)) – Geometric mean separation per bin [Mpc].

  • xi_dict (dict[int, ndarray]) – Keys are ell, values shape (n_bins,).

  • Performance

  • ———–

  • ~60 s/call (N_gal=5000, N_rand=50000, n_s=20, n_pi=50, CPU x86-64)

Parameters:
Return type:

tuple[ndarray, dict[int, ndarray]]