use integrator rng API

.github/workflows/Tests.yml

@@ -29,6 +29,7 @@ jobs:
         group:
           - InterfaceI
           - InterfaceII
+          - ThreadSafety
           - QA
         exclude:
           - version: "pre"

.github/workflows/ThreadSafety.yml

@@ -0,0 +1,65 @@
+name: "Thread Safety Tests"
+
+on:
+  push:
+    branches:
+      - master
+    paths-ignore:
+      - 'docs/**'
+  pull_request:
+    branches:
+      - master
+    paths-ignore:
+      - 'docs/**'
+
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: ${{ github.ref_name != github.event.repository.default_branch || github.ref != 'refs/tags/v*' }}
+
+jobs:
+  thread-safety:
+    name: "Thread Safety"
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - "1"
+          - "lts"
+          - "pre"
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: "Setup Julia ${{ matrix.version }}"
+        uses: julia-actions/setup-julia@v2
+        with:
+          version: "${{ matrix.version }}"
+
+      - name: "Cache Julia packages"
+        uses: julia-actions/cache@v2
+        with:
+          token: "${{ secrets.GITHUB_TOKEN }}"
+
+      - name: "Build package"
+        uses: julia-actions/julia-buildpkg@v1
+
+      - name: "Run thread safety tests (4 threads)"
+        run: |
+          julia --threads=4 --code-coverage=user --check-bounds=yes --compiled-modules=yes \
+                --project=@. --color=yes -e '
+            using Pkg
+            Pkg.test()
+          '
+        env:
+          GROUP: ThreadSafety
+
+      - name: "Process Coverage"
+        uses: julia-actions/julia-processcoverage@v1
+
+      - name: "Report Coverage"
+        uses: codecov/codecov-action@v5
+        if: ${{ github.event.pull_request.head.repo.full_name == github.repository }}
+        with:
+          files: lcov.info
+          token: "${{ secrets.CODECOV_TOKEN }}"
+          fail_ci_if_error: false

HISTORY.md

@@ -1,6 +1,26 @@
 # Breaking updates and feature summaries across releases
 
-## JumpProcesses unreleased (master branch)
+## 10.0 (Breaking)
+
+  - **Breaking**: The `rng` keyword argument has been removed from
+    `JumpProblem`. Pass `rng` to `solve` or `init` instead:
+    ```julia
+    # Before (no longer works):
+    jprob = JumpProblem(dprob, Direct(), jump; rng = Xoshiro(1234))
+    sol = solve(jprob, SSAStepper())
+
+    # After:
+    jprob = JumpProblem(dprob, Direct(), jump)
+    sol = solve(jprob, SSAStepper(); rng = Xoshiro(1234))
+    ```
+  - RNG state is now owned by the integrator, not the aggregator. This
+    eliminates data races when sharing a `JumpProblem` across threads and
+    ensures a single, consistent RNG priority across all solver pathways:
+    `rng` > `seed` > `Random.default_rng()`.
+  - `rng` and `seed` kwargs are fully supported on `solve`/`init` for all
+    solver pathways (SSAStepper, ODE, SDE, tau-leaping).
+  - `SSAIntegrator` now supports the `SciMLBase` RNG interface (`has_rng`,
+    `get_rng`, `set_rng!`).
 
 ## 9.14
 

Project.toml

@@ -45,22 +45,25 @@ KernelAbstractions = "0.9"
 LinearAlgebra = "1"
 LinearSolve = "3"
 OrdinaryDiffEq = "6"
-OrdinaryDiffEqCore = "1.32.0"
+OrdinaryDiffEqCore = "3.11"
 Pkg = "1"
 PoissonRandom = "0.4"
 Random = "1"
 RecursiveArrayTools = "3.35"
 Reexport = "1.2"
 SafeTestsets = "0.1"
-SciMLBase = "2.115"
+SciMLBase = "2.147"
 StableRNGs = "1"
 StaticArrays = "1.9.8"
 Statistics = "1"
-StochasticDiffEq = "6.82"
+StochasticDiffEq = "6.95"
 SymbolicIndexingInterface = "0.3.36"
 Test = "1"
 julia = "1.10"
 
+[sources]
+SciMLBase = {url = "https://github.com/isaacsas/SciMLBase.jl.git", rev = "ensemble_rng_redesign"}
+
 [extras]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"

docs/src/api.md

@@ -40,9 +40,91 @@ RSSACR
 SortingDirect
 ```
 
-# Private API Functions
+## Random Number Generator Control
+
+JumpProcesses supports controlling the random number generator (RNG) used for
+jump sampling via the `rng` and `seed` keyword arguments to `solve` or `init`.
+
+### `rng` keyword argument
+
+Pass any `AbstractRNG` to `solve` or `init`:
+
+```julia
+using Random, StableRNGs
+
+# Using a StableRNG for cross-version reproducibility
+sol = solve(jprob, SSAStepper(); rng = StableRNG(1234))
+
+# Using Julia's built-in Xoshiro
+sol = solve(jprob, Tsit5(); rng = Xoshiro(42))
+```
+
+### `seed` keyword argument
+
+As a shorthand, pass an integer `seed` to create a `Xoshiro` generator:
+
+```julia
+sol = solve(jprob, SSAStepper(); seed = 1234)
+# equivalent to: solve(jprob, SSAStepper(); rng = Xoshiro(1234))
+```
+
+### Resolution priority
+
+When both `rng` and `seed` are passed to the same `solve`/`init` call, `rng`
+takes priority:
+
+| User provides | Result |
+|---|---|
+| `rng` via `solve`/`init` | Uses that `rng` |
+| `seed` via `solve`/`init` | Creates `Xoshiro(seed)` |
+| Nothing | Uses `Random.default_rng()` (SSAStepper, ODE, tau-leaping) or a randomly-seeded `Xoshiro` (SDE) |
+
+### Behavior by solver pathway
+
+| Solver | Default RNG (nothing passed) | `rng` / `seed` support |
+|---|---|---|
+| `SSAStepper` | `Random.default_rng()` | Full support via `solve`/`init` kwargs |
+| ODE solvers (e.g., `Tsit5`) | `Random.default_rng()` | Full support via `solve`/`init` kwargs |
+| SDE solvers (e.g., `SRIW1`) | Randomly-seeded `Xoshiro` | Full support; `TaskLocalRNG` is auto-converted to `Xoshiro` |
+| `SimpleTauLeaping` | `Random.default_rng()` | Full support via `solve` kwargs |
+
+!!! note
+    For reproducible simulations, always pass an explicit `rng` or `seed`.
+    The default RNG is shared global state and may produce different results
+    depending on prior usage.
+
+# Private / Developer API
 
 ```@docs
 ExtendedJumpArray
 SSAIntegrator
 ```
+
+## Internal Dispatch Pathways
+
+The following table documents which code handles `solve`/`init` for each solver
+type. This is relevant for developers working on JumpProcesses or its solver
+backends.
+
+| Solver type | `__solve` handled by | `__init` handled by | Uses `__jump_init`? |
+|---|---|---|---|
+| `SSAStepper` | JumpProcesses (`solve.jl`) | JumpProcesses (`SSA_stepper.jl`) | No |
+| ODE (e.g., `Tsit5`) | JumpProcesses (`solve.jl`) | JumpProcesses (`solve.jl`) → OrdinaryDiffEq | Yes |
+| SDE (e.g., `SRIW1`) | StochasticDiffEq | StochasticDiffEq | No |
+| `SimpleTauLeaping` | JumpProcesses (`simple_regular_solve.jl`, custom `DiffEqBase.solve`) | N/A | No |
+
+For **SSAStepper**, `rng` is resolved via `resolve_rng` in `SSA_stepper.jl`'s
+`__init` and stored on the [`SSAIntegrator`](@ref).
+
+For **ODE solvers**, `rng` is resolved via `resolve_rng` in `__jump_init`
+(`solve.jl`) and forwarded to OrdinaryDiffEq's `init`, which stores it on the
+`ODEIntegrator`.
+
+For **SDE solvers**, StochasticDiffEq handles the full solve/init pathway
+directly (JumpProcesses' ambiguity-fix `__solve` method is never dispatched to).
+StochasticDiffEq has its own `_resolve_rng` that additionally handles
+`TaskLocalRNG` conversion and the problem's stored seed.
+
+For **tau-leaping**, JumpProcesses defines a custom `DiffEqBase.solve` that
+bypasses the standard `__solve`/`__init` pathway. It calls `resolve_rng`
+directly with the `rng` and `seed` kwargs from the `solve` call.

docs/src/applications/advanced_point_process.md

@@ -387,10 +387,10 @@ function Base.rand(rng::AbstractRNG,
     out = Array{History, 1}(undef, n)
     p = params(pp)
     dprob = DiscreteProblem([0], tspan, p)
-    jprob = JumpProblem(dprob, Coevolve(), jumps...; dep_graph = pp.g, save_positions, rng)
+    jprob = JumpProblem(dprob, Coevolve(), jumps...; dep_graph = pp.g, save_positions)
     for i in 1:n
         params!(pp, p)
-        solve(jprob, SSAStepper())
+        solve(jprob, SSAStepper(); rng)
         out[i] = deepcopy(p.h)
     end
     return out

docs/src/faq.md

@@ -55,21 +55,27 @@ jset = JumpSet(; constant_jumps = cjvec, variable_jumps = vjtuple,
 
 ## How can I set the random number generator used in the jump process sampling algorithms (SSAs)?
 
-Random number generators can be passed to `JumpProblem` via the `rng` keyword
+Random number generators can be passed to `solve` or `init` via the `rng` keyword
 argument. Continuing the previous example:
 
 ```julia
-#] add RandomNumbers
-using RandomNumbers
-jprob = JumpProblem(dprob, Direct(), maj,
-    rng = Xorshifts.Xoroshiro128Star(rand(UInt64)))
+using Random
+jprob = JumpProblem(dprob, Direct(), maj)
+sol = solve(jprob, SSAStepper(); rng = Xoshiro(1234))
+```
+
+Any `AbstractRNG` can be used. For example, to use a generator from
+[StableRNGs.jl](https://github.com/JuliaRandom/StableRNGs.jl):
+
+```julia
+using StableRNGs
+sol = solve(jprob, SSAStepper(); rng = StableRNG(1234))
 ```
 
-uses the `Xoroshiro128Star` generator from
-[RandomNumbers.jl](https://github.com/JuliaRandom/RandomNumbers.jl).
+A `seed` keyword argument is also supported as a shorthand for creating a `Xoshiro`
+generator: `solve(jprob, SSAStepper(); seed = 1234)`.
 
-On version 1.7 and up, JumpProcesses uses Julia's built-in random number generator by
-default. On versions below 1.7 it uses `Xoroshiro128Star`.
+By default, JumpProcesses uses Julia's built-in `Random.default_rng()`.
 
 ## What are these aggregators and aggregations in JumpProcesses?
 

docs/src/tutorials/simple_poisson_process.md

@@ -371,11 +371,10 @@ with ``N(t)`` a Poisson counting process with constant transition rate
 ``\lambda``, and the ``C_i`` independent and identical samples from a uniform
 distribution over ``\{-1,1\}``. We can simulate such a process as follows.
 
-We first ensure that we use the same random number generator as JumpProcesses. We
-can either pass one as an input to [`JumpProblem`](@ref) via the `rng` keyword
-argument, and make sure it is the same one we use in our `affect!` function, or
-we can just use the default generator chosen by JumpProcesses if one is not
-specified, `JumpProcesses.DEFAULT_RNG`. Let's do the latter
+We first ensure that we use the same random number generator as JumpProcesses.
+Custom RNGs can be passed to `solve` or `init` via the `rng` keyword argument.
+If no RNG is specified, JumpProcesses uses `Random.default_rng()`, which is also
+available as `JumpProcesses.DEFAULT_RNG`. Let's use the default
 
 ```@example tut1
 rng = JumpProcesses.DEFAULT_RNG

src/JumpProcesses.jl

@@ -40,7 +40,7 @@ using DiffEqBase: DiffEqBase, CallbackSet, ContinuousCallback, DAEFunction,
                   ODESolution, ReturnCode, SDEFunction, SDEProblem, add_tstop!,
                   deleteat!, isinplace, remake, savevalues!, step!,
                   u_modified!
-using SciMLBase: SciMLBase, DEIntegrator
+using SciMLBase: SciMLBase, DEIntegrator, has_rng, get_rng, set_rng!
 
 abstract type AbstractJump end
 abstract type AbstractMassActionJump <: AbstractJump end

src/SSA_stepper.jl

@@ -13,6 +13,9 @@ Highly efficient integrator for pure jump problems that involve only `ConstantRa
     `SSAStepper`.
   - Only supports a limited subset of the output controls from the common solver interface,
     specifically `save_start`, `save_end`, and `saveat`.
+  - Supports `rng` and `seed` keyword arguments in `solve`/`init` to control the random
+    number generator used for jump sampling. `rng` accepts any `AbstractRNG`, while `seed`
+    creates a `Xoshiro` generator. `rng` takes priority over `seed`.
   - As when using jumps with ODEs and SDEs, saving controls for whether to save each time a
     jump occurs are via the `save_positions` keyword argument to `JumpProblem`. Note that when
     choosing `SSAStepper` as the timestepper, `save_positions = (true,true)`, `(true,false)`,
@@ -58,17 +61,18 @@ for details.
 """
 struct SSAStepper <: DiffEqBase.DEAlgorithm end
 SciMLBase.allows_late_binding_tstops(::SSAStepper) = true
+SciMLBase.supports_solve_rng(::JumpProblem, ::SSAStepper) = true
 
 """
 $(TYPEDEF)
 
-Solution objects for pure jump problems solved via `SSAStepper`.
+Integrator for pure jump problems solved via `SSAStepper`.
 
 ## Fields
 
 $(FIELDS)
 """
-mutable struct SSAIntegrator{F, uType, tType, tdirType, P, S, CB, SA, OPT, TS} <:
+mutable struct SSAIntegrator{F, uType, tType, tdirType, P, S, CB, SA, OPT, TS, R} <:
                AbstractSSAIntegrator{SSAStepper, Nothing, uType, tType}
     """The underlying `prob.f` function. Not currently used."""
     f::F
@@ -108,6 +112,23 @@ mutable struct SSAIntegrator{F, uType, tType, tdirType, P, S, CB, SA, OPT, TS} <
     alias_tstops::Bool
     """If true indicates we have already allocated the tstops array"""
     copied_tstops::Bool
+    """The random number generator."""
+    rng::R
+end
+
+SciMLBase.has_rng(::SSAIntegrator) = true
+SciMLBase.get_rng(integrator::SSAIntegrator) = integrator.rng
+function SciMLBase.set_rng!(integrator::SSAIntegrator, rng)
+    R = typeof(integrator.rng)
+    if !isa(rng, R)
+        throw(ArgumentError(
+            "Cannot set RNG of type $(typeof(rng)) on an integrator " *
+            "whose RNG type parameter is $R. " *
+            "Construct a new integrator via `init(prob, alg; rng = your_rng)` instead."
+        ))
+    end
+    integrator.rng = rng
+    nothing
 end
 
 (integrator::SSAIntegrator)(t) = copy(integrator.u)
@@ -198,6 +219,7 @@ function DiffEqBase.__init(jump_prob::JumpProblem,
         save_start = true,
         save_end = true,
         seed = nothing,
+        rng = nothing,
         alias_jump = Threads.threadid() == 1,
         saveat = nothing,
         callback = nothing,
@@ -219,19 +241,13 @@ function DiffEqBase.__init(jump_prob::JumpProblem,
 
     # Check for continuous callbacks passed via kwargs (from JumpProblem constructor or solve)
     check_continuous_callback_error(callback)
+
+    _rng = resolve_rng(rng, seed)
+
     if alias_jump
         cb = jump_prob.jump_callback.discrete_callbacks[end]
-        if seed !== nothing
-            Random.seed!(cb.condition.rng, seed)
-        end
     else
         cb = deepcopy(jump_prob.jump_callback.discrete_callbacks[end])
-        # Only reseed if an explicit seed is provided. This respects the user's RNG choice
-        # and enables reproducibility. For EnsembleProblems, use prob_func to set unique seeds
-        # for each trajectory if different results are needed.
-        if seed !== nothing
-            Random.seed!(cb.condition.rng, seed)
-        end
     end
     opts = (callback = CallbackSet(callback),)
 
@@ -286,7 +302,7 @@ function DiffEqBase.__init(jump_prob::JumpProblem,
 
     integrator = SSAIntegrator(prob.f, copy(prob.u0), prob.tspan[1], prob.tspan[1], tdir,
         prob.p, sol, 1, prob.tspan[1], cb, _saveat, save_everystep,
-        save_end, cur_saveat, opts, _tstops, 1, false, true, alias_tstops, false)
+        save_end, cur_saveat, opts, _tstops, 1, false, true, alias_tstops, false, _rng)
     cb.initialize(cb, integrator.u, prob.tspan[1], integrator)
     DiffEqBase.initialize!(opts.callback, integrator.u, prob.tspan[1], integrator)
     if save_start