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

CHANGELOG.md

@@ -33,6 +33,8 @@ Attention: The newest changes should be on top -->
 ### Added
 
 - ENH: BUG/MNT: pre-release v1.13.0 review fixes [#1074](https://github.com/RocketPy-Team/RocketPy/pull/1074)
+- ENH: reproducible Monte Carlo runs via a random_seed argument [#1054](https://github.com/RocketPy-Team/RocketPy/pull/1054)
+
 ### Changed
 
 ### Fixed

rocketpy/simulation/monte_carlo.py

@@ -30,6 +30,7 @@
 from rocketpy.prints.monte_carlo_prints import _MonteCarloPrints
 from rocketpy.simulation.flight import Flight
 from rocketpy.tools import (
+    _seed_sequence_to_int,
     generate_monte_carlo_ellipses,
     generate_monte_carlo_ellipses_coordinates,
     import_optional_dependency,
@@ -170,6 +171,8 @@ def simulate(
         append=False,
         parallel=False,
         n_workers=None,
+        *,
+        random_seed=None,
         **kwargs,
     ):
         """
@@ -189,6 +192,22 @@ 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, numpy integer, sequence of ints, or 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: each simulation index derives its own
+            decorrelated child stream from this root, so index ``i`` receives the
+            same inputs no matter which worker runs it. A supplied ``SeedSequence``
+            is copied from its full state rather than consumed, so repeated calls
+            with the same seed reproduce the same inputs. Each model is reseeded
+            with a 128-bit integer -- the seed type a custom sampler's
+            ``reset_seed`` accepts. A stateful ``numpy.random.Generator`` or
+            ``BitGenerator`` is rejected (it is an RNG to draw from, not a fixed
+            seed); pass ``rng.bit_generator.seed_seq`` to seed from one. Default is
+            None, which draws fresh entropy on each run -- the previous,
+            non-reproducible default. This seeding is informed by Scientific Python
+            SPEC 7 but keeps immutable seed-snapshot semantics rather than sharing a
+            ``Generator``.
         kwargs : dict
             Custom arguments for simulation export of the ``inputs`` file. Options
             are:
@@ -223,6 +242,13 @@ def simulate(
         self.number_of_simulations = number_of_simulations
         self._initial_sim_idx = self.num_of_loaded_sims if append else 0
 
+        # Capture the small, picklable root seed state once per run (every
+        # simulation index derives its child seed from it, see __child_seed).
+        # This validates random_seed *before* __setup_files truncates any
+        # existing output, so an invalid seed cannot destroy prior results on
+        # the way to raising.
+        self.__capture_root_state(random_seed)
+
         print("Starting Monte Carlo analysis")
 
         self.__setup_files(append)
@@ -267,10 +293,84 @@ 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 __capture_root_state(self, random_seed):
+        """Capture the small, picklable root seed state for this run.
+
+        Stored once so serial mode and every parallel worker derive the same
+        per-index child seeds from it (see ``__child_seed``), instead of
+        materializing and pickling the full ``spawn(number_of_simulations)``
+        list to each process.
+        """
+        root = self.__root_seed_sequence(random_seed)
+        self.__root_state = (
+            root.entropy,
+            root.spawn_key,
+            root.pool_size,
+            root.n_children_spawned,
+        )
+
+    def __child_seed(self, sim_idx):
+        """Return the seed sequence for a single simulation index.
+
+        This equals ``root.spawn(number_of_simulations)[sim_idx]`` but is O(1)
+        in time and memory: ``SeedSequence.spawn`` derives child ``i`` by
+        appending ``n_children_spawned + i`` to the parent ``spawn_key``, so
+        rebuilding that one child directly reproduces it bit-for-bit while
+        letting a worker reconstruct any index from the small root state alone.
+        """
+        entropy, spawn_key, pool_size, base = self.__root_state
+        return np.random.SeedSequence(
+            entropy=entropy,
+            spawn_key=(*spawn_key, base + sim_idx),
+            pool_size=pool_size,
+        )
+
+    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. Each sub-stream is handed over as a 128-bit ``int`` (see
+        ``_seed_sequence_to_int``) so custom samplers keep working.
+        """
+        env_seed, rocket_seed, flight_seed = child_seed.spawn(3)
+        self.environment._set_stochastic(_seed_sequence_to_int(env_seed))
+        self.rocket._set_stochastic(_seed_sequence_to_int(rocket_seed))
+        self.flight._set_stochastic(_seed_sequence_to_int(flight_seed))
+
+    def __run_in_serial(self):  # pylint: disable=too-many-statements
         """
         Runs the monte carlo simulation in serial mode.
 
+        The root seed state is captured by ``simulate`` before this runs, so each
+        simulation index derives its child seed from ``self.__root_state``.
+
         Returns
         -------
         None
@@ -282,12 +382,13 @@ def __run_in_serial(self):
         )
         try:
             while sim_monitor.keep_simulating():
-                sim_monitor.increment()
+                sim_idx = sim_monitor.increment() - 1
                 inputs_json, outputs_json = "", ""
 
+                self.__seed_simulation(self.__child_seed(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)
@@ -313,6 +414,10 @@ def __run_in_parallel(self, n_workers=None):
         """
         Runs the monte carlo simulation in parallel.
 
+        The root seed state is captured by ``simulate`` before this runs and
+        travels with the pickled instance, so every worker derives the same
+        per-index child seed from ``self.__root_state``.
+
         Parameters
         ----------
         n_workers: int, optional
@@ -339,13 +444,16 @@ def __run_in_parallel(self, n_workers=None):
             )
 
             processes = []
-            seeds = np.random.SeedSequence().spawn(n_workers)
-
-            for seed in seeds:
+            # Each worker derives one independent child seed per simulation
+            # index (not per worker) from the shared root state: the counter
+            # assigns indices and index i always seeds from __child_seed(i), so
+            # the sampled inputs do not depend on the number of workers. The
+            # root state is small and travels with the pickled instance, so no
+            # per-index seed list is materialized or sent to each process.
+            for _ in range(n_workers):
                 sim_producer = multiprocess.Process(
                     target=self.__sim_producer,
                     args=(
-                        seed,
                         sim_monitor,
                         mutex,
                         simulation_error_event,
@@ -387,13 +495,11 @@ 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, 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.
         sim_monitor : _SimMonitor
             The simulation monitor object to keep track of the simulations.
         mutex : multiprocess.Lock
@@ -402,15 +508,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(self.__child_seed(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 +1712,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, running more simulations
+    than were requested (and duplicating a simulation index).
+    """
+    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.

rocketpy/stochastic/stochastic_model.py

@@ -3,8 +3,6 @@
 Stochastic classes.
 """
 
-from random import choice
-
 import numpy as np
 
 from rocketpy.mathutils.function import Function
@@ -508,6 +506,21 @@ def _validate_airfoil(self, airfoil):
                         "the first item"
                     )
 
+    def _random_choice(self, values):
+        """Choose one value from a list using this model's seeded generator.
+
+        The index is drawn from the seeded generator, not the stdlib global
+        ``random.choice`` (an unseeded shared instance), so the choice is
+        governed by ``random_seed``. Indexing rather than ``numpy.random.choice``
+        keeps a heterogeneous list -- ``Function`` objects, paths, arrays --
+        returned as itself instead of coerced to a common dtype. An empty
+        ``values`` is returned unchanged.
+        """
+        if not values:
+            return values
+        index = int(self.__random_number_generator.integers(len(values)))
+        return values[index]
+
     def dict_generator(self):
         """
         Generate a dictionary with randomly generated input arguments.
@@ -532,7 +545,7 @@ def dict_generator(self):
                 dist_sampler = value[-1]
                 generated_dict[arg] = dist_sampler(value[0], value[1])
             elif isinstance(value, list):
-                generated_dict[arg] = choice(value) if value else value
+                generated_dict[arg] = self._random_choice(value)
             elif isinstance(value, CustomSampler):
                 try:
                     generated_dict[arg] = value.sample(n_samples=1)[0]

rocketpy/stochastic/stochastic_rocket.py

@@ -1,7 +1,8 @@
 """Defines the StochasticRocket class."""
 
 import warnings
-from random import choice
+
+import numpy as np
 
 from rocketpy.control import _Controller
 from rocketpy.mathutils.vector_matrix import Vector
@@ -21,6 +22,7 @@
 from rocketpy.rocket.rocket import Rocket
 from rocketpy.stochastic.stochastic_generic_motor import StochasticGenericMotor
 from rocketpy.stochastic.stochastic_motor_model import StochasticMotorModel
+from rocketpy.tools import _seed_sequence_to_int
 
 from .stochastic_aero_surfaces import (
     StochasticAirBrakes,
@@ -176,21 +178,29 @@ def _set_stochastic(self, seed=None):
         """Set the stochastic attributes for Components, positions and
         inputs.
 
+        Every nested component -- the rocket body, each aerodynamic surface,
+        motor, rail button and parachute -- is reseeded from its own child of a
+        ``SeedSequence`` root, so components that sample the same distribution do
+        not draw identical values (a main and a drogue parachute get independent
+        ``cd_s`` and ``lag`` samples, not the same one). Children are spawned in a
+        fixed order, so the result stays reproducible under ``random_seed``.
+
         Parameters
         ----------
         seed : int, optional
             Seed for the random number generator.
         """
-        super()._set_stochastic(seed)
+        root = np.random.SeedSequence(seed)
+        super()._set_stochastic(_seed_sequence_to_int(root.spawn(1)[0]))
         self.aerodynamic_surfaces = self.__reset_components(
-            self.aerodynamic_surfaces, seed
+            self.aerodynamic_surfaces, root
         )
-        self.motors = self.__reset_components(self.motors, seed)
-        self.rail_buttons = self.__reset_components(self.rail_buttons, seed)
+        self.motors = self.__reset_components(self.motors, root)
+        self.rail_buttons = self.__reset_components(self.rail_buttons, root)
         for parachute in self.parachutes:
-            parachute._set_stochastic(seed)
+            parachute._set_stochastic(_seed_sequence_to_int(root.spawn(1)[0]))
 
-    def __reset_components(self, components, seed):
+    def __reset_components(self, components, root):
         """Creates a new Components whose stochastic structures
         and their positions are reset.
 
@@ -199,8 +209,9 @@ def __reset_components(self, components, seed):
         components : Components
             The components which contains the stochastic structure that
             will be used to create the new components.
-        seed : int, optional
-            Seed for the random number generator.
+        root : numpy.random.SeedSequence
+            The run's seed root. Each component is reseeded from its own spawned
+            child, so components sampling the same distribution stay decorrelated.
 
         Returns
         -------
@@ -212,7 +223,7 @@ def __reset_components(self, components, seed):
         new_components = Components()
         for stochastic_obj, _ in components:
             stochastic_obj_position_info = self.__components_map[stochastic_obj]
-            stochastic_obj._set_stochastic(seed)
+            stochastic_obj._set_stochastic(_seed_sequence_to_int(root.spawn(1)[0]))
             new_components.add(
                 stochastic_obj,
                 self._validate_position(stochastic_obj, stochastic_obj_position_info),
@@ -628,7 +639,7 @@ def _randomize_position(self, position):
                 return position[-1](position[0].z, position[1])
             return position[-1](position[0], position[1])
         elif isinstance(position, list):
-            return choice(position) if position else position
+            return self._random_choice(position)
 
     # pylint: disable=stop-iteration-return
     def dict_generator(self):
@@ -638,8 +649,8 @@ def dict_generator(self):
         all attributes of the class and generating a random value for each
         attribute. The random values are generated according to the format of
         each attribute. Tuples are generated using the distribution function
-        specified in the tuple. Lists are generated using the random.choice
-        function.
+        specified in the tuple. Lists are sampled through the model's seeded
+        generator so the choice is governed by ``random_seed``.
 
         Parameters
         ----------