cosmocore.basis module 

abstractmethod get_full_inverse(C_ell)[source]

Return the full n_pix x n_pix inverse covariance.

For HarmonicBasis: exact (N + S)^{-1} via the SMW formula. For PixelBasis in pixel-direct mode: also exact (U is the identity, so U @ get_inverse(C_ell) @ U.T equals the full inverse).

For PixelBasis on a truncated compressed basis: the returned matrix is U @ get_inverse(C_ell) @ U.T, which lives in the kept subspace and is zero on the truncated complement. It is not the inverse of the full pixel-space covariance; it is the inverse of the restricted operator, lifted back to n_pix dimensions. Use with care.

Callers who only need the basis-space inverse should use get_inverse() instead — cheaper and equivalent in basis-space algebra.

Return type:: ndarray

get_full_logdet(C_ell)[source]

Return log|N + S|, the full pixel-space log determinant.

HarmonicBasis overrides this to compute it exactly via the SMW formula. The default below applies to PixelBasis and is exact only in pixel-direct mode (U is the identity, so the basis-space logdet equals the full logdet).

On a truncated compressed pixel basis the full logdet cannot be recovered from the kept quantities: the discarded complement of N + S is needed and is not stored. The default falls through to get_logdet(), which returns the logdet of the restricted operator U^T (N + S) U — a different matrix, not an approximation to log|N + S|. Callers that need the exact full logdet on a truncated basis must arrange for it elsewhere.

Parameters:: C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
Return type:: float

get_inverse(C_ell)[source]

Compute the inverse of the basis-space covariance.

Parameters:: C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
Returns:: Basis-space inverse covariance of shape (dim, dim).
Return type:: numpy.ndarray

get_logdet(C_ell)[source]

Compute the log determinant of the basis-space covariance.

Parameters:: C_ell (numpy.ndarray) – Power spectrum values for ell = 2 to lmax.
Returns:: Log determinant of basis-space covariance.
Return type:: float

abstractmethod get_noise_for_bias()[source]

Basis-projected noise consumed by the QML noise-bias sandwich.

Spectra computes Cov(w | noise) = X get_noise_for_bias() X^T where X is the basis’s natural inverse-equivalent (A = (I + Lambda M)^{-T} for the harmonic basis; C_c^{-1} for the pixel basis). The form of the returned matrix is therefore basis-specific:

HarmonicBasis: V N_eff^{-1} N N_eff^{-1} V^T (the SMW intermediate T from N_eff = N + S_fixed when the switch optimisation is active; otherwise V N^{-1} N N^{-1} V^T).
PixelBasis: U^T N U (raw noise projected once into the eigenmode basis).

The two are not interchangeable. Callers must compose with the basis’s own inverse-equivalent.

Returns:: Symmetric (dim, dim) matrix.
Return type:: numpy.ndarray

abstractmethod get_projected_inverse(C_ell)[source]

Get the projected inverse for Fisher computation.

Parameters:: C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
Returns:: Projected inverse covariance matrix.
Return type:: numpy.ndarray

abstractmethod get_weighted_data(data, C_ell, C_c_inv=None, stable_inner_inv=None)[source]

Compute weighted basis-space data for QML estimation.

Parameters:

data (numpy.ndarray) – Pixel-space data vector of shape (n_pix,) or (n_pix, n_sims).
C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
C_c_inv (numpy.ndarray, optional) – Precomputed basis-space inverse. Used by PixelBasis only; ignored by HarmonicBasis.
stable_inner_inv (numpy.ndarray, optional) – Precomputed (I + Lambda M)^{-1}, as returned by HarmonicBasis.prepare_stable_inner_inv(). Used by HarmonicBasis only; ignored by PixelBasis.

Returns:

Weighted basis-space data of shape (dim,) or (dim, n_sims).

Return type:

abstract property method: str

“harmonic” or “pixel”.

Type:: Computation basis name

abstractmethod prepare_for_basis(C_ell_dict)[source]

Pre-bake per-C_ell quantities for fast likelihood reuse.

Returns a BasisPrepared carrying the basis-specific factor and the full log determinant. The factor is consumable by quadratic_form_from_prepared() on the same basis instance.

Parameters:: C_ell_dict (dict) – Multi-field power spectrum dict keyed by SpectrumKey.
Return type:: BasisPrepared

abstract property projector: ndarray

Get the projection matrix that maps pixel space to the basis.

The fundamental change-of-basis operator: - HarmonicBasis: V (n_modes × n_pix) - PixelBasis: U^T (dim × n_pix) when compression is active;

identity-equivalent in pixel-direct mode.

Returns:: Projection matrix of shape (dim, n_pix).
Return type:: numpy.ndarray

abstractmethod quadratic_form(data, C_ell)[source]

Compute d^T C^{-1} d for full pixel-space data.

The basis chooses its fastest path. Returns a single float for a 1-D data vector, or a length-n_sims array if data is shape (n_pix, n_sims).

Parameters:

data (numpy.ndarray) – Pixel-space data, shape (n_pix,) or (n_pix, n_sims).
C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).

Return type:

abstractmethod quadratic_form_from_prepared(data, factor)[source]

Compute d^T C^{-1} d using a pre-baked factor from prepare_for_basis().

Fast path for evaluating the quadratic form many times against the same C_ell: caller calls prepare_for_basis() once and reuses the resulting factor across simulations.

Parameters:

data (numpy.ndarray) – Pixel-space data, shape (n_pix,) or (n_pix, n_sims).
factor (numpy.ndarray) – Basis-specific pre-baked quantity: BasisPrepared.factor from prepare_for_basis().

Return type:

release_pixel_projector()[source]

No-op default. HarmonicBasis overrides to drop V post-SMW build.

Return type:: None

abstractmethod setup()[source]

Initialize basis-specific components.

Must be called after initialization before using any other methods.

Return type:: None

to_basis(data)[source]

Project pixel data into the basis: d_basis = P @ d.

Uses the basis-specific projector P: - HarmonicBasis: P = V (n_modes × n_pix) - PixelBasis: P = U^T (dim × n_pix) when compression is active;

identity short-circuit in pixel-direct mode.

Parameters:: data (numpy.ndarray) – Pixel-space data vector of shape (n_pix,) or (n_pix, n_sims).
Returns:: Basis-space data of shape (dim,) or (dim, n_sims).
Return type:: numpy.ndarray

class cosmocore.basis.HarmonicBasis(N, theta, phi, lmax_signal, beam=None, spins=None, lmin_signal=None, lmin=None, lmax=None, fiducial_C_ell=None, S_fixed=None, compress=False, delta_m=0)[source]

Direct harmonic space computation basis (Tegmark-like).

Transforms directly to n_modes dimensions via V. Fast and efficient when n_modes << n_pix. No additional eigenvalue truncation.

The key operations are: - Data projection: d̄ = P @ d where P = V - Projected covariance: C̄ = V N V^T + Λ - Projected inverse: V C^{-1} V^T = M - M K^{-1} M (via SMW)

Parameters:

N (numpy.ndarray) – Noise covariance matrix of shape (n_pix, n_pix). The basis takes ownership of this buffer; _factorise_noise() overwrites it in place during setup.
theta (numpy.ndarray) – Colatitude angles for active pixels in radians.
phi (numpy.ndarray) – Longitude angles for active pixels in radians.
lmax (int) – Maximum multipole for harmonic expansion.
beam (numpy.ndarray or None, optional) – Beam window function B_ℓ for ℓ=2 to lmax.
compress (bool, optional) – If True, use m-block compression for K inversion. Approximates K as block-diagonal in azimuthal quantum number |m|, giving ~lmax^2 speedup. Default is False.
delta_m (int, optional) – Bandwidth for m-block coupling. delta_m=0 means block-diagonal (no coupling between different |m|). delta_m=lmax recovers exact result. Default is 0.

Examples

>>> import numpy as np
>>> from cosmocore.basis import HarmonicBasis
>>> hc = HarmonicBasis(N, theta, phi, lmax=100)
>>> hc.setup()
>>> fisher_element = hc.compute_fisher_element(C_ell, ell_i=10, ell_j=10)

References

__init__(N, theta, phi, lmax_signal, beam=None, spins=None, lmin_signal=None, lmin=None, lmax=None, fiducial_C_ell=None, S_fixed=None, compress=False, delta_m=0)[source]

Initialize computation basis.

Parameters:

N (numpy.ndarray) – Noise covariance matrix (n_pix_total, n_pix_total). The basis takes ownership; _factorise_noise() overwrites this buffer in place during setup.
theta (numpy.ndarray or tuple of numpy.ndarray) – Colatitude angles for active pixels in radians. Single array for single-field, tuple of arrays for multi-field (one per component).
phi (numpy.ndarray or tuple of numpy.ndarray) – Longitude angles for active pixels in radians. Single array for single-field, tuple of arrays for multi-field (one per component).
lmax_signal (int) – Signal-cov ceiling. The basis represents multipoles up to and including this value.
beam (numpy.ndarray or None, optional) – Beam window function B_ℓ for ℓ=0..lmax_signal.
spins (list of int or None, optional) – Spin weight for each component (0 for scalar, 2 for polarization). Default is [0, …] for all components. For spin-2 components, theta/phi represent physical pixel locations but V is built for (Q, U) → (E, B) transformation with doubled dimensions.
lmin_signal (list of int or None, optional) – Per-component signal-cov floor. Defaults to [2]*n_components. Each entry must satisfy lmin_signal[i] >= |spins[i]|.
lmin (int or None, optional) – Inference window lower bound. Defaults to min(lmin_signal). Together with lmax enables the switch optimisation: S_fixed absorbs the signal contribution outside [lmin, lmax].
lmax (int or None, optional) – Inference window upper bound. Defaults to lmax_signal. Multipoles in (lmax, lmax_signal] are absorbed into S_fixed.
fiducial_C_ell (numpy.ndarray or None, optional) – Deprecated: Use S_fixed instead. Fiducial power spectrum for fixed multipoles (ℓ outside [lmin, lmax]).
S_fixed (numpy.ndarray or None, optional) – Precomputed signal matrix for fixed multipoles outside [lmin, lmax]. Recommended way to pass the fixed signal contribution. Shape (n_pix_total, n_pix_total).

compute_fisher_matrix(C_ell, spectra_list=None, ell_min=2, ell_max=None, *, symmetry_mode=None)[source]

Compute Fisher matrix.

Parameters:

C_ell (numpy.ndarray or dict) – Power spectrum. Single-field: numpy.ndarray (spectra_list=None). Multi-field: dict keyed by SpectrumKey (spectra_list required).
spectra_list (list[SpectrumKey] or None) – For multi-field, the list of SpectrumKey instances enumerating the spectra to include. None for single-field.
ell_min (int) – Minimum multipole.
ell_max (int or None) – Maximum multipole.
symmetry_mode (SymmetryMode or None) – Forwarded to the per-ℓ derivative builder. None keeps the SYMMETRIC default. Required as DIRECTIONAL when spectra_list contains a CG cross-pair entry — the builder otherwise rejects the mode=2 slot used by CG.

Returns:

Fisher matrix.

Return type:

get_covariance(C_ell)[source]

Compute harmonic-space covariance C̄ = V N V^T + Λ.

Parameters:: C_ell (numpy.ndarray or dict) – Power spectrum. Can be: - numpy.ndarray: C_ell values for ell = 2 to lmax (single-field) - dict: Multi-field with 2-tuple or 3-tuple keys
Returns:: Compressed covariance matrix.
Return type:: numpy.ndarray

get_derivative_diagonal(ell, key)[source]

Get derivative diagonal for an auto-pair diagonal SpectrumKey.

Returns a 1D vector — the diagonal of ∂S/∂C_ℓ — for use in the sparse-trace fast path. Valid only for key.comp_i == key.comp_j and diagonal kinds (SS, GG, CC) — i.e. TT, EE, or BB. Cross-component keys and GC/CG/SG/SC/GS/CS raise ValueError.

Return type:: ndarray

get_derivative_matrix(ell, key, *, symmetry_mode=None)[source]

Get the derivative matrix ∂S/∂C_ℓ in harmonic space.

Parameters:

ell (int) – Multipole.
key (SpectrumKey) – Identifies the spectrum (component pair + spin kind).
symmetry_mode (SymmetryMode or None) – None keeps the SYMMETRIC default.

Returns:

Derivative matrix of shape (n_modes, n_modes).

Return type:

get_full_inverse(C_ell)[source]

Compute full (N + S)^{-1} using the SMW formula.

Parameters:: C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
Returns:: Full inverse covariance matrix of shape (n_pix, n_pix).
Return type:: numpy.ndarray

get_full_logdet(C_ell)[source]

Compute log|N + S| using SMW formula.

Parameters:: C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
Returns:: Log determinant of the full covariance matrix.
Return type:: float

get_noise_for_bias()[source]

Return the SMW noise intermediate T = V N_eff^{-1} N N_eff^{-1} V^T.

See ComputationBasis.get_noise_for_bias() for the cross-basis contract. T is precomputed during basis setup.

Return type:: ndarray

get_projected_inverse(C_ell, field_groups=None)[source]

Compute V C^{-1} V^T efficiently using SMW formula.

Parameters:

C_ell (numpy.ndarray or dict) – Power spectrum. Can be array (single-field) or dict (multi-field).
field_groups (list of list of int or None, optional) – Independent field groups from _detect_field_blocks. If provided and contains more than one group, exploits field block-diagonal structure for faster K inversion.

Returns:

Projected inverse covariance V C^{-1} V^T.

Return type:

get_weighted_data(data, C_ell, C_c_inv=None, stable_inner_inv=None)[source]

Compute w = V C^{-1} d via the stable SMW algebra.

w = (I + Lambda M)^{-T} (V N^{-1} d) — algebraically equivalent to M K^{-1} d but never subtracts large nearly-equal matrices, so it retains precision in the cosmic-variance-limited regime where the legacy subtractive form y - M K^{-1} y cancelled.

Parameters:

data (numpy.ndarray) – Pixel-space data, shape (n_pix,) or (n_pix, n_sims).
C_ell (numpy.ndarray or dict) – Power spectrum (single-field array or multi-field dict).
C_c_inv (numpy.ndarray, optional) – Unused for harmonic basis; accepted for ABC signature parity.
stable_inner_inv (numpy.ndarray, optional) – Precomputed (I + Lambda M)^{-1}. If supplied, avoids rebuilding it internally. Callers that weight several data vectors with the same C_ell should build it once via prepare_stable_inner_inv() and pass it in.

Returns:

Weighted harmonic-space data w = V C^{-1} d.

Return type:

property method: str: Computation basis name.

prepare_for_basis(C_ell_dict)[source]

Precompute K Cholesky factor and log determinant for reuse across sims.

Return type:: BasisPrepared

prepare_stable_inner_inv(C_ell, lambda_matrix=None)[source]

Build (I + Lambda M)^{-1} for the stable QML algebra.

This matrix appears identically in three places of the QML algebra and replaces the unstable subtractive forms:

Projected inverse: V C^{-1} V^T = M @ (I + Lambda M)^{-1}
Data weighting: V C^{-1} d = (I + Lambda M)^{-1} @ V N^{-1} d
Noise bias matrix: A = I - M K^{-1} = (I + Lambda M)^{-1}

All three follow from the identity M - M K^{-1} M = M K^{-1} Lambda^{-1} = M (I + Lambda M)^{-1} and never subtract large nearly-equal matrices, so they remain accurate in the cosmic-variance-limited regime where the legacy SMW form loses precision (e.g. T at low ell with sub-uK polarization noise).

Parameters:

C_ell (np.ndarray or dict) – Power spectrum (single-field array or multi-field dict).
lambda_matrix (np.ndarray, optional) – Pre-built full Lambda matrix. If None, built from C_ell.

Returns:

(I + Lambda M)^{-1}, n_modes_total x n_modes_total.

Return type:

np.ndarray

property projector: ndarray: Projection matrix V (n_modes × n_pix).

quadratic_form(data, C_ell)[source]

Compute d^T C^{-1} d efficiently using SMW formula.

Parameters:

data (numpy.ndarray) – Pixel-space data vector of length n_pix.
C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).

Returns:

Quadratic form value d^T C^{-1} d.

Return type:

quadratic_form_from_prepared(data, K_chol)[source]

Compute d^T C^{-1} d using precomputed K Cholesky factor.

Return type:: float

release_pixel_projector()[source]

Drop V after SMW components are built.

QML hot paths (Fisher trace, Spectra weighted-data) read only _V_N_inv and _V_Ninv_VT. The remaining V consumers — get_covariance/get_inverse, to_basis, m-block compression, and PICSLike do_cross with the harmonic basis — must not be invoked after this call. Raises if m-block compression was requested at construction.

Return type:: None

setup()[source]

Build V and precompute SMW components.

Computes: - V: harmonic operator (n_modes × n_pix) - V N^{-1}: precomputed for SMW - V N^{-1} V^T: SMW kernel - V N V^T: harmonic-space noise covariance - log|N|: for determinant calculations - Derivative diagonals: E_ℓ for all ℓ

When switch optimization is enabled (lswitch_high < lmax): - V is built only for ℓ in [lswitch_low, lswitch_high] - S_fixed is computed for ℓ > lswitch_high using fiducial C_ell - N_eff = N + S_fixed is used instead of N This dramatically reduces the SMW dimension.

Return type:: None

class cosmocore.basis.PixelBasis(N, theta, phi, lmax_signal, beam=None, spins=None, compression_target='noise_weighted', C_ell=None, epsilon=None, mode_fraction=None, lmin_signal=None, lmin=None, lmax=None, S_fixed=None, fields=None)[source]

Pixel-space compression with projector (Gjerløw-like).

Stays in n_pix space with projector P_h, then uses eigenvalue decomposition to find optimal subspace. Compression: n_pix → dim.

The key operations are: - Projector: P_h = V^T V (harmonic projector, n_pix × n_pix but rank n_modes) - Eigendecomposition of compression matrix to find optimal modes - Data compression: d_c = P @ d where P = U^T (dim × n_pix)

This approach is more flexible than HarmonicBasis because it allows custom projectors to handle systematics (foreground deprojection, etc.).

Available compression bases (from Gjerløw et al. 2019):

harmonic: Pure harmonic projector P_h = V^T V. Selects modes based purely on harmonic content, ignoring noise properties.
noise_weighted: P_h N^{-1} P_h (default). Weights modes by inverse noise, preferring low-noise harmonic modes.
total_covariance: P_h C^{-1} P_h where C = N + S. Uses full signal+noise covariance for optimal compression.
snr: S^{1/2} N^{-1} S^{1/2}. Signal-to-noise ratio matrix, prioritizing modes with highest SNR.

Parameters:

N (numpy.ndarray) – Noise covariance matrix of shape (n_pix, n_pix). The basis takes ownership of this buffer; _factorise_noise() overwrites it in place during setup (V-based mode only; pixel-direct keeps it).
theta (numpy.ndarray) – Colatitude angles for active pixels in radians.
phi (numpy.ndarray) – Longitude angles for active pixels in radians.
lmax (int) – Maximum multipole for harmonic expansion.
beam (numpy.ndarray or None, optional) – Beam window function B_ℓ for ℓ=2 to lmax.
basis (str, default "noise_weighted") – Compression basis to use. Options: - “harmonic”: P_h = V^T V (pure harmonic projector) - “noise_weighted”: P_h N^{-1} P_h (inverse noise weighting) - “total_covariance”: P_h C^{-1} P_h where C = N + S - “snr”: S^{1/2} N^{-1} S^{1/2} (signal-to-noise ratio)
C_ell (numpy.ndarray or None, optional) – Power spectrum values for ell = 2 to lmax. Required for “total_covariance” and “snr” bases.
epsilon (float, optional) – Eigenvalue threshold relative to maximum.
mode_fraction (float, optional) – Fraction of modes to keep.

dim

Number of modes kept after compression (initially n_pix).

Type:: int

Examples

>>> import numpy as np
>>> from cosmocore.basis import PixelBasis
>>> N = np.diag(noise_variance)  # Noise covariance matrix
>>> # Compressed mode: pass epsilon at construction.
>>> ppc = PixelBasis(
...     N, theta, phi, lmax_signal=100,
...     epsilon=1e-4, compression_target="noise_weighted",
... )
>>> ppc.setup()
>>> fisher_element = ppc.compute_fisher_element(C_ell, ell_i=10, ell_j=10)

References

__init__(N, theta, phi, lmax_signal, beam=None, spins=None, compression_target='noise_weighted', C_ell=None, epsilon=None, mode_fraction=None, lmin_signal=None, lmin=None, lmax=None, S_fixed=None, fields=None)[source]

Initialize computation basis.

Parameters:

N (numpy.ndarray) – Noise covariance matrix (n_pix_total, n_pix_total). The basis takes ownership; _factorise_noise() overwrites this buffer in place during setup.
theta (numpy.ndarray or tuple of numpy.ndarray) – Colatitude angles for active pixels in radians. Single array for single-field, tuple of arrays for multi-field (one per component).
phi (numpy.ndarray or tuple of numpy.ndarray) – Longitude angles for active pixels in radians. Single array for single-field, tuple of arrays for multi-field (one per component).
lmax_signal (int) – Signal-cov ceiling. The basis represents multipoles up to and including this value.
beam (numpy.ndarray or None, optional) – Beam window function B_ℓ for ℓ=0..lmax_signal.
spins (list of int or None, optional) – Spin weight for each component (0 for scalar, 2 for polarization). Default is [0, …] for all components. For spin-2 components, theta/phi represent physical pixel locations but V is built for (Q, U) → (E, B) transformation with doubled dimensions.
lmin_signal (list of int or None, optional) – Per-component signal-cov floor. Defaults to [2]*n_components. Each entry must satisfy lmin_signal[i] >= |spins[i]|.
lmin (int or None, optional) – Inference window lower bound. Defaults to min(lmin_signal). Together with lmax enables the switch optimisation: S_fixed absorbs the signal contribution outside [lmin, lmax].
lmax (int or None, optional) – Inference window upper bound. Defaults to lmax_signal. Multipoles in (lmax, lmax_signal] are absorbed into S_fixed.
fiducial_C_ell (numpy.ndarray or None, optional) – Deprecated: Use S_fixed instead. Fiducial power spectrum for fixed multipoles (ℓ outside [lmin, lmax]).
S_fixed (numpy.ndarray or None, optional) – Precomputed signal matrix for fixed multipoles outside [lmin, lmax]. Recommended way to pass the fixed signal contribution. Shape (n_pix_total, n_pix_total).

classmethod available_bases()[source]

Get available compression bases and their descriptions.

Returns:: Dictionary mapping basis names to their descriptions.
Return type:: dict

Examples

>>> PixelBasis.available_bases()
{'harmonic': 'P_h = V^T V (pure harmonic projector)',
 'noise_weighted': 'P_h N^{-1} P_h (inverse noise weighting)',
 'total_covariance': 'P_h C^{-1} P_h (full covariance weighting, requires C_ell)',
 'snr': 'S^{1/2} N^{-1} S^{1/2} (signal-to-noise ratio, requires C_ell)'}

property compression_ratio: float

Ratio of kept modes to original pixels.

Returns:: Compression ratio (1.0 means no compression).
Return type:: float

property compression_target: str | None

The selector configured at construction (compression_target kwarg).

Names the matrix that _apply_compression eigendecomposes: "harmonic", "noise_weighted", "total_covariance", or "snr". Returns the configured value regardless of whether compression has actually been applied yet.

compute_eigenspectrum(basis='noise_weighted', C_ell=None)[source]

Compute eigenvalue spectrum for a given compression basis.

This method computes the eigendecomposition without applying compression, allowing inspection of the spectrum to choose an appropriate threshold.

Parameters:

basis (str, default "noise_weighted") – Compression basis: “harmonic”, “noise_weighted”, “total_covariance”, “snr”.
C_ell (numpy.ndarray, optional) – Power spectrum values for ell = 2 to lmax. Required for “total_covariance” and “snr” bases.

Return type:

tuple[ndarray, ndarray]

Returns:

eigenvalues (numpy.ndarray) – Eigenvalues sorted in descending order.
normalized_eigenvalues (numpy.ndarray) – Eigenvalues normalized by maximum value (for threshold selection).

compute_eigenspectrum_per_field(basis='noise_weighted', C_ell=None)[source]

Compute per-component eigenvalue spectra.

Returns one entry per component with eigenvalues and normalized eigenvalues. For spin-2 components the result additionally contains separate E- and B-mode eigenspectra.

Parameters:

basis (str, default "noise_weighted") – Compression basis.
C_ell (numpy.ndarray, dict, or None) – Power spectrum (required for “total_covariance” and “snr” bases).

Returns:

One dict per component with keys: component, spin, label, eigenvalues, normalized_eigenvalues. Spin-2 components additionally have E_eigenvalues, E_normalized, B_eigenvalues, B_normalized.

Return type:

list of dict

compute_fisher_matrix(C_ell, spectra_list=None, ell_min=2, ell_max=None, *, symmetry_mode=None)[source]

Compute Fisher matrix.

Parameters:

C_ell (numpy.ndarray or dict) – Power spectrum. Single-field: numpy.ndarray. Multi-field: dict keyed by SpectrumKey.
spectra_list (list[SpectrumKey] or None) – For multi-field, the list of SpectrumKey instances enumerating the spectra to include. None for single-field.
ell_min (int) – Minimum multipole.
ell_max (int or None) – Maximum multipole.
symmetry_mode (SymmetryMode or None) – Forwarded to the per-ℓ derivative builder. None keeps the SYMMETRIC default. Required as DIRECTIONAL when spectra_list contains a CG cross-pair entry — the builder otherwise rejects the mode=2 slot used by CG.

Returns:

Fisher matrix.

Return type:

property eigenvalues: ndarray | None

Eigenvalues of kept modes (sorted descending).

Returns:: Eigenvalues if compression has been applied, None otherwise.
Return type:: numpy.ndarray or None

get_binned_derivative_direct(bin_idx, bins, beam_smoothing, key)[source]

Build the binned derivative dC^b in a single Legendre/Wigner pass.

Mathematically equivalent to summing per-ℓ _get_derivative_direct results with weights beam²(ℓ) over ℓ ∈ [bins.lmins[bin_idx], bins.lmaxs[bin_idx]], but uses one compute_*_contribution call instead of N (where N = bin width). The contribution functions accumulate Σ_ℓ cl[ell] × kernel[ell] per pixel pair on ℓ-indexed arrays, so populating cl[ell] = beam²(ell) inside the bin (zero elsewhere) yields the exact binned derivative.

Parameters:: key (SpectrumKey) – Identifies the spectrum (component pair + spin kind). Decomposed into (comp_i, comp_j, mode) once at entry for the Numba kernels below.
Return type:: ndarray

get_covariance(C_ell)[source]

Compute the basis-space covariance.

In pixel-direct mode this is the exact N + S; on a truncated compressed basis it is U^T (N + S) U — the restriction of the full covariance to the kept subspace.

Parameters:: C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
Returns:: Basis-space covariance, shape (dim, dim).
Return type:: numpy.ndarray

get_derivative_matrix(ell, key, *, symmetry_mode=None)[source]

Get derivative matrix ∂C/∂C_ℓ in the basis space.

Parameters:

ell (int) – Multipole.
key (SpectrumKey) – Identifies the spectrum (component pair + spin kind).
symmetry_mode (SymmetryMode or None) – None keeps the SYMMETRIC default.

Returns:

Derivative matrix of shape (dim, dim).

Return type:

get_full_inverse(C_ell)[source]

Reconstruct full n_pix x n_pix inverse from the basis-space form.

Exact in pixel-direct mode (U is the identity). Lossy on a truncated compressed pixel basis — the result lives in the kept subspace and is zero on the truncated complement.

Return type:: ndarray

get_logdet(C_ell)[source]

Compute the log determinant of the basis-space covariance.

In pixel-direct mode this equals the exact full log|N + S|. On a truncated compressed basis it is the logdet of the restricted operator U^T (N + S) U — not the full log|N + S|; see get_full_logdet() for the ABC contract.

Parameters:: C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
Return type:: float

get_noise_for_bias()[source]

U^T N_raw U — raw noise projected into the eigenmode basis.

See ComputationBasis.get_noise_for_bias() for the cross-basis contract. Use this for noise-bias computations: Tr[C^{-1} N C^{-1} dS] requires the actual noise N, not the effective N_eff that includes the absorbed high-ℓ signal.

Without switch optimization, N_eff = N and this is identical to get_covariance(0).

Return type:: ndarray

get_projected_inverse(C_ell)[source]

Compute the basis-space inverse covariance.

In pixel-direct mode this is the exact (N + S)^{-1}; on a truncated compressed basis it is (U^T (N + S) U)^{-1} — the inverse of the restricted operator on the kept subspace.

Parameters:: C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
Returns:: Basis-space inverse covariance, shape (dim, dim).
Return type:: numpy.ndarray

get_weighted_data(data, C_ell, C_c_inv=None, stable_inner_inv=None)[source]

Compute C^{-1} @ d for QML estimation.

Parameters:

data (numpy.ndarray) – Pixel-space data vector of shape (n_pix,) or (n_pix, n_sims).
C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
C_c_inv (numpy.ndarray, optional) – Precomputed inverse covariance.
stable_inner_inv (numpy.ndarray, optional) – Unused for pixel basis; accepted for ABC signature parity.

Returns:

Weighted data of shape (dim,) or (dim, n_sims).

Return type:

property method: str: Compression method name.

plot_eigenvalue_comparison(bases=None, C_ell=None, axes=None, log_scale=True)[source]

Compare eigenvalue spectra across different compression bases.

Creates one subplot per component, overlaying the different bases.

Parameters:

bases (list of str or None) – Compression bases to compare. Default: all available (or just “harmonic”/”noise_weighted” if C_ell is not provided).
C_ell (numpy.ndarray, dict, or None) – Power spectrum (required for “total_covariance” and “snr” bases).
axes (numpy.ndarray of Axes or None) – Pre-created axes array (length n_components). If None, a new figure is created.
log_scale (bool, default True) – Whether to use logarithmic y-axis.

Return type:

tuple[Figure, ndarray]

Returns:

fig (matplotlib.figure.Figure) – The figure containing the plot.
axes (numpy.ndarray) – 1-D array of Axes (length n_components).

plot_eigenvalue_spectrum(basis='noise_weighted', C_ell=None, axes=None, log_scale=True, show_threshold_lines=True, threshold_values=None, show_eb_split=True)[source]

Plot eigenvalue spectrum for compression threshold selection.

Creates one subplot per component. The y-axis shows eigenvalues normalized by the maximum value, so values can be directly used as the epsilon constructor kwarg of PixelBasis. For spin-2 components the E and B sub-spectra are shown as dashed curves when show_eb_split is True.

Parameters:

basis (str, default "noise_weighted") – Compression basis.
C_ell (numpy.ndarray, dict, or None) – Power spectrum (required for “total_covariance” and “snr” bases).
axes (numpy.ndarray of Axes or None) – Pre-created axes array (length n_components). If None, a new figure is created.
log_scale (bool, default True) – Whether to use logarithmic y-axis.
show_threshold_lines (bool, default True) – Whether to show reference threshold lines.
threshold_values (list of float or None) – Custom threshold values to show. Default: [1e-2, 1e-4, 1e-6, 1e-8].
show_eb_split (bool, default True) – For spin-2 components, overlay E and B sub-spectra.

Return type:

tuple[Figure, ndarray]

Returns:

fig (matplotlib.figure.Figure) – The figure containing the plot.
axes (numpy.ndarray) – 1-D array of Axes (length n_components).

Examples

>>> # Pick threshold by plotting first, then construct a compressed basis.
>>> # The probe construction uses epsilon=0.0 to opt into the V-based
>>> # path (keeps all modes) so plot_eigenvalue_spectrum has V to read.
>>> ppc_probe = PixelBasis(
...     N, theta, phi, lmax_signal=100, epsilon=0.0,
... )
>>> ppc_probe.setup()
>>> fig, axes = ppc_probe.plot_eigenvalue_spectrum(basis="noise_weighted")
>>> # From the plot, decide threshold (e.g., 1e-4).
>>> ppc = PixelBasis(
...     N, theta, phi, lmax_signal=100,
...     epsilon=1e-4, compression_target="noise_weighted",
... )
>>> ppc.setup()

prepare_for_basis(C_ell_dict)[source]

Pre-bake the basis-space inverse and logdet for reuse across sims.

Return type:: BasisPrepared

property projector: ndarray

Get the projection matrix U^T (dim × n_pix).

Maps pixel space to the basis. In pixel-direct mode the projector is the identity; in compressed mode it is U^T, the transpose of the kept eigenvectors.

Raises:: RuntimeError – If compression has not been applied yet (V-based mode only).

quadratic_form(data, C_ell)[source]

Compute d^T C^{-1} d using the basis-space inverse.

In pixel-direct mode this equals the full pixel-space quadratic form d^T (N + S)^{-1} d. On a truncated compressed basis it is (U^T d)^T (U^T (N + S) U)^{-1} (U^T d) — the quadratic form of the restricted operator on the kept subspace.

Parameters:

data (numpy.ndarray) – Pixel-space data vector of length n_pix.
C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).

Returns:

Quadratic form value.

Return type:

quadratic_form_from_prepared(data, C_c_inv)[source]

Compute d^T C^{-1} d using a pre-baked basis-space inverse.

Return type:: float

setup()[source]

Initialize the basis.

In V-based mode: build V and P_h, optionally apply eigenvalue compression. In direct mode (n_pix small enough that V is overkill): skip V construction and rely on existing pixel-space machinery (compute_signal_matrix, do_derivative_step).

Return type:: None

to_basis(data)[source]

Project pixel-space data into the basis.

Pixel-direct mode short-circuits the multiplication against the identity projector (otherwise an n_pix × n_sims matmul against an n_pix × n_pix identity, allocated and run on every call).

Return type:: ndarray

class cosmocore.basis.BasisPrepared(factor: np.ndarray, logdet: float)[source]

Bases: NamedTuple

Per-C_ell quantities pre-computed for fast likelihood evaluation.

Returned by ComputationBasis.prepare_for_basis(). The stored factor is basis-specific in form but always usable as the second argument of ComputationBasis.quadratic_form_from_prepared().

Parameters:

factor (numpy.ndarray) – K Cholesky factor (harmonic) or basis-space inverse (pixel).
logdet (float) –
Log determinant. Semantics are basis-specific and intentionally match what quadratic_form_from_prepared consumes:
- HarmonicBasis: exact full pixel-space log|N + S| via SMW.
- PixelBasis in pixel-direct mode: exact full log|N + S| (U is the identity, so basis-space equals full-space).
- PixelBasis on a truncated compressed basis: basis-space log|U^T (N + S) U| — the logdet of the restricted operator, NOT the full log|N + S|. Combine with quadratic_form_from_prepared() (also basis-space) for an internally-consistent basis-space likelihood.

factor: ndarray: Alias for field number 0

logdet: float: Alias for field number 1

cosmocore.basis.create_computation_basis(method, N, theta, phi, lmax_signal, **kwargs)[source]

Factory function to create the appropriate computation basis implementation.

Parameters:

method (str) –
Computation basis method: - “harmonic”: Tegmark-like direct harmonic transformation - “pixel”: Gjerløw-like pixel-space projector with eigenvalue truncation - “auto”: Pick the cheaper path using a leading-order cost model.

Compares ~n_modes^3 (harmonic + sparse traces) vs ~(n_bins+1)*n_pix^3 (pixel-direct). Pass n_bins to refine the estimate; defaults to lmax-1 (unbinned, worst case for pixel-direct).
N (numpy.ndarray) – Noise covariance matrix. The basis takes ownership; the caller should drop its reference after construction.
theta (numpy.ndarray) – Colatitude angles for active pixels in radians.
phi (numpy.ndarray) – Longitude angles for active pixels in radians.
lmax (int) – Maximum multipole for harmonic expansion.
**kwargs – Additional keyword arguments passed to the basis constructor (beam, spins, basis, C_ell, epsilon, mode_fraction, fields, lmin_signal, lmin, lmax, etc.). Arguments not accepted by the chosen class are silently ignored.

Returns:

Configured basis instance (not yet set up — call .setup()).

Return type:

ComputationBasis

Submodules

cosmocore.basis.base

Base class for computation basis methods.

This module provides the abstract base class that defines the interface for all computation basis methods used in CMB Fisher matrix computation.

class cosmocore.basis.base.BasisPrepared(factor: np.ndarray, logdet: float)[source]

Bases: NamedTuple

Per-C_ell quantities pre-computed for fast likelihood evaluation.

Returned by ComputationBasis.prepare_for_basis(). The stored factor is basis-specific in form but always usable as the second argument of ComputationBasis.quadratic_form_from_prepared().

Parameters:

factor (numpy.ndarray) – K Cholesky factor (harmonic) or basis-space inverse (pixel).
logdet (float) –
Log determinant. Semantics are basis-specific and intentionally match what quadratic_form_from_prepared consumes:
- HarmonicBasis: exact full pixel-space log|N + S| via SMW.
- PixelBasis in pixel-direct mode: exact full log|N + S| (U is the identity, so basis-space equals full-space).
- PixelBasis on a truncated compressed basis: basis-space log|U^T (N + S) U| — the logdet of the restricted operator, NOT the full log|N + S|. Combine with quadratic_form_from_prepared() (also basis-space) for an internally-consistent basis-space likelihood.

factor: ndarray: Alias for field number 0

logdet: float: Alias for field number 1

class cosmocore.basis.base.ComputationBasis(N, theta, phi, lmax_signal, beam=None, spins=None, lmin_signal=None, lmin=None, lmax=None, fiducial_C_ell=None, S_fixed=None)[source]

Bases: ABC

Abstract base class for computation basis methods.

This class defines the interface for computation basis methods and provides shared implementations for spherical harmonic evaluation and ell-to-mode mapping.

Supports both single-field and multi-field configurations. For multi-field, theta and phi are passed as tuples of arrays (one per component), and the harmonic operator V is built as a block-diagonal matrix.

Parameters:

N (numpy.ndarray) – Noise covariance matrix of shape (n_pix_total, n_pix_total). The basis takes ownership of this buffer; the in-place Cholesky factor overwrites it during setup.
theta (numpy.ndarray or tuple of numpy.ndarray) – Colatitude angles for active pixels in radians. Single array for single-field, tuple of arrays for multi-field.
phi (numpy.ndarray or tuple of numpy.ndarray) – Longitude angles for active pixels in radians. Single array for single-field, tuple of arrays for multi-field.
lmax (int) – Maximum multipole for harmonic expansion.
beam (numpy.ndarray or None, optional) – Beam window function B_ℓ for ℓ=2 to lmax. Shape should be (lmax-1,). If provided, the harmonic operator V is multiplied by beam factors so that V_ℓm = B_ℓ × Y_ℓm. If None, no beam correction is applied.
spins (list of int or None, optional) – Spin weight for each component. Default is [0, 0, …] (all spin-0). Use spin=2 for polarization (Q/U) fields, which doubles pixel count and uses spin-weighted spherical harmonics for E/B decomposition.

n_pix

Total number of active pixels across all components.

Type:: int

n_modes

Number of harmonic modes per component (for ℓ=2 to lmax).

Type:: int

n_modes_total

Total harmonic modes across all components (n_components × n_modes).

Type:: int

dim

Number of modes kept after compression (if applicable).

Type:: int

n_components

Number of field components (1 for single-field, N for multi-field).

Type:: int

__init__(N, theta, phi, lmax_signal, beam=None, spins=None, lmin_signal=None, lmin=None, lmax=None, fiducial_C_ell=None, S_fixed=None)[source]

Initialize computation basis.

Parameters:

N (numpy.ndarray) – Noise covariance matrix (n_pix_total, n_pix_total). The basis takes ownership; _factorise_noise() overwrites this buffer in place during setup.
theta (numpy.ndarray or tuple of numpy.ndarray) – Colatitude angles for active pixels in radians. Single array for single-field, tuple of arrays for multi-field (one per component).
phi (numpy.ndarray or tuple of numpy.ndarray) – Longitude angles for active pixels in radians. Single array for single-field, tuple of arrays for multi-field (one per component).
lmax_signal (int) – Signal-cov ceiling. The basis represents multipoles up to and including this value.
beam (numpy.ndarray or None, optional) – Beam window function B_ℓ for ℓ=0..lmax_signal.
spins (list of int or None, optional) – Spin weight for each component (0 for scalar, 2 for polarization). Default is [0, …] for all components. For spin-2 components, theta/phi represent physical pixel locations but V is built for (Q, U) → (E, B) transformation with doubled dimensions.
lmin_signal (list of int or None, optional) – Per-component signal-cov floor. Defaults to [2]*n_components. Each entry must satisfy lmin_signal[i] >= |spins[i]|.
lmin (int or None, optional) – Inference window lower bound. Defaults to min(lmin_signal). Together with lmax enables the switch optimisation: S_fixed absorbs the signal contribution outside [lmin, lmax].
lmax (int or None, optional) – Inference window upper bound. Defaults to lmax_signal. Multipoles in (lmax, lmax_signal] are absorbed into S_fixed.
fiducial_C_ell (numpy.ndarray or None, optional) – Deprecated: Use S_fixed instead. Fiducial power spectrum for fixed multipoles (ℓ outside [lmin, lmax]).
S_fixed (numpy.ndarray or None, optional) – Precomputed signal matrix for fixed multipoles outside [lmin, lmax]. Recommended way to pass the fixed signal contribution. Shape (n_pix_total, n_pix_total).

abstractmethod setup()[source]

Initialize basis-specific components.

Must be called after initialization before using any other methods.

Return type:: None

property N_inv: ndarray | None

abstract property method: str

“harmonic” or “pixel”.

Type:: Computation basis name

abstract property projector: ndarray

Get the projection matrix that maps pixel space to the basis.

The fundamental change-of-basis operator: - HarmonicBasis: V (n_modes × n_pix) - PixelBasis: U^T (dim × n_pix) when compression is active;

identity-equivalent in pixel-direct mode.

Returns:: Projection matrix of shape (dim, n_pix).
Return type:: numpy.ndarray

abstractmethod get_projected_inverse(C_ell)[source]

Get the projected inverse for Fisher computation.

Parameters:: C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
Returns:: Projected inverse covariance matrix.
Return type:: numpy.ndarray

abstractmethod get_derivative_matrix(ell, key, *, symmetry_mode=None)[source]

Get the derivative matrix ∂C/∂C_ℓ in the basis space.

Parameters:

ell (int) – Multipole for which to compute the derivative.
key (SpectrumKey) – Identifies the spectrum (component pair + spin kind).
symmetry_mode (SymmetryMode or None) – None keeps the SYMMETRIC default.

Returns:

Derivative matrix of shape (dim, dim).

Return type:

abstractmethod get_covariance(C_ell)[source]

Compute the covariance matrix in the basis space.

HarmonicBasis: C̄ = V @ N @ V^T + Λ
PixelBasis: C_c = U^T @ C @ U

Parameters:: C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
Returns:: Basis-space covariance of shape (dim, dim).
Return type:: numpy.ndarray

abstractmethod get_weighted_data(data, C_ell, C_c_inv=None, stable_inner_inv=None)[source]

Compute weighted basis-space data for QML estimation.

Parameters:

data (numpy.ndarray) – Pixel-space data vector of shape (n_pix,) or (n_pix, n_sims).
C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
C_c_inv (numpy.ndarray, optional) – Precomputed basis-space inverse. Used by PixelBasis only; ignored by HarmonicBasis.
stable_inner_inv (numpy.ndarray, optional) – Precomputed (I + Lambda M)^{-1}, as returned by HarmonicBasis.prepare_stable_inner_inv(). Used by HarmonicBasis only; ignored by PixelBasis.

Returns:

Weighted basis-space data of shape (dim,) or (dim, n_sims).

Return type:

abstractmethod get_noise_for_bias()[source]

Basis-projected noise consumed by the QML noise-bias sandwich.

Spectra computes Cov(w | noise) = X get_noise_for_bias() X^T where X is the basis’s natural inverse-equivalent (A = (I + Lambda M)^{-T} for the harmonic basis; C_c^{-1} for the pixel basis). The form of the returned matrix is therefore basis-specific:

HarmonicBasis: V N_eff^{-1} N N_eff^{-1} V^T (the SMW intermediate T from N_eff = N + S_fixed when the switch optimisation is active; otherwise V N^{-1} N N^{-1} V^T).
PixelBasis: U^T N U (raw noise projected once into the eigenmode basis).

The two are not interchangeable. Callers must compose with the basis’s own inverse-equivalent.

Returns:: Symmetric (dim, dim) matrix.
Return type:: numpy.ndarray

get_inverse(C_ell)[source]

Compute the inverse of the basis-space covariance.

Parameters:: C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
Returns:: Basis-space inverse covariance of shape (dim, dim).
Return type:: numpy.ndarray

get_logdet(C_ell)[source]

Compute the log determinant of the basis-space covariance.

Parameters:: C_ell (numpy.ndarray) – Power spectrum values for ell = 2 to lmax.
Returns:: Log determinant of basis-space covariance.
Return type:: float

get_full_logdet(C_ell)[source]

Return log|N + S|, the full pixel-space log determinant.

HarmonicBasis overrides this to compute it exactly via the SMW formula. The default below applies to PixelBasis and is exact only in pixel-direct mode (U is the identity, so the basis-space logdet equals the full logdet).

On a truncated compressed pixel basis the full logdet cannot be recovered from the kept quantities: the discarded complement of N + S is needed and is not stored. The default falls through to get_logdet(), which returns the logdet of the restricted operator U^T (N + S) U — a different matrix, not an approximation to log|N + S|. Callers that need the exact full logdet on a truncated basis must arrange for it elsewhere.

Parameters:: C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
Return type:: float

abstractmethod get_full_inverse(C_ell)[source]

Return the full n_pix x n_pix inverse covariance.

For HarmonicBasis: exact (N + S)^{-1} via the SMW formula. For PixelBasis in pixel-direct mode: also exact (U is the identity, so U @ get_inverse(C_ell) @ U.T equals the full inverse).

For PixelBasis on a truncated compressed basis: the returned matrix is U @ get_inverse(C_ell) @ U.T, which lives in the kept subspace and is zero on the truncated complement. It is not the inverse of the full pixel-space covariance; it is the inverse of the restricted operator, lifted back to n_pix dimensions. Use with care.

Callers who only need the basis-space inverse should use get_inverse() instead — cheaper and equivalent in basis-space algebra.

Return type:: ndarray

abstractmethod quadratic_form(data, C_ell)[source]

Compute d^T C^{-1} d for full pixel-space data.

The basis chooses its fastest path. Returns a single float for a 1-D data vector, or a length-n_sims array if data is shape (n_pix, n_sims).

Parameters:

data (numpy.ndarray) – Pixel-space data, shape (n_pix,) or (n_pix, n_sims).
C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).

Return type:

abstractmethod quadratic_form_from_prepared(data, factor)[source]

Compute d^T C^{-1} d using a pre-baked factor from prepare_for_basis().

Fast path for evaluating the quadratic form many times against the same C_ell: caller calls prepare_for_basis() once and reuses the resulting factor across simulations.

Parameters:

data (numpy.ndarray) – Pixel-space data, shape (n_pix,) or (n_pix, n_sims).
factor (numpy.ndarray) – Basis-specific pre-baked quantity: BasisPrepared.factor from prepare_for_basis().

Return type:

abstractmethod prepare_for_basis(C_ell_dict)[source]

Pre-bake per-C_ell quantities for fast likelihood reuse.

Returns a BasisPrepared carrying the basis-specific factor and the full log determinant. The factor is consumable by quadratic_form_from_prepared() on the same basis instance.

Parameters:: C_ell_dict (dict) – Multi-field power spectrum dict keyed by SpectrumKey.
Return type:: BasisPrepared

release_pixel_projector()[source]

No-op default. HarmonicBasis overrides to drop V post-SMW build.

Return type:: None

to_basis(data)[source]

Project pixel data into the basis: d_basis = P @ d.

Uses the basis-specific projector P: - HarmonicBasis: P = V (n_modes × n_pix) - PixelBasis: P = U^T (dim × n_pix) when compression is active;

identity short-circuit in pixel-direct mode.

Parameters:: data (numpy.ndarray) – Pixel-space data vector of shape (n_pix,) or (n_pix, n_sims).
Returns:: Basis-space data of shape (dim,) or (dim, n_sims).
Return type:: numpy.ndarray

get_cls_vector()[source]

Get a placeholder C_ell vector for the configured lmax_signal.

Returns:: ℓ-indexed array of zeros with length lmax_signal + 1.
Return type:: numpy.ndarray

property compression_ratio: float

Ratio of kept modes to original modes.

Returns:: Compression ratio (1.0 means no compression).
Return type:: float

cosmocore.basis.harmonic

Harmonic computation basis (Tegmark-like) for CMB Fisher matrix computation.

This module implements direct harmonic space projection where data is transformed directly to n_modes dimensions via the harmonic operator V.

class cosmocore.basis.harmonic.HarmonicBasis(N, theta, phi, lmax_signal, beam=None, spins=None, lmin_signal=None, lmin=None, lmax=None, fiducial_C_ell=None, S_fixed=None, compress=False, delta_m=0)[source]

Direct harmonic space computation basis (Tegmark-like).

Transforms directly to n_modes dimensions via V. Fast and efficient when n_modes << n_pix. No additional eigenvalue truncation.

The key operations are: - Data projection: d̄ = P @ d where P = V - Projected covariance: C̄ = V N V^T + Λ - Projected inverse: V C^{-1} V^T = M - M K^{-1} M (via SMW)

Parameters:

N (numpy.ndarray) – Noise covariance matrix of shape (n_pix, n_pix). The basis takes ownership of this buffer; _factorise_noise() overwrites it in place during setup.
theta (numpy.ndarray) – Colatitude angles for active pixels in radians.
phi (numpy.ndarray) – Longitude angles for active pixels in radians.
lmax (int) – Maximum multipole for harmonic expansion.
beam (numpy.ndarray or None, optional) – Beam window function B_ℓ for ℓ=2 to lmax.
compress (bool, optional) – If True, use m-block compression for K inversion. Approximates K as block-diagonal in azimuthal quantum number |m|, giving ~lmax^2 speedup. Default is False.
delta_m (int, optional) – Bandwidth for m-block coupling. delta_m=0 means block-diagonal (no coupling between different |m|). delta_m=lmax recovers exact result. Default is 0.

Examples

>>> import numpy as np
>>> from cosmocore.basis import HarmonicBasis
>>> hc = HarmonicBasis(N, theta, phi, lmax=100)
>>> hc.setup()
>>> fisher_element = hc.compute_fisher_element(C_ell, ell_i=10, ell_j=10)

References

__init__(N, theta, phi, lmax_signal, beam=None, spins=None, lmin_signal=None, lmin=None, lmax=None, fiducial_C_ell=None, S_fixed=None, compress=False, delta_m=0)[source]

Initialize computation basis.

Parameters:

N (numpy.ndarray) – Noise covariance matrix (n_pix_total, n_pix_total). The basis takes ownership; _factorise_noise() overwrites this buffer in place during setup.
theta (numpy.ndarray or tuple of numpy.ndarray) – Colatitude angles for active pixels in radians. Single array for single-field, tuple of arrays for multi-field (one per component).
phi (numpy.ndarray or tuple of numpy.ndarray) – Longitude angles for active pixels in radians. Single array for single-field, tuple of arrays for multi-field (one per component).
lmax_signal (int) – Signal-cov ceiling. The basis represents multipoles up to and including this value.
beam (numpy.ndarray or None, optional) – Beam window function B_ℓ for ℓ=0..lmax_signal.
spins (list of int or None, optional) – Spin weight for each component (0 for scalar, 2 for polarization). Default is [0, …] for all components. For spin-2 components, theta/phi represent physical pixel locations but V is built for (Q, U) → (E, B) transformation with doubled dimensions.
lmin_signal (list of int or None, optional) – Per-component signal-cov floor. Defaults to [2]*n_components. Each entry must satisfy lmin_signal[i] >= |spins[i]|.
lmin (int or None, optional) – Inference window lower bound. Defaults to min(lmin_signal). Together with lmax enables the switch optimisation: S_fixed absorbs the signal contribution outside [lmin, lmax].
lmax (int or None, optional) – Inference window upper bound. Defaults to lmax_signal. Multipoles in (lmax, lmax_signal] are absorbed into S_fixed.
fiducial_C_ell (numpy.ndarray or None, optional) – Deprecated: Use S_fixed instead. Fiducial power spectrum for fixed multipoles (ℓ outside [lmin, lmax]).
S_fixed (numpy.ndarray or None, optional) – Precomputed signal matrix for fixed multipoles outside [lmin, lmax]. Recommended way to pass the fixed signal contribution. Shape (n_pix_total, n_pix_total).

property method: str: Computation basis name.

property projector: ndarray: Projection matrix V (n_modes × n_pix).

release_pixel_projector()[source]

Drop V after SMW components are built.

QML hot paths (Fisher trace, Spectra weighted-data) read only _V_N_inv and _V_Ninv_VT. The remaining V consumers — get_covariance/get_inverse, to_basis, m-block compression, and PICSLike do_cross with the harmonic basis — must not be invoked after this call. Raises if m-block compression was requested at construction.

Return type:: None

setup()[source]

Build V and precompute SMW components.

Computes: - V: harmonic operator (n_modes × n_pix) - V N^{-1}: precomputed for SMW - V N^{-1} V^T: SMW kernel - V N V^T: harmonic-space noise covariance - log|N|: for determinant calculations - Derivative diagonals: E_ℓ for all ℓ

When switch optimization is enabled (lswitch_high < lmax): - V is built only for ℓ in [lswitch_low, lswitch_high] - S_fixed is computed for ℓ > lswitch_high using fiducial C_ell - N_eff = N + S_fixed is used instead of N This dramatically reduces the SMW dimension.

Return type:: None

prepare_stable_inner_inv(C_ell, lambda_matrix=None)[source]

Build (I + Lambda M)^{-1} for the stable QML algebra.

This matrix appears identically in three places of the QML algebra and replaces the unstable subtractive forms:

Projected inverse: V C^{-1} V^T = M @ (I + Lambda M)^{-1}
Data weighting: V C^{-1} d = (I + Lambda M)^{-1} @ V N^{-1} d
Noise bias matrix: A = I - M K^{-1} = (I + Lambda M)^{-1}

All three follow from the identity M - M K^{-1} M = M K^{-1} Lambda^{-1} = M (I + Lambda M)^{-1} and never subtract large nearly-equal matrices, so they remain accurate in the cosmic-variance-limited regime where the legacy SMW form loses precision (e.g. T at low ell with sub-uK polarization noise).

Parameters:

C_ell (np.ndarray or dict) – Power spectrum (single-field array or multi-field dict).
lambda_matrix (np.ndarray, optional) – Pre-built full Lambda matrix. If None, built from C_ell.

Returns:

(I + Lambda M)^{-1}, n_modes_total x n_modes_total.

Return type:

np.ndarray

get_projected_inverse(C_ell, field_groups=None)[source]

Compute V C^{-1} V^T efficiently using SMW formula.

Parameters:

C_ell (numpy.ndarray or dict) – Power spectrum. Can be array (single-field) or dict (multi-field).
field_groups (list of list of int or None, optional) – Independent field groups from _detect_field_blocks. If provided and contains more than one group, exploits field block-diagonal structure for faster K inversion.

Returns:

Projected inverse covariance V C^{-1} V^T.

Return type:

get_derivative_matrix(ell, key, *, symmetry_mode=None)[source]

Get the derivative matrix ∂S/∂C_ℓ in harmonic space.

Parameters:

ell (int) – Multipole.
key (SpectrumKey) – Identifies the spectrum (component pair + spin kind).
symmetry_mode (SymmetryMode or None) – None keeps the SYMMETRIC default.

Returns:

Derivative matrix of shape (n_modes, n_modes).

Return type:

get_covariance(C_ell)[source]

Compute harmonic-space covariance C̄ = V N V^T + Λ.

Parameters:: C_ell (numpy.ndarray or dict) – Power spectrum. Can be: - numpy.ndarray: C_ell values for ell = 2 to lmax (single-field) - dict: Multi-field with 2-tuple or 3-tuple keys
Returns:: Compressed covariance matrix.
Return type:: numpy.ndarray

get_full_inverse(C_ell)[source]

Compute full (N + S)^{-1} using the SMW formula.

Parameters:: C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
Returns:: Full inverse covariance matrix of shape (n_pix, n_pix).
Return type:: numpy.ndarray

get_full_logdet(C_ell)[source]

Compute log|N + S| using SMW formula.

Parameters:: C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
Returns:: Log determinant of the full covariance matrix.
Return type:: float

get_noise_for_bias()[source]

Return the SMW noise intermediate T = V N_eff^{-1} N N_eff^{-1} V^T.

See ComputationBasis.get_noise_for_bias() for the cross-basis contract. T is precomputed during basis setup.

Return type:: ndarray

get_weighted_data(data, C_ell, C_c_inv=None, stable_inner_inv=None)[source]

Compute w = V C^{-1} d via the stable SMW algebra.

w = (I + Lambda M)^{-T} (V N^{-1} d) — algebraically equivalent to M K^{-1} d but never subtracts large nearly-equal matrices, so it retains precision in the cosmic-variance-limited regime where the legacy subtractive form y - M K^{-1} y cancelled.

Parameters:

data (numpy.ndarray) – Pixel-space data, shape (n_pix,) or (n_pix, n_sims).
C_ell (numpy.ndarray or dict) – Power spectrum (single-field array or multi-field dict).
C_c_inv (numpy.ndarray, optional) – Unused for harmonic basis; accepted for ABC signature parity.
stable_inner_inv (numpy.ndarray, optional) – Precomputed (I + Lambda M)^{-1}. If supplied, avoids rebuilding it internally. Callers that weight several data vectors with the same C_ell should build it once via prepare_stable_inner_inv() and pass it in.

Returns:

Weighted harmonic-space data w = V C^{-1} d.

Return type:

quadratic_form(data, C_ell)[source]

Compute d^T C^{-1} d efficiently using SMW formula.

Parameters:

data (numpy.ndarray) – Pixel-space data vector of length n_pix.
C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).

Returns:

Quadratic form value d^T C^{-1} d.

Return type:

prepare_for_basis(C_ell_dict)[source]

Precompute K Cholesky factor and log determinant for reuse across sims.

Return type:: BasisPrepared

quadratic_form_from_prepared(data, K_chol)[source]

Compute d^T C^{-1} d using precomputed K Cholesky factor.

Return type:: float

get_derivative_diagonal(ell, key)[source]

Get derivative diagonal for an auto-pair diagonal SpectrumKey.

Returns a 1D vector — the diagonal of ∂S/∂C_ℓ — for use in the sparse-trace fast path. Valid only for key.comp_i == key.comp_j and diagonal kinds (SS, GG, CC) — i.e. TT, EE, or BB. Cross-component keys and GC/CG/SG/SC/GS/CS raise ValueError.

Return type:: ndarray

compute_fisher_matrix(C_ell, spectra_list=None, ell_min=2, ell_max=None, *, symmetry_mode=None)[source]

Compute Fisher matrix.

Parameters:

C_ell (numpy.ndarray or dict) – Power spectrum. Single-field: numpy.ndarray (spectra_list=None). Multi-field: dict keyed by SpectrumKey (spectra_list required).
spectra_list (list[SpectrumKey] or None) – For multi-field, the list of SpectrumKey instances enumerating the spectra to include. None for single-field.
ell_min (int) – Minimum multipole.
ell_max (int or None) – Maximum multipole.
symmetry_mode (SymmetryMode or None) – Forwarded to the per-ℓ derivative builder. None keeps the SYMMETRIC default. Required as DIRECTIONAL when spectra_list contains a CG cross-pair entry — the builder otherwise rejects the mode=2 slot used by CG.

Returns:

Fisher matrix.

Return type:

cosmocore.basis.harmonic_basis

Harmonic basis construction for computation basis methods.

This module contains the HarmonicBasisBuilder class which encapsulates all spherical harmonic operator (V), Lambda matrix, and derivative matrix construction logic. It is an internal helper used by ComputationBasis and its subclasses.

class cosmocore.basis.harmonic_basis.HarmonicBasisBuilder(parent)[source]

Bases: object

Builds and caches harmonic operator V, Lambda matrices, and derivative matrices.

This is an internal helper class owned by ComputationBasis. It groups all the spherical harmonic basis construction code that was previously spread across ComputationBasis methods.

Parameters:: parent (ComputationBasis) – Parent computation basis instance providing configuration. The following attributes are read (all set during ComputationBasis.__init__): _theta_tuple, _phi_tuple, _spins, n_components, lmax, _lmin_smw, _lmax_smw, _n_modes_base, _n_modes_per_component, _n_modes_per_component_list, n_modes, n_modes_total, _mode_offsets, _pix_offsets, n_pix, _beam.

__init__(parent)[source]

build()[source]

Build harmonic operator V, mode mappings, and derivative diagonals.

Return type:: None

cosmocore.basis.pixel

Pixel-projected compression (Gjerløw-like) for CMB Fisher matrix computation.

This module implements pixel-space compression with a harmonic projector P_h, followed by eigenvalue decomposition to find the optimal subspace for compression.

Available compression bases (from Gjerløw et al. 2019): - “harmonic”: Pure harmonic projector P_h = V^T V - “noise_weighted”: P_h N^{-1} P_h (default) - inverse noise weighting - “total_covariance”: P_h C^{-1} P_h where C = N + S - full covariance weighting - “snr”: S^{1/2} N^{-1} S^{1/2} - signal-to-noise ratio matrix

class cosmocore.basis.pixel.PixelBasis(N, theta, phi, lmax_signal, beam=None, spins=None, compression_target='noise_weighted', C_ell=None, epsilon=None, mode_fraction=None, lmin_signal=None, lmin=None, lmax=None, S_fixed=None, fields=None)[source]

Pixel-space compression with projector (Gjerløw-like).

Stays in n_pix space with projector P_h, then uses eigenvalue decomposition to find optimal subspace. Compression: n_pix → dim.

The key operations are: - Projector: P_h = V^T V (harmonic projector, n_pix × n_pix but rank n_modes) - Eigendecomposition of compression matrix to find optimal modes - Data compression: d_c = P @ d where P = U^T (dim × n_pix)

This approach is more flexible than HarmonicBasis because it allows custom projectors to handle systematics (foreground deprojection, etc.).

Available compression bases (from Gjerløw et al. 2019):

harmonic: Pure harmonic projector P_h = V^T V. Selects modes based purely on harmonic content, ignoring noise properties.
noise_weighted: P_h N^{-1} P_h (default). Weights modes by inverse noise, preferring low-noise harmonic modes.
total_covariance: P_h C^{-1} P_h where C = N + S. Uses full signal+noise covariance for optimal compression.
snr: S^{1/2} N^{-1} S^{1/2}. Signal-to-noise ratio matrix, prioritizing modes with highest SNR.

Parameters:

N (numpy.ndarray) – Noise covariance matrix of shape (n_pix, n_pix). The basis takes ownership of this buffer; _factorise_noise() overwrites it in place during setup (V-based mode only; pixel-direct keeps it).
theta (numpy.ndarray) – Colatitude angles for active pixels in radians.
phi (numpy.ndarray) – Longitude angles for active pixels in radians.
lmax (int) – Maximum multipole for harmonic expansion.
beam (numpy.ndarray or None, optional) – Beam window function B_ℓ for ℓ=2 to lmax.
basis (str, default "noise_weighted") – Compression basis to use. Options: - “harmonic”: P_h = V^T V (pure harmonic projector) - “noise_weighted”: P_h N^{-1} P_h (inverse noise weighting) - “total_covariance”: P_h C^{-1} P_h where C = N + S - “snr”: S^{1/2} N^{-1} S^{1/2} (signal-to-noise ratio)
C_ell (numpy.ndarray or None, optional) – Power spectrum values for ell = 2 to lmax. Required for “total_covariance” and “snr” bases.
epsilon (float, optional) – Eigenvalue threshold relative to maximum.
mode_fraction (float, optional) – Fraction of modes to keep.

dim

Number of modes kept after compression (initially n_pix).

Type:: int

Examples

>>> import numpy as np
>>> from cosmocore.basis import PixelBasis
>>> N = np.diag(noise_variance)  # Noise covariance matrix
>>> # Compressed mode: pass epsilon at construction.
>>> ppc = PixelBasis(
...     N, theta, phi, lmax_signal=100,
...     epsilon=1e-4, compression_target="noise_weighted",
... )
>>> ppc.setup()
>>> fisher_element = ppc.compute_fisher_element(C_ell, ell_i=10, ell_j=10)

References

__init__(N, theta, phi, lmax_signal, beam=None, spins=None, compression_target='noise_weighted', C_ell=None, epsilon=None, mode_fraction=None, lmin_signal=None, lmin=None, lmax=None, S_fixed=None, fields=None)[source]

Initialize computation basis.

Parameters:

N (numpy.ndarray) – Noise covariance matrix (n_pix_total, n_pix_total). The basis takes ownership; _factorise_noise() overwrites this buffer in place during setup.
theta (numpy.ndarray or tuple of numpy.ndarray) – Colatitude angles for active pixels in radians. Single array for single-field, tuple of arrays for multi-field (one per component).
phi (numpy.ndarray or tuple of numpy.ndarray) – Longitude angles for active pixels in radians. Single array for single-field, tuple of arrays for multi-field (one per component).
lmax_signal (int) – Signal-cov ceiling. The basis represents multipoles up to and including this value.
beam (numpy.ndarray or None, optional) – Beam window function B_ℓ for ℓ=0..lmax_signal.
spins (list of int or None, optional) – Spin weight for each component (0 for scalar, 2 for polarization). Default is [0, …] for all components. For spin-2 components, theta/phi represent physical pixel locations but V is built for (Q, U) → (E, B) transformation with doubled dimensions.
lmin_signal (list of int or None, optional) – Per-component signal-cov floor. Defaults to [2]*n_components. Each entry must satisfy lmin_signal[i] >= |spins[i]|.
lmin (int or None, optional) – Inference window lower bound. Defaults to min(lmin_signal). Together with lmax enables the switch optimisation: S_fixed absorbs the signal contribution outside [lmin, lmax].
lmax (int or None, optional) – Inference window upper bound. Defaults to lmax_signal. Multipoles in (lmax, lmax_signal] are absorbed into S_fixed.
fiducial_C_ell (numpy.ndarray or None, optional) – Deprecated: Use S_fixed instead. Fiducial power spectrum for fixed multipoles (ℓ outside [lmin, lmax]).
S_fixed (numpy.ndarray or None, optional) – Precomputed signal matrix for fixed multipoles outside [lmin, lmax]. Recommended way to pass the fixed signal contribution. Shape (n_pix_total, n_pix_total).

property method: str: Compression method name.

property projector: ndarray

Get the projection matrix U^T (dim × n_pix).

Maps pixel space to the basis. In pixel-direct mode the projector is the identity; in compressed mode it is U^T, the transpose of the kept eigenvectors.

Raises:: RuntimeError – If compression has not been applied yet (V-based mode only).

to_basis(data)[source]

Project pixel-space data into the basis.

Pixel-direct mode short-circuits the multiplication against the identity projector (otherwise an n_pix × n_sims matmul against an n_pix × n_pix identity, allocated and run on every call).

Return type:: ndarray

setup()[source]

Initialize the basis.

In V-based mode: build V and P_h, optionally apply eigenvalue compression. In direct mode (n_pix small enough that V is overkill): skip V construction and rely on existing pixel-space machinery (compute_signal_matrix, do_derivative_step).

Return type:: None

get_noise_for_bias()[source]

U^T N_raw U — raw noise projected into the eigenmode basis.

See ComputationBasis.get_noise_for_bias() for the cross-basis contract. Use this for noise-bias computations: Tr[C^{-1} N C^{-1} dS] requires the actual noise N, not the effective N_eff that includes the absorbed high-ℓ signal.

Without switch optimization, N_eff = N and this is identical to get_covariance(0).

Return type:: ndarray

get_binned_derivative_direct(bin_idx, bins, beam_smoothing, key)[source]

Build the binned derivative dC^b in a single Legendre/Wigner pass.

Mathematically equivalent to summing per-ℓ _get_derivative_direct results with weights beam²(ℓ) over ℓ ∈ [bins.lmins[bin_idx], bins.lmaxs[bin_idx]], but uses one compute_*_contribution call instead of N (where N = bin width). The contribution functions accumulate Σ_ℓ cl[ell] × kernel[ell] per pixel pair on ℓ-indexed arrays, so populating cl[ell] = beam²(ell) inside the bin (zero elsewhere) yields the exact binned derivative.

Parameters:: key (SpectrumKey) – Identifies the spectrum (component pair + spin kind). Decomposed into (comp_i, comp_j, mode) once at entry for the Numba kernels below.
Return type:: ndarray

compute_eigenspectrum(basis='noise_weighted', C_ell=None)[source]

Compute eigenvalue spectrum for a given compression basis.

This method computes the eigendecomposition without applying compression, allowing inspection of the spectrum to choose an appropriate threshold.

Parameters:

basis (str, default "noise_weighted") – Compression basis: “harmonic”, “noise_weighted”, “total_covariance”, “snr”.
C_ell (numpy.ndarray, optional) – Power spectrum values for ell = 2 to lmax. Required for “total_covariance” and “snr” bases.

Return type:

tuple[ndarray, ndarray]

Returns:

eigenvalues (numpy.ndarray) – Eigenvalues sorted in descending order.
normalized_eigenvalues (numpy.ndarray) – Eigenvalues normalized by maximum value (for threshold selection).

compute_eigenspectrum_per_field(basis='noise_weighted', C_ell=None)[source]

Compute per-component eigenvalue spectra.

Returns one entry per component with eigenvalues and normalized eigenvalues. For spin-2 components the result additionally contains separate E- and B-mode eigenspectra.

Parameters:

basis (str, default "noise_weighted") – Compression basis.
C_ell (numpy.ndarray, dict, or None) – Power spectrum (required for “total_covariance” and “snr” bases).

Returns:

One dict per component with keys: component, spin, label, eigenvalues, normalized_eigenvalues. Spin-2 components additionally have E_eigenvalues, E_normalized, B_eigenvalues, B_normalized.

Return type:

list of dict

plot_eigenvalue_spectrum(basis='noise_weighted', C_ell=None, axes=None, log_scale=True, show_threshold_lines=True, threshold_values=None, show_eb_split=True)[source]

Plot eigenvalue spectrum for compression threshold selection.

Creates one subplot per component. The y-axis shows eigenvalues normalized by the maximum value, so values can be directly used as the epsilon constructor kwarg of PixelBasis. For spin-2 components the E and B sub-spectra are shown as dashed curves when show_eb_split is True.

Parameters:

basis (str, default "noise_weighted") – Compression basis.
C_ell (numpy.ndarray, dict, or None) – Power spectrum (required for “total_covariance” and “snr” bases).
axes (numpy.ndarray of Axes or None) – Pre-created axes array (length n_components). If None, a new figure is created.
log_scale (bool, default True) – Whether to use logarithmic y-axis.
show_threshold_lines (bool, default True) – Whether to show reference threshold lines.
threshold_values (list of float or None) – Custom threshold values to show. Default: [1e-2, 1e-4, 1e-6, 1e-8].
show_eb_split (bool, default True) – For spin-2 components, overlay E and B sub-spectra.

Return type:

tuple[Figure, ndarray]

Returns:

fig (matplotlib.figure.Figure) – The figure containing the plot.
axes (numpy.ndarray) – 1-D array of Axes (length n_components).

Examples

>>> # Pick threshold by plotting first, then construct a compressed basis.
>>> # The probe construction uses epsilon=0.0 to opt into the V-based
>>> # path (keeps all modes) so plot_eigenvalue_spectrum has V to read.
>>> ppc_probe = PixelBasis(
...     N, theta, phi, lmax_signal=100, epsilon=0.0,
... )
>>> ppc_probe.setup()
>>> fig, axes = ppc_probe.plot_eigenvalue_spectrum(basis="noise_weighted")
>>> # From the plot, decide threshold (e.g., 1e-4).
>>> ppc = PixelBasis(
...     N, theta, phi, lmax_signal=100,
...     epsilon=1e-4, compression_target="noise_weighted",
... )
>>> ppc.setup()

plot_eigenvalue_comparison(bases=None, C_ell=None, axes=None, log_scale=True)[source]

Compare eigenvalue spectra across different compression bases.

Creates one subplot per component, overlaying the different bases.

Parameters:

bases (list of str or None) – Compression bases to compare. Default: all available (or just “harmonic”/”noise_weighted” if C_ell is not provided).
C_ell (numpy.ndarray, dict, or None) – Power spectrum (required for “total_covariance” and “snr” bases).
axes (numpy.ndarray of Axes or None) – Pre-created axes array (length n_components). If None, a new figure is created.
log_scale (bool, default True) – Whether to use logarithmic y-axis.

Return type:

tuple[Figure, ndarray]

Returns:

fig (matplotlib.figure.Figure) – The figure containing the plot.
axes (numpy.ndarray) – 1-D array of Axes (length n_components).

get_projected_inverse(C_ell)[source]

Compute the basis-space inverse covariance.

In pixel-direct mode this is the exact (N + S)^{-1}; on a truncated compressed basis it is (U^T (N + S) U)^{-1} — the inverse of the restricted operator on the kept subspace.

Parameters:: C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
Returns:: Basis-space inverse covariance, shape (dim, dim).
Return type:: numpy.ndarray

get_derivative_matrix(ell, key, *, symmetry_mode=None)[source]

Get derivative matrix ∂C/∂C_ℓ in the basis space.

Parameters:

ell (int) – Multipole.
key (SpectrumKey) – Identifies the spectrum (component pair + spin kind).
symmetry_mode (SymmetryMode or None) – None keeps the SYMMETRIC default.

Returns:

Derivative matrix of shape (dim, dim).

Return type:

get_covariance(C_ell)[source]

Compute the basis-space covariance.

In pixel-direct mode this is the exact N + S; on a truncated compressed basis it is U^T (N + S) U — the restriction of the full covariance to the kept subspace.

Parameters:: C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
Returns:: Basis-space covariance, shape (dim, dim).
Return type:: numpy.ndarray

get_weighted_data(data, C_ell, C_c_inv=None, stable_inner_inv=None)[source]

Compute C^{-1} @ d for QML estimation.

Parameters:

data (numpy.ndarray) – Pixel-space data vector of shape (n_pix,) or (n_pix, n_sims).
C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).
C_c_inv (numpy.ndarray, optional) – Precomputed inverse covariance.
stable_inner_inv (numpy.ndarray, optional) – Unused for pixel basis; accepted for ABC signature parity.

Returns:

Weighted data of shape (dim,) or (dim, n_sims).

Return type:

quadratic_form(data, C_ell)[source]

Compute d^T C^{-1} d using the basis-space inverse.

In pixel-direct mode this equals the full pixel-space quadratic form d^T (N + S)^{-1} d. On a truncated compressed basis it is (U^T d)^T (U^T (N + S) U)^{-1} (U^T d) — the quadratic form of the restricted operator on the kept subspace.

Parameters:

data (numpy.ndarray) – Pixel-space data vector of length n_pix.
C_ell (numpy.ndarray or dict) – Power spectrum (array for single-field, dict for multi-field).

Returns:

Quadratic form value.

Return type:

compute_fisher_matrix(C_ell, spectra_list=None, ell_min=2, ell_max=None, *, symmetry_mode=None)[source]

Compute Fisher matrix.

Parameters:

C_ell (numpy.ndarray or dict) – Power spectrum. Single-field: numpy.ndarray. Multi-field: dict keyed by SpectrumKey.
spectra_list (list[SpectrumKey] or None) – For multi-field, the list of SpectrumKey instances enumerating the spectra to include. None for single-field.
ell_min (int) – Minimum multipole.
ell_max (int or None) – Maximum multipole.
symmetry_mode (SymmetryMode or None) – Forwarded to the per-ℓ derivative builder. None keeps the SYMMETRIC default. Required as DIRECTIONAL when spectra_list contains a CG cross-pair entry — the builder otherwise rejects the mode=2 slot used by CG.

Returns:

Fisher matrix.

Return type: