foscat.alm_loc_optim ==================== .. py:module:: foscat.alm_loc_optim .. autoapi-nested-parse:: alm_loc_optim.py ================ 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 ------- .. autoapisummary:: foscat.alm_loc_optim.alm_loc_optim Module Contents --------------- .. py:class:: alm_loc_optim(backend=None, lmax=24, limit_range=10000000000.0) Bases: :py:obj:`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 .. py:method:: analyze_domain(nside: int, cell_ids, nest: bool = False, lmax: int = None) Pre-compute the set of (l, m) modes that are effectively constrained by the partial patch described by cell_ids. :Parameters: * **nside** (:py:class:`HEALPix resolution`) * **cell_ids** (:py:class:`indices de pixels (ring ou nested)`) * **nest** (:py:class:`True si cell_ids utilise l'ordre NESTED`) * **lmax** (:py:class:`maximum multipole (default`: :py:class:`min(self.lmax`, :py:class:`3*nside-1))`) :returns: * **rings_used** (:py:class:`ndarray int32 [R] indices` of :py:class:`present rings`) * **counts** (:py:class:`ndarray int32 [R] nb de pixels par ring`) * **sizes** (:py:class:`ndarray int32 N_ring pour chaque ring du nside`) * **valid_idx_per_m** (:py:class:`dict m -> ndarray int32 indices dans rings_used`) * **lmax_eff_per_m** (:py:class:`dict m -> int lmax effectif pour ce m`) * **m_valid** (:py:class:`list[int] valeurs de m ayant au moins 1 ring valide`) .. py:method:: map2alm_loc_optim(im, nside: int, cell_ids, nest: bool = False, lmax: int = None) Compute spherical-harmonic coefficients on a partial patch, restricted to the (l, m) modes effectively constrained by the partial sky coverage. :Parameters: * **im** (:py:class:`[...`, :py:class:`n_pixels] map values on the patch`) * **nside** (:py:class:`HEALPix resolution`) * **cell_ids** (:py:class:`pixel indices in the patch`) * **nest** (:py:class:`True if NESTED ordering`) * **lmax** (:py:class:`maximum multipole`) :returns: * **alm_per_m** (:py:class:`list` of :py:class:`tensors`, :py:class:`one per valid m.`) -- alm_per_m[i] a la forme [..., lmax_eff(m)-m+1] pour m = m_valid[i] * **m_valid** (:py:class:`list[int] m values in the same order`) * **lmax_eff_per_m** (:py:class:`dict m -> int`) .. py:method:: sparse_to_dense(alm_per_m, m_valid, lmax_eff_per_m, lmax: int) 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** (:py:class:`liste de tenseurs (sortie de map2alm_loc_optim)`) * **m_valid** (:py:class:`list[int]`) * **lmax_eff_per_m** (:py:class:`dict m -> int`) * **lmax** (:py:class:`maximum multipole used during computation`) * **Retourne** * **--------** * **out** (:py:class:`tensor [...`, :py:class:`total_alm] same dtype` and :py:class:`device as input`) .. py:method:: anafast_loc_optim(im, nside: int, cell_ids, nest: bool = False, lmax: int = None) 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** (:py:class:`[...`, :py:class:`n_pixels]`) * **nside** (:py:class:`HEALPix resolution`) * **cell_ids** (:py:class:`pixel indices in the patch`) * **nest** (:py:class:`True if NESTED ordering`) * **lmax** (:py:class:`maximum multipole`) * **Retourne** * **--------** * **cl** (:py:class:`tenseur [...`, :py:class:`lmax+1]`) * **m_count** (:py:class:`ndarray int32 [lmax+1]`) -- number of m modes contributing to each l (flags poorly constrained multipoles: m_count[l] == 0 → Cl[l] undefined) .. py:method:: domain_summary(nside: int, cell_ids, nest: bool = False, lmax: int = None) 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.