From bf93b990790f069d9b415fdbff8e9cd11b4cce19 Mon Sep 17 00:00:00 2001 From: thc1006 <84045975+thc1006@users.noreply.github.com> Date: Tue, 7 Jul 2026 09:20:05 +0800 Subject: [PATCH 1/8] 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> --- rocketpy/simulation/monte_carlo.py | 94 +++++++-- .../test_monte_carlo_determinism.py | 193 ++++++++++++++++++ 2 files changed, 269 insertions(+), 18 deletions(-) create mode 100644 tests/unit/simulation/test_monte_carlo_determinism.py diff --git a/rocketpy/simulation/monte_carlo.py b/rocketpy/simulation/monte_carlo.py index 21c665d01..88461de52 100644 --- a/rocketpy/simulation/monte_carlo.py +++ b/rocketpy/simulation/monte_carlo.py @@ -170,6 +170,8 @@ def simulate( append=False, parallel=False, n_workers=None, + *, + random_seed=None, **kwargs, ): """ @@ -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) diff --git a/tests/unit/simulation/test_monte_carlo_determinism.py b/tests/unit/simulation/test_monte_carlo_determinism.py new file mode 100644 index 000000000..c3e0e47e9 --- /dev/null +++ b/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)) From 58cf1201179b25243280a818a271781297da37fc Mon Sep 17 00:00:00 2001 From: thc1006 <84045975+thc1006@users.noreply.github.com> Date: Thu, 9 Jul 2026 04:49:34 +0800 Subject: [PATCH 2/8] TST: cover Monte Carlo seeding helpers directly The reproducible-seeding change added __root_seed_sequence and __seed_simulation plus the per-index serial and parallel seeding, but the only tests that reached them ran a full Monte Carlo and were marked slow, so the coverage job (which does not pass --runslow) never executed them. Add fast unit tests that drive the two helpers directly: every supported random_seed type normalizes to the same root stream, None draws fresh entropy, existing SeedSequence/Generator/BitGenerator objects are reused rather than copied, and each child seed splits three ways so environment, rocket and flight get independent streams. Move the end-to-end simulate reproducibility tests into tests/integration, next to the existing Monte Carlo simulate test. The serial reproducibility run now lives in the non-slow suite; only the fork-based worker-invariance test stays slow, and it imports multiprocess lazily like the library does. Signed-off-by: thc1006 <84045975+thc1006@users.noreply.github.com> --- .../test_monte_carlo_determinism.py | 177 ++++++++++ .../test_monte_carlo_determinism.py | 307 ++++++++---------- 2 files changed, 304 insertions(+), 180 deletions(-) create mode 100644 tests/integration/simulation/test_monte_carlo_determinism.py diff --git a/tests/integration/simulation/test_monte_carlo_determinism.py b/tests/integration/simulation/test_monte_carlo_determinism.py new file mode 100644 index 000000000..b47629a55 --- /dev/null +++ b/tests/integration/simulation/test_monte_carlo_determinism.py @@ -0,0 +1,177 @@ +"""End-to-end determinism tests for ``MonteCarlo.simulate(random_seed=...)``. + +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 seed-handling helpers themselves +are unit tested in ``tests/unit/simulation/test_monte_carlo_determinism``.) + +The trajectory integration (``Flight``) is stubbed: worker invariance is a +property of the *input sampling*, which happens before ``Flight`` is built, so a +stub keeps the runs fast while still driving the real serial and parallel loops. +Stubbing the module-level ``Flight`` symbol reaches the parallel workers only +under the ``fork`` start method, so the worker-invariance test skips otherwise and +is marked ``slow`` to match the other Monte Carlo multiprocessing tests. + +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 pytest + +import rocketpy.simulation.monte_carlo as mc_module +from rocketpy.simulation import MonteCarlo +from rocketpy.stochastic import StochasticRocket, StochasticSolidMotor + + +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 seed yield byte-identical inputs per index. + + This drives the serial ``simulate`` path end to end; the flexible seed types + are covered by the unit test of ``__root_seed_sequence``. + """ + models = (stochastic_environment, stochastic_calisto_numpy_only, stochastic_flight) + run_a = _simulate_inputs( + monkeypatch, tmp_path, *models, "a", number_of_simulations=3, random_seed=7 + ) + run_b = _simulate_inputs( + monkeypatch, tmp_path, *models, "b", number_of_simulations=3, random_seed=7 + ) + assert sorted(run_a) == list(range(3)) + assert run_a == run_b + + +@pytest.mark.slow +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.""" + multiprocess = pytest.importorskip("multiprocess") + if multiprocess.get_start_method() != "fork": + pytest.skip( + "stub-based parallel determinism test requires the 'fork' start method" + ) + + 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}" diff --git a/tests/unit/simulation/test_monte_carlo_determinism.py b/tests/unit/simulation/test_monte_carlo_determinism.py index c3e0e47e9..939980508 100644 --- a/tests/unit/simulation/test_monte_carlo_determinism.py +++ b/tests/unit/simulation/test_monte_carlo_determinism.py @@ -1,193 +1,140 @@ -"""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. +"""Unit tests for the Monte Carlo seeding helpers. + +``MonteCarlo.simulate(random_seed=...)`` makes the sampled inputs reproducible by +turning the run's root seed into one independent child stream per simulation +index. Two private helpers do the work: + +* ``__root_seed_sequence`` normalizes the flexible ``random_seed`` argument (int, + ``SeedSequence``, ``Generator``, ``BitGenerator`` or None) into a + ``SeedSequence`` that can be spawned; +* ``__seed_simulation`` splits one per-index child seed three ways so the + environment, rocket and flight draw from independent streams. + +These tests exercise the helpers directly, with no fixtures and no simulation, so +they stay fast. The end-to-end reproducibility of ``simulate`` (serial and across +workers) is covered by ``tests/integration/simulation/test_monte_carlo_determinism``. + +Reaching a name-mangled member is an established pattern in this suite (see +``tests/unit/test_sensitivity.py`` and ``tests/unit/environment/test_environment.py``); +it lets the seeding invariants be asserted without running a Monte Carlo. """ -import json +from types import SimpleNamespace -import multiprocess +import numpy as np 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 +_root_seed_sequence = MonteCarlo._MonteCarlo__root_seed_sequence +_seed_simulation = MonteCarlo._MonteCarlo__seed_simulation + + +def _entropy(seed_sequence, n=4): + """A stable, comparable fingerprint of a ``SeedSequence``'s stream.""" + return tuple(int(x) for x in seed_sequence.generate_state(n)) + + +# --------------------------------------------------------------------------- # +# __root_seed_sequence: normalizing the flexible seed argument # +# --------------------------------------------------------------------------- # -requires_fork = pytest.mark.skipif( - multiprocess.get_start_method() != "fork", - reason="stub-based parallel determinism test requires the 'fork' start method", + +@pytest.mark.parametrize( + "make_seed", + [ + pytest.param(lambda: 12345, id="int"), + pytest.param(lambda: np.random.SeedSequence(12345), id="seedsequence"), + pytest.param(lambda: np.random.default_rng(12345), id="generator"), + pytest.param(lambda: np.random.PCG64(12345), id="bitgenerator"), + ], +) +def test_root_seed_sequence_accepts_supported_types(make_seed): + """int, SeedSequence, Generator and BitGenerator all normalize to the same + root SeedSequence stream for an equivalent seed value.""" + root = _root_seed_sequence(make_seed()) + assert isinstance(root, np.random.SeedSequence) + assert _entropy(root) == _entropy(_root_seed_sequence(12345)) + + +def test_root_seed_sequence_none_draws_fresh_entropy(): + """None yields a SeedSequence seeded from fresh OS entropy (not reproducible).""" + root = _root_seed_sequence(None) + assert isinstance(root, np.random.SeedSequence) + assert root.entropy is not None + + +@pytest.mark.parametrize( + "make_seed, resolve", + [ + pytest.param( + lambda: np.random.SeedSequence(999), + lambda seed: seed, + id="seedsequence", + ), + pytest.param( + lambda: np.random.default_rng(999), + lambda seed: seed.bit_generator.seed_seq, + id="generator", + ), + pytest.param( + lambda: np.random.PCG64(999), + lambda seed: seed.seed_seq, + id="bitgenerator", + ), + ], ) +def test_root_seed_sequence_reuses_existing_seed_sequence(make_seed, resolve): + """When given something that already carries a SeedSequence, the helper + reuses that object rather than copying it.""" + seed = make_seed() + assert _root_seed_sequence(seed) is resolve(seed) -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 - ) +# --------------------------------------------------------------------------- # +# __seed_simulation: splitting one child seed across the three models # +# --------------------------------------------------------------------------- # + + +class _RecordingModel: + """Stand-in stochastic model that records the seeds it is handed.""" + + def __init__(self): + self.seeds = [] + + def _set_stochastic(self, seed=None): + self.seeds.append(seed) - 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, + +def _split_seeds(child_seed): + """Run ``__seed_simulation`` against recording models; return the three seeds.""" + models = SimpleNamespace( + environment=_RecordingModel(), + rocket=_RecordingModel(), + flight=_RecordingModel(), ) - assert sorted(inputs) == list(range(5)) + _seed_simulation(models, child_seed) + return models.environment.seeds, models.rocket.seeds, models.flight.seeds + + +def test_seed_simulation_decorrelates_env_rocket_flight(): + """The per-index child seed is split three ways so environment, rocket and + flight draw from independent streams instead of sharing one.""" + env_seeds, rocket_seeds, flight_seeds = _split_seeds(np.random.SeedSequence(2024)) + assert [len(env_seeds), len(rocket_seeds), len(flight_seeds)] == [1, 1, 1] + fingerprints = { + _entropy(env_seeds[0]), + _entropy(rocket_seeds[0]), + _entropy(flight_seeds[0]), + } + assert len(fingerprints) == 3 + + +def test_seed_simulation_is_deterministic_per_child(): + """A given child seed reseeds the three models identically every time.""" + + def split(child): + env, rocket, flight = _split_seeds(child) + return [_entropy(env[0]), _entropy(rocket[0]), _entropy(flight[0])] + + assert split(np.random.SeedSequence(2024)) == split(np.random.SeedSequence(2024)) From 8eea1d85ae523e5259a3532fc9501c206864a5d6 Mon Sep 17 00:00:00 2001 From: thc1006 <84045975+thc1006@users.noreply.github.com> Date: Sat, 11 Jul 2026 09:14:37 +0800 Subject: [PATCH 3/8] BUG: fix Monte Carlo seeding race and non-reproducible SeedSequence Addresses review feedback on #1054. Parallel workers claimed the next index with an unlocked keep_simulating() + increment(), so near the end of a run two workers could both pass the count < n check and then claim sim_idx == n; the per-index child_seeds lookup turned that into an IndexError (before, it only wrote one extra record). Move the claim into a _claim_next_index helper that holds the shared mutex across the check and the increment, so each index is handed out once and the counter never overshoots. A deterministic unit test (a barrier plus a widened check-to-increment window) over-claims and fails if the lock is dropped. __root_seed_sequence returned the caller's SeedSequence, and spawn() advances its child counter, so passing the same object to simulate() twice produced different children. Copy it from its full state instead, which leaves the caller untouched and keeps repeated calls reproducible. Also drop Generator/BitGenerator from the accepted types: a stateful generator is not a seed, and reducing it to its underlying SeedSequence ignores how far it has been consumed. random_seed now takes an int, a sequence of ints, or a SeedSequence. Signed-off-by: thc1006 <84045975+thc1006@users.noreply.github.com> --- CHANGELOG.md | 1 + rocketpy/simulation/monte_carlo.py | 70 +++++++--- .../test_monte_carlo_determinism.py | 129 ++++++++++++++---- 3 files changed, 149 insertions(+), 51 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e68153f6e..fe6332fdd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -33,6 +33,7 @@ Attention: The newest changes should be on top --> ### Added - ENH: update master with develop [#1081](https://github.com/RocketPy-Team/RocketPy/pull/1081) +- ENH: reproducible Monte Carlo runs via a random_seed argument [#1054](https://github.com/RocketPy-Team/RocketPy/pull/1054) ### Changed diff --git a/rocketpy/simulation/monte_carlo.py b/rocketpy/simulation/monte_carlo.py index 88461de52..ff0fe20fc 100644 --- a/rocketpy/simulation/monte_carlo.py +++ b/rocketpy/simulation/monte_carlo.py @@ -191,15 +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 + 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 - (``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. + 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: @@ -280,19 +280,26 @@ def __setup_files(self, append): @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. + """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 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.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): @@ -315,7 +322,7 @@ def __run_in_serial(self, random_seed=None): # pylint: disable=too-many-stateme Parameters ---------- - random_seed : int, SeedSequence, Generator, optional + random_seed : int or SeedSequence, optional Root seed for the run. See ``simulate``. Returns @@ -369,7 +376,7 @@ def __run_in_parallel(self, n_workers=None, random_seed=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 + random_seed : int or SeedSequence, optional Root seed for the run. See ``simulate``. Returns @@ -464,8 +471,11 @@ def __sim_producer(self, child_seeds, sim_monitor, mutex, error_event): # pylin Event signaling an error occurred during the simulation. """ try: - while sim_monitor.keep_simulating(): - sim_idx = sim_monitor.increment() - 1 + while True: + sim_idx = _claim_next_index(sim_monitor, mutex) + if sim_idx is None: + break + inputs_json, outputs_json = "", "" self.__seed_simulation(child_seeds[sim_idx]) @@ -1665,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. diff --git a/tests/unit/simulation/test_monte_carlo_determinism.py b/tests/unit/simulation/test_monte_carlo_determinism.py index 939980508..526717182 100644 --- a/tests/unit/simulation/test_monte_carlo_determinism.py +++ b/tests/unit/simulation/test_monte_carlo_determinism.py @@ -4,9 +4,9 @@ turning the run's root seed into one independent child stream per simulation index. Two private helpers do the work: -* ``__root_seed_sequence`` normalizes the flexible ``random_seed`` argument (int, - ``SeedSequence``, ``Generator``, ``BitGenerator`` or None) into a - ``SeedSequence`` that can be spawned; +* ``__root_seed_sequence`` normalizes the ``random_seed`` argument (an int, a + sequence of ints, a ``SeedSequence`` or None) into a fresh ``SeedSequence`` + that can be spawned; * ``__seed_simulation`` splits one per-index child seed three ways so the environment, rocket and flight draw from independent streams. @@ -19,12 +19,15 @@ it lets the seeding invariants be asserted without running a Monte Carlo. """ +import threading +import time from types import SimpleNamespace import numpy as np import pytest from rocketpy.simulation import MonteCarlo +from rocketpy.simulation.monte_carlo import _SimMonitor, _claim_next_index _root_seed_sequence = MonteCarlo._MonteCarlo__root_seed_sequence _seed_simulation = MonteCarlo._MonteCarlo__seed_simulation @@ -44,17 +47,17 @@ def _entropy(seed_sequence, n=4): "make_seed", [ pytest.param(lambda: 12345, id="int"), + pytest.param(lambda: np.int64(12345), id="numpy-int"), + pytest.param(lambda: [1, 2, 3], id="sequence"), pytest.param(lambda: np.random.SeedSequence(12345), id="seedsequence"), - pytest.param(lambda: np.random.default_rng(12345), id="generator"), - pytest.param(lambda: np.random.PCG64(12345), id="bitgenerator"), ], ) -def test_root_seed_sequence_accepts_supported_types(make_seed): - """int, SeedSequence, Generator and BitGenerator all normalize to the same - root SeedSequence stream for an equivalent seed value.""" +def test_root_seed_sequence_accepts_seed_like_values(make_seed): + """An int, a numpy integer, a sequence of ints and a SeedSequence are all + accepted, normalize to a SeedSequence, and are reproducible.""" root = _root_seed_sequence(make_seed()) assert isinstance(root, np.random.SeedSequence) - assert _entropy(root) == _entropy(_root_seed_sequence(12345)) + assert _entropy(root) == _entropy(_root_seed_sequence(make_seed())) def test_root_seed_sequence_none_draws_fresh_entropy(): @@ -65,30 +68,42 @@ def test_root_seed_sequence_none_draws_fresh_entropy(): @pytest.mark.parametrize( - "make_seed, resolve", + "make_generator", [ - pytest.param( - lambda: np.random.SeedSequence(999), - lambda seed: seed, - id="seedsequence", - ), - pytest.param( - lambda: np.random.default_rng(999), - lambda seed: seed.bit_generator.seed_seq, - id="generator", - ), - pytest.param( - lambda: np.random.PCG64(999), - lambda seed: seed.seed_seq, - id="bitgenerator", - ), + pytest.param(lambda: np.random.default_rng(999), id="generator"), + pytest.param(lambda: np.random.PCG64(999), id="bitgenerator"), ], ) -def test_root_seed_sequence_reuses_existing_seed_sequence(make_seed, resolve): - """When given something that already carries a SeedSequence, the helper - reuses that object rather than copying it.""" - seed = make_seed() - assert _root_seed_sequence(seed) is resolve(seed) +def test_root_seed_sequence_rejects_stateful_generators(make_generator): + """A Generator/BitGenerator is a stateful RNG, not a seed value, so it is + rejected instead of being reduced to its underlying SeedSequence.""" + with pytest.raises(TypeError, match="SeedSequence"): + _root_seed_sequence(make_generator()) + + +def test_root_seed_sequence_copies_seedsequence_without_mutating_it(): + """A supplied SeedSequence is copied from its full state: repeated calls with + the same object reproduce the same children, the caller's spawn counter is + left untouched, and a spawned child (non-empty spawn_key) round-trips too.""" + + def children(seed_sequence): + return [ + _entropy(child) for child in _root_seed_sequence(seed_sequence).spawn(3) + ] + + # A SeedSequence that has already spawned children, so its counter is not 0. + seed = np.random.SeedSequence(2024) + seed.spawn(5) + counter_before = seed.n_children_spawned + + assert children(seed) == children(seed), "same object twice must reproduce" + assert seed.n_children_spawned == counter_before, "caller must not be mutated" + assert _root_seed_sequence(seed) is not seed, "must return a copy, not the caller" + + # A spawned child carries a non-empty spawn_key that the copy must preserve. + child = np.random.SeedSequence(2024).spawn(1)[0] + assert child.spawn_key != () + assert children(child) == children(child) # --------------------------------------------------------------------------- # @@ -138,3 +153,57 @@ def split(child): return [_entropy(env[0]), _entropy(rocket[0]), _entropy(flight[0])] assert split(np.random.SeedSequence(2024)) == split(np.random.SeedSequence(2024)) + + +# --------------------------------------------------------------------------- # +# _claim_next_index: atomic hand-out of the next simulation index # +# --------------------------------------------------------------------------- # + + +def test_claim_next_index_hands_out_each_index_once_under_contention(): + """Holding the mutex across keep_simulating() and increment() must hand out + each index exactly once, even when every worker reaches the claim together. + + A barrier releases all workers at once and a widened check-to-increment + window would let an unlocked claim run several workers past the count < n + check before any increments; the lock is what keeps the result to exactly + n_simulations indices (0..n-1, none repeated) and the counter from + overshooting. + """ + n_simulations = 5 + n_workers = 8 + monitor = _SimMonitor(initial_count=0, n_simulations=n_simulations, start_time=0.0) + + # Widen the window between the check and the increment so that, without the + # lock, several workers could pass count < n before any of them increments. + real_keep_simulating = monitor.keep_simulating + + def slow_keep_simulating(): + result = real_keep_simulating() + time.sleep(0.02) + return result + + monitor.keep_simulating = slow_keep_simulating + + mutex = threading.Lock() + barrier = threading.Barrier(n_workers) + claimed = [] + claimed_lock = threading.Lock() + + def worker(): + barrier.wait() + while True: + index = _claim_next_index(monitor, mutex) + if index is None: + break + with claimed_lock: + claimed.append(index) + + workers = [threading.Thread(target=worker) for _ in range(n_workers)] + for thread in workers: + thread.start() + for thread in workers: + thread.join() + + assert sorted(claimed) == list(range(n_simulations)) + assert monitor.count == n_simulations From 1446ee723ffd68f42bcae97f9ee672f6454f41c1 Mon Sep 17 00:00:00 2001 From: thc1006 <84045975+thc1006@users.noreply.github.com> Date: Mon, 20 Jul 2026 14:25:02 +0800 Subject: [PATCH 4/8] ENH: derive Monte Carlo per-index seeds in O(1) and seed models with 128-bit ints Each simulation index is seeded from its own child of the run's root seed. Building that child by extending the captured root spawn_key is bit-identical to root.spawn(number_of_simulations)[index] but O(1) in time and memory, so a worker reconstructs any index from a small root state instead of the full spawned list being materialized and pickled to every process. Hand each model a 128-bit int rather than a SeedSequence: a plain int is the seed type accepted alike by numpy.random.default_rng, RandomState and the stdlib random.Random (which rejects a SeedSequence with a TypeError from Python 3.11), so a custom sampler whose reset_seed documents an int keeps working; all four uint32 words are combined by value (not via tobytes) so the seed is byte-order independent and keeps the full 128-bit pool instead of collapsing to 32 bits. The random_seed docstring now lists the accepted types and notes the seeding is informed by SPEC 7 while keeping immutable seed-snapshot semantics. The unit tests assert the SeedSequence copy preserves full .state (an entropy-only copy would fail), the O(1) child equals spawn bit-for-bit -- including a root whose child counter has advanced and indices past 2**32 -- and each model receives a distinct 128-bit int. Signed-off-by: thc1006 <84045975+thc1006@users.noreply.github.com> --- rocketpy/simulation/monte_carlo.py | 121 ++++++++---- .../test_monte_carlo_determinism.py | 187 ++++++++++++++---- 2 files changed, 241 insertions(+), 67 deletions(-) diff --git a/rocketpy/simulation/monte_carlo.py b/rocketpy/simulation/monte_carlo.py index ff0fe20fc..dfb571b00 100644 --- a/rocketpy/simulation/monte_carlo.py +++ b/rocketpy/simulation/monte_carlo.py @@ -191,15 +191,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 or numpy.random.SeedSequence, optional + 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, 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. + 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: @@ -233,6 +240,9 @@ def simulate( self._export_config = kwargs self.number_of_simulations = number_of_simulations self._initial_sim_idx = self.num_of_loaded_sims if append else 0 + # Small, picklable root seed state captured once per run; every + # simulation index derives its child seed from it (see __child_seed). + self.__root_state = None print("Starting Monte Carlo analysis") @@ -302,6 +312,38 @@ def __root_seed_sequence(random_seed): ) 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. @@ -309,12 +351,13 @@ def __seed_simulation(self, child_seed): 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. + 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(env_seed) - self.rocket._set_stochastic(rocket_seed) - self.flight._set_stochastic(flight_seed) + 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, random_seed=None): # pylint: disable=too-many-statements """ @@ -334,15 +377,13 @@ def __run_in_serial(self, random_seed=None): # pylint: disable=too-many-stateme n_simulations=self.number_of_simulations, start_time=time(), ) - child_seeds = self.__root_seed_sequence(random_seed).spawn( - self.number_of_simulations - ) + self.__capture_root_state(random_seed) try: while sim_monitor.keep_simulating(): sim_idx = sim_monitor.increment() - 1 inputs_json, outputs_json = "", "" - self.__seed_simulation(child_seeds[sim_idx]) + 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) @@ -399,19 +440,18 @@ def __run_in_parallel(self, n_workers=None, random_seed=None): ) processes = [] - # 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 - ) + # 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. + self.__capture_root_state(random_seed) for _ in range(n_workers): sim_producer = multiprocess.Process( target=self.__sim_producer, args=( - child_seeds, sim_monitor, mutex, simulation_error_event, @@ -453,16 +493,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, child_seeds, 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 ---------- - 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 @@ -478,7 +513,7 @@ def __sim_producer(self, child_seeds, sim_monitor, mutex, error_event): # pylin inputs_json, outputs_json = "", "" - self.__seed_simulation(child_seeds[sim_idx]) + 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) @@ -1675,14 +1710,34 @@ def export_errors_to_json(self, filename): self._write_log_to_json(self.errors_log, filename) +def _seed_sequence_to_int(seed_sequence): + """Collapse a ``SeedSequence`` into a 128-bit Python ``int``. + + A plain ``int`` is the one seed type accepted alike by + ``numpy.random.default_rng``, ``numpy.random.RandomState`` and the stdlib + ``random.Random`` (which rejects a ``SeedSequence`` with a ``TypeError`` + since Python 3.11), so a custom sampler whose ``reset_seed`` documents an + ``int`` keeps working. All four ``uint32`` words are combined to keep the + full 128-bit pool, so the environment/rocket/flight sub-streams stay + decorrelated instead of collapsing to a single 32-bit word. + + The words are combined by value (little-endian word order), not via + ``tobytes()``, so the seed is the same on big- and little-endian machines + -- a byte-order-dependent seed would break the cross-platform + reproducibility this whole scheme exists to provide. + """ + words = seed_sequence.generate_state(4, dtype=np.uint32) + return sum(int(word) << (32 * position) for position, word in enumerate(words)) + + 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. + either increments, and both then claim an index, running more simulations + than were requested (and duplicating a simulation index). """ mutex.acquire() try: diff --git a/tests/unit/simulation/test_monte_carlo_determinism.py b/tests/unit/simulation/test_monte_carlo_determinism.py index 526717182..546b2c1d1 100644 --- a/tests/unit/simulation/test_monte_carlo_determinism.py +++ b/tests/unit/simulation/test_monte_carlo_determinism.py @@ -2,11 +2,16 @@ ``MonteCarlo.simulate(random_seed=...)`` makes the sampled inputs reproducible by turning the run's root seed into one independent child stream per simulation -index. Two private helpers do the work: +index. Four small helpers do the work: * ``__root_seed_sequence`` normalizes the ``random_seed`` argument (an int, a - sequence of ints, a ``SeedSequence`` or None) into a fresh ``SeedSequence`` - that can be spawned; + sequence of ints, a ``SeedSequence`` or None) into a fresh ``SeedSequence``; +* ``__child_seed`` derives the child seed for one simulation index in O(1) by + extending the captured root ``spawn_key`` -- bit-identical to + ``root.spawn(n)[index]`` but without materializing the whole spawned list, so + a worker can rebuild any index from the small root state alone; +* ``_seed_sequence_to_int`` collapses a child into a 128-bit ``int`` (the seed + type a documented ``CustomSampler.reset_seed`` accepts); * ``__seed_simulation`` splits one per-index child seed three ways so the environment, rocket and flight draw from independent streams. @@ -19,6 +24,8 @@ it lets the seeding invariants be asserted without running a Monte Carlo. """ +import random as stdlib_random +import sys import threading import time from types import SimpleNamespace @@ -27,9 +34,14 @@ import pytest from rocketpy.simulation import MonteCarlo -from rocketpy.simulation.monte_carlo import _SimMonitor, _claim_next_index +from rocketpy.simulation.monte_carlo import ( + _SimMonitor, + _claim_next_index, + _seed_sequence_to_int, +) _root_seed_sequence = MonteCarlo._MonteCarlo__root_seed_sequence +_child_seed = MonteCarlo._MonteCarlo__child_seed _seed_simulation = MonteCarlo._MonteCarlo__seed_simulation @@ -38,6 +50,26 @@ def _entropy(seed_sequence, n=4): return tuple(int(x) for x in seed_sequence.generate_state(n)) +def _plan(root): + """A stand-in ``self`` carrying only the root state ``__child_seed`` reads.""" + return SimpleNamespace( + _MonteCarlo__root_state=( + root.entropy, + root.spawn_key, + root.pool_size, + root.n_children_spawned, + ) + ) + + +def _advanced_root(seed, already_spawned): + """A root whose own child counter has advanced (n_children_spawned != 0), + the state a user's already-spawned SeedSequence would arrive in.""" + root = np.random.SeedSequence(seed) + root.spawn(already_spawned) + return root + + # --------------------------------------------------------------------------- # # __root_seed_sequence: normalizing the flexible seed argument # # --------------------------------------------------------------------------- # @@ -81,33 +113,123 @@ def test_root_seed_sequence_rejects_stateful_generators(make_generator): _root_seed_sequence(make_generator()) -def test_root_seed_sequence_copies_seedsequence_without_mutating_it(): - """A supplied SeedSequence is copied from its full state: repeated calls with - the same object reproduce the same children, the caller's spawn counter is - left untouched, and a spawned child (non-empty spawn_key) round-trips too.""" +def test_root_seed_sequence_copies_full_state_without_mutating_caller(): + """A supplied SeedSequence is copied from its FULL state -- entropy, spawn_key, + pool_size and n_children_spawned -- not just its entropy, and the caller object + is not mutated. Asserting on ``.state`` is what gives this teeth: an + entropy-only copy would silently drop spawn_key/n_children_spawned (making a + spawned-child seed collide with its parent) and fail the state comparison.""" + source = np.random.SeedSequence(2024).spawn(3)[2] # non-empty spawn_key + source.spawn(5) # advance its own child counter, so it is not 0 + assert source.spawn_key == (2,) + assert source.n_children_spawned == 5 + + state_before = dict(source.state) + clone = _root_seed_sequence(source) + + assert clone is not source, "must return a copy, not the caller" + assert clone.state == state_before, "copy must preserve the full seed state" + assert source.state == state_before, "caller must not be mutated" + # The copy reproduces exactly what an independent full-state rebuild produces. + rebuilt = np.random.SeedSequence(**state_before) + assert [_entropy(c) for c in clone.spawn(3)] == [ + _entropy(c) for c in rebuilt.spawn(3) + ] + + +# --------------------------------------------------------------------------- # +# __child_seed: O(1) per-index derivation, bit-identical to spawn(n)[index] # +# --------------------------------------------------------------------------- # + + +@pytest.mark.parametrize( + "make_root", + [ + pytest.param(lambda: np.random.SeedSequence(2024), id="int-root"), + pytest.param(lambda: np.random.SeedSequence([7, 8, 9]), id="sequence-root"), + pytest.param( + lambda: np.random.SeedSequence(2024).spawn(3)[2], id="spawned-root" + ), + pytest.param(lambda: _advanced_root(99, 4), id="advanced-counter-root"), + ], +) +def test_child_seed_matches_spawn_bit_for_bit(make_root): + """Deriving index i by extending the root spawn_key equals ``root.spawn(n)[i]`` + exactly, so the O(1) derivation changes no sampled inputs versus a full spawn. + A fresh, identical root is built on each side so neither run mutates the other. + The ``advanced-counter`` root (n_children_spawned != 0) covers a user passing a + SeedSequence they have already spawned from: the base offset must equal + n_children_spawned or the derived index would collide with those children. + """ + n = 6 + derived = [_entropy(_child_seed(_plan(make_root()), i)) for i in range(n)] + spawned = [_entropy(child) for child in make_root().spawn(n)] + assert derived == spawned + + +def test_child_seed_is_worker_order_independent(): + """Any index maps to the same child regardless of the order indices are asked + for -- the property that makes a run invariant to worker scheduling.""" + plan = _plan(np.random.SeedSequence(2024)) + forward = {i: _entropy(_child_seed(plan, i)) for i in range(5)} + backward = {i: _entropy(_child_seed(plan, i)) for i in reversed(range(5))} + assert forward == backward + + +def test_child_seed_supports_indices_beyond_32_bits(): + """A simulation index past 2**32 is not truncated: it derives a distinct child + from its neighbour and matches the direct spawn_key construction for it.""" + root = np.random.SeedSequence(11) + plan = _plan(root) + big = 2**32 + 5 + assert _entropy(_child_seed(plan, big)) != _entropy(_child_seed(plan, big + 1)) + expected = np.random.SeedSequence( + entropy=root.entropy, spawn_key=(big,), pool_size=root.pool_size + ) + assert _entropy(_child_seed(plan, big)) == _entropy(expected) + + +# --------------------------------------------------------------------------- # +# _seed_sequence_to_int: 128-bit int seed for the samplers # +# --------------------------------------------------------------------------- # + + +def test_seed_sequence_to_int_is_deterministic_128_bit_int(): + def child(): + return np.random.SeedSequence(42).spawn(1)[0] + + seed = _seed_sequence_to_int(child()) + assert isinstance(seed, int) + assert 0 <= seed < 2**128 + assert seed == _seed_sequence_to_int(child()), "must be deterministic" - def children(seed_sequence): - return [ - _entropy(child) for child in _root_seed_sequence(seed_sequence).spawn(3) - ] - # A SeedSequence that has already spawned children, so its counter is not 0. - seed = np.random.SeedSequence(2024) - seed.spawn(5) - counter_before = seed.n_children_spawned +def test_seed_sequence_to_int_uses_all_128_bits(): + """The int combines all four uint32 words, not a single 32-bit word, so it + keeps the full entropy pool rather than collapsing collision risk to n**2 / + 2**32. A single-word reduction would compare unequal here.""" + ss = np.random.SeedSequence(42).spawn(1)[0] + one_word = int(np.random.SeedSequence(42).spawn(1)[0].generate_state(1)[0]) + assert _seed_sequence_to_int(ss) != one_word + assert _seed_sequence_to_int(ss).bit_length() > 32 - assert children(seed) == children(seed), "same object twice must reproduce" - assert seed.n_children_spawned == counter_before, "caller must not be mutated" - assert _root_seed_sequence(seed) is not seed, "must return a copy, not the caller" - # A spawned child carries a non-empty spawn_key that the copy must preserve. - child = np.random.SeedSequence(2024).spawn(1)[0] - assert child.spawn_key != () - assert children(child) == children(child) +def test_seed_int_is_accepted_by_the_modern_rng_apis(): + """The 128-bit int a sampler receives works with random.Random and + numpy.random.default_rng -- the paths a CustomSampler uses. Passing a + SeedSequence there instead is unsafe: from Python 3.11 random.Random rejects + it with a TypeError, and before 3.11 it is silently hashed rather than used as + entropy. Either way an int is the right thing to hand a sampler.""" + seed = _seed_sequence_to_int(np.random.SeedSequence(1).spawn(1)[0]) + assert isinstance(stdlib_random.Random(seed).random(), float) + assert np.random.default_rng(seed).random() is not None + if sys.version_info >= (3, 11): + with pytest.raises(TypeError): + stdlib_random.Random(np.random.SeedSequence(1)) # --------------------------------------------------------------------------- # -# __seed_simulation: splitting one child seed across the three models # +# __seed_simulation: splitting one child seed across the three models # # --------------------------------------------------------------------------- # @@ -132,17 +254,14 @@ def _split_seeds(child_seed): return models.environment.seeds, models.rocket.seeds, models.flight.seeds -def test_seed_simulation_decorrelates_env_rocket_flight(): - """The per-index child seed is split three ways so environment, rocket and - flight draw from independent streams instead of sharing one.""" +def test_seed_simulation_hands_each_model_a_distinct_128_bit_int(): + """The per-index child seed is split three ways, and each model receives a + plain 128-bit int (not a SeedSequence) from an independent stream.""" env_seeds, rocket_seeds, flight_seeds = _split_seeds(np.random.SeedSequence(2024)) assert [len(env_seeds), len(rocket_seeds), len(flight_seeds)] == [1, 1, 1] - fingerprints = { - _entropy(env_seeds[0]), - _entropy(rocket_seeds[0]), - _entropy(flight_seeds[0]), - } - assert len(fingerprints) == 3 + seeds = [env_seeds[0], rocket_seeds[0], flight_seeds[0]] + assert all(isinstance(s, int) and 0 <= s < 2**128 for s in seeds) + assert len(set(seeds)) == 3, "env/rocket/flight must be decorrelated" def test_seed_simulation_is_deterministic_per_child(): @@ -150,7 +269,7 @@ def test_seed_simulation_is_deterministic_per_child(): def split(child): env, rocket, flight = _split_seeds(child) - return [_entropy(env[0]), _entropy(rocket[0]), _entropy(flight[0])] + return [env[0], rocket[0], flight[0]] assert split(np.random.SeedSequence(2024)) == split(np.random.SeedSequence(2024)) From d840d543872032867b91d49398d724745621c295 Mon Sep 17 00:00:00 2001 From: thc1006 <84045975+thc1006@users.noreply.github.com> Date: Mon, 20 Jul 2026 14:25:02 +0800 Subject: [PATCH 5/8] BUG: sample list-valued stochastic attributes through the seeded generator dict_generator drew list-valued attributes with the stdlib random.choice, which reads from an unseeded global Random instance, so random_seed did not make those attributes reproducible. Draw the index from this model's seeded numpy generator instead. Indexing (not numpy.random.choice) also avoids coercing a heterogeneous list -- Function objects, paths, arrays -- to a single dtype. Adds a unit test that a list-valued attribute is reproducible under a fixed seed and that heterogeneous entries are returned unchanged. Signed-off-by: thc1006 <84045975+thc1006@users.noreply.github.com> --- rocketpy/stochastic/stochastic_model.py | 13 ++++++-- .../unit/stochastic/test_stochastic_model.py | 33 +++++++++++++++++++ 2 files changed, 43 insertions(+), 3 deletions(-) diff --git a/rocketpy/stochastic/stochastic_model.py b/rocketpy/stochastic/stochastic_model.py index ca26f6578..fd96124cb 100644 --- a/rocketpy/stochastic/stochastic_model.py +++ b/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 @@ -532,7 +530,16 @@ 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 + # Draw the index from this model's seeded generator so a + # list-valued attribute is reproducible under random_seed. The + # stdlib random.choice draws from an unseeded global instance, + # and numpy's choice coerces a heterogeneous list (Function, + # paths, arrays) to a single dtype; indexing avoids both. + if value: + index = int(self.__random_number_generator.integers(len(value))) + generated_dict[arg] = value[index] + else: + generated_dict[arg] = value elif isinstance(value, CustomSampler): try: generated_dict[arg] = value.sample(n_samples=1)[0] diff --git a/tests/unit/stochastic/test_stochastic_model.py b/tests/unit/stochastic/test_stochastic_model.py index 9e35a5330..cf802f767 100644 --- a/tests/unit/stochastic/test_stochastic_model.py +++ b/tests/unit/stochastic/test_stochastic_model.py @@ -1,5 +1,38 @@ +from types import SimpleNamespace + import pytest +from rocketpy.stochastic.stochastic_model import StochasticModel + + +def _sampled_option(model): + """Return the value ``dict_generator`` picks for the ``options`` attribute.""" + return next(model.dict_generator())["options"] + + +def test_list_attribute_sampling_is_reproducible_under_seed(): + """A list-valued stochastic attribute is drawn through the model's own seeded + numpy generator, so a fixed seed reproduces the choice. It used to be drawn + with the stdlib ``random.choice`` (an unseeded global instance), which + ``random_seed`` could not govern. Heterogeneous entries (paths, callables, + lists) are returned unchanged rather than coerced to a numpy dtype the way + ``numpy.random.choice`` would. + """ + options = ["/motor/a.eng", "/motor/b.eng", (lambda t: t), [1, 2, 3]] + model = StochasticModel(obj=SimpleNamespace(), options=options) + + model._set_stochastic(42) + first = _sampled_option(model) + model._set_stochastic(42) + assert _sampled_option(model) == first, "same seed must reproduce the choice" + assert any(first is option for option in options), "object returned unchanged" + + chosen_ids = set() + for seed in range(16): + model._set_stochastic(seed) + chosen_ids.add(id(_sampled_option(model))) + assert len(chosen_ids) > 1, "different seeds must be able to pick differently" + @pytest.mark.parametrize( "fixture_name", From 1e4670d60562c7d519a99eb926ed052fcb829c1e Mon Sep 17 00:00:00 2001 From: thc1006 <84045975+thc1006@users.noreply.github.com> Date: Mon, 20 Jul 2026 14:25:02 +0800 Subject: [PATCH 6/8] TST: verify Monte Carlo seed derivation is start-method invariant The existing worker-invariance test stubs the module-level Flight and so reaches workers only under fork. Add a test that the per-index seed derived in a worker matches the main process under every available start method (fork, spawn, forkserver), using a top-level picklable target and small picklable arguments so it is valid under spawn/forkserver -- Python 3.14's POSIX default -- without relying on inherited parent state. It runs in ordinary CI. Signed-off-by: thc1006 <84045975+thc1006@users.noreply.github.com> --- .../test_monte_carlo_determinism.py | 77 +++++++++++++++++-- 1 file changed, 72 insertions(+), 5 deletions(-) diff --git a/tests/integration/simulation/test_monte_carlo_determinism.py b/tests/integration/simulation/test_monte_carlo_determinism.py index b47629a55..f0355ca71 100644 --- a/tests/integration/simulation/test_monte_carlo_determinism.py +++ b/tests/integration/simulation/test_monte_carlo_determinism.py @@ -14,21 +14,52 @@ under the ``fork`` start method, so the worker-invariance test skips otherwise and is marked ``slow`` to match the other Monte Carlo multiprocessing tests. -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. +A dedicated numpy-only rocket keeps the fork-based end-to-end test simple: it +gives the motor a single ``thrust_source`` so the run has no list-valued attribute +at all. List sampling is itself seeded now (it draws through the model generator, +not the stdlib ``random.choice``) and is covered directly in +``tests/unit/stochastic/test_stochastic_model``. + +Seed derivation being independent of the multiprocessing start method (fork, +spawn or forkserver) is verified separately by +``test_seed_derivation_is_start_method_invariant``, which uses a top-level +picklable target so it is safe under ``spawn``/``forkserver`` -- unlike the +``Flight``-stub test above, which reaches workers only under ``fork``. """ import json +import multiprocessing +from types import SimpleNamespace +import numpy as np import pytest import rocketpy.simulation.monte_carlo as mc_module from rocketpy.simulation import MonteCarlo +from rocketpy.simulation.monte_carlo import _seed_sequence_to_int from rocketpy.stochastic import StochasticRocket, StochasticSolidMotor +_child_seed = MonteCarlo._MonteCarlo__child_seed + + +def _available_start_methods(): + """The multiprocessing start methods this platform actually supports.""" + supported = multiprocessing.get_all_start_methods() + return [method for method in ("fork", "spawn", "forkserver") if method in supported] + + +def _derive_index_seeds(root_state, indices): + """Derive the per-index seed fingerprints from ``root_state``. + + Top-level and picklable (only a small tuple and a list of ints cross the + process boundary), so it runs unchanged under every start method -- including + ``spawn``/``forkserver``, which re-import this module rather than inheriting + the parent's memory. It calls the real production helpers (``__child_seed`` + and ``_seed_sequence_to_int``) so the test tracks the shipped derivation. + """ + plan = SimpleNamespace(_MonteCarlo__root_state=root_state) + return {index: _seed_sequence_to_int(_child_seed(plan, index)) for index in indices} + class _StubFlight: """Minimal stand-in for ``Flight`` that skips trajectory integration.""" @@ -175,3 +206,39 @@ def test_inputs_are_worker_invariant( 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}" + + +@pytest.mark.parametrize("start_method", _available_start_methods()) +def test_seed_derivation_is_start_method_invariant(start_method): + """Per-index seeds derived in a worker match the main process under every + available start method (fork, spawn, forkserver). + + The full worker-invariance test above stubs the module-level ``Flight`` and so + only reaches workers under ``fork``. This one instead checks the property that + actually has to hold cross-platform -- that a simulation index maps to the same + seed no matter which process derives it -- using a top-level picklable target + and small picklable arguments, so it is valid under ``spawn``/``forkserver`` + (Python 3.14's POSIX default) without relying on any inherited parent state. + Two workers split the indices; their combined result must equal the + single-process derivation. + """ + root = np.random.SeedSequence(2718281828) + root_state = ( + root.entropy, + root.spawn_key, + root.pool_size, + root.n_children_spawned, + ) + indices = list(range(6)) + expected = _derive_index_seeds(root_state, indices) + + context = multiprocessing.get_context(start_method) + chunks = [(root_state, indices[0::2]), (root_state, indices[1::2])] + with context.Pool(2) as pool: + results = pool.starmap(_derive_index_seeds, chunks) + + combined = {} + for result in results: + combined.update(result) + assert combined == expected + assert sorted(combined) == indices From 5ba7598c821dc8b9a69f15296c5f9e363e89e961 Mon Sep 17 00:00:00 2001 From: thc1006 <84045975+thc1006@users.noreply.github.com> Date: Mon, 20 Jul 2026 15:59:29 +0800 Subject: [PATCH 7/8] BUG: seed list/position sampling and decorrelate rocket components StochasticModel list-valued attributes were sampled with the stdlib random.choice (an unseeded global) and StochasticRocket._randomize_position did the same for list-valued component positions, so random_seed did not govern either. Both now draw the index through the model's seeded generator via a shared _random_choice helper -- indexing, not numpy.random.choice, so heterogeneous objects (Function, paths, arrays) stay intact. StochasticRocket._set_stochastic also handed the same seed to the rocket body and every surface, motor, rail button and parachute, so components sampling the same distribution drew identical values (a main and a drogue parachute got the same cd_s and lag quantiles). Each component is now reseeded from its own spawned child of a SeedSequence root, in a fixed order, so they stay independent and reproducible under random_seed. Signed-off-by: thc1006 <84045975+thc1006@users.noreply.github.com> --- rocketpy/stochastic/stochastic_model.py | 26 +++++++---- rocketpy/stochastic/stochastic_rocket.py | 37 +++++++++------ .../test_stochastic_rocket_seeding.py | 46 +++++++++++++++++++ 3 files changed, 86 insertions(+), 23 deletions(-) create mode 100644 tests/unit/stochastic/test_stochastic_rocket_seeding.py diff --git a/rocketpy/stochastic/stochastic_model.py b/rocketpy/stochastic/stochastic_model.py index fd96124cb..676bf2b73 100644 --- a/rocketpy/stochastic/stochastic_model.py +++ b/rocketpy/stochastic/stochastic_model.py @@ -506,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. @@ -530,16 +545,7 @@ def dict_generator(self): dist_sampler = value[-1] generated_dict[arg] = dist_sampler(value[0], value[1]) elif isinstance(value, list): - # Draw the index from this model's seeded generator so a - # list-valued attribute is reproducible under random_seed. The - # stdlib random.choice draws from an unseeded global instance, - # and numpy's choice coerces a heterogeneous list (Function, - # paths, arrays) to a single dtype; indexing avoids both. - if value: - index = int(self.__random_number_generator.integers(len(value))) - generated_dict[arg] = value[index] - else: - generated_dict[arg] = value + generated_dict[arg] = self._random_choice(value) elif isinstance(value, CustomSampler): try: generated_dict[arg] = value.sample(n_samples=1)[0] diff --git a/rocketpy/stochastic/stochastic_rocket.py b/rocketpy/stochastic/stochastic_rocket.py index 794a66c85..d349ec467 100644 --- a/rocketpy/stochastic/stochastic_rocket.py +++ b/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 ---------- diff --git a/tests/unit/stochastic/test_stochastic_rocket_seeding.py b/tests/unit/stochastic/test_stochastic_rocket_seeding.py new file mode 100644 index 000000000..83704eb90 --- /dev/null +++ b/tests/unit/stochastic/test_stochastic_rocket_seeding.py @@ -0,0 +1,46 @@ +"""Nested StochasticRocket components are reseeded from distinct SeedSequence +children, so components that sample the same distribution (a main and a drogue +parachute, for example) do not draw identical values. Reproducible under a fixed +seed. See the seeding design in ``StochasticRocket._set_stochastic``. +""" + +from rocketpy.stochastic.stochastic_model import StochasticModel + +# Captured once, before any patching, so wrapping it repeatedly in one test does +# not stack (each recorder wraps the real method, not a previous recorder). +_REAL_SET_STOCHASTIC = StochasticModel._set_stochastic + + +def _record_component_seeds(monkeypatch, rocket, seed): + """Return the seeds handed to every nested component for one reseed.""" + recorded = [] + + def recording(self, seed=None): + recorded.append(seed) + return _REAL_SET_STOCHASTIC(self, seed) + + monkeypatch.setattr(StochasticModel, "_set_stochastic", recording) + rocket._set_stochastic(seed) + return recorded + + +def test_rocket_components_receive_distinct_seeds(monkeypatch, stochastic_calisto): + """Every nested component (body, aerodynamic surfaces, motor, rail buttons and + the two parachutes) is reseeded from its own child, so none collide.""" + seeds = _record_component_seeds(monkeypatch, stochastic_calisto, 42) + + assert len(seeds) > 3, "expected the rocket body plus several components" + assert len(seeds) == len(set(seeds)), ( + "components share a seed -- they would draw perfectly correlated samples" + ) + + +def test_rocket_component_seeds_are_reproducible(monkeypatch, stochastic_calisto): + """The same root seed reseeds every component identically; a different root + seed changes them.""" + first = _record_component_seeds(monkeypatch, stochastic_calisto, 42) + again = _record_component_seeds(monkeypatch, stochastic_calisto, 42) + different = _record_component_seeds(monkeypatch, stochastic_calisto, 43) + + assert again == first, "same seed must reproduce every component seed" + assert different != first, "a different seed must change the component seeds" From e88f4df39102b93c595260df59d5df5d6dd5ebf8 Mon Sep 17 00:00:00 2001 From: thc1006 <84045975+thc1006@users.noreply.github.com> Date: Mon, 20 Jul 2026 15:59:29 +0800 Subject: [PATCH 8/8] BUG: validate the run seed before truncating output; tidy the seed helper simulate() set up (and, for append=False, truncated with w+) the input/output/error files before the seed was validated, so passing a rejected seed such as a Generator destroyed a previous run's results on the way to raising a TypeError. The seed is now captured and validated before __setup_files runs. Moved _seed_sequence_to_int to rocketpy.tools so the stochastic models can share it, and corrected its docstring: a 128-bit int is accepted by default_rng and random.Random, but the legacy RandomState caps a single-int seed at 2**32-1, so the earlier 'accepted by RandomState' claim was wrong. Signed-off-by: thc1006 <84045975+thc1006@users.noreply.github.com> --- rocketpy/simulation/monte_carlo.py | 54 +++++++------------ rocketpy/tools.py | 23 ++++++++ .../test_monte_carlo_determinism.py | 30 +++++++++++ 3 files changed, 71 insertions(+), 36 deletions(-) diff --git a/rocketpy/simulation/monte_carlo.py b/rocketpy/simulation/monte_carlo.py index dfb571b00..2e97bbd65 100644 --- a/rocketpy/simulation/monte_carlo.py +++ b/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, @@ -240,18 +241,22 @@ def simulate( self._export_config = kwargs self.number_of_simulations = number_of_simulations self._initial_sim_idx = self.num_of_loaded_sims if append else 0 - # Small, picklable root seed state captured once per run; every - # simulation index derives its child seed from it (see __child_seed). - self.__root_state = None + + # 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) if parallel: - self.__run_in_parallel(n_workers, random_seed) + self.__run_in_parallel(n_workers) else: - self.__run_in_serial(random_seed) + self.__run_in_serial() self.__terminate_simulation() @@ -359,14 +364,12 @@ def __seed_simulation(self, child_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, random_seed=None): # pylint: disable=too-many-statements + def __run_in_serial(self): # 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``. + The root seed state is captured by ``simulate`` before this runs, so each + simulation index derives its child seed from ``self.__root_state``. Returns ------- @@ -377,7 +380,6 @@ def __run_in_serial(self, random_seed=None): # pylint: disable=too-many-stateme n_simulations=self.number_of_simulations, start_time=time(), ) - self.__capture_root_state(random_seed) try: while sim_monitor.keep_simulating(): sim_idx = sim_monitor.increment() - 1 @@ -408,17 +410,19 @@ def __run_in_serial(self, random_seed=None): # pylint: disable=too-many-stateme f.write(inputs_json) raise error - def __run_in_parallel(self, n_workers=None, random_seed=None): + 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 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 ------- @@ -446,8 +450,6 @@ def __run_in_parallel(self, n_workers=None, random_seed=None): # 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. - self.__capture_root_state(random_seed) - for _ in range(n_workers): sim_producer = multiprocess.Process( target=self.__sim_producer, @@ -1710,26 +1712,6 @@ def export_errors_to_json(self, filename): self._write_log_to_json(self.errors_log, filename) -def _seed_sequence_to_int(seed_sequence): - """Collapse a ``SeedSequence`` into a 128-bit Python ``int``. - - A plain ``int`` is the one seed type accepted alike by - ``numpy.random.default_rng``, ``numpy.random.RandomState`` and the stdlib - ``random.Random`` (which rejects a ``SeedSequence`` with a ``TypeError`` - since Python 3.11), so a custom sampler whose ``reset_seed`` documents an - ``int`` keeps working. All four ``uint32`` words are combined to keep the - full 128-bit pool, so the environment/rocket/flight sub-streams stay - decorrelated instead of collapsing to a single 32-bit word. - - The words are combined by value (little-endian word order), not via - ``tobytes()``, so the seed is the same on big- and little-endian machines - -- a byte-order-dependent seed would break the cross-platform - reproducibility this whole scheme exists to provide. - """ - words = seed_sequence.generate_state(4, dtype=np.uint32) - return sum(int(word) << (32 * position) for position, word in enumerate(words)) - - def _claim_next_index(sim_monitor, mutex): """Atomically claim the next 0-based simulation index, or ``None`` if done. diff --git a/rocketpy/tools.py b/rocketpy/tools.py index 0d7f1a74e..9df900eb5 100644 --- a/rocketpy/tools.py +++ b/rocketpy/tools.py @@ -1467,6 +1467,29 @@ def find_obj_from_hash(obj, hash_, depth_limit=None): return None +def _seed_sequence_to_int(seed_sequence): + """Collapse a ``SeedSequence`` into a 128-bit Python ``int``. + + A plain ``int`` is what ``numpy.random.default_rng`` and the stdlib + ``random.Random`` both accept (``random.Random`` rejects a ``SeedSequence`` + with a ``TypeError`` since Python 3.11), so a custom sampler whose + ``reset_seed`` documents an ``int`` and builds a modern generator keeps + working. The legacy ``numpy.random.RandomState`` is the exception: it caps a + single-integer seed at ``2**32 - 1``, so a sampler still built on it would + have to reduce the value (``RandomState`` is a frozen legacy API NumPy steers + new code away from). All four ``uint32`` words are combined to keep the full + 128-bit pool, so sub-streams stay decorrelated instead of collapsing to a + single 32-bit word. + + The words are combined by value (little-endian word order), not via + ``tobytes()``, so the seed is the same on big- and little-endian machines -- + a byte-order-dependent seed would break the cross-platform reproducibility + this exists to provide. + """ + words = seed_sequence.generate_state(4, dtype=np.uint32) + return sum(int(word) << (32 * position) for position, word in enumerate(words)) + + if __name__ == "__main__": # pragma: no cover import doctest diff --git a/tests/integration/simulation/test_monte_carlo_determinism.py b/tests/integration/simulation/test_monte_carlo_determinism.py index f0355ca71..2d2cfe756 100644 --- a/tests/integration/simulation/test_monte_carlo_determinism.py +++ b/tests/integration/simulation/test_monte_carlo_determinism.py @@ -150,6 +150,36 @@ def _simulate_inputs( return _read_inputs_by_index(montecarlo.input_file) +def test_invalid_seed_does_not_truncate_existing_output( + monkeypatch, + tmp_path, + stochastic_environment, + stochastic_calisto_numpy_only, + stochastic_flight, +): + """A rejected seed must fail before any output file is truncated, so passing + an invalid seed cannot destroy the results of a previous run.""" + monkeypatch.setattr(mc_module, "Flight", _StubFlight) + montecarlo = MonteCarlo( + filename=str(tmp_path / "keep"), + environment=stochastic_environment, + rocket=stochastic_calisto_numpy_only, + flight=stochastic_flight, + ) + with open(montecarlo.input_file, "w", encoding="utf-8") as existing: + existing.write("previous results\n") + + # A Generator is not a seed and is rejected; the run must raise before the + # ``w+`` file setup truncates anything. + with pytest.raises(TypeError): + montecarlo.simulate( + number_of_simulations=3, random_seed=np.random.default_rng(0) + ) + + with open(montecarlo.input_file, encoding="utf-8") as kept: + assert kept.read() == "previous results\n" + + def test_serial_inputs_are_reproducible( monkeypatch, tmp_path,