Enable the use of different dtypes in observations (#79)

The method `Simulation.create_observation` now accepts a `dtype_tod` keyword that can be used to specify the numeric type to use for the TOD, like in the following example:

```python
obs_list = sim.create_observations(
    detectors=detector_list, num_of_obs_per_detector=1, dtype_tod=numpy.float128,
)
```

CHANGELOG.md

@@ -1,5 +1,7 @@
 # HEAD
 
+- Add the ability to specify the size of the floating-point type used in `Observation` objects [PR#79](https://github.com/litebird/litebird_sim/pull/79)
+
 - Simple bin map-maker [PR#73](https://github.com/litebird/litebird_sim/pull/73)
 
 - Use SI units in class `SpinningScanningStrategy` (**breaking change**) [PR#69](https://github.com/litebird/litebird_sim/pull/69)

docs/source/observations.rst

@@ -48,12 +48,31 @@ object::
   )
 
 Note that the 2-D array ``obs.tod`` is created for you. Its shape is
-``(n_detectors, n_samples)``. In full scale simulations it may get too large to
-fit in memory. You can chunk it along the time or detector dimension
-(or both) using ``n_blocks_det, n_blocks_time, comm`` at construction time or
-the `set_n_blocks` method. The same chunking is applied also to any detector
-information that you add with `detector_global_info` or `detector_info`. Note
-note ``det_idx`` is added for you at construction.
+``(n_detectors, n_samples)``, and the default type is
+``numpy.float32``: the choice of a 32-bit type is usually good enough
+for the purposes of LiteBIRD simulations, but if you need less/more
+precision, you are free to use any of the floating-point types
+provided by NumPy (``float16``, ``float32``, ``float64``,
+``float128``) using the keyword ``dtype_tod``::
+
+  import numpy as np
+  
+  obs = lbs.Observation(
+      detectors=2,
+      start_time=0.0,
+      sampling_rate_hz=5.0,
+      n_samples=5,
+      dtype_tod=np.float64,  # Use a 64-bit floating point type
+  )
+
+  
+In full scale simulations the TOD may get too large to fit in memory.
+You can chunk it along the time or detector dimension (or both) using
+``n_blocks_det, n_blocks_time, comm`` at construction time or the
+`set_n_blocks` method. The same chunking is applied also to any
+detector information that you add with `detector_global_info` or
+`detector_info`. Note note ``det_idx`` is added for you at
+construction.
 
 When you distribute the observation the first
 ``n_blocks_det x n_blocks_time``  MPI ranks are organized in a row-major grid

litebird_sim/simulations.py

@@ -25,6 +25,7 @@
 import astropy.time
 import astropy.units
 import markdown
+import numpy as np
 import jinja2
 import tomlkit
 
@@ -606,6 +607,7 @@ def create_observations(
         n_blocks_det=1,
         n_blocks_time=1,
         root=0,
+        dtype_tod=np.float32,
     ):
         "Create a set of Observation objects"
 
@@ -637,6 +639,7 @@ def create_observations(
                 n_blocks_time=n_blocks_time,
                 comm=(None if distribute else self.mpi_comm),
                 root=0,
+                dtype_tod=dtype_tod,
             )
             observations.append(cur_obs)
 

test/test_simulations.py

@@ -224,15 +224,21 @@ def test_duration_units_in_parameter_file():
 
 
 def test_distribute_observation(tmp_path):
-    sim = lbs.Simulation(
-        base_path=tmp_path / "simulation_dir", start_time=1.0, duration_s=11.0
-    )
-    det = lbs.DetectorInfo("dummy", sampling_rate_hz=15)
-    obs_list = sim.create_observations(detectors=[det], num_of_obs_per_detector=5)
+    for dtype in (np.float16, np.float32, np.float64, np.float128):
+        sim = lbs.Simulation(
+            base_path=tmp_path / "simulation_dir", start_time=1.0, duration_s=11.0
+        )
+        det = lbs.DetectorInfo("dummy", sampling_rate_hz=15)
+        obs_list = sim.create_observations(
+            detectors=[det], num_of_obs_per_detector=5, dtype_tod=dtype
+        )
 
-    assert len(obs_list) == 5
-    assert int(obs_list[-1].get_times()[-1] - obs_list[0].get_times()[0]) == 10
-    assert sum([o.n_samples for o in obs_list]) == sim.duration_s * det.sampling_rate_hz
+        assert len(obs_list) == 5
+        assert int(obs_list[-1].get_times()[-1] - obs_list[0].get_times()[0]) == 10
+        assert (
+            sum([o.n_samples for o in obs_list])
+            == sim.duration_s * det.sampling_rate_hz
+        )
 
 
 def test_distribute_observation_astropy(tmp_path):