sys_mapping.utils

Angular two-point function measurement and covariance utilities.

Wrappers around TreeCorr and Corrfunc:

Utility functions for covariance estimation and two-point measurements.

Implements:

  • compute_covariance_matrix (Eq. 24): template auto-correlations

  • compute_amplitude_bias (Eq. 25-26): expected bias on clustering amplitude

  • measure_two_point_function: TreeCorr Landy-Szalay estimator wrapper

  • measure_two_point_function_corrfunc: Corrfunc Landy-Szalay estimator wrapper

  • measure_kk_correlation_treecorr: TreeCorr KK scalar-field correlator

  • measure_kk_correlation_corrfunc: Corrfunc KK scalar-field correlator

sys_mapping.utils.compute_covariance_matrix(delta_t)[source]

Compute the template covariance matrix (Eq. 24).

C_{ij} = (1/N_pix) Σ_p δ_{ti,p} δ_{tj,p}

Parameters:

delta_t ((n_sys, n_pix) template values at unmasked pixels)

Returns:

  • (n_sys, n_sys) symmetric positive-semi-definite covariance matrix

  • Performance

  • ———–

  • Measured on CPU (n_sys=5, n_pix=10_000) (*~120 μs/call.*)

  • Uses numpy matrix multiply; scales as O(n_sys² × n_pix).

  • Precision

  • ———

  • Matrix is symmetric to atol < 1e-14 (exact in float64).

  • All eigenvalues ≥ −1e-12 (positive semi-definite).

  • For unit-variance uncorrelated templates, result is close to the

  • identity matrix (deviation < 0.01 for n_pix ≥ 200_000).

Return type:

ndarray

Examples

>>> import numpy as np
>>> from sys_mapping import compute_covariance_matrix
>>> rng = np.random.default_rng(0)
>>> delta_t = rng.standard_normal((3, 10_000))
>>> C = compute_covariance_matrix(delta_t)
>>> C.shape
(3, 3)
>>> np.allclose(C, C.T)  # symmetric
True
sys_mapping.utils.compute_amplitude_bias(a_sq, b_sq, cov_matrix)[source]

Estimate bias on the galaxy clustering amplitude (Eq. 25-26).

Returns the additive offset and multiplicative factor that contamination imprints on the angle-averaged angular power spectrum.

Additive bias: Δw = Σ_i a²_i C_{ii} Multiplicative: f  = 1 + Σ_i b²_i C_{ii}

Parameters:

  • a_sq ((n_sys,) debiased squared additive parameters)

  • b_sq ((n_sys,) debiased squared multiplicative parameters)

  • cov_matrix ((n_sys, n_sys) template covariance from compute_covariance_matrix())

Returns:

  • additive_bias (float Δw added to the measured correlation function)

  • mult_factor (float factor by which w is inflated (1 = no bias))

  • Performance

  • ———–

  • Measured on CPU (n_sys=5) (*~11 μs/call. Two dot products.*)

  • Precision

  • ———

  • Matches analytic formula to rtol < 1e-10 (float64 dot products).

Return type:

tuple[float, float]

Examples

>>> import numpy as np
>>> from sys_mapping import compute_covariance_matrix, compute_amplitude_bias
>>> rng = np.random.default_rng(0)
>>> delta_t = rng.standard_normal((3, 10_000))
>>> C = compute_covariance_matrix(delta_t)
>>> a_sq = np.array([0.01, 0.02, 0.005])
>>> b_sq = np.array([0.005, 0.01, 0.002])
>>> add_bias, mult_factor = compute_amplitude_bias(a_sq, b_sq, C)
>>> add_bias     # additive offset on w(θ)
>>> mult_factor  # typically very close to 1.0
sys_mapping.utils.measure_two_point_function(ra_gal, dec_gal, ra_rand, dec_rand, min_sep=0.06, max_sep=30.0, nbins=15, sep_units='arcmin')[source]

Measure the angular two-point correlation function via Landy-Szalay.

Uses TreeCorr with log-spaced angular bins and the standard estimator w(θ) = (DD − 2DR + RR) / RR.

Parameters:

  • ra_gal ((n_gal,) galaxy positions in degrees)

  • dec_gal ((n_gal,) galaxy positions in degrees)

  • ra_rand ((n_rand,) random positions in degrees)

  • dec_rand ((n_rand,) random positions in degrees)

  • min_sep (float angular bin limits in sep_units)

  • max_sep (float angular bin limits in sep_units)

  • nbins (int number of log-spaced bins)

  • sep_units (str angular units for min_sep/max_sep (default 'arcmin'))

Returns:

  • theta ((nbins,) bin centers in sep_units)

  • w ((nbins,) Landy-Szalay angular correlation function estimate)

  • Performance

  • ———–

  • Dominated by TreeCorr pair-counting; scales as O(n_gal log n_gal).

  • Typical wall time (*seconds to minutes for survey-scale catalogs.*)

  • Precision

  • ———

  • Landy-Szalay estimator is unbiased to O(1/N). Statistical error is

  • approximately 1/sqrt(DD_pairs) per bin.

Return type:

tuple[ndarray, ndarray]

Examples

>>> import numpy as np
>>> from sys_mapping import measure_two_point_function
>>> rng = np.random.default_rng(0)
>>> ra_g   = rng.uniform(30, 60, 10_000)  # small patch
>>> dec_g  = rng.uniform(-10, 10, 10_000)
>>> ra_r   = rng.uniform(30, 60, 50_000)
>>> dec_r  = rng.uniform(-10, 10, 50_000)
>>> theta, w = measure_two_point_function(
...     ra_g, dec_g, ra_r, dec_r,
...     min_sep=1.0, max_sep=100.0, nbins=10, sep_units="arcmin",
... )
>>> theta.shape, w.shape
((10,), (10,))
sys_mapping.utils.measure_cross_two_point_function(ra1, dec1, ra2, dec2, ra_rand, dec_rand, w1=None, w2=None, min_sep=0.06, max_sep=30.0, nbins=15, sep_units='arcmin', metric='Arc', bin_slop=None, dr=None, rr=None)[source]

Measure the angular cross-correlation of two point samples (Landy-Szalay).

Estimates w_x(θ) between sample 1 (e.g. galaxies) and sample 2 (e.g. stars) that share a common footprint, using a single random catalog drawn uniformly over that footprint for both samples:

w_x(θ) = (D1D2 − D1R − R D2 + RR) / RR

computed with TreeCorr’s NNCorrelation.calculateXi(rr, dr, rd)() where dr = D1R and rd = R D2. Unlike measure_two_point_function() (auto-correlation), this counts distinct D1×D2 pairs, so it is the correct estimator for a galaxy×star signal.

The D1R (dr) and RR (rr) pair counts depend only on sample 1 and the randoms — not on sample 2 — so when scanning many sample-2 catalogs (e.g. magnitude-binned star samples) against a fixed sample 1, pass the returned dr/rr back in to skip recomputing them.

Parameters:

  • ra1 ((n1,) sample-1 positions in degrees)

  • dec1 ((n1,) sample-1 positions in degrees)

  • ra2 ((n2,) sample-2 positions in degrees)

  • dec2 ((n2,) sample-2 positions in degrees)

  • ra_rand ((n_rand,) random positions in degrees (common footprint))

  • dec_rand ((n_rand,) random positions in degrees (common footprint))

  • w1 (optional per-object weights for samples 1 and 2)

  • w2 (optional per-object weights for samples 1 and 2)

  • min_sep (float angular bin limits in sep_units)

  • max_sep (float angular bin limits in sep_units)

  • nbins (int number of log-spaced bins)

  • sep_units (str angular units for min_sep/max_sep (default 'arcmin'))

  • metric (str TreeCorr metric (default 'Arc', great-circle separation))

  • bin_slop (float or None TreeCorr bin_slop (None → TreeCorr default))

  • dr (optional pre-computed NNCorrelation objects (D1R and RR)) – to reuse across many sample-2 catalogs sharing the same sample 1 and randoms. When None they are computed and returned.

  • rr (optional pre-computed NNCorrelation objects (D1R and RR)) – to reuse across many sample-2 catalogs sharing the same sample 1 and randoms. When None they are computed and returned.

Returns:

  • theta ((nbins,) bin centers in sep_units)

  • w ((nbins,) Landy-Szalay cross-correlation estimate)

  • varxi ((nbins,) TreeCorr shot-noise variance per bin)

  • dr (NNCorrelation the D1R counts (reusable across sample-2 catalogs))

  • rr (NNCorrelation the RR counts (reusable across sample-2 catalogs))

Notes

For two statistically independent samples over the same footprint the estimator is consistent with w_x(θ) = 0 within shot noise — the basis for the GLASS null test.

Examples

>>> import numpy as np
>>> from sys_mapping import measure_cross_two_point_function
>>> rng = np.random.default_rng(0)
>>> ra1  = rng.uniform(30, 60, 5_000);  dec1 = rng.uniform(-10, 10, 5_000)
>>> ra2  = rng.uniform(30, 60, 5_000);  dec2 = rng.uniform(-10, 10, 5_000)
>>> ra_r = rng.uniform(30, 60, 50_000); dec_r = rng.uniform(-10, 10, 50_000)
>>> theta, w, var, dr, rr = measure_cross_two_point_function(
...     ra1, dec1, ra2, dec2, ra_r, dec_r,
...     min_sep=1.0, max_sep=100.0, nbins=8, sep_units="arcmin")
>>> theta.shape, w.shape
((8,), (8,))
sys_mapping.utils.measure_two_point_function_corrfunc(ra_gal, dec_gal, ra_rand, dec_rand, min_sep=0.06, max_sep=30.0, nbins=15, sep_units='arcmin', nthreads=1)[source]

Measure the angular two-point correlation function via Landy-Szalay.

Drop-in replacement for measure_two_point_function() using Corrfunc instead of TreeCorr. Pair-counting is exact (no bin_slop).

Parameters:

  • ra_gal ((n_gal,) galaxy positions in degrees)

  • dec_gal ((n_gal,) galaxy positions in degrees)

  • ra_rand ((n_rand,) random positions in degrees)

  • dec_rand ((n_rand,) random positions in degrees)

  • min_sep (float angular bin limits in sep_units)

  • max_sep (float angular bin limits in sep_units)

  • nbins (int number of log-spaced bins)

  • sep_units (str angular units (default 'arcmin'))

  • nthreads (int OpenMP threads passed to Corrfunc (default 1))

Returns:

  • theta ((nbins,) mean pair separation in sep_units)

  • w ((nbins,) Landy-Szalay angular correlation function estimate)

Return type:

tuple[ndarray, ndarray]

sys_mapping.utils.measure_kk_correlation_treecorr(ra, dec, k, w=None, *, ra2=None, dec2=None, k2=None, w2=None, min_sep=10.0, max_sep=250.0, nbins=20, sep_units='arcmin', bin_slop=0.01)[source]

Angular scalar-field correlation (KK) using TreeCorr.

Computes xi(theta) = <k1_i k2_j>_w for auto- or cross-correlation.

Parameters:

  • ra ((n,) positions in degrees)

  • dec ((n,) positions in degrees)

  • k ((n,) scalar field values (e.g. galaxy overdensity))

  • w ((n,) optional pair weights; ones if None)

  • ra2 (second catalog for cross-correlation; if None, auto-correlation)

  • dec2 (second catalog for cross-correlation; if None, auto-correlation)

  • k2 (second catalog for cross-correlation; if None, auto-correlation)

  • w2 (second catalog for cross-correlation; if None, auto-correlation)

  • min_sep (float bin limits in sep_units)

  • max_sep (float bin limits in sep_units)

  • nbins (int number of log-spaced bins)

  • sep_units (str angular units (default 'arcmin'))

  • bin_slop (float TreeCorr bin_slop parameter (default 0.01))

Returns:

  • theta ((nbins,) bin centres in sep_units)

  • xi ((nbins,) correlation function)

Return type:

tuple[ndarray, ndarray]

sys_mapping.utils.measure_kk_correlation_corrfunc(ra, dec, k, w=None, *, ra2=None, dec2=None, k2=None, w2=None, min_sep=10.0, max_sep=250.0, nbins=20, sep_units='arcmin', nthreads=1)[source]

Angular scalar-field correlation (KK) using Corrfunc.

Drop-in replacement for measure_kk_correlation_treecorr(). Computes xi(theta) = sum(w_i w_j k_i k_j) / sum(w_i w_j) per bin via two weighted pair-count passes with weight_type='pair_product'.

Parameters:

  • ra ((n,) positions in degrees)

  • dec ((n,) positions in degrees)

  • k ((n,) scalar field values)

  • w ((n,) optional pair weights; ones if None)

  • ra2 (second catalog for cross-correlation; if None, auto-correlation)

  • dec2 (second catalog for cross-correlation; if None, auto-correlation)

  • k2 (second catalog for cross-correlation; if None, auto-correlation)

  • w2 (second catalog for cross-correlation; if None, auto-correlation)

  • min_sep (float bin limits in sep_units)

  • max_sep (float bin limits in sep_units)

  • nbins (int number of log-spaced bins)

  • sep_units (str angular units (default 'arcmin'))

  • nthreads (int OpenMP threads (default 1))

Returns:

  • theta ((nbins,) mean pair separation in sep_units)

  • xi ((nbins,) correlation function)

Return type:

tuple[ndarray, ndarray]