diff --git a/CHANGELOG.md b/CHANGELOG.md index e3881573..e3ae8887 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/Project.toml b/Project.toml index e8a57200..e4b70148 100644 --- a/Project.toml +++ b/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" diff --git a/ext/StochasticSystemsBase.jl b/ext/StochasticSystemsBase.jl index 28adc3df..ac30d1ca 100644 --- a/ext/StochasticSystemsBase.jl +++ b/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") diff --git a/ext/src/CoupledSDEs.jl b/ext/src/CoupledSDEs.jl index 4001fbdb..bc0f0bc9 100644 --- a/ext/src/CoupledSDEs.jl +++ b/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 diff --git a/src/core_systems/continuous_time_sde.jl b/src/core_systems/continuous_time_sde.jl index fdac6a6b..152e8882 100644 --- a/src/core_systems/continuous_time_sde.jl +++ b/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 diff --git a/test/stochastic.jl b/test/stochastic.jl index ba4c67b8..3a628024 100644 --- a/test/stochastic.jl +++ b/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