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

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, numpy.random.SeedSequence, numpy.random.Generator, 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
+            (``SeedSequence(random_seed).spawn(number_of_simulations)``). It
+            accepts an int, a ``SeedSequence`` or a ``Generator``. 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,46 @@ 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):
+        """Return a ``SeedSequence`` root from a flexible seed argument.
+
+        Accepts what the scientific-Python SPEC 7 seeding convention accepts
+        (an int, a ``SeedSequence``, a ``Generator`` or ``BitGenerator``, or
+        None for fresh entropy) and returns a ``SeedSequence`` so it can be
+        spawned into one independent child stream per simulation index.
+        """
+        if isinstance(random_seed, np.random.SeedSequence):
+            return random_seed
+        if isinstance(random_seed, np.random.Generator):
+            return random_seed.bit_generator.seed_seq
+        if isinstance(random_seed, np.random.BitGenerator):
+            return random_seed.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, SeedSequence, Generator, optional
+            Root seed for the run. See ``simulate``.
+
         Returns
         -------
         None
@@ -280,14 +327,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 +360,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 +369,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, SeedSequence, Generator, optional
+            Root seed for the run. See ``simulate``.
 
         Returns
         -------
@@ -339,13 +392,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 +446,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 +464,11 @@ 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 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)

tests/unit/simulation/test_monte_carlo_determinism.py

@@ -0,0 +1,193 @@
+"""Determinism tests for the ``random_seed`` argument of ``MonteCarlo.simulate``.
+
+With a fixed ``random_seed`` the generated random *inputs* are reproducible and
+identical across serial and parallel execution and across any number of workers.
+Each simulation index draws from its own child stream spawned from the run's root
+seed, and ``SeedSequence.spawn`` is prefix-stable, so index ``i`` maps to the same
+seed regardless of the worker that runs it.
+
+The trajectory integration (``Flight``) is stubbed so the tests stay fast: worker
+invariance is a property of the input sampling, which happens before ``Flight`` is
+built. Stubbing the module-level ``Flight`` symbol reaches the parallel workers
+only under the ``fork`` start method, so those tests are guarded accordingly.
+
+A dedicated numpy-only rocket is used so *all* randomness flows through the seeded
+numpy generator. List-valued stochastic attributes are sampled with the standard
+library ``random.choice`` (an unseeded global generator) which ``random_seed``
+does not govern; the fixture drops the only such attribute (a multi-element
+``thrust_source``) so the inputs are byte-for-byte reproducible from the seed.
+"""
+
+import json
+
+import multiprocess
+import pytest
+
+import rocketpy.simulation.monte_carlo as mc_module
+from rocketpy.simulation import MonteCarlo
+from rocketpy.stochastic import StochasticRocket, StochasticSolidMotor
+
+pytestmark = pytest.mark.slow
+
+requires_fork = pytest.mark.skipif(
+    multiprocess.get_start_method() != "fork",
+    reason="stub-based parallel determinism test requires the 'fork' start method",
+)
+
+
+class _StubFlight:
+    """Minimal stand-in for ``Flight`` that skips trajectory integration."""
+
+    def __init__(self, **kwargs):  # accepts and ignores MonteCarlo's Flight kwargs
+        pass
+
+    def __getattr__(self, name):
+        return 0.0
+
+
+@pytest.fixture
+def stochastic_calisto_numpy_only(
+    cesaroni_m1670,
+    calisto_robust,
+    stochastic_nose_cone,
+    stochastic_trapezoidal_fins,
+    stochastic_tail,
+    stochastic_rail_buttons,
+    stochastic_main_parachute,
+    stochastic_drogue_parachute,
+):
+    """A ``StochasticRocket`` whose randomness flows entirely through numpy.
+
+    Mirrors the shared ``stochastic_calisto`` fixture but gives the solid motor a
+    single ``thrust_source`` instead of a multi-element list, so no attribute is
+    sampled through the unseeded standard-library ``random.choice``.
+    """
+    motor = StochasticSolidMotor(
+        solid_motor=cesaroni_m1670,
+        burn_out_time=(4, 0.1),
+        grains_center_of_mass_position=0.001,
+        grain_density=50,
+        grain_separation=1 / 1000,
+        grain_initial_height=1 / 1000,
+        grain_initial_inner_radius=0.375 / 1000,
+        grain_outer_radius=0.375 / 1000,
+        total_impulse=(6500, 1000),
+        throat_radius=0.5 / 1000,
+        nozzle_radius=0.5 / 1000,
+        nozzle_position=0.001,
+    )
+    rocket = StochasticRocket(
+        rocket=calisto_robust,
+        radius=0.0127 / 2000,
+        mass=(15.426, 0.5, "normal"),
+        inertia_11=(6.321, 0),
+        inertia_22=0.01,
+        inertia_33=0.01,
+        center_of_mass_without_motor=0,
+    )
+    rocket.add_motor(motor, position=0.001)
+    rocket.add_nose(stochastic_nose_cone, position=(1.134, 0.001))
+    rocket.add_trapezoidal_fins(stochastic_trapezoidal_fins, position=(0.001, "normal"))
+    rocket.add_tail(stochastic_tail)
+    rocket.set_rail_buttons(
+        stochastic_rail_buttons, lower_button_position=(-0.618, 0.001, "normal")
+    )
+    rocket.add_parachute(parachute=stochastic_main_parachute)
+    rocket.add_parachute(parachute=stochastic_drogue_parachute)
+    return rocket
+
+
+def _read_inputs_by_index(input_file):
+    """Read a ``.inputs.txt`` file into ``{index: raw_json_line}``."""
+    by_index = {}
+    with open(input_file, mode="r", encoding="utf-8") as rows:
+        for line in rows:
+            line = line.strip()
+            if not line:
+                continue
+            by_index[json.loads(line)["index"]] = line
+    return by_index
+
+
+def _simulate_inputs(
+    monkeypatch, tmp_path, environment, rocket, flight, tag, **simulate_kwargs
+):
+    """Run a Monte Carlo with a stubbed ``Flight`` and return inputs by index."""
+    monkeypatch.setattr(mc_module, "Flight", _StubFlight)
+    montecarlo = MonteCarlo(
+        filename=str(tmp_path / tag),
+        environment=environment,
+        rocket=rocket,
+        flight=flight,
+    )
+    montecarlo.simulate(**simulate_kwargs)
+    return _read_inputs_by_index(montecarlo.input_file)
+
+
+def test_serial_inputs_are_reproducible(
+    monkeypatch,
+    tmp_path,
+    stochastic_environment,
+    stochastic_calisto_numpy_only,
+    stochastic_flight,
+):
+    """Two serial runs with the same random_seed yield identical inputs."""
+    models = (stochastic_environment, stochastic_calisto_numpy_only, stochastic_flight)
+    run_a = _simulate_inputs(
+        monkeypatch, tmp_path, *models, "a", number_of_simulations=6, random_seed=7
+    )
+    run_b = _simulate_inputs(
+        monkeypatch, tmp_path, *models, "b", number_of_simulations=6, random_seed=7
+    )
+    assert sorted(run_a) == list(range(6))
+    assert run_a == run_b
+
+
+@requires_fork
+def test_inputs_are_worker_invariant(
+    monkeypatch,
+    tmp_path,
+    stochastic_environment,
+    stochastic_calisto_numpy_only,
+    stochastic_flight,
+):
+    """serial == parallel(2) == parallel(4): inputs are bit-identical per index."""
+    models = (stochastic_environment, stochastic_calisto_numpy_only, stochastic_flight)
+    common = {"number_of_simulations": 8, "random_seed": 314159}
+
+    serial = _simulate_inputs(monkeypatch, tmp_path, *models, "serial", **common)
+    par2 = _simulate_inputs(
+        monkeypatch, tmp_path, *models, "par2", parallel=True, n_workers=2, **common
+    )
+    par4 = _simulate_inputs(
+        monkeypatch, tmp_path, *models, "par4", parallel=True, n_workers=4, **common
+    )
+
+    expected = list(range(8))
+    assert sorted(serial) == expected
+    assert sorted(par2) == expected
+    assert sorted(par4) == expected
+    for index in expected:
+        assert serial[index] == par2[index], f"serial vs parallel(2) differ at {index}"
+        assert serial[index] == par4[index], f"serial vs parallel(4) differ at {index}"
+
+
+def test_none_seed_still_runs(
+    monkeypatch,
+    tmp_path,
+    stochastic_environment,
+    stochastic_calisto_numpy_only,
+    stochastic_flight,
+):
+    """random_seed=None draws fresh entropy but still exports one record per index."""
+    inputs = _simulate_inputs(
+        monkeypatch,
+        tmp_path,
+        stochastic_environment,
+        stochastic_calisto_numpy_only,
+        stochastic_flight,
+        "none",
+        number_of_simulations=5,
+        random_seed=None,
+    )
+    assert sorted(inputs) == list(range(5))