Informative error on nonexisting symbolic parameter (#232)

* Informative error on nonexisting symbolic parameter

* mention init cond as parameters

* correct typo

* add a new function for named current parameters

* organize MTK tests better

* add reinit_dae = false

See SciML/ModelingToolkit.jl#3451

* remove `named_parametersd`

Needs MTK to be defined

* allow kw propagation in reinit

* fix reinit for SDE

* simplify ODE test dependency

* continuous tests pass

* all tests pass?

* don't forget array for SDEs

* abvoid method overwrites in tests

* remove RODAS usage

* fix jacobian (?)

* fix jacian?

* add back the type check

only way to ensure the symbolic jacobian was extracted.

* correct version

Author: Datseris

CHANGELOG.md

@@ -1,3 +1,8 @@
+# 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
 
 - `jacobian` will now re-use symbolically generated Jacobian functions if the dynamical system is made through a ModelingToolkit.jl-generated `DEProblem`.

Project.toml

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

src/core/dynamicalsystem_interface.jl

@@ -56,14 +56,12 @@ The referenced MTK model corresponding to the dynamical system can be obtained w
 
 See also the DynamicalSystems.jl tutorial online for an example.
 
-!!! warn "ModelingToolkit.jl v9"
-    In ModelingToolkit.jl v9 the default `split` behavior of the parameter container
-    is `true`. This means that the parameter container is no longer a `Vector{Float64}`
-    by default, which means that you cannot use integers to access parameters.
-    It is recommended to keep `split = true` (default) and only access
-    parameters via their symbolic parameter binding.
-    Use `structural_simplify(sys; split = false)` to allow accessing parameters
-    with integers again.
+!!! warn "Initial conditions and parameters"
+    ModelingToolkit.jl treats initial conditions of all variables as additional parameters.
+    This is because its terminology is aimed primarily at initial value problems rather
+    than dynamical systems. It is strongly recommended when making a dynamical system from
+    an MTK model to only access parameters via symbols, and to not use the `split = false`
+    keyword when creating the problem type from the MTK model.
 
 ## API
 
@@ -79,8 +77,8 @@ unless when developing new algorithm implementations that use dynamical systems.
 ### API - obtain information
 
 - `ds(t)` with `ds` an instance of `DynamicalSystem`: return the state of `ds` at time `t`.
-  For continuous time systems this interpolates and extrapolates,
-  while for discrete time systems it only works if `t` is the current time.
+  For continuous time systems this interpolates current and previous time step and is very accurate;
+  for discrete time systems it only works if `t` is the current time.
 - [`current_state`](@ref)
 - [`initial_state`](@ref)
 - [`observe_state`](@ref)
@@ -194,7 +192,7 @@ Return the state `u` of `ds` _observed_ at "index" `i`. Possibilities are:
   symbolic variables.
 
 For [`ProjectedDynamicalSystem`](@ref), this function assumes that the
-state of the system is the full state space state, not the projected one
+state of the system is the original state space state, not the projected one
 (this makes the most sense for allowing MTK-based indexing).
 
 Use [`state_name`](@ref) for an accompanying name.
@@ -247,7 +245,10 @@ function current_parameter(ds::DynamicalSystem, index, p = current_parameters(ds
     prob = referrenced_sciml_prob(ds)
     if !has_referrenced_model(prob)
         return _get_parameter(p, index)
-    else # symbolic dispatch
+    else # symbolic parameter
+        if !SymbolicIndexingInterface.is_parameter(prob, index)
+            error("Symbolic parameter with name $(index) does not exist in the system.")
+        end
         i = SymbolicIndexingInterface.getp(prob, index)
         return i(p)
     end
@@ -431,7 +432,10 @@ function _set_parameter!(ds::DynamicalSystem, index, value, p = current_paramete
         else
             setproperty!(p, index, value)
         end
-    else
+    else # symbolic parameter
+        if !SymbolicIndexingInterface.is_parameter(prob, index)
+            error("Symbolic parameter with name $(index) does not exist in the system.")
+        end
         set! = SymbolicIndexingInterface.setp(prob, index)
         set!(p, value)
     end

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; reset_dt = true, t0)
+    reinit!(ds.integ, u; reinit_dae = false, reset_dt = true, t0, kw...)
     return ds
 end
 

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

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

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"

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

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