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

[thc1006]

Scope, up front, because three reviews have now asked me not to claim things this does not claim.

This makes the sampled inputs reproducible: for a given root seed, simulation index i draws the same stochastic parameters whether the run is serial or parallel and however many workers it uses. That is what .inputs.txt records and what the tests compare.

It does not make the flown inputs or the trajectory reproducible, and it does not close #1053. What is still open is filed rather than folded in:

• MonteCarlo draws the flight dictionary three times, so the logged inputs are not the ones flown #1090 the flight dictionary is drawn three times, so the row written to .inputs.txt is not the row that was flown: logged inclination 85.60, flown 84.46 on the shared fixture. Every draw is deterministic under a fixed seed, so this costs the log as a transcript rather than the seeding. Not fixed here: BUG: sample StochasticFlight inputs once per simulation (#1090) #1126 is open against the same method, adds the shared _sample_flight_inputs the fix wants, and also closes the second draw inside StochasticFlight.create_object that a Monte-Carlo-only fix leaves alone. I had written my own and took it back out rather than have two of them.
• MonteCarlo append=True does not continue the same seeded run #1075 append=True is a resume, but not a repair. A run records the root, the count and which logs it wrote, an append continues that root or is refused, and the whole manifest is validated first, so a checkpoint from another lineage or another release no longer reads as continuable. What is left of MonteCarlo append=True does not continue the same seeded run #1075 is filling holes: the next index still comes from a row count rather than from a plan, so a checkpoint with a gap is refused rather than repaired.
• simulate_convergence() cannot reproduce a study from a seed #1077 simulate_convergence seeds neither its batches nor its bootstrap resampling, so a convergence study is outside this guarantee.

#1091 (Parachute pressure noise on the process-global np.random) was the other one. It is closed: #1134 gave Parachute a per-instance RNG and StochasticParachute now derives its seed from this tree.

#1109 was the third one, and it is fixed on develop now: #1122 made dict_generator walk the inputs a model declares rather than every attribute on it, so a tuple initial_solution is no longer read as a distribution. The issue is still open because a pull request into develop cannot close one. This branch dropped its own narrower guard, which walked the whole instance and skipped the component collections, since walking the declared inputs never saw them.

---

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.

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

Follow-up review round

A closer pass after the first reviews turned up four more fixes, all pushed here:

• StochasticRocket._set_stochastic gave the same seed to the rocket body and to every surface, motor, rail button and parachute, so components that sample the same distribution (a main and a drogue parachute, for instance) drew identical cd_s and lag quantiles. Each component now gets its own child of the run's seed, in a fixed order, so they stay independent and reproducible.
• dict_generator and StochasticRocket._randomize_position sampled list-valued attributes (component positions included) with the stdlib random.choice, which random_seed did not govern. Both now draw the index through the model's seeded generator, via a shared _random_choice helper.
• simulate() set up (and, for append=False, truncated) the output files before the seed was validated, so passing a rejected seed destroyed a previous run's results. The seed is validated first now.
• Corrected the seed helper's docstring about RandomState, and moved it to rocketpy.tools so the stochastic models can share it.

Known limitations and follow-ups

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.

Larger items from the review are better handled on their own, so they are filed separately rather than growing this PR:

• MonteCarlo append=True does not continue the same seeded run #1075: the manifest carries the root, the count and the log names now, and an append checks all of them. What is left is claiming indices from a plan instead of counting on from the end, so a checkpoint with a hole is refused rather than filled.
• simulate_convergence() cannot reproduce a study from a seed #1077: simulate_convergence() has no seed, so a convergence study is not reproducible.

Failure safety and log integrity

A later review round found several paths where a run that went wrong could still be reported as a success. Those are fixed here too, with tests:

• The parent waited for every worker with an unbounded join(), in start order. One worker stuck in a native call held it there while another had already set the error event, so neither the error nor the cleanup after it was reached, and Ctrl-C hung on the same join a second time. The wait is bounded now and returns as soon as the event is set. Shutdown signals the whole fleet before waiting on any of it, with a kill fallback.
• The completeness check accepted a corrupt file: unreadable rows were skipped, rows with no index or an index outside the run were ignored, and JSON true or 1.0 passed for index 1 because both compare equal to it. Every row now has to be an object with a plain non-negative int index, the two files have to agree on the exact set, and an interrupted run may be short but not corrupt.
• Both run paths cleared the current payload after the interruptible call rather than before it. Serial Ctrl-C on the first lap surfaced as an UnboundLocalError over the interrupt, and between laps the handler still held the row that had just been written. In the worker, a claim that failed on a later lap reported the simulation that had just succeeded.
• The normal write path released the mutex in finally whether or not acquire() had returned, so a manager that died during acquire raised a second error over the first.
• The error record kept either the inputs or the traceback, never both.
• n_workers was validated after the logs were opened "w+", so asking for a worker count the run cannot use destroyed the previous results on the way to raising.

Tests

Seed handling is unit tested in tests/unit/simulation/test_monte_carlo_determinism.py: accepted seed types, the SeedSequence copy preserving the full .state, 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 width, and the parallel index claim.

tests/integration/simulation/test_monte_carlo_determinism.py runs the real parallel path, no stub, under fork, spawn and forkserver, comparing serial against parallel(2) and parallel(4) per index. The fixtures are built so the properties can fail: the shared stochastic environment has zero wind at every altitude and zero times any factor is zero, so a compounding baseline cannot show up in it, and a bare StochasticAirBrakes gives every parameter a standard deviation of zero. The assertions check the eccentricities and the air brake are among the compared fields, or stripping object identity could quietly empty the comparison.

tests/unit/simulation/test_monte_carlo_log_integrity.py covers what the run is allowed to call a success and how the fleet comes down when it is not.

The thorough version of that test is marked slow and pull-request CI skips those, so it gated nothing. I said earlier that the weekly run would still cover it; that was wrong. Scheduled Tests triggers on schedule and on push to master, with no pull_request, so it only ever runs the default branch and never sees a test that is still on a PR. There is a small version of it now that is not marked, and because it takes the platform default each job ends up gating the start method it actually runs: spawn on Windows and macOS, forkserver on Python 3.14's POSIX default, fork below that. It costs about six seconds.

That turned out to be worth more than the argument for it. While fixing the shutdown timing below I reused a helper that sets the error event, and every successful parallel run started reporting itself as failed. The new gate caught it within seconds of being written, on a path the slow-marked version would not have run in CI at all.

Later review round

Three more from a closer look at the head, all fixed here:

• The completeness check judged the whole file, so a duplicate, a torn row or a pair that disagrees anywhere in it failed the run. The documented way to resume is to interrupt a run and carry on with append=True, and an interrupted run is exactly what leaves that damage, so a file could be damaged once and never resumed. It judges only the indices this run claimed now, and reports rather than raises on what an earlier run left.
• The parent terminated the fleet the moment a worker reported an error, so a worker part way through a write was cut off. Measured at 0.0 ms against the 5000 ms the interrupt path already gave. The shutdown was producing the torn rows the check then reported. Both paths give the same window now.
• A torn row also makes the two files disagree, and the cross-file check ran first, so the error named the symptom rather than the cause. Reordered.

Checkpoint validation

The resume point came from a line count rather than from the indices on disk, and the completeness check trusted it. A blank line makes the two disagree:

two rows                  -> num_of_loaded_sims = 2
two rows + one blank line -> 3
two rows + two blanks     -> 4      (the file holds only 0 and 1 either way)

So the next run starts at 2, index 1 is never written, and a check scoped to the new range reports success. Appending now reads both logs first and refuses unless what it finds is the run it is being asked to continue: every row readable, no index twice, both files holding the same set, and the indices forming exactly the range below the resume point. Nothing is opened for writing until that passes. A run that was not interrupted is then held to the whole range rather than to its own share of it.

A file numbered from 1 is named rather than reported as an off-by-one, since serial runs used to be numbered that way and the answer is to re-baseline.

A file with a hole in it is refused rather than repaired. Filling holes needs workers to claim from a plan instead of counting on from the end, which is
#1075. Until then, refusing loudly beats resuming in the wrong place quietly. The manifest settles which generation wrote the rows and how many there were; it does not decide which index to run next.

The test that used to empty both logs and then assert a four-simulation result holding only indices 2 and 3 was a success is gone. That was the shape of the bug rather than a guard against it.

Two other ways a previous run could be lost:

• multiprocess is an optional extra and was imported inside the parallel path, which runs after both logs have been emptied. An install without rocketpy[monte-carlo] lost its results on the way to the ImportError.
• The Generator rejection advised rng.bit_generator.seed_seq, which NumPy grew in 1.25 while this package declared numpy>=1.13, so the advice raised AttributeError on versions it claimed to support. The floor moves to 1.17, which default_rng has needed all along.

Also: the custom sampler fixture built a Generator in reset_seed and dropped it while sample() drew from the process-global np.random, so nothing in it answered to a seed and the 128-bit path went untested.

Round after the review above

Everything here came out of a reading of the append protocol, and each one was reproduced before it was fixed.

A seed that is an array crashed the lineage check. numpy.random.SeedSequence takes array_like[ints], so entropy can be an ndarray, and comparing two roots holding one answers with an array rather than a verdict: ValueError: The truth value of an array with more than one element is ambiguous. [1, 2, 3] and (1, 2, 3) are the same seed and compared unequal, which warned about a lineage nobody had left. Roots are compared by a canonical fingerprint of the words they generate; n_children_spawned stays in it because __child_seed counts from there. Only int, None and a different int had been covered.

A warning was not enough, and the lineage moved before the files did. One manifest holds one root, so appending from another left the file with two lineages while the manifest named only the newer half: the provenance it reported was wrong rather than incomplete, and my own test wrote that down as correct. An append now continues the root the logs already hold. Leaving random_seed out means continue, a matching seed is allowed, and a different one is refused before anything opens a file, so a file can no longer be given two lineages at all.

A checkpoint cannot be recognised by its indices. The previous release numbered parallel runs from 0 as well, so a clean one of those passed every structural check while its rows came from per-worker entropy, shared component seeds and a different sampling order. Each run now writes a manifest beside its output log naming the schema version, the sampling scheme, the root, a run_id, a committed_count and the two log filenames. It is staged, fsynced and moved into place, and an append validates the whole document: a non-empty checkpoint with no manifest, one that parses but describes no rebuildable root, or one whose seed_chosen is the string "false" rather than a boolean, are all refused. The manifest outlives the object, so the check survives a fresh interpreter.

Three atomic renames are not one transaction. Staging the empty logs first only covered a staging failure. os.replace is atomic per file, so a failure on the second left the first already emptied, the untouched temporaries behind, and the destination narrowed from 0644 to the 0600 a staged file opens at. Each destination is moved aside before its replacement goes in, anything already installed is put back if a later one fails, and the mode is carried across. Every one of the six renames is driven to fail in the tests, plus a KeyboardInterrupt. It is still not a filesystem transaction and the docstring says so.

A bool is an int. simulate refuses one through _is_whole_number; simulate_convergence checked with a plain isinstance, so batch_size=True ran one simulation. tolerance had the same hole.

A row that stopped said nothing about how. A simulation that raised carried error; one dropped because a peer crashed or the user interrupted carried nothing, so it read as a simulation that simply had no error. Those now carry a status.

An early failure lost the inputs. The row is only assembled once the flight is, so anything raising inside create_object or Flight itself left a traceback with nothing about the inputs behind it. dict_generator builds a local dictionary and binds last_rnd_dict once, after its whole loop, so recovery is per model rather than per field: a model that finished publishes, one that raised part way through its own draw does not. Both failure paths recover what was published and mark the row partial. Telling a fresh publication from the previous simulation's is done by keeping a reference to what each model held at seeding, so nothing clears last_rnd_dict, which dict_generator documents as the last generated dictionary.

The manifest recorded a target, not an outcome. It moved on the number of simulations asked for, so a Ctrl-C before the first new row still moved it while the logs stayed where they were, and a target that was never reached claimed rows that are not there. A run opens a generation once its logs exist, so a fresh set belongs to its root even if the first row never lands, and an append stays in the generation it continues. The count is taken from the rows themselves, after the completeness check, and best effort: the run has finished and the rows are on disk, so a count that cannot be taken leaves the previous one for the next append to refuse on.

Breaking change

• Yes

The exact numbers a run produces change (per-index seeding, the env/rocket/flight split and the per-component split within a rocket, the 128-bit int seeds, the serial index now counting from 0 to match parallel, and list-valued attributes and positions 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.

Partially addresses #1053. The per-index seeding for serial and parallel runs is here. The log not matching the flight is #1090 and belongs to #1126, simulate_convergence is #1077, and what is left of #1075 is claiming indices from a plan. I would rather leave #1053 open until those land than close it on a guarantee that does not cover them.

What this PR does and does not promise

The guarantee here is over the sampled inputs: for a given root seed, simulation index i draws the same stochastic parameters whether the run is serial or parallel and however many workers it uses. That is what .inputs.txt records and what the tests compare.

It is deliberately not a guarantee about the whole trajectory yet. The gaps are filed rather than grown into this PR:

• MonteCarlo draws the flight dictionary three times, so the logged inputs are not the ones flown #1090: the flight dictionary is drawn once per argument rather than once per simulation, so the row written to .inputs.txt is not the row that was flown. BUG: sample StochasticFlight inputs once per simulation (#1090) #1126 is open against the same method and fixes it, including the second draw inside StochasticFlight.create_object that a Monte-Carlo-only fix leaves alone. I wrote my own, found theirs, and took mine back out.
• MonteCarlo append=True does not continue the same seeded run #1075: the root, the count and the log identity are persisted and checked. What is left is claiming indices from a plan rather than counting on from the end.
• simulate_convergence() cannot reproduce a study from a seed #1077: simulate_convergence passes no seed to its batches and no random_state to the bootstrap, so a study is not reproducible.

#1091, Parachute pressure noise drawn from the process-global np.random, was
open when this PR was written. It is closed now: #1134 gave Parachute a per-instance RNG and StochasticParachute feeds it a seed derived from this tree, so parachute deployment is inside the guarantee.

Until those land, two runs agreeing on .inputs.txt does not prove they flew the same thing. Once they do, the promise can be restated in terms of results.

[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 96.40719% with 18 lines in your changes missing coverage. Please review.
✅ Project coverage is 85.76%. Comparing base (6c15649) to head (cf99e87).

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

Additional details and impacted files

@@             Coverage Diff             @@
##           develop    #1054      +/-   ##
===========================================
+ Coverage    84.33%   85.76%   +1.43%     
===========================================
  Files          130      130              
  Lines        17258    17673     +415     
===========================================
+ Hits         14554    15158     +604     
+ Misses        2704     2515     -189     

☔ 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?

[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.

[thc1006]

@MateusStano @phmbressan a follow-up pass turned up a few more things worth fixing, so I have pushed them and would appreciate another look when you have time.

Since your reviews:

• Component seeds: StochasticRocket._set_stochastic handed the same seed to the body and to every surface, motor, rail button and parachute, so a main and a drogue parachute drew the same cd_s and lag quantiles. Each component now gets its own child of the run's seed.
• List sampling: dict_generator and _randomize_position sampled list-valued attributes (component positions among them) with the stdlib random.choice, which random_seed did not control. Both now draw through the model's seeded generator.
• File safety: simulate() truncated the output files before the seed was validated, so a rejected seed destroyed a previous run's results. Validation runs first now.
• Docs: corrected the seed helper's note about RandomState, since a 128-bit int does not fit its 32-bit seed.

I also marked the earlier threads resolved. The race and the SeedSequence copy are both fixed in the current code, and the dangling-files question checked out: the run writes only under tmp_path.

A few larger items from the same review are better as their own issues, so I opened #1075 (append continuation), #1076 (a full parallel test under spawn and forkserver) and #1077 (a seed for simulate_convergence), and linked them from the description. Thanks for the careful reviews.

[thc1006]

A quick note on the red CI here, so it isn't mistaken for a regression from this change: the failing jobs crash in test_flight_animation_export_gif (a VTK/PyVista off-screen GIF export) with Fatal Python error: Bus error. The same crash, at the same file and line, also hit develop's own Tests run a day ago (https://github.com/RocketPy-Team/RocketPy/actions/runs/29697106388), so it looks like a pre-existing flake in the off-screen rendering tests rather than anything this PR introduced. I checked the dependency set too: the green and red runs installed identical vtk, pyvista, matplotlib and pillow versions, and nothing in this branch touches the plotting or animation code.

Re-running usually clears it. Happy to help look at the flaky animation tests on their own if that would be useful.

Filed #1078 to track the flaky animation tests.