moments — Power spectra and higher-order statistics

moments is the primary validation module. It computes the summary statistics used throughout the paper to compare DDPM-generated maps against AGORA training maps and a Gaussian baseline.

Power spectra

mean_cls() and mean_cross_cls() average power spectra over an ensemble of maps. Both functions cache the ℓ-bin grid and support joblib parallelism via n_jobs=-1:

from foregrounds_diffusion.moments import mean_cls, mean_cross_cls

mapparams = [256, 256, 1.40625, 1.40625]

# Auto-spectrum of CIB maps
el, cl_cib, err_cib = mean_cls(cib_maps, mapparams, lmin=300, lmax=4000,
                                binsize=60, n_jobs=-1)

# Cross-spectrum between CIB and tSZ
el, cl_cross, err_cross = mean_cross_cls(cib_maps, tsz_maps, mapparams,
                                          lmin=300, lmax=4000, binsize=60)

Higher-order moments

compute_cross_moments() computes 12 cross-moments of bandpass-filtered CIB (a) and tSZ (b) fields per ℓ-band — variance, skewness, and kurtosis terms up to 4th order:

from foregrounds_diffusion.flatmaps import get_lpf_hpf
from foregrounds_diffusion.moments import compute_cross_moments

# Build one bandpass filter per ℓ-band
bp_edges = [(300 + i*720, 300 + (i+1)*720) for i in range(8)]
bp_filters = [get_lpf_hpf(mapparams, edges, filter_type=2) for edges in bp_edges]

moments, labels = compute_cross_moments(cib_maps, tsz_maps, bp_filters, n_jobs=-1)
# moments: (N, 8, 12),  labels: ['S2aa', 'S2bb', 'S2ab', ...]

See 07 — Higher-Order Statistics (Bispectrum and Trispectrum) for the full comparison between AGORA, DDPM, and Gaussian maps.

API

Power-spectrum and higher-order moment statistics.

This module computes the summary statistics used to validate DDPM samples against the AGORA training maps (paper §4.1 and §4.2). All functions accept stacks of flat-sky maps and return statistics averaged over the ensemble, making them suitable for comparing large sets of DDPM samples.

Power spectra

mean_cls() — mean auto-power spectrum C_ℓ over N maps. mean_cross_cls() — mean cross-power spectrum between two map stacks.

Both functions accept n_jobs=-1 to parallelise over maps using joblib, and share a pre-computed ℓ-bin cache to amortise the FFT grid construction.

Higher-order moments

compute_summed_moments() — variance (S2), skewness (S3), and excess kurtosis (S4) of the bandpass-filtered sum field CIB + tSZ. Equivalent to the collapsed bispectrum and trispectrum proxies used in the paper.

compute_cross_moments() — the full set of 12 cross-moments between bandpass-filtered CIB (a) and tSZ (b) fields:

S2^{aa}, S2^{bb}, S2^{ab},
S3^{aaa}, S3^{bbb}, S3^{aab}, S3^{abb},
S4^{aaaa}, S4^{bbbb}, S4^{aaab}, S4^{aabb}, S4^{abbb}

These are computed per ℓ-band by applying 2D bandpass filters from get_lpf_hpf() before taking moments. The output shape is (N, n_bands, 12) where N is the number of maps.

See tutorial docs/tutorials/07_higher_order_stats.ipynb for the full comparison between AGORA, DDPM samples, and a Gaussian baseline.

foregrounds_diffusion.moments.mean_cls(maps_nhw, mapparams, lmin, lmax, binsize, n_jobs=1)[source]

Compute mean auto-power spectrum over a stack of maps.

Parameters:

maps_nhw (ndarray, shape (N, H, W)) – Stack of flat-sky maps.
mapparams (list) – [nx, ny, dx, dy] — see get_lxly().
lmin, lmax (float) – Multipole range.
binsize (float) – Bin width in ℓ.
n_jobs (int) – Number of parallel workers (joblib). 1 = serial (default). -1 = all cores.

Returns:

el (ndarray) – Bin centres.
mean_cl (ndarray) – Mean power spectrum across maps.
std_cl (ndarray) – Standard deviation across maps.

foregrounds_diffusion.moments.mean_cross_cls(maps1, maps2, mapparams, lmin, lmax, binsize, n_jobs=1)[source]

Compute mean cross-power spectrum between two stacks of maps.

Parameters:

maps1, maps2 (ndarray, shape (N, H, W)) – Two stacks of flat-sky maps.
mapparams (list) – [nx, ny, dx, dy].
lmin, lmax (float) – Multipole range.
binsize (float) – Bin width in ℓ.
n_jobs (int) – Number of parallel workers (joblib). 1 = serial (default). -1 = all cores.

Returns:

el (ndarray)
mean_cl (ndarray)
std_cl (ndarray)

foregrounds_diffusion.moments.compute_summed_moments(cib_arr, tsz_arr, bp_filters, n_jobs=1)[source]

Compute S2, S3, S4 of the summed CIB+tSZ field per ℓ-band.

Parameters:

cib_arr (ndarray, shape (N, H, W))
tsz_arr (ndarray, shape (N, H, W))
bp_filters (list of ndarray) – 2D bandpass filters from get_lpf_hpf().
n_jobs (int) – Number of parallel workers (joblib). 1 = serial (default). -1 = all cores.

Returns:

ndarray, shape (N, len(bp_filters), 3) – Columns: variance (S2), skewness (S3), excess kurtosis (S4).

foregrounds_diffusion.moments.compute_cross_moments(cib_arr, tsz_arr, bp_filters, n_jobs=1)[source]

Compute all 12 cross-moments per ℓ-band (a=CIB, b=tSZ).

Moments: S2^{aa}, S2^{bb}, S2^{ab},: S3^{aaa}, S3^{bbb}, S3^{aab}, S3^{abb}, S4^{aaaa}, S4^{bbbb}, S4^{aaab}, S4^{aabb}, S4^{abbb}.

Parameters:

cib_arr (ndarray, shape (N, H, W))
tsz_arr (ndarray, shape (N, H, W))
bp_filters (list of ndarray)
n_jobs (int) – Number of parallel workers (joblib). 1 = serial (default). −1 = use all cores.

Returns:

moments (ndarray, shape (N, len(bp_filters), 12))
labels (list of str)