API Reference

cmbstack.StackingPipeline

class cmbstack.StackingPipeline(sky_map, nside)

Bases: object

End-to-end stacking pipeline.

Construct from a power spectrum (from_cl()) or from an existing map (from_map()), then call run().

Parameters:

  • sky_map (numpy.ndarray) – The HEALPix map to stack on.

  • nside (int) – Resolution parameter of the map.

normalized

Set after run(); the normalized map.

Type:

numpy.ndarray or None

positions

Set after run(); the peak positions in (theta, phi) radians.

Type:

numpy.ndarray or None

patches

Set after run(); the extracted patches.

Type:

numpy.ndarray or None

stacked

Set after run(); the stacked image.

Type:

numpy.ndarray or None

radius

Set after run(); the radial distances of the profile in arcmin.

Type:

numpy.ndarray or None

profile

Set after run(); the radial profile of the stacked image.

Type:

numpy.ndarray or None

classmethod from_cl(cl_path, nside=128, seed=None)[source]

Build a pipeline by simulating a map from a power-spectrum file.

classmethod from_fits(path, field=0)[source]

Build a pipeline from any HEALPix FITS file on disk.

classmethod from_map(sky_map)[source]

Build a pipeline from a HEALPix map array already in memory. nside is inferred from the map length, so the caller doesn’t have to supply it.

run(size_deg=10.0, reso_arcmin=3.0, profile=True, threshold=3.0, n_peaks=None)[source]

Run the full stacking loop.

Parameters:

  • size_deg (float) – Patch geometry.

  • reso_arcmin (float) – Patch geometry.

  • profile (bool) – Whether to also compute the radial profile.

  • threshold (float) – Peak-finding threshold in units of the map std.

  • n_peaks (int or None) – Maximum number of peaks to use; None means use all.

Returns:

self – The pipeline object itself, with all intermediate products stored as attributes.

Return type:

StackingPipeline

cmbstack.maps

Map simulation and preprocessing utilities for HEALPix CMB maps.

Provides the building blocks for turning a raw power spectrum into a normalised HEALPix map ready for stacking:

  1. load_cl() — read a D_ell spectrum file and convert to C_ell

  2. simulate_map() — draw a Gaussian random realisation with healpy.synfast

  3. normalize_map()— subtract the monopole and divide by the std

These functions are intentionally field-agnostic: they work on any scalar power spectrum (temperature TT, lensing convergence, y-map, …).

cmbstack.maps.dl_to_cl(ell, dl, lmax)[source]

Convert D_ell = ell(ell+1) C_ell / (2 pi) to C_ell.

Parameters:

  • ell (array_like) – Multipole values. May start at 0; the ell=0 and ell=1 entries are set to zero in the output to avoid division by zero (they carry no usable power for this purpose).

  • dl (array_like) – D_ell values in the same units you want C_ell returned in (e.g. uK^2).

Returns:

cl – The angular power spectrum C_ell, same shape as dl.

Return type:

numpy.ndarray

Notes

The inverse normalization is C_ell = D_ell * 2*pi / (ell*(ell+1)).

cmbstack.maps.load_cl(path)[source]

Load a power-spectrum file and return C_ell for the requested spectrum.

The expected file columns are ell, Dl_TT, Dl_TE, Dl_EE, Dl_BB, Dl_dd, with D_ell in uK^2. The chosen column is converted from D_ell to C_ell via dl_to_cl() before being returned.

Parameters:

  • path (str) – Path to the whitespace-delimited spectrum file.

  • column (str, optional) – Which spectrum to return: one of “TT”, “TE”, “EE”, “BB”, “dd”. Default “TT”.

Returns:

cl – C_ell array indexed from ell=0, suitable for passing to simulate_map().

Return type:

numpy.ndarray

cmbstack.maps.load_map(path, field=0)[source]

Wraps hp.read_map

cmbstack.maps.normalize_map(m)[source]

Subtract the monopole and divide by the standard deviation.

After this, peak thresholds can be expressed in units of sigma, which is the natural convention for peak statistics.

Parameters:

m (array_like) – Input HEALPix map. May contain UNSEEN/NaN pixels, which are ignored in the mean and standard deviation.

Returns:

m_norm – The normalized map, with mean ~0 and std ~1.

Return type:

numpy.ndarray

cmbstack.maps.save_map(path, sky_map, overwrite=True)[source]

Wraps hp.write_map

cmbstack.maps.simulate_map(cl, nside=128, seed=None)[source]

Simulate a Gaussian random HEALPix map from a power spectrum.

Thin wrapper over healpy.synfast with an optional seed so results are reproducible in tests.

Parameters:

  • cl (array_like) – Angular power spectrum C_ell (not D_ell).

  • nside (int, optional) – HEALPix resolution parameter. Default 128.

  • seed (int or None, optional) – Seed for the random number generator. If None, the draw is random.

Returns:

m – A HEALPix map (RING ordering) of length 12 * nside**2.

Return type:

numpy.ndarray

cmbstack.stacking

Stacking of patches around positions on a HEALPix map.

The pipeline is field-agnostic: it operates on any scalar HEALPix map (CMB temperature, lensing convergence, a y-map, galaxy density, …). Positions to stack on can either be auto-detected peaks (local maxima) or supplied as an external catalogue, so the same code serves peak stacking and stacking-on-catalogue (clusters, voids, filaments, …).

Typical use

>>> peaks = find_peaks(m, nside, threshold=3.0)      # positions in (theta, phi)
>>> patches = extract_patches(m, peaks)              # fixed-grid 2D cutouts
>>> stacked = stack_patches(patches)                 # mean 2D image
>>> r, profile = radial_profile(stacked, reso_arcmin=3.0)
cmbstack.stacking.extract_patches(sky_map, positions, size_deg=10.0, reso_arcmin=3.0)[source]

Extract fixed-grid gnomonic patches centred on each position.

Each patch is a square 2D array produced by a gnomonic (tangent-plane) projection centred on the position, so every patch shares the same grid and the centre pixel always corresponds to the position itself.

Parameters:

  • sky_map (numpy.ndarray) – Input HEALPix map (RING ordering). Any scalar field.

  • positions (array_like, shape (N, 2)) – Sky positions as (theta, phi) in radians (e.g. the output of find_peaks(), or an external catalogue converted to this format).

  • size_deg (float, optional) – Full side length of the square patch in degrees. Default 10.0.

  • reso_arcmin (float, optional) – Pixel size of the projected patch in arcminutes. Default 3.0.

Returns:

patches – One square 2D array per position, all of identical shape (xsize, xsize) with xsize = size_deg * 60 / reso_arcmin.

Return type:

list of numpy.ndarray

cmbstack.stacking.find_peaks(sky_map, nside, threshold=None, n_peaks=None)[source]

Find local maxima of a HEALPix map and return their sky positions.

A pixel is a local maximum if its value is strictly greater than all of its immediate HEALPix neighbours. Peaks can be filtered by a significance threshold and/or capped at the n_peaks highest.

Parameters:

  • sky_map (numpy.ndarray) – Input HEALPix map (RING ordering). For a normalised map, values are in units of sigma, so threshold is a significance nu.

  • nside (int) – HEALPix resolution parameter of sky_map.

  • threshold (float, optional) – If given, keep only peaks with value greater than this (e.g. 3.0 for 3-sigma peaks on a normalised map). Default None (no threshold).

  • n_peaks (int, optional) – If given, keep only the n_peaks highest peaks (applied after the threshold). Default None (keep all).

Returns:

positions – Sky positions of the selected peaks as (theta, phi) in radians, the same format accepted by extract_patches().

Return type:

numpy.ndarray, shape (N, 2)

cmbstack.stacking.radial_profile(stacked, reso_arcmin=3.0, n_bins=None)[source]

Azimuthally average a stacked patch into a 1D radial profile.

Collapses the 2D stacked image to value-versus-radius by averaging in concentric annuli about the centre. This 1D profile is the characterisation of the mean peak: a central maximum, and (for CMB temperature) a faint acoustic ring further out.

Parameters:

  • stacked (numpy.ndarray) – Square 2D stacked patch from stack_patches().

  • reso_arcmin (float, optional) – Pixel size in arcminutes, so the returned radius is in physical angular units. Default 3.0.

  • n_bins (int, optional) – Number of radial bins. Default is half the patch side length.

Returns:

  • radius_arcmin (numpy.ndarray) – Bin-centre radii in arcminutes.

  • profile (numpy.ndarray) – Mean value in each annulus.

cmbstack.stacking.stack_patches(patches)[source]

Average patches pixel-by-pixel into a single stacked image.

Because every patch shares the same fixed grid (see extract_patches()), the mean is a genuine stacked image: incoherent noise averages towards zero while the coherent central profile survives. NaN/UNSEEN pixels (patch edges that fall off the map near the poles) are ignored.

Parameters:

patches (sequence of numpy.ndarray) – Patches of identical shape, e.g. the output of extract_patches().

Returns:

stacked – The mean 2D patch. Display with plt.imshow(stacked).

Return type:

numpy.ndarray