diff --git a/CHANGELOG.md b/CHANGELOG.md index 3c3e09bc..8269a559 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,7 @@ # HEAD +- **Breaking change**: Reworking of the IO, `write_observations` and `read_observations` are now part of the class simulation [#293](https://github.com/litebird/litebird_sim/pull/293) + - Support for numpy.float128 made optional, this fixes importing issue on ARM architectures [#286](https://github.com/litebird/litebird_sim/pull/286) - Improve the documentation about noise simulations [#283](https://github.com/litebird/litebird_sim/pull/283) diff --git a/docs/source/mapmaking.rst b/docs/source/mapmaking.rst index f41f55ca..f832f382 100644 --- a/docs/source/mapmaking.rst +++ b/docs/source/mapmaking.rst @@ -1,3 +1,5 @@ +.. _mapmaking: + Map-making ========== @@ -117,7 +119,7 @@ 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 :func:`.make_binned_map` has a high level interface in the class +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 The syntax is identical to :func:`.make_binned_map`:: @@ -335,6 +337,14 @@ instance of the class :class:`.DestriperParameters`:: The result is an instance of the class :class:`.DestriperResult`, which is similar to :class:`.BinnerResult` but it contains much more information. +The function :func:`.make_destriped_map` has a high level interface in the class +:class:`.Simulation` that applies the destriper algorithm to all the +observations in the simulation. +The syntax is identical to :func:`.make_destriped_map`:: + + result = sim.make_destriped_map(nside=nside) + healpy.mollview(result.destriped_map) + We will now explain how a destriper works and what is the meaning of each parameter in the classes :class:`.DestriperParameters` and :class:`.DestriperResult`. Apart from KS2009, another source of information @@ -580,7 +590,6 @@ You can save the results of the destriper using the function using MPI, you should call both functions on *all* the MPI processes, and the number of processes should be the same between the two calls. - How the N_obs matrix is stored ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/docs/source/observations.rst b/docs/source/observations.rst index 4832333b..c7004a4f 100644 --- a/docs/source/observations.rst +++ b/docs/source/observations.rst @@ -1,3 +1,5 @@ +.. _observations: + Observations ============ @@ -225,12 +227,12 @@ following information are saved and restored: - Global and local flags saved in ``.global_flags`` and ``.local_flags`` (see below). -The function used to save observations is :func:`.write_observations`, -which works with a :class:`.Simulation` object; if you prefer to +The function used to save observations is :func:`.Simulation.write_observations`, +which acts on a :class:`.Simulation` object; if you prefer to operate without a :class:`.Simulation` object, you can call :func:`.write_list_of_observations`. -To read observations, you can use :func:`.read_observations` and +To read observations, you can use :func:`.Simulation.read_observations` and :func:`.read_list_of_observations`. diff --git a/docs/source/simulations.rst b/docs/source/simulations.rst index fde86f63..a57bb4a1 100644 --- a/docs/source/simulations.rst +++ b/docs/source/simulations.rst @@ -487,6 +487,38 @@ contains the size of the TOD array, the names of the detectors, and other useful information. Refer to the documentation of each class to know what is inside. +High level interface +-------------------- + +The class :class:`.Simulation` has a powerfull high level +interface that allows to quickly generate a scanning strategy +allocate the observations, generate simulated timelines +cointaing signal, noise and dipole, build maps, and +save(read) the entire simulation object. The syntax is +straightforward:: + + sim = lbs.Simulation(...) + + sim.set_scanning_strategy(...) + sim.set_instrument(...) + sim.create_observations(...) + sim.compute_pointings() + + sim.compute_pos_and_vel() + sim.add_dipole() + + sim.fill_tods(...) + sim.add_noise(...) + + result = sim.make_destriped_map(nside=nside) + healpy.mollview(result.destriped_map) + + sim.write_observations(...) + sim.read_observations(...) + +See the documentation in :ref:`observations`, :ref:`scanning-strategy` +:ref:`dipole-anisotropy`, :ref:`timeordered`, :ref:`mapmaking` for +details of the single functions. API reference ------------- diff --git a/litebird_sim/io.py b/litebird_sim/io.py index 98bc3309..2b9bb701 100644 --- a/litebird_sim/io.py +++ b/litebird_sim/io.py @@ -12,11 +12,13 @@ import h5py import numpy as np +from deprecation import deprecated + +from .version import __version__ as litebird_sim_version from .compress import rle_compress, rle_decompress from .detectors import DetectorInfo from .mpi import MPI_ENABLED, MPI_COMM_WORLD from .observations import Observation, TodDescription -from .simulations import Simulation __NUMPY_INT_TYPES = [ np.int8, @@ -361,59 +363,26 @@ def write_list_of_observations( return file_list +@deprecated( + deprecated_in="0.11", + current_version=litebird_sim_version, + details="Use Simulation.write_observations", +) def write_observations( - sim: Simulation, + sim, subdir_name: Union[None, str] = "tod", include_in_report: bool = True, *args, **kwargs, ) -> List[Path]: - """Write a set of observations as HDF5 - - This function is a wrapper to :func:`.write_list_of_observations` that saves - the observations associated with the simulation to a subdirectory within the - output path of the simulation. The subdirectory is named `subdir_name`; if - you want to avoid creating a subdirectory, just pass an empty string or None. - - This function only writes HDF5 for the observations that belong to the current - MPI process. If you have distributed the observations over several processes, - you must call this function on each MPI process. - - For a full explanation of the available parameters, see the documentation for - :func:`.write_list_of_observations`. - """ - - if subdir_name: - tod_path = sim.base_path / subdir_name - # Ensure that the subdirectory exists - tod_path.mkdir(exist_ok=True) - else: - tod_path = sim.base_path - - file_list = write_list_of_observations( - obs=sim.observations, path=tod_path, *args, **kwargs + # Here we call the method moved inside Simulation + return sim.write_observations( + subdir_name, + include_in_report, + *args, + **kwargs, ) - if include_in_report: - sim.append_to_report( - """ -## TOD files - -{% if file_list %} -The following files containing Time-Ordered Data (TOD) have been written: - -{% for file in file_list %} -- {{ file }} -{% endfor %} -{% else %} -No TOD files have been written to disk. -{% endif %} -""", - file_list=file_list, - ) - - return file_list - def __find_flags(inpf, expected_num_of_dets: int, expected_num_of_samples: int): flags_matches = [__FLAGS_GROUP_NAME_REGEXP.fullmatch(x) for x in inpf] @@ -661,28 +630,17 @@ def read_list_of_observations( return observations +@deprecated( + deprecated_in="0.11", + current_version=litebird_sim_version, + details="Use Simulation.read_observations", +) def read_observations( - sim: Simulation, + sim, path: Union[str, Path] = None, subdir_name: Union[None, str] = "tod", *args, **kwargs, ): - """Read a list of observations from a set of files in a simulation - - This function is a wrapper around the function :func:`.read_list_of_observations`. - It reads all the HDF5 files that are present in the directory whose name is - `subdir_name` and is a child of `path`, and it stores them in the - :class:`.Simulation` object `sim`. - - If `path` is not specified, the default is to use the value of ``sim.base_path``; - this is meaningful if you are trying to read back HDF5 files that have been saved - earlier in the same session. - """ - if path is None: - path = sim.base_path - - obs = read_list_of_observations( - file_name_list=list((path / subdir_name).glob("**/*.h5")), *args, **kwargs - ) - sim.observations = obs + # Here we call the method moved inside Simulation + sim.read_observations(path, subdir_name, *args, **kwargs) diff --git a/litebird_sim/simulations.py b/litebird_sim/simulations.py index 454851a9..c00821a3 100644 --- a/litebird_sim/simulations.py +++ b/litebird_sim/simulations.py @@ -23,6 +23,8 @@ from .healpix import write_healpix_map_to_file, npix_to_nside from .imo.imo import Imo from .mapmaking import ( + make_binned_map, + BinnerResult, make_destriped_map, DestriperParameters, DestriperResult, @@ -39,7 +41,7 @@ from .scan_map import scan_map_in_observations from .spacecraft import SpacecraftOrbit, spacecraft_pos_and_vel from .noise import add_noise_to_observations -from litebird_sim.mapmaking.binner import make_binned_map +from .io import write_list_of_observations, read_list_of_observations from .gaindrifts import GainDriftType, GainDriftParams, apply_gaindrift_to_observations import astropy.time @@ -1426,8 +1428,9 @@ def make_binned_map( self, nside: int, output_coordinate_system: CoordinateSystem = CoordinateSystem.Galactic, + components: Optional[List[str]] = None, append_to_report: bool = True, - ): + ) -> BinnerResult: """ Bins the tods of `sim.observations` into maps. The syntax mimics the one of :meth:`litebird_sim.make_binned_map` @@ -1447,8 +1450,152 @@ def make_binned_map( obs=self.observations, nside=nside, output_coordinate_system=output_coordinate_system, + components=components, + ) + + def make_destriped_map( + self, + nside: int, + params: DestriperParameters = DestriperParameters(), + components: Optional[List[str]] = None, + keep_weights: bool = False, + keep_pixel_idx: bool = False, + keep_pol_angle_rad: bool = False, + callback: Any = destriper_log_callback, + callback_kwargs: Optional[Dict[Any, Any]] = None, + append_to_report: bool = True, + ) -> DestriperResult: + """ + Bins the tods of `sim.observations` into maps. + The syntax mimics the one of :meth:`litebird_sim.make_binned_map` + """ + + results = make_destriped_map( + nside=nside, + obs=self.observations, + pointings=None, + params=params, + components=components, + keep_weights=keep_weights, + keep_pixel_idx=keep_pixel_idx, + keep_pol_angle_rad=keep_pol_angle_rad, + callback=callback, + callback_kwargs=callback_kwargs, + ) + + if append_to_report: + fig, ax = plt.subplots() + ax.set_xlabel("Iteration number") + ax.set_ylabel("Residual [K]") + ax.set_title("CG convergence of the destriper") + ax.semilogy( + np.arange(len(results.history_of_stopping_factors)), + results.history_of_stopping_factors, + "ko-", + ) + + template_file_path = get_template_file_path("report_destriper.md") + with template_file_path.open("rt") as inpf: + markdown_template = "".join(inpf.readlines()) + + cg_plot_filename = f"destriper-cg-convergence-{uuid4()}.png" + + self.append_to_report( + markdown_text=markdown_template, + results=results, + history_of_stopping_factors=[ + float(x) for x in results.history_of_stopping_factors + ], + bytes_in_cholesky_matrices=results.nobs_matrix_cholesky.nbytes, + cg_plot_filename=cg_plot_filename, + figures=[ + # Using uuid4() we can have more than one section + # about “destriping” in the report + (fig, cg_plot_filename), + ], + ) + + return results + + def write_observations( + self, + subdir_name: Union[None, str] = "tod", + append_to_report: bool = True, + *args, + **kwargs, + ) -> List[Path]: + """Write a set of observations as HDF5 + + This function is a wrapper to :func:`.write_list_of_observations` that saves + the observations associated with the simulation to a subdirectory within the + output path of the simulation. The subdirectory is named `subdir_name`; if + you want to avoid creating a subdirectory, just pass an empty string or None. + + This function only writes HDF5 for the observations that belong to the current + MPI process. If you have distributed the observations over several processes, + you must call this function on each MPI process. + + For a full explanation of the available parameters, see the documentation for + :func:`.write_list_of_observations`. + """ + + if subdir_name: + tod_path = self.base_path / subdir_name + # Ensure that the subdirectory exists + tod_path.mkdir(exist_ok=True) + else: + tod_path = self.base_path + + file_list = write_list_of_observations( + obs=self.observations, path=tod_path, *args, **kwargs ) + if append_to_report: + self.append_to_report( + """ + ## TOD files + + {% if file_list %} + The following files containing Time-Ordered Data (TOD) have been written: + + {% for file in file_list %} + - {{ file }} + {% endfor %} + {% else %} + No TOD files have been written to disk. + {% endif %} + """, + file_list=file_list, + ) + + return file_list + + def read_observations( + self, + path: Union[str, Path] = None, + subdir_name: Union[None, str] = "tod", + *args, + **kwargs, + ): + """Read a list of observations from a set of files in a simulation + + This function is a wrapper around the function :func:`.read_list_of_observations`. + It reads all the HDF5 files that are present in the directory whose name is + `subdir_name` and is a child of `path`, and it stores them in the + :class:`.Simulation` object `sim`. + + If `path` is not specified, the default is to use the value of ``sim.base_path``; + this is meaningful if you are trying to read back HDF5 files that have been saved + earlier in the same session. + """ + if path is None: + path = self.base_path + + obs = read_list_of_observations( + file_name_list=list((path / subdir_name).glob("**/*.h5")), *args, **kwargs + ) + self.observations = obs + def apply_gaindrift( self, drift_params: GainDriftParams = None, @@ -1528,62 +1675,3 @@ def apply_gaindrift( user_seed=user_seed, **dictionary, ) - - def make_destriped_map( - self, - nside: int, - params: DestriperParameters = DestriperParameters(), - components: Optional[List[str]] = None, - keep_weights: bool = False, - keep_pixel_idx: bool = False, - keep_pol_angle_rad: bool = False, - callback: Any = destriper_log_callback, - callback_kwargs: Optional[Dict[Any, Any]] = None, - append_to_report: bool = True, - ) -> DestriperResult: - results = make_destriped_map( - nside=nside, - obs=self.observations, - pointings=None, - params=params, - components=components, - keep_weights=keep_weights, - keep_pixel_idx=keep_pixel_idx, - keep_pol_angle_rad=keep_pol_angle_rad, - callback=callback, - callback_kwargs=callback_kwargs, - ) - - if append_to_report: - fig, ax = plt.subplots() - ax.set_xlabel("Iteration number") - ax.set_ylabel("Residual [K]") - ax.set_title("CG convergence of the destriper") - ax.semilogy( - np.arange(len(results.history_of_stopping_factors)), - results.history_of_stopping_factors, - "ko-", - ) - - template_file_path = get_template_file_path("report_destriper.md") - with template_file_path.open("rt") as inpf: - markdown_template = "".join(inpf.readlines()) - - cg_plot_filename = f"destriper-cg-convergence-{uuid4()}.png" - - self.append_to_report( - markdown_text=markdown_template, - results=results, - history_of_stopping_factors=[ - float(x) for x in results.history_of_stopping_factors - ], - bytes_in_cholesky_matrices=results.nobs_matrix_cholesky.nbytes, - cg_plot_filename=cg_plot_filename, - figures=[ - # Using uuid4() we can have more than one section - # about “destriping” in the report - (fig, cg_plot_filename), - ], - ) - - return results diff --git a/test/test_io.py b/test/test_io.py index 30b2866e..7f8e6d43 100644 --- a/test/test_io.py +++ b/test/test_io.py @@ -115,8 +115,7 @@ def __write_complex_observation( return ( obs, det, - lbs.write_observations( - sim=sim, + sim.write_observations( subdir_name="", gzip_compression=gzip_compression, tod_fields=["tod1", "tod2"], diff --git a/test/test_mpi.py b/test/test_mpi.py index 4ea02975..4f6d5deb 100644 --- a/test/test_mpi.py +++ b/test/test_mpi.py @@ -414,8 +414,8 @@ def test_write_hdf5_mpi(tmp_path): num_of_obs = 12 sim.create_observations(detectors=[det], num_of_obs_per_detector=num_of_obs) - file_names = lbs.write_observations( - sim, subdir_name="tod", file_name_mask="litebird_tod{global_index:04d}.h5" + file_names = sim.write_observations( + subdir_name="tod", file_name_mask="litebird_tod{global_index:04d}.h5" ) assert len(file_names) == len(sim.observations)