feat: reinit! seed kwarg (#254)

* feat: reinit! seed kwarg

* format

Author: oameye

CHANGELOG.md

@@ -1,5 +1,7 @@
-# main
+# v3.16
 
+- `reinit!(ds::CoupledSDEs)` now accepts a `seed::UInt64` keyword that re-seeds the noise process random number generator. The default is `rand(UInt64)`, so each `reinit!` produces a fresh, independent noise realization unless an explicit seed is given.
+- The `CoupledSDEs` constructor's default `seed` is now `rand(UInt64)` instead of `UInt64(0)`. Behavior is unchanged (the SDE problem already used a random seed when given `0`), but `prob.seed` now stores the actual seed used, so a realization can be reproduced by reading it back and passing it to `reinit!`.
 - Support StochasticDiffEq v6.96+, where the diffusion function is no longer reachable via `integrator.g`. The extension now reads it from the `SDEProblem` (`prob.g`), which is a stable SciMLBase accessor that works across both old and new StochasticDiffEq versions.
 - Run the multiplicative and non-diagonal noise examples on the `CoupledSDEs` docs page so they actually integrate a trajectory.
 

Project.toml

@@ -1,12 +1,13 @@
 name = "DynamicalSystemsBase"
 uuid = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
 repo = "https://github.com/JuliaDynamics/DynamicalSystemsBase.jl.git"
-version = "3.16.1"
+version = "3.16.0"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"

ext/StochasticSystemsBase.jl

@@ -2,13 +2,14 @@ module StochasticSystemsBase
 
 using DynamicalSystemsBase: DynamicalSystemsBase, SciMLBase, correct_state, CoupledODEs,
     CoupledSDEs, StateSpaceSets, isinplace, _delete, set_parameter!,
-    set_state!, dynamic_rule, isdeterministic, current_state,
+    set_parameters!, set_state!, dynamic_rule, isdeterministic, current_state,
     DynamicalSystemsBase, _set_parameter!, u_modified!,
     additional_details, referrenced_sciml_prob, DEFAULT_DIFFEQ,
-    SVector, SMatrix, current_parameters
+    SVector, SMatrix, current_parameters, initial_state, initial_time
 using SciMLBase: SDEProblem, AbstractSDEIntegrator, __init, SDEFunction, step!
-using StochasticDiffEq: SOSRA
+using StochasticDiffEq: SOSRA, setup_next_step!
 using LinearAlgebra
+import Random
 
 include("src/CoupledSDEs.jl")
 include("src/classification.jl")

ext/src/CoupledSDEs.jl

@@ -35,7 +35,7 @@ function DynamicalSystemsBase.CoupledSDEs(
         diffeq = DEFAULT_STOCH_DIFFEQ,
         noise_prototype = nothing,
         noise_process = nothing,
-        seed = UInt64(0)
+        seed = rand(UInt64)
     )
     IIP = isinplace(f, 4) # from SciMLBase
     if !isnothing(g)
@@ -108,6 +108,32 @@ end
 # This preserves the referrenced MTK system and the originally passed diffeq kwargs
 CoupledSDEs(ds::CoupledSDEs, diffeq) = CoupledSDEs(SDEProblem(ds), merge(ds.diffeq, diffeq))
 
+"""
+    reinit!(ds::CoupledSDEs, u = initial_state(ds); kwargs...) → ds
+
+Re-initialize a [`CoupledSDEs`](@ref). In addition to the keywords accepted by
+`reinit!` for any [`DynamicalSystem`](@ref), the following is supported:
+
+- `seed::UInt64`: re-seed the noise process random number generator. Defaults to
+  `rand(UInt64)`, so by default every `reinit!` produces a fresh, independent noise
+  realization. Pass an explicit `seed` to obtain a reproducible noise stream.
+"""
+function SciMLBase.reinit!(
+        ds::CoupledSDEs, u::AbstractArray = DynamicalSystemsBase.initial_state(ds);
+        p = DynamicalSystemsBase.current_parameters(ds),
+        t0 = DynamicalSystemsBase.initial_time(ds),
+        seed::UInt64 = rand(UInt64), kw...
+    )
+    DynamicalSystemsBase.set_parameters!(ds, p)
+    # Adaptive solvers carry the previous run's last `dt`; resetting it ensures
+    # `seed` reproducibility. Non-adaptive solvers must keep their user-specified `dt`.
+    reset_dt = ds.integ.opts.adaptive
+    SciMLBase.reinit!(ds.integ, u; reinit_dae = false, reset_dt, t0, kw...)
+    Random.seed!(ds.integ.W.rng, seed)
+    setup_next_step!(ds.integ.W, ds.integ.u, ds.integ.p)
+    return ds
+end
+
 """
     CoupledSDEs(ds::CoupledODEs, p = current_parameters(ds); kwargs...)
 
@@ -124,7 +150,7 @@ function DynamicalSystemsBase.CoupledSDEs(
         diffeq = DEFAULT_STOCH_DIFFEQ,
         noise_prototype = nothing,
         noise_process = nothing,
-        seed = UInt64(0)
+        seed = rand(UInt64)
     )
     prob = referrenced_sciml_prob(ds)
     # we want the symbolic jacobian to be transfered over

src/core_systems/continuous_time_sde.jl

@@ -76,7 +76,7 @@ found in the [online docs for `CoupledSDEs`](@ref defining-stochastic-dynamics).
 - `noise_process`: stochastic process as provided by [DiffEqNoiseProcess.jl](https://docs.sciml.ai/DiffEqNoiseProcess/stable/) (default `nothing`, i.e. standard Wiener process)
 - `t0`: initial time (default `0.0`)
 - `diffeq`: DiffEq solver settings (see below)
-- `seed`: random number seed (default `UInt64(0)`)
+- `seed`: random number seed (default `rand(UInt64)`, so each `CoupledSDEs` produces a different noise realization unless a seed is given explicitly)
 
 ## DifferentialEquations.jl interfacing
 
@@ -112,12 +112,3 @@ struct CoupledSDEs{IIP, D, I, P} <: ContinuousTimeDynamicalSystem
     diffeq # isn't parameterized because it is only used for display
     noise_type::NamedTuple
 end
-
-function SciMLBase.reinit!(
-        ds::CoupledSDEs, u::AbstractArray = initial_state(ds);
-        p = current_parameters(ds), t0 = initial_time(ds), kw...
-    )
-    set_parameters!(ds, p)
-    reinit!(ds.integ, u; t0, kw...)
-    return ds
-end

test/stochastic.jl

@@ -168,6 +168,47 @@ end
     end
 end
 
+@testset "seed" begin
+    f(u, p, t) = -u
+
+    @testset "constructor default is randomized" begin
+        # Two CoupledSDEs without explicit seed get different seeds and different trajectories.
+        ds1 = CoupledSDEs(f, [1.0])
+        ds2 = CoupledSDEs(f, [1.0])
+        @test ds1.integ.sol.prob.seed != 0
+        @test ds1.integ.sol.prob.seed != ds2.integ.sol.prob.seed
+        step!(ds1, 1.0); step!(ds2, 1.0)
+        @test current_state(ds1) != current_state(ds2)
+    end
+
+    @testset "constructor explicit seed is reproducible" begin
+        seed = UInt64(20250502)
+        ds1 = CoupledSDEs(f, [1.0]; seed = seed)
+        ds2 = CoupledSDEs(f, [1.0]; seed = seed)
+        step!(ds1, 1.0); step!(ds2, 1.0)
+        @test current_state(ds1) ≈ current_state(ds2)
+    end
+
+    @testset "reinit! default reseeds with fresh randomness" begin
+        ds = CoupledSDEs(f, [1.0]; seed = UInt64(1))
+        reinit!(ds); step!(ds, 1.0); ua = copy(current_state(ds))
+        reinit!(ds); step!(ds, 1.0); ub = copy(current_state(ds))
+        # Different default seeds → different trajectories
+        @test ua != ub
+    end
+
+    @testset "reinit! accepts explicit seed kwarg" begin
+        ds = CoupledSDEs(f, [1.0])
+        # Same explicit seed → same noise stream after reinit
+        reinit!(ds; seed = UInt64(42)); step!(ds, 1.0); ua = copy(current_state(ds))
+        reinit!(ds; seed = UInt64(42)); step!(ds, 1.0); ub = copy(current_state(ds))
+        @test ua ≈ ub
+        # Different explicit seed → different trajectory
+        reinit!(ds; seed = UInt64(43)); step!(ds, 1.0); uc = copy(current_state(ds))
+        @test ua != uc
+    end
+end
+
 # Regression test for https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/issues/251:
 # the auto-generated diffusion closure used to recompute its (constant) output on every
 # call, allocating ~1 KB per `step!` and dominating long integrations. The closure must