07 — Higher-Order Statistics (Bispectrum and Trispectrum)
Purpose: Compute and compare collapsed equilateral bispectrum (S3) and trispectrum (S4) statistics for Agora, DDPM, and Gaussian samples.
Following Lee et al., the bispectrum and trispectrum are estimated from the skewness and kurtosis of harmonic band-filtered maps. This notebook covers both levels of analysis:
Summed channel — adds CIB and tSZ into a single map and computes S2, S3, S4 per ℓ-band (8 bands of width 720, centred from ℓ ≈ 300 to ℓ ≈ 5700). This is the primary result shown in the paper (Figure 7).
Joint cross-moments — computes all 12 channel combinations (S2^{aa}, S2^{bb}, S2^{ab}, S3^{aaa}, …, S4^{abbb}) separately for CIB and tSZ. This is the extended analysis shown in Appendix C (Figures 10–11).
Both analyses add three tiers of ILC residual noise (SPT-3G, S4-Wide, S4-Ultra Deep) before computing moments, and save intermediate results as .npy arrays so the slow computation only needs to run once.
Note: Computing moments for 800+ samples across 8 bands is slow (~hours). Run the moment computation cells once and reload from the saved arrays for plotting.
Inputs:
Test maps, DDPM samples, Gaussian baseline:
data/low_pass/2mJy/*.npyILC noise spectra:
data/ilc/ilc_weights_residuals_agora_fg_model.npy
Outputs:
Moment arrays:
data/moments_train_*.npy,data/moments_samples_*.npy,data/moments_gaussian_*.npy(shape: N × 8 × 12)Bispectrum/trispectrum comparison plots (Figures 7, 10, 11)
Key module functions:
foregrounds_diffusion.flatmaps.get_lpf_hpf(bandpass filter construction)foregrounds_diffusion.preprocessing.load_all_moments
Paper reference: §3.2 (moment statistics definition), §4.6 (Figure 7), Appendix C (Figures 10–11).
[10]:
!pip install -e ~/cmb_foregrounds_diffusion/
Obtaining file:///home/apb86/cmb_foregrounds_diffusion
Installing build dependencies ... done
Checking if build backend supports build_editable ... done
Getting requirements to build editable ... done
Preparing editable metadata (pyproject.toml) ... done
Requirement already satisfied: numpy>=1.26 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from foregrounds_diffusion==0.1.0) (1.26.4)
Requirement already satisfied: scipy>=1.11 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from foregrounds_diffusion==0.1.0) (1.17.0)
Requirement already satisfied: torch>=2.10 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from foregrounds_diffusion==0.1.0) (2.10.0)
Requirement already satisfied: torchvision>=0.25 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from foregrounds_diffusion==0.1.0) (0.25.0)
Requirement already satisfied: accelerate>=1.12 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from foregrounds_diffusion==0.1.0) (1.12.0)
Requirement already satisfied: denoising-diffusion-pytorch>=2.2.5 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from foregrounds_diffusion==0.1.0) (2.2.6)
Requirement already satisfied: healpy>=1.19 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from foregrounds_diffusion==0.1.0) (1.19.0)
Requirement already satisfied: astropy>=7.2 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from foregrounds_diffusion==0.1.0) (7.2.0)
Requirement already satisfied: einops>=0.8 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from foregrounds_diffusion==0.1.0) (0.8.2)
Requirement already satisfied: ema-pytorch>=0.7 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from foregrounds_diffusion==0.1.0) (0.7.9)
Requirement already satisfied: pillow>=12.0 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from foregrounds_diffusion==0.1.0) (12.1.1)
Requirement already satisfied: pytorch-fid>=0.3 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from foregrounds_diffusion==0.1.0) (0.3.0)
Requirement already satisfied: tqdm>=4.67 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from foregrounds_diffusion==0.1.0) (4.67.3)
Requirement already satisfied: packaging>=20.0 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from accelerate>=1.12->foregrounds_diffusion==0.1.0) (26.0)
Requirement already satisfied: psutil in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from accelerate>=1.12->foregrounds_diffusion==0.1.0) (7.2.2)
Requirement already satisfied: pyyaml in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from accelerate>=1.12->foregrounds_diffusion==0.1.0) (6.0.3)
Requirement already satisfied: huggingface_hub>=0.21.0 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from accelerate>=1.12->foregrounds_diffusion==0.1.0) (1.4.1)
Requirement already satisfied: safetensors>=0.4.3 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from accelerate>=1.12->foregrounds_diffusion==0.1.0) (0.7.0)
Requirement already satisfied: pyerfa>=2.0.1.1 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from astropy>=7.2->foregrounds_diffusion==0.1.0) (2.0.1.5)
Requirement already satisfied: astropy-iers-data>=0.2025.10.27.0.39.10 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from astropy>=7.2->foregrounds_diffusion==0.1.0) (0.2026.2.9.0.50.33)
Requirement already satisfied: filelock in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from huggingface_hub>=0.21.0->accelerate>=1.12->foregrounds_diffusion==0.1.0) (3.24.0)
Requirement already satisfied: fsspec>=2023.5.0 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from huggingface_hub>=0.21.0->accelerate>=1.12->foregrounds_diffusion==0.1.0) (2026.2.0)
Requirement already satisfied: hf-xet<2.0.0,>=1.2.0 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from huggingface_hub>=0.21.0->accelerate>=1.12->foregrounds_diffusion==0.1.0) (1.2.0)
Requirement already satisfied: httpx<1,>=0.23.0 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from huggingface_hub>=0.21.0->accelerate>=1.12->foregrounds_diffusion==0.1.0) (0.28.1)
Requirement already satisfied: shellingham in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from huggingface_hub>=0.21.0->accelerate>=1.12->foregrounds_diffusion==0.1.0) (1.5.4)
Requirement already satisfied: typer-slim in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from huggingface_hub>=0.21.0->accelerate>=1.12->foregrounds_diffusion==0.1.0) (0.23.1)
Requirement already satisfied: typing-extensions>=4.1.0 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from huggingface_hub>=0.21.0->accelerate>=1.12->foregrounds_diffusion==0.1.0) (4.15.0)
Requirement already satisfied: anyio in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from httpx<1,>=0.23.0->huggingface_hub>=0.21.0->accelerate>=1.12->foregrounds_diffusion==0.1.0) (4.12.1)
Requirement already satisfied: certifi in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from httpx<1,>=0.23.0->huggingface_hub>=0.21.0->accelerate>=1.12->foregrounds_diffusion==0.1.0) (2026.1.4)
Requirement already satisfied: httpcore==1.* in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from httpx<1,>=0.23.0->huggingface_hub>=0.21.0->accelerate>=1.12->foregrounds_diffusion==0.1.0) (1.0.9)
Requirement already satisfied: idna in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from httpx<1,>=0.23.0->huggingface_hub>=0.21.0->accelerate>=1.12->foregrounds_diffusion==0.1.0) (3.11)
Requirement already satisfied: h11>=0.16 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from httpcore==1.*->httpx<1,>=0.23.0->huggingface_hub>=0.21.0->accelerate>=1.12->foregrounds_diffusion==0.1.0) (0.16.0)
Requirement already satisfied: sympy>=1.13.3 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (1.14.0)
Requirement already satisfied: networkx>=2.5.1 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (3.6.1)
Requirement already satisfied: jinja2 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (3.1.6)
Requirement already satisfied: cuda-bindings==12.9.4 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (12.9.4)
Requirement already satisfied: nvidia-cuda-nvrtc-cu12==12.8.93 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (12.8.93)
Requirement already satisfied: nvidia-cuda-runtime-cu12==12.8.90 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (12.8.90)
Requirement already satisfied: nvidia-cuda-cupti-cu12==12.8.90 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (12.8.90)
Requirement already satisfied: nvidia-cudnn-cu12==9.10.2.21 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (9.10.2.21)
Requirement already satisfied: nvidia-cublas-cu12==12.8.4.1 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (12.8.4.1)
Requirement already satisfied: nvidia-cufft-cu12==11.3.3.83 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (11.3.3.83)
Requirement already satisfied: nvidia-curand-cu12==10.3.9.90 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (10.3.9.90)
Requirement already satisfied: nvidia-cusolver-cu12==11.7.3.90 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (11.7.3.90)
Requirement already satisfied: nvidia-cusparse-cu12==12.5.8.93 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (12.5.8.93)
Requirement already satisfied: nvidia-cusparselt-cu12==0.7.1 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (0.7.1)
Requirement already satisfied: nvidia-nccl-cu12==2.27.5 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (2.27.5)
Requirement already satisfied: nvidia-nvshmem-cu12==3.4.5 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (3.4.5)
Requirement already satisfied: nvidia-nvtx-cu12==12.8.90 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (12.8.90)
Requirement already satisfied: nvidia-nvjitlink-cu12==12.8.93 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (12.8.93)
Requirement already satisfied: nvidia-cufile-cu12==1.13.1.3 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (1.13.1.3)
Requirement already satisfied: triton==3.6.0 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from torch>=2.10->foregrounds_diffusion==0.1.0) (3.6.0)
Requirement already satisfied: cuda-pathfinder~=1.1 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from cuda-bindings==12.9.4->torch>=2.10->foregrounds_diffusion==0.1.0) (1.3.4)
Requirement already satisfied: mpmath<1.4,>=1.1.0 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from sympy>=1.13.3->torch>=2.10->foregrounds_diffusion==0.1.0) (1.3.0)
Requirement already satisfied: MarkupSafe>=2.0 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from jinja2->torch>=2.10->foregrounds_diffusion==0.1.0) (3.0.3)
Requirement already satisfied: typer>=0.23.1 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from typer-slim->huggingface_hub>=0.21.0->accelerate>=1.12->foregrounds_diffusion==0.1.0) (0.23.1)
Requirement already satisfied: click>=8.0.0 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from typer>=0.23.1->typer-slim->huggingface_hub>=0.21.0->accelerate>=1.12->foregrounds_diffusion==0.1.0) (8.3.1)
Requirement already satisfied: rich>=10.11.0 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from typer>=0.23.1->typer-slim->huggingface_hub>=0.21.0->accelerate>=1.12->foregrounds_diffusion==0.1.0) (14.3.2)
Requirement already satisfied: annotated-doc>=0.0.2 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from typer>=0.23.1->typer-slim->huggingface_hub>=0.21.0->accelerate>=1.12->foregrounds_diffusion==0.1.0) (0.0.4)
Requirement already satisfied: markdown-it-py>=2.2.0 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from rich>=10.11.0->typer>=0.23.1->typer-slim->huggingface_hub>=0.21.0->accelerate>=1.12->foregrounds_diffusion==0.1.0) (4.0.0)
Requirement already satisfied: pygments<3.0.0,>=2.13.0 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from rich>=10.11.0->typer>=0.23.1->typer-slim->huggingface_hub>=0.21.0->accelerate>=1.12->foregrounds_diffusion==0.1.0) (2.19.2)
Requirement already satisfied: mdurl~=0.1 in /home/apb86/diffusion_project_env/lib64/python3.11/site-packages (from markdown-it-py>=2.2.0->rich>=10.11.0->typer>=0.23.1->typer-slim->huggingface_hub>=0.21.0->accelerate>=1.12->foregrounds_diffusion==0.1.0) (0.1.2)
Building wheels for collected packages: foregrounds_diffusion
Building editable for foregrounds_diffusion (pyproject.toml) ... done
Created wheel for foregrounds_diffusion: filename=foregrounds_diffusion-0.1.0-0.editable-py3-none-any.whl size=4647 sha256=9d3bb3d6aa21f816caa1c61ecba54c7ac9d3d800d32df0d942fb8567eb370553
Stored in directory: /tmp/pip-ephem-wheel-cache-hqaxvf3s/wheels/ad/05/47/d622ea03cc2c607997f65b7643dab9c4f27fa5f37d82097115
Successfully built foregrounds_diffusion
Installing collected packages: foregrounds_diffusion
Attempting uninstall: foregrounds_diffusion
Found existing installation: foregrounds_diffusion 0.1.0
Uninstalling foregrounds_diffusion-0.1.0:
Successfully uninstalled foregrounds_diffusion-0.1.0
Successfully installed foregrounds_diffusion-0.1.0
[notice] A new release of pip is available: 26.0.1 -> 26.1.2
[notice] To update, run: pip install --upgrade pip
1 Configuration
Band-pass filter parameters. We tile [ℓ_min, ℓ_max] = [300, 6060] into N_BANDS = 8 uniform bands of width Δℓ = 720, centred at ℓ ≈ 660, 1380, …,
These choices follow Lee et al. and the paper’s Figure 7. The Gaussian difference filter (
filter_type=2) is a band-pass with smooth roll-off that avoids the Gibbs ringing of a hard top-hat.
[ ]:
# compute_summed_moments(cib_arr + tsz_arr, bp_filters):
# For each band b and map i, band-filter (cib_i + tsz_i) and compute:
# col 0 = S2 = var(filtered_map) — proportional to C_ell
# col 1 = S3 = mean(f^3) / var(f)^1.5 — skewness
# col 2 = S4 = mean(f^4) / var(f)^2 - 3 — excess kurtosis
# Returns ndarray(N, B, 3).
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
from foregrounds_diffusion.flatmaps import get_lpf_hpf, bandpass_filter
from foregrounds_diffusion.preprocessing import load_all_moments, renormalize_dm_maps, denormalize_dm_maps
from foregrounds_diffusion.moments import compute_summed_moments, compute_cross_moments
PTSRC = 2
N_BANDS = 8
LMIN_BAND = 300
LMAX_BAND = 6060
BAND_WIDTH = (LMAX_BAND - LMIN_BAND) // N_BANDS # 720
flatskymapparams = [256, 256, 1.40625, 1.40625]
PROJECT_ROOT = Path("/home/apb86/cmb_foregrounds_diffusion")
PATCHES_DIR = Path(f"data/low_pass/{PTSRC}mJy")
bp_edges = [(LMIN_BAND + i * BAND_WIDTH, LMIN_BAND + (i + 1) * BAND_WIDTH)
for i in range(N_BANDS)]
bandpass_centers = np.array([0.5 * (lo + hi) for lo, hi in bp_edges])
print("Bandpass centres:", bandpass_centers.astype(int))
2 Pre-compute band-pass filters
get_lpf_hpf with filter_type=2 returns a 2D Fourier-space filter array of shape (ny, nx) whose values lie in [0, 1]. Pre-computing these once and passing them as a list to compute_summed_moments / compute_cross_moments avoids redundant FFT frequency grid construction during the per-map loop. Store the filter list as bp_filters — it is reused for all three sample populations (Agora, DDPM, Gaussian).
[ ]:
# Pre-compute the 2D band-pass filters once. Each filter is a (ny, nx) float
# array in Fourier space with values in [0,1]. filter_type=2 gives a Gaussian
# difference (DoG) filter, which has a smooth roll-off that avoids ringing
# artefacts at the band edges. Storing the list avoids repeated FFT-grid
# construction inside the moment-computation loops.
# Pre-compute all bandpass filters (stored as 2D arrays)
bp_filters = [get_lpf_hpf(flatskymapparams, lminmax, filter_type=2)
for lminmax in bp_edges]
3 Optional: ILC noise spectra
Instrument noise from the SPT-3G and CMB-S4 experiments can be added to each sample before computing statistics, to simulate realistic map-level signal-to-noise. The ILC residual spectra are pre-computed from a Wiener-filtered combination of frequency channels. This cell loads those spectra but does not apply them in the main analysis — uncomment the example line to add a noise realisation to any map.
[13]:
# Load ILC noise spectra (SPT-3G, S4-Wide, S4-Ultra Deep)
ilc_noise = np.load(
PROJECT_ROOT / "data" / "ilc" / "ilc_weights_residuals_agora_fg_model.npy", allow_pickle=True
).item()
# Example: add SPT-3G noise realisation to a map
# noise_map = cl2map(flatskymapparams, ilc_noise['spt3g_150'], el=ilc_noise['el'])
4 Load and denormalise maps
Identical loading pattern to notebook 06. The Agora test set is defined by the same rng.permutation(seed=42) split used during training, so there is no data leakage.
[ ]:
cib_maps = np.load(PATCHES_DIR / f"CIB_map_150GHz_256_st6_minmax_{PTSRC}mJy_zero_lp.npy")
tsz_maps = np.load(PATCHES_DIR / f"tSZ3_map_150GHz_256_st6_minmax_{PTSRC}mJy_norm_lp.npy")
ddpm_raw = np.load(PROJECT_ROOT / "data" / "low_pass" / f"{PTSRC}mJy" / f"new_samples_cib_tsz_{PTSRC}mJy_zero_norm_6x6_w_au_lp.npy")
gauss_maps = np.load(PATCHES_DIR / f"gaussian_cib_tsz_{PTSRC}mJy_lp.npy")
agora_cib = cib_maps[:, :, :, 0]
agora_tsz = tsz_maps[:, :, :, 0]
train_maps = np.concatenate([cib_maps, tsz_maps], axis=-1) # (N, H, W, 2)
print(f"train_maps shape: {train_maps.shape}")
print(f"ddpm_raw shape: {ddpm_raw.shape}")
# Load saved normalisation parameters from patch extraction
norm_params = np.load(PATCHES_DIR / f"norm_params_{PTSRC}mJy.npy")
cib_mean, cib_std, tsz_mean, tsz_std = norm_params
print(f"CIB — mean: {cib_mean:.4f} std: {cib_std:.4f}")
print(f"tSZ — mean: {tsz_mean:.4f} std: {tsz_std:.4f}")
ddpm_renorm = denormalize_dm_maps(ddpm_raw, cib_mean, cib_std, tsz_mean, tsz_std)
print(f"Denormalised DDPM — CIB range: [{ddpm_renorm[:,0].min():.4f}, {ddpm_renorm[:,0].max():.4f}]")
print(f"Denormalised DDPM — tSZ range: [{ddpm_renorm[:,1].min():.4f}, {ddpm_renorm[:,1].max():.4f}]")
# DDPM samples are (N, 2, H, W) channels-first; Agora/Gaussian are (N, H, W, 1)
ddpm_cib = ddpm_renorm[:, 0] # (N, H, W)
ddpm_tsz = ddpm_renorm[:, 1]
if gauss_maps.shape[1] == 2:
gauss_cib = gauss_maps[:, 0]; gauss_tsz = gauss_maps[:, 1]
else:
gauss_cib = gauss_maps[:, :, :, 0]; gauss_tsz = gauss_maps[:, :, :, 1]
5 Compute summed-channel moments
compute_summed_moments(cib + tsz, bp_filters) band-filters the sum CIB + tSZ, then computes the variance S2, skewness S3, and excess kurtosis S4 of the filtered field. A Gaussian field has S3 = S4 = 0 by definition, so non-zero values indicate non-Gaussianity. The expected ordering is |S3_Agora| > |S3_DDPM| > |S3_Gaussian| if the DDPM captures some but not all non-Gaussian structure. This computation takes ≈ 30 min for N = 200 maps; reduce N_TEST for a quick check.
[ ]:
# compute_summed_moments(cib_arr + tsz_arr, bp_filters):
# For each band b and map i, band-filter (cib_i + tsz_i) and compute:
# col 0 = S2 = var(filtered_map) — proportional to C_ell
# col 1 = S3 = mean(f^3) / var(f)^1.5 — skewness
# col 2 = S4 = mean(f^4) / var(f)^2 - 3 — excess kurtosis
# Returns ndarray(N, B, 3).
# NOTE: this takes ~30 min for full datasets — compute once and save
print("Computing Agora moments … (set N=50 for a quick test)")
N_TEST = 200 # replace with -1 for full run
agora_moments = compute_summed_moments(
agora_cib[:N_TEST], agora_tsz[:N_TEST], bp_filters)
ddpm_moments = compute_summed_moments(
ddpm_cib[:N_TEST], ddpm_tsz[:N_TEST], bp_filters)
gauss_moments = compute_summed_moments(
gauss_cib[:N_TEST], gauss_tsz[:N_TEST], bp_filters)
6 Compute cross-channel moments
compute_cross_moments computes all 12 mixed moments between the CIB (a) and tSZ (b) channels per ℓ-band: three second-order (S2^{aa}, S2^{bb}, S2^{ab}), four third-order, and five fourth-order terms. The second-order terms correspond to the auto- and cross-power spectra after band-filtering. The higher-order terms capture non-linear coupling between the two foreground components that cannot be recovered from the power spectrum alone.
[ ]:
# compute_cross_moments(cib_arr, tsz_arr, bp_filters) → (ndarray(N, B, 12), labels)
# The 12 labels are: S2aa, S2bb, S2ab, S3aaa, S3bbb, S3aab, S3abb,
# S4aaaa, S4bbbb, S4aaab, S4aabb, S4abbb
# where a = CIB channel, b = tSZ channel.
# NOTE: this takes ~30-60 min for full datasets — compute once and save
print("Computing Agora cross-moments …")
agora_cross, cross_labels = compute_cross_moments(
agora_cib[:N_TEST], agora_tsz[:N_TEST], bp_filters)
print("Computing DDPM cross-moments …")
ddpm_cross, _ = compute_cross_moments(
ddpm_cib[:N_TEST], ddpm_tsz[:N_TEST], bp_filters)
print("Computing Gaussian cross-moments …")
gauss_cross, _ = compute_cross_moments(
gauss_cib[:N_TEST], gauss_tsz[:N_TEST], bp_filters)
print(f"Done. Shape: {agora_cross.shape} — labels: {cross_labels}")
7 Save and reload
Moment computation is the most expensive step (≈ 30–60 min for N = 200 maps on a single CPU core). Save results to disk immediately after computation so that subsequent plotting cells can be re-run quickly without recomputing. The commented np.load lines below illustrate how to reload from disk.
[ ]:
# Save and reload — run the computation cells once, then reload from disk
OUT_DIR = Path("data")
np.save(OUT_DIR / f"moments_train_summed_{PTSRC}mJy.npy", agora_moments)
np.save(OUT_DIR / f"moments_samples_summed_{PTSRC}mJy.npy", ddpm_moments)
np.save(OUT_DIR / f"moments_gaussian_summed_{PTSRC}mJy.npy", gauss_moments)
np.save(OUT_DIR / f"moments_train_cross_{PTSRC}mJy.npy", agora_cross)
np.save(OUT_DIR / f"moments_samples_cross_{PTSRC}mJy.npy", ddpm_cross)
np.save(OUT_DIR / f"moments_gaussian_cross_{PTSRC}mJy.npy", gauss_cross)
# Reload (uncomment to skip recomputation)
# agora_cross = np.load(OUT_DIR / f"moments_train_cross_{PTSRC}mJy.npy")
# ddpm_cross = np.load(OUT_DIR / f"moments_samples_cross_{PTSRC}mJy.npy")
# gauss_cross = np.load(OUT_DIR / f"moments_gaussian_cross_{PTSRC}mJy.npy")
# cross_labels = ['S2aa', 'S2bb', 'S2ab',
# 'S3aaa', 'S3bbb', 'S3aab', 'S3abb',
# 'S4aaaa', 'S4bbbb', 'S4aaab', 'S4aabb', 'S4abbb']
8 Plot summed-channel moments
[18]:
MOMENT_LABELS = ['S2', 'S3 (skewness)', 'S4 (excess kurtosis)']
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for ax, col, label in zip(axes, range(3), MOMENT_LABELS):
ax.plot(bandpass_centers, agora_moments[:, :, col].mean(axis=0), label="Agora", color="C0")
ax.plot(bandpass_centers, ddpm_moments[:, :, col].mean(axis=0), label="DDPM", color="C1")
ax.plot(bandpass_centers, gauss_moments[:, :, col].mean(axis=0), label="Gaussian", color="C2", ls="--")
ax.set_xlabel(r"$\ell$ (band centre)"); ax.set_title(label); ax.legend()
plt.tight_layout(); plt.show()
Cross-channel bi- and trispectrum (all CIB × tSZ combinations)
The 9 higher-order cross-moments between the CIB channel (a) and tSZ channel (b), plotted per ℓ-band. Bispectra (S3, top two rows) and trispectra (S4, bottom row). Gaussian baseline is dashed — any deviation indicates non-Gaussianity.
[ ]:
# Indices 3-11 in the cross-moments array are the higher-order terms
# (indices 0-2 are S2aa, S2bb, S2ab — second-order, plotted separately if needed)
HIGHER_ORDER_IDX = list(range(3, 12))
LATEX_LABELS = [
r'$S_3^{aaa}$', r'$S_3^{bbb}$', r'$S_3^{aab}$', r'$S_3^{abb}$',
r'$S_4^{aaaa}$', r'$S_4^{bbbb}$', r'$S_4^{aaab}$', r'$S_4^{aabb}$', r'$S_4^{abbb}$',
]
sources = [
('Agora', agora_cross, 'C0', '-'),
('DDPM', ddpm_cross, 'C1', '-'),
('Gaussian', gauss_cross, 'C2', '--'),
]
fig, axes = plt.subplots(3, 3, figsize=(15, 11), sharex=True)
axes = axes.flatten()
for ax, idx, latex in zip(axes, HIGHER_ORDER_IDX, LATEX_LABELS):
for label, data, color, ls in sources:
vals = data[:, :, idx] # (N, n_bands)
mean = vals.mean(axis=0)
std = vals.std(axis=0) / np.sqrt(len(vals))
ax.plot(bandpass_centers, mean, color=color, ls=ls, lw=1.5, label=label)
ax.fill_between(bandpass_centers, mean - std, mean + std, alpha=0.2, color=color)
ax.axhline(0, color='k', lw=0.6, ls=':')
ax.set_title(latex, fontsize=13)
ax.grid(True, ls=':', lw=0.4, alpha=0.5)
for ax in axes[6:]:
ax.set_xlabel(r'$\ell_c$ (band centre)', fontsize=11)
for ax in axes[::3]:
ax.set_ylabel('Moment value', fontsize=11)
handles, lbls = axes[0].get_legend_handles_labels()
fig.legend(handles, lbls, loc='lower center', ncol=3, fontsize=12,
bbox_to_anchor=(0.5, -0.02), frameon=False)
# Label the channels
fig.text(0.01, 0.97, r'$a$ = CIB, $b$ = tSZ', fontsize=11, va='top')
plt.suptitle('Higher-order cross-moments: bispectrum (S3) and trispectrum (S4)', fontsize=14)
plt.tight_layout(rect=[0, 0.04, 1, 0.97])
Path('plots').mkdir(exist_ok=True)
plt.savefig('plots/cross_moments_bispec_trispec.pdf', dpi=200, bbox_inches='tight')
plt.show()
[ ]: