Reference

Core functions

class pywiggle.core.Wiggle(lmax, bin_edges=None, verbose=True, xlmax=2, xgllmax=2)

Position-space Power Spectrum Mode-Decoupler

This class provides the numerical backend for computing mode-decoupled angular power spectra using Gauss-Legendre quadrature methods. It is particularly suited for pseudo-spectrum estimation in the presence of incomplete sky coverage (e.g., masked observations in CMB and LSS analyses).

param lmax:

Maximum multipole to consider in the analysis.

type lmax:

int

param bin_edges:

Array of bin edges in multipole space for binning power spectra. If provided, the spectrum will be binned accordingly. Binned calculations are significantly faster. Note, these are of the form low_edge <= ℓ < upper_edge.

type bin_edges:

array_like, optional

verbosebool, default=True

Whether to print verbose information during computation.

xlmaxint, default=2

Controls how far in multipole space to use mask pseudo-spectra. Must be ≥2 for accurate decoupling, i.e. mask spectra need to provided for at least 2 x lmax.

xgllmaxint, default=2

Multiplier controlling the number of Gauss-Legendre quadrature points used, relative to lmax. Should typically be ≥2 for accurate integrals.

Notes

• The object precomputes the Legendre matrix used for spin-0 correlations (cd00), as well as a view truncated to lmax (ud00).

• The internal memory use increases with the number of masks or beams added.

• All internal 2D matrices (e.g., cd00) are stored in (theta, ell) order.

get_coupling_matrix_from_ids(mask_id1, mask_id2, spintype, bin_weight_id=None, beam_id1=None, beam_id2=None, pure_E=False, pure_B=False)

Compute the mode-coupling matrix for a given pair of masks and spin configuration.

This method retrieves the mode-coupling matrix used in power spectrum estimation, based on the specified spin configuration (spintype) and mask identifiers. The method supports scalar (spin-0) and spin-2 field combinations, e.g. for use in CMB or large-scale structure analyses.

Parameters:

• mask_id1 (str or int) – Identifier for the first mask to use in the coupling matrix calculation.

• mask_id2 (str or int) – Identifier for the second mask to use in the coupling matrix calculation.

• spintype (str) –

  The spin configuration to compute the coupling matrix for. Supported values are:

  • ’00’ : scalar-scalar (spin-0 × spin-0)

  • ’+’ : spin-2 auto (E-mode like)

  • ’-’ : spin-2 cross (B-mode like)

  • ’20’ : spin-2 × spin-0 mixed term

• bin_weight_id (str or int, optional) – Identifier for the binning weights, if applicable.

• beam_id1 (str or int, optional) – Identifier for the beam function in the first map.

• beam_id2 (str or int, optional) – Identifier for the beam function in the second map.

Returns:

The computed coupling matrix corresponding to the specified spin type and masks.

Return type:

numpy.ndarray

Raises:

ValueError – If an unsupported spintype is provided.

get_coupling_matrix_from_mask_cls(mask_cls, spintype, bin_weight_id=None, beam_id1=None, beam_id2=None, pure_E=False, pure_B=False)

Compute the mode-coupling matrix for a given mask spectrum and spin configuration.

This method retrieves the mode-coupling matrix used in power spectrum estimation, based on the specified spin configuration (spintype). The method supports scalar (spin-0) and spin-2 field combinations, e.g. for use in CMB or large-scale structure analyses.

Parameters:

• mask_cls (array_like (nells,) or (nmask,nells)) – Pseudo-Cl power spectrum of the mask, covering multipoles up to at least 2 * lmax, starting at 0. This should be precomputed externally.

• spintype (str) –

  The spin configuration to compute the coupling matrix for. Supported values are:

  • ’00’ : scalar-scalar (spin-0 × spin-0)

  • ’+’ : spin-2 auto (E-mode like)

  • ’-’ : spin-2 cross (B-mode like)

  • ’20’ : spin-2 × spin-0 mixed term

• bin_weight_id (str or int, optional) – Identifier for the binning weights, if applicable.

• beam_id1 (str or int, optional) – Identifier for the beam function in the first map.

• beam_id2 (str or int, optional) – Identifier for the beam function in the second map.

Returns:

The computed coupling matrix corresponding to the specified spin type and masks.

Return type:

numpy.ndarray

Raises:

ValueError – If an unsupported spintype is provided.

get_theory_filter(mask_id1, mask_id2=None, spintype='TT', bin_weight_id=None, pure_E=False, pure_B=False)

Construct the theoretical bandpower filter \(\mathcal{F}^{s_as_b}_{q\ell}\), as defined in arxiv:1809.09603.

This method computes the filter matrix that transforms the theoretical full-sky power spectrum \(\mathsf{C}^{ab,\mathrm{th}}_\ell\) into the corresponding prediction for the binned, decoupled bandpowers \(\mathsf{B}^{ab,\mathrm{th}}_q\) in the presence of mode coupling and nontrivial binning. The filter is given by (see arxiv:1809.09603) :

\[\mathrm{vec}\left[\mathsf{B}^{ab,\mathrm{th}}_q\right] = \sum_{\ell} \mathcal{F}^{s_as_b}_{q\ell} \cdot \mathrm{vec}\left[\mathsf{C}^{ab,\mathrm{th}}_\ell\right],\]

where:

\[\mathcal{F}^{s_as_b}_{q\ell} = \sum_{q'} \left(\mathcal{M}^{s_as_b}\right)^{-1}_{qq'} \sum_{\ell' \in \vec{\ell}_{q'}} w_{q'}^{\ell'} \mathsf{M}^{s_as_b}_{\ell'\ell}.\]

The theory spectrum this filter is dotted with must not contain the beam, even if beams were deconvolved during the Wiggle analysis.

Parameters:

• mask_id1 (str) – Identifier for the first mask used in computing the coupling matrices, previously provided through self.add_mask.

• mask_id2 (str or None, optional) – Identifier for the second mask (e.g., in cross-spectra), previously provided through self.add_mask. If None, defaults to mask_id1.

• spintype (str, default='TT') – Specifies the spin combination of the fields: - ‘TT’ for scalar × scalar (e.g., temperature or κ) - ‘+’, ‘-’ for E/B-mode combinations (spin-2 × spin-2) - ‘20’ for spin-2 × spin-0 cross spectra (e.g., shear × κ)

• bin_weight_id (str or None, optional) – ID of the binning weights to use. If not provided, defaults to uniform binning.

Returns:

thfilt – A matrix of shape (n_bins, lmax + 1) representing the filter \(\mathcal{F}^{s_as_b}_{q\ell}\) to apply to theory spectra for direct comparison with bandpower estimates.

Return type:

ndarray

add_mask(mask_id, mask_alms=None, mask_cls=None)

Register a mask for use in mode-coupling and decoupling calculations.

This method adds a sky mask to the internal cache, either in spherical harmonic (alm) form or as a precomputed pseudo-power spectrum (Cl). The mask is used to compute mode-coupling matrices that correct for incomplete sky coverage.

Parameters:

• mask_id (str) – A unique string identifier for the mask. This ID will be used in subsequent calls that reference masks for pseudo-Cl computation or coupling matrix generation.

• mask_alms (array_like, optional) – Spherical harmonic coefficients (1D array) of the mask map. Must have sufficient resolution, i.e., support at least xlmax * lmax in multipole space. Required if mask_cls is not provided.

• mask_cls (array_like, optional) – Pseudo-Cl spectrum of the mask, used as a shortcut to avoid computing mask_alms. Must cover multipoles up to xlmax * lmax. Cannot be provided at the same time as mask_alms.

Raises:

ValueError –

• If both mask_alms and mask_cls are provided. - If mask_alms is multidimensional (should be a 1D array). - If the resolution of mask_alms or mask_cls is insufficient for the configured xlmax * lmax.

Notes

• If mask_cls is provided, it is stored directly and the harmonic coefficients are not needed.

• If mask_alms is provided, its lmax is checked against the required cutoff for accurate mode-coupling computation.

• The mask is cached internally using the specified mask_id and can be reused in auto- and cross-spectrum computations.

add_bin_weights(weight_id, bin_weights)

Register custom binning weights for multipole binning.

This method allows the user to provide a weighting scheme when projecting multipole spectra into bandpowers. This is commonly used to apply inverse-variance weighting or to mimic a specific theoretical spectrum shape during the binning operation. The weights are normalized within each bin and stored under a user-defined ID.

Parameters:

• weight_id (str) – A unique string identifier for the binning weights. Must not use reserved internal IDs.

• bin_weights (array_like or None) – An array of weights of shape (lmax + 1,). Each element corresponds to a weight for the multipole \(\ell\), starting at 0. If None is passed, the method falls back to using uniform weights within each bin.

Raises:

ValueError – If weight_id is reserved for internal use. If bin_weights is provided but does not cover at least up to lmax.

Notes

• The weights are automatically normalized within each bin.

• Only the first lmax + 1 values are used.

• Once registered, the weights can be referred to by their weight_id in calls to power spectrum estimation methods.

add_beam(beam_id, beam_fl)

Register a beam transfer function for use in power spectrum beam deconvolution.

This method adds a 1D multiplicative beam transfer function \(B_\ell\) that can be deconvolved from power spectra.

Parameters:

• beam_id (str) – A unique string identifier for the beam. This ID will be referenced in calls to decoupling or filtering methods that support beam correction.

• beam_fl (array_like) – A 1D array of shape (lmax + 1,) or larger, specifying the multiplicative filter to apply to each multipole \(\ell\), starting at zero. Only the first lmax + 1 values will be retained.

Raises:

ValueError – If the length of beam_fl is less than lmax + 1, indicating insufficient multipole support for the configured maximum multipole.

Notes

• Multiple beams can be registered simultaneously under different IDs and used in cross-spectrum configurations.

get_powers(alms1, alms2, mask_ids1, mask_ids2=None, bin_weight_id=None, beam_id1=None, beam_id2=None, pure_E=False, pure_B=False, return_theory_filter=False)

Compute decoupled angular power spectra (\(C_{\ell}\)) from the spherical harmonics of maps that have already been masked.

This method estimates the angular power spectrum between two input fields in harmonic space (alm`s), accounting for mode coupling due to partial sky coverage, beam smoothing, and bandpower binning. The result is a debiased, decoupled bandpower estimate suitable for direct comparison with theoretical predictions. Note that a theory filter of shape (nbins,nells) needs to be applied to (nells,) shaped theory spectra (no beam) if using bandpower binning. This filter can be obtained from this function call, but is the most expensive part of the calculation, so consider obtaining it from `self.get_theory_filter once.

Parameters:

• alms1 (ndarray) –

  Spherical-harmonic coefficients of the first and second map, both with shape

  • (N_alm,) – a single scalar field;

  • (1, N_alm) – same as above (explicit 1-pol);

  • (2, N_alm) – spin-2 field (E, B);

  • (3, N_alm) – scalar + spin-2 (T, E, B).

• alms2 (ndarray) –

  Spherical-harmonic coefficients of the first and second map, both with shape

  • (N_alm,) – a single scalar field;

  • (1, N_alm) – same as above (explicit 1-pol);

  • (2, N_alm) – spin-2 field (E, B);

  • (3, N_alm) – scalar + spin-2 (T, E, B).

• the (The two arrays must be identical in shape. If you pass)

• auto-spectra. (*same* object twice the function returns)

mask_ids1str or sequence(str)

Mask identifier(s) previously registered via add_mask() for the first field. If you pass a single string it is applied to all polarisation components; otherwise give a two-element list/tuple with one entry for the spin-0 field and the second for the spin-2 field.

mask_ids2str or sequence(str), optional

Mask identifier(s) for the second field. If None (default) the same mask(s) as mask_ids1 are used.

bin_weight_idstr or None, optional

Identifier for custom multipole binning weights. If not specified, unity weights will be used by default.

beam_id1str or None, optional

Beam ID for the first map to deconvolve a beam previously provided through self.add_beam. If None, no beam is deconvolved for the first map.

beam_id2str or None, optional

Beam ID for the second map to deconvolve a beam previously provided through self.add_beam. If None, no beam is deconvolved for the second map.

return_theory_filterbool, default=False

If True, also return the theory bandpower filter \(\mathcal{F}^{s_as_b}_{q\ell}\) for use in model comparison.

• scalar–scalar -> {'TT': {'Cls': (nbin | lmax+1,)}}

• spin-2 only -> {'EE', 'EB', 'BE', 'BB'}

• scalar+spin-2 -> ``{‘TT’, ‘TE’, ‘ET’, ‘TB’, ‘BT’,

  ‘EE’, ‘EB’, ‘BE’, ‘BB’}``

Each entry holds a dict with the field 'Cls' containing the decoupled (and possibly binned) spectrum.

If return_theory_filter is True, a second key 'Th' is added containing the corresponding theory filters with shape (n_bins, lmax+1).

For the pure-polarization part, the theory filter combines EE, EB, BE and BB and is in spectra[‘ThPol’]. For most cases involving null EB, BE and BB, you only need the EE portion of this matrix.

• ———-

Raises:

ValueError –

• The two alm arrays differ in shape. * Unsupported number of polarisation components. * Inconsistent mask list length. * Internal binning/mixing dimensions do not match.

Notes

• This method is the recommended interface for evaluating both auto- and cross-spectra in masked sky analyses, but other convenience wrappers are provided elsewhere.

decoupled_cl(pcls, mask_id1, mask_id2=None, spintype='TT', bin_weight_id=None, beam_id1=None, beam_id2=None, return_theory_filter=False, pure_E=False, pure_B=False)

pywiggle.core.get_binned_theory(theory_filter_dict, cl_dict)

Compute binned theoretical power spectra using provided filter matrices and input C_ells.

This function takes in a dictionary of theory filter matrices (typically returned by Wiggle functions) and a dictionary of angular power spectra (cl_dict), and returns a dictionary of binned theoretical spectra after matrix multiplication with the appropriate filters.

Supports scalar spectra (TT, TE, TB) and polarization spectra (EE, EB, BE, BB) via a block matrix labeled ThPol.

Parameters:

• theory_filter_dict (dict) –

  Dictionary containing filter matrices:

  • Keys may include ‘TT’, ‘TE’, ‘TB’, each mapping to a dict with key ‘Th’ (2D array) each of shape (nbins, lmax).

  • May also include ‘ThPol’, a 2D array of shape (4 * nbins, lmax) for polarization.

• cl_dict (dict) – Dictionary mapping spectrum keys (‘TT’, ‘TE’, ‘TB’, ‘EE’, ‘BB’, ‘EB’) to 1D numpy arrays of C_ell values, starting at ell=0.

Returns:

ret_th – Dictionary containing the binned theory spectra with keys corresponding to input types: - ‘TT’, ‘TE’, ‘TB’ (if present in inputs) - ‘EE’, ‘EB’, ‘BE’, ‘BB’ (if ‘ThPol’ is provided)

Return type:

dict

Notes

• Missing TB, BB, or EB spectra in cl_dict are assumed to be zero with a warning.

• The maximum multipole lmax is inferred from the second dimension of each ‘Th’ matrix.

pywiggle.core.get_coupling_matrix_from_mask_cls(mask_cls, lmax, spintype='TT', bin_edges=None, bin_weights=None, beam_fl1=None, beam_fl2=None, pure_E=False, pure_B=False, return_obj=False, verbose=True, xlmax=2, xgllmax=2)

Compute the (optionally, binned) mode-coupling matrix from the pseudo-Cl of a sky mask.

This function is a high-level wrapper to generate the (optionally, binned) mode-coupling matrix \(\mathcal{M}^{s_as_b}_{q q'}\) using the pseudo-power spectrum of a mask. This matrix describes how true sky power at one multipole leaks into others due to incomplete sky coverage, beam smoothing, and binning.

Parameters:

• mask_cls (array_like (nells,) or (nmask,nells)) – Pseudo-Cl power spectrum of the mask, covering multipoles up to at least 2 * lmax, starting at 0. This should be precomputed externally.

• lmax (int) – Maximum multipole for the output coupling matrix.

• spintype (str, default='TT', or list of spins.) – Spin combination of the fields: - ‘TT’: scalar × scalar (e.g., T × T or κ × κ) - ‘+’, ‘-’: spin-2 × spin-2 (e.g., E × E or B × B) - ‘TE’: spin-2 E × spin-0 (e.g., E × κ or γ_E × κ) - ‘TB’: spin-2 E × spin-0 (e.g., B × κ or γ_B × κ)

• bin_edges (array_like, optional) – Array of bin edges defining bandpowers. If not provided, no binning is applied. Note, these are of the form low_edge <= ℓ < upper_edge.

• bin_weights (array_like, optional) – Weights for each multipole used in binning. Must have at least lmax + 1 entries, starting at 0. If not provided, uniform weights will be assumed.

• beam_fl1 (array_like, optional) – Beam transfer function for the first field (length ≥ lmax + 1). Optional.

• beam_fl2 (array_like, optional) – Beam transfer function for the second field. Optional.

• return_obj (bool, default=False) – If True, also return the internal Wiggle object for further manipulation.

• xlmax (int, default=2) – Controls how far in multipole space to use mask pseudo-spectra. Must be ≥2 for accurate decoupling, i.e. mask spectra need to provided for at least 2 x lmax.

• xgllmax (int, default=2) – Multiplier controlling the number of Gauss-Legendre quadrature points used, relative to lmax. Should typically be ≥2 for accurate integrals.

Returns:

• m (ndarray) – The (optionally, binned) mode-coupling matrix of shape (nmask, nspin, n_bins, n_bins) if binned or (nmask, nspin, lmax+1, lmax+1) if unbinned.

• g (Wiggle, optional) – The Wiggle object used to generate the matrix, returned only if return_obj=True.

pywiggle.core.get_pure_EB_alms(Qmap, Umap, mask, masked_on_input=False, return_mask=False, lmax=None, eps=0.0001)

Perform pure E/B mode decomposition on a masked polarization map. Based on code from Irene Abril-Cabezas and Will Coulton.

Parameters:

• Qmap (ndarray or enmap) – Stokes Q polarization map. If masked_on_input=True, this is assumed to have already been multiplied by mask.

• Umap (ndarray or enmap) – Stokes U polarization map. Same geometry as Qmap.

• mask (ndarray or enmap) – Scalar apodized mask in the same geometry as Qmap and Umap.

• masked_on_input (bool, optional) – If True, assumes Q and U have already been multiplied by mask and divides it out (with threshold protection) internally. Default is False.

• return_mask (bool, optional) – If True, returns the spin-1 and spin-2 derivatives of the mask used in the purification, instead of the pure E/B multipoles. Default is False.

• lmax (int or None, optional) – Maximum multipole to compute. If None, it defaults to pi/pixel_size/2 based on the input map geometry.

• eps (float, optional) – Threshold below which mask values are treated as zero when dividing out the mask. Default is 1e-4.

Returns:

(pureE_alms, pureB_alms) – Spherical harmonic coefficients of the purified E- and B-mode signals. If return_mask=True, returns (mask1, mask2) instead, where mask1 and mask2 are the spin-1 and spin-2 derivatives of the mask.

Return type:

tuple of ndarrays

Notes

• Requires the pixell library for spherical harmonic transforms and map handling.

• Implements the Smith & Zaldarriaga (2007) pure E/B formalism to suppress E->B leakage induced by masking.

class pywiggle.core.EBPurifier(mask, lmax, eps=0.0001)

Precomputes mask-dependent operators for pure E/B projection.

This class builds the spherical-harmonic mask derivatives and l-space weights once. Subsequent calls to project only perform the per-map transforms.

Parameters:

• mask (ndarray or enmap) – Scalar apodized mask, same geometry as the input maps.

• lmax (int) – Maximum multipole.

• eps (float) – Threshold for mask safety when dividing.

project(Qmap, Umap, masked_on_input=False)

Compute pure E/B multipoles for the provided Q/U maps using the precomputed mask operators.

Parameters:

• Qmap (ndarray or enmap) – Input Stokes Q/U maps. If masked_on_input=True, they are assumed to have already been multiplied by mask.

• Umap (ndarray or enmap) – Input Stokes Q/U maps. If masked_on_input=True, they are assumed to have already been multiplied by mask.

• masked_on_input (bool, optional) – If True, divides out mask (with threshold protection) before building pure estimators. Default is False.

Returns:

pureE_alms, pureB_alms – Purified E- and B-mode alm arrays.

Return type:

ndarrays

pywiggle.core.get_powers(alms1, alms2, mask_alm1, mask_alm2=None, bin_weights=None, beam_fl1=None, beam_fl2=None, pure_E=False, pure_B=False, return_theory_filter=False, lmax=None, bin_edges=None, verbose=False, xlmax=2, xgllmax=2)

Compute decoupled angular power spectra (\(C_{\ell}\)) from the spherical harmonics of maps that have already been masked.

This method estimates the angular power spectrum between two input fields in harmonic space (alm`s), accounting for mode coupling due to partial sky coverage, beam smoothing, and bandpower binning. The result is a debiased, decoupled bandpower estimate suitable for direct comparison with theoretical predictions. Note that a theory filter of shape (nbins,nells) needs to be applied to (nells,) shaped theory spectra (no beam) if using bandpower binning. This filter can be obtained from this function call, but is the most expensive part of the calculation, so consider obtaining it from `self.get_theory_filter once.

Parameters:

• alms1 (ndarray) –

  Spherical-harmonic coefficients of the first and second map, both with shape

  • (N_alm,) – a single scalar (spin-0) field;

  • (1, N_alm) – same as above (explicit 1-pol);

  • (2, N_alm) – spin-2 field (E, B);

  • (3, N_alm) – scalar + spin-2 (T, E, B).

• alms2 (ndarray) –

  Spherical-harmonic coefficients of the first and second map, both with shape

  • (N_alm,) – a single scalar (spin-0) field;

  • (1, N_alm) – same as above (explicit 1-pol);

  • (2, N_alm) – spin-2 field (E, B);

  • (3, N_alm) – scalar + spin-2 (T, E, B).

• the (The two arrays must be identical in shape. If you pass)

• auto-spectra. (*same* object twice the function returns)

• mask_alm1 (array_like, optional) – Spherical harmonic coefficients of the mask that was applied. Must have sufficient resolution, i.e., support at least xlmax * lmax in multipole space, where xlmax is 2 by default. If 1d, or (1,nalm) shaped, the same mask is used for spin-0 and spin-2. If (2,nalm) shaped, the first mask is used for the spin-0 field and the second for the spin-2 fields.

• mask_alm2 (array_like, optional) – Second mask in alm space. If None, uses mask_alm1 for both fields.

• bin_weights (array_like, optional) – Weights to apply when binning the resulting power spectra.

• beam_fl1 (array_like, optional) – Beam transfer function for the first field.

• beam_fl2 (array_like, optional) – Beam transfer function for the second field.

• pure_E (bool, default=False) – If True, apply pure-E mode estimation.

• pure_B (bool, default=False) – If True, apply pure-B mode estimation.

• return_theory_filter (bool, default=False) – If True, also return the theory bandpower filter \(\mathcal{F}^{s_as_b}_{q\ell}\) for use in model comparison.

• lmax (int) – Maximum multipole to consider in the analysis.

• bin_edges (array_like, optional) – Array of bin edges in multipole space for binning power spectra. If provided, the spectrum will be binned accordingly. Binned calculations are significantly faster. Note, these are of the form low_edge <= ℓ < upper_edge.

verbosebool, default=True

Whether to print verbose information during computation.

xlmaxint, default=2

Controls how far in multipole space to use mask pseudo-spectra. Must be ≥2 for accurate decoupling, i.e. mask spectra need to provided for at least 2 x lmax.

xgllmaxint, default=2

Multiplier controlling the number of Gauss-Legendre quadrature points used, relative to lmax. Should typically be ≥2 for accurate integrals.

Returns:

spectra – Keys and shapes depend on the input polarisation:

• scalar–scalar -> {'TT': {'Cls': (nbin | lmax+1,)}}

• spin-2 only -> {'EE', 'EB', 'BE', 'BB'}

• scalar+spin-2 -> ``{‘TT’, ‘TE’, ‘ET’, ‘TB’, ‘BT’,

  ’EE’, ‘EB’, ‘BE’, ‘BB’}``

Each entry holds a dict with the field 'Cls' containing the decoupled (and possibly binned) spectrum.

If return_theory_filter is True, a second key 'Th' is added containing the corresponding theory filters with shape (n_bins, lmax+1).

For the pure-polarization part, the theory filter combines EE, EB, BE and BB and is in spectra[‘ThPol’]. For most cases involving null EB, BE and BB, you only need the EE portion of this matrix.

Return type:

dict

Raises:

ValueError –

• If input alms have mismatched shapes or unsupported dimensionality. * Unsupported number of polarisation components. * Inconsistent mask list length. * Internal binning/mixing dimensions do not match.

Utilities

pywiggle.utils.change_alm_lmax(alms, lmax, mmax_out=None)

pywiggle.utils.multipole_to_bin_indices(lmax, bin_edges)

Returns an array mapping multipole index l (0 to lmax) to bin index, based on provided bin_edges.

Parameters: - lmax: int, maximum multipole - bin_edges: 1D array-like, sorted bin edges of length nbins + 1

Returns: - bin_indices: np.ndarray of shape (lmax + 1,), where bin_indices[l] gives the bin index for multipole l

pywiggle.utils.bin_array(arr, bin_indices, nbins, weights=None)

Bin an array using bin indices and optional per-element weights.

Parameters:

• arr (ndarray of shape (lmax + 1,)) – Array defined over multipole ℓ values.

• bin_indices (ndarray of shape (lmax + 1,)) – Per-ℓ bin index (-1 for entries to ignore).

• weights (ndarray of shape (lmax + 1,), optional) – Optional per-element weights w_ell. If not provided, defaults to 1.

Returns:

binned – Array with weighted sum of values in each bin.

Return type:

ndarray of shape (nbins,)

pywiggle.utils.normalize_weights_per_bin(nbins, bin_indices, weights)

Normalize weights so that values in each bin sum to 1.

Parameters:

• nbins (int) – Number of bins.

• bin_indices (ndarray of shape (N,)) – Bin index for each element. Use -1 for elements to be ignored.

• weights (ndarray of shape (N,)) – Weights to normalize.

Returns:

normalized_weights – Weights normalized per bin.

Return type:

ndarray of shape (N,)

pywiggle.utils.get_wigner_d(lmax, s1, s2, cos_thetas, bin_edges=None, bin_weights=None, bin_weights2=None)

Compute Wigner-d matrices for given spins and angles, with optional weighted binning over multipoles.

This function evaluates the Wigner small-d matrices, ( d^{ell}_{s_1,s_2}(theta) ), for multipoles from 0 up to lmax at the specified cosines of angles (cos_thetas). It supports optional binning of the multipole axis using specified bin edges and weighting schemes.

Parameters:

• lmax (int) – Maximum multipole (ℓ) for which to compute the Wigner-d matrices.

• s1 (int) – First spin index for the Wigner-d matrix.

• s2 (int) – Second spin index for the Wigner-d matrix.

• cos_thetas (array_like) – Array of cosines of the angles θ at which to evaluate the Wigner-d matrices.

• bin_edges (array_like, optional) – Bin edges defining multipole bins for aggregating the Wigner-d matrices. Bins are interpreted as intervals: low_edge ≤ ℓ < upper_edge. If not provided, no binning is performed and the full multipole range is returned.

• bin_weights (array_like, optional) – Multipole weights used for binning along the ℓ axis. If provided, enables weighted binning. If bin_weights2 is not provided, single weighted binning is performed.

• bin_weights2 (array_like, optional) – Optional second set of multipole weights for double-weighted binning. Only used if bin_weights is also provided. If bin_weights2 is provided, double weighted binning is performed.

Returns:

If no binning is requested (bin_edges is None), returns a 2D array of shape (len(cos_thetas), lmax + 1) containing the Wigner-d matrices evaluated at each angle and multipole.

If binning is requested, returns a binned (or double-binned) Wigner-d matrix array with shape depending on the number of bins and angles.

Return type:

ndarray

Raises:

ValueError – If bin_weights is provided without bin_edges.

Notes

• This function is intended primarily for testing purposes and is not used in wiggle.

pywiggle.utils.bin_square_matrix(matrix, bin_edges, lmax, bin_weights=None, bin_weights2=None)

Bin a square matrix in along both directions using specified multipole binning and weights.

This function reduces a 2D square matrix (e.g., a coupling matrix) from (lmax+1) × (lmax+1) to a binned form of shape (n_bins × n_bins), where binning is defined by multipole bin_edges and optional weights. It supports different weighting schemes along each axis of the matrix.

Parameters:

• matrix (2D array_like) – Square matrix of shape (lmax+1, lmax+1) to be binned. Typically a function of multipole indices (ℓ, ℓ’).

• bin_edges (array_like) – Array defining the edges of the multipole bins. Must span the relevant ℓ-range up to lmax. Note, these are of the form low_edge <= ℓ < upper_edge.

• lmax (int) – Maximum multipole to consider. Defines the size of the unbinned matrix as (lmax+1) × (lmax+1).

• bin_weights (array_like or None, optional) – Weights to apply within each bin along the first axis (rows). If None, unit weights are used and internally normalized.

• bin_weights2 (array_like or None, optional) – Weights to apply within each bin along the second axis (columns). If None, unit weights are used but not normalized—useful for asymmetric binning or theoretical filters.

Returns:

binned_matrix – Binned square matrix of shape (n_bins, n_bins), where n_bins = len(bin_edges) - 1.

Return type:

ndarray

Notes

• Internally uses normalized weights along the first axis, and unnormalized unity for the second axis if weights for the second axis are not provided. This correctly transforms an unbinned coupling matrix to a binned coupling matrix.

pywiggle.utils.wfactor(n, mask, sht=True, pmap=None, equal_area=False)

Copied from msyriac/orphics/maps.py

Approximate correction to an n-point function for the loss of power due to the application of a mask.

For an n-point function using SHTs, this is the ratio of area weighted by the nth power of the mask to the full sky area 4 pi. This simplifies to mean(mask**n) for equal area pixelizations like healpix. For SHTs on CAR, it is sum(mask**n * pixel_area_map) / 4pi. When using FFTs, it is the area weighted by the nth power normalized to the area of the map. This also simplifies to mean(mask**n) for equal area pixels. For CAR, it is sum(mask**n * pixel_area_map) / sum(pixel_area_map).

If not, it does an expensive calculation of the map of pixel areas. If this has been pre-calculated, it can be provided as the pmap argument.

pywiggle.utils.get_camb_spectra(lmax=512, tensor=True, ns=0.965, As=2e-09, r=0.05)

Returns CMB Cls [TT, EE, BB, TE] from CAMB with low accuracy for fast testing.

pywiggle.utils.get_cldict(ps)

pywiggle.utils.cosine_apodize(bmask, width_deg)

pywiggle.utils.cprint(string, color=None, bold=False, uline=False)

class pywiggle.utils.bcolors

Colored output for print commands

HEADER = '\x1b[95m'

OKBLUE = '\x1b[94m'

OKGREEN = '\x1b[92m'

WARNING = '\x1b[93m'

FAIL = '\x1b[91m'

ENDC = '\x1b[0m'

BOLD = '\x1b[1m'

UNDERLINE = '\x1b[4m'

pywiggle.utils.hplot(img, savename=None, verbose=True, grid=False, **kwargs)

pywiggle.utils.mollview(hp_map, filename=None, lim=None, coord='C', verbose=True, return_projected_map=False, xsize=1200, grat_deg=None, dpi=None, grat_color='gray', grat_alpha=0.5, **kwargs)

mollview plot for healpix wrapper

class pywiggle.utils.Benchmark(lmax, bin_edges=None, nthreads=None, verbose=False)

get_mcm(code, spin=0, bin_edges=None, bin_weights=None)