.. _random-numbers: Random numbers ============== .. contents:: Table of Contents :depth: 2 :local: The LiteBIRD Simulation Framework can be used to generate random numbers, used for example for producing noise timelines. In order to do so, a seed and a Hierarchy of Random Number Generators (RNGs) are used. Seed ---- The ``random_seed`` is used to control the behavior of the RNG. The seed can be ``None`` or an integer number. If you are **not** interested in the reproducibility of your results, you can set ``random_seed`` to ``None``. However, this is not recommended, as if you run many times a function or method that uses an RNG, e.g., for producing noise timelines, then the outputs will be **different**, and you will not be able to re-obtain these results again. If, instead, you are interested in the reproducibility of your results, you can set ``random_seed`` to an integer number. With this choice, if you run a function or method multiple times using an RNG, then the outputs will be **the same**, and you will obtain these results again by re-running your script, as the random numbers produced by the generator will be the same. How should you set the ``random_seed``? This parameter **must** be passed to the constructor of a :class:`.Simulation` class. If you create a :class:`.Simulation` instance without passing the seed, the constructor will raise an error, and your script will fail. We implemented this feature to avoid automatic settings of the seed and unclear behavior of the generation of random numbers. If you run a parallel script using MPI, you do not have to worry about setting a different seed for different MPI ranks. The random number generator ensures that the random numbers generated by the ranks will produce orthogonal sequences as the first layer of the RNG hierarchy below the `random_seed`. Furthermore, a second layer of per-detector RNGs is created, ensuring that each detector of each MPI task is independent. This structure is stored in a :class:`.RNGHierarchy` instance, which handles this under the hood. We want this behavior as we do not want repetitions of, e.g., the same noise TOD if we split their computation on different MPI ranks. For example, in this way, if you split the TOD matrix of an :class:`.Observation` class by the time, you will not encounter the same noise after the samples generated by a certain rank; if you split the TOD matrix of an :class:`.Observation` class by the detectors, each one will have a different noise timestream. Regarding the reproducibility of the results in a parallel code, there is an **important thing** to bear in mind. If you set the seed to an integer number but run your script with a different number of MPI ranks, the outputs will be **different**! Think about a noise time stream of 4 samples. If you use 2 MPI ranks, then the first two samples will be generated by one RNG (rank 0), while the last two samples will be generated by a different RNG (rank 1). If you then run the same script with the same seed but with 4 MPI ranks, each of the samples will be generated by a different RNG, and only the first sample will be the same for the two runs, as it is always the first one generated by rank 0’s RNG. Despite this, the closest thing you can do towards **complete** reproducibility among different runs is the possibility to save to disk the instance of :class:`.RNGHierarchy`. This can be done at any moment calling :meth:`.RNGHierarchy.save`, which saves a frozen image of the state of the random generators, that can be reloaded through :meth:`.RNGHierarchy.load`. The setting of the ``random_seed`` is as simple as this:: sim = lbs.Simulation( base_path="/storage/output", start_time=astropy.time.Time("2020-02-01T10:30:00"), duration_s=3600.0, name="My noise simulation", description="I want to generate some noise and be able to reproduce my results", random_seed=12345, # Here the seed for the random number generator is set ) The :class:`.Simulation` constructor runs the :meth:`.Simulation.init_rng_hierarchy` method, which builds the MPI-layer of the RNG hierarchy and seeds it with ``random_seed``. When observations are created with :meth:`.Simulation.create_observations`, the detector-layer of the RNG hierarchy is initialized with :meth:`.Simulation.init_detectors_random`. You can then access these RNGs by calling ``sim.dets_random``, which is a list of RNGs. There is also the possibility to change the random seed after creating a :class:`.Simulation` instance. You can achieve this by calling :meth:`.Simulation.init_random`, which will regenerate the full hierarchy:: sim = lbs.Simulation(random_seed=12345) [...] sim.init_random(random_seed=42, num_detectors=3) Random number generators ------------------------ The random number generator used by the :class:`.RNGHierarchy` is `PCG64 `_. After creating this RNG hierarchy from the :class:`.Simulation` constructor, for example, you can access the detector-level geenrators via `sim.dets_random`:: sim = lbs.Simulation(random_seed=12345) [...] sim.add_noise(noise_type='white', dets_random=sim.dets_random) You can also use your own RNGs with the functions and methods of ``litebird_sim``:: sim = lbs.Simulation(random_seed=12345) [...] my_rngs = ... # New list of RNG definition (length must match number of detectors) sim.add_noise(noise_type='white', dets_random=my_rngs) You should just make sure that your custom RNG implements the ``normal`` method, so it can be used for white noise generation. Details on the RNG hierarchy ============================ As mentioned above, the RNG hierarchy system in ``litebird_sim`` provides a reproducible and flexible structure for managing random number generators across multiple levels of a simulation. This is especially useful in MPI-parallel contexts where independent and deterministic streams of randomness are required at different granularity levels—such as per-task or per-detector. The key classes and methods provided by this module are: 1. :class:`.RNGHierarchy`: A container managing a nested RNG tree from a single base seed. 2. :meth:`.RNGHierarchy.build_mpi_layer`: Method to build the first level of the RNG hierarchy for MPI tasks. 3. :meth:`.RNGHierarchy.build_detector_layer`: Method to build the second level of the RNG hierarchy for detectors. 4. :meth:`.RNGHierarchy.get_generator`: Method to get a speficic RNG generator from the hierarchy. 5. :meth:`.RNGHierarchy.get_detector_level_generators_on_rank`: Method to extract the full list of generators for a specific MPI rank. 6. :meth:`.RNGHierarchy.regenerate_or_check_detector_generators`: Method to check validity of a set of detector-level generators and eventually regenerate them from a new base seed. The following example demonstrates how you can create a layered RNG hierarchy: .. code-block:: python import litebird_sim as lbs base_seed = 12345 num_mpi_tasks = 2 num_detectors = 3 # Create a hierarchy of RNGs with two levels: task and detector rng_hierarchy = lbs.RNGHierarchy(base_seed) rng_hierarchy.build_mpi_layer(num_mpi_tasks) rng_hierarchy.build_detector_layer(num_detectors) # or directly with rng = lbs.RNGHierarchy(base_seed, num_mpi_tasks, num_detectors) Once the :class:`.RNGHierarchy` is instanciated and RNGs are created, it is sufficient to use :meth:`.RNGHierarchy.get_detector_level_generators_on_rank` to obtain all the generators of a specific MPI rank. If needed, :meth:`.RNGHierarchy.get_generator` allows you to retrieve a specific generator from the hierarchy. .. code-block:: python rank = 0 generators_of_first_rank = rng_hierarchy.get_detector_level_generators_on_rank(rank) indeces = (1, 2) generetor_second_rank_third_det = rng_hierarchy.get_generator(indeces) Furthermore, :meth:`.RNGHierarchy.add_extra_layer` allows you to define a new level of the hierarchy from the current leaf layer. Finally, the :meth:`.RNGHierarchy.regenerate_or_check_detector_generators` allows you to quickly regenerate a set of new detector-level generators, or to check the validity of the existing one. This is useful if you passes a random seed specific for each systematics. This needs a list of :class:`.Observation` to work. Saving and loading ------------------ RNG hierarchies can be serialized and deserialized via: - :meth:`.RNGHierarchy.save` - :meth:`.RNGHierarchy.load` This is useful for ensuring reproducibility even across different simulation runs. .. code-block:: python rng_hierarchy.save("my_rng_state.pkl") restored_rng_hierarchy = lbs.RNGHierarchy.load("my_rng_state.pkl") assert rng == restored_rng API reference ------------- .. automodule:: litebird_sim.seeding :members: :undoc-members: :show-inheritance: