lenspyx.experimental

Experimental spherical harmonic transforms with optimizations for spherical caps.

This module provides wrapper functions and classes for spherical harmonic transforms (SHTs) that automatically select optimized implementations based on the geometry of the data.

Overview

The module wraps functions from ducc0 and capsht, providing:

• Automatic selection between general SHT and optimized capped SHT

• Unified interface for different pixel location types (HEALPix, arbitrary locations)

• Support for both synthesis (alm → map) and adjoint synthesis (map → alm)

• Optimized implementations for data confined to spherical caps

When data is confined to a spherical cap (less than half the sky), the capped transforms can provide significant speedups compared to full-sky transforms.

Note

This module requires capsht to be installed for the optimized capped transforms. If capsht is not available, functions will fall back to standard ducc0 implementations.

API Reference

Synthesis Functions

lenspyx.experimental.synthesis_general(alm: ndarray, spin: int, lmax: int, loc: ndarray, epsilon: float, thtcap: float | None = None, eps_apo: float | None = None, tht_min: float | None = None, tht_max: float | None = None, verbose: bool = False, **kwargs)

Spherical harmonic synthesis with automatic optimization for spherical caps.

This function automatically selects the most efficient implementation based on the geometry of the target locations. If data is confined to a spherical cap, it uses an optimized capped transform; otherwise it falls back to the general transform.

Parameters:

• alm (array_like) – Spherical harmonic coefficients. Shape (ncomp, nalm) where ncomp is 1 for spin-0 fields or 2 for spin-s fields with s > 0.

• spin (int) – Spin weight of the field (0 for temperature/scalar, ±2 for polarization)

• lmax (int) – Maximum multipole moment

• loc (array_like) – Evaluation locations, shape (npix, 2) with (colatitude, longitude) in radians

• epsilon (float) – Target accuracy (e.g., 1e-7 for typical CMB applications)

• thtcap (float, optional) – Cap radius in radians. If provided, assumes all locations lie within this cap and uses optimized capped transform. If None, attempts automatic detection.

• eps_apo (float, optional) – Apodization parameter. If None, computed automatically from epsilon and lmax.

• tht_min (float, optional) – Minimum colatitude in the data (currently experimental)

• tht_max (float, optional) – Maximum colatitude in the data (currently experimental)

• verbose (bool, optional) – Print diagnostic information about which transform is used

• **kwargs (dict) –

  Additional keyword arguments passed to the underlying transform:

  • mode : str, optional (e.g., ‘STANDARD’, ‘GRAD_ONLY’)

  • map : array_like, optional (pre-allocated output array)

  • mmax : int, optional (maximum azimuthal mode number)

  • nthreads : int, optional (number of threads)

Returns:

Synthesized map(s). Shape (ncomp, npix) where ncomp matches the input alm.

Return type:

array_like

Notes

The function automatically determines which implementation to use:

• If thtcap is provided and all locations satisfy \(\theta \leq \theta_{\text{cap}}\), uses capsht.synthesis_general_cap() for optimal performance

• Otherwise, uses ducc0.sht.synthesis_general() for arbitrary locations

The capped transform can provide significant speedups for data confined to caps covering less than half the sky.

Examples

>>> import numpy as np
>>> from lenspyx.experimental import synthesis_general
>>> # Synthesize on arbitrary locations
>>> alm = np.random.randn(1, nalm) + 1j * np.random.randn(1, nalm)
>>> loc = np.random.rand(1000, 2) * [np.pi, 2*np.pi]  # Random locations
>>> maps = synthesis_general(alm, spin=0, lmax=100, loc=loc, epsilon=1e-7)

>>> # Use capped transform for polar cap
>>> theta_cap = 30 * np.pi / 180  # 30 degree cap
>>> loc_cap = np.random.rand(1000, 2) * [theta_cap, 2*np.pi]
>>> maps_cap = synthesis_general(alm, spin=0, lmax=100, loc=loc_cap,
...                              epsilon=1e-7, thtcap=theta_cap, verbose=True)

See also

adjoint_synthesis_general

Adjoint operation (map → alm)

ducc0.sht.synthesis_general

Underlying general transform

capsht.synthesis_general_cap

Underlying capped transform

lenspyx.experimental.adjoint_synthesis_general(map: ndarray, spin: int, lmax: int, loc: ndarray, epsilon: float, thtcap: float | None = None, eps_apo: float | None = None, tht_min: float | None = None, tht_max: float | None = None, verbose: bool = False, **kwargs)

Adjoint spherical harmonic synthesis with automatic optimization for spherical caps.

This is the adjoint operation to synthesis_general(), computing spherical harmonic coefficients from map data. It automatically selects the most efficient implementation based on the geometry.

Parameters:

• map (array_like) – Input map(s). Shape (ncomp, npix) where ncomp is 1 for spin-0 or 2 for spin-s.

• spin (int) – Spin weight of the field (0 for temperature/scalar, ±2 for polarization)

• lmax (int) – Maximum multipole moment

• loc (array_like) – Map pixel locations, shape (npix, 2) with (colatitude, longitude) in radians

• epsilon (float) – Target accuracy (e.g., 1e-7 for typical CMB applications)

• thtcap (float, optional) – Cap radius in radians. If provided, uses optimized capped transform.

• eps_apo (float, optional) – Apodization parameter. If None, computed automatically.

• tht_min (float, optional) – Minimum colatitude in the data (currently experimental)

• tht_max (float, optional) – Maximum colatitude in the data (currently experimental)

• verbose (bool, optional) – Print diagnostic information about which transform is used

• **kwargs (dict) –

  Additional keyword arguments:

  • mode : str, optional (e.g., ‘STANDARD’, ‘GRAD_ONLY’)

  • alm : array_like, optional (pre-allocated output array)

  • mmax : int, optional (maximum azimuthal mode number)

  • nthreads : int, optional (number of threads)

Returns:

Spherical harmonic coefficients. Shape (ncomp, nalm) where nalm depends on lmax and mmax.

Return type:

array_like

Notes

This function computes the adjoint operation (not the inverse!) of the synthesis transform. For a forward transform \(\mathbf{m} = \mathbf{S} \mathbf{a}\), this computes \(\mathbf{a}' = \mathbf{S}^{\dagger} \mathbf{m}\) where \(\dagger\) denotes the adjoint (conjugate transpose).

The adjoint is used in:

• Quadratic estimators for CMB lensing and other fields

• Iterative map-making and Wiener filtering

• Maximum likelihood power spectrum estimation

Like synthesis_general(), this automatically selects between general and capped implementations based on the data geometry.

Examples

>>> import numpy as np
>>> from lenspyx.experimental import adjoint_synthesis_general
>>> from lenspyx.utils_hp import Alm
>>> # Compute alm from map on arbitrary locations
>>> maps = np.random.randn(1, 1000)  # Random map
>>> loc = np.random.rand(1000, 2) * [np.pi, 2*np.pi]
>>> lmax = 100
>>> alm = np.zeros((1, Alm.getsize(lmax, lmax)), dtype=complex)
>>> alm = adjoint_synthesis_general(maps, spin=0, lmax=lmax, loc=loc,
...                                 epsilon=1e-7, alm=alm)

See also

synthesis_general

Forward transform (alm → map)

ducc0.sht.adjoint_synthesis_general

Underlying general adjoint transform

capsht.adjoint_synthesis_general_cap

Underlying capped adjoint transform

Locations Class

class lenspyx.experimental.Locations(loc: ndarray | None = None, geom=None, epsilon=1e-07)

Bases: object

Unified interface for locations (or pixelizations) in spherical harmonic transforms.

This class provides a common API for working with both structured (iso-latitude rings) and unstructured (arbitrary locations) pixelizations. It automatically selects the appropriate transform implementation based on the location type.

Parameters:

• loc (array_like, optional) – Arbitrary pixel locations, shape (npix, 2) with (colatitude, longitude) in radians. Use this for unstructured grids or partial sky coverage.

• geom (lenspyx.remapping.utils_geom.Geom, optional) – Iso-latitude ring geometry (e.g., HEALPix, Gauss-Legendre). Use this for structured, symmetric pixelizations.

• epsilon (float, optional) – Target accuracy for transforms (default: 1e-7). Only relevant when using loc.

Notes

Exactly one of loc or geom must be provided.

This class serves as a wrapper that:

• Provides consistent synthesis() and adjoint_synthesis() methods

• Automatically routes to optimized implementations (capped vs general, structured vs unstructured)

• Handles differences in calling conventions between backends

loc

Pixel locations if using unstructured grid

Type:

array_like or None

geom

Geometry object if using structured grid

Type:

Geom or None

thtrange

(min_theta, max_theta) colatitude range

Type:

tuple

epsilon

Transform accuracy parameter

Type:

float

Examples

Using arbitrary locations:

>>> import numpy as np
>>> from lenspyx.experimental import Locations
>>> # Create random pixel locations in a polar cap
>>> npix = 10000
>>> loc = np.random.rand(npix, 2) * [np.pi/4, 2*np.pi]
>>> locs = Locations(loc=loc, epsilon=1e-7)
>>> # Synthesize a map
>>> alm = np.random.randn(1, nalm) + 1j * np.random.randn(1, nalm)
>>> maps = locs.synthesis(alm, spin=0, lmax=100, mmax=100, nthreads=4)

Using HEALPix geometry:

>>> from lenspyx.remapping.utils_geom import Geom
>>> geom = Geom.get_healpix_geometry(nside=512)
>>> locs = Locations(geom=geom)
>>> maps = locs.synthesis(alm, spin=0, lmax=1500, mmax=1500, nthreads=4)

See also

lenspyx.remapping.utils_geom.Geom

Geometry objects for structured grids

synthesis_general

Underlying synthesis function

adjoint_synthesis_general

Underlying adjoint synthesis function

adjoint_synthesis(alm: ndarray, spin: int, lmax: int, mmax: int, nthreads: int, m: ndarray, **kwargs)

Compute adjoint spherical harmonic synthesis (map → alm).

Transforms a map to spherical harmonic coefficients using the adjoint operation. This is NOT the inverse transform, but rather the conjugate transpose.

Parameters:

• alm (array_like) – Pre-allocated output array for coefficients, shape (ncomp, nalm)

• spin (int) – Spin weight of the field

• lmax (int) – Maximum multipole moment

• mmax (int) – Maximum azimuthal mode number

• nthreads (int) – Number of threads to use for the transform

• m (array_like) –

  Input map array, shape (ncomp, npix) where:

  • ncomp = 1 for spin-0 fields

  • ncomp = 2 for spin-s fields with s > 0

• **kwargs (dict) – Additional arguments passed to the underlying transform

Returns:

Spherical harmonic coefficients, shape (ncomp, nalm)

Return type:

array_like

Notes

For a forward transform \(\mathbf{m} = \mathbf{S} \mathbf{a}\), this computes \(\mathbf{a}' = \mathbf{S}^{\dagger} \mathbf{m}\) where \(\dagger\) is the adjoint (conjugate transpose), not the inverse.

The adjoint is used in quadratic estimators, iterative algorithms, and maximum likelihood estimation.

Examples

>>> locs = Locations(loc=my_locations)
>>> alm_out = np.zeros((1, nalm), dtype=complex)
>>> alm_out = locs.adjoint_synthesis(alm_out, spin=0, lmax=2000, mmax=2000,
...                                  nthreads=8, m=input_map)

npix()

synthesis(alm: ndarray, spin: int, lmax: int, mmax: int, nthreads: int, m: ndarray | None = None, **kwargs)

Compute spherical harmonic synthesis (alm → map).

Transforms spherical harmonic coefficients to a map on the pixelization.

Parameters:

• alm (array_like) –

  Spherical harmonic coefficients, shape (ncomp, nalm) where:

  • ncomp = 1 for spin-0 fields (temperature, scalar)

  • ncomp = 2 for spin-s fields with s > 0 (gradient and curl components)

• spin (int) – Spin weight of the field

• lmax (int) – Maximum multipole moment in the alm array

• mmax (int) – Maximum azimuthal mode number in the alm array

• nthreads (int) – Number of threads to use for the transform

• m (array_like, optional) – Pre-allocated output array, shape (ncomp, npix). If None, allocated automatically.

• **kwargs (dict) – Additional arguments passed to the underlying transform (e.g., mode=’GRAD_ONLY’)

Returns:

Map array of shape (ncomp, npix)

Return type:

array_like

Examples

>>> locs = Locations(loc=my_locations)
>>> maps = locs.synthesis(alm, spin=2, lmax=2000, mmax=2000, nthreads=8)

Examples

Using synthesis_general with automatic cap detection

import numpy as np
from lenspyx.experimental import synthesis_general
from lenspyx.utils_hp import Alm

# Create random alm coefficients
lmax = 2000
nalm = Alm.getsize(lmax, lmax)
alm = np.random.randn(1, nalm) + 1j * np.random.randn(1, nalm)

# Define locations in a polar cap (30 degrees)
npix = 10000
theta_cap = 30 * np.pi / 180
loc = np.random.rand(npix, 2) * [theta_cap, 2*np.pi]

# Synthesize with automatic cap detection
maps = synthesis_general(alm, spin=0, lmax=lmax, loc=loc,
                        epsilon=1e-7, thtcap=theta_cap, verbose=True)

Using Locations class with arbitrary locations

import numpy as np
from lenspyx.experimental import Locations
from lenspyx.utils_hp import Alm

# Create Locations object from arbitrary pixel positions
npix = 10000
loc = np.random.rand(npix, 2) * [np.pi, 2*np.pi]
locs = Locations(loc=loc, epsilon=1e-7)

# Forward transform (alm → map)
lmax = 1000
nalm = Alm.getsize(lmax, lmax)
alm = np.random.randn(1, nalm) + 1j * np.random.randn(1, nalm)
maps = locs.synthesis(alm, spin=0, lmax=lmax, mmax=lmax, nthreads=4)

# Adjoint transform (map → alm)
alm_out = np.zeros((1, nalm), dtype=complex)
alm_out = locs.adjoint_synthesis(alm_out, spin=0, lmax=lmax,
                                 mmax=lmax, nthreads=4, m=maps)

Using Locations with HEALPix geometry

from lenspyx.experimental import Locations
from lenspyx.remapping.utils_geom import Geom
from lenspyx.utils_hp import Alm
import numpy as np

# Create HEALPix geometry
nside = 512
geom = Geom.get_healpix_geometry(nside)
locs = Locations(geom=geom)

# Synthesize on HEALPix grid
lmax = 1500
nalm = Alm.getsize(lmax, lmax)
alm = np.random.randn(1, nalm) + 1j * np.random.randn(1, nalm)
maps = locs.synthesis(alm, spin=0, lmax=lmax, mmax=lmax, nthreads=8)

Polarization (spin-2) example

import numpy as np
from lenspyx.experimental import synthesis_general
from lenspyx.utils_hp import Alm

# Create E and B mode coefficients
lmax = 2000
nalm = Alm.getsize(lmax, lmax)
alm_eb = np.random.randn(2, nalm) + 1j * np.random.randn(2, nalm)

# Define locations
npix = 5000
loc = np.random.rand(npix, 2) * [np.pi/2, 2*np.pi]

# Synthesize Q and U maps (spin=2)
maps_qu = synthesis_general(alm_eb, spin=2, lmax=lmax, loc=loc,
                            epsilon=1e-7, nthreads=8)
# maps_qu has shape (2, npix) containing [Q, U]

Notes

• The capped transforms are most beneficial when data covers much less than ~50% of the sky

• The Locations class provides a convenient unified interface regardless of geometry type

• Use verbose=True to see which transform implementation is selected

References

[1]

Carron, J. and Reinecke, M., 2026. “Fast partial-sky spherical harmonic transforms.” arXiv:2603.17166. https://arxiv.org/abs/2603.17166

[2]

Reinecke, M., Belkner, S., and Carron, J., 2023. “Improved cosmic microwave background (de-)lensing using general spherical harmonic transforms.” arXiv:2304.10431. https://arxiv.org/abs/2304.10431

See Also

ducc0.sht : Standard spherical harmonic transforms capsht : Capped spherical harmonic transforms library lenspyx.remapping.utils_geom : Geometry objects for structured grids