03 — Patch Extraction and Normalisation

Purpose: Convert masked full-sky HEALPix maps into the flat-sky training patches consumed by train.py.

This notebook completes the preprocessing pipeline in three steps:

  1. Low-pass filtering — applies a sharp spherical-harmonic cutoff at ℓ = 7000 to remove aliasing artifacts, then zeros any negative pixels introduced by the filter (Gibbs ringing around masked regions).

  2. Flat-sky patch extraction — uses FlatCutter to project 6°×6° patches onto 256×256 Cartesian grids (1.41 arcmin/pixel) at centres computed by get_patch_centers with a 6° step size.

  3. Normalisation and saving — min-max normalises all patches jointly to [0, 1] and saves them as (N, 256, 256, 1) arrays. Also generates correlated Gaussian realisations from the measured CIB–tSZ cross-power spectrum for use as a non-Gaussian baseline in validation.

Inputs:

  • Masked CIB map: data/cib_150_masked.fits (from 02_masking.ipynb)

  • Masked tSZ map: data/tsz_150_masked.fits (from 02_masking.ipynb)

Outputs:

  • CIB patches: data/low_pass/2mJy/CIB_map_150GHz_256_st6_minmax_2mJy_zero_lp.npy

  • tSZ patches: data/low_pass/2mJy/tSZ3_map_150GHz_256_st6_minmax_2mJy_norm_lp.npy

  • Gaussian baseline patches: data/low_pass/2mJy/gaussian_cib_tsz_2mJy_lp.npy

Key module functions:

  • foregrounds_diffusion.preprocessing.get_patch_centers

  • foregrounds_diffusion.preprocessing.FlatCutter

  • foregrounds_diffusion.preprocessing.apply_maxmin_normalization

  • foregrounds_diffusion.flatmaps.make_gaussian_realisation

Paper reference: §2 (Figure 1 pipeline), §3.2 (Gaussian baseline construction).

[29]:

!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=b8e1b76c417a91c237d8f2edc311370937bfd9657f345b665d3c60625383ff07
  Stored in directory: /tmp/pip-ephem-wheel-cache-15415n7t/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.1
[notice] To update, run: pip install --upgrade pip

[30]:

# Remove all-zero patches: these arise when the patch footprint falls entirely
# within a masked region (Galactic plane, deep cluster mask).  Keeping them
# would introduce a large fraction of uninformative training samples.
# get_patch_centers returns a list of (RA, Dec) astropy SkyCoord objects on a
# quasi-uniform grid avoiding the Galactic plane (|b| < GAL_CUT) and the poles.
# The grid spacing STEP_DEG controls the patch density and overlap.
import numpy as np
import healpy as hp
import astropy.units as u
from pathlib import Path

from foregrounds_diffusion.preprocessing import (
    get_patch_centers, FlatCutter, apply_maxmin_normalization, apply_stdnorm,
)
from foregrounds_diffusion.flatmaps import map2cl, make_gaussian_realisation

PATCH_DEG  = 6.0   # patch side length in degrees
NRES       = 256   # pixels per side
STEP_DEG   = 6.0   # patch centre spacing in degrees
GAL_CUT    = 20.0  # Galactic-plane exclusion half-width in degrees
LMAX_FILTER = 7000  # low-pass cut-off multipole
PTSRC      = 2     # point-source threshold label (mJy)

flatskymapparams = [NRES, NRES, PATCH_DEG * 60. / NRES, PATCH_DEG * 60. / NRES]
# => [256, 256, 1.40625, 1.40625] arcmin/pixel
print(f"Pixel resolution: {flatskymapparams[2]:.4f} arcmin")

Pixel resolution: 1.4062 arcmin

1 Configuration

Patch geometry: NRES = 256 pixels per side at 1.40625′/px gives a 6° × 6° flat-sky projection per patch. A gnomonic (tangential) projection is used so that pixel area distortion is < 0.5% across the patch.

GAL_CUT = 20° removes patches centred within 20° of the Galactic plane, where the CIB and tSZ maps contain strong Galactic dust contamination not modelled by AGORA. STEP_DEG = 3 gives a 3° grid of patch centres, producing mild overlap between neighbouring patches and maximising sky coverage while keeping the training set large.

[31]:

PROJECT_ROOT = Path("/home/apb86/cmb_foregrounds_diffusion")
HPC_WORK = Path("/home/apb86/rds/hpc-work")

CIB_MASKED_FITS = PROJECT_ROOT / "data" / "cib_150_masked.fits"
TSZ_MASKED_FITS = PROJECT_ROOT / "data" / "tsz_150_masked.fits"

cib_map = hp.read_map(CIB_MASKED_FITS)
print(f"Loaded CIB — dtype: {cib_map.dtype}  min: {cib_map.min():.4e}  max: {cib_map.max():.4e}")

tsz_map = hp.read_map(TSZ_MASKED_FITS)
nside   = hp.get_nside(cib_map)

# Low-pass filter at ell = LMAX_FILTER via spherical harmonic transform
print(f"Applying low-pass filter at ell = {LMAX_FILTER} …")
cib_alm = hp.map2alm(cib_map, lmax=LMAX_FILTER)
print(f"ALMs — max abs: {np.abs(cib_alm).max():.4e}  any nan: {np.isnan(cib_alm).any()}")

tsz_alm = hp.map2alm(tsz_map, lmax=LMAX_FILTER)

cib_lp  = hp.alm2map(cib_alm, nside=nside)
print(f"Filtered CIB — min: {cib_lp.min():.4e}  max: {cib_lp.max():.4e}")

tsz_lp  = hp.alm2map(tsz_alm, nside=nside)

# Zero negative pixels (Gibbs ringing artefacts around masked regions)
n_neg_cib = (cib_lp < 0).sum()
print(f"Zeroing {n_neg_cib:,} negative CIB pixels ({100 * n_neg_cib / len(cib_lp):.2f}%)")
cib_lp[cib_lp < 0] = 0.

Loaded CIB — dtype: >f4  min: -9.8152e+00  max: 1.2193e+02
Applying low-pass filter at ell = 7000 …
ALMs — max abs: 7.4134e+01  any nan: False
Filtered CIB — min: -3.9952e+00  max: 1.0336e+02
Zeroing 39 negative CIB pixels (0.00%)

2 Generate patch centres

get_patch_centers returns a list of (RA, Dec) pairs on a quasi-uniform grid obeying the Galactic and polar cuts. The centre list is deterministic (no randomness), so the train/val/test split in the next steps is reproducible by controlling only the random seed passed to rng.permutation.

[32]:

# Remove all-zero patches: these arise when the patch footprint falls entirely
# within a masked region (Galactic plane, deep cluster mask).  Keeping them
# would introduce a large fraction of uninformative training samples.
# get_patch_centers returns a list of (RA, Dec) astropy SkyCoord objects on a
# quasi-uniform grid avoiding the Galactic plane (|b| < GAL_CUT) and the poles.
# The grid spacing STEP_DEG controls the patch density and overlap.
centers = get_patch_centers(gal_cut=GAL_CUT * u.deg, step_size=STEP_DEG * u.deg, pole_cut=PATCH_DEG * u.deg)
print(f"Patch centres (|b| > {GAL_CUT}°, step {STEP_DEG}°): {len(centers)}")

cutter = FlatCutter(ang_x=PATCH_DEG * u.deg, ang_y=PATCH_DEG * u.deg,
                    xres=NRES, yres=NRES)

cib_maps, tsz_maps = [], []
for lon, lat in centers:
    patch = cutter.rotate_to_pole_and_interpolate(lon, lat, [cib_lp, tsz_lp])
    cib_maps.append(patch[:, :, 0:1])
    tsz_maps.append(patch[:, :, 1:2])

cib_maps = np.stack(cib_maps)  # (N, 256, 256, 1)
tsz_maps = np.stack(tsz_maps)  # (N, 256, 256, 1)
print(f"Extracted {len(cib_maps)} patches, shape {cib_maps.shape}")

Patch centres (|b| > 20.0°, step 6.0°): 701
Extracted 701 patches, shape (701, 256, 256, 1)

3 Extract patches and filter empty regions

FlatCutter applies the gnomonic projection centred at each sky position and reads off a NRES × NRES pixel grid from the HEALPix map. Patches that land entirely in the Galactic plane or in the cluster-masked region will be all-zero after masking; we remove these to avoid training on empty patches.

[33]:

# Remove all-zero patches: these arise when the patch footprint falls entirely
# within a masked region (Galactic plane, deep cluster mask).  Keeping them
# would introduce a large fraction of uninformative training samples.
valid = np.array([p.max() > 0 for p in cib_maps])
n_removed = (~valid).sum()
print(f"Removing {n_removed} all-zero patches ({100 * n_removed / len(valid):.1f}%)")
cib_maps = cib_maps[valid]
tsz_maps = tsz_maps[valid]

# CIB: global Z-score
cib_mean = cib_maps.mean()
cib_std  = cib_maps.std()
cib_norm = (cib_maps - cib_mean) / cib_std

# tSZ: global standard normalisation across all patches
tsz_mean = tsz_maps.mean()
tsz_std  = tsz_maps.std()
tsz_norm = (tsz_maps - tsz_mean) / tsz_std

# Save normalisation parameters for denormalisation later
np.save(OUT_DIR / f"norm_params_{PTSRC}mJy.npy", np.array([cib_mean, cib_std, tsz_mean, tsz_std]))

print(id(tsz_norm))
print(tsz_norm.mean())

Removing 0 all-zero patches (0.0%)
139633529606896
-7.730383142932082e-16

4 Normalise and measure power spectra

CIB is min-max normalised to [0, 1] (the _zero suffix in filenames indicates the minimum is mapped to zero). tSZ is z-score normalised (zero mean, unit variance; the _norm suffix). The choice of different normalisation schemes for the two channels matches the paper: CIB pixels are non-negative by construction, making min-max natural; tSZ spans negative and positive values, making z-score more appropriate.

Power spectra are measured on a subsample of 200 patches to check that normalisation has not altered the spectral shape.

[34]:

# Measure mean auto- and cross-power spectra from normalised patches
LMIN, LMAX, LBIN = 300, 7000, 100
N_FOR_SPECTRA = min(200, len(cib_norm))

cl_cib_arr, cl_tsz_arr, cl_cross_arr = [], [], []
for i in range(N_FOR_SPECTRA):
    cib_m = cib_norm[i, :, :, 0]
    tsz_m = tsz_norm[i, :, :, 0]
    el, cl_c  = map2cl(flatskymapparams, cib_m, binsize=LBIN, minbin=LMIN, maxbin=LMAX)
    _,  cl_t  = map2cl(flatskymapparams, tsz_m, binsize=LBIN, minbin=LMIN, maxbin=LMAX)
    _,  cl_x  = map2cl(flatskymapparams, cib_m, tsz_m, binsize=LBIN, minbin=LMIN, maxbin=LMAX)
    cl_cib_arr.append(cl_c);  cl_tsz_arr.append(cl_t);  cl_cross_arr.append(cl_x)

cl_cib_mean  = np.mean(cl_cib_arr,  axis=0)
cl_tsz_mean  = np.mean(cl_tsz_arr,  axis=0)
cl_cross_mean = np.mean(cl_cross_arr, axis=0)
print(f"Mean CIB  C_ell range: {cl_cib_mean.min():.3e} – {cl_cib_mean.max():.3e}")
print(f"Mean tSZ  C_ell range: {cl_tsz_mean.min():.3e} – {cl_tsz_mean.max():.3e}")

# Generate correlated Gaussian realisations from the measured spectra
N_GAUSSIAN = len(cib_norm)
gaussian_maps = []
for _ in range(N_GAUSSIAN):
    sim = make_gaussian_realisation(
        flatskymapparams, el,
        cl=cl_cib_mean, cl2=cl_tsz_mean, cl12=cl_cross_mean,
    )
    gaussian_maps.append(sim[:2])  # keep only the two correlated fields

gaussian_arr = np.array(gaussian_maps)  # (N, 2, H, W)  channels-first
print(f"Gaussian baseline shape: {gaussian_arr.shape}")

Mean CIB  C_ell range: 3.256e-08 – 4.858e-06
Mean tSZ  C_ell range: 3.904e-09 – 6.381e-06
Gaussian baseline shape: (701, 2, 256, 256)

5 Save training arrays

Save four .npy files — normalised CIB patches, normalised tSZ patches, and the normalisation parameters — to data/low_pass/{ptsrc}mJy/. The arrays have shape (N, H, W, 1) (channels-last, one channel per file) so they can be stacked in train.py into (N, H, W, 2) before transposing to channels-first for PyTorch.

[35]:

OUT_DIR = Path("data/low_pass/2mJy")  # relative to where notebook runs
OUT_DIR.mkdir(parents=True, exist_ok=True)

np.save(OUT_DIR / f"CIB_map_150GHz_{NRES}_st6_zscore_{PTSRC}mJy_lp.npy", cib_norm)
np.save(OUT_DIR / f"tSZ3_map_150GHz_{NRES}_st6_zscore_{PTSRC}mJy_lp.npy", tsz_norm)
np.save(OUT_DIR / f"norm_params_{PTSRC}mJy.npy", np.array([cib_mean, cib_std, tsz_mean, tsz_std]))
np.save(OUT_DIR / f"gaussian_cib_tsz_{PTSRC}mJy_lp.npy", gaussian_arr)

[36]:

import matplotlib.pyplot as plt

N_SHOW = 9
fig, axes = plt.subplots(2, N_SHOW, figsize=(N_SHOW * 3, 6))

# Pick N_SHOW random patches
rng_vis = np.random.default_rng(seed=42)
idx = rng_vis.choice(len(cib_norm), size=N_SHOW, replace=False)

VMIN_CIB, VMAX_CIB = -3, 6
VMIN_TSZ, VMAX_TSZ = -5, 5

for col, i in enumerate(idx):
    # CIB
    ax = axes[0, col]
    im = ax.imshow(cib_norm[i, :, :, 0], origin='upper', cmap='magma',vmin=VMIN_CIB, vmax=VMAX_CIB)
    plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
    ax.set_title(f"CIB patch {i}", fontsize=8)
    ax.axis('off')

    # tSZ
    ax = axes[1, col]
    im = ax.imshow(tsz_norm[i, :, :, 0], origin='upper', cmap='magma',vmin=VMIN_TSZ, vmax=VMAX_TSZ)
    plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
    ax.set_title(f"tSZ patch {i}", fontsize=8)
    ax.axis('off')

axes[0, 0].set_ylabel("CIB", fontsize=10)
axes[1, 0].set_ylabel("tSZ", fontsize=10)

plt.suptitle("Sample CIB and tSZ training patches (Z-score normalised)", fontsize=12)
plt.tight_layout()
Path("plots").mkdir(parents=True, exist_ok=True)
plt.savefig("plots/sample_patches.png", dpi=150, bbox_inches='tight')
plt.show()

# Print basic statistics
print(f"CIB patches — shape: {cib_norm.shape}  min: {cib_norm.min():.4f}  "
      f"max: {cib_norm.max():.4f}  mean: {cib_norm.mean():.4f}  std: {cib_norm.std():.4f}")
print(f"tSZ patches — shape: {tsz_norm.shape}  min: {tsz_norm.min():.4f}  "
      f"max: {tsz_norm.max():.4f}  mean: {tsz_norm.mean():.4f}  std: {tsz_norm.std():.4f}")

# Flag any suspicious patches
masked_frac_cib = [(p == 0).mean() for p in cib_norm[:, :, :, 0]]
masked_frac_tsz = [(p == 0).mean() for p in tsz_norm[:, :, :, 0]]
print(f"\nPatches with >10% zero pixels:")
print(f"  CIB: {(np.array(masked_frac_cib) > 0.1).sum()}")
print(f"  tSZ: {(np.array(masked_frac_tsz) > 0.1).sum()}")
print(f"Mean zero fraction — CIB: {np.mean(masked_frac_cib):.3f}  "
      f"tSZ: {np.mean(masked_frac_tsz):.3f}")
../_images/tutorials_03_patch_extraction_13_0.png 
CIB patches — shape: (701, 256, 256, 1)  min: -19.0688  max: 10.5208  mean: -0.0000  std: 1.0000
tSZ patches — shape: (701, 256, 256, 1)  min: -36.3356  max: 17.3672  mean: -0.0000  std: 1.0000

Patches with >10% zero pixels:
  CIB: 0
  tSZ: 0
Mean zero fraction — CIB: 0.000  tSZ: 0.000

[37]:

print(f"tsz_maps raw — mean: {tsz_maps.mean():.4e}  std: {tsz_maps.std():.4e}  min: {tsz_maps.min():.4e}  max: {tsz_maps.max():.4e}")

tsz_patches raw — mean: -5.1154e+00  std: 3.7389e+00  min: -1.4097e+02  max: 5.9819e+01

[38]:

print(f"tsz_norm — mean: {tsz_norm.mean():.4f}  std: {tsz_norm.std():.4f}  min: {tsz_norm.min():.4f}  max: {tsz_norm.max():.4f}")

tsz_norm — mean: -0.0000  std: 1.0000  min: -36.3356  max: 17.3672

[ ]: