ENH: reproducible Monte Carlo via per-simulation-index seeding

CHANGELOG.md

@@ -32,7 +32,9 @@ Attention: The newest changes should be on top -->
 
 ### Added
 
+- ENH: reproducible Monte Carlo runs via a random_seed argument [#1054](https://github.com/RocketPy-Team/RocketPy/pull/1054)
 - ENH: MNT: remove duplicate and merge-sync entries from Unreleased changelog [#1062](https://github.com/RocketPy-Team/RocketPy/pull/1062)
+
 ### Changed
 
 ### Fixed

rocketpy/simulation/monte_carlo.py

@@ -170,6 +170,8 @@ def simulate(
         append=False,
         parallel=False,
         n_workers=None,
+        *,
+        random_seed=None,
         **kwargs,
     ):  # pylint: disable=too-many-statements
         """
@@ -189,6 +191,15 @@ def simulate(
             number of workers will be equal to the number of CPUs available.
             A minimum of 2 workers is required for parallel mode.
             Default is None.
+        random_seed : int or numpy.random.SeedSequence, optional
+            Root seed for the run. When provided, the sampled inputs are
+            reproducible and identical across serial and parallel execution
+            and across any number of workers, because each simulation index
+            draws from its own child stream spawned from this root. It accepts
+            an int or a ``SeedSequence``; a supplied ``SeedSequence`` is copied
+            rather than consumed, so repeated calls with the same seed produce
+            the same inputs. Default is None, which draws fresh entropy (not
+            reproducible), preserving the previous behavior.
         kwargs : dict
             Custom arguments for simulation export of the ``inputs`` file. Options
             are:
@@ -228,9 +239,9 @@ def simulate(
         self.__setup_files(append)
 
         if parallel:
-            self.__run_in_parallel(n_workers)
+            self.__run_in_parallel(n_workers, random_seed)
         else:
-            self.__run_in_serial()
+            self.__run_in_serial(random_seed)
 
         self.__terminate_simulation()
 
@@ -267,10 +278,53 @@ def __setup_files(self, append):
         except OSError as error:
             raise OSError(f"Error creating files: {error}") from error
 
-    def __run_in_serial(self):
+    @staticmethod
+    def __root_seed_sequence(random_seed):
+        """Build a fresh ``SeedSequence`` root from ``random_seed``.
+
+        ``random_seed`` may be an int (or any entropy ``numpy.random.SeedSequence``
+        accepts), an existing ``SeedSequence``, or ``None`` for fresh entropy. A
+        supplied ``SeedSequence`` is copied from its full ``state``, so the
+        spawning below neither mutates the caller's object nor advances a shared
+        child counter between calls; repeated ``simulate`` calls with the same
+        seed then stay reproducible. A stateful ``Generator``/``BitGenerator`` is
+        not accepted, since using it as an immutable seed would contradict its
+        consume-on-use semantics; pass ``rng.bit_generator.seed_seq`` to seed
+        from an existing generator's stream.
+        """
+        if isinstance(random_seed, np.random.SeedSequence):
+            return np.random.SeedSequence(**random_seed.state)
+        if isinstance(random_seed, (np.random.Generator, np.random.BitGenerator)):
+            raise TypeError(
+                "random_seed must be an int or a numpy.random.SeedSequence, not "
+                f"a {type(random_seed).__name__}; to seed from an existing "
+                "generator pass rng.bit_generator.seed_seq."
+            )
+        return np.random.SeedSequence(random_seed)
+
+    def __seed_simulation(self, child_seed):
+        """Reseed the stochastic models for a single simulation index.
+
+        The per-index child seed is split three ways so the environment,
+        rocket and flight draw from independent streams instead of sharing
+        one. Seeding per simulation index (not per worker) is what makes the
+        sampled inputs invariant to the execution mode and to the number of
+        workers.
+        """
+        env_seed, rocket_seed, flight_seed = child_seed.spawn(3)
+        self.environment._set_stochastic(env_seed)
+        self.rocket._set_stochastic(rocket_seed)
+        self.flight._set_stochastic(flight_seed)
+
+    def __run_in_serial(self, random_seed=None):  # pylint: disable=too-many-statements
         """
         Runs the monte carlo simulation in serial mode.
 
+        Parameters
+        ----------
+        random_seed : int or SeedSequence, optional
+            Root seed for the run. See ``simulate``.
+
         Returns
         -------
         None
@@ -280,14 +334,18 @@ def __run_in_serial(self):
             n_simulations=self.number_of_simulations,
             start_time=time(),
         )
+        child_seeds = self.__root_seed_sequence(random_seed).spawn(
+            self.number_of_simulations
+        )
         try:
             while sim_monitor.keep_simulating():
-                sim_monitor.increment()
+                sim_idx = sim_monitor.increment() - 1
                 inputs_json, outputs_json = "", ""
 
+                self.__seed_simulation(child_seeds[sim_idx])
                 flight = self.__run_single_simulation()
-                inputs_json = self.__evaluate_flight_inputs(sim_monitor.count)
-                outputs_json = self.__evaluate_flight_outputs(flight, sim_monitor.count)
+                inputs_json = self.__evaluate_flight_inputs(sim_idx)
+                outputs_json = self.__evaluate_flight_outputs(flight, sim_idx)
 
                 with open(self.input_file, "a", encoding="utf-8") as f:
                     f.write(inputs_json)
@@ -309,7 +367,7 @@ def __run_in_serial(self):
                 f.write(inputs_json)
             raise error
 
-    def __run_in_parallel(self, n_workers=None):
+    def __run_in_parallel(self, n_workers=None, random_seed=None):
         """
         Runs the monte carlo simulation in parallel.
 
@@ -318,6 +376,8 @@ def __run_in_parallel(self, n_workers=None):
         n_workers: int, optional
             Number of workers to be used. If None, the number of workers
             will be equal to the number of CPUs available. Default is None.
+        random_seed : int or SeedSequence, optional
+            Root seed for the run. See ``simulate``.
 
         Returns
         -------
@@ -339,13 +399,19 @@ def __run_in_parallel(self, n_workers=None):
             )
 
             processes = []
-            seeds = np.random.SeedSequence().spawn(n_workers)
+            # One independent child seed per simulation index (not per
+            # worker), shared with every worker. The shared counter assigns
+            # indices, and index i always seeds from child_seeds[i], so the
+            # sampled inputs do not depend on the number of workers.
+            child_seeds = self.__root_seed_sequence(random_seed).spawn(
+                self.number_of_simulations
+            )
 
-            for seed in seeds:
+            for _ in range(n_workers):
                 sim_producer = multiprocess.Process(
                     target=self.__sim_producer,
                     args=(
-                        seed,
+                        child_seeds,
                         sim_monitor,
                         mutex,
                         simulation_error_event,
@@ -387,13 +453,16 @@ def __validate_number_of_workers(self, n_workers):
             raise ValueError("Number of workers must be at least 2 for parallel mode.")
         return n_workers
 
-    def __sim_producer(self, seed, sim_monitor, mutex, error_event):  # pylint: disable=too-many-statements
+    def __sim_producer(self, child_seeds, sim_monitor, mutex, error_event):  # pylint: disable=too-many-statements
         """Simulation producer to be used in parallel by multiprocessing.
 
         Parameters
         ----------
-        seed : int
-            The seed to set the random number generator.
+        child_seeds : list[numpy.random.SeedSequence]
+            One seed sequence per simulation index. Before each simulation
+            the worker seeds the stochastic models from
+            ``child_seeds[sim_idx]``, where ``sim_idx`` comes from the shared
+            counter, so the inputs are invariant to the number of workers.
         sim_monitor : _SimMonitor
             The simulation monitor object to keep track of the simulations.
         mutex : multiprocess.Lock
@@ -402,15 +471,14 @@ def __sim_producer(self, seed, sim_monitor, mutex, error_event):  # pylint: disa
             Event signaling an error occurred during the simulation.
         """
         try:
-            # Ensure Processes generate different random numbers
-            self.environment._set_stochastic(seed)
-            self.rocket._set_stochastic(seed)
-            self.flight._set_stochastic(seed)
+            while True:
+                sim_idx = _claim_next_index(sim_monitor, mutex)
+                if sim_idx is None:
+                    break
 
-            while sim_monitor.keep_simulating():
-                sim_idx = sim_monitor.increment() - 1
                 inputs_json, outputs_json = "", ""
 
+                self.__seed_simulation(child_seeds[sim_idx])
                 flight = self.__run_single_simulation()
                 inputs_json = self.__evaluate_flight_inputs(sim_idx)
                 outputs_json = self.__evaluate_flight_outputs(flight, sim_idx)
@@ -1607,6 +1675,24 @@ def export_errors_to_json(self, filename):
         self._write_log_to_json(self.errors_log, filename)
 
 
+def _claim_next_index(sim_monitor, mutex):
+    """Atomically claim the next 0-based simulation index, or ``None`` if done.
+
+    ``keep_simulating()`` and ``increment()`` are two separate manager calls, so
+    the shared ``mutex`` has to be held across both. Without it, two workers can
+    each pass the ``count < number_of_simulations`` check at the tail before
+    either increments, and both then claim an index, overrunning the requested
+    number of simulations and indexing past the per-index seed list.
+    """
+    mutex.acquire()
+    try:
+        if not sim_monitor.keep_simulating():
+            return None
+        return sim_monitor.increment() - 1
+    finally:
+        mutex.release()
+
+
 def _import_multiprocess():
     """Import the necessary modules and submodules for the
     multiprocess library.