Implementation of h_n maps

CHANGELOG.md

@@ -1,5 +1,11 @@
 # HEAD
 
+-   Ducc0 1/f noise generation module included [#490](https://github.com/litebird/litebird_sim/pull/490)
+
+-   Small optimizations in beam convolver [#492](https://github.com/litebird/litebird_sim/pull/492)
+
+-   Fixed a TypeError in Observation when allocate_tod=False in MPI jobs [#491](https://github.com/litebird/litebird_sim/pull/491)
+
 -   Save memory in pointing generation [#488](https://github.com/litebird/litebird_sim/pull/488)
 
 -   **Breaking change**: Major reworking of the interfaces and handling of inputs across the framework [#479](https://github.com/litebird/litebird_sim/pull/479), in detail:

docs/source/beam_convolution.rst

@@ -80,6 +80,7 @@ as demonstrated in the following example:
         values=hp.synalm(np.ones((4,lmax+1)),lmax=lmax, mmax=lmax,),
         lmax=lmax,
         mmax=lmax,
+        coordinates=lbs.CoordinateSystem.Galactic,
     )
 
     Convparams = lbs.BeamConvolutionParameters(
@@ -156,6 +157,19 @@ dataclass. Allowed parameters:
 
 If convolution parameters are omitted, defaults are inferred from sky and beam alms, triggering a warning.
 
+.. note::
+   **Type Coherence and Memory Efficiency**
+
+   To optimize performance and minimize memory footprint, the convolution engine (based on ``ducc0``) 
+   requires specific floating-point precision. If the data types of the input ``sky_alms`` and 
+   ``beam_alms`` (instances of :class:`SphericalHarmonics`) do not match the precision 
+   specified in ``convolution_params`` (e.g., ``complex64`` for ``single_precision=True``), 
+   the code will perform a type conversion using ``astype(..., copy=False)``. 
+
+   This means that if the types are already compatible, no copy is made, and the data is 
+   processed in-place. However, if a conversion is necessary, a new array might be 
+   created or the existing one casted, but the framework is designed to avoid 
+   unnecessary deep copies whenever possible.
 
 Container for Spherical Harmonics
 ---------------------------------
@@ -180,7 +194,8 @@ It supports:
 - Resizing via zero-padding or truncation
 - I/O through Healpy-compatible FITS files
 
-
+Full documentation can be found in :ref: `maps_and_harmonics`.
+
 Example usage:
 
 .. testcode::
@@ -389,9 +404,6 @@ API reference
     :undoc-members:
     :show-inheritance:
 
-.. automodule:: litebird_sim.maps_and_harmonics
-    :members:
-
 .. automodule:: litebird_sim.beam_synthesis
     :members:
 

docs/source/index.rst

@@ -19,7 +19,6 @@ This is the User’s Manual of the LiteBIRD Simulation Framework.
    part5.rst
    part6.rst
    appendix.rst
-   installation
 
 Indices and tables
 ==================

docs/source/map_scanning.rst

@@ -210,7 +210,7 @@ transparent:
   sim.prepare_pointings()
 
   values_tuple = (np.ones(12*nside*nside)*1e-4, np.ones(12*nside*nside)*1e-4, np.ones(12*nside*nside)*1e-4)
-  sky_signal = lbs.HealpixMap(values=values_tuple, nside=nside)
+  sky_signal = lbs.HealpixMap(values=values_tuple, nside=nside, coordinates=lbs.CoordinateSystem.Galactic)
 
   sim.fill_tods(sky_signal)
 

docs/source/noise.rst

@@ -8,7 +8,10 @@ function :func:`.add_noise_to_observations` and the low-level versions
 :func:`.add_noise`, :func:`.add_white_noise`, and
 :func:`.add_one_over_f_noise`.
 
-Here is a short example that shows how to add noise:
+White noise
+-----------
+
+Here is a short example that shows how to add white noise to timelines:
 
 .. testcode::
 
@@ -90,7 +93,42 @@ call the low level function directly:
    custom_sigma_uk = 1234
    lbs.noise.add_white_noise(obs[0].tod[0], custom_sigma_uk, random=sim.dets_random[0])
 
-We can also add 1/f noise using a very similar call to the above:
+
+1/f Models and Engines
+----------------------
+
+The framework supports also the generatio of 1/f noise. Here you can choose the computational **engine** 
+(how it is calculated) and the physical **model** (the shape of the power spectrum).
+
+Engines
+^^^^^^^
+
+The engine is selected via the ``engine`` parameter:
+
+1.  **"fft" (Default)**: Generates noise in the Fourier domain.
+2.  **"ducc"**: Uses time-domain infinite impulse response (IIR) filtering provided by the `ducc0` library. It only supports the "keshner" model.
+
+Models
+^^^^^^
+
+The physical shape of the Power Spectral Density (PSD) is selected via the ``model`` parameter:
+
+1.  **"toast" (Default)**:
+    The classic power-law ratio, also implemented in https://github.com/hpc4cmb/toast/blob/372fa7642bbe61a5f01d239e707c04b80ad4bf46/src/toast/tod/sim_noise.py#L74. The PSD is proportional to:
+
+    .. math::
+
+        P(f) \propto \frac{f^\alpha + f_{knee}^\alpha}{f^\alpha + f_{min}^\alpha}
+
+2.  **"keshner"**:
+    Corresponds to a sum of relaxation processes. This is the native model of the `ducc` engine. The PSD is proportional to:
+
+    .. math::
+
+        P(f) \propto \left( \frac{f^2 + f_{knee}^2}{f^2 + f_{min}^2} \right)^{\alpha/2}
+
+
+This call allows to add 1/f noise:
 
 .. testcode::
 

litebird_sim/__init__.py

@@ -9,6 +9,9 @@
     DestriperParameters,
     DestriperResult,
     ExternalDestriperParameters,
+    HnMapResult,
+    load_h_map_from_file,
+    make_h_maps,
 )
 from .bandpasses import BandPassInfo
 from .beam_convolution import (

litebird_sim/beam_convolution.py

@@ -94,6 +94,11 @@ def add_convolved_sky_to_one_detector(
     -----
     - If no HWP is present, a standard 4π convolution is performed.
     - If HWP is present, the Mueller matrix is used to properly handle polarization.
+    - **Memory Management**: To ensure maximum efficiency, the sky and beam coefficients
+      are passed to the underlying engine using `copy=False` during type casting.
+      If the input coefficients' precision (e.g., complex128) differs from the
+      requested `convolution_params.single_precision`, a conversion will occur.
+      Users should ensure types are coherent to avoid even temporary memory overhead.
     - The function modifies `tod_det` in place by adding the convolved signal.
     """
 
@@ -130,6 +135,13 @@ def add_convolved_sky_to_one_detector(
         _slm = sky_alms_det.resize_alm(
             convolution_params.lmax, convolution_params.lmax, inplace=False
         ).values
+
+        logging.warning(
+            ("Input sky alms resized, using ℓ_max={lmax}.").format(
+                lmax=convolution_params.lmax
+            )
+        )
+
     else:
         _slm = sky_alms_det.values
 
@@ -139,6 +151,13 @@ def add_convolved_sky_to_one_detector(
         _blm = beam_alms_det.resize_alm(
             convolution_params.lmax, convolution_params.mmax, inplace=False
         ).values
+
+        logging.warning(
+            ("Input beam alms resized, using ℓ_max={lmax} and m_max={mmax}.").format(
+                lmax=convolution_params.lmax, mmax=convolution_params.mmax
+            )
+        )
+
     else:
         _blm = beam_alms_det.values
 
@@ -154,8 +173,8 @@ def add_convolved_sky_to_one_detector(
             intertype = Interpolator
 
         inter = intertype(
-            sky=_slm.astype(complex_type),
-            beam=_blm.astype(complex_type),
+            sky=_slm.astype(complex_type, copy=False),
+            beam=_blm.astype(complex_type, copy=False),
             separate=False,
             lmax=convolution_params.lmax,
             kmax=convolution_params.mmax,

litebird_sim/beam_synthesis.py

@@ -286,6 +286,7 @@ def gauss_bl(lmax: int, fwhm_rad: float, pol: bool = True) -> np.ndarray:
           - Index 0: Temperature beam b_l^T
           - Index 1: E-mode polarization beam b_l^E
           - Index 2: B-mode polarization beam b_l^B
+
         If False, returns a 1D array of shape (lmax+1,) (Temperature only).
 
     Returns

litebird_sim/hwp_diff_emiss.py

@@ -68,13 +68,13 @@ def add_2f_to_observations(
     By default, the TOD is added to ``Observation.tod``. If you want to add it to some
     other field of the :class:`.Observation` class, use `component`::
 
-    for cur_obs in sim.observations:
-        # Allocate a new TOD for the 2f alone
-        cur_obs.2f_tod = np.zeros_like(cur_obs.tod)
+        for cur_obs in sim.observations:
+            # Allocate a new TOD for the 2f alone
+            cur_obs.hwp_2f_tod = np.zeros_like(cur_obs.tod)
 
-        # Ask `add_2f_to_observations` to store the 2f
-        # in `observations.2f_tod`
-        add_2f_to_observations(sim.observations, component="2f_tod")
+            # Ask `add_2f_to_observations` to store the 2f
+            # in `observations.hwp_2f_tod`
+            add_2f_to_observations(sim.observations, component="hwp_2f_tod")
     """
     if isinstance(observations, Observation):
         obs_list = [observations]

litebird_sim/hwp_harmonics.py

@@ -486,7 +486,7 @@ def fill_tod(
             ``allocate_tod=False`` in :class:`.Observation`.
 
         maps : HealpixMap | SphericalHarmonics | dict[str, HealpixMap] |
-               dict[str, SphericalHarmonics]
+            dict[str, SphericalHarmonics]
             Sky model to be scanned. In dictionary form, the keys must match the
             entries of `input_names`.
 

litebird_sim/mapmaking/__init__.py

@@ -1,5 +1,6 @@
 from .common import ExternalDestriperParameters
 from .binner import make_binned_map, check_valid_splits, BinnerResult
+from .h_maps import HnMapResult, make_h_maps,load_h_map_from_file
 from .brahmap_gls import make_brahmap_gls_map
 from .destriper import (
     make_destriped_map,
@@ -19,6 +20,10 @@
     "BinnerResult",
     "make_binned_map",
     "check_valid_splits",
+    # h_n.py
+    "HnMapResult",
+    "make_h_maps",
+    "load_h_map_from_files",
     # brahmap_gls
     "make_brahmap_gls_map",
     # destriper.py

litebird_sim/mapmaking/binner.py

@@ -11,10 +11,10 @@
 from dataclasses import dataclass
 from typing import Any
 
-import healpy as hp
 import numpy as np
 import numpy.typing as npt
 from ducc0.healpix import Healpix_Base
+from litebird_sim.healpix import UNSEEN_PIXEL_VALUE
 from numba import njit
 
 from litebird_sim import mpi
@@ -86,8 +86,8 @@ def _solve_binning(nobs_matrix, atd):
             atd[ipix] = np.linalg.solve(nobs_matrix[ipix], atd[ipix])
             nobs_matrix[ipix] = np.linalg.inv(nobs_matrix[ipix])
         else:
-            nobs_matrix[ipix].fill(hp.UNSEEN)
-            atd[ipix].fill(hp.UNSEEN)
+            nobs_matrix[ipix].fill(UNSEEN_PIXEL_VALUE)
+            atd[ipix].fill(UNSEEN_PIXEL_VALUE)
 
 
 @njit

litebird_sim/mapmaking/common.py

@@ -131,6 +131,7 @@ def _compute_pixel_indices(
     hwp_angle: npt.NDArray | None,
     output_coordinate_system: CoordinateSystem,
     pointings_dtype=np.float64,
+    hmap_generation: bool = False,
 ) -> tuple[npt.NDArray, npt.NDArray]:
     """Compute the index of each pixel and its attack angle
 
@@ -144,12 +145,15 @@ def _compute_pixel_indices(
     ``N_d`` the number of detectors and ``N_t`` the number of samples in the TOD,
     and the last rank represents the θ and φ angles (in radians) expressed in the
     Ecliptic reference frame.
+
+    The option `hmap_generation` returns only the telescope orientation instead of
+    the polarization angle.
     """
 
     assert len(pol_angle_detectors) == num_of_detectors
 
-    pixidx_all = np.empty((num_of_detectors, num_of_samples), dtype=int)
-    polang_all = np.empty((num_of_detectors, num_of_samples), dtype=np.float64)
+    pixidx_all = np.empty((num_of_detectors, num_of_samples), dtype=np.int32)
+    polang_all = np.empty((num_of_detectors, num_of_samples), dtype=pointings_dtype)
 
     for idet in range(num_of_detectors):
         curr_pointings_det, hwp_angle = _get_pointings_array(
@@ -159,12 +163,14 @@ def _compute_pixel_indices(
             output_coordinate_system=output_coordinate_system,
             pointings_dtype=pointings_dtype,
         )
-
-        polang_all[idet] = _get_pol_angle(
-            curr_pointings_det=curr_pointings_det,
-            hwp_angle=hwp_angle,
-            pol_angle_detectors=pol_angle_detectors[idet],
-        )
+        if hmap_generation:
+            polang_all[idet] = curr_pointings_det[:, 2]
+        else:
+            polang_all[idet] = _get_pol_angle(
+                curr_pointings_det=curr_pointings_det,
+                hwp_angle=hwp_angle,
+                pol_angle_detectors=pol_angle_detectors[idet],
+            )
 
         pixidx_all[idet] = hpx.ang2pix(curr_pointings_det[:, :2])
 
@@ -177,6 +183,49 @@ def _compute_pixel_indices(
 
     return pixidx_all, polang_all
 
+def _compute_pixel_indices_single_detector(hpx: Healpix_Base,
+            pointings: npt.NDArray | Callable,
+            pol_angle_detector: float,
+            num_of_samples: int,
+            detector_index: int,
+            hwp_angle: npt.NDArray | None,
+            output_coordinate_system: CoordinateSystem,
+            pointings_dtype=np.float64,
+            hmap_generation: bool = False,
+            )-> tuple[npt.NDArray, npt.NDArray]:
+    """
+   Same as _compute_pixel_indices but for a single detector, thus returning only the pixel indices and polarization angles for that detector.
+    """
+    pixidx = np.empty((num_of_samples), dtype=np.int32)
+    polang = np.empty((num_of_samples), dtype=pointings_dtype)
+
+    curr_pointings_det, hwp_angle = _get_pointings_array(
+            detector_idx=detector_index,
+            pointings=pointings,
+            hwp_angle=hwp_angle,
+            output_coordinate_system=output_coordinate_system,
+            pointings_dtype=pointings_dtype,
+        )
+
+    if hmap_generation:
+        polang = curr_pointings_det[:, 2]
+    else:
+        polang = _get_pol_angle(
+            curr_pointings_det=curr_pointings_det,
+            hwp_angle=hwp_angle,
+            pol_angle_detectors=pol_angle_detector,
+        )
+    pixidx= hpx.ang2pix(curr_pointings_det[:, :2])
+
+    if output_coordinate_system == CoordinateSystem.Galactic:
+        # Free curr_pointings_det if the output map is already in Galactic coordinates
+        try:
+            del curr_pointings_det
+        except UnboundLocalError:
+            pass
+
+    return pixidx, polang
+
 
 def _cholesky_plain(A: npt.NDArray, dest_L: npt.NDArray) -> None:
     "Store a lower-triangular matrix in L such that A = L·L†"

litebird_sim/mapmaking/destriper.py

@@ -6,10 +6,11 @@
 from typing import Any, cast
 from collections.abc import Callable
 
-import healpy as hp
 import numpy as np
 import numpy.typing as npt
 from ducc0.healpix import Healpix_Base
+from litebird_sim.healpix import UNSEEN_PIXEL_VALUE
+
 from numba import njit, prange
 
 from litebird_sim.coordinates import CoordinateSystem, coord_sys_to_healpix_string
@@ -20,6 +21,8 @@
     _get_hwp_angle,
     _normalize_observations_and_pointings,
 )
+
+from litebird_sim.maps_and_harmonics import HealpixMap
 from .common import (
     _compute_pixel_indices,
     COND_THRESHOLD,
@@ -2142,7 +2145,7 @@ def _load_rank0_destriper_results(file_path: Path) -> DestriperResult:
             valid_pixel=np.array(inpf["MASK"].data.field("VALID"), dtype=bool),
             is_cholesky=bool(inpf["NOBSMATR"].header["ISCHOL"]),
         )
-        nside = hp.npix2nside(len(nobs_matrix.valid_pixel))
+        nside = HealpixMap.npix_to_nside(len(nobs_matrix.valid_pixel))
 
         if "GAL" in str(inpf["HITMAP"].header["COORDSYS"]).upper():
             coord_sys = CoordinateSystem.Galactic