ENH: reproducible Monte Carlo via per-simulation-index seeding#1054
ENH: reproducible Monte Carlo via per-simulation-index seeding#1054thc1006 wants to merge 8 commits into
Conversation
Codecov Report❌ Patch coverage is
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. 🚀 New features to boost your workflow:
|
|
The seeding logic is unit-tested in The lines codecov still shows uncovered are all in the parallel path: |
0d37ed6 to
761c092
Compare
phmbressan
left a comment
There was a problem hiding this comment.
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?
Addresses review feedback on RocketPy-Team#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>
Addresses review feedback on RocketPy-Team#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>
9c020b6 to
3e22729
Compare
Addresses review feedback on RocketPy-Team#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>
3e22729 to
6bf8bb6
Compare
|
@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 |
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>
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>
Addresses review feedback on RocketPy-Team#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>
6bf8bb6 to
c529d0a
Compare
|
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:
Both earlier concerns are still handled: the parallel claim holds the mutex across the check and the increment, and a supplied |
…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>
…rator 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>
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>
c529d0a to
e5ba940
Compare
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>
…lper 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>
|
@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:
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 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 |
5a9c119 to
da2ba5c
Compare
|
A quick note on the red CI here, so it isn't mistaken for a regression from this change: the failing jobs crash in 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. |
Pull request type
Current behavior
MonteCarlo.simulate()seeds the stochastic models per worker in parallel mode (from a fresh, unseedednp.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_seedtosimulate(). From that root, simulation indexiis seeded from its own child of the root seed, derived before simulationiruns, so indeximaps 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:
iis built by extending the root'sspawn_key, which is exactly howSeedSequence.spawnderives it, sochild(i)is bit-identical toroot.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.int, not aSeedSequence. An int is the seed typenumpy.random.default_rngand the stdlibrandom.Randomboth accept (aSeedSequenceraisesTypeErrorinrandom.Randomsince Python 3.11), so a custom sampler whosereset_seeddocuments an int keeps working. All fouruint32words are combined by value, so the seed is byte-order independent and keeps the full 128-bit pool rather than collapsing to 32 bits.StochasticModel.dict_generatordrew list attributes with the stdlibrandom.choice(an unseeded global instance), sorandom_seeddid not govern them. It now draws the index from the model's own seeded generator, which also avoidsnumpy.random.choicecoercing a heterogeneous list (Function, paths, arrays) to a single dtype.random_seedis a seed, not a live RNG: it takes an int, a numpy integer, a sequence of ints, or aSeedSequence, withNone= fresh entropy so existing behavior is unchanged unless you pass a seed. A suppliedSeedSequenceis 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 statefulrng: aGenerator/BitGeneratoris not accepted, because reducing it to its underlyingSeedSequencewould ignore how far it has been consumed. Passrng.bit_generator.seed_seqto 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.pyunit-tests the seed handling: the accepted seed types; theSeedSequencecopy preserving the full.state(an entropy-only copy fails the test); the O(1) child equal tospawnbit-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.pycovers end-to-end reproducibility (serial, and serial == parallel(2) == parallel(4)) with a stubbedFlight, plus a start-method-invariance test that runs under fork, spawn and forkserver in ordinary CI.tests/unit/stochastic/test_stochastic_model.pycovers the list-sampling fix.Notes from review
keep_simulating()+increment(). Near the end of a run two workers could both passcount < nand 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.SeedSequencewas returned as-is, andspawn()advances its child counter, so passing the same object twice was not reproducible. It is now copied from its full state, andGenerator/BitGeneratorare 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_stochasticgave 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 identicalcd_sandlagquantiles. Each component now gets its own child of the run's seed, in a fixed order, so they stay independent and reproducible.dict_generatorandStochasticRocket._randomize_positionsampled list-valued attributes (component positions included) with the stdlibrandom.choice, whichrandom_seeddid not govern. Both now draw the index through the model's seeded generator, via a shared_random_choicehelper.simulate()set up (and, forappend=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.RandomState, and moved it torocketpy.toolsso 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 moderndefault_rng(or the stdlibrandom.Random) takes it fine; one built onRandomStatewould need to reduce it. SinceRandomStateis 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.Three larger items from the review are better handled on their own, so I filed them separately rather than growing this PR:
append=Truedoes not continue the same seeded run #1075:append=Truedoes not yet continue the same seeded run (the root is not persisted, and the next index comes from a row count).simulate_convergence()cannot reproduce a study from a seed #1077:simulate_convergence()has no seed, so a convergence study is not reproducible.Breaking change
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_simulatechecks apogee and impact velocity within a tolerance and still passes), andrandom_seedis opt-in.Closes #1053