\(h\)-maps#

Overview#

The h_maps module provides tools to generate complex harmonic maps \(h_{n,m}\) from LiteBIRD observations. These maps allow to perform a spin-based mapmaking, they encode the scanning strategy information. The spin-based mapmaking can be performed with the SMARTIES package, it has an interface with LBS \(h\) maps through the functions read_lbs_h_maps() and build_mask_from_lbs_h_maps().

This implementation follows the formalism introduced in McCallum et al (2021) (arXiv:2109.05038), extended to the HWP case in Takase et al (2024) (arXiv:2408.03040).

Mathematical formalism#

For a scanning strategy with HWP modulation, the harmonic map of spin order \(n,m\) at pixel \(p\) is defined as:

(1)#\[h_{n,m}(p) = \frac{1}{N_p} \sum_{t \in p} e^{i n \psi_t + i m \chi_t}\]

where:

  • \(N_p\) is the number of observations (hits) in pixel \(p\)

  • \(\psi_t\) is the parralactic angle of the detector at time sample \(t\)

  • \(\chi_t\) is the HWP rotation angle at time \(t\)

  • The sum runs over all time samples \(t\) that fall within pixel \(p\)

These maps capture how the detector orientation is distributed across the sky during the scanning strategy. Setting \(m=0\) is equivalent to the definition of h maps in McCallum et al (2021) without HWP modulation.

The maps can be generated by using the make_h_maps() function with a list of observations, or by using the observations of a simulation through the method make_h_maps() of simulation

Example#

import litebird_sim as lbs
import numpy as np
from astropy.time import Time

# load an IMO
imo = lbs.Imo()
# Create a simple simulation
sim = lbs.Simulation(
 base_path="./output_sim",
 start_time=0,
 duration_s=24 * 3600.0,
 random_seed=12345,
 imo=imo,
)

# ... (set up detectors, scanning strategy, etc.) ...

#Create an observation

sim.create_observations(
 detectors=dets,
 n_blocks_det=1,
 n_blocks_time=1,
 tod_dtype=precision,
 split_list_over_processes=False,
 allocate_tod=False,    #We only need the pointing information to compute h maps, so we can save memory by not allocating the full TOD
)
sim.prepare_pointings()

result = lbs.make_h_maps(
    observations=sim.observations,
    nside=64,
    n_m_couples=np.array([[0, 0], [2, 0], [4, 0]]),
    output_coordinate_system=lbs.CoordinateSystem.Galactic,
    save_to_file=True,
    output_directory="./h_n_maps",
)

The n_m_couples parameter#

The spin orders to compute are specified as a 2D array of (n, m) pairs:

# Compute h_{0,0}, h_{2,0} and h_{4,0}
n_m_couples = np.array([
    [0, 0],
    [2, 0],
    [4, 0],
])

Note: for the couple (0,0) the resulting map is not the one expected from the definition of eq (1). This definition yield \(h_{0,0} = \frac{1}{N_p} \sum_{t \in p} e^{i 0 \psi_t + i 0 \chi_t} = 1\) for all observed pixels, which is not very useful. Instead the \(h_{0,0}\) is an hitmap, i.e. \(h_{0,0}(p) = N_p\).

Output: HMapsResult#

make_h_maps() returns a HMapsResult object, which contains:

  • h_maps: a dictionary indexed by detector name, then by (n, m) tuple, each entry being a h_map_Re_and_Im object.

  • coordinate_system: the output coordinate system (a CoordinateSystem object).

  • detector_split: the detector split used.

  • time_split: the time split used.

  • duration_s: the duration of the observation in seconds.

  • sampling_rate_Hz: the sampling rate of the observation in Hz.

Accessing individual maps#

Each entry in h_maps is a h_map_Re_and_Im object with:

  • real: the real part of \(h_{n,m}\) (NumPy array of shape (Npix,))

  • imag: the imaginary part (NumPy array of shape (Npix,))

  • norm: property returning \(\sqrt{\mathrm{Re}^2 + \mathrm{Im}^2}\), with unobserved pixels set to UNSEEN_PIXEL_VALUE

  • n, m: the spin orders

  • det_info: the detector name

Example:

h_2_0 = result.h_maps["detector_name"][2, 0]
print(h_2_0.real)   # real part
print(h_2_0.imag)   # imaginary part
print(h_2_0.norm)   # amplitude

Saving and loading maps#

Maps are saved in HDF5 format, one file per detector:

output_directory/
└── h_maps_det_{detector_name}.h5
    ├── attrs: coordinate_system, det, detector_split,
    │          time_split, duration, sampling_rate
    ├── 0,0/
    │   ├── Re
    │   └── Im
    ├── 2,0/
    │   ├── Re
    │   └── Im
    └── ...

To reload maps from disk, use load_h_maps_from_file():

from litebird_sim.mapmaking.h_maps import load_h_maps_from_file

result = load_h_maps_from_file("./h_n_maps/h_maps_det_mydetector.h5")

Detector and time splits#

The detector_split and time_split parameters follow the same conventions as the rest of the litebird_sim mapmaking framework (see Map-making). Use "full" (default) to include all detectors and all time samples.

MPI support#

The current implementation of make_h_maps() only allows to distribute observation by detector, i.e. each MPI process computes the h maps for a subset of detectors, but using all time samples. Example:

sim.create_observations(
   detectors=dets,
   n_blocks_det=2, # The detectors can be split over several MPI tasks
   n_blocks_time=1,  # But the time samples cannot be split, all are used for each detector
   tod_dtype=precision,
   split_list_over_processes=False,
   allocate_tod=False,
)
sim.prepare_pointings()

result = lbs.make_h_maps(
      observations=sim.observations,
      nside=64,
      n_m_couples=np.array([[0, 0], [2, 0], [4, 0]]),
      output_coordinate_system=lbs.CoordinateSystem.Galactic,
      save_to_file=True,
      output_directory="./h_n_maps",
   )

API reference#

class litebird_sim.mapmaking.h_maps.h_map_Re_and_Im(real: ndarray[tuple[Any, ...], dtype[_ScalarT]], imag: ndarray[tuple[Any, ...], dtype[_ScalarT]], n: int, m: int, det_info: str)#

Single \(h_{n,m}\) map component for one detector.

Variables:

  • real – Real part of \(h_{n,m}\).

  • imag – Imaginary part of \(h_{n,m}\).

  • n – Spin index \(n\).

  • m – Spin index \(m\).

  • det_info – Detector name associated with this map.

property norm#

Return \(|h_{n,m}|\) per pixel.

Returns:

Pixel-wise magnitude computed from real and imag.

Return type:

npt.NDArray

class litebird_sim.mapmaking.h_maps.HMapsResult(h_maps: dict[str, dict[tuple[int, int], h_map_Re_and_Im]], duration_s: float, sampling_rate_Hz: float, coordinate_system: CoordinateSystem = CoordinateSystem.Ecliptic, detector_split: str = 'full', time_split: str = 'full')#

Result of make_h_maps().

Variables:

  • h_maps – Mapping detector -> {(n, m): h_map_Re_and_Im} containing all computed \(h_{n,m}\) maps.

  • duration_s – Duration of the simulated observation in seconds.

  • sampling_rate_Hz – Sampling rate in hertz.

  • coordinate_system – Coordinate system of the output maps.

  • detector_split – Detector split used to build the maps.

  • time_split – Time split used to build the maps.

litebird_sim.mapmaking.h_maps.make_h_maps(observations: ~litebird_sim.observations.Observation | list[~litebird_sim.observations.Observation], nside: int, pointings: ~numpy.ndarray | list[~numpy.ndarray] | None = None, n_m_couples: ~numpy.ndarray = array([[0, 0],        [2, 0],        [4, 0]]), hwp: ~litebird_sim.hwp.HWP | None = None, output_coordinate_system: ~litebird_sim.coordinates.CoordinateSystem = CoordinateSystem.Galactic, detector_split: str = 'full', time_split: str = 'full', pointings_dtype=<class 'numpy.float64'>, save_to_file: bool = True, output_directory: str = './h_n_maps') → HMapsResult#

Generate complex harmonic maps \(h_{n,m}\) from observations.

The map for (n, m) = (0, 0) contains hit counts per pixel.

Parameters:

  • observations (Observation | list[Observation]) – Observation object or list of observations used to build the maps.

  • nside (int) – HEALPix nside of the output maps.

  • pointings (np.ndarray | list[np.ndarray] | None) – Optional external pointings. Expected shape is (n_detectors, n_samples, 3) per observation.

  • n_m_couples (np.ndarray) – Array of (n, m) pairs with shape (N, 2).

  • hwp (HWP | None) – Half-wave plate model.

  • output_coordinate_system (CoordinateSystem) – Coordinate system of the output maps.

  • detector_split (str) – Detector split to apply.

  • time_split (str) – Time split to apply.

  • pointings_dtype (type) – Data type used when pointings are computed on the fly.

  • save_to_file (bool) – If True, save maps to HDF5 files.

  • output_directory (str) – Output directory for generated HDF5 files.

Returns:

Result container with all computed maps and metadata.

Return type:

HnMapResult

litebird_sim.mapmaking.h_maps.save_h_maps(result, output_directory: str) → None#

Save \(h_n\) maps to HDF5 files.

One file is produced per detector in output_directory. The directory is created if it does not exist.

Parameters:

  • result (HnMapResult) – Container with maps and metadata to save.

  • output_directory (str) – Destination directory for output HDF5 files.

Returns:

None

Return type:

None

litebird_sim.mapmaking.h_maps.load_h_maps_from_file(filepath: str) → HMapsResult#

Load \(h_{n,m}\) maps from an HDF5 file.

Parameters:

filepath (str) – Path to an HDF5 file produced by save_hn_maps().

Returns:

Loaded maps and metadata.

Return type:

HnMapResult