test: fix RNG-stream degeneracy in saveat_regression under threading

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 <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Author: ChrisRackauckas

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)