cosmocore.basics module

Basic mathematical utilities for cosmological computations.

This package provides fundamental mathematical functions optimized for cosmological analysis, including Legendre polynomials, Wigner d-matrices, rotation calculations, and matrix operations.

cosmocore.basics.legendre_00(scalar_prod, legendre)

In-place computation of normalized Legendre polynomials (2ℓ+1)/(4π) P_l(x).

Parameters:

scalar_prod (float) – Argument x for the Legendre polynomials.
legendre (numpy.ndarray) – Pre-allocated buffer of length lmax + 1. On output, legendre[ell] contains the normalized P_ℓ(x) for ℓ = 0..lmax. legendre[0] = 1/(4π) (the monopole base case) is now populated so lmin_signal=0 callers (e.g. foreground templates) sum the P_0 contribution correctly.

Notes

Computes the recurrence for P_l(x) and absorbs the (2ℓ+1)/(4π) normalization factor, so callers receive the fully normalized basis functions ready for signal matrix construction.

cosmocore.basics.legendre_22(scalar_prod, legendre, f1, f2)

In-place computation of normalized spin-2 Legendre functions.

Computes (2ℓ+1)/(4π) × factor2 × P_l^{22}(x) and related f1, f2 arrays, where factor2 = 1/((ℓ+2)(ℓ+1)ℓ(ℓ-1)).

Parameters:

scalar_prod (float) – Argument x for the associated Legendre polynomials.
legendre (numpy.ndarray) – Pre-allocated array to fill with normalized values.
f1 (numpy.ndarray) – Pre-allocated arrays for the EE/BB decomposition factors.
f2 (numpy.ndarray) – Pre-allocated arrays for the EE/BB decomposition factors.

Notes

Memory-efficient version for spin-2 computations. Fills the provided arrays in-place to avoid allocations in performance-critical loops. Used for polarization auto-correlation calculations. Buffers are ℓ-indexed: legendre[ell], f1[ell], f2[ell] correspond to multipole ℓ. Indices 0 and 1 are zero by spin-2 physics.

cosmocore.basics.legendre_02(scalar_prod, legendre)

In-place computation of normalized spin-0×2 Legendre functions.

Computes (2ℓ+1)/(4π) × factor × P_l^{02}(x), where factor = 1/sqrt((ℓ+2)(ℓ+1)ℓ(ℓ-1)).

Parameters:

scalar_prod (float) – Argument x for the associated Legendre polynomials.
legendre (numpy.ndarray) – Pre-allocated array to fill with normalized values (will be zeroed first).

Notes

Memory-efficient version that avoids allocations in hot loops. Used for temperature-polarization cross-correlations. Buffer is ℓ-indexed: legendre[ell] corresponds to multipole ℓ. Indices 0 and 1 are zero (spin-0×spin-2 requires ℓ ≥ 2).

cosmocore.basics.legendre_plm(cos_theta, sin_theta, plm)

In-place computation of fully normalized associated Legendre polynomials.

Computes sqrt((2ℓ+1)/(4π)) × N_l^m(x) for all l = 0..lmax and m = 0..l, where N_l^m = sqrt((l-m)!/(l+m)!) × P_l^m(x) and x = cos(theta).

Parameters:

cos_theta (float) – Cosine of colatitude angle.
sin_theta (float) – Sine of colatitude angle (must be non-negative).
plm (numpy.ndarray) – Pre-allocated array of shape (lmax+1, lmax+1) to fill with values. On output, plm[ell, m] contains the fully normalized value for m <= l. Values for m > l are set to zero.

cosmocore.basics.wigner_d_small(ell, m, s, cos_theta, sin_theta)

Compute Wigner d-matrix element d^l_{m,s}(theta) for spin-weighted spherical harmonics.

Uses stable recurrence in l derived from the Jacobi polynomial three-term recurrence.

Parameters:

ell (int) – Multipole degree (l >= max(|m|, |s|)).
m (int) – Azimuthal order (-l <= m <= l).
s (int) – Spin weight (typically +/-2 for polarization).
cos_theta (float) – Cosine of colatitude angle.
sin_theta (float) – Sine of colatitude angle (must be >= 0).

Returns:

Value of d^l_{m,s}(theta).

Return type:

float

cosmocore.basics.wigner_d_matrix(ell, s, cos_theta, sin_theta, d_out)

Compute all d^l_{m,s}(theta) for m = -l to l for a fixed l and s.

Parameters:

ell (int) – Multipole degree.
s (int) – Spin weight (typically +/-2 for polarization).
cos_theta (float) – Cosine of colatitude angle.
sin_theta (float) – Sine of colatitude angle.
d_out (numpy.ndarray) – Pre-allocated output array of length (2*l+1) to store d^l_{m,s} for m = -l, -l+1, …, l-1, l. Index i corresponds to m = i - l.

Return type:

None

cosmocore.basics.spec2idx(i, j, nfields)[source]

Convert field indices to spectrum index for compressed storage.

Parameters:

i (int) – Field indices for the spectrum.
j (int) – Field indices for the spectrum.
nfields (int) – Total number of fields.

Returns:

Linear index for spectrum storage in compressed format.

Return type:

int

Notes

Auto-spectra (i==j) are stored first, followed by cross-spectra in upper triangular order. Uses LRU cache for performance.

cosmocore.basics.idx2spec(idx, nfields)[source]

Convert spectrum index back to field indices.

Parameters:

idx (int) – Linear spectrum index in compressed storage.
nfields (int) – Total number of fields.

Returns:

Field indices (i, j) corresponding to the spectrum index.

Return type:

tuple of int

Raises:

ValueError – If index is out of bounds for the given number of fields.

Notes

Inverse operation of spec2idx. Uses LRU cache for performance.

cosmocore.basics.get_rotation_angle(r1, r2)

Compute rotation angles between two 3D vectors on the sphere.

Parameters:

r1 (numpy.ndarray) – 3D unit vectors pointing to different positions on the sphere.
r2 (numpy.ndarray) – 3D unit vectors pointing to different positions on the sphere.

Returns:

Two rotation angles (a12, a21) needed for coordinate transformations between the local coordinate systems at the two positions.

Return type:

tuple of float

Notes

This function computes the rotation angles needed for transforming spin-2 quantities (like polarization) between different coordinate systems on the sphere.

cosmocore.basics.matrix_mult(A, B)[source]

Matrix multiplication wrapper.

Parameters:

A (numpy.ndarray) – 2D matrices to multiply.
B (numpy.ndarray) – 2D matrices to multiply.

Returns:

Result of matrix multiplication A @ B.

Return type:

numpy.ndarray

cosmocore.basics.add_diagonal(M, d, out=None)

Add diagonal vector to matrix: result = M + diag(d).

Parameters:

M (numpy.ndarray) – 2D square matrix.
d (numpy.ndarray) – 1D vector to add to diagonal.
out (numpy.ndarray, optional) – Output array. If None, a new array is created. If provided, must have same shape as M.

Returns:

M + diag(d), computed without creating full diagonal matrix.

Return type:

numpy.ndarray

cosmocore.basics.matrix_trace(A, B)

Compute trace of matrix product A @ B.

Parameters:

A (numpy.ndarray) – 2D square matrices of the same size.
B (numpy.ndarray) – 2D square matrices of the same size.

Returns:

Trace of the matrix product: tr(A @ B) = sum_i sum_j A_ij B_ji.

Return type:

float

cosmocore.basics.matrix_inverse_symm(M, overwrite=False)[source]

Compute inverse of symmetric positive definite matrix.

Parameters:

M (numpy.ndarray) – 2D square symmetric positive definite matrix to invert.
overwrite (bool, optional) – If True, the input matrix M may be overwritten with intermediate results for better performance. Default is False.

Returns:

Inverse of the input matrix.

Return type:

numpy.ndarray

Raises:

ValueError – If matrix is not square or if Cholesky decomposition fails even after symmetric-scaling preconditioning.

Notes

Uses LAPACK’s dpotrf (Cholesky) and dpotri (inverse) for efficient and numerically stable inversion of symmetric positive definite matrices.

The matrix is always symmetrically preconditioned before Cholesky as D M D with D = diag(1/sqrt(|M_ii|)), then unscaled back after inversion via F^{-1} = D (D M D)^{-1} D. This recovers PD-ness for matrices that are mathematically positive definite but numerically ill-conditioned because of large dynamic range across the diagonal — the typical case of multi-spectrum QML Fisher matrices that mix cosmic-variance-limited (huge signal, tiny Fisher entry) and noise- limited (large Fisher entry) bandpowers. The preconditioning is cheap (two diagonal scalings) and well-conditioned matrices are unaffected in finite precision.

cosmocore.basics.matrix_slogdet(M)[source]

Compute sign and logarithm of the determinant of a matrix.

Parameters:: M (numpy.ndarray) – 2D square matrix for which to compute the signed log determinant.
Returns:: (sign, logdet) where sign is +/-1 and logdet is log(|det(M)|). If det(M) = 0, returns (0.0, -inf).
Return type:: tuple of float
Raises:: ValueError – If matrix is not square or if LU decomposition fails.

Notes

Uses LAPACK’s dgetrf (LU decomposition with partial pivoting).

cosmocore.basics.cholesky_decomposition(M)[source]

Perform Cholesky decomposition of a symmetric positive definite matrix.

Parameters:: M (numpy.ndarray) – 2D square symmetric positive definite matrix to decompose.
Returns:: Lower triangular matrix L such that M = L @ L.T.
Return type:: numpy.ndarray
Raises:: ValueError – If matrix is not square or if Cholesky decomposition fails.

cosmocore.basics.cholesky_factor(M, overwrite_a=False, clean=True)[source]

Cholesky factor (L, lower=True) of a symmetric positive definite matrix.

Pairs with cholesky_solve(). With clean=True (default) the upper triangle is zeroed so consumers that treat the factor as a dense matrix (e.g., L @ L.T) get correct results without reading garbage.

Parameters:

M (numpy.ndarray) – 2D square symmetric positive definite matrix.
overwrite_a (bool, optional) – If True, M is overwritten with the factor in place (no extra n×n allocation). Default is False.
clean (bool, optional) – If True, zero the upper triangle of the result. Default is True.

Returns:

(c, True) where c is the n×n array holding L in its lower triangle. Suitable as the first argument to cholesky_solve().

Return type:

tuple

Raises:

numpy.linalg.LinAlgError – If M is not positive definite.
ValueError – If overwrite_a=True but M is not F-contiguous (LAPACK would silently ignore overwrite_a and copy, defeating the in-place intent).

cosmocore.basics.cholesky_solve(c_and_lower, b, overwrite_b=False)[source]

Solve A x = b given the Cholesky factor of A from cholesky_factor().

Parameters:

c_and_lower (tuple) – (c, lower) tuple as returned by cholesky_factor().
b (numpy.ndarray) – Right-hand side, shape (n,) or (n, k).
overwrite_b (bool, optional) – If True, b may be overwritten with the solution. Default is False.

Returns:

Solution x with the same shape as b.

Return type:

numpy.ndarray

cosmocore.basics.matrix_slogdet_symm(M)[source]

Compute sign and logarithm of determinant for symmetric positive definite matrix.

Parameters:: M (numpy.ndarray) – 2D square symmetric positive definite matrix.
Returns:: (sign, logdet) where sign is +1 and logdet is log(det(M)). For positive definite matrices, sign is always +1.
Return type:: tuple of float
Raises:: ValueError – If matrix is not square or if Cholesky decomposition fails.

Notes

Uses Cholesky decomposition: det(M) = det(L)^2 = (prod L_ii)^2.

cosmocore.basics.symmetrize_inplace(M)

Symmetrize M in place: M[i,j] = M[j,i] = (M[i,j] + M[j,i]) / 2.

Equivalent to M = 0.5 * (M + M.T) but without the two n² temporary allocations of the broadcast form. Important on the QML hot path where M is the SMW kernel V N⁻¹ Vᵀ (n_modes × n_modes); at production scale the temporaries can dominate the basis-setup memory peak.

Parameters:: M (numpy.ndarray) – 2D square matrix.
Returns:: The same array, now symmetric (lower and upper triangles agree).
Return type:: numpy.ndarray

cosmocore.basics.inv(A)[source]

General LU-based matrix inverse.

For symmetric positive-definite matrices, prefer matrix_inverse_symm() (Cholesky, ~2× faster and more stable).

cosmocore.basics.solve_linear(A, b)[source]

General LU-based linear solve A x = b.

For symmetric positive-definite systems, prefer the cholesky_factor() + cholesky_solve() pair.

cosmocore.basics.eigh(A)[source]

Eigendecomposition of a symmetric / Hermitian matrix.

Returns (eigenvalues, eigenvectors) with eigenvalues in ascending order.

cosmocore.basics.eigvalsh(A)[source]: Eigenvalues of a symmetric / Hermitian matrix, ascending order.

cosmocore.basics.svd(A, full_matrices=False)[source]

Singular value decomposition with full_matrices=False by default.

The thin SVD is the right default for the orthogonalization use cases in this codebase; the wrapper exists so the choice is uniform.

cosmocore.basics.smw_inverse(N_inv, V_N_inv, V_Ninv_VT, Lambda_diag, threshold=1e-30)[source]

Sherman-Morrison-Woodbury inverse: (N + V^T Lambda V)^{-1}.

Parameters:

N_inv (numpy.ndarray) – Precomputed inverse of N, shape (n, n).
V_N_inv (numpy.ndarray) – Precomputed V @ N^{-1}, shape (k, n).
V_Ninv_VT (numpy.ndarray) – Precomputed V @ N^{-1} @ V^T, shape (k, k).
Lambda_diag (numpy.ndarray) – Diagonal of Lambda, shape (k,).
threshold (float, optional) – Minimum value for Lambda diagonal elements. Default is 1e-30.

Returns:

The inverse (N + V^T Lambda V)^{-1}, shape (n, n).

Return type:

numpy.ndarray

cosmocore.basics.smw_logdet(log_det_N, V_Ninv_VT, Lambda_diag, threshold=1e-30)[source]

Sherman-Morrison-Woodbury log determinant: log|N + V^T Lambda V|.

Parameters:

log_det_N (float) – Precomputed log|N|.
V_Ninv_VT (numpy.ndarray) – Precomputed V @ N^{-1} @ V^T, shape (k, k).
Lambda_diag (numpy.ndarray) – Diagonal of Lambda, shape (k,).
threshold (float, optional) – Minimum value for Lambda diagonal elements. Default is 1e-30.

Returns:

The log determinant log|N + V^T Lambda V|.

Return type:

float

cosmocore.basics.smw_kernel(V_Ninv_VT, Lambda_diag, threshold=1e-30)[source]

Build the SMW kernel matrix K = Lambda^{-1} + V N^{-1} V^T.

Parameters:

V_Ninv_VT (numpy.ndarray) – Precomputed V @ N^{-1} @ V^T, shape (k, k).
Lambda_diag (numpy.ndarray) – Diagonal of Lambda, shape (k,).
threshold (float, optional) – Minimum value for Lambda diagonal elements. Default is 1e-30.

Returns:

The kernel matrix K, shape (k, k), in Fortran order.

Return type:

numpy.ndarray

cosmocore.basics.smw_quadratic_form(data, N_chol, V_N_inv, V_Ninv_VT, Lambda_diag, threshold=1e-30)[source]

SMW quadratic form: d^T (N + V^T Lambda V)^{-1} d.

Parameters:

data (numpy.ndarray) – Data vector d, shape (n,).
N_chol (tuple) – Cholesky factor (L, lower=True) of N as returned by cholesky_factor().
V_N_inv (numpy.ndarray) – Precomputed V @ N^{-1}, shape (k, n).
V_Ninv_VT (numpy.ndarray) – Precomputed V @ N^{-1} @ V^T, shape (k, k).
Lambda_diag (numpy.ndarray) – Diagonal of Lambda, shape (k,).
threshold (float, optional) – Minimum value for Lambda diagonal elements. Default is 1e-30.

Returns:

The quadratic form d^T C^{-1} d.

Return type:

float