Implement time-dependent quaternions and save memory when computing p…

…ointings (#319) * Fix wrong terminology in the manual * Fix wrong docstring for _compute_pixel_indices * Add a few more comments and improve docstrings * Fix the formatting of a docstring * Improve the docstring for Simulation.set_scanning_strategy * Improve the clarity of a few tests * Improve the pointing generation benchmark * Ignore ASDF files * Add missing @njit * Remove a comment that's no longer relevant * Rename Spin2EclipticQuaternions into TimeDependentQuaternion * Rename get_detector_quats() into slerp() and add __mul__() * TimeDependentQuaternion → RotQuaternion, implement pointing calculation * Update benchmarks/*.py to the new API for get_pointings * Add `HWP.apply_hwp_to_pointings` * Implement the class `PointingProvider` * Update the benckmark * Implement PointingProvider.has_hwp() * Export PointingProvider * Remove deprecated function `read_observations` * Remove declaration for deprecated function `read_observations` * Add I/O capabilities to the HWP class * Remove "dtype_pointing" keyword from Simulation.compute_pointings * Load/save quaternions and HWP information in HDF5 files * Update the documentation about HDF5 files * Remove get_pointing_buffer_shape, get_pointings, get_pointings_for_observations * Test the ability to save the pointings * Make test_coordinates() work with the new API * Rename Simulation.compute_pointings() to Simulation.prepare_pointings() * Support multiple detectors in Observation.get_pointings() * Fix the order of imports * Create a low-level prepare_pointings function * Support multiple detectors and lbs.prepare_pointings() in Simulation class * Update the signature of `rotate_coordinates_e2g` * Remove unused private functions * Update the parameter names used with rotate_coordinates_e2g * Update formatting * Add 'pointings_dtype=' parameter to prepare_pointings() and get_pointings() * Extract `precompute_pointings()` from prepare_pointings() * Remove wrong reference to Observation.pointing_coords * Remove wrong reference to Observation.pointing_coords * Update the documentation * Improve docstrings for HWP * updated to new pointings hwp_sys.py and test_hwp_sys.py * reformatted hwp_sys.py * update example notebook hwp_sys, still missing the pointing update regarding lbs.make_binned_map() * minor adjustment to _compute_pixel_indices() * hwp_sys example notebook update (make_binned_map still missing computation of pointing on the fly) * choosing dtype for pointings computation in hwp_sys * reformatting * dipole.py updated * test fixed * binner fixed * fix test_scan_map.py * test_binner fixed * test_scanning fixed * fixed hwp_sys.py * test_coordinates fixed * madam interface fixed * updated test_hwp_sys.py * updated hwp_sys.py with coord in gal choice * update hwp_sys/examples/simple_scan.ipynb with pointings onthefly for make_binned_map * simulation.prepare_pointings cleaned * standardize variable names * fix syntax * Update `destriper` to new scheme * New naming convention (obs -> observations) * Fix tests destriper * Propagate new naming * Docs * Correct typing * Correct docs * Avoid printing polangle in docs * Print all (θ, φ, ψ) in docs * update hwp_sys.py * update hwp_angle in hwp_sys.py and relative test * some documentation added * get_pointings has "all" as default * add quaternion function in `quaternions.py` * add `pointing_sys.py` * fix a document error in `quaternions.py` * add `__init__.py` in `pointing_sys` * add functions of `pointing_sys` * add `test_quat_rotations` in `test` * add `test_pointing_sys.py` in `test` * fix to change github action macos image from arm64 to x86_64 architecture to allow installation of toast2 on macos * add some documents in `pointing_sys.py` * simplify: `get_detector_orientation()` * debug: `_ecl2focalplane()` * Changed all tests to compare with reference file. * add script `gen_mock_ficalplane.py` which creates a toml file includes `DetectorInfo`. * add `mock_focalplane.toml` which is used by `test_pointing_sys.py` * add reference file `pointing_sys_reference.json.gz` for test * refactoring `pointing_sys.py` * add a notebook `Tutorial_PointingSys.ipynb` * add comments to `Tutorial_PointingSys.ipynb` * update changelog and status * fix Tutorial_PointingSys.ipynb * Adapt data splits notebook to new API * Change `dtype_tod` into `tod_dtype` for consistency * update notebook * put `atol` for every tests * Test revert macos change * revert * fix pointing_matrix returned by get_pointings * test fixed * fix syntax * again syntax fix * changed the github action image for macos to macos-13; fixed syntax issues * updated poetry.lock to sync markdown-katex version * Fix poetry.lock * [skip ci] Update CHANGELOG --------- Co-authored-by: Nicolò Raffuzzi <[email protected]> Co-authored-by: sgiardie <[email protected]> Co-authored-by: Luca Pagano <[email protected]> Co-authored-by: Giacomo Galloni <[email protected]> Co-authored-by: yusuke-takase <[email protected]> Co-authored-by: Avinash Anand <[email protected]>
litebird · Jun 25, 2024 · d6c6386 · d6c6386
1 parent 058a12d
commit d6c6386
Show file tree

Hide file tree

Showing 68 changed files with 8,222 additions and 4,050 deletions.
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
@@ -8,7 +8,7 @@ jobs:
     timeout-minutes: 30
     strategy:
       matrix:
-        os: [ubuntu-latest, macos-latest]
+        os: [ubuntu-latest, macos-13]
         mpi: ["none", openmpi] # mpich
         python: [3.9, 3.11]
         exclude:

diff --git a/.gitignore b/.gitignore
@@ -1,3 +1,6 @@
+# This file is created by ASDF
+.tool-versions
+
 ### PyCharm ###
 # Covers JetBrains IDEs: IntelliJ, RubyMine, PhpStorm, AppCode, PyCharm, CLion, Android Studio, WebStorm and Rider
 # Reference: https://intellij-support.jetbrains.com/hc/en-us/articles/206544839

diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -1,5 +1,23 @@
 # HEAD
 
+-   **Breaking change**: new API for pointing computation [#319](https://github.com/litebird/litebird_sim/pull/319). Here is a in-depth list of all the breaking changes in this PR:
+
+    1.  Quaternions describing the orientation of the detectors must now be encoded using a `RotQuaternion` object; plain NumPy arrays are no longer supported.
+
+    2.  Quaternions are now computed using the function `prepare_pointings()` (low-level) and the method `Simulation.prepare_pointings()` (high-level, you should use this). Pointings are no longer kept in memory until you retrieve them using `Observation.get_pointings()`.
+
+    3.  Pointings are no longer accessible using the field `pointings` in the `Observation` class. (Not 100% true, see below.) They are computed on the fly by the method `Observation.get_pointings()`.
+
+    4.  The way pointings are returned differs from how they were stored before. The result of a call to `Observation.get_pointings()` is a 2-element tuple: the first element contains a `(N, 3)` NumPy array containing the colatitude θ, the longitude φ, and the orientation ψ, while the second element is an array of the angles of the HWP. Thus, the orientation angle ψ is now stored together with θ and φ.
+
+    5.  If you want to pre-compute all the pointings instead of computing them on the fly each time you call `Observation.get_pointings()`, you can use the function `precompute_pointings()` (low-level) and the method `Simulation.precompute_pointings()` (high-level). This initializes a number of fields in each `Observation` object, but they are shaped as described in the previous point, i.e., ψ is kept in the same matrix as θ and φ.
+
+    6.  The argument `dtype_tod` of the method `Simulation.create_observations` has become `tod_type` for consistency with other similar parameters.
+
+    7.  The format of the HDF5 files has been slightly changed to let additional information about pointings to be stored.
+
+    See the comments in [PR#319](https://github.com/litebird/litebird_sim/pull/319) and discussion [#312](https://github.com/litebird/litebird_sim/discussions/312) for more details.
+
 -   Add data splits in time and detector space to destriped maps [#309](https://github.com/litebird/litebird_sim/pull/309)
 
 -   Fix issue [#317](https://github.com/litebird/litebird_sim/issues/317)

diff --git a/STATUS.md b/STATUS.md
@@ -11,15 +11,18 @@ Simulation Team plans to implement in `litebird_sim`.
 | Interface with Ephemerides       | Complete |          | Through AstroPy                 | [#48](https://github.com/litebird/litebird_sim/pull/48)                                                              |
 | Synthetic sky map generation     | Complete |          | Based on PySM3                  | [#76](https://github.com/litebird/litebird_sim/pull/76)                                                              |
 | Binning map-maker                | Complete |          |                                 | [#73](https://github.com/litebird/litebird_sim/pull/76)                                                              |
-| Destriping+calibration map-maker | Complete |          | Provided by TOAST+Madam         | [#86](https://github.com/litebird/litebird_sim/pull/86)                                                              |
+| Destriping+calibration map-maker | Complete |          | Internal destriper +            | [#260](https://github.com/litebird/litebird_sim/pull/260)                                                            |
+|                                  |          |          | interface with TOAST and Madam  | [#86](https://github.com/litebird/litebird_sim/pull/86)                                                              |
 |                                  |          |          |                                 | [#186](https://github.com/litebird/litebird_sim/pull/186)                                                            |
+| Splits in map-makers             | Complete |          |                                 | [#291](https://github.com/litebird/litebird_sim/pull/291)                                                            |
 | Spacecraft simulator             | Complete |          |                                 | [#122](https://github.com/litebird/litebird_sim/pull/122)                                                            |
 | Dipole calculation               | Complete |          |                                 | [#122](https://github.com/litebird/litebird_sim/pull/122)                                                            |
 | Map scanning                     | Complete |          |                                 | [#131](https://github.com/litebird/litebird_sim/pull/131)                                                            |
 | White+1/f noise generation       | Complete |          |                                 | [#100](https://github.com/litebird/litebird_sim/pull/100)                                                            |
 | Gain drift simulation            | Complete |          |                                 | [#243](https://github.com/litebird/litebird_sim/pull/243)                                                            |
 | Synthetic bandpass generation    | Complete |          |                                 | [#160](https://github.com/litebird/litebird_sim/pull/160), [#200](https://github.com/litebird/litebird_sim/pull/200) |
 | Calibration non-idealities       | Complete |          |                                 | [#243](https://github.com/litebird/litebird_sim/pull/243)                                                            |
+| Pointing systematics             | Complete |          |                                 | [#319](https://github.com/litebird/litebird_sim/pull/319)                                                            |
 | Beam convolution                 | Partial  |          | Through ducc0                   | [ducc.totalconvolve](https://gitlab.mpcdf.mpg.de/mtr/ducc/-/tree/ducc0/)                                             |
 | Cosmic-ray glitch generation     | Partial  |          |                                 | No PRs yet                                                                                                           |
 | HWP simulation                   | Partial  |          |                                 | [#117](https://github.com/litebird/litebird_sim/pull/117)                                                            |
@@ -29,14 +32,13 @@ Simulation Team plans to implement in `litebird_sim`.
 
 ## Beam convolution
 
--   `ducc0` already provides a 4π convolution code, and it is already
-    available within `litebird_sim`
+-   `ducc0` already provides a 4π convolution code, and it is already available within `litebird_sim`
 -   A high-level interface to `ducc0` is still missing
 
 ## Destriping+calibration map-maker
 
+-   Internal destriper, [#260](https://github.com/litebird/litebird_sim/pull/260)
 -   Provided by TOAST2, [PR#86](https://github.com/litebird/litebird_sim/pull/86)
-
 -   [PR#186](https://github.com/litebird/litebird_sim/pull/186) adds the possibility to interface Madam
 
 ## Calibration non-idealities
@@ -55,17 +57,14 @@ Simulation Team plans to implement in `litebird_sim`.
 
 -   A mathematical model is already available, based on [Giardiello et al. (2021)](https://ui.adsabs.harvard.edu/abs/2021arXiv210608031G/abstract).
 -   [PR#117](https://github.com/litebird/litebird_sim/pull/117).
+-   This supports both Jones and Mueller formalisms, and it's interfaced with the bandpass generation module
 
 ## ADC simulation
 
 -   Need to simulate the following effects:
-
     -   Signal quantization
-
     -   Clipping of the signal outside the dynamic range
-
     -   Non-linearity effects
-
 -   Signal clipping is already available in the Cosmic-ray glitch
     generator (see above)
 

diff --git a/benchmarks/.gitignore b/benchmarks/.gitignore
@@ -0,0 +1,2 @@
+pointings.npy
+results/
diff --git a/benchmarks/pointing_generation.py b/benchmarks/pointing_generation.py
@@ -1,13 +1,15 @@
+#!/usr/bin/env python3
 # -*- encoding: utf-8 -*-
 
 """
 This program tests the speed of the code that generates quaternions
 and pointings.
 """
 
-import time
 import sys
+import time
 from pathlib import Path
+
 import numpy as np
 
 # Add the `..` directory to PYTHONPATH, so that we can import "litebird_sim"
@@ -28,7 +30,12 @@
 )
 
 sim.create_observations(
-    detectors=[lbs.DetectorInfo(name="dummy", sampling_rate_hz=50.0)],
+    detectors=[
+        lbs.DetectorInfo(name="dummy1", sampling_rate_hz=30.0),
+        lbs.DetectorInfo(name="dummy2", sampling_rate_hz=30.0),
+        lbs.DetectorInfo(name="dummy3", sampling_rate_hz=30.0),
+        lbs.DetectorInfo(name="dummy4", sampling_rate_hz=30.0),
+    ],
     num_of_obs_per_detector=1,
     split_list_over_processes=False,
 )
@@ -58,15 +65,22 @@
 )
 
 instr = lbs.InstrumentInfo(spin_boresight_angle_rad=np.deg2rad(15.0))
+sim.set_instrument(instr)
 
 # Compute the pointings by running a "slerp" operation
+sim.prepare_pointings()
+
 start = time.perf_counter_ns()
-pointings_and_orientation = lbs.get_pointings(
-    obs,
-    spin2ecliptic_quats=sim.spin2ecliptic_quats,
-    detector_quats=np.array([[0.0, 0.0, 0.0, 1.0]]),
-    bore2spin_quat=instr.bore2spin_quat,
+
+pointings_and_orientation = np.empty(
+    shape=(len(sim.detectors), sim.observations[0].n_samples, 3),
+    dtype=np.float64,
 )
+for cur_obs in sim.observations:
+    for det_idx in range(cur_obs.n_detectors):
+        (cur_pointings, hwp_angle) = cur_obs.get_pointings(det_idx)
+        pointings_and_orientation[det_idx, :, :] = cur_pointings
+
 stop = time.perf_counter_ns()
 elapsed_time = (stop - start) * 1.0e-9
 
@@ -77,3 +91,16 @@
         pointings_and_orientation.shape[1] / elapsed_time
     ),
 )
+
+array_file = Path("pointings.npy")
+
+if array_file.exists():
+    with array_file.open("rb") as inp_f:
+        reference = np.load(inp_f)
+        np.save(file=Path("difference.npy"), arr=reference - pointings_and_orientation)
+        np.testing.assert_array_almost_equal(reference, pointings_and_orientation)
+        print(f'The array looks the same as the one in "{array_file}"')
+else:
+    with array_file.open("wb") as out_f:
+        np.save(out_f, pointings_and_orientation)
+    print(f'Array saved for reference in "{array_file}"')
diff --git a/docs/source/dipole.rst b/docs/source/dipole.rst
@@ -235,7 +235,7 @@ simulation :func:`.Simulation.add_dipole`.
 
   sim.create_observations(detectors=det)
 
-  sim.compute_pointings()
+  sim.prepare_pointings()
 
   sim.compute_pos_and_vel()
 

diff --git a/docs/source/gaindrifts.rst b/docs/source/gaindrifts.rst
@@ -67,7 +67,7 @@ For any kind of gain drift, one can use either the method of :class:`.Simulation
 
     # Applying gain drift on the given TOD component of an `Observation` object
     lbs.apply_gaindrift_to_observations(
-        obs=sim1.observations,
+        observations=sim1.observations,
         drift_params=drift_params,
         component="gain_2_obs",
     )

diff --git a/docs/source/images/tutorial-maps.png b/docs/source/images/tutorial-maps.png
diff --git a/docs/source/index.rst b/docs/source/index.rst
@@ -13,20 +13,20 @@ Welcome to litebird_sim's documentation!
    installation
    tutorial
    simulations
+   imo
    detectors
    observations
    data_layout
-   mapmaking
-   sky_maps
    scanning
-   bandpasses
-   dipole
-   imo
+   sky_maps
    timeordered
-   reports
+   dipole
+   mapmaking
+   bandpasses
+   random_numbers
    mpi
+   reports
    gaindrifts
-   random_numbers
    external_modules
    profiling
    integrating

diff --git a/docs/source/mapmaking.rst b/docs/source/mapmaking.rst
@@ -75,7 +75,7 @@ detectors::
                 name="0B", sampling_rate_hz=10, quat=lbs.quat_rotation_z(np.pi / 2)
             ),
         ],
-        dtype_tod=np.float64,  # Needed if you use the TOAST destriper
+        tod_dtype=np.float64,  # Needed if you use the TOAST destriper
         n_blocks_time=lbs.MPI_COMM_WORLD.size,
         split_list_over_processes=False,
     )
@@ -111,20 +111,21 @@ The output map is in Galactic coordinates, but you can specify the
 coordinate system you want via the parameter `output_coordinates`.
 This is how it should be called::
 
-    result = lbs.make_binned_map(nside=128, obs=obs)
-    healpy.mollview(result.binned_map)
+    result = lbs.make_binned_map(nside=128, observations=observations)
+    healpy.mollview(result.binned_map[0])
 
-(The pointing information is included in the :class:`.Observation`,
-alternatively pointings can be provided as a list of numpy arrays.)
+The pointing information is obtained from the :class:`.Observation` (either 
+using the precomputed pointing or computed on the fly), alternatively 
+pointings can be provided as a list of numpy arrays.
 The result is an instance of the class :class:`.BinnerResult`
 and contains both the I/Q/U maps and the covariance matrix.
 
 The function :func:`.make_binned_map` has a high level interface in the class
-:class:`.Simulation` that bins the content of the observations into maps
+:class:`.Simulation` that bins the content of the observations into maps.
 The syntax is identical to :func:`.make_binned_map`::
 
     result = sim.make_binned_map(nside=nside)
-    healpy.mollview(result.binned_map)
+    healpy.mollview(result.binned_map[0])
 
 
 The :class:`.BinnerResult` contains the field ``binned_map``, which
@@ -367,9 +368,13 @@ instance of the class :class:`.DestriperParameters`::
         ...
     )
 
-    result = lbs.make_destriped_map(nside=nside, obs=obs, params=params)
+    result = lbs.make_destriped_map(nside=nside,
+                                    observations=observations,
+                                    params=params,
+                                    )
     healpy.mollview(result.destriped_map)
 
+The pointing information is handled identically to :func:`.make_binned_map`.
 The result is an instance of the class :class:`.DestriperResult`, which
 is similar to :class:`.BinnerResult` but it contains much more information.
 
@@ -563,7 +568,7 @@ following fields in the :class:`.DestriperParameters` class:
   then the CG algorithm stops.
 - ``samples_per_baseline``: this can either be an integer, in which case it will
   be used for *all* the baselines, or a list of 1D arrays, each containing the
-  length of each baseline for each observation passed through the parameter ``obs``.
+  length of each baseline for each observation passed through the parameter ``observations``.
   Note that if you provide an integer, it might be that not all baselines will
   have exactly that length: it depends whether the number :math:`N_t` of samples
   in the TOD is evenly divisibile by ``samples_per_baseline`` or not. The
@@ -587,7 +592,7 @@ went:
 
 The baselines are saved in the field ``baselines`` of the :class:`.DestriperResult`
 class; this is a list of 2D arrays, where each element in the list is
-associated with one of the observations passed in the parameter ``obs``. The
+associated with one of the observations passed in the parameter ``observations``. The
 shape of each 2D arrays is :math:`(N_d, N_b),` where
 :math:`N_d` is the number of detectors for the observation and :math:`N_b` is
 the number of baselines. A visual representation of the memory layout of
@@ -632,7 +637,7 @@ the destriper will skip the CG iterations and proceed directly to the
 map-making step.
 
 How the N_obs matrix is stored
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 The destriper uses a different method to store the matrix :math:`M` in
 memory. As the 3×3 sub-blocks of this matrix need to be inverted often
@@ -725,7 +730,7 @@ would be produced by the *binner*, see above), and the *destriped map*.
    floating point numbers. As the default data type for timelines
    created by ``sim.create_observations`` is a 32-bit float, if you
    plan to run the destriper you should pass the flag
-   ``dtype_tod=np.float64`` to ``sim.create_observations`` (see the
+   ``tod_dtype=np.float64`` to ``sim.create_observations`` (see the
    code above), otherwise ``destripe`` will create an internal copy of
    the TOD converted in 64-bit floating-point numbers, which is
    usually a waste of space.
@@ -780,7 +785,7 @@ so you should call ``mpiexec``, ``mpirun``, or something similar.
 
 
 Creating several maps with Madam
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 There are cases where you want to create several maps out of
 one simulation. A common case is when you simulate several