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

[thc1006]

Pull request type

• Code changes (bugfix, features)

Current behavior

MonteCarlo.simulate() seeds the stochastic models per worker in parallel mode (from a fresh, unseeded np.random.SeedSequence().spawn(n_workers)) and once at construction in serial mode. So the sampled inputs depend on the execution mode and the number of workers, and parallel runs are not reproducible run to run. This is #1053.

New behavior

Adds a keyword-only random_seed to simulate(). From that root, simulation index i is seeded from its own child of the root seed, derived before simulation i runs, so index i maps to the same seed no matter which worker runs it. The sampled inputs come out identical across serial, parallel(2) and parallel(N), and reproducible from the seed.

A few specifics:

• O(1) per-index derivation. The child for index i is built by extending the root's spawn_key, which is exactly how SeedSequence.spawn derives it, so child(i) is bit-identical to root.spawn(number_of_simulations)[i]. Nothing pre-spawns a full list: a worker reconstructs any index from a small root state (entropy, spawn_key, pool_size, counter) that travels with the pickled instance, so nothing O(N) is sent to each process.
• 128-bit int seeds. Each model is reseeded with a plain 128-bit int, not a SeedSequence. An int is the seed type numpy.random.default_rng and the stdlib random.Random both accept (a SeedSequence raises TypeError in random.Random since Python 3.11), so a custom sampler whose reset_seed documents an int keeps working. All four uint32 words are combined by value, so the seed is byte-order independent and keeps the full 128-bit pool rather than collapsing to 32 bits.
• List-valued attributes are now seeded too. StochasticModel.dict_generator drew list attributes with the stdlib random.choice (an unseeded global instance), so random_seed did not govern them. It now draws the index from the model's own seeded generator, which also avoids numpy.random.choice coercing a heterogeneous list (Function, paths, arrays) to a single dtype.

random_seed is a seed, not a live RNG: it takes an int, a numpy integer, a sequence of ints, or a SeedSequence, with None = fresh entropy so existing behavior is unchanged unless you pass a seed. A supplied SeedSequence is copied from its full state before use, so it is never mutated and repeated calls with the same object reproduce the same run. This is informed by SPEC 7 and NumPy's parallel idiom, but keeps immutable seed-snapshot semantics rather than SPEC 7's stateful rng: a Generator/BitGenerator is not accepted, because reducing it to its underlying SeedSequence would ignore how far it has been consumed. Pass rng.bit_generator.seed_seq to seed from an existing generator.

Relation to #1071

#1071 targets the same issue. This PR takes the two ideas it got right, deriving each index's seed on demand instead of pre-spawning a list, and handing the samplers a plain int, and combines them with the parallel-claim lock below, the full 128-bit width (a single 32-bit word collides near 2**16 streams), the list-sampling fix, and cross-platform tests. Happy to reconcile the two however the maintainers prefer.

Tests

tests/unit/simulation/test_monte_carlo_determinism.py unit-tests the seed handling: the accepted seed types; the SeedSequence copy preserving the full .state (an entropy-only copy fails the test); the O(1) child equal to spawn bit-for-bit, including a root whose counter has advanced and indices past 2**32; the 128-bit int width; and the parallel index claim (a deterministic barrier-based test that fails if the claim's lock is dropped). tests/integration/simulation/test_monte_carlo_determinism.py covers end-to-end reproducibility (serial, and serial == parallel(2) == parallel(4)) with a stubbed Flight, plus a start-method-invariance test that runs under fork, spawn and forkserver in ordinary CI. tests/unit/stochastic/test_stochastic_model.py covers the list-sampling fix.

Notes from review

• Parallel workers claimed the next index with an unlocked keep_simulating() + increment(). Near the end of a run two workers could both pass count < n and then both claim an index, running past the requested count. The claim now holds the shared mutex across the check and the increment, so each index is handed out once.
• A supplied SeedSequence was returned as-is, and spawn() advances its child counter, so passing the same object twice was not reproducible. It is now copied from its full state, and Generator/BitGenerator are no longer accepted (see above).

Known limitation

The 128-bit int does not fit the legacy numpy.random.RandomState, which caps seeds at 2**32 - 1. A custom sampler built on the modern default_rng (or the stdlib random.Random) takes it fine; one built on RandomState would need to reduce it. Since RandomState is the discouraged legacy path this felt like the right trade for keeping the full 128-bit decorrelation, but I am happy to revisit if you would rather cap the width.

Breaking change

• Yes

The exact numbers a run produces change (per-index seeding, the env/rocket/flight decorrelation, the 128-bit int seeds, the serial index now counting from 0 to match parallel, and list attributes now sampled through the seeded generator), so external code that pinned exact Monte Carlo samples would need to re-baseline. The in-repo Monte Carlo tests do not pin exact values (test_monte_carlo_simulate checks apogee and impact velocity within a tolerance and still passes), and random_seed is opt-in.

Closes #1053

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[codecov]

Codecov Report

❌ Patch coverage is 84.00000% with 8 lines in your changes missing coverage. Please review.
✅ Project coverage is 82.89%. Comparing base (e0ff281) to head (e5ba940).
⚠️ Report is 7 commits behind head on develop.

Files with missing lines	Patch %	Lines
rocketpy/simulation/monte_carlo.py	82.60%	8 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           develop    #1054      +/-   ##
===========================================
+ Coverage    82.18%   82.89%   +0.71%     
===========================================
  Files          122      122              
  Lines        16355    16410      +55     
===========================================
+ Hits         13441    13603     +162     
+ Misses        2914     2807     -107     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[thc1006]

The seeding logic is unit-tested in tests/unit/simulation/test_monte_carlo_determinism.py: __root_seed_sequence accepts an int, a sequence of ints, or a SeedSequence (copied from its full state so the caller's object is not mutated and repeated calls reproduce) and rejects a stateful Generator/BitGenerator; __seed_simulation splits each child seed three ways; and _claim_next_index (the atomic index claim from the race fix) has a deterministic barrier-based test that over-claims and fails if the lock is removed. End-to-end reproducibility (serial, and serial == parallel) lives in tests/integration/, with only the fork-based worker-invariance test marked slow.

The lines codecov still shows uncovered are all in the parallel path: simulate's parallel=True dispatch, the worker setup in __run_in_parallel, and the __sim_producer loop. The coverage jobs cannot reach them because parallel=True is only exercised by the slow worker-invariance test (the jobs do not pass --runslow), and the producer body runs in forked worker processes that coverage.py does not instrument without concurrency = multiprocessing. The behavior is covered by the slow determinism and test_monte_carlo_simulate[parallel] tests, and the claim logic by the fast unit test above. Glad to set up multiprocessing coverage separately if you want the parallel path counted, but that felt out of scope for this PR.

[phmbressan]

The implementation is very clear and throughout, nice work.

The explanation on the concepts behind per index seeding (both in the issue and PR description) were rather helpful. I agree having reproducible results was an issue with the parallel per worker seeding.

Regarding the decisions on parameter naming, I agree with most of the decisions taken here. Moreover, the rng attribute is well docstringed, so it shouldn't be a matter of confusion to the user.

@MateusStano could you give your two cents on the changes here before we proceed with a merge?

[phmbressan]

Could you double check if after running this test, no files are left dangling on the test file system. If that is the case, using a try-finally block to remove those (as done in a few other MonteCarlo tests) is an option.

[thc1006]

Thanks for the review. On the dangling files: I checked, and nothing gets left behind. The test passes filename=tmp_path / tag into MonteCarlo, so the .inputs.txt, .outputs.txt and .errors.txt all land under pytest's per-test tmp_path and get cleaned up with it. I reran it and confirmed the repo and working directory stay clean, so there's nothing here for a try/finally to remove.

The other MonteCarlo tests need that cleanup because their fixture uses a fixed filename="monte_carlo_test", which writes into the cwd. This one keeps everything inside tmp_path instead. I'm happy to match the try/finally style if you'd prefer consistency across the file, but the tmp_path route seemed cleaner to me. Your call.

[MateusStano]

I am not sure about this so @thc1006 could you verify, please?

But I believe adding this child_seeds[sim_idx] can break at the very end of the simulation. If you are running 8 simulations, and we are at step 7. We could have one worker check 7<8 -> yes, and a second worker also check 7<8 -> yes, because keep_simulating() and increment() are two separate steps and nothing locks them together. Then both workers increment: the first gets sim_idx = 7 (ok), but the second gets sim_idx = 8. Since child_seeds only has 8 item, child_seeds[8] does not exist

I think this race was already there before, but it didnt crash the simulation. Could you check and fix, please?

[thc1006]

You're right, this is a real race. keep_simulating() and increment() are separate manager calls, so near the end multiple workers can pass the count < n check before any of them increments and then claim sim_idx == n, which the new child_seeds[sim_idx] lookup turned into an IndexError (before, it just wrote an extra record). Fixed in 9c020b6: the claim now goes through a small _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 runs past number_of_simulations. There's a deterministic unit test for it (a barrier plus a widened check-to-increment window) that over-claims and fails if the lock is removed. Thanks for catching it.

[MateusStano]

This .spawn() call seems to make the case of using SeedSequence/Generator not repeatable

From what I understand if a user builds a SeedSequence (or Generator) once and passes the same object to simulate() twice, the first will spawns a set of children and the second run will spawn another set

[thc1006]

Agreed, the current code mutated the caller's seed lineage. SeedSequence.spawn() advances n_children_spawned, and the helper returned the original object, so a second simulate() with the same object spawned a different set. Fixed in 9c020b6 with two changes:

1. A supplied SeedSequence is now copied from its full state (entropy, spawn_key, pool_size, n_children_spawned) before spawning, so the caller's object is untouched and repeated calls reproduce the same children. Copying only entropy would drop spawn_key, which matters for a spawned child.
2. I dropped Generator/BitGenerator from the accepted types. A generator is a stateful RNG, not a seed: its bit_generator.seed_seq is only the original lineage, so a generator advanced by a million draws would still seed the run identically to a fresh one. SPEC 7's rng semantics would instead consume it, the opposite of a reproducible seed. random_seed now takes an int, a sequence of ints, or a SeedSequence; a generator user can pass rng.bit_generator.seed_seq.

Tests cover the same-object reproducibility, n_children_spawned staying put, a spawned child with a non-empty spawn_key, and the generator rejection.

[thc1006]

@MateusStano friendly ping when you have a moment. Both points from your last pass are addressed: the parallel index claim now holds the shared mutex across the check-and-increment (with a deterministic test that goes red if the lock is removed), and a supplied SeedSequence is copied from its full state before spawning, so repeated calls reproduce and the caller is left untouched. I replied inline on both threads. The test matrix and lint pass on the current head; the only red is the soft codecov patch check, which is the parallel-only lines I covered in the thread above. Whenever you get a chance to take another look, I would appreciate it.

[thc1006]

Heads up that this changed enough since the last look to be worth a fresh pass rather than merging on the earlier approval. @MateusStano @phmbressan when you have a moment.

What is new since the review:

• Per-index child seeds are now derived in O(1) by extending the root spawn_key (bit-identical to spawn(n)[i]) instead of pre-spawning the whole list, so nothing O(N) is pickled to each worker.
• The samplers get a plain 128-bit int rather than a SeedSequence, so a custom sampler's int-typed reset_seed keeps working (a SeedSequence raises TypeError in random.Random since 3.11). The int combines all four words by value, so it is byte-order independent.
• List-valued stochastic attributes now draw from the model's seeded generator, so random_seed governs them too. That closes the gap the previous description called out as a known limitation.
• Added a start-method-invariance test that runs under fork, spawn and forkserver in ordinary CI, since 3.14 moved the POSIX default to forkserver.

Both earlier concerns are still handled: the parallel claim holds the mutex across the check and the increment, and a supplied SeedSequence is copied from its full state. #1071 opened for the same issue in the meantime; the description notes how this relates and what it borrows. A re-review whenever you get the chance would be appreciated.