Merge branch 'main' into sde-seed

Author: oameye

CHANGELOG.md

@@ -1,7 +1,9 @@
-# v3.17
+# 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.
 
 # v3.15
 

Project.toml

@@ -1,7 +1,7 @@
 name = "DynamicalSystemsBase"
 uuid = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
 repo = "https://github.com/JuliaDynamics/DynamicalSystemsBase.jl.git"
-version = "3.17.0"
+version = "3.16.0"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"

docs/src/CoupledSDEs.md

@@ -84,6 +84,8 @@ function g!(du, u, p, t)
     return nothing
 end
 sde = CoupledSDEs(f!, rand(2)./10; g=g!)
+tr = trajectory(sde, 1.0)
+plot_trajectory(tr...)
 ```
 
 #### Non-diagonal noise
@@ -105,6 +107,8 @@ function g!(du, u, p, t)
 end
 diffeq = (alg = RKMilCommute(), reltol = 1e-3, abstol = 1e-3, dt=0.1)
 sde = CoupledSDEs(f!, rand(2)./10; g=g!, noise_prototype = zeros(2, 2), diffeq = diffeq)
+tr = trajectory(sde, 1.0)
+plot_trajectory(tr...)
 ```
 
 !!! warning

ext/src/CoupledSDEs.jl

@@ -85,6 +85,10 @@ function DynamicalSystemsBase.CoupledSDEs(
     end
 
     solver, remaining = _decompose_into_sde_solver_and_remaining(diffeq)
+    # The default `dtmin` from SciML scales with `tspan`. With our open-ended
+    # `tspan = (0, 1e11)` it becomes ~1e-5, which is too coarse for the SDE
+    # adaptive controller and causes spurious `DtLessThanMin` aborts.
+    remaining = haskey(remaining, :dtmin) ? remaining : merge((dtmin = 0.0,), remaining)
     integ = __init(
         prob,
         solver;
@@ -280,21 +284,29 @@ function construct_diffusion_function(
         A = sqrt(cov)
         if IIP
             if isdiag(cov)
-                g = (du, u, p, t) -> du .= noise_strength .* diag(A)
+                diag_const = collect(noise_strength .* diag(A))
+                g = let diag_const = diag_const
+                    (du, u, p, t) -> (du .= diag_const; nothing)
+                end
             else
-                g = (du, u, p, t) -> du .= noise_strength .* A
+                A_const = collect(noise_strength .* A)
+                g = let A_const = A_const
+                    (du, u, p, t) -> (du .= A_const; nothing)
+                end
                 noise_prototype = zeros(size(A))
                 # ^ we could make this sparse to make it more performant
             end
         else
             if isdiag(cov)
-                g = (u, p, t) -> SVector{length(diag(A)), eltype(A)}(
-                    diag(noise_strength .* A)
-                )
+                diag_const = SVector{D, eltype(A)}(diag(noise_strength .* A))
+                g = let diag_const = diag_const
+                    (u, p, t) -> diag_const
+                end
             else
-                g = (u, p, t) -> SMatrix{size(A)..., eltype(A)}(
-                    noise_strength .* A
-                )
+                A_const = SMatrix{size(A)..., eltype(A)}(noise_strength .* A)
+                g = let A_const = A_const
+                    (u, p, t) -> A_const
+                end
                 noise_prototype = zeros(size(A))
             end
         end

ext/src/classification.jl

@@ -50,7 +50,8 @@ function diffusion_function(g, IIP, noise_prototype)
 end
 
 function diffusion_function(ds::CoupledSDEs{IIP}) where {IIP}
-    return diffusion_function(ds.integ.g, IIP, referrenced_sciml_prob(ds).noise_rate_prototype)
+    prob = referrenced_sciml_prob(ds)
+    return diffusion_function(prob.g, IIP, prob.noise_rate_prototype)
 end
 
 """

src/core_systems/continuous_time_ode.jl

@@ -52,7 +52,7 @@ to access the solvers. The default `diffeq` is:
 $(DynamicalSystemsBase.DEFAULT_DIFFEQ)
 
 `diffeq` keywords can also include `callback` for [event handling
-](http://docs.juliadiffeq.org/latest/features/callback_functions.html).
+](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/).
 
 The convenience constructors `CoupledODEs(prob::ODEProblem [, diffeq])` and
 `CoupledODEs(ds::CoupledODEs [, diffeq])` are also available.

test/continuous.jl

@@ -41,6 +41,7 @@ end
     prob = lorenz_oop.integ.sol.prob
     ds = CoupledODEs(prob, (alg = Vern9(), abstol = 0.0, reltol = 1.0e-6, verbose = false))
     @test ds.integ.alg isa Vern9
-    @test ds.integ.opts.verbose == false
+    # SciML moved from Bool verbose to a `DEVerbosity` struct of per-toggle verbosities.
+    @test nameof(typeof(ds.integ.opts.verbose.linear_verbosity)) == :None
 
 end

test/jacobian.jl

@@ -48,23 +48,23 @@ end
 
     jac = jacobian(ds)
     @test jac isa GeneratedFunctionWrapper
-    @test jac([1.0, 1.0], [], 0.0) == [-3 0;0 3]
+    @test jac([1.0, 1.0], prob.p, 0.0) == [-3 0;0 3]
 
     @testset "CoupledSDEs" begin
         # just to check if MTK @brownian does not give any problems
         using StochasticDiffEq
-        @brownian β
+        @brownians β
         eqs = [
             D(u[1]) ~ 3.0 * u[1] + β,
             D(u[2]) ~ -3.0 * u[2] + β,
         ]
-        @mtkbuild sys = System(eqs, t)
+        @mtkcompile sys = System(eqs, t)
 
         prob = SDEProblem(sys, [1.0, 1.0], (0.0, 1.0), jac = true)
         sde = CoupledSDEs(prob)
 
         jac = jacobian(sde)
         @test jac isa GeneratedFunctionWrapper
-        @test jac([1.0, 1.0], [], 0.0) == [-3 0; 0 3]
+        @test jac([1.0, 1.0], prob.p, 0.0) == [-3 0; 0 3]
     end
 end

test/stochastic.jl

@@ -73,15 +73,16 @@ end
         diffeq = (alg = SRA(), abstol = 1.0e-3, reltol = 1.0e-3, verbose = false)
     )
     @test lorenz_SRA.integ.alg isa SRA
-    @test lorenz_SRA.integ.opts.verbose == false
+    # SciML moved from Bool verbose to a `DEVerbosity` struct of per-toggle verbosities.
+    @test nameof(typeof(lorenz_SRA.integ.opts.verbose.linear_verbosity)) == :None
 
     # also test SDEproblem creation
     prob = lorenz_SRA.integ.sol.prob
 
     ds = CoupledSDEs(prob, (alg = SRA(), abstol = 0.0, reltol = 1.0e-3, verbose = false))
 
     @test ds.integ.alg isa SRA
-    @test ds.integ.opts.verbose == false
+    @test nameof(typeof(ds.integ.opts.verbose.linear_verbosity)) == :None
 
     @test_throws ArgumentError CoupledSDEs(prob; diffeq = (alg = SRA(),))
 
@@ -104,7 +105,10 @@ end
         corr = CoupledSDEs(f, zeros(2); covariance = [1 0.3; 0.3 1])
         corr_alt = CoupledSDEs(f, zeros(2); g = g, noise_prototype = zeros(2, 2))
         @test corr.noise_type == corr_alt.noise_type
-        @test all(corr.integ.g(zeros(2), (), 0.0) .== corr_alt.integ.g(zeros(2), (), 0.0))
+        @test all(
+            DynamicalSystemsBase.referrenced_sciml_prob(corr).g(zeros(2), (), 0.0) .==
+                DynamicalSystemsBase.referrenced_sciml_prob(corr_alt).g(zeros(2), (), 0.0)
+        )
     end
 
     @testset "ArgumentError" begin
@@ -202,5 +206,73 @@ end
         # Different explicit seed → different trajectory
         reinit!(ds; seed = UInt64(43)); step!(ds, 1.0); uc = copy(current_state(ds))
         @test ua != uc
+# 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
+# now return a precomputed constant with no allocations.
+@testset "auto-diffusion closure is allocation-free (#251)" begin
+    f_oop(u, p, t) = SVector{2}(0.0, 0.0)
+    f_iip(du, u, p, t) = (du .= 0; nothing)
+    Γ = [1.0 0.3; 0.3 1.0]
+
+    function call_oop(g, n)
+        s = SVector(0.1, 0.2)
+        out = SVector(0.0, 0.0)
+        for _ in 1:n
+            out = g(s, nothing, 0.0)
+        end
+        return out
+    end
+    function call_iip!(g, du, n)
+        u = [0.0, 0.0]
+        for _ in 1:n
+            g(du, u, nothing, 0.0)
+        end
+        return du
+    end
+
+    @testset "OOP diagonal" begin
+        ds = CoupledSDEs(f_oop, SVector(0.0, 0.0); noise_strength = 2.5)
+        g = ds.integ.sol.prob.f.g
+        @test g(SVector(0.0, 0.0), nothing, 0.0) === g(SVector(1.0, 1.0), nothing, 5.0)
+        @test g(SVector(0.0, 0.0), nothing, 0.0) == SVector(2.5, 2.5)
+        call_oop(g, 10) # warmup
+        @test (@allocated call_oop(g, 10_000)) < 100
+    end
+
+    @testset "OOP non-diagonal" begin
+        ds = CoupledSDEs(f_oop, SVector(0.0, 0.0); covariance = Γ, noise_strength = 1.5)
+        g = ds.integ.sol.prob.f.g
+        @test g(SVector(0.0, 0.0), nothing, 0.0) === g(SVector(1.0, 1.0), nothing, 5.0)
+        @test g(SVector(0.0, 0.0), nothing, 0.0) ≈ 1.5 .* sqrt(Γ)
+        call_oop(g, 10)
+        @test (@allocated call_oop(g, 10_000)) < 100
+    end
+
+    @testset "IIP diagonal" begin
+        ds = CoupledSDEs(f_iip, [0.0, 0.0]; noise_strength = 2.5)
+        g = ds.integ.sol.prob.f.g
+        du = zeros(2)
+        g(du, [0.0, 0.0], nothing, 0.0)
+        @test du == [2.5, 2.5]
+        call_iip!(g, du, 10)
+        @test (@allocated call_iip!(g, du, 10_000)) < 100
+    end
+
+    @testset "IIP non-diagonal" begin
+        ds = CoupledSDEs(f_iip, [0.0, 0.0]; covariance = Γ, noise_strength = 1.5)
+        g = ds.integ.sol.prob.f.g
+        DU = zeros(2, 2)
+        g(DU, [0.0, 0.0], nothing, 0.0)
+        @test DU ≈ 1.5 .* sqrt(Γ)
+        function call_iip_mat!(g, DU, n)
+            u = [0.0, 0.0]
+            for _ in 1:n
+                g(DU, u, nothing, 0.0)
+            end
+            return DU
+        end
+        call_iip_mat!(g, DU, 10)
+        @test (@allocated call_iip_mat!(g, DU, 10_000)) < 100
     end
 end