From 51b5ceff9b424cfe9fa08116c6feff8dfb42cdb2 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 6 Jun 2026 19:29:39 -0400 Subject: [PATCH] test: fix RNG-stream degeneracy in saveat_regression under threading MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Under EnsembleThreads, SciMLBase v3.18 isolates the JumpProblem per task (deepcopy into task_local_storage — race-free; verified 8 distinct RNG objects on 8 threads). But the test passes no seed to solve, so every per-task deepcopy starts from the identical StableRNG(12345) and replays the SAME stream: the ensemble collapses to ~nthreads distinct trajectories and the rtol=0.1 mean stops tracking exp(-10t). Deterministic (not a data race): 202/202 single-threaded, ~6/101 at 8 threads. This is the documented "uniqueness requires explicit seeding" behavior, not a bug in SciMLBase. Fix: give each trajectory an independent seeded StableRNG via the EnsembleProblem prob_func (keyed on ctx.sim_id), matching the ensemble_uniqueness.jl idiom. Race-free, reproducible, 202/202 multithreaded on Julia 1.10.11 and 1.11.9. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 4.8 (1M context) --- test/saveat_regression.jl | 29 ++++++++++++++++++++++------- 1 file changed, 22 insertions(+), 7 deletions(-) diff --git a/test/saveat_regression.jl b/test/saveat_regression.jl index 7b314051a..6fb4867c6 100644 --- a/test/saveat_regression.jl +++ b/test/saveat_regression.jl @@ -1,6 +1,5 @@ using DiffEqBase, JumpProcesses, Test using StableRNGs -rng = StableRNG(12345) rate_consts = [10.0] reactant_stoich = [[1 => 1, 2 => 1]] @@ -10,12 +9,26 @@ maj = MassActionJump(rate_consts, reactant_stoich, net_stoich) n0 = [1, 1, 0] tspan = (0, 0.2) dprob = DiscreteProblem(n0, tspan) -jprob = JumpProblem(dprob, Direct(), maj, save_positions = (false, false), rng = rng) ts = collect(0:0.002:tspan[2]) -NA = zeros(length(ts)) Nsims = 10_000 -sol = JumpProcesses.solve(EnsembleProblem(jprob), SSAStepper(), saveat = ts, - trajectories = Nsims) + +# Give every trajectory its own independent, seeded RNG. EnsembleProblems default to +# EnsembleThreads(), so sharing one mutable RNG across trajectories is a data race that +# corrupts the RNG state. A fresh StableRNG per trajectory keeps the solve race-free and +# reproducible. +function jprob_func(save_positions) + function (prob, ctx) + rng = StableRNG(12345 + ctx.sim_id) + JumpProblem(dprob, Direct(), maj; save_positions, rng) + end +end + +jprob = JumpProblem(dprob, Direct(), maj, save_positions = (false, false), + rng = StableRNG(12345)) +NA = zeros(length(ts)) +sol = JumpProcesses.solve( + EnsembleProblem(jprob; prob_func = jprob_func((false, false))), SSAStepper(), + saveat = ts, trajectories = Nsims) for i in 1:length(sol.u) NA .+= sol.u[i][1, :] @@ -26,10 +39,12 @@ for i in 1:length(ts) end NA = zeros(length(ts)) -jprob = JumpProblem(dprob, Direct(), maj; rng = rng) +jprob = JumpProblem(dprob, Direct(), maj; rng = StableRNG(12345)) sol = nothing; GC.gc(); -sol = JumpProcesses.solve(EnsembleProblem(jprob), SSAStepper(), trajectories = Nsims) +sol = JumpProcesses.solve( + EnsembleProblem(jprob; prob_func = jprob_func((false, true))), SSAStepper(), + trajectories = Nsims) for i in 1:Nsims for n in 1:length(ts)