Merge branch 'main' into named_param_error

Datseris · web-flow · commit 83a8da3ec092 · 2025-09-09T15:10:03.000+01:00
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -29,7 +29,7 @@ jobs:
           access_token: ${{ github.token }}
 
       # Do tests
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}
diff --git a/.github/workflows/doccleanup.yml b/.github/workflows/doccleanup.yml
@@ -9,7 +9,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: Checkout gh-pages branch
-        uses: actions/checkout@v4
+        uses: actions/checkout@v5
         with:
           ref: gh-pages
 
diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml
@@ -13,7 +13,7 @@ jobs:
       contents: write
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v1
         with:
           version: '1'
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -6,6 +6,7 @@
 # v3.14.0
 
 - `set_parameter!` and `current_parameter` will now throw an informative error message explicitly naming the parameter if the user tries to get/set a symbolic parameter that does not exist in the MTK-generated dynamical system.
+- Fixed a plethora of bugs and test failures associated with updates in the SciML ecosystem.
 
 # v3.13.0
 
diff --git a/Project.toml b/Project.toml
@@ -15,8 +15,14 @@ StateSpaceSets = "40b095a5-5852-4c12-98c7-d43bf788e795"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 
+[weakdeps]
+StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
+
+[extensions]
+StochasticSystemsBase = "StochasticDiffEq"
+
 [compat]
-ForwardDiff = "0.10"
+ForwardDiff = "0.10, 1"
 OrdinaryDiffEqTsit5 = "1.1"
 Reexport = "1"
 Roots = "1, 2"
@@ -26,9 +32,3 @@ Statistics = "1"
 StochasticDiffEq = "6.66.0"
 SymbolicIndexingInterface = "0.3.4"
 julia = "1.9"
-
-[weakdeps]
-StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
-
-[extensions]
-StochasticSystemsBase = "StochasticDiffEq"
diff --git a/src/core_systems/continuous_time_ode.jl b/src/core_systems/continuous_time_ode.jl
@@ -142,10 +142,10 @@ set_state!(ds::CoupledODEs, u::AbstractArray) = (set_state!(ds.integ, u); ds)
 SciMLBase.step!(ds::CoupledODEs, args...) = (step!(ds.integ, args...); ds)
 
 function SciMLBase.reinit!(ds::ContinuousTimeDynamicalSystem, u::AbstractArray = initial_state(ds);
-        p = current_parameters(ds), t0 = initial_time(ds)
+        p = current_parameters(ds), t0 = initial_time(ds), kw...
     )
     set_parameters!(ds, p)
-    reinit!(ds.integ, u; reinit_dae = false, reset_dt = true, t0)
+    reinit!(ds.integ, u; reinit_dae = false, reset_dt = true, t0, kw...)
     return ds
 end
 
diff --git a/src/core_systems/continuous_time_sde.jl b/src/core_systems/continuous_time_sde.jl
@@ -13,7 +13,7 @@ where
 ``\\mathbf{u}(t)`` is the state vector at time ``t``, ``\\mathbf{f}`` describes the
 deterministic dynamics, and the noise term
 ``\\mathbf{g}(\\mathbf{u}, p, t) \\text{d}\\mathcal{N}_t`` describes
-the stochastic forcing in terms of a noise function (or *diffusion function*) 
+the stochastic forcing in terms of a noise function (or *diffusion function*)
 ``\\mathbf{g}`` and a noise process ``\\mathcal{N}_t``. The parameters of the functions
 ``\\mathcal{f}`` and ``\\mathcal{g}`` are contained in the vector ``p``.
 
@@ -55,7 +55,7 @@ Commonly, this is a matrix or sparse matrix. If this is not given, it
 defaults to `nothing`, which means the `g` should be interpreted as being diagonal.
 
 The noise process can be specified via `noise_process`. It defaults to a standard Wiener
-process (Gaussian white noise). 
+process (Gaussian white noise).
 For details on defining noise processes, see the docs of [DiffEqNoiseProcess.jl
 ](https://docs.sciml.ai/DiffEqNoiseProcess/stable/). A complete list of the pre-defined
 processes can be found [here](https://docs.sciml.ai/DiffEqNoiseProcess/stable/noise_processes/).
@@ -111,4 +111,12 @@ struct CoupledSDEs{IIP,D,I,P} <: ContinuousTimeDynamicalSystem
     p0::P
     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/src/core_systems/jacobian.jl b/src/core_systems/jacobian.jl
@@ -6,7 +6,7 @@ import ForwardDiff
     jacobian(ds::CoreDynamicalSystem)
 
 Construct the Jacobian rule for the dynamical system `ds`.
-If the system already has a Jacobian rule constructed via ModelingToolkit it returns this,
+If the system already has a Jacobian rule constructed via ModelingToolkit.jl it returns this,
 otherwise it constructs the Jacobian rule with automatic differentiation using module
 [`ForwardDiff`](https://github.com/JuliaDiff/ForwardDiff.jl).
 
@@ -21,6 +21,17 @@ at the state `u`, parameters `p` and time `t` and save the result in `J0`.
 """
 function jacobian(ds::CoreDynamicalSystem{IIP}) where {IIP}
     if ds isa ContinuousTimeDynamicalSystem
+        # TODO: This is the correct API way to obtain the Jacobian,
+        # however it relies on MTK dependency, so we can't do it.
+        # if has_referrenced_model(ds)
+        #     model = referrenced_sciml_model(ds)
+        #     Joop, Jiip = generate_jacobian(model; expression = Val{false})
+        #     if IIP
+        #         jac = Jiip
+        #     else
+        #         jac = Joop
+        #     end
+        # end
         prob = referrenced_sciml_prob(ds)
         if prob.f isa SciMLBase.AbstractDiffEqFunction && !isnothing(prob.f.jac)
             jac = prob.f.jac
diff --git a/test/Project.toml b/test/Project.toml
@@ -3,7 +3,8 @@ DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503"
 DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
+OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/continuous.jl b/test/continuous.jl
@@ -1,6 +1,7 @@
 using DynamicalSystemsBase, Test
 
-using OrdinaryDiffEq: Vern9, ODEProblem, Rodas5, Tsit5
+using OrdinaryDiffEqTsit5: ODEProblem, Tsit5
+using OrdinaryDiffEqVerner: Vern9
 
 include("test_system_function.jl")
 
@@ -25,14 +26,9 @@ p0 = [10, 28, 8/3]
 
 lorenz_oop = CoupledODEs(lorenz_rule, u0, p0)
 lorenz_iip = CoupledODEs(ODEProblem(lorenz_rule_iip, copy(u0), (0.0, Inf), p0))
-lorenz_vern = CoupledODEs(lorenz_rule, u0, p0;
-    diffeq = (alg = Vern9(), abstol = 1e-9, reltol = 1e-9)
-)
 
-for (ds, iip) in zip((lorenz_oop, lorenz_iip, lorenz_vern), (false, true, false))
-
-    name = (ds === lorenz_vern) ? "lorvern" : "lorenz"
-    @testset "$(name) IIP=$(iip)" begin
+for (ds, iip) in zip((lorenz_oop, lorenz_iip), (false, true))
+    @testset "IIP=$(iip)" begin
         @test dynamic_rule(ds) == (iip ? lorenz_rule_iip : lorenz_rule)
         test_dynamical_system(ds, u0, p0; idt = false, iip=iip)
     end
@@ -42,20 +38,9 @@ end
     lorenz_oop = CoupledODEs(lorenz_rule, u0, p0)
     @test lorenz_oop.integ.alg isa Tsit5
 
-    lorenz_vern = CoupledODEs(lorenz_rule, u0, p0;
-        diffeq = (alg = Vern9(), verbose = false, abstol = 1e-9, reltol = 1e-9)
-    )
-    @test lorenz_vern.integ.alg isa Vern9
-    @test lorenz_vern.integ.opts.verbose == false
-
-    # also test ODEproblem creation
-    prob = lorenz_vern.integ.sol.prob
-
-    ds = CoupledODEs(prob, (alg=Rodas5(autodiff=false), abstol=0.0, reltol=1e-6, verbose=false))
-
-    @test ds.integ.alg isa Rodas5
+    prob = lorenz_oop.integ.sol.prob
+    ds = CoupledODEs(prob, (alg=Vern9(), abstol=0.0, reltol=1e-6, verbose=false))
+    @test ds.integ.alg isa Vern9
     @test ds.integ.opts.verbose == false
 
-    @test_throws ArgumentError CoupledODEs(prob; diffeq = (alg=Rodas5(autodiff=false), ))
-
 end
diff --git a/test/jacobian.jl b/test/jacobian.jl
@@ -30,7 +30,7 @@ end
 
 @testset "MTK Jacobian" begin
     using ModelingToolkit
-    using ModelingToolkit: Num, RuntimeGeneratedFunctions.RuntimeGeneratedFunction
+    using ModelingToolkit: Num, GeneratedFunctionWrapper
     using DynamicalSystemsBase: SciMLBase
     @independent_variables t
     @variables u(t)[1:2]
@@ -42,10 +42,10 @@ end
     sys = structural_simplify(sys)
 
     prob = ODEProblem(sys, [1.0, 1.0], (0.0, 1.0); jac=true)
-    ode = CoupledODEs(prob)
+    ds = CoupledODEs(prob)
 
-    jac = jacobian(ode)
-    @test jac.jac_oop isa RuntimeGeneratedFunction
+    jac = jacobian(ds)
+    @test jac isa GeneratedFunctionWrapper
     @test jac([1.0, 1.0], [], 0.0) == [3 0;0 -3]
 
     @testset "CoupledSDEs" begin
@@ -60,7 +60,7 @@ end
         sde = CoupledSDEs(prob)
 
         jac = jacobian(sde)
-        @test jac.jac_oop isa RuntimeGeneratedFunction
-        @test jac([1.0, 1.0], [], 0.0) == [3 0;0 -3]
+        @test jac isa GeneratedFunctionWrapper
+        @test jac([1.0, 1.0], [], 0.0) == [3 0; 0 -3]
     end
 end
diff --git a/test/stochastic.jl b/test/stochastic.jl
@@ -1,5 +1,5 @@
 using DynamicalSystemsBase, Test
-using OrdinaryDiffEq: Tsit5
+using OrdinaryDiffEqTsit5: Tsit5
 using StochasticDiffEq: SDEProblem, SRA, SOSRA, LambaEM, CorrelatedWienerProcess, EM
 
 StochasticSystemsBase = Base.get_extension(DynamicalSystemsBase, :StochasticSystemsBase)
diff --git a/test/stochastic_characterisation.jl b/test/stochastic_characterisation.jl
@@ -25,7 +25,7 @@ f!(du, u, p, t) = du .= 1.01u
     @test CoupledSDEs(sciml_prob(corr)).noise_type == types
 
     # invertible correlated Wiener noise
-    g!(du, u, p, t) = (du .= [1 0.3; 0.3 1]; return nothing)
+    g! = (du, u, p, t) -> (du .= [1 0.3; 0.3 1]; return nothing)
     corr_alt = CoupledSDEs(f!, zeros(2); g = g!, noise_prototype = zeros(2, 2))
     types = corr_alt.noise_type
     @test values(types) == (true, true, true, true)
@@ -38,14 +38,14 @@ f!(du, u, p, t) = du .= 1.01u
     @test CoupledSDEs(sciml_prob(corr_noninv)).noise_type == types
 
     # non-invertible correlated Wiener noise
-    g!(du, u, p, t) = (du .= [1 0.3 1; 0.3 1 1]; return nothing)
+    g! = (du, u, p, t) -> (du .= [1 0.3 1; 0.3 1 1]; return nothing)
     corr_alt = CoupledSDEs(f!, zeros(2); g = g!, noise_prototype = zeros(2, 3))
     types = corr_alt.noise_type
     @test values(types) == (true, true, true, false)
     @test CoupledSDEs(sciml_prob(corr_alt)).noise_type == types
 
     # non-autonomous Wiener noise
-    g!(du, u, p, t) = (du .= 1 / (1 + t); return nothing)
+    g! = (du, u, p, t) -> (du .= 1 / (1 + t); return nothing)
     addit_non_autom = CoupledSDEs(f!, zeros(2); g = g!)
     types = addit_non_autom.noise_type
     @test values(types) == (true, false, true, false)
@@ -54,7 +54,7 @@ end
 
 @testset "Multiplicative noise" begin
     # multiplicative linear Wiener noise
-    g!(du, u, p, t) = (du .= u; return nothing)
+    g! = (du, u, p, t) -> (du .= u; return nothing)
     linear_multipli = CoupledSDEs(f!, rand(2) ./ 10; g = g!)
     types = linear_multipli.noise_type
     @test values(types) == (false, true, true, false)
@@ -67,21 +67,21 @@ end
     @test CoupledSDEs(sciml_prob(lin_multipli_alt)).noise_type == types
 
     # multiplicative nonlinear Wiener noise
-    g!(du, u, p, t) = (du .= u .^ 2; return nothing)
+    g! = (du, u, p, t) -> (du .= u .^ 2; return nothing)
     nonlinear_multiplicative = CoupledSDEs(f!, rand(2); g = g!)
     types = nonlinear_multiplicative.noise_type
     @test values(types) == (false, true, false, false)
     @test CoupledSDEs(sciml_prob(nonlinear_multiplicative)).noise_type == types
 
     # non-autonomous linear multiplicative Wiener noise
-    g!(du, u, p, t) = (du .= u ./ (1 + t); return nothing)
+    g! = (du, u, p, t) -> (du .= u ./ (1 + t); return nothing)
     addit_non_autom = CoupledSDEs(f!, zeros(2); g = g!)
     types = addit_non_autom.noise_type
     @test values(types) == (false, false, true, false)
     @test CoupledSDEs(sciml_prob(addit_non_autom)).noise_type == types
 
     # non-autonomous nonlinear multiplicative Wiener noise
-    g!(du, u, p, t) = (du .= u .^ 2 ./ (1 + t); return nothing)
+    g! = (du, u, p, t) -> (du .= u .^ 2 ./ (1 + t); return nothing)
     addit_non_autom = CoupledSDEs(f!, zeros(2); g = g!)
     types = addit_non_autom.noise_type
     @test values(types) == (false, false, false, false)
diff --git a/test/stroboscopic.jl b/test/stroboscopic.jl
@@ -1,6 +1,6 @@
 using DynamicalSystemsBase, Test
 
-using OrdinaryDiffEq: Vern9
+using OrdinaryDiffEqVerner: Vern9
 
 include("test_system_function.jl")
 
