foscat.alm_loc_optim#

Optimised variant of alm_loc for spherical-harmonic computation on a restricted HEALPix domain.

Two key optimisations relative to alm_loc#

  1. Restriction in m (longitude axis) For each ring r with cnt pixels out of N_ring total:

    • full ring : mode m valid if (m % N_ring) <= N_ring // 2

    • partial ring : mode m valid if (m % N_ring) <= cnt // 2 (partial Nyquist)

    For each m, only rings satisfying this criterion feed the Legendre projection.

  2. Restriction in l (latitude axis) With R_m valid rings for mode m, the Legendre system has R_m equations. At most R_m ell degrees can be constrained. Hence:

    lmax_eff(m) = min(lmax, m + R_m - 1)

    and Legendre polynomials are computed only up to this ceiling.

Resulting savings#

La projection de Legendre passe de O(R × (lmax - m + 1)) to O(R_m × (lmax_eff(m) - m + 1)), which is quadratically more favorable pour les petits domaines.

Public API (mirroring alm_loc)#

analyze_domain(nside, cell_ids, nest, lmax)

→ rings_used, counts, sizes, valid_idx_per_m, lmax_eff_per_m, m_valid

map2alm_loc_optim(im, nside, cell_ids, nest, lmax)

→ alm_per_m, m_valid, lmax_eff_per_m (sparse alm)

anafast_loc_optim(im, nside, cell_ids, nest, lmax)

→ cl [lmax+1], m_count [lmax+1]

sparse_to_dense(alm_per_m, m_valid, lmax_eff_per_m, lmax)
→ dense alm vector (zeros for uncalculated modes),

compatible avec la mise en page de alm_loc.map2alm_loc

Classes#

alm_loc_optim

Local/partial-sky variant of foscat.alm.alm.

Module Contents#

class foscat.alm_loc_optim.alm_loc_optim(backend=None, lmax=24, limit_range=10000000000.0)[source]#

Bases: foscat.alm_loc.alm_loc

Local/partial-sky variant of foscat.alm.alm.

Key design choice (to match alm.py exactly when full-sky is provided): - Reuse all Legendre/normalization machinery from the parent class (alm),

i.e. shift_ph(), compute_legendre_m(), ratio_mm, A/B recurrences, etc. This is critical for matching alm.map2alm() numerically.

Differences vs alm.map2alm(): - Input map is […, n] with explicit (nside, cell_ids) - Only rings touched by cell_ids are processed. - For rings with full coverage, we run the exact same FFT+tiling logic as alm.comp_tf()

(but only for those rings) -> bitwise comparable up to backend FFT differences.

  • For rings with partial coverage, we compute a partial DFT for m=0..mmax, using the same phase convention as alm.comp_tf():

    FFT kernel uses exp(-i 2pi (m mod Nring) j / Nring) then apply the per-ring shift exp(-i m phi0) via self.matrix_shift_ph

analyze_domain(nside: int, cell_ids, nest: bool = False, lmax: int = None)[source]#

Pre-compute the set of (l, m) modes that are effectively constrained by the partial patch described by cell_ids.

Parameters:

  • nside (HEALPix resolution)

  • cell_ids (indices de pixels (ring ou nested))

  • nest (True si cell_ids utilise l'ordre NESTED)

  • lmax (maximum multipole (default: min(self.lmax, 3*nside-1)))

Returns:

  • rings_used (ndarray int32 [R]    indices of present rings)

  • counts (ndarray int32 [R]    nb de pixels par ring)

  • sizes (ndarray int32         N_ring pour chaque ring du nside)

  • valid_idx_per_m (dict  m -> ndarray int32  indices dans rings_used)

  • lmax_eff_per_m (dict  m -> int   lmax effectif pour ce m)

  • m_valid (list[int]  valeurs de m ayant au moins 1 ring valide)

map2alm_loc_optim(im, nside: int, cell_ids, nest: bool = False, lmax: int = None)[source]#

Compute spherical-harmonic coefficients on a partial patch, restricted to the (l, m) modes effectively constrained by the partial sky coverage.

Parameters:

  • im ([..., n_pixels]  map values on the patch)

  • nside (HEALPix resolution)

  • cell_ids (pixel indices in the patch)

  • nest (True if NESTED ordering)

  • lmax (maximum multipole)

Returns:

  • alm_per_m (list of tensors, one per valid m.) – alm_per_m[i] a la forme […, lmax_eff(m)-m+1] pour m = m_valid[i]

  • m_valid (list[int]   m values in the same order)

  • lmax_eff_per_m (dict  m -> int)

sparse_to_dense(alm_per_m, m_valid, lmax_eff_per_m, lmax: int)[source]#

Convertit l’alm creux (sortie de map2alm_loc_optim) vers le vecteur flat dense array used by alm_loc.map2alm_loc.

Format dense : [m=0: l=0..lmax | m=1: l=1..lmax | …] Uncalculated modes are filled with zeros.

Parameters:

  • alm_per_m (liste de tenseurs (sortie de map2alm_loc_optim))

  • m_valid (list[int])

  • lmax_eff_per_m (dict m -> int)

  • lmax (maximum multipole used during computation)

  • Retourne

  • ——–

  • out (tensor [..., total_alm]  same dtype and device as input)

anafast_loc_optim(im, nside: int, cell_ids, nest: bool = False, lmax: int = None)[source]#

Estimate the angular power spectrum Cl on a partial patch.

Only the (l, m) modes effectively constrained by the patch contribute to the estimate. For each l, Cl is normalised by the number of available m modes (graceful degradation rather than dilution by zero modes).

Parameters:

  • im ([..., n_pixels])

  • nside (HEALPix resolution)

  • cell_ids (pixel indices in the patch)

  • nest (True if NESTED ordering)

  • lmax (maximum multipole)

  • Retourne

  • ——–

  • cl (tenseur [..., lmax+1])

  • m_count (ndarray int32 [lmax+1]) – number of m modes contributing to each l (flags poorly constrained multipoles:

    m_count[l] == 0 → Cl[l] undefined)

domain_summary(nside: int, cell_ids, nest: bool = False, lmax: int = None)[source]#

Print a human-readable summary of the partial domain: sky fraction covered, number of rings, effective (m, l) mode count vs brute-force, and estimated computation gain.