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

MonteCarlo seeded the stochastic models per worker in parallel mode (from a
fresh, unseeded SeedSequence) and once at construction in serial mode, so the
sampled inputs depended on the execution mode and the worker count, and parallel
runs were not reproducible run to run.

Add a keyword-only random_seed to simulate() (SPEC 7 style: accepts an int, a
SeedSequence, or a Generator; None keeps the previous fresh-entropy behavior).
Spawn one child seed per simulation index from that root and reseed the
stochastic models from child_seeds[i] before simulation i. SeedSequence.spawn is
prefix-stable, so index i maps to the same seed regardless of which worker runs
it, making the inputs identical across serial, parallel(2) and parallel(N). Each
index seed is split three ways so the environment, rocket and flight draw from
independent streams rather than sharing one.

The serial index field now counts from 0 to match the parallel path. Both changes
alter the numbers a fixed seed produces, so stored baselines regenerate.

Adds tests/unit/simulation/test_monte_carlo_determinism.py: serial
reproducibility, worker invariance (serial == parallel(2) == parallel(4)), and
the None-seed path.

Signed-off-by: thc1006 <84045975+thc1006@users.noreply.github.com>

Author: thc1006

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))