peak_counts — Peak and minima counting

peak_counts implements the peak and minima counting statistics from Sabyr, Hill & Haiman (2024, arXiv:2410.21247), adapted for flat-sky patches.

Peak counts measure how many pixels exceed a given signal-to-noise threshold ν = T/σ, and minima counts measure how many pixels fall below -ν. Together they characterise the one-point distribution and local topology of the field in a way that is sensitive to non-Gaussianity — a Gaussian field has a symmetric peak/minima distribution, while tSZ maps are skewed toward positive peaks due to massive clusters.

Requires only numpy and scipy — no LensTools dependency.

Usage

from foregrounds_diffusion.peak_counts import compute_peak_minima_counts

nu_bins = np.linspace(-4, 4, 17)   # ν = T/σ thresholds
smoothing_arcmin = 1.0             # Gaussian smoothing scale

result = compute_peak_minima_counts(
    maps_nhw, nu_bins, smoothing_arcmin, dx_arcmin=1.40625
)
# result["peak_counts"]: (N, len(nu_bins)-1)
# result["minima_counts"]: (N, len(nu_bins)-1)
# result["nu_centres"]: bin centres

To smooth at multiple angular scales before counting:

from foregrounds_diffusion.peak_counts import smooth_map, find_peaks, count_peaks_binned

smoothed = smooth_map(maps_nhw, sigma_arcmin=2.0, dx_arcmin=1.40625)
peaks = find_peaks(smoothed)                           # local maxima coordinates
counts = count_peaks_binned(smoothed, peaks, nu_bins) # (N, n_bins)

See 10 — Peak and Minima Counts for comparisons between AGORA, DDPM, and Gaussian maps across multiple smoothing scales.

API

Peak and minima counts for flat-sky CMB foreground maps.

Implements the peak and minima counting statistics from Sabyr, Hill & Haiman (2024), arXiv:2410.21247, adapted for flat-sky patches rather than full-sky maps.

The pipeline:

Smooth each patch with a Gaussian kernel at one or more angular scales.
Identify local maxima (peaks) and local minima using scipy’s maximum_filter / minimum_filter.
Bin the peak/minima pixel values by their amplitude in units of the map standard deviation ν = T / σ.
Return counts as a function of ν for comparison across map ensembles.

No external packages beyond numpy and scipy are required — the LensTools dependency used by Sabyr et al. is not needed for flat-sky patches.

foregrounds_diffusion.peak_counts.smooth_map(patch, fwhm_arcmin, pixel_res_arcmin=1.40625)[source]

Smooth a 2D patch with a Gaussian kernel.

Parameters:

patch (ndarray, shape (H, W)) – Input flat-sky map in any units.
fwhm_arcmin (float) – Full-width at half-maximum of the Gaussian smoothing kernel in arcmin.
pixel_res_arcmin (float) – Pixel resolution in arcmin/pixel (default: 6°/256px ≈ 1.406 arcmin).

Returns:

ndarray, shape (H, W) – Smoothed map.

foregrounds_diffusion.peak_counts.find_peaks(patch, filter_size=3)[source]

Find local maxima (peaks) in a 2D map.

A pixel is a peak if it equals the maximum of its neighbourhood defined by filter_size × filter_size.

Parameters:

patch (ndarray, shape (H, W)) – Smoothed 2D map.
filter_size (int) – Side length of the sliding window for local maximum detection.

Returns:

ndarray of float – Pixel values at all peak locations.

foregrounds_diffusion.peak_counts.find_minima(patch, filter_size=3)[source]

Find local minima in a 2D map.

Parameters:

patch (ndarray, shape (H, W)) – Smoothed 2D map.
filter_size (int) – Side length of the sliding window for local minimum detection.

Returns:

ndarray of float – Pixel values at all minima locations.

foregrounds_diffusion.peak_counts.count_peaks_binned(patches_nhw, thresholds, fwhm_arcmin, pixel_res_arcmin=1.40625, filter_size=3)[source]

Compute mean peak counts per map as a function of threshold ν.

Following Sabyr et al. (2024), thresholds are defined in units of the per-map standard deviation: ν = T / σ. Peaks with ν > threshold are counted, giving a cumulative count curve.

Parameters:

patches_nhw (ndarray, shape (N, H, W)) – Stack of flat-sky patches.
thresholds (array_like) – Threshold values in units of σ (e.g. np.linspace(-4, 4, 40)).
fwhm_arcmin (float) – Gaussian smoothing scale in arcmin applied before peak finding.
pixel_res_arcmin (float) – Pixel resolution in arcmin.
filter_size (int) – Neighbourhood size for local maximum detection.

Returns:

counts (ndarray, shape (N, len(thresholds))) – Peak counts per map per threshold bin.

foregrounds_diffusion.peak_counts.count_minima_binned(patches_nhw, thresholds, fwhm_arcmin, pixel_res_arcmin=1.40625, filter_size=3)[source]

Compute mean minima counts per map as a function of threshold ν.

Minima with ν < threshold are counted (threshold should be negative).

Parameters:

patches_nhw (ndarray, shape (N, H, W)) – Stack of flat-sky patches.
thresholds (array_like) – Threshold values in units of σ (e.g. np.linspace(-4, 0, 20)).
fwhm_arcmin (float) – Gaussian smoothing scale in arcmin.
pixel_res_arcmin (float) – Pixel resolution in arcmin.
filter_size (int) – Neighbourhood size for local minimum detection.

Returns:

counts (ndarray, shape (N, len(thresholds))) – Minima counts per map per threshold bin.

foregrounds_diffusion.peak_counts.compute_peak_minima_counts(patches_nhw, thresholds_peaks, thresholds_minima, smoothing_scales_arcmin=(1.0, 2.5, 5.0), pixel_res_arcmin=1.40625, filter_size=3, n_jobs=1)[source]

Compute peak and minima counts at multiple smoothing scales.

Parameters:

patches_nhw (ndarray, shape (N, H, W)) – Stack of flat-sky patches in physical units (e.g. µK).
thresholds_peaks (array_like) – Threshold values ν for peak counting (e.g. np.linspace(-1, 5, 30)).
thresholds_minima (array_like) – Threshold values ν for minima counting (e.g. np.linspace(-5, 1, 30)).
smoothing_scales_arcmin (tuple of float) – Gaussian FWHM smoothing scales in arcmin.
pixel_res_arcmin (float) – Pixel resolution in arcmin.
filter_size (int) – Neighbourhood size for local extremum detection.
n_jobs (int) – Number of parallel workers (joblib). 1 = serial (default). -1 = all cores. When n_jobs != 1, all smoothing scales are processed together per map, eliminating redundant smoothing passes.

Returns:

results (dict) – Nested dictionary keyed by smoothing scale (arcmin), then 'peaks' and 'minima', each containing an ndarray of shape (N, len(thresholds)).

Examples

>>> results = compute_peak_minima_counts(
...     tsz_patches, thresholds_peaks, thresholds_minima)
>>> mean_peaks = results[1.0]['peaks'].mean(axis=0)