diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index 14bea532b3..e696b42380 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -4,6 +4,7 @@ on: push: branches: - master + - v10 tags: '*' pull_request: diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml index 6185015c44..0d3052b969 100644 --- a/.github/workflows/FormatCheck.yml +++ b/.github/workflows/FormatCheck.yml @@ -4,6 +4,7 @@ on: push: branches: - 'master' + - v10 tags: '*' pull_request: diff --git a/.github/workflows/Tests.yml b/.github/workflows/Tests.yml index 52c5482970..a62f7f272d 100644 --- a/.github/workflows/Tests.yml +++ b/.github/workflows/Tests.yml @@ -5,6 +5,7 @@ on: branches: - master - 'release-' + - v10 paths-ignore: - 'docs/**' push: diff --git a/NEWS.md b/NEWS.md index 038b1d79f6..d316ac23fb 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,22 @@ +# ModelingToolkit v10 Release Notes + +### Callbacks + +Callback semantics have changed. + + - There is a new `Pre` operator that is used to specify which values are before the callback. + For example, the affect `A ~ A + 1` should now be written as `A ~ Pre(A) + 1`. This is + **required** to be specified - `A ~ A + 1` will now be interpreted as an equation to be + satisfied after the callback (and will thus error since it is unsatisfiable). + + - All parameters that are changed by a callback must be declared as discrete parameters to + the callback constructor, using the `discrete_parameters` keyword argument. + +```julia +event = SymbolicDiscreteCallback( + [t == 1] => [p ~ Pre(p) + 1], discrete_parameters = [p]) +``` + # ModelingToolkit v9 Release Notes ### Upgrade guide diff --git a/Project.toml b/Project.toml index b661dc402d..1afe6ab76d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ModelingToolkit" uuid = "961ee093-0014-501f-94e3-6117800e7a78" authors = ["Yingbo Ma ", "Chris Rackauckas and contributors"] -version = "9.78.0" +version = "10.0.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" @@ -30,6 +30,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e" FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" +ImplicitDiscreteSolve = "3263718b-31ed-49cf-8a0f-35a466e8af96" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5" @@ -42,6 +43,7 @@ NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" @@ -113,6 +115,7 @@ ForwardDiff = "0.10.3" FunctionWrappers = "1.1" FunctionWrappersWrappers = "0.1" Graphs = "1.5.2" +ImplicitDiscreteSolve = "0.1.2" InfiniteOpt = "0.5" InteractiveUtils = "1" JuliaFormatter = "1.0.47, 2" @@ -174,6 +177,7 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" +OptimizationBase = "bca83a33-5cc9-4baa-983d-23429ab6bcbb" OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1" OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" @@ -193,5 +197,10 @@ StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +[sources] +ModelingToolkitStandardLibrary = { url = "https://github.com/SciML/ModelingToolkitStandardLibrary.jl/", rev = "mtk-v10" } +OptimizationBase = { url = "https://github.com/AayushSabharwal/OptimizationBase.jl", rev = "as/mtk-v10" } +OptimizationMOI = { url = "https://github.com/AayushSabharwal/Optimization.jl", subdir = "lib/OptimizationMOI", rev = "as/mtk-v10" } + [targets] -test = ["AmplNLWriter", "BenchmarkTools", "BoundaryValueDiffEqMIRK", "BoundaryValueDiffEqAscher", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "OrdinaryDiffEqDefault", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET", "OrdinaryDiffEqNonlinearSolve", "Logging"] +test = ["AmplNLWriter", "BenchmarkTools", "BoundaryValueDiffEqMIRK", "BoundaryValueDiffEqAscher", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "OrdinaryDiffEqDefault", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET", "OrdinaryDiffEqNonlinearSolve", "Logging", "OptimizationBase"] diff --git a/TestEnv/Project.toml b/TestEnv/Project.toml new file mode 100644 index 0000000000..4c82d62efa --- /dev/null +++ b/TestEnv/Project.toml @@ -0,0 +1,8 @@ +[deps] +ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739" +NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" +SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" +Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" diff --git a/docs/Project.toml b/docs/Project.toml index ca13773492..ae0d89dd3a 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -40,7 +40,7 @@ Documenter = "1" DynamicQuantities = "^0.11.2, 0.12, 1" FMI = "0.14" FMIZoo = "1" -ModelingToolkit = "8.33, 9" +ModelingToolkit = "10" ModelingToolkitStandardLibrary = "2.19" NonlinearSolve = "3, 4" Optim = "1.7" diff --git a/docs/src/basics/Composition.md b/docs/src/basics/Composition.md index 2e5d4be831..d3de71d696 100644 --- a/docs/src/basics/Composition.md +++ b/docs/src/basics/Composition.md @@ -135,16 +135,14 @@ sys.y = u * 1.1 In a hierarchical system, variables of the subsystem get namespaced by the name of the system they are in. This prevents naming clashes, but also enforces that every unknown and parameter is local to the subsystem it is used in. In some cases it might be desirable to have variables and parameters that are shared between subsystems, or even global. This can be accomplished as follows. ```julia -@parameters a b c d e f +@parameters a b c d # a is a local variable b = ParentScope(b) # b is a variable that belongs to one level up in the hierarchy c = ParentScope(ParentScope(c)) # ParentScope can be nested -d = DelayParentScope(d) # skips one level before applying ParentScope -e = DelayParentScope(e, 2) # second argument allows skipping N levels -f = GlobalScope(f) +d = GlobalScope(d) -p = [a, b, c, d, e, f] +p = [a, b, c, d] level0 = ODESystem(Equation[], t, [], p; name = :level0) level1 = ODESystem(Equation[], t, [], []; name = :level1) ∘ level0 @@ -152,25 +150,19 @@ parameters(level1) #level0₊a #b #c -#level0₊d -#level0₊e -#f +#d level2 = ODESystem(Equation[], t, [], []; name = :level2) ∘ level1 parameters(level2) #level1₊level0₊a #level1₊b #c -#level0₊d -#level1₊level0₊e -#f +#d level3 = ODESystem(Equation[], t, [], []; name = :level3) ∘ level2 parameters(level3) #level2₊level1₊level0₊a #level2₊level1₊b #level2₊c -#level2₊level0₊d -#level1₊level0₊e -#f +#d ``` ## Structural Simplify diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index 3a76f478f1..2808204d61 100644 --- a/docs/src/basics/Events.md +++ b/docs/src/basics/Events.md @@ -25,6 +25,67 @@ the event occurs). These can both be specified symbolically, but a more [general functional affect](@ref func_affects) representation is also allowed, as described below. +## Symbolic Callback Semantics + +In callbacks, there is a distinction between values of the unknowns and parameters +*before* the callback, and the desired values *after* the callback. In MTK, this +is provided by the `Pre` operator. For example, if we would like to add 1 to an +unknown `x` in a callback, the equation would look like the following: + +```julia +x ~ Pre(x) + 1 +``` + +Non `Pre`-d values will be interpreted as values *after* the callback. As such, +writing + +```julia +x ~ x + 1 +``` + +will be interpreted as an algebraic equation to be satisfied after the callback. +Since this equation obviously cannot be satisfied, an error will result. + +Callbacks must maintain the consistency of DAEs, meaning that they must satisfy +all the algebraic equations of the system after their update. However, the affect +equations often do not fully specify which unknowns/parameters should be modified +to maintain consistency. To make this clear, MTK uses the following rules: + + 1. All unknowns are treated as modifiable by the callback. In order to enforce that an unknown `x` remains the same, one can add `x ~ Pre(x)` to the affect equations. + 2. All parameters are treated as un-modifiable, *unless* they are declared as `discrete_parameters` to the callback. In order to be a discrete parameter, the parameter must be time-dependent (the terminology *discretes* here means [discrete variables](@ref save_discretes)). + +For example, consider the following system. + +```julia +@variables x(t) y(t) +@parameters p(t) +@mtkbuild sys = ODESystem([x * y ~ p, D(x) ~ 0], t) +event = [t == 1] => [x ~ Pre(x) + 1] +``` + +By default what will happen is that `x` will increase by 1, `p` will remain constant, +and `y` will change in order to compensate the increase in `x`. But what if we +wanted to keep `y` constant and change `p` instead? We could use the callback +constructor as follows: + +```julia +event = SymbolicDiscreteCallback( + [t == 1] => [x ~ Pre(x) + 1, y ~ Pre(y)], discrete_parameters = [p]) +``` + +This way, we enforce that `y` will remain the same, and `p` will change. + +!!! warning + + Symbolic affects come with the guarantee that the state after the callback + will be consistent. However, when using [general functional affects](@ref func_affects) + or [imperative affects](@ref imp_affects) one must be more careful. In + particular, one can pass in `reinitializealg` as a keyword arg to the + callback constructor to re-initialize the system. This will default to + `SciMLBase.NoInit()` in the case of symbolic affects and `SciMLBase.CheckInit()` + in the case of functional affects. This keyword should *not* be provided + if the affect is purely symbolic. + ## Continuous Events The basic purely symbolic continuous event interface to encode *one* continuous @@ -91,7 +152,7 @@ like this @variables x(t)=1 v(t)=0 root_eqs = [x ~ 0] # the event happens at the ground x(t) = 0 -affect = [v ~ -v] # the effect is that the velocity changes sign +affect = [v ~ -Pre(v)] # the effect is that the velocity changes sign @mtkbuild ball = ODESystem([D(x) ~ v D(v) ~ -9.8], t; continuous_events = root_eqs => affect) # equation => affect @@ -110,8 +171,8 @@ Multiple events? No problem! This example models a bouncing ball in 2D that is e ```@example events @variables x(t)=1 y(t)=0 vx(t)=0 vy(t)=2 -continuous_events = [[x ~ 0] => [vx ~ -vx] - [y ~ -1.5, y ~ 1.5] => [vy ~ -vy]] +continuous_events = [[x ~ 0] => [vx ~ -Pre(vx)] + [y ~ -1.5, y ~ 1.5] => [vy ~ -Pre(vy)]] @mtkbuild ball = ODESystem( [ @@ -204,7 +265,7 @@ bb_sol = solve(bb_prob, Tsit5()) plot(bb_sol) ``` -## Discrete events support +## Discrete Events In addition to continuous events, discrete events are also supported. The general interface to represent a collection of discrete events is @@ -227,13 +288,13 @@ Suppose we have a population of `N(t)` cells that can grow and die, and at time `t1` we want to inject `M` more cells into the population. We can model this by ```@example events -@parameters M tinject α +@parameters M tinject α(t) @variables N(t) Dₜ = Differential(t) eqs = [Dₜ(N) ~ α - N] # at time tinject we inject M cells -injection = (t == tinject) => [N ~ N + M] +injection = (t == tinject) => [N ~ Pre(N) + M] u0 = [N => 0.0] tspan = (0.0, 20.0) @@ -255,7 +316,7 @@ its steady-state value (which is 100). We can encode this by modifying the event to ```@example events -injection = ((t == tinject) & (N < 50)) => [N ~ N + M] +injection = ((t == tinject) & (N < 50)) => [N ~ Pre(N) + M] @mtkbuild osys = ODESystem(eqs, t, [N], [M, tinject, α]; discrete_events = injection) oprob = ODEProblem(osys, u0, tspan, p) @@ -269,16 +330,18 @@ event time, the event condition now returns false. Here we used logical and, cannot be used within symbolic expressions. Let's now also add a drug at time `tkill` that turns off production of new -cells, modeled by setting `α = 0.0` +cells, modeled by setting `α = 0.0`. Since this is a parameter we must explicitly +set it as `discrete_parameters`. ```@example events @parameters tkill # we reset the first event to just occur at tinject -injection = (t == tinject) => [N ~ N + M] +injection = (t == tinject) => [N ~ Pre(N) + M] # at time tkill we turn off production of cells -killing = (t == tkill) => [α ~ 0.0] +killing = ModelingToolkit.SymbolicDiscreteCallback( + (t == tkill) => [α ~ 0.0]; discrete_parameters = α, iv = t) tspan = (0.0, 30.0) p = [α => 100.0, tinject => 10.0, M => 50, tkill => 20.0] @@ -298,7 +361,7 @@ A preset-time event is triggered at specific set times, which can be passed in a vector like ```julia -discrete_events = [[1.0, 4.0] => [v ~ -v]] +discrete_events = [[1.0, 4.0] => [v ~ -Pre(v)]] ``` This will change the sign of `v` *only* at `t = 1.0` and `t = 4.0`. @@ -306,8 +369,9 @@ This will change the sign of `v` *only* at `t = 1.0` and `t = 4.0`. As such, our last example with treatment and killing could instead be modeled by ```@example events -injection = [10.0] => [N ~ N + M] -killing = [20.0] => [α ~ 0.0] +injection = [10.0] => [N ~ Pre(N) + M] +killing = ModelingToolkit.SymbolicDiscreteCallback( + [20.0] => [α ~ 0.0], discrete_parameters = α, iv = t) p = [α => 100.0, M => 50] @mtkbuild osys = ODESystem(eqs, t, [N], [α, M]; @@ -325,7 +389,7 @@ specify a periodic interval, pass the interval as the condition for the event. For example, ```julia -discrete_events = [1.0 => [v ~ -v]] +discrete_events = [1.0 => [v ~ -Pre(v)]] ``` will change the sign of `v` at `t = 1.0`, `2.0`, ... @@ -334,10 +398,10 @@ Finally, we note that to specify an event at precisely one time, say 2.0 below, one must still use a vector ```julia -discrete_events = [[2.0] => [v ~ -v]] +discrete_events = [[2.0] => [v ~ -Pre(v)]] ``` -## Saving discrete values +## [Saving discrete values](@id save_discretes) Time-dependent parameters which are updated in callbacks are termed as discrete variables. ModelingToolkit enables automatically saving the timeseries of these discrete variables, @@ -348,8 +412,10 @@ example: @variables x(t) @parameters c(t) +ev = ModelingToolkit.SymbolicDiscreteCallback( + 1.0 => [c ~ Pre(c) + 1], discrete_parameters = c, iv = t) @mtkbuild sys = ODESystem( - D(x) ~ c * cos(x), t, [x], [c]; discrete_events = [1.0 => [c ~ c + 1]]) + D(x) ~ c * cos(x), t, [x], [c]; discrete_events = [ev]) prob = ODEProblem(sys, [x => 0.0], (0.0, 2pi), [c => 1.0]) sol = solve(prob, Tsit5()) @@ -362,15 +428,15 @@ The solution object can also be interpolated with the discrete variables sol([1.0, 2.0], idxs = [c, c * cos(x)]) ``` -Note that only time-dependent parameters will be saved. If we repeat the above example with -this change: +Note that only time-dependent parameters that are explicitly passed as `discrete_parameters` +will be saved. If we repeat the above example with `c` not a `discrete_parameter`: ```@example events @variables x(t) -@parameters c +@parameters c(t) @mtkbuild sys = ODESystem( - D(x) ~ c * cos(x), t, [x], [c]; discrete_events = [1.0 => [c ~ c + 1]]) + D(x) ~ c * cos(x), t, [x], [c]; discrete_events = [1.0 => [c ~ Pre(c) + 1]]) prob = ODEProblem(sys, [x => 0.0], (0.0, 2pi), [c => 1.0]) sol = solve(prob, Tsit5()) diff --git a/docs/src/basics/MTKLanguage.md b/docs/src/basics/MTKLanguage.md index e91f2bcb67..ba6d2c34b5 100644 --- a/docs/src/basics/MTKLanguage.md +++ b/docs/src/basics/MTKLanguage.md @@ -203,6 +203,7 @@ getdefault(model_c3.model_a.k_array[2]) - Defining continuous events as described [here](https://docs.sciml.ai/ModelingToolkit/stable/basics/Events/#Continuous-Events). - If this block is not defined in the model, no continuous events will be added. + - Discrete parameters and other keyword arguments should be specified in a vector, as seen below. ```@example mtkmodel-example using ModelingToolkit @@ -210,7 +211,7 @@ using ModelingToolkit: t @mtkmodel M begin @parameters begin - k + k(t) end @variables begin x(t) @@ -223,6 +224,7 @@ using ModelingToolkit: t @continuous_events begin [x ~ 1.5] => [x ~ 5, y ~ 5] [t ~ 4] => [x ~ 10] + [t ~ 5] => [k ~ 3], [discrete_parameters = k] end end ``` @@ -231,13 +233,14 @@ end - Defining discrete events as described [here](https://docs.sciml.ai/ModelingToolkit/stable/basics/Events/#Discrete-events-support). - If this block is not defined in the model, no discrete events will be added. + - Discrete parameters and other keyword arguments should be specified in a vector, as seen below. ```@example mtkmodel-example using ModelingToolkit @mtkmodel M begin @parameters begin - k + k(t) end @variables begin x(t) @@ -248,7 +251,8 @@ using ModelingToolkit D(y) ~ -k end @discrete_events begin - (t == 1.5) => [x ~ x + 5, y ~ 5] + (t == 1.5) => [x ~ Pre(x) + 5, y ~ 5] + (t == 2.5) => [k ~ Pre(k) * 2], [discrete_parameters = k] end end ``` diff --git a/docs/src/systems/DiscreteSystem.md b/docs/src/systems/DiscreteSystem.md index f8a71043ab..55a02e5714 100644 --- a/docs/src/systems/DiscreteSystem.md +++ b/docs/src/systems/DiscreteSystem.md @@ -12,7 +12,6 @@ DiscreteSystem - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns in the discrete system. - `get_ps(sys)` or `parameters(sys)`: The parameters of the discrete system. - `get_iv(sys)`: The independent variable of the discrete system - - `discrete_events(sys)`: The set of discrete events in the discrete system. ## Transformations diff --git a/docs/src/systems/ImplicitDiscreteSystem.md b/docs/src/systems/ImplicitDiscreteSystem.md index d69f88f106..d687502b49 100644 --- a/docs/src/systems/ImplicitDiscreteSystem.md +++ b/docs/src/systems/ImplicitDiscreteSystem.md @@ -12,7 +12,6 @@ ImplicitDiscreteSystem - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns in the implicit discrete system. - `get_ps(sys)` or `parameters(sys)`: The parameters of the implicit discrete system. - `get_iv(sys)`: The independent variable of the implicit discrete system - - `discrete_events(sys)`: The set of discrete events in the implicit discrete system. ## Transformations diff --git a/docs/src/tutorials/fmi.md b/docs/src/tutorials/fmi.md index ef00477c78..0e01393652 100644 --- a/docs/src/tutorials/fmi.md +++ b/docs/src/tutorials/fmi.md @@ -94,7 +94,8 @@ we will create a model from a CoSimulation FMU. ```@example fmi fmu = loadFMU("SpringPendulum1D", "Dymola", "2023x", "3.0"; type = :CS) @named inner = ModelingToolkit.FMIComponent( - Val(3); fmu, communication_step_size = 0.001, type = :CS) + Val(3); fmu, communication_step_size = 0.001, type = :CS, + reinitializealg = BrownFullBasicInit()) ``` This FMU has fewer equations, partly due to missing aliasing variables and partly due to being a CS FMU. @@ -170,7 +171,8 @@ end `a` and `b` are inputs, `c` is a state, and `out` and `out2` are outputs of the component. ```@repl fmi -@named adder = ModelingToolkit.FMIComponent(Val(2); fmu, type = :ME); +@named adder = ModelingToolkit.FMIComponent( + Val(2); fmu, type = :ME, reinitializealg = BrownFullBasicInit()); isinput(adder.a) isinput(adder.b) isoutput(adder.out) diff --git a/ext/MTKFMIExt.jl b/ext/MTKFMIExt.jl index 5cfe9a82ef..32023d0749 100644 --- a/ext/MTKFMIExt.jl +++ b/ext/MTKFMIExt.jl @@ -93,7 +93,7 @@ with the name `namespace__variable`. - `name`: The name of the system. """ function MTK.FMIComponent(::Val{Ver}; fmu = nothing, tolerance = 1e-6, - communication_step_size = nothing, reinitializealg = SciMLBase.NoInit(), type, name) where {Ver} + communication_step_size = nothing, reinitializealg = nothing, type, name) where {Ver} if Ver != 2 && Ver != 3 throw(ArgumentError("FMI Version must be `2` or `3`")) end @@ -238,7 +238,7 @@ function MTK.FMIComponent(::Val{Ver}; fmu = nothing, tolerance = 1e-6, finalize_affect = MTK.FunctionalAffect(fmiFinalize!, [], [wrapper], []) step_affect = MTK.FunctionalAffect(Returns(nothing), [], [], []) instance_management_callback = MTK.SymbolicDiscreteCallback( - (t != t - 1), step_affect; finalize = finalize_affect, reinitializealg = reinitializealg) + (t != t - 1), step_affect; finalize = finalize_affect, reinitializealg) push!(params, wrapper) append!(observed, der_observed) @@ -279,8 +279,7 @@ function MTK.FMIComponent(::Val{Ver}; fmu = nothing, tolerance = 1e-6, fmiCSStep!; observed = cb_observed, modified = cb_modified, ctx = _functor) instance_management_callback = MTK.SymbolicDiscreteCallback( communication_step_size, step_affect; initialize = initialize_affect, - finalize = finalize_affect, reinitializealg = reinitializealg - ) + finalize = finalize_affect, reinitializealg) # guarded in case there are no outputs/states and the variable is `[]`. symbolic_type(__mtk_internal_o) == NotSymbolic() || push!(params, __mtk_internal_o) diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index d2d20f902e..ec29b73c6b 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -54,6 +54,7 @@ import Moshi using Moshi.Data: @data using NonlinearSolve import SCCNonlinearSolve +using ImplicitDiscreteSolve using Reexport using RecursiveArrayTools import Graphs: SimpleDiGraph, add_edge!, incidence_matrix @@ -158,7 +159,6 @@ include("systems/model_parsing.jl") include("systems/connectors.jl") include("systems/analysis_points.jl") include("systems/imperative_affect.jl") -include("systems/callbacks.jl") include("systems/codegen_utils.jl") include("systems/problem_utils.jl") include("linearization.jl") @@ -168,19 +168,20 @@ include("systems/optimization/optimizationsystem.jl") include("systems/optimization/modelingtoolkitize.jl") include("systems/nonlinear/nonlinearsystem.jl") -include("systems/nonlinear/homotopy_continuation.jl") +include("systems/discrete_system/discrete_system.jl") +include("systems/discrete_system/implicit_discrete_system.jl") +include("systems/callbacks.jl") + include("systems/diffeqs/odesystem.jl") include("systems/diffeqs/sdesystem.jl") include("systems/diffeqs/abstractodesystem.jl") +include("systems/nonlinear/homotopy_continuation.jl") include("systems/nonlinear/modelingtoolkitize.jl") include("systems/nonlinear/initializesystem.jl") include("systems/diffeqs/first_order_transform.jl") include("systems/diffeqs/modelingtoolkitize.jl") include("systems/diffeqs/basic_transformations.jl") -include("systems/discrete_system/discrete_system.jl") -include("systems/discrete_system/implicit_discrete_system.jl") - include("systems/jumps/jumpsystem.jl") include("systems/pde/pdesystem.jl") @@ -297,12 +298,13 @@ export PDESystem export Differential, expand_derivatives, @derivatives export Equation, ConstrainedEquation export Term, Sym -export SymScope, LocalScope, ParentScope, DelayParentScope, GlobalScope +export SymScope, LocalScope, ParentScope, GlobalScope export independent_variable, equations, controls, observed, full_equations export initialization_equations, guesses, defaults, parameter_dependencies, hierarchy export structural_simplify, expand_connections, linearize, linearization_function, LinearizationProblem export solve +export Pre export calculate_jacobian, generate_jacobian, generate_function, generate_custom_function, generate_W diff --git a/src/inputoutput.jl b/src/inputoutput.jl index 8c1a084c9b..18d7830cba 100644 --- a/src/inputoutput.jl +++ b/src/inputoutput.jl @@ -163,7 +163,7 @@ has_var(ex, x) = x ∈ Set(get_variables(ex)) (f_oop, f_ip), x_sym, p_sym, io_sys = generate_control_function( sys::AbstractODESystem, inputs = unbound_inputs(sys), - disturbance_inputs = nothing; + disturbance_inputs = disturbances(sys); implicit_dae = false, simplify = false, ) @@ -179,9 +179,6 @@ The return values also include the chosen state-realization (the remaining unkno If `disturbance_inputs` is an array of variables, the generated dynamics function will preserve any state and dynamics associated with disturbance inputs, but the disturbance inputs themselves will (by default) not be included as inputs to the generated function. The use case for this is to generate dynamics for state observers that estimate the influence of unmeasured disturbances, and thus require unknown variables for the disturbance model, but without disturbance inputs since the disturbances are not available for measurement. To add an input argument corresponding to the disturbance inputs, either include the disturbance inputs among the control inputs, or set `disturbance_argument=true`, in which case an additional input argument `w` is added to the generated function `(x,u,p,t,w)->rhs`. -!!! note "Un-simplified system" - This function expects `sys` to be un-simplified, i.e., `structural_simplify` or `@mtkbuild` should not be called on the system before passing it into this function. `generate_control_function` calls a special version of `structural_simplify` internally. - # Example ``` @@ -200,18 +197,18 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu simplify = false, eval_expression = false, eval_module = @__MODULE__, + check_simplified = true, kwargs...) + # Remove this when the ControlFunction gets merged. + if check_simplified && !iscomplete(sys) + error("A completed `ODESystem` is required. Call `complete` or `structural_simplify` on the system before creating the control function.") + end isempty(inputs) && @warn("No unbound inputs were found in system.") - if disturbance_inputs !== nothing # add to inputs for the purposes of io processing inputs = [inputs; disturbance_inputs] end - if !iscomplete(sys) - sys, _ = io_preprocessing(sys, inputs, []; simplify, kwargs...) - end - dvs = unknowns(sys) ps = parameters(sys; initial_parameters = true) ps = setdiff(ps, inputs) @@ -220,7 +217,7 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu inputs = setdiff(inputs, disturbance_inputs) # ps = [ps; disturbance_inputs] end - inputs = map(x -> time_varying_as_func(value(x), sys), inputs) + inputs = map(value, inputs) disturbance_inputs = unwrap.(disturbance_inputs) eqs = [eq for eq in full_equations(sys)] @@ -257,8 +254,11 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu (; f, dvs, ps, io_sys = sys) end -function inputs_to_parameters!(state::TransformationState, io) - check_bound = io === nothing +""" +Turn input variables into parameters of the system. +""" +function inputs_to_parameters!(state::TransformationState, inputsyms) + check_bound = inputsyms === nothing @unpack structure, fullvars, sys = state @unpack var_to_diff, graph, solvable_graph = structure @assert solvable_graph === nothing @@ -287,7 +287,7 @@ function inputs_to_parameters!(state::TransformationState, io) push!(new_fullvars, v) end end - ninputs == 0 && return (state, 1:0) + ninputs == 0 && return state nvars = ndsts(graph) - ninputs new_graph = BipartiteGraph(nsrcs(graph), nvars, Val(false)) @@ -316,24 +316,11 @@ function inputs_to_parameters!(state::TransformationState, io) @set! sys.unknowns = setdiff(unknowns(sys), keys(input_to_parameters)) ps = parameters(sys) - if io !== nothing - inputs, = io - # Change order of new parameters to correspond to user-provided order in argument `inputs` - d = Dict{Any, Int}() - for (i, inp) in enumerate(new_parameters) - d[inp] = i - end - permutation = [d[i] for i in inputs] - new_parameters = new_parameters[permutation] - end - @set! sys.ps = [ps; new_parameters] - @set! state.sys = sys @set! state.fullvars = new_fullvars @set! state.structure = structure - base_params = length(ps) - return state, (base_params + 1):(base_params + length(new_parameters)) # (1:length(new_parameters)) .+ base_params + return state end """ @@ -359,7 +346,7 @@ function get_disturbance_system(dist::DisturbanceModel{<:ODESystem}) end """ - (f_oop, f_ip), augmented_sys, dvs, p = add_input_disturbance(sys, dist::DisturbanceModel, inputs = nothing) + (f_oop, f_ip), augmented_sys, dvs, p = add_input_disturbance(sys, dist::DisturbanceModel, inputs = Any[]) Add a model of an unmeasured disturbance to `sys`. The disturbance model is an instance of [`DisturbanceModel`](@ref). @@ -408,13 +395,13 @@ model_outputs = [model.inertia1.w, model.inertia2.w, model.inertia1.phi, model.i `f_oop` will have an extra state corresponding to the integrator in the disturbance model. This state will not be affected by any input, but will affect the dynamics from where it enters, in this case it will affect additively from `model.torque.tau.u`. """ -function add_input_disturbance(sys, dist::DisturbanceModel, inputs = nothing; kwargs...) +function add_input_disturbance(sys, dist::DisturbanceModel, inputs = Any[]; kwargs...) t = get_iv(sys) @variables d(t)=0 [disturbance = true] @variables u(t)=0 [input = true] # New system input dsys = get_disturbance_system(dist) - if inputs === nothing + if isempty(inputs) all_inputs = [u] else i = findfirst(isequal(dist.input), inputs) @@ -429,8 +416,9 @@ function add_input_disturbance(sys, dist::DisturbanceModel, inputs = nothing; kw dist.input ~ u + dsys.output.u[1]] augmented_sys = ODESystem(eqs, t, systems = [dsys], name = gensym(:outer)) augmented_sys = extend(augmented_sys, sys) + ssys = structural_simplify(augmented_sys, inputs = all_inputs, disturbance_inputs = [d]) - f, dvs, p, io_sys = generate_control_function(augmented_sys, all_inputs, + f, dvs, p, io_sys = generate_control_function(ssys, all_inputs, [d]; kwargs...) f, augmented_sys, dvs, p, io_sys end diff --git a/src/linearization.jl b/src/linearization.jl index 77f4422b63..12340435b8 100644 --- a/src/linearization.jl +++ b/src/linearization.jl @@ -58,9 +58,8 @@ function linearization_function(sys::AbstractSystem, inputs, outputs = mapreduce(vcat, outputs; init = []) do var symbolic_type(var) == ArraySymbolic() ? collect(var) : [var] end - ssys, diff_idxs, alge_idxs, input_idxs = io_preprocessing(sys, inputs, outputs; - simplify, - kwargs...) + ssys = structural_simplify(sys; inputs, outputs, simplify, kwargs...) + diff_idxs, alge_idxs = eq_idxs(ssys) if zero_dummy_der dummyder = setdiff(unknowns(ssys), unknowns(sys)) defs = Dict(x => 0.0 for x in dummyder) @@ -87,9 +86,9 @@ function linearization_function(sys::AbstractSystem, inputs, p = parameter_values(prob) t0 = current_time(prob) - inputvals = [p[idx] for idx in input_idxs] + inputvals = [prob.ps[i] for i in inputs] - hp_fun = let fun = h, setter = setp_oop(sys, input_idxs) + hp_fun = let fun = h, setter = setp_oop(sys, inputs) function hpf(du, input, u, p, t) p = setter(p, input) fun(du, u, p, t) @@ -113,7 +112,7 @@ function linearization_function(sys::AbstractSystem, inputs, # observed function is a `GeneratedFunctionWrapper` with iip component h_jac = PreparedJacobian{true}(h, similar(prob.u0, size(outputs)), autodiff, prob.u0, DI.Constant(p), DI.Constant(t0)) - pf_fun = let fun = prob.f, setter = setp_oop(sys, input_idxs) + pf_fun = let fun = prob.f, setter = setp_oop(sys, inputs) function pff(du, input, u, p, t) p = setter(p, input) SciMLBase.ParamJacobianWrapper(fun, t, u)(du, p) @@ -127,12 +126,23 @@ function linearization_function(sys::AbstractSystem, inputs, end lin_fun = LinearizationFunction( - diff_idxs, alge_idxs, input_idxs, length(unknowns(sys)), + diff_idxs, alge_idxs, inputs, length(unknowns(sys)), prob, h, u0 === nothing ? nothing : similar(u0), uf_jac, h_jac, pf_jac, hp_jac, initializealg, initialization_kwargs) return lin_fun, sys end +""" +Return the set of indexes of differential equations and algebraic equations in the simplified system. +""" +function eq_idxs(sys::AbstractSystem) + eqs = equations(sys) + alge_idxs = findall(!isdiffeq, eqs) + diff_idxs = setdiff(1:length(eqs), alge_idxs) + + diff_idxs, alge_idxs +end + """ $(TYPEDEF) @@ -192,7 +202,7 @@ A callable struct which linearizes a system. $(TYPEDFIELDS) """ struct LinearizationFunction{ - DI <: AbstractVector{Int}, AI <: AbstractVector{Int}, II, P <: ODEProblem, + DI <: AbstractVector{Int}, AI <: AbstractVector{Int}, I, P <: ODEProblem, H, C, J1, J2, J3, J4, IA <: SciMLBase.DAEInitializationAlgorithm, IK} """ The indexes of differential equations in the linearized system. @@ -206,7 +216,7 @@ struct LinearizationFunction{ The indexes of parameters in the linearized system which represent input variables. """ - input_idxs::II + inputs::I """ The number of unknowns in the linearized system. """ @@ -281,6 +291,7 @@ function (linfun::LinearizationFunction)(u, p, t) end fun = linfun.prob.f + input_vals = [linfun.prob.ps[i] for i in linfun.inputs] if u !== nothing # Handle systems without unknowns linfun.num_states == length(u) || error("Number of unknown variables ($(linfun.num_states)) does not match the number of input unknowns ($(length(u)))") @@ -294,15 +305,15 @@ function (linfun::LinearizationFunction)(u, p, t) end fg_xz = linfun.uf_jac(u, DI.Constant(p), DI.Constant(t)) h_xz = linfun.h_jac(u, DI.Constant(p), DI.Constant(t)) - fg_u = linfun.pf_jac([p[idx] for idx in linfun.input_idxs], + fg_u = linfun.pf_jac(input_vals, DI.Constant(u), DI.Constant(p), DI.Constant(t)) else linfun.num_states == 0 || error("Number of unknown variables (0) does not match the number of input unknowns ($(length(u)))") fg_xz = zeros(0, 0) - h_xz = fg_u = zeros(0, length(linfun.input_idxs)) + h_xz = fg_u = zeros(0, length(linfun.inputs)) end - h_u = linfun.hp_jac([p[idx] for idx in linfun.input_idxs], + h_u = linfun.hp_jac(input_vals, DI.Constant(u), DI.Constant(p), DI.Constant(t)) (f_x = fg_xz[linfun.diff_idxs, linfun.diff_idxs], f_z = fg_xz[linfun.diff_idxs, linfun.alge_idxs], @@ -482,9 +493,8 @@ function linearize_symbolic(sys::AbstractSystem, inputs, outputs; simplify = false, allow_input_derivatives = false, eval_expression = false, eval_module = @__MODULE__, kwargs...) - sys, diff_idxs, alge_idxs, input_idxs = io_preprocessing( - sys, inputs, outputs; simplify, - kwargs...) + sys = structural_simplify(sys; inputs, outputs, simplify, kwargs...) + diff_idxs, alge_idxs = eq_idxs(sys) sts = unknowns(sys) t = get_iv(sys) ps = parameters(sys; initial_parameters = true) @@ -545,10 +555,14 @@ function linearize_symbolic(sys::AbstractSystem, inputs, (; A, B, C, D, f_x, f_z, g_x, g_z, f_u, g_u, h_x, h_z, h_u), sys end -function markio!(state, orig_inputs, inputs, outputs; check = true) +""" +Modify the variable metadata of system variables to indicate which ones are inputs, outputs, and disturbances. Needed for `inputs`, `outputs`, `disturbances`, `unbound_inputs`, `unbound_outputs` to return the proper subsets. +""" +function markio!(state, orig_inputs, inputs, outputs, disturbances; check = true) fullvars = get_fullvars(state) inputset = Dict{Any, Bool}(i => false for i in inputs) outputset = Dict{Any, Bool}(o => false for o in outputs) + disturbanceset = Dict{Any, Bool}(d => false for d in disturbances) for (i, v) in enumerate(fullvars) if v in keys(inputset) if v in keys(outputset) @@ -570,6 +584,13 @@ function markio!(state, orig_inputs, inputs, outputs; check = true) v = setio(v, false, false) fullvars[i] = v end + + if v in keys(disturbanceset) + v = setio(v, true, false) + v = setdisturbance(v, true) + disturbanceset[v] = true + fullvars[i] = v + end end if check ikeys = keys(filter(!last, inputset)) @@ -578,11 +599,16 @@ function markio!(state, orig_inputs, inputs, outputs; check = true) "Some specified inputs were not found in system. The following variables were not found ", ikeys) end - end - check && (all(values(outputset)) || - error( - "Some specified outputs were not found in system. The following Dict indicates the found variables ", + dkeys = keys(filter(!last, disturbanceset)) + if !isempty(dkeys) + error( + "Specified disturbance inputs were not found in system. The following variables were not found ", + ikeys) + end + (all(values(outputset)) || error( + "Some specified outputs were not found in system. The following Dict indicates the found variables ", outputset)) + end state, orig_inputs end diff --git a/src/parameters.jl b/src/parameters.jl index 91121b7cbb..de7a722af3 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -63,7 +63,7 @@ toparam(s::Num) = wrap(toparam(value(s))) Maps the variable to an unknown. """ tovar(s::Symbolic) = setmetadata(s, MTKVariableTypeCtx, VARIABLE) -tovar(s::Num) = Num(tovar(value(s))) +tovar(s::Union{Num, Symbolics.Arr}) = wrap(tovar(unwrap(s))) """ $(SIGNATURES) diff --git a/src/structural_transformation/bipartite_tearing/modia_tearing.jl b/src/structural_transformation/bipartite_tearing/modia_tearing.jl index 5da873afdf..59b32abd56 100644 --- a/src/structural_transformation/bipartite_tearing/modia_tearing.jl +++ b/src/structural_transformation/bipartite_tearing/modia_tearing.jl @@ -96,7 +96,7 @@ function tear_graph_modia(structure::SystemStructure, isder::F = nothing, ieqs = Int[] filtered_vars = BitSet() free_eqs = free_equations(graph, var_sccs, var_eq_matching, varfilter) - is_overdetemined = !isempty(free_eqs) + is_overdetermined = !isempty(free_eqs) for vars in var_sccs for var in vars if varfilter(var) @@ -112,7 +112,7 @@ function tear_graph_modia(structure::SystemStructure, isder::F = nothing, filtered_vars, isder) # If the systems is overdetemined, we cannot assume the free equations # will not form algebraic loops with equations in the sccs. - if !is_overdetemined + if !is_overdetermined vargraph.ne = 0 for var in vars vargraph.matching[var] = unassigned @@ -121,7 +121,7 @@ function tear_graph_modia(structure::SystemStructure, isder::F = nothing, empty!(ieqs) empty!(filtered_vars) end - if is_overdetemined + if is_overdetermined free_vars = findall(x -> !(x isa Int), var_eq_matching) tear_graph_block_modia!(var_eq_matching, ict, solvable_graph, free_eqs, BitSet(free_vars), isder) diff --git a/src/structural_transformation/utils.jl b/src/structural_transformation/utils.jl index 14628f2958..7834e33b61 100644 --- a/src/structural_transformation/utils.jl +++ b/src/structural_transformation/utils.jl @@ -218,6 +218,7 @@ function find_eq_solvables!(state::TearingState, ieq, to_rm = Int[], coeffs = no all_int_vars = true coeffs === nothing || empty!(coeffs) empty!(to_rm) + for j in 𝑠neighbors(graph, ieq) var = fullvars[j] isirreducible(var) && (all_int_vars = false; continue) @@ -550,6 +551,7 @@ end function _distribute_shift(expr, shift) if iscall(expr) op = operation(expr) + (op isa Pre || op isa Initial) && return expr args = arguments(expr) if ModelingToolkit.isvariable(expr) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 6c00806b7d..1da2a8ff73 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -1190,55 +1190,6 @@ function ParentScope(sym::Union{Num, Symbolic, Symbolics.Arr{Num}}) end end -""" - $(TYPEDEF) - -Denotes that a variable belongs to a system that is at least `N + 1` levels up in the -hierarchy from the system whose equations it is involved in. It is namespaced by the -first `N` parents and not namespaced by the `N+1`th parent in the hierarchy. The scope -of the variable after this point is given by `parent`. - -In other words, this scope delays applying `ParentScope` by `N` levels, and applies -`LocalScope` in the meantime. - -# Fields - -$(TYPEDFIELDS) -""" -struct DelayParentScope <: SymScope - parent::SymScope - N::Int -end - -""" - $(TYPEDSIGNATURES) - -Apply `DelayParentScope` to `sym`, with a delay of `N` and `parent` being `LocalScope`. -""" -function DelayParentScope(sym::Union{Num, Symbolic, Symbolics.Arr{Num}}, N) - Base.depwarn( - "`DelayParentScope` is deprecated and will be removed soon", :DelayParentScope) - apply_to_variables(sym) do sym - if iscall(sym) && operation(sym) == getindex - args = arguments(sym) - a1 = setmetadata(args[1], SymScope, - DelayParentScope(getmetadata(value(args[1]), SymScope, LocalScope()), N)) - maketerm(typeof(sym), operation(sym), [a1, args[2:end]...], - metadata(sym)) - else - setmetadata(sym, SymScope, - DelayParentScope(getmetadata(value(sym), SymScope, LocalScope()), N)) - end - end -end - -""" - $(TYPEDSIGNATURES) - -Apply `DelayParentScope` to `sym`, with a delay of `1` and `parent` being `LocalScope`. -""" -DelayParentScope(sym::Union{Num, Symbolic, Symbolics.Arr{Num}}) = DelayParentScope(sym, 1) - """ $(TYPEDEF) @@ -1291,15 +1242,6 @@ function renamespace(sys, x) rename(x, renamespace(getname(sys), getname(x)))::T elseif scope isa ParentScope setmetadata(x, SymScope, scope.parent)::T - elseif scope isa DelayParentScope - if scope.N > 0 - x = setmetadata(x, SymScope, - DelayParentScope(scope.parent, scope.N - 1)) - rename(x, renamespace(getname(sys), getname(x)))::T - else - #rename(x, renamespace(getname(sys), getname(x))) - setmetadata(x, SymScope, scope.parent)::T - end else # GlobalScope x::T end @@ -1849,20 +1791,6 @@ function isaffine(sys::AbstractSystem) all(isaffine(r, unknowns(sys)) for r in rhs) end -function time_varying_as_func(x, sys::AbstractTimeDependentSystem) - # if something is not x(t) (the current unknown) - # but is `x(t-1)` or something like that, pass in `x` as a callable function rather - # than pass in a value in place of x(t). - # - # This is done by just making `x` the argument of the function. - if iscall(x) && - issym(operation(x)) && - !(length(arguments(x)) == 1 && isequal(arguments(x)[1], get_iv(sys))) - return operation(x) - end - return x -end - """ $(SIGNATURES) @@ -2400,11 +2328,7 @@ function default_to_parentscope(v) uv = unwrap(v) uv isa Symbolic || return v apply_to_variables(v) do sym - if !hasmetadata(uv, SymScope) - ParentScope(sym) - else - sym - end + ParentScope(sym) end end @@ -2562,21 +2486,6 @@ function eliminate_constants(sys::AbstractSystem) return sys end -function io_preprocessing(sys::AbstractSystem, inputs, - outputs; simplify = false, kwargs...) - sys, input_idxs = structural_simplify(sys, (inputs, outputs); simplify, kwargs...) - - eqs = equations(sys) - alg_start_idx = findfirst(!isdiffeq, eqs) - if alg_start_idx === nothing - alg_start_idx = length(eqs) + 1 - end - diff_idxs = 1:(alg_start_idx - 1) - alge_idxs = alg_start_idx:length(eqs) - - sys, diff_idxs, alge_idxs, input_idxs -end - @latexrecipe function f(sys::AbstractSystem) return latexify(equations(sys)) end diff --git a/src/systems/analysis_points.jl b/src/systems/analysis_points.jl index 0d1a2830cf..29f1e777ac 100644 --- a/src/systems/analysis_points.jl +++ b/src/systems/analysis_points.jl @@ -958,7 +958,6 @@ function linearization_function(sys::AbstractSystem, end sys = handle_loop_openings(sys, map(AnalysisPoint, collect(loop_openings))) - return linearization_function(system_modifier(sys), input_vars, output_vars; kwargs...) end diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index 07809bf611..785dc49338 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -1,15 +1,4 @@ -#################################### system operations ##################################### -has_continuous_events(sys::AbstractSystem) = isdefined(sys, :continuous_events) -function get_continuous_events(sys::AbstractSystem) - has_continuous_events(sys) || return SymbolicContinuousCallback[] - getfield(sys, :continuous_events) -end - -has_discrete_events(sys::AbstractSystem) = isdefined(sys, :discrete_events) -function get_discrete_events(sys::AbstractSystem) - has_discrete_events(sys) || return SymbolicDiscreteCallback[] - getfield(sys, :discrete_events) -end +abstract type AbstractCallback end struct FunctionalAffect f::Any @@ -38,7 +27,7 @@ function FunctionalAffect(; f, sts, pars, discretes, ctx = nothing) FunctionalAffect(f, sts, pars, discretes, ctx) end -func(f::FunctionalAffect) = f.f +func(a::FunctionalAffect) = a.f context(a::FunctionalAffect) = a.ctx parameters(a::FunctionalAffect) = a.pars parameters_syms(a::FunctionalAffect) = a.pars_syms @@ -62,33 +51,118 @@ function Base.hash(a::FunctionalAffect, s::UInt) hash(a.ctx, s) end -namespace_affect(affect, s) = namespace_equation(affect, s) -function namespace_affect(affect::FunctionalAffect, s) - FunctionalAffect(func(affect), - renamespace.((s,), unknowns(affect)), - unknowns_syms(affect), - renamespace.((s,), parameters(affect)), - parameters_syms(affect), - renamespace.((s,), discretes(affect)), - context(affect)) -end - function has_functional_affect(cb) (affects(cb) isa FunctionalAffect || affects(cb) isa ImperativeAffect) end -function vars!(vars, aff::FunctionalAffect; op = Differential) +struct AffectSystem + """The internal implicit discrete system whose equations are solved to obtain values after the affect.""" + system::ImplicitDiscreteSystem + """Unknowns of the parent ODESystem whose values are modified or accessed by the affect.""" + unknowns::Vector + """Parameters of the parent ODESystem whose values are accessed by the affect.""" + parameters::Vector + """Parameters of the parent ODESystem whose values are modified by the affect.""" + discretes::Vector + """Maps the symbols of unknowns/observed in the ImplicitDiscreteSystem to its corresponding unknown/parameter in the parent system.""" + aff_to_sys::Dict +end + +system(a::AffectSystem) = a.system +discretes(a::AffectSystem) = a.discretes +unknowns(a::AffectSystem) = a.unknowns +parameters(a::AffectSystem) = a.parameters +aff_to_sys(a::AffectSystem) = a.aff_to_sys +all_equations(a::AffectSystem) = vcat(equations(system(a)), observed(system(a))) + +function Base.show(iio::IO, aff::AffectSystem) + println(iio, "Affect system defined by equations:") + eqs = all_equations(aff) + show(iio, eqs) +end + +function Base.:(==)(a1::AffectSystem, a2::AffectSystem) + isequal(system(a1), system(a2)) && + isequal(discretes(a1), discretes(a2)) && + isequal(unknowns(a1), unknowns(a2)) && + isequal(parameters(a1), parameters(a2)) && + isequal(aff_to_sys(a1), aff_to_sys(a2)) +end + +function Base.hash(a::AffectSystem, s::UInt) + s = hash(system(a), s) + s = hash(unknowns(a), s) + s = hash(parameters(a), s) + s = hash(discretes(a), s) + hash(aff_to_sys(a), s) +end + +function vars!(vars, aff::Union{FunctionalAffect, AffectSystem}; op = Differential) for var in Iterators.flatten((unknowns(aff), parameters(aff), discretes(aff))) vars!(vars, var) end - return vars + vars end -#################################### continuous events ##################################### +""" + Pre(x) -const NULL_AFFECT = Equation[] +The `Pre` operator. Used by the callback system to indicate the value of a parameter or variable +before the callback is triggered. """ - SymbolicContinuousCallback(eqs::Vector{Equation}, affect, affect_neg, rootfind) +struct Pre <: Symbolics.Operator end +Pre(x) = Pre()(x) +SymbolicUtils.promote_symtype(::Type{Pre}, T) = T +SymbolicUtils.isbinop(::Pre) = false +Base.nameof(::Pre) = :Pre +Base.show(io::IO, x::Pre) = print(io, "Pre") +input_timedomain(::Pre, _ = nothing) = ContinuousClock() +output_timedomain(::Pre, _ = nothing) = ContinuousClock() +unPre(x::Num) = unPre(unwrap(x)) +unPre(x::BasicSymbolic) = (iscall(x) && operation(x) isa Pre) ? only(arguments(x)) : x + +function (p::Pre)(x) + iw = Symbolics.iswrapped(x) + x = unwrap(x) + # non-symbolic values don't change + if symbolic_type(x) == NotSymbolic() + return x + end + # differential variables are default-toterm-ed + if iscall(x) && operation(x) isa Differential + x = default_toterm(x) + end + # don't double wrap + iscall(x) && operation(x) isa Pre && return x + result = if symbolic_type(x) == ArraySymbolic() + # create an array for `Pre(array)` + Symbolics.array_term(p, x) + elseif iscall(x) && operation(x) == getindex + # instead of `Pre(x[1])` create `Pre(x)[1]` + # which allows parameter indexing to handle this case automatically. + arr = arguments(x)[1] + term(getindex, p(arr), arguments(x)[2:end]...) + else + term(p, x) + end + # the result should be a parameter + result = toparam(result) + if iw + result = wrap(result) + end + return result +end +haspre(eq::Equation) = haspre(eq.lhs) || haspre(eq.rhs) +haspre(O) = recursive_hasoperator(Pre, O) + +############################### +###### Continuous events ###### +############################### +const Affect = Union{AffectSystem, FunctionalAffect, ImperativeAffect} + +""" + SymbolicContinuousCallback(eqs::Vector{Equation}, affect = nothing, iv = nothing; + affect_neg = affect, initialize = nothing, finalize = nothing, rootfind = SciMLBase.LeftRootFind, alg_eqs = Equation[]) A [`ContinuousCallback`](@ref SciMLBase.ContinuousCallback) specified symbolically. Takes a vector of equations `eq` as well as the positive-edge `affect` and negative-edge `affect_neg` that apply when *any* of `eq` are satisfied. @@ -128,995 +202,765 @@ Affects (i.e. `affect` and `affect_neg`) can be specified as either: + `ctx` is a user-defined context object passed to `f!` when invoked. This value is aliased for each problem. * A [`ImperativeAffect`](@ref); refer to its documentation for details. -DAEs will be reinitialized using `reinitializealg` (which defaults to `SciMLBase.CheckInit`) after callbacks are applied. -This reinitialization algorithm ensures that the DAE is satisfied after the callback runs. The default value of `CheckInit` will simply validate -that the newly-assigned values indeed satisfy the algebraic system; see the documentation on DAE initialization for a more detailed discussion of -initialization. +`reinitializealg` is used to set how the system will be reinitialized after the callback. +- Symbolic affects have reinitialization built in. In this case the algorithm will default to SciMLBase.NoInit(), and should **not** be provided. +- Functional and imperative affects will default to SciMLBase.CheckInit(), which will error if the system is not properly reinitialized after the callback. If your system is a DAE, pass in an algorithm like SciMLBase.BrownBasicFullInit() to properly re-initialize. -Initial and final affects can also be specified with SCC, which are specified identically to positive and negative edge affects. Initialization affects +Initial and final affects can also be specified identically to positive and negative edge affects. Initialization affects will run as soon as the solver starts, while finalization affects will be executed after termination. """ -struct SymbolicContinuousCallback - eqs::Vector{Equation} - initialize::Union{Vector{Equation}, FunctionalAffect, ImperativeAffect} - finalize::Union{Vector{Equation}, FunctionalAffect, ImperativeAffect} - affect::Union{Vector{Equation}, FunctionalAffect, ImperativeAffect} - affect_neg::Union{Vector{Equation}, FunctionalAffect, ImperativeAffect, Nothing} - rootfind::SciMLBase.RootfindOpt +struct SymbolicContinuousCallback <: AbstractCallback + conditions::Vector{Equation} + affect::Union{Affect, Nothing} + affect_neg::Union{Affect, Nothing} + initialize::Union{Affect, Nothing} + finalize::Union{Affect, Nothing} + rootfind::Union{Nothing, SciMLBase.RootfindOpt} reinitializealg::SciMLBase.DAEInitializationAlgorithm - function SymbolicContinuousCallback(; - eqs::Vector{Equation}, - affect = NULL_AFFECT, + + function SymbolicContinuousCallback( + conditions::Union{Equation, Vector{Equation}}, + affect = nothing; affect_neg = affect, - initialize = NULL_AFFECT, - finalize = NULL_AFFECT, + initialize = nothing, + finalize = nothing, rootfind = SciMLBase.LeftRootFind, - reinitializealg = SciMLBase.CheckInit()) - new(eqs, initialize, finalize, make_affect(affect), - make_affect(affect_neg), rootfind, reinitializealg) + reinitializealg = nothing, + kwargs...) + conditions = (conditions isa AbstractVector) ? conditions : [conditions] + + if isnothing(reinitializealg) + if any(a -> (a isa FunctionalAffect || a isa ImperativeAffect), + [affect, affect_neg, initialize, finalize]) + reinitializealg = SciMLBase.CheckInit() + else + reinitializealg = SciMLBase.NoInit() + end + end + + new(conditions, make_affect(affect; kwargs...), + make_affect(affect_neg; kwargs...), + make_affect(initialize; kwargs...), make_affect( + finalize; kwargs...), + rootfind, reinitializealg) end # Default affect to nothing end -make_affect(affect) = affect -make_affect(affect::Tuple) = FunctionalAffect(affect...) -make_affect(affect::NamedTuple) = FunctionalAffect(; affect...) - -function Base.:(==)(e1::SymbolicContinuousCallback, e2::SymbolicContinuousCallback) - isequal(e1.eqs, e2.eqs) && isequal(e1.affect, e2.affect) && - isequal(e1.initialize, e2.initialize) && isequal(e1.finalize, e2.finalize) && - isequal(e1.affect_neg, e2.affect_neg) && isequal(e1.rootfind, e2.rootfind) + +function SymbolicContinuousCallback(p::Pair, args...; kwargs...) + SymbolicContinuousCallback(p[1], p[2], args...; kwargs...) end -Base.isempty(cb::SymbolicContinuousCallback) = isempty(cb.eqs) -function Base.hash(cb::SymbolicContinuousCallback, s::UInt) - hash_affect(affect::AbstractVector, s) = foldr(hash, affect, init = s) - hash_affect(affect, s) = hash(affect, s) - s = foldr(hash, cb.eqs, init = s) - s = hash_affect(cb.affect, s) - s = hash_affect(cb.affect_neg, s) - s = hash_affect(cb.initialize, s) - s = hash_affect(cb.finalize, s) - s = hash(cb.reinitializealg, s) - hash(cb.rootfind, s) +SymbolicContinuousCallback(cb::SymbolicContinuousCallback, args...; kwargs...) = cb +SymbolicContinuousCallback(cb::Nothing, args...; kwargs...) = nothing +function SymbolicContinuousCallback(cb::Tuple, args...; kwargs...) + if length(cb) == 2 + SymbolicContinuousCallback(cb[1]; kwargs..., cb[2]...) + else + error("Malformed tuple specifying callback. Should be a condition => affect pair, followed by a vector of kwargs.") + end end -function Base.show(io::IO, cb::SymbolicContinuousCallback) +make_affect(affect::Nothing; kwargs...) = nothing +make_affect(affect::Tuple; kwargs...) = FunctionalAffect(affect...) +make_affect(affect::NamedTuple; kwargs...) = FunctionalAffect(; affect...) +make_affect(affect::Affect; kwargs...) = affect + +function make_affect(affect::Vector{Equation}; discrete_parameters = Any[], + iv = nothing, alg_eqs::Vector{Equation} = Equation[], warn_no_algebraic = true, kwargs...) + isempty(affect) && return nothing + if isnothing(iv) + iv = t_nounits + @warn "No independent variable specified. Defaulting to t_nounits." + end + + discrete_parameters isa AbstractVector || (discrete_parameters = [discrete_parameters]) + for p in discrete_parameters + occursin(unwrap(iv), unwrap(p)) || + error("Non-time dependent parameter $p passed in as a discrete. Must be declared as @parameters $p(t).") + end + + dvs = OrderedSet() + params = OrderedSet() + for eq in affect + if !haspre(eq) && !(symbolic_type(eq.rhs) === NotSymbolic() || + symbolic_type(eq.lhs) === NotSymbolic()) + @warn "Affect equation $eq has no `Pre` operator. As such it will be interpreted as an algebraic equation to be satisfied after the callback. If you intended to use the value of a variable x before the affect, use Pre(x). Errors may be thrown if there is no `Pre` and the algebraic equation is unsatisfiable, such as X ~ X + 1." + end + collect_vars!(dvs, params, eq, iv; op = Pre) + diffvs = collect_applied_operators(eq, Differential) + union!(dvs, diffvs) + end + for eq in alg_eqs + collect_vars!(dvs, params, eq, iv) + end + + pre_params = filter(haspre ∘ value, params) + sys_params = collect(setdiff(params, union(discrete_parameters, pre_params))) + discretes = map(tovar, discrete_parameters) + dvs = collect(dvs) + _dvs = map(default_toterm, dvs) + + aff_map = Dict(zip(discretes, discrete_parameters)) + rev_map = Dict(zip(discrete_parameters, discretes)) + subs = merge(rev_map, Dict(zip(dvs, _dvs))) + affect = Symbolics.fast_substitute(affect, subs) + alg_eqs = Symbolics.fast_substitute(alg_eqs, subs) + + @named affectsys = ImplicitDiscreteSystem( + vcat(affect, alg_eqs), iv, collect(union(_dvs, discretes)), + collect(union(pre_params, sys_params))) + affectsys = structural_simplify(affectsys; fully_determined = nothing) + # get accessed parameters p from Pre(p) in the callback parameters + accessed_params = filter(isparameter, map(unPre, collect(pre_params))) + union!(accessed_params, sys_params) + + # add scalarized unknowns to the map. + _dvs = reduce(vcat, map(scalarize, _dvs), init = Any[]) + for u in _dvs + aff_map[u] = u + end + + AffectSystem(affectsys, collect(_dvs), collect(accessed_params), + collect(discrete_parameters), aff_map) +end + +function make_affect(affect; kwargs...) + error("Malformed affect $(affect). This should be a vector of equations or a tuple specifying a functional affect.") +end + +function Base.show(io::IO, cb::AbstractCallback) indent = get(io, :indent, 0) iio = IOContext(io, :indent => indent + 1) + is_discrete(cb) ? print(io, "SymbolicDiscreteCallback(") : print(io, "SymbolicContinuousCallback(") - print(iio, "Equations:") + print(iio, "Conditions:") show(iio, equations(cb)) print(iio, "; ") - if affects(cb) != NULL_AFFECT + if affects(cb) != nothing print(iio, "Affect:") show(iio, affects(cb)) print(iio, ", ") end - if affect_negs(cb) != NULL_AFFECT + if !is_discrete(cb) && affect_negs(cb) != nothing print(iio, "Negative-edge affect:") show(iio, affect_negs(cb)) print(iio, ", ") end - if initialize_affects(cb) != NULL_AFFECT + if initialize_affects(cb) != nothing print(iio, "Initialization affect:") show(iio, initialize_affects(cb)) print(iio, ", ") end - if finalize_affects(cb) != NULL_AFFECT + if finalize_affects(cb) != nothing print(iio, "Finalization affect:") show(iio, finalize_affects(cb)) end print(iio, ")") end -function Base.show(io::IO, mime::MIME"text/plain", cb::SymbolicContinuousCallback) +function Base.show(io::IO, mime::MIME"text/plain", cb::AbstractCallback) indent = get(io, :indent, 0) iio = IOContext(io, :indent => indent + 1) + is_discrete(cb) ? println(io, "SymbolicDiscreteCallback:") : println(io, "SymbolicContinuousCallback:") - println(iio, "Equations:") + println(iio, "Conditions:") show(iio, mime, equations(cb)) print(iio, "\n") - if affects(cb) != NULL_AFFECT + if affects(cb) != nothing println(iio, "Affect:") show(iio, mime, affects(cb)) print(iio, "\n") end - if affect_negs(cb) != NULL_AFFECT - println(iio, "Negative-edge affect:") + if !is_discrete(cb) && affect_negs(cb) != nothing + print(iio, "Negative-edge affect:\n") show(iio, mime, affect_negs(cb)) print(iio, "\n") end - if initialize_affects(cb) != NULL_AFFECT + if initialize_affects(cb) != nothing println(iio, "Initialization affect:") show(iio, mime, initialize_affects(cb)) print(iio, "\n") end - if finalize_affects(cb) != NULL_AFFECT + if finalize_affects(cb) != nothing println(iio, "Finalization affect:") show(iio, mime, finalize_affects(cb)) print(iio, "\n") end end -to_equation_vector(eq::Equation) = [eq] -to_equation_vector(eqs::Vector{Equation}) = eqs -function to_equation_vector(eqs::Vector{Any}) - isempty(eqs) || error("This should never happen") - Equation[] -end - -function SymbolicContinuousCallback(args...) - SymbolicContinuousCallback(to_equation_vector.(args)...) -end # wrap eq in vector -SymbolicContinuousCallback(p::Pair) = SymbolicContinuousCallback(p[1], p[2]) -SymbolicContinuousCallback(cb::SymbolicContinuousCallback) = cb # passthrough -function SymbolicContinuousCallback(eqs::Equation, affect = NULL_AFFECT; - initialize = NULL_AFFECT, finalize = NULL_AFFECT, - affect_neg = affect, rootfind = SciMLBase.LeftRootFind) - SymbolicContinuousCallback( - eqs = [eqs], affect = affect, affect_neg = affect_neg, - initialize = initialize, finalize = finalize, rootfind = rootfind) -end -function SymbolicContinuousCallback(eqs::Vector{Equation}, affect = NULL_AFFECT; - affect_neg = affect, initialize = NULL_AFFECT, finalize = NULL_AFFECT, - rootfind = SciMLBase.LeftRootFind) - SymbolicContinuousCallback( - eqs = eqs, affect = affect, affect_neg = affect_neg, - initialize = initialize, finalize = finalize, rootfind = rootfind) -end - -SymbolicContinuousCallbacks(cb::SymbolicContinuousCallback) = [cb] -SymbolicContinuousCallbacks(cbs::Vector{<:SymbolicContinuousCallback}) = cbs -SymbolicContinuousCallbacks(cbs::Vector) = SymbolicContinuousCallback.(cbs) -function SymbolicContinuousCallbacks(ve::Vector{Equation}) - SymbolicContinuousCallbacks(SymbolicContinuousCallback(ve)) -end -function SymbolicContinuousCallbacks(others) - SymbolicContinuousCallbacks(SymbolicContinuousCallback(others)) -end -SymbolicContinuousCallbacks(::Nothing) = SymbolicContinuousCallback[] - -equations(cb::SymbolicContinuousCallback) = cb.eqs -function equations(cbs::Vector{<:SymbolicContinuousCallback}) - mapreduce(equations, vcat, cbs, init = Equation[]) -end - -affects(cb::SymbolicContinuousCallback) = cb.affect -function affects(cbs::Vector{SymbolicContinuousCallback}) - mapreduce(affects, vcat, cbs, init = Equation[]) -end - -affect_negs(cb::SymbolicContinuousCallback) = cb.affect_neg -function affect_negs(cbs::Vector{SymbolicContinuousCallback}) - mapreduce(affect_negs, vcat, cbs, init = Equation[]) -end - -reinitialization_alg(cb::SymbolicContinuousCallback) = cb.reinitializealg -function reinitialization_algs(cbs::Vector{SymbolicContinuousCallback}) - mapreduce( - reinitialization_alg, vcat, cbs, init = SciMLBase.DAEInitializationAlgorithm[]) -end - -initialize_affects(cb::SymbolicContinuousCallback) = cb.initialize -function initialize_affects(cbs::Vector{SymbolicContinuousCallback}) - mapreduce(initialize_affects, vcat, cbs, init = Equation[]) -end - -finalize_affects(cb::SymbolicContinuousCallback) = cb.finalize -function finalize_affects(cbs::Vector{SymbolicContinuousCallback}) - mapreduce(finalize_affects, vcat, cbs, init = Equation[]) -end - -namespace_affects(af::Vector, s) = Equation[namespace_affect(a, s) for a in af] -namespace_affects(af::FunctionalAffect, s) = namespace_affect(af, s) -namespace_affects(::Nothing, s) = nothing - -function namespace_callback(cb::SymbolicContinuousCallback, s)::SymbolicContinuousCallback - SymbolicContinuousCallback(; - eqs = namespace_equation.(equations(cb), (s,)), - affect = namespace_affects(affects(cb), s), - affect_neg = namespace_affects(affect_negs(cb), s), - initialize = namespace_affects(initialize_affects(cb), s), - finalize = namespace_affects(finalize_affects(cb), s), - rootfind = cb.rootfind) -end - -""" - continuous_events(sys::AbstractSystem)::Vector{SymbolicContinuousCallback} - -Returns a vector of all the `continuous_events` in an abstract system and its component subsystems. -The `SymbolicContinuousCallback`s in the returned vector are structs with two fields: `eqs` and -`affect` which correspond to the first and second elements of a `Pair` used to define an event, i.e. -`eqs => affect`. -""" -function continuous_events(sys::AbstractSystem) - cbs = get_continuous_events(sys) - filter(!isempty, cbs) - - systems = get_systems(sys) - cbs = [cbs; - reduce(vcat, - (map(cb -> namespace_callback(cb, s), continuous_events(s)) - for s in systems), - init = SymbolicContinuousCallback[])] - filter(!isempty, cbs) -end - -function vars!(vars, cb::SymbolicContinuousCallback; op = Differential) - for eq in equations(cb) - vars!(vars, eq; op) - end - for aff in (affects(cb), affect_negs(cb), initialize_affects(cb), finalize_affects(cb)) - if aff isa Vector{Equation} - for eq in aff +function vars!(vars, cb::AbstractCallback; op = Differential) + if symbolic_type(conditions(cb)) == NotSymbolic + if conditions(cb) isa AbstractArray + for eq in conditions(cb) vars!(vars, eq; op) end - elseif aff !== nothing - vars!(vars, aff; op) end + else + vars!(vars, conditions(cb); op) + end + for aff in (affects(cb), initialize_affects(cb), finalize_affects(cb)) + isnothing(aff) || vars!(vars, aff; op) end + !is_discrete(cb) && vars!(vars, affect_negs(cb); op) return vars end +################################ +######## Discrete events ####### +################################ """ - continuous_events_toplevel(sys::AbstractSystem) + SymbolicDiscreteCallback(conditions::Vector{Equation}, affect = nothing, iv = nothing; + initialize = nothing, finalize = nothing, alg_eqs = Equation[]) -Replicates the behaviour of `continuous_events`, but ignores events of subsystems. +A callback that triggers at the first timestep that the conditions are satisfied. -Notes: -- Cannot be applied to non-complete systems. -""" -function continuous_events_toplevel(sys::AbstractSystem) - if has_parent(sys) && (parent = get_parent(sys)) !== nothing - return continuous_events_toplevel(parent) - end - return get_continuous_events(sys) -end +The condition can be one of: +- Δt::Real - periodic events with period Δt +- ts::Vector{Real} - events trigger at these preset times given by `ts` +- eqs::Vector{Symbolic} - events trigger when the condition evaluates to true -#################################### discrete events ##################################### - -struct SymbolicDiscreteCallback - # condition can be one of: - # Δt::Real - Periodic with period Δt - # Δts::Vector{Real} - events trigger in this times (Preset) - # condition::Vector{Equation} - event triggered when condition is true - # TODO: Iterative - condition::Any - affects::Any - initialize::Any - finalize::Any +Arguments: +- iv: The independent variable of the system. This must be specified if the independent variable appears in one of the equations explicitly, as in x ~ t + 1. +- alg_eqs: Algebraic equations of the system that must be satisfied after the callback occurs. +""" +struct SymbolicDiscreteCallback <: AbstractCallback + conditions::Union{Number, Vector{<:Number}, Symbolic{Bool}} + affect::Union{Affect, Nothing} + initialize::Union{Affect, Nothing} + finalize::Union{Affect, Nothing} reinitializealg::SciMLBase.DAEInitializationAlgorithm function SymbolicDiscreteCallback( - condition, affects = NULL_AFFECT; reinitializealg = SciMLBase.CheckInit(), - initialize = NULL_AFFECT, finalize = NULL_AFFECT) - c = scalarize_condition(condition) - a = scalarize_affects(affects) - new(c, a, scalarize_affects(initialize), - scalarize_affects(finalize), reinitializealg) + condition::Union{Symbolic{Bool}, Number, Vector{<:Number}}, affect = nothing; + initialize = nothing, finalize = nothing, + reinitializealg = nothing, kwargs...) + c = is_timed_condition(condition) ? condition : value(scalarize(condition)) + + if isnothing(reinitializealg) + any(a -> (a isa FunctionalAffect || a isa ImperativeAffect), + [affect, initialize, finalize]) ? + reinitializealg = SciMLBase.CheckInit() : + reinitializealg = SciMLBase.NoInit() + end + new(c, make_affect(affect; kwargs...), + make_affect(initialize; kwargs...), + make_affect(finalize; kwargs...), reinitializealg) end # Default affect to nothing end -is_timed_condition(cb) = false -is_timed_condition(::R) where {R <: Real} = true -is_timed_condition(::V) where {V <: AbstractVector} = eltype(V) <: Real -is_timed_condition(::Num) = false -is_timed_condition(cb::SymbolicDiscreteCallback) = is_timed_condition(condition(cb)) - -function scalarize_condition(condition) - is_timed_condition(condition) ? condition : value(scalarize(condition)) +function SymbolicDiscreteCallback(p::Pair, args...; kwargs...) + SymbolicDiscreteCallback(p[1], p[2], args...; kwargs...) end -function namespace_condition(condition, s) - is_timed_condition(condition) ? condition : namespace_expr(condition, s) -end - -scalarize_affects(affects) = scalarize(affects) -scalarize_affects(affects::Tuple) = FunctionalAffect(affects...) -scalarize_affects(affects::NamedTuple) = FunctionalAffect(; affects...) -scalarize_affects(affects::FunctionalAffect) = affects - -SymbolicDiscreteCallback(p::Pair) = SymbolicDiscreteCallback(p[1], p[2]) -SymbolicDiscreteCallback(cb::SymbolicDiscreteCallback) = cb # passthrough - -function Base.show(io::IO, db::SymbolicDiscreteCallback) - println(io, "condition: ", db.condition) - println(io, "affects:") - if db.affects isa FunctionalAffect || db.affects isa ImperativeAffect - # TODO - println(io, " ", db.affects) +SymbolicDiscreteCallback(cb::SymbolicDiscreteCallback, args...; kwargs...) = cb +SymbolicDiscreteCallback(cb::Nothing, args...; kwargs...) = nothing +function SymbolicDiscreteCallback(cb::Tuple, args...; kwargs...) + if length(cb) == 2 + SymbolicDiscreteCallback(cb[1]; cb[2]...) else - for affect in db.affects - println(io, " ", affect) - end + error("Malformed tuple specifying callback. Should be a condition => affect pair, followed by a vector of kwargs.") end end -function Base.:(==)(e1::SymbolicDiscreteCallback, e2::SymbolicDiscreteCallback) - isequal(e1.condition, e2.condition) && isequal(e1.affects, e2.affects) && - isequal(e1.initialize, e2.initialize) && isequal(e1.finalize, e2.finalize) -end -function Base.hash(cb::SymbolicDiscreteCallback, s::UInt) - s = hash(cb.condition, s) - s = cb.affects isa AbstractVector ? foldr(hash, cb.affects, init = s) : - hash(cb.affects, s) - s = cb.initialize isa AbstractVector ? foldr(hash, cb.initialize, init = s) : - hash(cb.initialize, s) - s = cb.finalize isa AbstractVector ? foldr(hash, cb.finalize, init = s) : - hash(cb.finalize, s) - s = hash(cb.reinitializealg, s) - return s +function is_timed_condition(condition::T) where {T} + if T === Num + false + elseif T <: Real + true + elseif T <: AbstractVector + eltype(condition) <: Real + else + false + end end -condition(cb::SymbolicDiscreteCallback) = cb.condition -function conditions(cbs::Vector{<:SymbolicDiscreteCallback}) - reduce(vcat, condition(cb) for cb in cbs) +to_cb_vector(cbs::Vector{<:AbstractCallback}; kwargs...) = cbs +to_cb_vector(cbs::Union{Nothing, Vector{Nothing}}; kwargs...) = AbstractCallback[] +to_cb_vector(cb::AbstractCallback; kwargs...) = [cb] +function to_cb_vector(cbs; CB_TYPE = SymbolicContinuousCallback, kwargs...) + if cbs isa Pair + [CB_TYPE(cbs; kwargs...)] + else + Vector{CB_TYPE}([CB_TYPE(cb; kwargs...) for cb in cbs]) + end end -affects(cb::SymbolicDiscreteCallback) = cb.affects - -function affects(cbs::Vector{SymbolicDiscreteCallback}) - reduce(vcat, affects(cb) for cb in cbs; init = []) +############################################ +########## Namespacing Utilities ########### +############################################ +function namespace_affects(affect::FunctionalAffect, s) + FunctionalAffect(func(affect), + renamespace.((s,), unknowns(affect)), + unknowns_syms(affect), + renamespace.((s,), parameters(affect)), + parameters_syms(affect), + renamespace.((s,), discretes(affect)), + context(affect)) end -reinitialization_alg(cb::SymbolicDiscreteCallback) = cb.reinitializealg -function reinitialization_algs(cbs::Vector{SymbolicDiscreteCallback}) - mapreduce( - reinitialization_alg, vcat, cbs, init = SciMLBase.DAEInitializationAlgorithm[]) +function namespace_affects(affect::AffectSystem, s) + AffectSystem(renamespace(s, system(affect)), + renamespace.((s,), unknowns(affect)), + renamespace.((s,), parameters(affect)), + renamespace.((s,), discretes(affect)), + Dict([k => renamespace(s, v) for (k, v) in aff_to_sys(affect)])) end +namespace_affects(af::Nothing, s) = nothing -initialize_affects(cb::SymbolicDiscreteCallback) = cb.initialize -function initialize_affects(cbs::Vector{SymbolicDiscreteCallback}) - mapreduce(initialize_affects, vcat, cbs, init = Equation[]) +function namespace_callback(cb::SymbolicContinuousCallback, s)::SymbolicContinuousCallback + SymbolicContinuousCallback( + namespace_equation.(equations(cb), (s,)), + namespace_affects(affects(cb), s), + affect_neg = namespace_affects(affect_negs(cb), s), + initialize = namespace_affects(initialize_affects(cb), s), + finalize = namespace_affects(finalize_affects(cb), s), + rootfind = cb.rootfind, reinitializealg = cb.reinitializealg) end -finalize_affects(cb::SymbolicDiscreteCallback) = cb.finalize -function finalize_affects(cbs::Vector{SymbolicDiscreteCallback}) - mapreduce(finalize_affects, vcat, cbs, init = Equation[]) +function namespace_conditions(condition, s) + is_timed_condition(condition) ? condition : namespace_expr(condition, s) end function namespace_callback(cb::SymbolicDiscreteCallback, s)::SymbolicDiscreteCallback - function namespace_affects(af) - return af isa AbstractVector ? namespace_affect.(af, Ref(s)) : - namespace_affect(af, s) - end SymbolicDiscreteCallback( - namespace_condition(condition(cb), s), namespace_affects(affects(cb)), - reinitializealg = cb.reinitializealg, initialize = namespace_affects(initialize_affects(cb)), - finalize = namespace_affects(finalize_affects(cb))) + namespace_conditions(conditions(cb), s), + namespace_affects(affects(cb), s), + initialize = namespace_affects(initialize_affects(cb), s), + finalize = namespace_affects(finalize_affects(cb), s), reinitializealg = cb.reinitializealg) end -SymbolicDiscreteCallbacks(cb::Pair) = SymbolicDiscreteCallback[SymbolicDiscreteCallback(cb)] -SymbolicDiscreteCallbacks(cbs::Vector) = SymbolicDiscreteCallback.(cbs) -SymbolicDiscreteCallbacks(cb::SymbolicDiscreteCallback) = [cb] -SymbolicDiscreteCallbacks(cbs::Vector{<:SymbolicDiscreteCallback}) = cbs -SymbolicDiscreteCallbacks(::Nothing) = SymbolicDiscreteCallback[] +function Base.hash(cb::AbstractCallback, s::UInt) + s = conditions(cb) isa AbstractVector ? foldr(hash, conditions(cb), init = s) : + hash(conditions(cb), s) + s = hash(affects(cb), s) + !is_discrete(cb) && (s = hash(affect_negs(cb), s)) + s = hash(initialize_affects(cb), s) + s = hash(finalize_affects(cb), s) + !is_discrete(cb) && (s = hash(cb.rootfind, s)) + hash(cb.reinitializealg, s) +end -""" - discrete_events(sys::AbstractSystem) :: Vector{SymbolicDiscreteCallback} +########################### +######### Helpers ######### +########################### -Returns a vector of all the `discrete_events` in an abstract system and its component subsystems. -The `SymbolicDiscreteCallback`s in the returned vector are structs with two fields: `condition` and -`affect` which correspond to the first and second elements of a `Pair` used to define an event, i.e. -`condition => affect`. -""" -function discrete_events(sys::AbstractSystem) - cbs = get_discrete_events(sys) - systems = get_systems(sys) - cbs = [cbs; - reduce(vcat, - (map(cb -> namespace_callback(cb, s), discrete_events(s)) for s in systems), - init = SymbolicDiscreteCallback[])] - cbs +conditions(cb::AbstractCallback) = cb.conditions +function conditions(cbs::Vector{<:AbstractCallback}) + reduce(vcat, conditions(cb) for cb in cbs; init = []) end +equations(cb::AbstractCallback) = conditions(cb) +equations(cb::Vector{<:AbstractCallback}) = conditions(cb) -function vars!(vars, cb::SymbolicDiscreteCallback; op = Differential) - if symbolic_type(cb.condition) == NotSymbolic - if cb.condition isa AbstractArray - for eq in cb.condition - vars!(vars, eq; op) - end - end - else - vars!(vars, cb.condition; op) - end - for aff in (cb.affects, cb.initialize, cb.finalize) - if aff isa Vector{Equation} - for eq in aff - vars!(vars, eq; op) - end - elseif aff !== nothing - vars!(vars, aff; op) - end - end - return vars +affects(cb::AbstractCallback) = cb.affect +function affects(cbs::Vector{<:AbstractCallback}) + reduce(vcat, affects(cb) for cb in cbs; init = []) end -""" - discrete_events_toplevel(sys::AbstractSystem) - -Replicates the behaviour of `discrete_events`, but ignores events of subsystems. - -Notes: -- Cannot be applied to non-complete systems. -""" -function discrete_events_toplevel(sys::AbstractSystem) - if has_parent(sys) && (parent = get_parent(sys)) !== nothing - return discrete_events_toplevel(parent) - end - return get_discrete_events(sys) +affect_negs(cb::SymbolicContinuousCallback) = cb.affect_neg +function affect_negs(cbs::Vector{SymbolicContinuousCallback}) + reduce(vcat, affect_negs(cb) for cb in cbs; init = []) end -################################# compilation functions #################################### - -# handles ensuring that affect! functions work with integrator arguments -function add_integrator_header( - sys::AbstractSystem, integrator = gensym(:MTKIntegrator), out = :u) - expr -> Func([DestructuredArgs(expr.args, integrator, inds = [:u, :p, :t])], [], - expr.body), - expr -> Func( - [DestructuredArgs(expr.args, integrator, inds = [out, :u, :p, :t])], [], - expr.body) +initialize_affects(cb::AbstractCallback) = cb.initialize +function initialize_affects(cbs::Vector{<:AbstractCallback}) + reduce(initialize_affects, vcat, cbs; init = []) end -function condition_header(sys::AbstractSystem, integrator = gensym(:MTKIntegrator)) - expr -> Func( - [expr.args[1], expr.args[2], - DestructuredArgs(expr.args[3:end], integrator, inds = [:p])], - [], - expr.body) +finalize_affects(cb::AbstractCallback) = cb.finalize +function finalize_affects(cbs::Vector{<:AbstractCallback}) + reduce(finalize_affects, vcat, cbs; init = []) end -function callback_save_header(sys::AbstractSystem, cb) - if !(has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing) - return (identity, identity) - end - save_idxs = get(ic.callback_to_clocks, cb, Int[]) - isempty(save_idxs) && return (identity, identity) - - wrapper = function (expr) - return Func(expr.args, [], - LiteralExpr(quote - $(expr.body) - save_idxs = $(save_idxs) - for idx in save_idxs - $(SciMLBase.save_discretes!)($(expr.args[1]), idx) - end - end)) - end - - return wrapper, wrapper +function Base.:(==)(e1::AbstractCallback, e2::AbstractCallback) + (is_discrete(e1) === is_discrete(e2)) || return false + (isequal(e1.conditions, e2.conditions) && isequal(e1.affect, e2.affect) && + isequal(e1.initialize, e2.initialize) && isequal(e1.finalize, e2.finalize)) && + isequal(e1.reinitializealg, e2.reinitializealg) || + return false + is_discrete(e1) || + (isequal(e1.affect_neg, e2.affect_neg) && isequal(e1.rootfind, e2.rootfind)) end -""" - compile_condition(cb::SymbolicDiscreteCallback, sys, dvs, ps; expression, kwargs...) - -Returns a function `condition(u,t,integrator)` returning the `condition(cb)`. +Base.isempty(cb::AbstractCallback) = isempty(cb.conditions) -Notes +#################################### +####### Compilation functions ###### +#################################### +""" + compile_condition(cb::AbstractCallback, sys, dvs, ps; expression, kwargs...) - - `expression = Val{true}`, causes the generated function to be returned as an expression. - If set to `Val{false}` a `RuntimeGeneratedFunction` will be returned. - - `kwargs` are passed through to `Symbolics.build_function`. +Returns a function `condition(u,t,integrator)`, condition(out,u,t,integrator)` returning the `condition(cb)`. """ -function compile_condition(cb::SymbolicDiscreteCallback, sys, dvs, ps; - expression = Val{true}, eval_expression = false, eval_module = @__MODULE__, kwargs...) - u = map(x -> time_varying_as_func(value(x), sys), dvs) - p = map.(x -> time_varying_as_func(value(x), sys), reorder_parameters(sys, ps)) +function compile_condition(cbs::Union{AbstractCallback, Vector{<:AbstractCallback}}, sys, dvs, ps; + eval_expression = false, eval_module = @__MODULE__, kwargs...) + u = map(value, dvs) + p = map.(value, reorder_parameters(sys, ps)) t = get_iv(sys) - condit = condition(cb) + condit = conditions(cbs) cs = collect_constants(condit) if !isempty(cs) cmap = map(x -> x => getdefault(x), cs) - condit = substitute(condit, cmap) - end - expr = build_function_wrapper(sys, - condit, u, t, p...; expression = Val{true}, - p_start = 3, p_end = length(p) + 2, - wrap_code = condition_header(sys), - kwargs...) - if expression == Val{true} - return expr + condit = substitute(condit, Dict(cmap)) end - return eval_or_rgf(expr; eval_expression, eval_module) -end - -function compile_affect(cb::SymbolicContinuousCallback, args...; kwargs...) - compile_affect(affects(cb), cb, args...; kwargs...) -end -""" - compile_affect(eqs::Vector{Equation}, sys, dvs, ps; expression, outputidxs, kwargs...) - compile_affect(cb::SymbolicContinuousCallback, args...; kwargs...) + if !is_discrete(cbs) + condit = reduce(vcat, flatten_equations(condit)) + condit = condit isa AbstractVector ? [c.lhs - c.rhs for c in condit] : + [condit.lhs - condit.rhs] + end -Returns a function that takes an integrator as argument and modifies the state with the -affect. The generated function has the signature `affect!(integrator)`. + fs = build_function_wrapper(sys, condit, u, p..., t; kwargs..., expression = Val{false}, cse = false) + (f_oop, f_iip) = is_discrete(cbs) ? (fs, nothing) : fs -Notes - - - `expression = Val{true}`, causes the generated function to be returned as an expression. - If set to `Val{false}` a `RuntimeGeneratedFunction` will be returned. - - `outputidxs`, a vector of indices of the output variables which should correspond to - `unknowns(sys)`. If provided, checks that the LHS of affect equations are variables are - dropped, i.e. it is assumed these indices are correct and affect equations are - well-formed. - - `kwargs` are passed through to `Symbolics.build_function`. -""" -function compile_affect(eqs::Vector{Equation}, cb, sys, dvs, ps; outputidxs = nothing, - expression = Val{true}, checkvars = true, eval_expression = false, - eval_module = @__MODULE__, - postprocess_affect_expr! = nothing, kwargs...) - if isempty(eqs) - if expression == Val{true} - return :((args...) -> ()) - else - return (args...) -> () # We don't do anything in the callback, we're just after the event - end + cond = if cbs isa AbstractVector + (out, u, t, integ) -> f_iip(out, u, parameter_values(integ), t) + elseif is_discrete(cbs) + (u, t, integ) -> f_oop(u, parameter_values(integ), t) else - eqs = flatten_equations(eqs) - rhss = map(x -> x.rhs, eqs) - outvar = :u - if outputidxs === nothing - lhss = map(x -> x.lhs, eqs) - all(isvariable, lhss) || - error("Non-variable symbolic expression found on the left hand side of an affect equation. Such equations must be of the form variable ~ symbolic expression for the new value of the variable.") - update_vars = collect(Iterators.flatten(map(ModelingToolkit.vars, lhss))) # these are the ones we're changing - length(update_vars) == length(unique(update_vars)) == length(eqs) || - error("affected variables not unique, each unknown can only be affected by one equation for a single `root_eqs => affects` pair.") - alleq = all(isequal(isparameter(first(update_vars))), - Iterators.map(isparameter, update_vars)) - if !isparameter(first(lhss)) && alleq - unknownind = Dict(reverse(en) for en in enumerate(dvs)) - update_inds = map(sym -> unknownind[sym], update_vars) - elseif isparameter(first(lhss)) && alleq - if has_index_cache(sys) && get_index_cache(sys) !== nothing - update_inds = map(update_vars) do sym - return parameter_index(sys, sym) - end - else - psind = Dict(reverse(en) for en in enumerate(ps)) - update_inds = map(sym -> psind[sym], update_vars) - end - outvar = :p + function (u, t, integ) + if DiffEqBase.isinplace(SciMLBase.get_sol(integ).prob) + tmp, = DiffEqBase.get_tmp_cache(integ) + f_iip(tmp, u, parameter_values(integ), t) + tmp[1] else - error("Error, building an affect function for a callback that wants to modify both parameters and unknowns. This is not currently allowed in one individual callback.") + f_oop(u, parameter_values(integ), t) end - else - update_inds = outputidxs end - - _ps = ps - ps = reorder_parameters(sys, ps) - if checkvars - u = map(x -> time_varying_as_func(value(x), sys), dvs) - p = map.(x -> time_varying_as_func(value(x), sys), ps) - else - u = dvs - p = ps - end - t = get_iv(sys) - integ = gensym(:MTKIntegrator) - rf_oop, rf_ip = build_function_wrapper( - sys, rhss, u, p..., t; expression = Val{true}, - wrap_code = callback_save_header(sys, cb) .∘ - add_integrator_header(sys, integ, outvar), - outputidxs = update_inds, - create_bindings = false, - kwargs..., cse = false) - # applied user-provided function to the generated expression - if postprocess_affect_expr! !== nothing - postprocess_affect_expr!(rf_ip, integ) - end - if expression == Val{false} - return eval_or_rgf(rf_ip; eval_expression, eval_module) - end - return rf_ip end end -function generate_rootfinding_callback(sys::AbstractTimeDependentSystem, - dvs = unknowns(sys), ps = parameters(sys; initial_parameters = true); kwargs...) - cbs = continuous_events(sys) - isempty(cbs) && return nothing - generate_rootfinding_callback(cbs, sys, dvs, ps; kwargs...) -end """ -Generate a single rootfinding callback; this happens if there is only one equation in `cbs` passed to -generate_rootfinding_callback and thus we can produce a ContinuousCallback instead of a VectorContinuousCallback. +Compile user-defined functional affect. """ -function generate_single_rootfinding_callback( - eq, cb, sys::AbstractTimeDependentSystem, dvs = unknowns(sys), - ps = parameters(sys; initial_parameters = true); kwargs...) - if !isequal(eq.lhs, 0) - eq = 0 ~ eq.lhs - eq.rhs - end +function compile_functional_affect(affect::FunctionalAffect, sys; kwargs...) + dvs = unknowns(sys) + ps = parameters(sys) + dvs_ind = Dict(reverse(en) for en in enumerate(dvs)) + v_inds = map(sym -> dvs_ind[sym], unknowns(affect)) - rf_oop, rf_ip = generate_custom_function( - sys, [eq.rhs], dvs, ps; expression = Val{false}, kwargs..., cse = false) - affect_function = compile_affect_fn(cb, sys, dvs, ps, kwargs) - cond = function (u, t, integ) - if DiffEqBase.isinplace(integ.sol.prob) - tmp, = DiffEqBase.get_tmp_cache(integ) - rf_ip(tmp, u, parameter_values(integ), t) - tmp[1] - else - rf_oop(u, parameter_values(integ), t) - end - end - user_initfun = isnothing(affect_function.initialize) ? SciMLBase.INITIALIZE_DEFAULT : - (c, u, t, i) -> affect_function.initialize(i) - if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing && - (save_idxs = get(ic.callback_to_clocks, cb, nothing)) !== nothing - initfn = let save_idxs = save_idxs - function (cb, u, t, integrator) - user_initfun(cb, u, t, integrator) - for idx in save_idxs - SciMLBase.save_discretes!(integrator, idx) - end - end - end + if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing + p_inds = [(pind = parameter_index(sys, sym)) === nothing ? sym : pind + for sym in parameters(affect)] else - initfn = user_initfun + ps_ind = Dict(reverse(en) for en in enumerate(ps)) + p_inds = map(sym -> get(ps_ind, sym, sym), parameters(affect)) end + # HACK: filter out eliminated symbols. Not clear this is the right thing to do + # (MTK should keep these symbols) + u = filter(x -> !isnothing(x[2]), collect(zip(unknowns_syms(affect), v_inds))) |> + NamedTuple + p = filter(x -> !isnothing(x[2]), collect(zip(parameters_syms(affect), p_inds))) |> + NamedTuple - return ContinuousCallback( - cond, affect_function.affect, affect_function.affect_neg, rootfind = cb.rootfind, - initialize = initfn, - finalize = isnothing(affect_function.finalize) ? SciMLBase.FINALIZE_DEFAULT : - (c, u, t, i) -> affect_function.finalize(i), - initializealg = reinitialization_alg(cb)) + let u = u, p = p, user_affect = func(affect), ctx = context(affect) + (integ) -> begin + user_affect(integ, u, p, ctx) + end + end end -function generate_vector_rootfinding_callback( - cbs, sys::AbstractTimeDependentSystem, dvs = unknowns(sys), - ps = parameters(sys; initial_parameters = true); rootfind = SciMLBase.RightRootFind, - reinitialization = SciMLBase.CheckInit(), kwargs...) - eqs = map(cb -> flatten_equations(cb.eqs), cbs) - num_eqs = length.(eqs) - # fuse equations to create VectorContinuousCallback - eqs = reduce(vcat, eqs) - # rewrite all equations as 0 ~ interesting stuff - eqs = map(eqs) do eq - isequal(eq.lhs, 0) && return eq - 0 ~ eq.lhs - eq.rhs - end +is_discrete(cb::AbstractCallback) = cb isa SymbolicDiscreteCallback +is_discrete(cb::Vector{<:AbstractCallback}) = eltype(cb) isa SymbolicDiscreteCallback - rhss = map(x -> x.rhs, eqs) - _, rf_ip = generate_custom_function( - sys, rhss, dvs, ps; expression = Val{false}, kwargs..., cse = false) - - affect_functions = @NamedTuple{ - affect::Function, - affect_neg::Union{Function, Nothing}, - initialize::Union{Function, Nothing}, - finalize::Union{Function, Nothing}}[ - compile_affect_fn(cb, sys, dvs, ps, kwargs) - for cb in cbs] - cond = function (out, u, t, integ) - rf_ip(out, u, parameter_values(integ), t) - end +function generate_continuous_callbacks(sys::AbstractSystem, dvs = unknowns(sys), + ps = parameters(sys; initial_parameters = true); kwargs...) + cbs = continuous_events(sys) + isempty(cbs) && return nothing + cb_classes = Dict{Tuple{SciMLBase.RootfindOpt, SciMLBase.DAEInitializationAlgorithm}, + Vector{SymbolicContinuousCallback}}() - # since there may be different number of conditions and affects, - # we build a map that translates the condition eq. number to the affect number - eq_ind2affect = reduce(vcat, - [fill(i, num_eqs[i]) for i in eachindex(affect_functions)]) - @assert length(eq_ind2affect) == length(eqs) - @assert maximum(eq_ind2affect) == length(affect_functions) - - affect = let affect_functions = affect_functions, eq_ind2affect = eq_ind2affect - function (integ, eq_ind) # eq_ind refers to the equation index that triggered the event, each event has num_eqs[i] equations - affect_functions[eq_ind2affect[eq_ind]].affect(integ) - end - end - affect_neg = let affect_functions = affect_functions, eq_ind2affect = eq_ind2affect - function (integ, eq_ind) # eq_ind refers to the equation index that triggered the event, each event has num_eqs[i] equations - affect_neg = affect_functions[eq_ind2affect[eq_ind]].affect_neg - if isnothing(affect_neg) - return # skip if the neg function doesn't exist - don't want to split this into a separate VCC because that'd break ordering - end - affect_neg(integ) - end - end - function handle_optional_setup_fn(funs, default) - if all(isnothing, funs) - return default - else - return let funs = funs - function (cb, u, t, integ) - for func in funs - if isnothing(func) - continue - else - func(integ) - end - end - end - end - end + # Sort the callbacks by their rootfinding method + for cb in cbs + _cbs = get!(() -> SymbolicContinuousCallback[], + cb_classes, (cb.rootfind, cb.reinitializealg)) + push!(_cbs, cb) end - initialize = nothing - if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing - initialize = handle_optional_setup_fn( - map(cbs, affect_functions) do cb, fn - if (save_idxs = get(ic.callback_to_clocks, cb, nothing)) !== nothing - let save_idxs = save_idxs - custom_init = fn.initialize - (i) -> begin - !isnothing(custom_init) && custom_init(i) - for idx in save_idxs - SciMLBase.save_discretes!(i, idx) - end - end - end - else - fn.initialize - end - end, - SciMLBase.INITIALIZE_DEFAULT) - + sort!(OrderedDict(cb_classes), by = cb -> cb[1]) + compiled_callbacks = [generate_callback(cb, sys; kwargs...) + for ((rf, reinit), cb) in cb_classes] + if length(compiled_callbacks) == 1 + return only(compiled_callbacks) else - initialize = handle_optional_setup_fn( - map(fn -> fn.initialize, affect_functions), SciMLBase.INITIALIZE_DEFAULT) + return CallbackSet(compiled_callbacks...) end +end - finalize = handle_optional_setup_fn( - map(fn -> fn.finalize, affect_functions), SciMLBase.FINALIZE_DEFAULT) - return VectorContinuousCallback( - cond, affect, affect_neg, length(eqs), rootfind = rootfind, - initialize = initialize, finalize = finalize, initializealg = reinitialization) +function generate_discrete_callbacks(sys::AbstractSystem, dvs = unknowns(sys), + ps = parameters(sys; initial_parameters = true); kwargs...) + dbs = discrete_events(sys) + isempty(dbs) && return nothing + [generate_callback(db, sys; kwargs...) for db in dbs] end +EMPTY_AFFECT(args...) = nothing + """ -Compile a single continuous callback affect function(s). +Codegen a DifferentialEquations callback. A (set of) continuous callback with multiple equations becomes a VectorContinuousCallback. +Continuous callbacks with only one equation will become a ContinuousCallback. +Individual discrete callbacks become DiscreteCallback, PresetTimeCallback, PeriodicCallback depending on the case. """ -function compile_affect_fn(cb, sys::AbstractTimeDependentSystem, dvs, ps, kwargs) - eq_aff = affects(cb) - eq_neg_aff = affect_negs(cb) - affect = compile_affect(eq_aff, cb, sys, dvs, ps; expression = Val{false}, kwargs...) - function compile_optional_affect(aff, default = nothing) - if isnothing(aff) || aff == default - return nothing - else - return compile_affect(aff, cb, sys, dvs, ps; expression = Val{false}, kwargs...) - end - end - if eq_neg_aff === eq_aff - affect_neg = affect - else - affect_neg = _compile_optional_affect( - NULL_AFFECT, eq_neg_aff, cb, sys, dvs, ps; kwargs...) - end - initialize = _compile_optional_affect( - NULL_AFFECT, initialize_affects(cb), cb, sys, dvs, ps; kwargs...) - finalize = _compile_optional_affect( - NULL_AFFECT, finalize_affects(cb), cb, sys, dvs, ps; kwargs...) - (affect = affect, affect_neg = affect_neg, initialize = initialize, finalize = finalize) -end - -function generate_rootfinding_callback(cbs, sys::AbstractTimeDependentSystem, - dvs = unknowns(sys), ps = parameters(sys; initial_parameters = true); kwargs...) - eqs = map(cb -> flatten_equations(cb.eqs), cbs) +function generate_callback(cbs::Vector{SymbolicContinuousCallback}, sys; kwargs...) + eqs = map(cb -> flatten_equations(equations(cb)), cbs) num_eqs = length.(eqs) - total_eqs = sum(num_eqs) - (isempty(eqs) || total_eqs == 0) && return nothing - if total_eqs == 1 - # find the callback with only one eq + (isempty(eqs) || sum(num_eqs) == 0) && return nothing + if sum(num_eqs) == 1 cb_ind = findfirst(>(0), num_eqs) - if isnothing(cb_ind) - error("Inconsistent state in affect compilation; one equation but no callback with equations?") - end - cb = cbs[cb_ind] - return generate_single_rootfinding_callback(cb.eqs[], cb, sys, dvs, ps; kwargs...) + return generate_callback(cbs[cb_ind], sys; kwargs...) end - # group the cbs by what rootfind op they use - # groupby would be very useful here, but alas - cb_classes = Dict{ - @NamedTuple{ - rootfind::SciMLBase.RootfindOpt, - reinitialization::SciMLBase.DAEInitializationAlgorithm}, Vector{SymbolicContinuousCallback}}() + trigger = compile_condition( + cbs, sys, unknowns(sys), parameters(sys; initial_parameters = true); kwargs...) + affects = [] + affect_negs = [] + inits = [] + finals = [] for cb in cbs - push!( - get!(() -> SymbolicContinuousCallback[], cb_classes, - ( - rootfind = cb.rootfind, - reinitialization = reinitialization_alg(cb))), - cb) + affect = compile_affect(cb.affect, cb, sys; default = EMPTY_AFFECT, kwargs...) + push!(affects, affect) + affect_neg = (cb.affect_neg === cb.affect) ? affect : + compile_affect( + cb.affect_neg, cb, sys; default = EMPTY_AFFECT, kwargs...) + push!(affect_negs, affect_neg) + push!(inits, + compile_affect( + cb.initialize, cb, sys; default = nothing, is_init = true, kwargs...)) + push!(finals, compile_affect(cb.finalize, cb, sys; default = nothing, kwargs...)) end - # generate the callbacks out; we sort by the equivalence class to ensure a deterministic preference order - compiled_callbacks = map(collect(pairs(sort!( - OrderedDict(cb_classes); by = p -> p.rootfind)))) do (equiv_class, cbs_in_class) - return generate_vector_rootfinding_callback( - cbs_in_class, sys, dvs, ps; rootfind = equiv_class.rootfind, - reinitialization = equiv_class.reinitialization, kwargs...) + # Since there may be different number of conditions and affects, + # we build a map that translates the condition eq. number to the affect number + eq2affect = reduce(vcat, + [fill(i, num_eqs[i]) for i in eachindex(affects)]) + eqs = reduce(vcat, eqs) + + affect = let eq2affect = eq2affect, affects = affects + function (integ, idx) + affects[eq2affect[idx]](integ) + end end - if length(compiled_callbacks) == 1 - return compiled_callbacks[] - else - return CallbackSet(compiled_callbacks...) + affect_neg = let eq2affect = eq2affect, affect_negs = affect_negs + function (integ, idx) + f = affect_negs[eq2affect[idx]] + isnothing(f) && return + f(integ) + end end -end + initialize = wrap_vector_optional_affect(inits, SciMLBase.INITIALIZE_DEFAULT) + finalize = wrap_vector_optional_affect(finals, SciMLBase.FINALIZE_DEFAULT) -function compile_user_affect(affect::FunctionalAffect, cb, sys, dvs, ps; kwargs...) - dvs_ind = Dict(reverse(en) for en in enumerate(dvs)) - v_inds = map(sym -> dvs_ind[sym], unknowns(affect)) + return VectorContinuousCallback( + trigger, affect, affect_neg, length(eqs); initialize, finalize, + rootfind = cbs[1].rootfind, initializealg = cbs[1].reinitializealg) +end - if has_index_cache(sys) && get_index_cache(sys) !== nothing - p_inds = [if (pind = parameter_index(sys, sym)) === nothing - sym - else - pind - end - for sym in parameters(affect)] - else - ps_ind = Dict(reverse(en) for en in enumerate(ps)) - p_inds = map(sym -> get(ps_ind, sym, sym), parameters(affect)) - end - # HACK: filter out eliminated symbols. Not clear this is the right thing to do - # (MTK should keep these symbols) - u = filter(x -> !isnothing(x[2]), collect(zip(unknowns_syms(affect), v_inds))) |> - NamedTuple - p = filter(x -> !isnothing(x[2]), collect(zip(parameters_syms(affect), p_inds))) |> - NamedTuple +function generate_callback(cb, sys; kwargs...) + is_timed = is_timed_condition(conditions(cb)) + dvs = unknowns(sys) + ps = parameters(sys; initial_parameters = true) - if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing - save_idxs = get(ic.callback_to_clocks, cb, Int[]) + trigger = is_timed ? conditions(cb) : compile_condition(cb, sys, dvs, ps; kwargs...) + affect = compile_affect(cb.affect, cb, sys; default = EMPTY_AFFECT, kwargs...) + affect_neg = if is_discrete(cb) + nothing else - save_idxs = Int[] + (cb.affect === cb.affect_neg) ? affect : + compile_affect(cb.affect_neg, cb, sys; default = EMPTY_AFFECT, kwargs...) end - let u = u, p = p, user_affect = func(affect), ctx = context(affect), - save_idxs = save_idxs - - function (integ) - user_affect(integ, u, p, ctx) - for idx in save_idxs - SciMLBase.save_discretes!(integ, idx) - end + init = compile_affect(cb.initialize, cb, sys; default = SciMLBase.INITIALIZE_DEFAULT, + is_init = true, kwargs...) + final = compile_affect( + cb.finalize, cb, sys; default = SciMLBase.FINALIZE_DEFAULT, kwargs...) + + initialize = isnothing(cb.initialize) ? init : ((c, u, t, i) -> init(i)) + finalize = isnothing(cb.finalize) ? final : ((c, u, t, i) -> final(i)) + + if is_discrete(cb) + if is_timed && conditions(cb) isa AbstractVector + return PresetTimeCallback(trigger, affect; initialize, + finalize, initializealg = cb.reinitializealg) + elseif is_timed + return PeriodicCallback( + affect, trigger; initialize, finalize, initializealg = cb.reinitializealg) + else + return DiscreteCallback(trigger, affect; initialize, + finalize, initializealg = cb.reinitializealg) end + else + return ContinuousCallback(trigger, affect, affect_neg; initialize, finalize, + rootfind = cb.rootfind, initializealg = cb.reinitializealg) end end -function invalid_variables(sys, expr) - filter(x -> !any(isequal(x), all_symbols(sys)), reduce(vcat, vars(expr); init = [])) -end -function unassignable_variables(sys, expr) - assignable_syms = reduce( - vcat, Symbolics.scalarize.(vcat( - unknowns(sys), parameters(sys; initial_parameters = true))); - init = []) - written = reduce(vcat, Symbolics.scalarize.(vars(expr)); init = []) - return filter( - x -> !any(isequal(x), assignable_syms), written) -end +""" + compile_affect(cb::AbstractCallback, sys::AbstractSystem, dvs, ps; expression, outputidxs, kwargs...) -@generated function _generated_writeback(integ, setters::NamedTuple{NS1, <:Tuple}, - values::NamedTuple{NS2, <:Tuple}) where {NS1, NS2} - setter_exprs = [] - for name in NS2 - if !(name in NS1) - missing_name = "Tried to write back to $name from affect; only declared states ($NS1) may be written to." - error(missing_name) - end - push!(setter_exprs, :(setters.$name(integ, values.$name))) - end - return :(begin - $(setter_exprs...) - end) -end +Returns a function that takes an integrator as argument and modifies the state with the +affect. The generated function has the signature `affect!(integrator)`. -function check_assignable(sys, sym) - if symbolic_type(sym) == ScalarSymbolic() - is_variable(sys, sym) || is_parameter(sys, sym) - elseif symbolic_type(sym) == ArraySymbolic() - is_variable(sys, sym) || is_parameter(sys, sym) || - all(x -> check_assignable(sys, x), collect(sym)) - elseif sym isa Union{AbstractArray, Tuple} - all(x -> check_assignable(sys, x), sym) +Notes + - `kwargs` are passed through to `Symbolics.build_function`. +""" +function compile_affect( + aff::Union{Nothing, Affect}, cb::AbstractCallback, sys::AbstractSystem; + default = nothing, is_init = false, kwargs...) + save_idxs = if !(has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing) + Int[] else - false + get(ic.callback_to_clocks, cb, Int[]) end -end -function compile_affect(affect::FunctionalAffect, cb, sys, dvs, ps; kwargs...) - compile_user_affect(affect, cb, sys, dvs, ps; kwargs...) -end -function _compile_optional_affect(default, aff, cb, sys, dvs, ps; kwargs...) - if isnothing(aff) || aff == default - return nothing - else - return compile_affect(aff, cb, sys, dvs, ps; expression = Val{false}, kwargs...) + if isnothing(aff) + is_init ? wrap_save_discretes(default, save_idxs) : default + elseif aff isa AffectSystem + f = compile_equational_affect(aff, sys; kwargs...) + wrap_save_discretes(f, save_idxs) + elseif aff isa FunctionalAffect || aff isa ImperativeAffect + f = compile_functional_affect(aff, sys; kwargs...) + wrap_save_discretes(f, save_idxs) end end -function generate_timed_callback(cb, sys, dvs, ps; postprocess_affect_expr! = nothing, - kwargs...) - cond = condition(cb) - as = compile_affect(affects(cb), cb, sys, dvs, ps; expression = Val{false}, - postprocess_affect_expr!, kwargs...) - - user_initfun = _compile_optional_affect( - NULL_AFFECT, initialize_affects(cb), cb, sys, dvs, ps; kwargs...) - user_finfun = _compile_optional_affect( - NULL_AFFECT, finalize_affects(cb), cb, sys, dvs, ps; kwargs...) - if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing && - (save_idxs = get(ic.callback_to_clocks, cb, nothing)) !== nothing - initfn = let - save_idxs = save_idxs - initfun = user_initfun - function (cb, u, t, integrator) - if !isnothing(initfun) - initfun(integrator) + +function wrap_save_discretes(f, save_idxs) + let save_idxs = save_idxs, f = f + if f === SciMLBase.INITIALIZE_DEFAULT + (c, u, t, i) -> begin + f(c, u, t, i) + for idx in save_idxs + SciMLBase.save_discretes!(i, idx) end + end + else + (i) -> begin + isnothing(f) || f(i) for idx in save_idxs - SciMLBase.save_discretes!(integrator, idx) + SciMLBase.save_discretes!(i, idx) end end end - else - initfn = isnothing(user_initfun) ? SciMLBase.INITIALIZE_DEFAULT : - (_, _, _, i) -> user_initfun(i) - end - finfun = isnothing(user_finfun) ? SciMLBase.FINALIZE_DEFAULT : - (_, _, _, i) -> user_finfun(i) - if cond isa AbstractVector - # Preset Time - return PresetTimeCallback( - cond, as; initialize = initfn, finalize = finfun, - initializealg = reinitialization_alg(cb)) - else - # Periodic - return PeriodicCallback( - as, cond; initialize = initfn, finalize = finfun, - initializealg = reinitialization_alg(cb)) end end -function generate_discrete_callback(cb, sys, dvs, ps; postprocess_affect_expr! = nothing, - kwargs...) - if is_timed_condition(cb) - return generate_timed_callback(cb, sys, dvs, ps; postprocess_affect_expr!, - kwargs...) - else - c = compile_condition(cb, sys, dvs, ps; expression = Val{false}, kwargs...) - as = compile_affect(affects(cb), cb, sys, dvs, ps; expression = Val{false}, - postprocess_affect_expr!, kwargs...) - - user_initfun = _compile_optional_affect( - NULL_AFFECT, initialize_affects(cb), cb, sys, dvs, ps; kwargs...) - user_finfun = _compile_optional_affect( - NULL_AFFECT, finalize_affects(cb), cb, sys, dvs, ps; kwargs...) - if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing && - (save_idxs = get(ic.callback_to_clocks, cb, nothing)) !== nothing - initfn = let save_idxs = save_idxs, initfun = user_initfun - function (cb, u, t, integrator) - if !isnothing(initfun) - initfun(integrator) - end - for idx in save_idxs - SciMLBase.save_discretes!(integrator, idx) - end - end +""" +Initialize and finalize for VectorContinuousCallback. +""" +function wrap_vector_optional_affect(funs, default) + all(isnothing, funs) && return default + return let funs = funs + function (cb, u, t, integ) + for func in funs + isnothing(func) ? continue : func(integ) end - else - initfn = isnothing(user_initfun) ? SciMLBase.INITIALIZE_DEFAULT : - (_, _, _, i) -> user_initfun(i) end - finfun = isnothing(user_finfun) ? SciMLBase.FINALIZE_DEFAULT : - (_, _, _, i) -> user_finfun(i) - return DiscreteCallback( - c, as; initialize = initfn, finalize = finfun, - initializealg = reinitialization_alg(cb)) end end -function generate_discrete_callbacks(sys::AbstractSystem, dvs = unknowns(sys), - ps = parameters(sys; initial_parameters = true); kwargs...) - has_discrete_events(sys) || return nothing - symcbs = discrete_events(sys) - isempty(symcbs) && return nothing +function add_integrator_header( + sys::AbstractSystem, integrator = gensym(:MTKIntegrator), out = :u) + expr -> Func([DestructuredArgs(expr.args, integrator, inds = [:u, :p, :t])], [], + expr.body), + expr -> Func( + [DestructuredArgs(expr.args, integrator, inds = [out, :u, :p, :t])], [], + expr.body) +end + +""" +Compile an affect defined by a set of equations. Systems with algebraic equations will solve implicit discrete problems to obtain their next state. Systems without will generate functions that perform explicit updates. +""" +function compile_equational_affect( + aff::Union{AffectSystem, Vector{Equation}}, sys; reset_jumps = false, kwargs...) + if aff isa AbstractVector + aff = make_affect( + aff; iv = get_iv(sys), warn_no_algebraic = false) + end + affsys = system(aff) + ps_to_update = discretes(aff) + dvs_to_update = setdiff(unknowns(aff), getfield.(observed(sys), :lhs)) + aff_map = aff_to_sys(aff) + sys_map = Dict([v => k for (k, v) in aff_map]) + + if isempty(equations(affsys)) + update_eqs = Symbolics.fast_substitute( + observed(affsys), Dict([p => unPre(p) for p in parameters(affsys)])) + rhss = map(x -> x.rhs, update_eqs) + lhss = map(x -> aff_map[x.lhs], update_eqs) + is_p = [lhs ∈ Set(ps_to_update) for lhs in lhss] + is_u = [lhs ∈ Set(dvs_to_update) for lhs in lhss] + dvs = unknowns(sys) + ps = parameters(sys) + t = get_iv(sys) + + u_idxs = indexin((@view lhss[is_u]), dvs) + + wrap_mtkparameters = has_index_cache(sys) && (get_index_cache(sys) !== nothing) + p_idxs = if wrap_mtkparameters + [parameter_index(sys, p) for (i, p) in enumerate(lhss) + if is_p[i]] + else + indexin((@view lhss[is_p]), ps) + end + _ps = reorder_parameters(sys, ps) + integ = gensym(:MTKIntegrator) + + u_up, u_up! = build_function_wrapper(sys, (@view rhss[is_u]), dvs, _ps..., t; + wrap_code = add_integrator_header(sys, integ, :u), + expression = Val{false}, outputidxs = u_idxs, wrap_mtkparameters, cse = false) + p_up, p_up! = build_function_wrapper(sys, (@view rhss[is_p]), dvs, _ps..., t; + wrap_code = add_integrator_header(sys, integ, :p), + expression = Val{false}, outputidxs = p_idxs, wrap_mtkparameters, cse = false) + + return let dvs_to_update = dvs_to_update, ps_to_update = ps_to_update, + reset_jumps = reset_jumps, u_up! = u_up!, p_up! = p_up! - dbs = map(symcbs) do cb - generate_discrete_callback(cb, sys, dvs, ps; kwargs...) + function explicit_affect!(integ) + isempty(dvs_to_update) || u_up!(integ) + isempty(ps_to_update) || p_up!(integ) + reset_jumps && reset_aggregated_jumps!(integ) + end + end + else + return let dvs_to_update = dvs_to_update, aff_map = aff_map, sys_map = sys_map, + affsys = affsys, ps_to_update = ps_to_update, aff = aff, sys = sys + + dvs_to_access = [aff_map[u] for u in unknowns(affsys)] + ps_to_access = [unPre(p) for p in parameters(affsys)] + + affu_getter = getsym(sys, dvs_to_access) + affp_getter = getsym(sys, ps_to_access) + affu_setter! = setsym(affsys, unknowns(affsys)) + affp_setter! = setsym(affsys, parameters(affsys)) + u_setter! = setsym(sys, dvs_to_update) + p_setter! = setsym(sys, ps_to_update) + u_getter = getsym(affsys, [sys_map[u] for u in dvs_to_update]) + p_getter = getsym(affsys, [sys_map[p] for p in ps_to_update]) + + affprob = ImplicitDiscreteProblem(affsys, [dv => 0 for dv in unknowns(affsys)], + (0, 0), [p => 0. for p in parameters(affsys)]; + build_initializeprob = false, check_length = false) + + function implicit_affect!(integ) + new_u0 = affu_getter(integ) + affu_setter!(affprob, new_u0) + new_ps = affp_getter(integ) + affp_setter!(affprob, new_ps) + + affprob = remake( + affprob, tspan = (integ.t, integ.t)) + affsol = init(affprob, IDSolve()) + (check_error(affsol) === ReturnCode.InitialFailure) && + throw(UnsolvableCallbackError(all_equations(aff))) + + u_setter!(integ, u_getter(affsol)) + p_setter!(integ, p_getter(affsol)) + end + end end +end - dbs +struct UnsolvableCallbackError + eqs::Vector{Equation} +end + +function Base.showerror(io::IO, err::UnsolvableCallbackError) + println(io, + "The callback defined by the following equations:\n\n$(join(err.eqs, "\n"))\n\nis not solvable. Please check that the algebraic equations and affect equations are correct, and that all parameters intended to be changed are passed in as `discrete_parameters`.") end merge_cb(::Nothing, ::Nothing) = nothing @@ -1124,18 +968,109 @@ merge_cb(::Nothing, x) = merge_cb(x, nothing) merge_cb(x, ::Nothing) = x merge_cb(x, y) = CallbackSet(x, y) +""" +Generate the CallbackSet for a ODESystem or SDESystem. +""" function process_events(sys; callback = nothing, kwargs...) - if has_continuous_events(sys) && !isempty(continuous_events(sys)) - contin_cb = generate_rootfinding_callback(sys; kwargs...) - else - contin_cb = nothing + contin_cbs = generate_continuous_callbacks(sys; kwargs...) + discrete_cbs = generate_discrete_callbacks(sys; kwargs...) + cb = merge_cb(contin_cbs, callback) + (discrete_cbs === nothing) ? cb : CallbackSet(contin_cbs, discrete_cbs...) +end + +""" + discrete_events(sys::AbstractSystem) :: Vector{SymbolicDiscreteCallback} + +Returns a vector of all the `discrete_events` in an abstract system and its component subsystems. +The `SymbolicDiscreteCallback`s in the returned vector are structs with two fields: `condition` and +`affect` which correspond to the first and second elements of a `Pair` used to define an event, i.e. +`condition => affect`. + +See also `get_discrete_events`, which only returns the events of the top-level system. +""" +function discrete_events(sys::AbstractSystem) + obs = get_discrete_events(sys) + systems = get_systems(sys) + cbs = [obs; + reduce(vcat, + (map(cb -> namespace_callback(cb, s), discrete_events(s)) for s in systems), + init = SymbolicDiscreteCallback[])] + cbs +end + +has_discrete_events(sys::AbstractSystem) = isdefined(sys, :discrete_events) +function get_discrete_events(sys::AbstractSystem) + has_discrete_events(sys) || return SymbolicDiscreteCallback[] + getfield(sys, :discrete_events) +end + +""" + discrete_events_toplevel(sys::AbstractSystem) + +Replicates the behaviour of `discrete_events`, but ignores events of subsystems. + +Notes: +- Cannot be applied to non-complete systems. +""" +function discrete_events_toplevel(sys::AbstractSystem) + if has_parent(sys) && (parent = get_parent(sys)) !== nothing + return discrete_events_toplevel(parent) end - if has_discrete_events(sys) && !isempty(discrete_events(sys)) - discrete_cb = generate_discrete_callbacks(sys; kwargs...) - else - discrete_cb = nothing + return get_discrete_events(sys) +end + +""" + continuous_events(sys::AbstractSystem)::Vector{SymbolicContinuousCallback} + +Returns a vector of all the `continuous_events` in an abstract system and its component subsystems. +The `SymbolicContinuousCallback`s in the returned vector are structs with two fields: `eqs` and +`affect` which correspond to the first and second elements of a `Pair` used to define an event, i.e. +`eqs => affect`. + +See also `get_continuous_events`, which only returns the events of the top-level system. +""" +function continuous_events(sys::AbstractSystem) + obs = get_continuous_events(sys) + filter(!isempty, obs) + + systems = get_systems(sys) + cbs = [obs; + reduce(vcat, + (map(o -> namespace_callback(o, s), continuous_events(s)) for s in systems), + init = SymbolicContinuousCallback[])] + filter(!isempty, cbs) +end + +has_continuous_events(sys::AbstractSystem) = isdefined(sys, :continuous_events) +function get_continuous_events(sys::AbstractSystem) + has_continuous_events(sys) || return SymbolicContinuousCallback[] + getfield(sys, :continuous_events) +end + +""" + continuous_events_toplevel(sys::AbstractSystem) + +Replicates the behaviour of `continuous_events`, but ignores events of subsystems. + +Notes: +- Cannot be applied to non-complete systems. +""" +function continuous_events_toplevel(sys::AbstractSystem) + if has_parent(sys) && (parent = get_parent(sys)) !== nothing + return continuous_events_toplevel(parent) end + return get_continuous_events(sys) +end - cb = merge_cb(contin_cb, callback) - (discrete_cb === nothing) ? cb : CallbackSet(contin_cb, discrete_cb...) +""" +Process the symbolic events of a system. +""" +function create_symbolic_events(cont_events, disc_events, sys_eqs, iv) + alg_eqs = filter(eq -> eq.lhs isa Union{Symbolic, Number} && !is_diff_equation(eq), + sys_eqs) + cont_callbacks = to_cb_vector(cont_events; CB_TYPE = SymbolicContinuousCallback, + iv = iv, alg_eqs = alg_eqs, warn_no_algebraic = false) + disc_callbacks = to_cb_vector(disc_events; CB_TYPE = SymbolicDiscreteCallback, + iv = iv, alg_eqs = alg_eqs, warn_no_algebraic = false) + cont_callbacks, disc_callbacks end diff --git a/src/systems/clock_inference.jl b/src/systems/clock_inference.jl index b535773061..42fe28f7c7 100644 --- a/src/systems/clock_inference.jl +++ b/src/systems/clock_inference.jl @@ -1,7 +1,11 @@ struct ClockInference{S} + """Tearing state.""" ts::S + """The time domain (discrete clock, continuous) of each equation.""" eq_domain::Vector{TimeDomain} + """The output time domain (discrete clock, continuous) of each variable.""" var_domain::Vector{TimeDomain} + """The set of variables with concrete domains.""" inferred::BitSet end @@ -67,6 +71,9 @@ function substitute_sample_time(ex, dt) end end +""" +Update the equation-to-time domain mapping by inferring the time domain from the variables. +""" function infer_clocks!(ci::ClockInference) @unpack ts, eq_domain, var_domain, inferred = ci @unpack var_to_diff, graph = ts.structure @@ -132,6 +139,9 @@ function is_time_domain_conversion(v) input_timedomain(o) != output_timedomain(o) end +""" +For multi-clock systems, create a separate system for each clock in the system, along with associated equations. Return the updated tearing state, and the sets of clocked variables associated with each time domain. +""" function split_system(ci::ClockInference{S}) where {S} @unpack ts, eq_domain, var_domain, inferred = ci fullvars = get_fullvars(ts) @@ -143,11 +153,14 @@ function split_system(ci::ClockInference{S}) where {S} cid_to_eq = Vector{Int}[] var_to_cid = Vector{Int}(undef, ndsts(graph)) cid_to_var = Vector{Int}[] + # cid_counter = number of clocks cid_counter = Ref(0) for (i, d) in enumerate(eq_domain) cid = let cid_counter = cid_counter, id_to_clock = id_to_clock, continuous_id = continuous_id + # Fill the clock_to_id dict as you go, + # ContinuousClock() => 1, ... get!(clock_to_id, d) do cid = (cid_counter[] += 1) push!(id_to_clock, d) diff --git a/src/systems/codegen_utils.jl b/src/systems/codegen_utils.jl index a3fe53b95d..bedfdbcc37 100644 --- a/src/systems/codegen_utils.jl +++ b/src/systems/codegen_utils.jl @@ -179,17 +179,6 @@ function build_function_wrapper(sys::AbstractSystem, expr, args...; p_start = 2, args = ntuple(Val(length(args))) do i arg = args[i] - # for time-dependent systems, all arguments are passed through `time_varying_as_func` - # TODO: This is legacy behavior and a candidate for removal in v10 since we have callable - # parameters now. - if is_time_dependent(sys) - arg = if symbolic_type(arg) == NotSymbolic() - arg isa AbstractArray ? - map(x -> time_varying_as_func(unwrap(x), sys), arg) : arg - else - time_varying_as_func(unwrap(arg), sys) - end - end # Make sure to use the proper names for arguments if symbolic_type(arg) == NotSymbolic() && arg isa AbstractArray DestructuredArgs(arg, generated_argument_name(i); create_bindings) diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index 4a13d7ccf1..e061706416 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -335,9 +335,9 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; if length(unique(sysnames)) != length(sysnames) throw(ArgumentError("System names must be unique.")) end - cont_callbacks = SymbolicContinuousCallbacks(continuous_events) - disc_callbacks = SymbolicDiscreteCallbacks(discrete_events) + cont_callbacks, disc_callbacks = create_symbolic_events( + continuous_events, discrete_events, deqs, iv) if is_dde === nothing is_dde = _check_if_dde(deqs, iv′, systems) end @@ -470,7 +470,7 @@ Generates a function that computes the observed value(s) `ts` in the system `sys - `eval_expression = false`: If true and `expression = false`, evaluates the returned function in the module `eval_module` - `output_type = Array` the type of the array generated by a out-of-place vector-valued function - `param_only = false` if true, only allow the generated function to access system parameters -- `inputs = nothing` additinoal symbolic variables that should be provided to the generated function +- `inputs = Any[]` additional symbolic variables that should be provided to the generated function - `checkbounds = true` checks bounds if true when destructuring parameters - `op = Operator` sets the recursion terminator for the walk done by `vars` to identify the variables that appear in `ts`. See the documentation for `vars` for more detail. - `throw = true` if true, throw an error when generating a function for `ts` that reference variables that do not exist. @@ -500,8 +500,8 @@ For example, a function `g(op, unknowns, p..., inputs, t)` will be the in-place an array of inputs `inputs` is given, and `param_only` is false for a time-dependent system. """ function build_explicit_observed_function(sys, ts; - inputs = nothing, - disturbance_inputs = nothing, + inputs = Any[], + disturbance_inputs = Any[], disturbance_argument = false, expression = false, eval_expression = false, @@ -574,13 +574,13 @@ function build_explicit_observed_function(sys, ts; else (unknowns(sys),) end - if inputs === nothing + if isempty(inputs) inputs = () else ps = setdiff(ps, inputs) # Inputs have been converted to parameters by io_preprocessing, remove those from the parameter list inputs = (inputs,) end - if disturbance_inputs !== nothing + if !isempty(disturbance_inputs) # Disturbance inputs may or may not be included as inputs, depending on disturbance_argument ps = setdiff(ps, disturbance_inputs) end diff --git a/src/systems/diffeqs/sdesystem.jl b/src/systems/diffeqs/sdesystem.jl index c5299c28be..7f92853008 100644 --- a/src/systems/diffeqs/sdesystem.jl +++ b/src/systems/diffeqs/sdesystem.jl @@ -269,8 +269,10 @@ function SDESystem(deqs::AbstractVector{<:Equation}, neqs::AbstractArray, iv, dv ctrl_jac = RefValue{Any}(EMPTY_JAC) Wfact = RefValue(EMPTY_JAC) Wfact_t = RefValue(EMPTY_JAC) - cont_callbacks = SymbolicContinuousCallbacks(continuous_events) - disc_callbacks = SymbolicDiscreteCallbacks(discrete_events) + + cont_callbacks, disc_callbacks = create_symbolic_events( + continuous_events, discrete_events, deqs, iv) + if is_dde === nothing is_dde = _check_if_dde(deqs, iv′, systems) end diff --git a/src/systems/discrete_system/discrete_system.jl b/src/systems/discrete_system/discrete_system.jl index 075aa27e4d..4dfed954e5 100644 --- a/src/systems/discrete_system/discrete_system.jl +++ b/src/systems/discrete_system/discrete_system.jl @@ -428,4 +428,14 @@ function DiscreteFunctionExpr(sys::DiscreteSystem, args...; kwargs...) DiscreteFunctionExpr{true}(sys, args...; kwargs...) end +function Base.:(==)(sys1::DiscreteSystem, sys2::DiscreteSystem) + sys1 === sys2 && return true + isequal(nameof(sys1), nameof(sys2)) && + isequal(get_iv(sys1), get_iv(sys2)) && + _eq_unordered(get_eqs(sys1), get_eqs(sys2)) && + _eq_unordered(get_unknowns(sys1), get_unknowns(sys2)) && + _eq_unordered(get_ps(sys1), get_ps(sys2)) && + all(s1 == s2 for (s1, s2) in zip(get_systems(sys1), get_systems(sys2))) +end + supports_initialization(::DiscreteSystem) = false diff --git a/src/systems/discrete_system/implicit_discrete_system.jl b/src/systems/discrete_system/implicit_discrete_system.jl index 3956c089d4..b818ffcc8d 100644 --- a/src/systems/discrete_system/implicit_discrete_system.jl +++ b/src/systems/discrete_system/implicit_discrete_system.jl @@ -270,7 +270,8 @@ function flatten(sys::ImplicitDiscreteSystem, noeqs = false) end function generate_function( - sys::ImplicitDiscreteSystem, dvs = unknowns(sys), ps = parameters(sys); wrap_code = identity, kwargs...) + sys::ImplicitDiscreteSystem, dvs = unknowns(sys), ps = parameters(sys); + wrap_code = identity, cachesyms::Tuple = (), kwargs...) iv = get_iv(sys) # Algebraic equations get shifted forward 1, to match with differential equations exprs = map(equations(sys)) do eq @@ -286,8 +287,9 @@ function generate_function( u_next = map(Shift(iv, 1), dvs) u = dvs + p = (reorder_parameters(sys, unwrap.(ps))..., cachesyms...) build_function_wrapper( - sys, exprs, u_next, u, ps..., iv; p_start = 3, extra_assignments, kwargs...) + sys, exprs, u_next, u, p..., iv; p_start = 3, extra_assignments, kwargs...) end function shift_u0map_forward(sys::ImplicitDiscreteSystem, u0map, defs) @@ -296,12 +298,9 @@ function shift_u0map_forward(sys::ImplicitDiscreteSystem, u0map, defs) for k in collect(keys(u0map)) v = u0map[k] if !((op = operation(k)) isa Shift) - isnothing(getunshifted(k)) && - error("Initial conditions must be for the past state of the unknowns. Instead of providing the condition for $k, provide the condition for $(Shift(iv, -1)(k)).") - updated[k] = v elseif op.steps > 0 - error("Initial conditions must be for the past state of the unknowns. Instead of providing the condition for $k, provide the condition for $(Shift(iv, -1)(only(arguments(k)))).") + error("Initial conditions must be for the current or past state of the unknowns. Instead of providing the condition for $k, provide the condition for $(Shift(iv, -1)(only(arguments(k)))).") else updated[k] = v end @@ -380,6 +379,12 @@ function SciMLBase.ImplicitDiscreteFunction{iip, specialize}( f(u_next, u, p, t) = f_oop(u_next, u, p, t) f(resid, u_next, u, p, t) = f_iip(resid, u_next, u, p, t) + if length(dvs) == length(equations(sys)) + resid_prototype = nothing + else + resid_prototype = calculate_resid_prototype(length(equations(sys)), u0, p) + end + if specialize === SciMLBase.FunctionWrapperSpecialize && iip if u0 === nothing || p === nothing || t === nothing error("u0, p, and t must be specified for FunctionWrapperSpecialize on ImplicitDiscreteFunction.") @@ -394,6 +399,7 @@ function SciMLBase.ImplicitDiscreteFunction{iip, specialize}( sys = sys, observed = observedfun, analytic = analytic, + resid_prototype = resid_prototype, kwargs...) end @@ -441,3 +447,13 @@ end function ImplicitDiscreteFunctionExpr(sys::ImplicitDiscreteSystem, args...; kwargs...) ImplicitDiscreteFunctionExpr{true}(sys, args...; kwargs...) end + +function Base.:(==)(sys1::ImplicitDiscreteSystem, sys2::ImplicitDiscreteSystem) + sys1 === sys2 && return true + isequal(nameof(sys1), nameof(sys2)) && + isequal(get_iv(sys1), get_iv(sys2)) && + _eq_unordered(get_eqs(sys1), get_eqs(sys2)) && + _eq_unordered(get_unknowns(sys1), get_unknowns(sys2)) && + _eq_unordered(get_ps(sys1), get_ps(sys2)) && + all(s1 == s2 for (s1, s2) in zip(get_systems(sys1), get_systems(sys2))) +end diff --git a/src/systems/imperative_affect.jl b/src/systems/imperative_affect.jl index 4c9ff3d248..f01682deb1 100644 --- a/src/systems/imperative_affect.jl +++ b/src/systems/imperative_affect.jl @@ -99,8 +99,7 @@ function Base.hash(a::ImperativeAffect, s::UInt) hash(a.ctx, s) end -namespace_affects(af::ImperativeAffect, s) = namespace_affect(af, s) -function namespace_affect(affect::ImperativeAffect, s) +function namespace_affects(affect::ImperativeAffect, s) ImperativeAffect(func(affect), namespace_expr.(observed(affect), (s,)), observed_syms(affect), @@ -110,11 +109,49 @@ function namespace_affect(affect::ImperativeAffect, s) affect.skip_checks) end -function compile_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs...) - compile_user_affect(affect, cb, sys, dvs, ps; kwargs...) +function invalid_variables(sys, expr) + filter(x -> !any(isequal(x), all_symbols(sys)), reduce(vcat, vars(expr); init = [])) end -function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs...) +function unassignable_variables(sys, expr) + assignable_syms = reduce( + vcat, Symbolics.scalarize.(vcat( + unknowns(sys), parameters(sys; initial_parameters = true))); + init = []) + written = reduce(vcat, Symbolics.scalarize.(vars(expr)); init = []) + return filter( + x -> !any(isequal(x), assignable_syms), written) +end + +@generated function _generated_writeback(integ, setters::NamedTuple{NS1, <:Tuple}, + values::NamedTuple{NS2, <:Tuple}) where {NS1, NS2} + setter_exprs = [] + for name in NS2 + if !(name in NS1) + missing_name = "Tried to write back to $name from affect; only declared states ($NS1) may be written to." + error(missing_name) + end + push!(setter_exprs, :(setters.$name(integ, values.$name))) + end + return :(begin + $(setter_exprs...) + end) +end + +function check_assignable(sys, sym) + if symbolic_type(sym) == ScalarSymbolic() + is_variable(sys, sym) || is_parameter(sys, sym) + elseif symbolic_type(sym) == ArraySymbolic() + is_variable(sys, sym) || is_parameter(sys, sym) || + all(x -> check_assignable(sys, x), collect(sym)) + elseif sym isa Union{AbstractArray, Tuple} + all(x -> check_assignable(sys, x), sym) + else + false + end +end + +function compile_functional_affect(affect::ImperativeAffect, sys; kwargs...) #= Implementation sketch: generate observed function (oop), should save to a component array under obs_syms @@ -138,6 +175,9 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs. return (syms_dedup, exprs_dedup) end + dvs = unknowns(sys) + ps = parameters(sys) + obs_exprs = observed(affect) if !affect.skip_checks for oexpr in obs_exprs @@ -191,14 +231,8 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs. upd_funs = NamedTuple{mod_names}((setu.((sys,), first.(mod_pairs))...,)) - if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing - save_idxs = get(ic.callback_to_clocks, cb, Int[]) - else - save_idxs = Int[] - end - let user_affect = func(affect), ctx = context(affect) - function (integ) + @inline function (integ) # update the to-be-mutated values; this ensures that if you do a no-op then nothing happens modvals = mod_og_val_fun(integ.u, integ.p, integ.t) upd_component_array = NamedTuple{mod_names}(modvals) @@ -212,10 +246,6 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs. # write the new values back to the integrator _generated_writeback(integ, upd_funs, upd_vals) - - for idx in save_idxs - SciMLBase.save_discretes!(integ, idx) - end end end end diff --git a/src/systems/index_cache.jl b/src/systems/index_cache.jl index d0b687c212..948af4dfa5 100644 --- a/src/systems/index_cache.jl +++ b/src/systems/index_cache.jl @@ -117,10 +117,11 @@ function IndexCache(sys::AbstractSystem) affs = [affs] end for affect in affs - if affect isa Equation - is_parameter(sys, affect.lhs) && push!(discs, affect.lhs) - elseif affect isa FunctionalAffect || affect isa ImperativeAffect + if affect isa AffectSystem || affect isa FunctionalAffect || + affect isa ImperativeAffect union!(discs, unwrap.(discretes(affect))) + elseif isnothing(affect) + continue else error("Unhandled affect type $(typeof(affect))") end diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl index 06f5e1b623..5cfabefe1e 100644 --- a/src/systems/jumps/jumpsystem.jl +++ b/src/systems/jumps/jumpsystem.jl @@ -1,19 +1,5 @@ const JumpType = Union{VariableRateJump, ConstantRateJump, MassActionJump} -# modifies the expression representing an affect function to -# call reset_aggregated_jumps!(integrator). -# assumes iip -function _reset_aggregator!(expr, integrator) - @assert Meta.isexpr(expr, :function) - body = expr.args[end] - body = quote - $body - $reset_aggregated_jumps!($integrator) - end - expr.args[end] = body - return nothing -end - """ $(TYPEDEF) @@ -230,8 +216,10 @@ function JumpSystem(eqs, iv, unknowns, ps; end end - cont_callbacks = SymbolicContinuousCallbacks(continuous_events) - disc_callbacks = SymbolicDiscreteCallbacks(discrete_events) + cont_callbacks = to_cb_vector(continuous_events; CB_TYPE = SymbolicContinuousCallback, + iv = iv, warn_no_algebraic = false) + disc_callbacks = to_cb_vector(discrete_events; CB_TYPE = SymbolicDiscreteCallback, + iv = iv, warn_no_algebraic = false) JumpSystem{typeof(ap)}(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)), ap, iv′, us′, ps′, var_to_name, observed, name, description, systems, @@ -282,15 +270,13 @@ function generate_rate_function(js::JumpSystem, rate) expression = Val{true}) end -function generate_affect_function(js::JumpSystem, affect, outputidxs) +function generate_affect_function(js::JumpSystem, affect) consts = collect_constants(affect) if !isempty(consts) # The SymbolicUtils._build_function method of this case doesn't support postprocess_fbody csubs = Dict(c => getdefault(c) for c in consts) affect = substitute(affect, csubs) end - compile_affect( - affect, nothing, js, unknowns(js), parameters(js); outputidxs = outputidxs, - expression = Val{true}, checkvars = false) + compile_equational_affect(affect, js; expression = Val{true}, checkvars = false) end function assemble_vrj( @@ -299,8 +285,7 @@ function assemble_vrj( rate = GeneratedFunctionWrapper{(2, 3, is_split(js))}(rate, nothing) outputvars = (value(affect.lhs) for affect in vrj.affect!) outputidxs = [unknowntoid[var] for var in outputvars] - affect = eval_or_rgf(generate_affect_function(js, vrj.affect!, outputidxs); - eval_expression, eval_module) + affect = generate_affect_function(js, vrj.affect!) VariableRateJump(rate, affect; save_positions = vrj.save_positions) end @@ -308,10 +293,9 @@ function assemble_vrj_expr(js, vrj, unknowntoid) rate = generate_rate_function(js, vrj.rate) outputvars = (value(affect.lhs) for affect in vrj.affect!) outputidxs = ((unknowntoid[var] for var in outputvars)...,) - affect = generate_affect_function(js, vrj.affect!, outputidxs) + affect = generate_affect_function(js, vrj.affect!) quote rate = $rate - affect = $affect VariableRateJump(rate, affect) end @@ -323,8 +307,7 @@ function assemble_crj( rate = GeneratedFunctionWrapper{(2, 3, is_split(js))}(rate, nothing) outputvars = (value(affect.lhs) for affect in crj.affect!) outputidxs = [unknowntoid[var] for var in outputvars] - affect = eval_or_rgf(generate_affect_function(js, crj.affect!, outputidxs); - eval_expression, eval_module) + affect = generate_affect_function(js, crj.affect!) ConstantRateJump(rate, affect) end @@ -332,10 +315,9 @@ function assemble_crj_expr(js, crj, unknowntoid) rate = generate_rate_function(js, crj.rate) outputvars = (value(affect.lhs) for affect in crj.affect!) outputidxs = ((unknowntoid[var] for var in outputvars)...,) - affect = generate_affect_function(js, crj.affect!, outputidxs) + affect = generate_affect_function(js, crj.affect!) quote rate = $rate - affect = $affect ConstantRateJump(rate, affect) end @@ -574,8 +556,7 @@ function JumpProcesses.JumpProblem(js::JumpSystem, prob, end # handle events, making sure to reset aggregators in the generated affect functions - cbs = process_events(js; callback, eval_expression, eval_module, - postprocess_affect_expr! = _reset_aggregator!) + cbs = process_events(js; callback, eval_expression, eval_module, reset_jumps = true) JumpProblem(prob, aggregator, jset; dep_graph = jtoj, vartojumps_map = vtoj, jumptovars_map = jtov, scale_rates = false, nocopy = true, diff --git a/src/systems/model_parsing.jl b/src/systems/model_parsing.jl index 024c249363..2a852813a2 100644 --- a/src/systems/model_parsing.jl +++ b/src/systems/model_parsing.jl @@ -129,8 +129,9 @@ function _model_macro(mod, fullname::Union{Expr, Symbol}, expr, isconnector) Ref(dict), [:constants, :defaults, :kwargs, :structural_parameters]) sys = :($type($(flatten_equations)(equations), $iv, variables, parameters; - name, description = $description, systems, gui_metadata = $gui_metadata, defaults, - costs = [$(costs...)], constraints = [$(cons...)], consolidate = $consolidate)) + name, description = $description, systems, gui_metadata = $gui_metadata, + continuous_events = [$(c_evts...)], discrete_events = [$(d_evts...)], + defaults, costs = [$(costs...)], constraints = [$(cons...)], consolidate = $consolidate)) if length(ext) == 0 push!(exprs.args, :(var"#___sys___" = $sys)) @@ -141,16 +142,6 @@ function _model_macro(mod, fullname::Union{Expr, Symbol}, expr, isconnector) isconnector && push!(exprs.args, :($Setfield.@set!(var"#___sys___".connector_type=$connector_type(var"#___sys___")))) - !isempty(c_evts) && push!(exprs.args, - :($Setfield.@set!(var"#___sys___".continuous_events=$SymbolicContinuousCallback.([ - $(c_evts...) - ])))) - - !isempty(d_evts) && push!(exprs.args, - :($Setfield.@set!(var"#___sys___".discrete_events=$SymbolicDiscreteCallback.([ - $(d_evts...) - ])))) - f = if length(where_types) == 0 :($(Symbol(:__, name, :__))(; name, $(kwargs...)) = $exprs) else @@ -1144,8 +1135,16 @@ end function parse_continuous_events!(c_evts, dict, body) dict[:continuous_events] = [] Base.remove_linenums!(body) - for arg in body.args - push!(c_evts, arg) + for line in body.args + if length(line.args) == 3 && line.args[1] == :(=>) + push!(c_evts, :(($line, ()))) + elseif length(line.args) == 2 + event = line.args[1] + kwargs = parse_event_kwargs(line.args[2]) + push!(c_evts, :(($event, $kwargs))) + else + error("Malformed continuous event $line.") + end push!(dict[:continuous_events], readable_code.(c_evts)...) end end @@ -1153,12 +1152,30 @@ end function parse_discrete_events!(d_evts, dict, body) dict[:discrete_events] = [] Base.remove_linenums!(body) - for arg in body.args - push!(d_evts, arg) + for line in body.args + if length(line.args) == 3 && line.args[1] == :(=>) + push!(d_evts, :(($line, ()))) + elseif length(line.args) == 2 + event = line.args[1] + kwargs = parse_event_kwargs(line.args[2]) + push!(d_evts, :(($event, $kwargs))) + else + error("Malformed discrete event $line.") + end push!(dict[:discrete_events], readable_code.(d_evts)...) end end +function parse_event_kwargs(disc_expr) + kwargs = :([]) + for arg in disc_expr.args + (arg.head != :(=)) && error("Malformed event kwarg $arg.") + (arg.args[1] isa Symbol) || error("Invalid keyword argument name $(arg.args[1]).") + push!(kwargs.args, :($(QuoteNode(arg.args[1])) => $(arg.args[2]))) + end + kwargs +end + function parse_constraints!(cons, dict, body) dict[:constraints] = [] Base.remove_linenums!(body) diff --git a/src/systems/systems.jl b/src/systems/systems.jl index 52f93afb9b..4e4a7f7ba4 100644 --- a/src/systems/systems.jl +++ b/src/systems/systems.jl @@ -17,22 +17,23 @@ $(SIGNATURES) Structurally simplify algebraic equations in a system and compute the topological sort of the observed equations in `sys`. -### Optional Arguments: -+ optional argument `io` may take a tuple `(inputs, outputs)`. This will convert all `inputs` to parameters and allow them to be unconnected, i.e., simplification will allow models where `n_unknowns = n_equations - n_inputs`. - ### Optional Keyword Arguments: + When `simplify=true`, the `simplify` function will be applied during the tearing process. + `allow_symbolic=false`, `allow_parameter=true`, and `conservative=false` limit the coefficient types during tearing. In particular, `conservative=true` limits tearing to only solve for trivial linear systems where the coefficient has the absolute value of ``1``. + `fully_determined=true` controls whether or not an error will be thrown if the number of equations don't match the number of inputs, outputs, and equations. ++ `inputs`, `outputs` and `disturbance_inputs` are passed as keyword arguments.` All inputs` get converted to parameters and are allowed to be unconnected, allowing models where `n_unknowns = n_equations - n_inputs`. + `sort_eqs=true` controls whether equations are sorted lexicographically before simplification or not. """ function structural_simplify( - sys::AbstractSystem, io = nothing; additional_passes = [], simplify = false, split = true, + sys::AbstractSystem; additional_passes = [], simplify = false, split = true, allow_symbolic = false, allow_parameter = true, conservative = false, fully_determined = true, + inputs = Any[], outputs = Any[], + disturbance_inputs = Any[], kwargs...) isscheduled(sys) && throw(RepeatedStructuralSimplificationError()) - newsys′ = __structural_simplify(sys, io; simplify, + newsys′ = __structural_simplify(sys; simplify, allow_symbolic, allow_parameter, conservative, fully_determined, + inputs, outputs, disturbance_inputs, kwargs...) if newsys′ isa Tuple @assert length(newsys′) == 2 @@ -42,10 +43,6 @@ function structural_simplify( end if newsys isa DiscreteSystem && any(eq -> symbolic_type(eq.lhs) == NotSymbolic(), equations(newsys)) - error(""" - Encountered algebraic equations when simplifying discrete system. Please construct \ - an ImplicitDiscreteSystem instead. - """) end for pass in additional_passes newsys = pass(newsys) @@ -70,8 +67,10 @@ function __structural_simplify(sys::SDESystem, args...; kwargs...) return __structural_simplify(ODESystem(sys), args...; kwargs...) end -function __structural_simplify( - sys::AbstractSystem, io = nothing; simplify = false, sort_eqs = true, +function __structural_simplify(sys::AbstractSystem; simplify = false, + inputs = Any[], outputs = Any[], + disturbance_inputs = Any[], + sort_eqs = true, kwargs...) sys = expand_connections(sys) state = TearingState(sys; sort_eqs) @@ -90,7 +89,8 @@ function __structural_simplify( end end if isempty(brown_vars) - return structural_simplify!(state, io; simplify, kwargs...) + return structural_simplify!( + state; simplify, inputs, outputs, disturbance_inputs, kwargs...) else Is = Int[] Js = Int[] @@ -122,8 +122,8 @@ function __structural_simplify( for (i, v) in enumerate(fullvars) if !iszero(new_idxs[i]) && invview(var_to_diff)[i] === nothing] - # TODO: IO is not handled. - ode_sys = structural_simplify(sys, io; simplify, kwargs...) + ode_sys = structural_simplify( + sys; simplify, inputs, outputs, disturbance_inputs, kwargs...) eqs = equations(ode_sys) sorted_g_rows = zeros(Num, length(eqs), size(g, 2)) for (i, eq) in enumerate(eqs) @@ -157,9 +157,9 @@ function __structural_simplify( get_iv(ode_sys), unknowns(ode_sys), parameters(ode_sys); name = nameof(ode_sys), is_scalar_noise, observed = observed(ode_sys), defaults = defaults(sys), parameter_dependencies = parameter_dependencies(sys), assertions = assertions(sys), - guesses = guesses(sys), initialization_eqs = initialization_equations(sys)) - @set! ssys.tearing_state = get_tearing_state(ode_sys) - return ssys + guesses = guesses(sys), initialization_eqs = initialization_equations(sys), + continuous_events = continuous_events(sys), + discrete_events = discrete_events(sys)) end end diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index e0feb0d34d..1b561f97b4 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -657,29 +657,21 @@ function Base.show(io::IO, mime::MIME"text/plain", ms::MatchedSystemStructure) printstyled(io, " SelectedState") end -# TODO: clean up -function merge_io(io, inputs) - isempty(inputs) && return io - if io === nothing - io = (inputs, []) - else - io = ([inputs; io[1]], io[2]) - end - return io -end - -function structural_simplify!(state::TearingState, io = nothing; simplify = false, +function structural_simplify!(state::TearingState; simplify = false, check_consistency = true, fully_determined = true, warn_initialize_determined = true, + inputs = Any[], outputs = Any[], + disturbance_inputs = Any[], kwargs...) if state.sys isa ODESystem ci = ModelingToolkit.ClockInference(state) ci = ModelingToolkit.infer_clocks!(ci) time_domains = merge(Dict(state.fullvars .=> ci.var_domain), Dict(default_toterm.(state.fullvars) .=> ci.var_domain)) - tss, inputs, continuous_id, id_to_clock = ModelingToolkit.split_system(ci) - cont_io = merge_io(io, inputs[continuous_id]) - sys, input_idxs = _structural_simplify!(tss[continuous_id], cont_io; simplify, + tss, clocked_inputs, continuous_id, id_to_clock = ModelingToolkit.split_system(ci) + cont_inputs = [inputs; clocked_inputs[continuous_id]] + sys = _structural_simplify!(tss[continuous_id]; simplify, check_consistency, fully_determined, + inputs = cont_inputs, outputs, disturbance_inputs, kwargs...) if length(tss) > 1 if continuous_id > 0 @@ -695,8 +687,9 @@ function structural_simplify!(state::TearingState, io = nothing; simplify = fals discrete_subsystems[i] = sys continue end - dist_io = merge_io(io, inputs[i]) - ss, = _structural_simplify!(state, dist_io; simplify, check_consistency, + disc_inputs = [inputs; clocked_inputs[i]] + ss, = _structural_simplify!(state; simplify, check_consistency, + inputs = disc_inputs, outputs, disturbance_inputs, fully_determined, kwargs...) append!(appended_parameters, inputs[i], unknowns(ss)) discrete_subsystems[i] = ss @@ -713,31 +706,30 @@ function structural_simplify!(state::TearingState, io = nothing; simplify = fals for sym in get_ps(sys)] @set! sys.ps = ps else - sys, input_idxs = _structural_simplify!(state, io; simplify, check_consistency, + sys = _structural_simplify!(state; simplify, check_consistency, + inputs, outputs, disturbance_inputs, fully_determined, kwargs...) end - has_io = io !== nothing - return has_io ? (sys, input_idxs) : sys + return sys end -function _structural_simplify!(state::TearingState, io; simplify = false, +function _structural_simplify!(state::TearingState; simplify = false, check_consistency = true, fully_determined = true, warn_initialize_determined = false, dummy_derivative = true, + inputs = Any[], outputs = Any[], + disturbance_inputs = Any[], kwargs...) if fully_determined isa Bool check_consistency &= fully_determined else check_consistency = true end - has_io = io !== nothing + has_io = !isempty(inputs) || !isempty(outputs) !== nothing || + !isempty(disturbance_inputs) orig_inputs = Set() if has_io - ModelingToolkit.markio!(state, orig_inputs, io...) - end - if io !== nothing - state, input_idxs = ModelingToolkit.inputs_to_parameters!(state, io) - else - input_idxs = 0:-1 # Empty range + ModelingToolkit.markio!(state, orig_inputs, inputs, outputs, disturbance_inputs) + state = ModelingToolkit.inputs_to_parameters!(state, [inputs; disturbance_inputs]) end sys, mm = ModelingToolkit.alias_elimination!(state; kwargs...) if check_consistency @@ -761,5 +753,5 @@ function _structural_simplify!(state::TearingState, io; simplify = false, fullunknowns = [observables(sys); unknowns(sys)] @set! sys.observed = ModelingToolkit.topsort_equations(observed(sys), fullunknowns) - ModelingToolkit.invalidate_cache!(sys), input_idxs + ModelingToolkit.invalidate_cache!(sys) end diff --git a/src/utils.jl b/src/utils.jl index 1884a91c19..c0b9517419 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -513,18 +513,6 @@ function collect_applied_operators(x, op) end end -function find_derivatives!(vars, expr::Equation, f = identity) - (find_derivatives!(vars, expr.lhs, f); find_derivatives!(vars, expr.rhs, f); vars) -end -function find_derivatives!(vars, expr, f) - !iscall(O) && return vars - operation(O) isa Differential && push!(vars, f(O)) - for arg in arguments(O) - vars!(vars, arg) - end - return vars -end - """ $(TYPEDSIGNATURES) @@ -603,6 +591,7 @@ end function collect_var!(unknowns, parameters, var, iv; depth = 0) isequal(var, iv) && return nothing + var = unwrap(var) check_scope_depth(getmetadata(var, SymScope, LocalScope()), depth) || return nothing var = setmetadata(var, SymScope, LocalScope()) if iscalledparameter(var) @@ -632,8 +621,6 @@ function check_scope_depth(scope, depth) return depth == 0 elseif scope isa ParentScope return depth > 0 && check_scope_depth(scope.parent, depth - 1) - elseif scope isa DelayParentScope - return depth >= scope.N && check_scope_depth(scope.parent, depth - scope.N - 1) elseif scope isa GlobalScope return depth == -1 end diff --git a/src/variables.jl b/src/variables.jl index 510bd5c28d..eb2e54a268 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -354,6 +354,8 @@ function isdisturbance(x) Symbolics.getmetadata(x, VariableDisturbance, false) end +setdisturbance(x, v) = setmetadata(x, VariableDisturbance, v) + function disturbances(sys) [filter(isdisturbance, unknowns(sys)); filter(isdisturbance, parameters(sys))] end diff --git a/test/accessor_functions.jl b/test/accessor_functions.jl index 7ce477155b..4136736a8b 100644 --- a/test/accessor_functions.jl +++ b/test/accessor_functions.jl @@ -54,8 +54,8 @@ let D(Y) ~ -Y^3, O ~ (p_bot + d) * X_bot + Y ] - cevs = [[t ~ 1.0] => [Y ~ Y + 2.0]] - devs = [(t == 2.0) => [Y ~ Y + 2.0]] + cevs = [[t ~ 1.0] => [Y ~ Pre(Y) + 2.0]] + devs = [(t == 2.0) => [Y ~ Pre(Y) + 2.0]] @named sys_bot = ODESystem( eqs_bot, t; systems = [], continuous_events = cevs, discrete_events = devs) @named sys_mid2 = ODESystem( @@ -149,20 +149,19 @@ let for sys in [sys_bot, sys_mid2, sys_mid1, sys_top]) # Checks `continuous_events_toplevel` and `discrete_events_toplevel` (straightforward - # as I stored the same singe event in all systems). Don't check for non-toplevel cases as + # as I stored the same single event in all systems). Don't check for non-toplevel cases as # technically not needed for these tests and name spacing the events is a mess. - mtk_cev = ModelingToolkit.SymbolicContinuousCallback.(cevs)[1] - mtk_dev = ModelingToolkit.SymbolicDiscreteCallback.(devs)[1] + bot_cev = ModelingToolkit.SymbolicContinuousCallback( + cevs[1], alg_eqs = [O ~ (d + p_bot) * X_bot + Y]) + mid_dev = ModelingToolkit.SymbolicDiscreteCallback( + devs[1], alg_eqs = [O ~ (d + p_mid1) * X_mid1 + Y]) @test all_sets_equal( - continuous_events_toplevel.( - [sys_bot, sys_bot_comp, sys_bot_ss, sys_mid1, sys_mid1_comp, sys_mid1_ss, - sys_mid2, sys_mid2_comp, sys_mid2_ss, sys_top, sys_top_comp, sys_top_ss])..., - [mtk_cev]) + continuous_events_toplevel.([sys_bot, sys_bot_comp, sys_bot_ss])..., + [bot_cev]) @test all_sets_equal( discrete_events_toplevel.( - [sys_bot, sys_bot_comp, sys_bot_ss, sys_mid1, sys_mid1_comp, sys_mid1_ss, - sys_mid2, sys_mid2_comp, sys_mid2_ss, sys_top, sys_top_comp, sys_top_ss])..., - [mtk_dev]) + [sys_mid1, sys_mid1_comp, sys_mid1_ss])..., + [mid_dev]) @test all(sym_issubset( continuous_events_toplevel(sys), get_continuous_events(sys)) for sys in [sys_bot, sys_mid2, sys_mid1, sys_top]) diff --git a/test/clock.jl b/test/clock.jl index c6051a52a8..c4c64dbf90 100644 --- a/test/clock.jl +++ b/test/clock.jl @@ -65,10 +65,11 @@ By inference: ci, varmap = infer_clocks(sys) eqmap = ci.eq_domain tss, inputs, continuous_id = ModelingToolkit.split_system(deepcopy(ci)) -sss, = ModelingToolkit._structural_simplify!( - deepcopy(tss[continuous_id]), (inputs[continuous_id], ())) +sss = ModelingToolkit._structural_simplify!( + deepcopy(tss[continuous_id]), inputs = inputs[continuous_id], outputs = []) @test equations(sss) == [D(x) ~ u - x] -sss, = ModelingToolkit._structural_simplify!(deepcopy(tss[1]), (inputs[1], ())) +sss = ModelingToolkit._structural_simplify!( + deepcopy(tss[1]), inputs = inputs[1], outputs = []) @test isempty(equations(sss)) d = Clock(dt) k = ShiftIndex(d) diff --git a/test/code_generation.jl b/test/code_generation.jl index cf3d660b81..2fbf8f13d7 100644 --- a/test/code_generation.jl +++ b/test/code_generation.jl @@ -70,7 +70,7 @@ end @parameters p[1:2] (f::Function)(..) @named sys = ODESystem( [D(x[0]) ~ p[1] * x[0] + x[2], D(x[1]) ~ p[2] * f(x) + x[2]], t) - sys, = structural_simplify(sys, ([x[2]], [])) + sys = structural_simplify(sys, inputs = [x[2]], outputs = []) @test is_parameter(sys, x[2]) prob = ODEProblem(sys, [x[0] => 1.0, x[1] => 1.0], (0.0, 1.0), [p => ones(2), f => sum, x[2] => 2.0]) diff --git a/test/discrete_system.jl b/test/discrete_system.jl index b0e2481e56..2c515dda1a 100644 --- a/test/discrete_system.jl +++ b/test/discrete_system.jl @@ -254,7 +254,6 @@ end @variables x(t) y(t) k = ShiftIndex(t) @named sys = DiscreteSystem([x ~ x^2 + y^2, y ~ x(k - 1) + y(k - 1)], t) -@test_throws ["algebraic equations", "ImplicitDiscreteSystem"] structural_simplify(sys) @testset "Passing `nothing` to `u0`" begin @variables x(t) = 1 diff --git a/test/downstream/jump_control.jl b/test/downstream/jump_control.jl index a93146be10..dd0edb080d 100644 --- a/test/downstream/jump_control.jl +++ b/test/downstream/jump_control.jl @@ -98,7 +98,7 @@ end cost = [-x(1.0)] # Maximize the final distance. @named block = ODESystem( [D(x(t)) ~ v(t), D(v(t)) ~ u(t)], t; costs = cost, constraints = constr) - block, input_idxs = structural_simplify(block, ([u(t)], [])) + block = structural_simplify(block; inputs = [u(t)]) u0map = [x(t) => 0.0, v(t) => 0.0] tspan = (0.0, 1.0) @@ -137,7 +137,7 @@ end costs = [-q(tspan[2])] @named beesys = ODESystem(eqs, t; costs) - beesys, input_idxs = structural_simplify(beesys, ([α], [])) + beesys = structural_simplify(beesys; inputs = [α]) u0map = [w(t) => 40, q(t) => 2] pmap = [b => 1, c => 1, μ => 1, s => 1, ν => 1, α => 1] @@ -180,7 +180,7 @@ end costs = [-h(te)] cons = [T(te) ~ 0, m(te) ~ m_c] @named rocket = ODESystem(eqs, t; costs, constraints = cons) - rocket, input_idxs = structural_simplify(rocket, ([T(t)], [])) + rocket = structural_simplify(rocket; inputs = [T(t)]) u0map = [h(t) => h₀, m(t) => m₀, v(t) => 0] pmap = [ @@ -223,7 +223,7 @@ end costs = [-Symbolics.Integral(t in (0, tf))(x(t) - u(t)), -x(tf)] consolidate(u) = u[1] + u[2] @named rocket = ODESystem(eqs, t; costs, consolidate) - rocket, input_idxs = structural_simplify(rocket, ([u(t)], [])) + rocket = structural_simplify(rocket; inputs = [u(t)]) u0map = [x(t) => 17.5] pmap = [u(t) => 0.0, tf => 8] @@ -243,7 +243,7 @@ end @named block = ODESystem( [D(x(t)) ~ v(t), D(v(t)) ~ u(t)], t; costs = cost, constraints = constr) - block, input_idxs = structural_simplify(block, ([u(t)], [])) + block = structural_simplify(block, inputs = [u(t)]) u0map = [x(t) => 1.0, v(t) => 0.0] tspan = (0.0, tf) @@ -282,7 +282,7 @@ end tspan = (0, tf) @named cartpole = ODESystem(eqs, t; costs, constraints = cons) - cartpole, input_idxs = structural_simplify(cartpole, ([u], [])) + cartpole, input_idxs = structural_simplify(cartpole; inputs = [u]) u0map = [D(x(t)) => 0.0, D(θ(t)) => 0.0, θ(t) => 0.0, x(t) => 0.0] pmap = [mₖ => 1.0, mₚ => 0.2, l => 0.5, g => 9.81, u => 0] diff --git a/test/downstream/linearize.jl b/test/downstream/linearize.jl index d82bd696dd..20a9d317e8 100644 --- a/test/downstream/linearize.jl +++ b/test/downstream/linearize.jl @@ -87,10 +87,10 @@ connections = [f.y ~ c.r # filtered reference to controller reference @named cl = ODESystem(connections, t, systems = [f, c, p]) -lsys0, ssys = linearize(cl, [f.u], [p.x]) +lsys0, ssys = linearize(cl) desired_order = [f.x, p.x] lsys = ModelingToolkit.reorder_unknowns(lsys0, unknowns(ssys), desired_order) -lsys1, ssys = linearize(cl, [f.u], [p.x]; autodiff = AutoFiniteDiff()) +lsys1, ssys = linearize(cl; autodiff = AutoFiniteDiff()) lsys2 = ModelingToolkit.reorder_unknowns(lsys1, unknowns(ssys), desired_order) @test lsys.A == lsys2.A == [-2 0; 1 -2] @@ -266,7 +266,7 @@ closed_loop = ODESystem(connections, t, systems = [model, pid, filt, sensor, r, filt.xd => 0.0 ]) -@test_nowarn linearize(closed_loop, :r, :y; warn_empty_op = false) +@test_nowarn linearize(closed_loop; warn_empty_op = false) # https://discourse.julialang.org/t/mtk-change-in-linearize/115760/3 @mtkmodel Tank_noi begin diff --git a/test/extensions/ad.jl b/test/extensions/ad.jl index 14649b6bb6..73b1710812 100644 --- a/test/extensions/ad.jl +++ b/test/extensions/ad.jl @@ -60,7 +60,8 @@ end @parameters a b[1:3] c(t) d::Integer e[1:3] f[1:3, 1:3]::Int g::Vector{AbstractFloat} h::String @named sys = ODESystem( Equation[], t, [], [a, b, c, d, e, f, g, h], - continuous_events = [[a ~ 0] => [c ~ 0]]) + continuous_events = [ModelingToolkit.SymbolicContinuousCallback( + [a ~ 0] => [c ~ 0], discrete_parameters = c)]) sys = complete(sys) ivs = Dict(c => 3a, b => ones(3), a => 1.0, d => 4, e => [5.0, 6.0, 7.0], diff --git a/test/extensions/test_infiniteopt.jl b/test/extensions/test_infiniteopt.jl index 833b9f3275..a891311b46 100644 --- a/test/extensions/test_infiniteopt.jl +++ b/test/extensions/test_infiniteopt.jl @@ -24,8 +24,7 @@ end model = complete(model) inputs = [model.τ] outputs = [model.y] -model, _ = structural_simplify(model, (inputs, outputs)) - +model = structural_simplify(model; inputs, outputs) f, dvs, psym, io_sys = ModelingToolkit.generate_control_function( model, split = false) diff --git a/test/fmi/fmi.jl b/test/fmi/fmi.jl index 98c93398ff..edbbb312d6 100644 --- a/test/fmi/fmi.jl +++ b/test/fmi/fmi.jl @@ -158,7 +158,7 @@ end fmu = loadFMU(joinpath(FMU_DIR, "SimpleAdder.fmu"); type = :CS) @named adder = MTK.FMIComponent( Val(2); fmu, type = :CS, communication_step_size = 1e-6, - reinitializealg = BrownFullBasicInit()) + reinitializealg = BrownFullBasicInit(abstol = 1e-7)) @test MTK.isinput(adder.a) @test MTK.isinput(adder.b) @test MTK.isoutput(adder.out) @@ -211,7 +211,7 @@ end fmu = loadFMU(joinpath(FMU_DIR, "StateSpace.fmu"); type = :CS) @named sspace = MTK.FMIComponent( Val(3); fmu, communication_step_size = 1e-6, type = :CS, - reinitializealg = BrownFullBasicInit()) + reinitializealg = BrownFullBasicInit(abstol = 1e-7)) @test MTK.isinput(sspace.u) @test MTK.isoutput(sspace.y) @test !MTK.isinput(sspace.x) && !MTK.isoutput(sspace.x) diff --git a/test/funcaffect.jl b/test/funcaffect.jl index 05543c9161..1e7c66f39a 100644 --- a/test/funcaffect.jl +++ b/test/funcaffect.jl @@ -26,8 +26,7 @@ cb1 = ModelingToolkit.SymbolicDiscreteCallback(t == zr, (affect1!, [], [], [], [ @test cb == cb1 @test ModelingToolkit.SymbolicDiscreteCallback(cb) === cb # passthrough @test hash(cb) == hash(cb1) -ModelingToolkit.generate_discrete_callback(cb, sys, ModelingToolkit.get_variables(sys), - ModelingToolkit.get_ps(sys)); +ModelingToolkit.generate_callback(cb, sys); cb = ModelingToolkit.SymbolicContinuousCallback([t ~ zr], (f = affect1!, sts = [], pars = [], discretes = [], @@ -48,7 +47,7 @@ sys1 = ODESystem(eqs, t, [u], [], name = :sys, de = ModelingToolkit.get_discrete_events(sys1) @test length(de) == 1 de = de[1] -@test ModelingToolkit.condition(de) == [4.0] +@test ModelingToolkit.conditions(de) == [4.0] @test ModelingToolkit.has_functional_affect(de) sys2 = ODESystem(eqs, t, [u], [], name = :sys, diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl index 1ccc70607a..1a8952e80f 100644 --- a/test/initializationsystem.jl +++ b/test/initializationsystem.jl @@ -1290,9 +1290,9 @@ end @parameters β γ S0 @variables S(t)=S0 I(t) R(t) rate₁ = β * S * I - affect₁ = [S ~ S - 1, I ~ I + 1] + affect₁ = [S ~ Pre(S) - 1, I ~ Pre(I) + 1] rate₂ = γ * I - affect₂ = [I ~ I - 1, R ~ R + 1] + affect₂ = [I ~ Pre(I) - 1, R ~ Pre(R) + 1] j₁ = ConstantRateJump(rate₁, affect₁) j₂ = ConstantRateJump(rate₂, affect₂) j₃ = MassActionJump(2 * β + γ, [R => 1], [S => 1, R => -1]) diff --git a/test/input_output_handling.jl b/test/input_output_handling.jl index 2c16c89320..70ac772ca6 100644 --- a/test/input_output_handling.jl +++ b/test/input_output_handling.jl @@ -7,10 +7,10 @@ using ModelingToolkit: t_nounits as t, D_nounits as D @variables xx(t) some_input(t) [input = true] eqs = [D(xx) ~ some_input] @named model = ODESystem(eqs, t) -@test_throws ExtraVariablesSystemException structural_simplify(model, ((), ())) +@test_throws ExtraVariablesSystemException structural_simplify(model) if VERSION >= v"1.8" err = "In particular, the unset input(s) are:\n some_input(t)" - @test_throws err structural_simplify(model, ((), ())) + @test_throws err structural_simplify(model) end # Test input handling @@ -50,7 +50,7 @@ end @test !is_bound(sys31, sys1.v[2]) # simplification turns input variables into parameters -ssys, _ = structural_simplify(sys, ([u], [])) +ssys = structural_simplify(sys, inputs = [u], outputs = []) @test ModelingToolkit.isparameter(unbound_inputs(ssys)[]) @test !is_bound(ssys, u) @test u ∈ Set(unbound_inputs(ssys)) @@ -88,7 +88,7 @@ fsys4 = flatten(sys4) @variables x(t) y(t) [output = true] @test isoutput(y) @named sys = ODESystem([D(x) ~ -x, y ~ x], t) # both y and x are unbound -syss = structural_simplify(sys) # This makes y an observed variable +syss = structural_simplify(sys, outputs = [y]) # This makes y an observed variable @named sys2 = ODESystem([D(x) ~ -sys.x, y ~ sys.y], t, systems = [sys]) @@ -106,7 +106,7 @@ syss = structural_simplify(sys) # This makes y an observed variable @test isequal(unbound_outputs(sys2), [y]) @test isequal(bound_outputs(sys2), [sys.y]) -syss = structural_simplify(sys2) +syss = structural_simplify(sys2, outputs = [sys.y]) @test !is_bound(syss, y) @test !is_bound(syss, x) @@ -165,6 +165,7 @@ end ] @named sys = ODESystem(eqs, t) + sys = structural_simplify(sys, inputs = [u]) f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(sys; simplify, split) @test isequal(dvs[], x) @@ -182,8 +183,8 @@ end ] @named sys = ODESystem(eqs, t) - f, dvs, ps, io_sys = ModelingToolkit.generate_control_function( - sys, [u], [d]; simplify, split) + sys = structural_simplify(sys, inputs = [u], disturbance_inputs = [d]) + f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(sys; simplify, split) @test isequal(dvs[], x) @test isempty(ps) @@ -200,8 +201,9 @@ end ] @named sys = ODESystem(eqs, t) + sys = structural_simplify(sys, inputs = [u], disturbance_inputs = [d]) f, dvs, ps, io_sys = ModelingToolkit.generate_control_function( - sys, [u], [d]; simplify, split, disturbance_argument = true) + sys; simplify, split, disturbance_argument = true) @test isequal(dvs[], x) @test isempty(ps) @@ -265,9 +267,9 @@ eqs = [connect_sd(sd, mass1, mass2) @named _model = ODESystem(eqs, t) @named model = compose(_model, mass1, mass2, sd); +model = structural_simplify(model, inputs = [u]) f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(model, simplify = true) @test length(dvs) == 4 -@test length(ps) == length(parameters(model)) p = MTKParameters(io_sys, [io_sys.u => NaN]) x = ModelingToolkit.varmap_to_vars( merge(ModelingToolkit.defaults(model), @@ -281,7 +283,7 @@ i = findfirst(isequal(u[1]), out) @variables x(t) u(t) [input = true] eqs = [D(x) ~ u] @named sys = ODESystem(eqs, t) -@test_nowarn structural_simplify(sys, ([u], [])) +@test_nowarn structural_simplify(sys, inputs = [u], outputs = []) #= ## Disturbance input handling @@ -366,9 +368,9 @@ eqs = [D(y₁) ~ -k₁ * y₁ + k₃ * y₂ * y₃ + u1 @named sys = ODESystem(eqs, t) m_inputs = [u[1], u[2]] m_outputs = [y₂] -sys_simp, input_idxs = structural_simplify(sys, (; inputs = m_inputs, outputs = m_outputs)) +sys_simp = structural_simplify(sys, inputs = m_inputs, outputs = m_outputs) @test isequal(unknowns(sys_simp), collect(x[1:2])) -@test length(input_idxs) == 2 +@test length(inputs(sys_simp)) == 2 # https://github.com/SciML/ModelingToolkit.jl/issues/1577 @named c = Constant(; k = 2) @@ -389,7 +391,7 @@ sys = structural_simplify(model) ## Disturbance models when plant has multiple inputs using ModelingToolkit, LinearAlgebra -using ModelingToolkit: DisturbanceModel, io_preprocessing, get_iv, get_disturbance_system +using ModelingToolkit: DisturbanceModel, get_iv, get_disturbance_system using ModelingToolkitStandardLibrary.Blocks A, C = [randn(2, 2) for i in 1:2] B = [1.0 0; 0 1.0] @@ -433,6 +435,7 @@ matrices = ModelingToolkit.reorder_unknowns( ] @named sys = ODESystem(eqs, t) + sys = structural_simplify(sys, inputs = [u]) (; io_sys,) = ModelingToolkit.generate_control_function(sys, simplify = true) obsfn = ModelingToolkit.build_explicit_observed_function( io_sys, [x + u * t]; inputs = [u]) @@ -444,9 +447,9 @@ end @constants c = 2.0 @variables x(t) eqs = [D(x) ~ c * x] - @named sys = ODESystem(eqs, t, [x], []) + @mtkbuild sys = ODESystem(eqs, t, [x], []) - f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(sys, simplify = true) + f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(sys) @test f([0.5], nothing, MTKParameters(io_sys, []), 0.0) ≈ [1.0] end @@ -455,7 +458,8 @@ end @parameters p(::Real) = (x -> 2x) eqs = [D(x) ~ -x + p(u)] @named sys = ODESystem(eqs, t) - f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(sys, simplify = true) + sys = structural_simplify(sys, inputs = [u]) + f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(sys) p = MTKParameters(io_sys, []) u = [1.0] x = [1.0] diff --git a/test/jumpsystem.jl b/test/jumpsystem.jl index 1ee0408758..c5dbe1c56c 100644 --- a/test/jumpsystem.jl +++ b/test/jumpsystem.jl @@ -2,6 +2,7 @@ using ModelingToolkit, DiffEqBase, JumpProcesses, Test, LinearAlgebra using Random, StableRNGs, NonlinearSolve using OrdinaryDiffEq using ModelingToolkit: t_nounits as t, D_nounits as D +using BenchmarkTools MT = ModelingToolkit rng = StableRNG(12345) @@ -11,9 +12,9 @@ rng = StableRNG(12345) @constants h = 1 @variables S(t) I(t) R(t) rate₁ = β * S * I * h -affect₁ = [S ~ S - 1 * h, I ~ I + 1] +affect₁ = [S ~ Pre(S) - 1 * h, I ~ Pre(I) + 1] rate₂ = γ * I + t -affect₂ = [I ~ I - 1, R ~ R + 1] +affect₂ = [I ~ Pre(I) - 1, R ~ Pre(R) + 1] j₁ = ConstantRateJump(rate₁, affect₁) j₂ = VariableRateJump(rate₂, affect₂) @named js = JumpSystem([j₁, j₂], t, [S, I, R], [β, γ]) @@ -59,7 +60,7 @@ jump2.affect!(integrator) # test MT can make and solve a jump problem rate₃ = γ * I * h -affect₃ = [I ~ I * h - 1, R ~ R + 1] +affect₃ = [I ~ Pre(I) * h - 1, R ~ Pre(R) + 1] j₃ = ConstantRateJump(rate₃, affect₃) @named js2 = JumpSystem([j₁, j₃], t, [S, I, R], [β, γ]) js2 = complete(js2) @@ -247,7 +248,7 @@ end rate = k affect = [X ~ X - 1] -crj = ConstantRateJump(1.0, [X ~ X - 1]) +crj = ConstantRateJump(1.0, [X ~ Pre(X) - 1]) js1 = complete(JumpSystem([crj], t, [X], [k]; name = :js1)) js2 = complete(JumpSystem([crj], t, [X], []; name = :js2)) @@ -274,9 +275,9 @@ dp4 = DiscreteProblem(js4, u0, tspan) @parameters k @variables X(t) rate = k -affect = [X ~ X - 1] +affect = [X ~ Pre(X) - 1] -j1 = ConstantRateJump(k, [X ~ X - 1]) +j1 = ConstantRateJump(k, [X ~ Pre(X) - 1]) @test_nowarn @mtkbuild js1 = JumpSystem([j1], t, [X], [k]) # test correct autosolver is selected, which implies appropriate dep graphs are available @@ -284,8 +285,8 @@ let @parameters k @variables X(t) rate = k - affect = [X ~ X - 1] - j1 = ConstantRateJump(k, [X ~ X - 1]) + affect = [X ~ Pre(X) - 1] + j1 = ConstantRateJump(k, [X ~ Pre(X) - 1]) Nv = [1, JumpProcesses.USE_DIRECT_THRESHOLD + 1, JumpProcesses.USE_RSSA_THRESHOLD + 1] algtypes = [Direct, RSSA, RSSACR] @@ -304,7 +305,7 @@ let Random.seed!(rng, 1111) @variables A(t) B(t) C(t) @parameters k - vrj = VariableRateJump(k * (sin(t) + 1), [A ~ A + 1, C ~ C + 2]) + vrj = VariableRateJump(k * (sin(t) + 1), [A ~ Pre(A) + 1, C ~ Pre(C) + 2]) js = complete(JumpSystem([vrj], t, [A, C], [k]; name = :js, observed = [B ~ C * A])) oprob = ODEProblem(js, [A => 0, C => 0], (0.0, 10.0), [k => 1.0]) jprob = JumpProblem(js, oprob, Direct(); rng) @@ -345,9 +346,9 @@ end let @variables x1(t) x2(t) x3(t) x4(t) x5(t) @parameters p1 p2 p3 p4 p5 - j1 = ConstantRateJump(p1, [x1 ~ x1 + 1]) + j1 = ConstantRateJump(p1, [x1 ~ Pre(x1) + 1]) j2 = MassActionJump(p2, [x2 => 1], [x3 => -1]) - j3 = VariableRateJump(p3, [x3 ~ x3 + 1, x4 ~ x4 + 1]) + j3 = VariableRateJump(p3, [x3 ~ Pre(x3) + 1, x4 ~ Pre(x4) + 1]) j4 = MassActionJump(p4 * p5, [x1 => 1, x5 => 1], [x1 => -1, x5 => -1, x2 => 1]) us = Set() ps = Set() @@ -378,22 +379,20 @@ end # scoping tests let - @variables x1(t) x2(t) x3(t) x4(t) x5(t) + @variables x1(t) x2(t) x3(t) x4(t) x2 = ParentScope(x2) x3 = ParentScope(ParentScope(x3)) - x4 = DelayParentScope(x4) - x5 = GlobalScope(x5) - @parameters p1 p2 p3 p4 p5 + x4 = GlobalScope(x4) + @parameters p1 p2 p3 p4 p2 = ParentScope(p2) p3 = ParentScope(ParentScope(p3)) - p4 = DelayParentScope(p4) - p5 = GlobalScope(p5) + p4 = GlobalScope(p4) - j1 = ConstantRateJump(p1, [x1 ~ x1 + 1]) + j1 = ConstantRateJump(p1, [x1 ~ Pre(x1) + 1]) j2 = MassActionJump(p2, [x2 => 1], [x3 => -1]) - j3 = VariableRateJump(p3, [x3 ~ x3 + 1, x4 ~ x4 + 1]) - j4 = MassActionJump(p4 * p5, [x1 => 1, x5 => 1], [x1 => -1, x5 => -1, x2 => 1]) - @named js = JumpSystem([j1, j2, j3, j4], t, [x1, x2, x3, x4, x5], [p1, p2, p3, p4, p5]) + j3 = VariableRateJump(p3, [x3 ~ Pre(x3) + 1, x4 ~ Pre(x4) + 1]) + j4 = MassActionJump(p4 * p4, [x1 => 1, x4 => 1], [x1 => -1, x4 => -1, x2 => 1]) + @named js = JumpSystem([j1, j2, j3, j4], t, [x1, x2, x3, x4], [p1, p2, p3, p4]) us = Set() ps = Set() @@ -414,13 +413,13 @@ let empty!.((us, ps)) MT.collect_scoped_vars!(us, ps, js, iv; depth = 2) - @test issetequal(us, [x3, x4]) - @test issetequal(ps, [p3, p4]) + @test issetequal(us, [x3]) + @test issetequal(ps, [p3]) empty!.((us, ps)) MT.collect_scoped_vars!(us, ps, js, iv; depth = -1) - @test issetequal(us, [x5]) - @test issetequal(ps, [p5]) + @test issetequal(us, [x4]) + @test issetequal(ps, [p4]) end # PDMP test @@ -429,8 +428,8 @@ let Random.seed!(rng, seed) @variables X(t) Y(t) @parameters k1 k2 - vrj1 = VariableRateJump(k1 * X, [X ~ X - 1]; save_positions = (false, false)) - vrj2 = VariableRateJump(k1, [Y ~ Y + 1]; save_positions = (false, false)) + vrj1 = VariableRateJump(k1 * X, [X ~ Pre(X) - 1]; save_positions = (false, false)) + vrj2 = VariableRateJump(k1, [Y ~ Pre(Y) + 1]; save_positions = (false, false)) eqs = [D(X) ~ k2, D(Y) ~ -k2 / 10 * Y] @named jsys = JumpSystem([vrj1, vrj2, eqs[1], eqs[2]], t, [X, Y], [k1, k2]) jsys = complete(jsys) @@ -471,8 +470,8 @@ let Random.seed!(rng, seed) @variables X(t) Y(t) @parameters α β - vrj = VariableRateJump(β * X, [X ~ X - 1]; save_positions = (false, false)) - crj = ConstantRateJump(β * Y, [Y ~ Y - 1]) + vrj = VariableRateJump(β * X, [X ~ Pre(X) - 1]; save_positions = (false, false)) + crj = ConstantRateJump(β * Y, [Y ~ Pre(Y) - 1]) maj = MassActionJump(α, [0 => 1], [Y => 1]) eqs = [D(X) ~ α * (1 + Y)] @named jsys = JumpSystem([maj, crj, vrj, eqs[1]], t, [X, Y], [α, β]) @@ -539,8 +538,8 @@ end @variables X(t) rate1 = p rate2 = X * d - affect1 = [X ~ X + 1] - affect2 = [X ~ X - 1] + affect1 = [X ~ Pre(X) + 1] + affect2 = [X ~ Pre(X) - 1] j1 = ConstantRateJump(rate1, affect1) j2 = ConstantRateJump(rate2, affect2) @@ -557,7 +556,7 @@ end @parameters a b eq = D(X) ~ a rate = b * X - affect = [X ~ X - 1] + affect = [X ~ Pre(X) - 1] crj = ConstantRateJump(rate, affect) @named jsys = JumpSystem([crj, eq], t, [X], [a, b]) jsys = complete(jsys) diff --git a/test/model_parsing.jl b/test/model_parsing.jl index fe2bcbfca6..ad457a0ba3 100644 --- a/test/model_parsing.jl +++ b/test/model_parsing.jl @@ -484,7 +484,7 @@ using ModelingToolkit: D_nounits [x ~ 1.5] => [x ~ 5, y ~ 1] end @discrete_events begin - (t == 1.5) => [x ~ x + 5, z ~ 2] + (t == 1.5) => [x ~ Pre(x) + 5, z ~ 2] end end diff --git a/test/mtkparameters.jl b/test/mtkparameters.jl index 55de0768e0..5cc26ba79a 100644 --- a/test/mtkparameters.jl +++ b/test/mtkparameters.jl @@ -10,7 +10,8 @@ using JET @parameters a b c(t) d::Integer e[1:3] f[1:3, 1:3]::Int g::Vector{AbstractFloat} h::String @named sys = ODESystem( Equation[], t, [], [a, c, d, e, f, g, h], parameter_dependencies = [b ~ 2a], - continuous_events = [[a ~ 0] => [c ~ 0]], defaults = Dict(a => 0.0)) + continuous_events = [ModelingToolkit.SymbolicContinuousCallback( + [a ~ 0] => [c ~ 0], discrete_parameters = c)], defaults = Dict(a => 0.0)) sys = complete(sys) ivs = Dict(c => 3a, d => 4, e => [5.0, 6.0, 7.0], diff --git a/test/odesystem.jl b/test/odesystem.jl index 9219fca5fb..68a62500f8 100644 --- a/test/odesystem.jl +++ b/test/odesystem.jl @@ -1,5 +1,6 @@ using ModelingToolkit, StaticArrays, LinearAlgebra -using ModelingToolkit: get_metadata, MTKParameters +using ModelingToolkit: get_metadata, MTKParameters, SymbolicDiscreteCallback, + SymbolicContinuousCallback using SymbolicIndexingInterface using OrdinaryDiffEq, Sundials using DiffEqBase, SparseArrays @@ -9,6 +10,7 @@ using SymbolicUtils: issym using ForwardDiff using ModelingToolkit: value using ModelingToolkit: t_nounits as t, D_nounits as D +using Symbolics: unwrap # Define some variables @parameters σ ρ β @@ -135,19 +137,7 @@ du = zeros(3) tgrad_iip(du, u, p, t) @test du == [0.0, -u[2], 0.0] -@parameters σ′(t - 1) -eqs = [D(x) ~ σ′ * (y - x), - D(y) ~ x * (ρ - z) - y, - D(z) ~ x * y - β * z * κ] -@named de = ODESystem(eqs, t) -test_diffeq_inference("global iv-varying", de, t, (x, y, z), (σ′, ρ, β)) - -f = generate_function(de, [x, y, z], [σ′, ρ, β], expression = Val{false})[2] -du = [0.0, 0.0, 0.0] -f(du, [1.0, 2.0, 3.0], [x -> x + 7, 2, 3], 5.0) -@test du ≈ [11, -3, -7] - -@parameters σ(..) +@parameters (σ::Function)(..) eqs = [D(x) ~ σ(t - 1) * (y - x), D(y) ~ x * (ρ - z) - y, D(z) ~ x * y - β * z * κ] @@ -477,7 +467,7 @@ let eqs = [D(x) ~ ẋ, D(ẋ) ~ f - k * x - d * ẋ] @named sys = ODESystem(eqs, t, [x, ẋ], [d, k]) - sys, _ = structural_simplify(sys, ([f], [])) + sys = structural_simplify(sys; inputs = [f]) @test isequal(calculate_control_jacobian(sys), reshape(Num[0, 1], 2, 1)) @@ -1026,24 +1016,27 @@ prob = ODEProblem(sys, [x => 1.0], (0.0, 10.0)) @test_nowarn solve(prob, Tsit5()) # Issue#2383 -@variables x(t)[1:3] -@parameters p[1:3, 1:3] -eqs = [ - D(x) ~ p * x -] -@mtkbuild sys = ODESystem(eqs, t; continuous_events = [[norm(x) ~ 3.0] => [x ~ ones(3)]]) -# array affect equations used to not work -prob1 = @test_nowarn ODEProblem(sys, [x => ones(3)], (0.0, 10.0), [p => ones(3, 3)]) -sol1 = @test_nowarn solve(prob1, Tsit5()) +@testset "Arrays in affect/condition equations" begin + @variables x(t)[1:3] + @parameters p[1:3, 1:3] + eqs = [ + D(x) ~ p * x + ] + @mtkbuild sys = ODESystem( + eqs, t; continuous_events = [[norm(x) ~ 3.0] => [x ~ ones(3)]]) + # array affect equations used to not work + prob1 = @test_nowarn ODEProblem(sys, [x => ones(3)], (0.0, 10.0), [p => ones(3, 3)]) + sol1 = @test_nowarn solve(prob1, Tsit5()) -# array condition equations also used to not work -@mtkbuild sys = ODESystem( - eqs, t; continuous_events = [[x ~ sqrt(3) * ones(3)] => [x ~ ones(3)]]) -# array affect equations used to not work -prob2 = @test_nowarn ODEProblem(sys, [x => ones(3)], (0.0, 10.0), [p => ones(3, 3)]) -sol2 = @test_nowarn solve(prob2, Tsit5()) + # array condition equations also used to not work + @mtkbuild sys = ODESystem( + eqs, t; continuous_events = [[x ~ sqrt(3) * ones(3)] => [x ~ ones(3)]]) + # array affect equations used to not work + prob2 = @test_nowarn ODEProblem(sys, [x => ones(3)], (0.0, 10.0), [p => ones(3, 3)]) + sol2 = @test_nowarn solve(prob2, Tsit5()) -@test sol1 ≈ sol2 + @test sol1.u ≈ sol2.u[2:end] +end # Requires fix in symbolics for `linear_expansion(p * x, D(y))` @test_skip begin @@ -1190,10 +1183,12 @@ end end # Namespacing of array variables -@variables x(t)[1:2] -@named sys = ODESystem(Equation[], t) -@test getname(unknowns(sys, x)) == :sys₊x -@test size(unknowns(sys, x)) == size(x) +@testset "Namespacing of array variables" begin + @variables x(t)[1:2] + @named sys = ODESystem(Equation[], t) + @test getname(unknowns(sys, x)) == :sys₊x + @test size(unknowns(sys, x)) == size(x) +end # Issue#2667 and Issue#2953 @testset "ForwardDiff through ODEProblem constructor" begin @@ -1294,11 +1289,11 @@ end @named sys = ODESystem( [D(u) ~ (sum(u) + sum(x) + sum(p) + sum(o)) * x, o ~ prod(u) * x], t, [u..., x..., o...], [p...]) - sys1, = structural_simplify(sys, ([x...], [])) + sys1 = structural_simplify(sys, inputs = [x...], outputs = []) fn1, = ModelingToolkit.generate_function(sys1; expression = Val{false}) ps = MTKParameters(sys1, [x => 2ones(2), p => 3ones(2, 2)]) @test_nowarn fn1(ones(4), ps, 4.0) - sys2, = structural_simplify(sys, ([x...], []); split = false) + sys2 = structural_simplify(sys, inputs = [x...], outputs = [], split = false) fn2, = ModelingToolkit.generate_function(sys2; expression = Val{false}) ps = zeros(8) setp(sys2, x)(ps, 2ones(2)) @@ -1407,7 +1402,7 @@ end o[2] ~ sum(p) * sum(x)] @named sys = ODESystem(eqs, t, [u..., x..., o], [p...]) - sys1, = structural_simplify(sys, ([x...], [o...]), split = false) + sys1 = structural_simplify(sys, inputs = [x...], outputs = [o...], split = false) @test_nowarn ModelingToolkit.build_explicit_observed_function(sys1, u; inputs = [x...]) @@ -1531,8 +1526,12 @@ end @testset "Observed variables dependent on discrete parameters" begin @variables x(t) obs(t) @parameters c(t) - @mtkbuild sys = ODESystem( - [D(x) ~ c * cos(x), obs ~ c], t, [x], [c]; discrete_events = [1.0 => [c ~ c + 1]]) + @mtkbuild sys = ODESystem([D(x) ~ c * cos(x), obs ~ c], + t, + [x], + [c]; + discrete_events = [SymbolicDiscreteCallback( + 1.0 => [c ~ Pre(c) + 1], discrete_parameters = [c])]) prob = ODEProblem(sys, [x => 0.0], (0.0, 2pi), [c => 1.0]) sol = solve(prob, Tsit5()) @test sol[obs] ≈ 1:7 @@ -1592,15 +1591,16 @@ end # Test `isequal` @testset "`isequal`" begin @variables X(t) - @parameters p d + @parameters p d(t) eq = D(X) ~ p - d * X osys1 = complete(ODESystem([eq], t; name = :osys)) osys2 = complete(ODESystem([eq], t; name = :osys)) @test osys1 == osys2 # true - continuous_events = [[X ~ 1.0] => [X ~ X + 5.0]] - discrete_events = [5.0 => [d ~ d / 2.0]] + continuous_events = [[X ~ 1.0] => [X ~ Pre(X) + 5.0]] + discrete_events = [SymbolicDiscreteCallback( + 5.0 => [d ~ d / 2.0], discrete_parameters = [d])] osys1 = complete(ODESystem([eq], t; name = :osys, continuous_events)) osys2 = complete(ODESystem([eq], t; name = :osys)) @@ -1730,3 +1730,26 @@ end @test obsfn_expr_oop isa Expr @test obsfn_expr_iip isa Expr end + +@testset "`@named` always wraps in `ParentScope`" begin + function SysA(; name, var1) + @variables x(t) + scope = ModelingToolkit.getmetadata(unwrap(var1), ModelingToolkit.SymScope, nothing) + @test scope isa ParentScope + @test scope.parent isa ParentScope + @test scope.parent.parent isa LocalScope + return ODESystem(D(x) ~ var1, t; name) + end + function SysB(; name, var1) + @variables x(t) + @named subsys = SysA(; var1) + return ODESystem(D(x) ~ x, t; systems = [subsys], name) + end + function SysC(; name) + @variables x(t) + @named subsys = SysB(; var1 = x) + return ODESystem(D(x) ~ x, t; systems = [subsys], name) + end + @mtkbuild sys = SysC() + @test length(unknowns(sys)) == 3 +end diff --git a/test/parameter_dependencies.jl b/test/parameter_dependencies.jl index 31881e1ca8..bdaa2a391b 100644 --- a/test/parameter_dependencies.jl +++ b/test/parameter_dependencies.jl @@ -1,6 +1,7 @@ using ModelingToolkit using Test -using ModelingToolkit: t_nounits as t, D_nounits as D +using ModelingToolkit: t_nounits as t, D_nounits as D, SymbolicDiscreteCallback, + SymbolicContinuousCallback using OrdinaryDiffEq using StochasticDiffEq using JumpProcesses @@ -10,14 +11,14 @@ using SymbolicIndexingInterface using NonlinearSolve @testset "ODESystem with callbacks" begin - @parameters p1=1.0 p2 + @parameters p1(t)=1.0 p2 @variables x(t) - cb1 = [x ~ 2.0] => [p1 ~ 2.0] # triggers at t=-2+√6 + cb1 = SymbolicContinuousCallback([x ~ 2.0] => [p1 ~ 2.0], discrete_parameters = [p1]) # triggers at t=-2+√6 function affect1!(integ, u, p, ctx) integ.ps[p[1]] = integ.ps[p[2]] end cb2 = [x ~ 4.0] => (affect1!, [], [p1, p2], [p1]) # triggers at t=-2+√7 - cb3 = [1.0] => [p1 ~ 5.0] + cb3 = SymbolicDiscreteCallback([1.0] => [p1 ~ 5.0], discrete_parameters = [p1]) @mtkbuild sys = ODESystem( [D(x) ~ p1 * t + p2], @@ -203,7 +204,7 @@ end @testset "Clock system" begin dt = 0.1 @variables x(t) y(t) u(t) yd(t) ud(t) r(t) z(t) - @parameters kp kq + @parameters kp(t) kq d = Clock(dt) k = ShiftIndex(d) @@ -225,7 +226,8 @@ end @test_nowarn solve(prob, Tsit5()) @mtkbuild sys = ODESystem(eqs, t; parameter_dependencies = [kq => 2kp], - discrete_events = [[0.5] => [kp ~ 2.0]]) + discrete_events = [SymbolicDiscreteCallback( + [0.5] => [kp ~ 2.0], discrete_parameters = [kp])]) prob = ODEProblem(sys, [x => 0.0, y => 0.0], (0.0, Tf), [kp => 1.0; z(k - 1) => 3.0; yd(k - 1) => 0.0; z(k - 2) => 4.0; yd(k - 2) => 2.0]) @@ -245,7 +247,7 @@ end end @testset "SDESystem" begin - @parameters σ ρ β + @parameters σ(t) ρ β @variables x(t) y(t) z(t) eqs = [D(x) ~ σ * (y - x), @@ -269,7 +271,8 @@ end @named sys = ODESystem(eqs, t) @named sdesys = SDESystem(sys, noiseeqs; parameter_dependencies = [ρ => 2σ], - discrete_events = [[10.0] => [σ ~ 15.0]]) + discrete_events = [SymbolicDiscreteCallback( + [10.0] => [σ ~ 15.0], discrete_parameters = [σ])]) sdesys = complete(sdesys) prob = SDEProblem( sdesys, [x => 1.0, y => 0.0, z => 0.0], (0.0, 100.0), [σ => 10.0, β => 2.33]) @@ -283,17 +286,17 @@ end @testset "JumpSystem" begin rng = StableRNG(12345) - @parameters β γ + @parameters β γ(t) @constants h = 1 @variables S(t) I(t) R(t) rate₁ = β * S * I * h - affect₁ = [S ~ S - 1 * h, I ~ I + 1] + affect₁ = [S ~ Pre(S) - 1 * h, I ~ Pre(I) + 1] rate₃ = γ * I * h - affect₃ = [I ~ I * h - 1, R ~ R + 1] + affect₃ = [I ~ Pre(I) * h - 1, R ~ Pre(R) + 1] j₁ = ConstantRateJump(rate₁, affect₁) j₃ = ConstantRateJump(rate₃, affect₃) @named js2 = JumpSystem( - [j₁, j₃], t, [S, I, R], [γ]; parameter_dependencies = [β => 0.01γ]) + [j₃], t, [S, I, R], [γ]; parameter_dependencies = [β => 0.01γ]) @test isequal(only(parameters(js2)), γ) @test Set(full_parameters(js2)) == Set([γ, β]) js2 = complete(js2) @@ -308,7 +311,8 @@ end @named js2 = JumpSystem( [j₁, j₃], t, [S, I, R], [γ]; parameter_dependencies = [β => 0.01γ], - discrete_events = [[10.0] => [γ ~ 0.02]]) + discrete_events = [SymbolicDiscreteCallback( + [10.0] => [γ ~ 0.02], discrete_parameters = [γ])]) js2 = complete(js2) dprob = DiscreteProblem(js2, u₀map, tspan, parammap) jprob = JumpProblem(js2, dprob, Direct(), save_positions = (false, false), rng = rng) diff --git a/test/reduction.jl b/test/reduction.jl index fa9029a652..9a515614d5 100644 --- a/test/reduction.jl +++ b/test/reduction.jl @@ -233,7 +233,7 @@ eqs = [D(x) ~ σ * (y - x) u ~ z + a] lorenz1 = ODESystem(eqs, t, name = :lorenz1) -lorenz1_reduced, _ = structural_simplify(lorenz1, ([z], [])) +lorenz1_reduced = structural_simplify(lorenz1, inputs = [z], outputs = []) @test z in Set(parameters(lorenz1_reduced)) # #2064 diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index 9099d32d14..e208a74dd5 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -1,10 +1,11 @@ using ModelingToolkit, OrdinaryDiffEq, StochasticDiffEq, JumpProcesses, Test using SciMLStructures: canonicalize, Discrete using ModelingToolkit: SymbolicContinuousCallback, - SymbolicContinuousCallbacks, NULL_AFFECT, + SymbolicDiscreteCallback, get_callback, t_nounits as t, - D_nounits as D + D_nounits as D, + affects, affect_negs, system, observed, AffectSystem using StableRNGs import SciMLBase using SymbolicIndexingInterface @@ -17,215 +18,74 @@ eqs = [D(x) ~ 1] affect = [x ~ 0] affect_neg = [x ~ 1] -## Test SymbolicContinuousCallback @testset "SymbolicContinuousCallback constructors" begin e = SymbolicContinuousCallback(eqs[]) @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == NULL_AFFECT - @test e.affect_neg == NULL_AFFECT + @test isequal(equations(e), eqs) + @test e.affect === nothing + @test e.affect_neg === nothing @test e.rootfind == SciMLBase.LeftRootFind e = SymbolicContinuousCallback(eqs) @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == NULL_AFFECT - @test e.affect_neg == NULL_AFFECT + @test isequal(equations(e), eqs) + @test e.affect === nothing + @test e.affect_neg === nothing @test e.rootfind == SciMLBase.LeftRootFind - e = SymbolicContinuousCallback(eqs, NULL_AFFECT) + e = SymbolicContinuousCallback(eqs, nothing) @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == NULL_AFFECT - @test e.affect_neg == NULL_AFFECT + @test isequal(equations(e), eqs) + @test e.affect === nothing + @test e.affect_neg === nothing @test e.rootfind == SciMLBase.LeftRootFind - e = SymbolicContinuousCallback(eqs[], NULL_AFFECT) + e = SymbolicContinuousCallback(eqs[], nothing) @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == NULL_AFFECT - @test e.affect_neg == NULL_AFFECT + @test isequal(equations(e), eqs) + @test e.affect === nothing + @test e.affect_neg === nothing @test e.rootfind == SciMLBase.LeftRootFind - e = SymbolicContinuousCallback(eqs => NULL_AFFECT) + e = SymbolicContinuousCallback(eqs => nothing) @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == NULL_AFFECT - @test e.affect_neg == NULL_AFFECT + @test isequal(equations(e), eqs) + @test e.affect === nothing + @test e.affect_neg === nothing @test e.rootfind == SciMLBase.LeftRootFind - e = SymbolicContinuousCallback(eqs[] => NULL_AFFECT) + e = SymbolicContinuousCallback(eqs[] => nothing) @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == NULL_AFFECT - @test e.affect_neg == NULL_AFFECT + @test isequal(equations(e), eqs) + @test e.affect === nothing + @test e.affect_neg === nothing @test e.rootfind == SciMLBase.LeftRootFind ## With affect - - e = SymbolicContinuousCallback(eqs[], affect) - @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == affect - @test e.affect_neg == affect - @test e.rootfind == SciMLBase.LeftRootFind - - e = SymbolicContinuousCallback(eqs, affect) - @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == affect - @test e.affect_neg == affect - @test e.rootfind == SciMLBase.LeftRootFind - - e = SymbolicContinuousCallback(eqs, affect) - @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == affect - @test e.affect_neg == affect - @test e.rootfind == SciMLBase.LeftRootFind - e = SymbolicContinuousCallback(eqs[], affect) @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == affect - @test e.affect_neg == affect - @test e.rootfind == SciMLBase.LeftRootFind - - e = SymbolicContinuousCallback(eqs => affect) - @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == affect - @test e.affect_neg == affect - @test e.rootfind == SciMLBase.LeftRootFind - - e = SymbolicContinuousCallback(eqs[] => affect) - @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == affect - @test e.affect_neg == affect + @test isequal(equations(e), eqs) @test e.rootfind == SciMLBase.LeftRootFind # with only positive edge affect - e = SymbolicContinuousCallback(eqs[], affect, affect_neg = nothing) @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == affect - @test isnothing(e.affect_neg) - @test e.rootfind == SciMLBase.LeftRootFind - - e = SymbolicContinuousCallback(eqs, affect, affect_neg = nothing) - @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == affect - @test isnothing(e.affect_neg) - @test e.rootfind == SciMLBase.LeftRootFind - - e = SymbolicContinuousCallback(eqs, affect, affect_neg = nothing) - @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == affect - @test isnothing(e.affect_neg) - @test e.rootfind == SciMLBase.LeftRootFind - - e = SymbolicContinuousCallback(eqs[], affect, affect_neg = nothing) - @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == affect + @test isequal(equations(e), eqs) @test isnothing(e.affect_neg) @test e.rootfind == SciMLBase.LeftRootFind # with explicit edge affects - - e = SymbolicContinuousCallback(eqs[], affect, affect_neg = affect_neg) - @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == affect - @test e.affect_neg == affect_neg - @test e.rootfind == SciMLBase.LeftRootFind - - e = SymbolicContinuousCallback(eqs, affect, affect_neg = affect_neg) - @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == affect - @test e.affect_neg == affect_neg - @test e.rootfind == SciMLBase.LeftRootFind - - e = SymbolicContinuousCallback(eqs, affect, affect_neg = affect_neg) - @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == affect - @test e.affect_neg == affect_neg - @test e.rootfind == SciMLBase.LeftRootFind - e = SymbolicContinuousCallback(eqs[], affect, affect_neg = affect_neg) @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == affect - @test e.affect_neg == affect_neg + @test isequal(equations(e), eqs) @test e.rootfind == SciMLBase.LeftRootFind # with different root finding ops - e = SymbolicContinuousCallback( eqs[], affect, affect_neg = affect_neg, rootfind = SciMLBase.LeftRootFind) @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == affect - @test e.affect_neg == affect_neg + @test isequal(equations(e), eqs) @test e.rootfind == SciMLBase.LeftRootFind - - e = SymbolicContinuousCallback( - eqs[], affect, affect_neg = affect_neg, rootfind = SciMLBase.RightRootFind) - @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == affect - @test e.affect_neg == affect_neg - @test e.rootfind == SciMLBase.RightRootFind - - e = SymbolicContinuousCallback( - eqs[], affect, affect_neg = affect_neg, rootfind = SciMLBase.NoRootFind) - @test e isa SymbolicContinuousCallback - @test isequal(e.eqs, eqs) - @test e.affect == affect - @test e.affect_neg == affect_neg - @test e.rootfind == SciMLBase.NoRootFind - # test plural constructor - - e = SymbolicContinuousCallbacks(eqs[]) - @test e isa Vector{SymbolicContinuousCallback} - @test isequal(e[].eqs, eqs) - @test e[].affect == NULL_AFFECT - - e = SymbolicContinuousCallbacks(eqs) - @test e isa Vector{SymbolicContinuousCallback} - @test isequal(e[].eqs, eqs) - @test e[].affect == NULL_AFFECT - - e = SymbolicContinuousCallbacks(eqs[] => affect) - @test e isa Vector{SymbolicContinuousCallback} - @test isequal(e[].eqs, eqs) - @test e[].affect == affect - - e = SymbolicContinuousCallbacks(eqs => affect) - @test e isa Vector{SymbolicContinuousCallback} - @test isequal(e[].eqs, eqs) - @test e[].affect == affect - - e = SymbolicContinuousCallbacks([eqs[] => affect]) - @test e isa Vector{SymbolicContinuousCallback} - @test isequal(e[].eqs, eqs) - @test e[].affect == affect - - e = SymbolicContinuousCallbacks([eqs => affect]) - @test e isa Vector{SymbolicContinuousCallback} - @test isequal(e[].eqs, eqs) - @test e[].affect == affect - - e = SymbolicContinuousCallbacks(SymbolicContinuousCallbacks([eqs => affect])) - @test e isa Vector{SymbolicContinuousCallback} - @test isequal(e[].eqs, eqs) - @test e[].affect == affect end @testset "ImperativeAffect constructors" begin @@ -341,409 +201,357 @@ end @test m.ctx === 3 end -## - -@named sys = ODESystem(eqs, t, continuous_events = [x ~ 1]) -@test getfield(sys, :continuous_events)[] == - SymbolicContinuousCallback(Equation[x ~ 1], NULL_AFFECT) -@test isequal(equations(getfield(sys, :continuous_events))[], x ~ 1) -fsys = flatten(sys) -@test isequal(equations(getfield(fsys, :continuous_events))[], x ~ 1) - -@named sys2 = ODESystem([D(x) ~ 1], t, continuous_events = [x ~ 2], systems = [sys]) -@test getfield(sys2, :continuous_events)[] == - SymbolicContinuousCallback(Equation[x ~ 2], NULL_AFFECT) -@test all(ModelingToolkit.continuous_events(sys2) .== [ - SymbolicContinuousCallback(Equation[x ~ 2], NULL_AFFECT), - SymbolicContinuousCallback(Equation[sys.x ~ 1], NULL_AFFECT) -]) - -@test isequal(equations(getfield(sys2, :continuous_events))[1], x ~ 2) -@test length(ModelingToolkit.continuous_events(sys2)) == 2 -@test isequal(ModelingToolkit.continuous_events(sys2)[1].eqs[], x ~ 2) -@test isequal(ModelingToolkit.continuous_events(sys2)[2].eqs[], sys.x ~ 1) - -sys = complete(sys) -sys_nosplit = complete(sys; split = false) -sys2 = complete(sys2) -# Functions should be generated for root-finding equations -prob = ODEProblem(sys, Pair[], (0.0, 2.0)) -p0 = 0 -t0 = 0 -@test get_callback(prob) isa ModelingToolkit.DiffEqCallbacks.ContinuousCallback -cb = ModelingToolkit.generate_rootfinding_callback(sys) -cond = cb.condition -out = [0.0] -cond.rf_ip(out, [0], p0, t0) -@test out[] ≈ -1 # signature is u,p,t -cond.rf_ip(out, [1], p0, t0) -@test out[] ≈ 0 # signature is u,p,t -cond.rf_ip(out, [2], p0, t0) -@test out[] ≈ 1 # signature is u,p,t - -prob = ODEProblem(sys, Pair[], (0.0, 2.0)) -prob_nosplit = ODEProblem(sys_nosplit, Pair[], (0.0, 2.0)) -sol = solve(prob, Tsit5()) -sol_nosplit = solve(prob_nosplit, Tsit5()) -@test minimum(t -> abs(t - 1), sol.t) < 1e-10 # test that the solver stepped at the root -@test minimum(t -> abs(t - 1), sol_nosplit.t) < 1e-10 # test that the solver stepped at the root - -# Test that a user provided callback is respected -test_callback = DiscreteCallback(x -> x, x -> x) -prob = ODEProblem(sys, Pair[], (0.0, 2.0), callback = test_callback) -prob_nosplit = ODEProblem(sys_nosplit, Pair[], (0.0, 2.0), callback = test_callback) -cbs = get_callback(prob) -cbs_nosplit = get_callback(prob_nosplit) -@test cbs isa CallbackSet -@test cbs.discrete_callbacks[1] == test_callback -@test cbs_nosplit isa CallbackSet -@test cbs_nosplit.discrete_callbacks[1] == test_callback - -prob = ODEProblem(sys2, Pair[], (0.0, 3.0)) -cb = get_callback(prob) -@test cb isa ModelingToolkit.DiffEqCallbacks.VectorContinuousCallback - -cond = cb.condition -out = [0.0, 0.0] -# the root to find is 2 -cond.rf_ip(out, [0, 0], p0, t0) -@test out[1] ≈ -2 # signature is u,p,t -cond.rf_ip(out, [1, 0], p0, t0) -@test out[1] ≈ -1 # signature is u,p,t -cond.rf_ip(out, [2, 0], p0, t0) # this should return 0 -@test out[1] ≈ 0 # signature is u,p,t - -# the root to find is 1 -out = [0.0, 0.0] -cond.rf_ip(out, [0, 0], p0, t0) -@test out[2] ≈ -1 # signature is u,p,t -cond.rf_ip(out, [0, 1], p0, t0) # this should return 0 -@test out[2] ≈ 0 # signature is u,p,t -cond.rf_ip(out, [0, 2], p0, t0) -@test out[2] ≈ 1 # signature is u,p,t - -sol = solve(prob, Tsit5(); abstol = 1e-14, reltol = 1e-14) -@test minimum(t -> abs(t - 1), sol.t) < 1e-10 # test that the solver stepped at the first root -@test minimum(t -> abs(t - 2), sol.t) < 1e-10 # test that the solver stepped at the second root - -@named sys = ODESystem(eqs, t, continuous_events = [x ~ 1, x ~ 2]) # two root eqs using the same unknown -sys = complete(sys) -prob = ODEProblem(sys, Pair[], (0.0, 3.0)) -@test get_callback(prob) isa ModelingToolkit.DiffEqCallbacks.VectorContinuousCallback -sol = solve(prob, Tsit5(); abstol = 1e-14, reltol = 1e-14) -@test minimum(t -> abs(t - 1), sol.t) < 1e-10 # test that the solver stepped at the first root -@test minimum(t -> abs(t - 2), sol.t) < 1e-10 # test that the solver stepped at the second root - -## Test bouncing ball with equation affect -@variables x(t)=1 v(t)=0 - -root_eqs = [x ~ 0] -affect = [v ~ -v] - -@named ball = ODESystem([D(x) ~ v - D(v) ~ -9.8], t, continuous_events = root_eqs => affect) - -@test getfield(ball, :continuous_events)[] == - SymbolicContinuousCallback(Equation[x ~ 0], Equation[v ~ -v]) -ball = structural_simplify(ball) - -@test length(ModelingToolkit.continuous_events(ball)) == 1 - -tspan = (0.0, 5.0) -prob = ODEProblem(ball, Pair[], tspan) -sol = solve(prob, Tsit5()) -@test 0 <= minimum(sol[x]) <= 1e-10 # the ball never went through the floor but got very close -# plot(sol) - -## Test bouncing ball in 2D with walls -@variables x(t)=1 y(t)=0 vx(t)=0 vy(t)=1 - -continuous_events = [[x ~ 0] => [vx ~ -vx] - [y ~ -1.5, y ~ 1.5] => [vy ~ -vy]] - -@named ball = ODESystem( - [D(x) ~ vx - D(y) ~ vy - D(vx) ~ -9.8 - D(vy) ~ -0.01vy], t; continuous_events) - -_ball = ball -ball = structural_simplify(_ball) -ball_nosplit = structural_simplify(_ball; split = false) - -tspan = (0.0, 5.0) -prob = ODEProblem(ball, Pair[], tspan) -prob_nosplit = ODEProblem(ball_nosplit, Pair[], tspan) - -cb = get_callback(prob) -@test cb isa ModelingToolkit.DiffEqCallbacks.VectorContinuousCallback -@test getfield(ball, :continuous_events)[1] == - SymbolicContinuousCallback(Equation[x ~ 0], Equation[vx ~ -vx]) -@test getfield(ball, :continuous_events)[2] == - SymbolicContinuousCallback(Equation[y ~ -1.5, y ~ 1.5], Equation[vy ~ -vy]) -cond = cb.condition -out = [0.0, 0.0, 0.0] -cond.rf_ip(out, [0, 0, 0, 0], p0, t0) -@test out ≈ [0, 1.5, -1.5] - -sol = solve(prob, Tsit5()) -sol_nosplit = solve(prob_nosplit, Tsit5()) -@test 0 <= minimum(sol[x]) <= 1e-10 # the ball never went through the floor but got very close -@test minimum(sol[y]) ≈ -1.5 # check wall conditions -@test maximum(sol[y]) ≈ 1.5 # check wall conditions -@test 0 <= minimum(sol_nosplit[x]) <= 1e-10 # the ball never went through the floor but got very close -@test minimum(sol_nosplit[y]) ≈ -1.5 # check wall conditions -@test maximum(sol_nosplit[y]) ≈ 1.5 # check wall conditions - -# tv = sort([LinRange(0, 5, 200); sol.t]) -# plot(sol(tv)[y], sol(tv)[x], line_z=tv) -# vline!([-1.5, 1.5], l=(:black, 5), primary=false) -# hline!([0], l=(:black, 5), primary=false) - -## Test multi-variable affect -# in this test, there are two variables affected by a single event. -continuous_events = [ - [x ~ 0] => [vx ~ -vx, vy ~ -vy] -] - -@named ball = ODESystem([D(x) ~ vx - D(y) ~ vy - D(vx) ~ -1 - D(vy) ~ 0], t; continuous_events) - -ball_nosplit = structural_simplify(ball) -ball = structural_simplify(ball) - -tspan = (0.0, 5.0) -prob = ODEProblem(ball, Pair[], tspan) -prob_nosplit = ODEProblem(ball_nosplit, Pair[], tspan) -sol = solve(prob, Tsit5()) -sol_nosplit = solve(prob_nosplit, Tsit5()) -@test 0 <= minimum(sol[x]) <= 1e-10 # the ball never went through the floor but got very close -@test -minimum(sol[y]) ≈ maximum(sol[y]) ≈ sqrt(2) # the ball will never go further than √2 in either direction (gravity was changed to 1 to get this particular number) -@test 0 <= minimum(sol_nosplit[x]) <= 1e-10 # the ball never went through the floor but got very close -@test -minimum(sol_nosplit[y]) ≈ maximum(sol_nosplit[y]) ≈ sqrt(2) # the ball will never go further than √2 in either direction (gravity was changed to 1 to get this particular number) - -# tv = sort([LinRange(0, 5, 200); sol.t]) -# plot(sol(tv)[y], sol(tv)[x], line_z=tv) -# vline!([-1.5, 1.5], l=(:black, 5), primary=false) -# hline!([0], l=(:black, 5), primary=false) - -# issue https://github.com/SciML/ModelingToolkit.jl/issues/1386 -# tests that it works for ODAESystem -@variables vs(t) v(t) vmeasured(t) -eq = [vs ~ sin(2pi * t) - D(v) ~ vs - v - D(vmeasured) ~ 0.0] -ev = [sin(20pi * t) ~ 0.0] => [vmeasured ~ v] -@named sys = ODESystem(eq, t, continuous_events = ev) -sys = structural_simplify(sys) -prob = ODEProblem(sys, zeros(2), (0.0, 5.1)) -sol = solve(prob, Tsit5()) -@test all(minimum((0:0.05:5) .- sol.t', dims = 2) .< 0.0001) # test that the solver stepped every 0.05s as dictated by event -@test sol([0.25 - eps()])[vmeasured][] == sol([0.23])[vmeasured][] # test the hold property - -## https://github.com/SciML/ModelingToolkit.jl/issues/1528 -Dₜ = D +@testset "Condition Compilation" begin + @named sys = ODESystem(eqs, t, continuous_events = [x ~ 1]) + @test getfield(sys, :continuous_events)[] == + SymbolicContinuousCallback(Equation[x ~ 1], nothing) + @test isequal(equations(getfield(sys, :continuous_events))[], x ~ 1) + fsys = flatten(sys) + @test isequal(equations(getfield(fsys, :continuous_events))[], x ~ 1) + + @named sys2 = ODESystem([D(x) ~ 1], t, continuous_events = [x ~ 2], systems = [sys]) + @test getfield(sys2, :continuous_events)[] == + SymbolicContinuousCallback(Equation[x ~ 2], nothing) + @test all(ModelingToolkit.continuous_events(sys2) .== [ + SymbolicContinuousCallback(Equation[x ~ 2], nothing), + SymbolicContinuousCallback(Equation[sys.x ~ 1], nothing) + ]) + + @test isequal(equations(getfield(sys2, :continuous_events))[1], x ~ 2) + @test length(ModelingToolkit.continuous_events(sys2)) == 2 + @test isequal(equations(ModelingToolkit.continuous_events(sys2)[1])[], x ~ 2) + @test isequal(equations(ModelingToolkit.continuous_events(sys2)[2])[], sys.x ~ 1) + + sys = complete(sys) + sys_nosplit = complete(sys; split = false) + sys2 = complete(sys2) + + # Test proper rootfinding + prob = ODEProblem(sys, Pair[], (0.0, 2.0)) + p0 = 0 + t0 = 0 + @test get_callback(prob) isa ModelingToolkit.DiffEqCallbacks.ContinuousCallback + cb = ModelingToolkit.generate_continuous_callbacks(sys) + cond = cb.condition + out = [0.0] + cond.f_iip(out, [0], p0, t0) + @test out[] ≈ -1 # signature is u,p,t + cond.f_iip(out, [1], p0, t0) + @test out[] ≈ 0 # signature is u,p,t + cond.f_iip(out, [2], p0, t0) + @test out[] ≈ 1 # signature is u,p,t + + prob = ODEProblem(sys, Pair[], (0.0, 2.0)) + prob_nosplit = ODEProblem(sys_nosplit, Pair[], (0.0, 2.0)) + sol = solve(prob, Tsit5()) + sol_nosplit = solve(prob_nosplit, Tsit5()) + @test minimum(t -> abs(t - 1), sol.t) < 1e-10 # test that the solver stepped at the root + @test minimum(t -> abs(t - 1), sol_nosplit.t) < 1e-10 # test that the solver stepped at the root + + # Test user-provided callback is respected + test_callback = DiscreteCallback(x -> x, x -> x) + prob = ODEProblem(sys, Pair[], (0.0, 2.0), callback = test_callback) + prob_nosplit = ODEProblem(sys_nosplit, Pair[], (0.0, 2.0), callback = test_callback) + cbs = get_callback(prob) + cbs_nosplit = get_callback(prob_nosplit) + @test cbs isa CallbackSet + @test cbs.discrete_callbacks[1] == test_callback + @test cbs_nosplit isa CallbackSet + @test cbs_nosplit.discrete_callbacks[1] == test_callback + + prob = ODEProblem(sys2, Pair[], (0.0, 3.0)) + cb = get_callback(prob) + @test cb isa ModelingToolkit.DiffEqCallbacks.VectorContinuousCallback + + cond = cb.condition + out = [0.0, 0.0] + # the root to find is 2 + cond.f_iip(out, [0, 0], p0, t0) + @test out[1] ≈ -2 # signature is u,p,t + cond.f_iip(out, [1, 0], p0, t0) + @test out[1] ≈ -1 # signature is u,p,t + cond.f_iip(out, [2, 0], p0, t0) # this should return 0 + @test out[1] ≈ 0 # signature is u,p,t + + # the root to find is 1 + out = [0.0, 0.0] + cond.f_iip(out, [0, 0], p0, t0) + @test out[2] ≈ -1 # signature is u,p,t + cond.f_iip(out, [0, 1], p0, t0) # this should return 0 + @test out[2] ≈ 0 # signature is u,p,t + cond.f_iip(out, [0, 2], p0, t0) + @test out[2] ≈ 1 # signature is u,p,t -@parameters u(t) [input = true] # Indicate that this is a controlled input -@parameters y(t) [output = true] # Indicate that this is a measured output + sol = solve(prob, Tsit5()) + @test minimum(t -> abs(t - 1), sol.t) < 1e-9 # test that the solver stepped at the first root + @test minimum(t -> abs(t - 2), sol.t) < 1e-9 # test that the solver stepped at the second root -function Mass(; name, m = 1.0, p = 0, v = 0) - ps = @parameters m = m - sts = @variables pos(t)=p vel(t)=v - eqs = Dₜ(pos) ~ vel - ODESystem(eqs, t, [pos, vel], ps; name) -end -function Spring(; name, k = 1e4) - ps = @parameters k = k - @variables x(t) = 0 # Spring deflection - ODESystem(Equation[], t, [x], ps; name) -end -function Damper(; name, c = 10) - ps = @parameters c = c - @variables vel(t) = 0 - ODESystem(Equation[], t, [vel], ps; name) -end -function SpringDamper(; name, k = false, c = false) - spring = Spring(; name = :spring, k) - damper = Damper(; name = :damper, c) - compose(ODESystem(Equation[], t; name), - spring, damper) -end -connect_sd(sd, m1, m2) = [sd.spring.x ~ m1.pos - m2.pos, sd.damper.vel ~ m1.vel - m2.vel] -sd_force(sd) = -sd.spring.k * sd.spring.x - sd.damper.c * sd.damper.vel -@named mass1 = Mass(; m = 1) -@named mass2 = Mass(; m = 1) -@named sd = SpringDamper(; k = 1000, c = 10) -function Model(u, d = 0) - eqs = [connect_sd(sd, mass1, mass2) - Dₜ(mass1.vel) ~ (sd_force(sd) + u) / mass1.m - Dₜ(mass2.vel) ~ (-sd_force(sd) + d) / mass2.m] - @named _model = ODESystem(eqs, t; observed = [y ~ mass2.pos]) - @named model = compose(_model, mass1, mass2, sd) + @named sys = ODESystem(eqs, t, continuous_events = [x ~ 1, x ~ 2]) # two root eqs using the same unknown + sys = complete(sys) + prob = ODEProblem(sys, Pair[], (0.0, 3.0)) + @test get_callback(prob) isa ModelingToolkit.DiffEqCallbacks.VectorContinuousCallback + sol = solve(prob, Tsit5()) + @test minimum(t -> abs(t - 1), sol.t) < 1e-9 # test that the solver stepped at the first root + @test minimum(t -> abs(t - 2), sol.t) < 1e-9 # test that the solver stepped at the second root end -model = Model(sin(30t)) -sys = structural_simplify(model) -@test isempty(ModelingToolkit.continuous_events(sys)) -let - function testsol(osys, u0, p, tspan; tstops = Float64[], paramtotest = nothing, - kwargs...) - oprob = ODEProblem(complete(osys), u0, tspan, p; kwargs...) - sol = solve(oprob, Tsit5(); tstops = tstops, abstol = 1e-10, reltol = 1e-10) - @test isapprox(sol(1.0000000001)[1] - sol(0.999999999)[1], 1.0; rtol = 1e-6) - paramtotest === nothing || (@test sol.ps[paramtotest] == 1.0) - @test isapprox(sol(4.0)[1], 2 * exp(-2.0)) - sol - end +@testset "Bouncing Ball" begin + ###### 1D Bounce + @variables x(t)=1 v(t)=0 - @parameters k t1 t2 - @variables A(t) B(t) + root_eqs = [x ~ 0] + affect = [v ~ -Pre(v)] - cond1 = (t == t1) - affect1 = [A ~ A + 1] - cb1 = cond1 => affect1 - cond2 = (t == t2) - affect2 = [k ~ 1.0] - cb2 = cond2 => affect2 + @named ball = ODESystem( + [D(x) ~ v + D(v) ~ -9.8], t, continuous_events = root_eqs => affect) - ∂ₜ = D - eqs = [∂ₜ(A) ~ -k * A] - @named osys = ODESystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb1, cb2]) - u0 = [A => 1.0] - p = [k => 0.0, t1 => 1.0, t2 => 2.0] - tspan = (0.0, 4.0) - testsol(osys, u0, p, tspan; tstops = [1.0, 2.0], paramtotest = k) + @test only(continuous_events(ball)) == + SymbolicContinuousCallback(Equation[x ~ 0], Equation[v ~ -Pre(v)]) + ball = structural_simplify(ball) - cond1a = (t == t1) - affect1a = [A ~ A + 1, B ~ A] - cb1a = cond1a => affect1a - @named osys1 = ODESystem(eqs, t, [A, B], [k, t1, t2], discrete_events = [cb1a, cb2]) - u0′ = [A => 1.0, B => 0.0] - sol = testsol( - osys1, u0′, p, tspan; tstops = [1.0, 2.0], check_length = false, paramtotest = k) - @test sol(1.0000001, idxs = B) == 2.0 + @test length(ModelingToolkit.continuous_events(ball)) == 1 - # same as above - but with set-time event syntax - cb1‵ = [1.0] => affect1 # needs to be a Vector for the event to happen only once - cb2‵ = [2.0] => affect2 - @named osys‵ = ODESystem(eqs, t, [A], [k], discrete_events = [cb1‵, cb2‵]) - testsol(osys‵, u0, p, tspan; paramtotest = k) + tspan = (0.0, 5.0) + prob = ODEProblem(ball, Pair[], tspan) + sol = solve(prob, Tsit5()) + @test 0 <= minimum(sol[x]) <= 1e-10 # the ball never went through the floor but got very close + + ###### 2D bouncing ball + @variables x(t)=1 y(t)=0 vx(t)=0 vy(t)=1 + + events = [[x ~ 0] => [vx ~ -Pre(vx)] + [y ~ -1.5, y ~ 1.5] => [vy ~ -Pre(vy)]] + + @named ball = ODESystem( + [D(x) ~ vx + D(y) ~ vy + D(vx) ~ -9.8 + D(vy) ~ -0.01vy], t; continuous_events = events) + + _ball = ball + ball = structural_simplify(_ball) + ball_nosplit = structural_simplify(_ball; split = false) + + tspan = (0.0, 5.0) + prob = ODEProblem(ball, Pair[], tspan) + prob_nosplit = ODEProblem(ball_nosplit, Pair[], tspan) + + cb = get_callback(prob) + @test cb isa ModelingToolkit.DiffEqCallbacks.VectorContinuousCallback + @test getfield(ball, :continuous_events)[1] == + SymbolicContinuousCallback(Equation[x ~ 0], Equation[vx ~ -Pre(vx)]) + @test getfield(ball, :continuous_events)[2] == + SymbolicContinuousCallback(Equation[y ~ -1.5, y ~ 1.5], Equation[vy ~ -Pre(vy)]) + cond = cb.condition + out = [0.0, 0.0, 0.0] + p0 = 0.0 + t0 = 0.0 + cond.f_iip(out, [0, 0, 0, 0], p0, t0) + @test out ≈ [0, 1.5, -1.5] - # mixing discrete affects - @named osys3 = ODESystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb1, cb2‵]) - testsol(osys3, u0, p, tspan; tstops = [1.0], paramtotest = k) + sol = solve(prob, Tsit5()) + sol_nosplit = solve(prob_nosplit, Tsit5()) + @test 0 <= minimum(sol[x]) <= 1e-10 # the ball never went through the floor but got very close + @test minimum(sol[y]) ≈ -1.5 # check wall conditions + @test maximum(sol[y]) ≈ 1.5 # check wall conditions + @test 0 <= minimum(sol_nosplit[x]) <= 1e-10 # the ball never went through the floor but got very close + @test minimum(sol_nosplit[y]) ≈ -1.5 # check wall conditions + @test maximum(sol_nosplit[y]) ≈ 1.5 # check wall conditions + + ## Test multi-variable affect + # in this test, there are two variables affected by a single event. + events = [[x ~ 0] => [vx ~ -Pre(vx), vy ~ -Pre(vy)]] + + @named ball = ODESystem( + [D(x) ~ vx + D(y) ~ vy + D(vx) ~ -1 + D(vy) ~ 0], t; continuous_events = events) + + ball_nosplit = structural_simplify(ball) + ball = structural_simplify(ball) + + tspan = (0.0, 5.0) + prob = ODEProblem(ball, Pair[], tspan) + prob_nosplit = ODEProblem(ball_nosplit, Pair[], tspan) + sol = solve(prob, Tsit5()) + sol_nosplit = solve(prob_nosplit, Tsit5()) + @test 0 <= minimum(sol[x]) <= 1e-10 # the ball never went through the floor but got very close + @test -minimum(sol[y]) ≈ maximum(sol[y]) ≈ sqrt(2) # the ball will never go further than √2 in either direction (gravity was changed to 1 to get this particular number) + @test 0 <= minimum(sol_nosplit[x]) <= 1e-10 # the ball never went through the floor but got very close + @test -minimum(sol_nosplit[y]) ≈ maximum(sol_nosplit[y]) ≈ sqrt(2) # the ball will never go further than √2 in either direction (gravity was changed to 1 to get this particular number) +end - # mixing with a func affect - function affect!(integrator, u, p, ctx) - integrator.ps[p.k] = 1.0 - nothing - end - cb2‵‵ = [2.0] => (affect!, [], [k], [k], nothing) - @named osys4 = ODESystem(eqs, t, [A], [k, t1], discrete_events = [cb1, cb2‵‵]) - oprob4 = ODEProblem(complete(osys4), u0, tspan, p) - testsol(osys4, u0, p, tspan; tstops = [1.0], paramtotest = k) +# issue https://github.com/SciML/ModelingToolkit.jl/issues/1386 +# tests that it works for ODAESystem +@testset "ODAESystem" begin + @variables vs(t) v(t) vmeasured(t) + eq = [vs ~ sin(2pi * t) + D(v) ~ vs - v + D(vmeasured) ~ 0.0] + ev = [sin(20pi * t) ~ 0.0] => [vmeasured ~ Pre(v)] + @named sys = ODESystem(eq, t, continuous_events = ev) + sys = structural_simplify(sys) + prob = ODEProblem(sys, zeros(2), (0.0, 5.1)) + sol = solve(prob, Tsit5()) + @test all(minimum((0:0.1:5) .- sol.t', dims = 2) .< 0.0001) # test that the solver stepped every 0.1s as dictated by event + @test sol([0.25 - eps()])[vmeasured][] == sol([0.23])[vmeasured][] # test the hold property +end - # mixing with symbolic condition in the func affect - cb2‵‵‵ = (t == t2) => (affect!, [], [k], [k], nothing) - @named osys5 = ODESystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb1, cb2‵‵‵]) - testsol(osys5, u0, p, tspan; tstops = [1.0, 2.0]) - @named osys6 = ODESystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb2‵‵‵, cb1]) - testsol(osys6, u0, p, tspan; tstops = [1.0, 2.0], paramtotest = k) +## https://github.com/SciML/ModelingToolkit.jl/issues/1528 +@testset "Handle Empty Events" begin + Dₜ = D - # mix a continuous event too - cond3 = A ~ 0.1 - affect3 = [k ~ 0.0] - cb3 = cond3 => affect3 - @named osys7 = ODESystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb1, cb2‵‵‵], - continuous_events = [cb3]) - sol = testsol(osys7, u0, p, (0.0, 10.0); tstops = [1.0, 2.0]) - @test isapprox(sol(10.0)[1], 0.1; atol = 1e-10, rtol = 1e-10) + @parameters u(t) [input = true] # Indicate that this is a controlled input + @parameters y(t) [output = true] # Indicate that this is a measured output + + function Mass(; name, m = 1.0, p = 0, v = 0) + ps = @parameters m = m + sts = @variables pos(t)=p vel(t)=v + eqs = Dₜ(pos) ~ vel + ODESystem(eqs, t, [pos, vel], ps; name) + end + function Spring(; name, k = 1e4) + ps = @parameters k = k + @variables x(t) = 0 # Spring deflection + ODESystem(Equation[], t, [x], ps; name) + end + function Damper(; name, c = 10) + ps = @parameters c = c + @variables vel(t) = 0 + ODESystem(Equation[], t, [vel], ps; name) + end + function SpringDamper(; name, k = false, c = false) + spring = Spring(; name = :spring, k) + damper = Damper(; name = :damper, c) + compose(ODESystem(Equation[], t; name), + spring, damper) + end + connect_sd(sd, m1, m2) = [ + sd.spring.x ~ m1.pos - m2.pos, sd.damper.vel ~ m1.vel - m2.vel] + sd_force(sd) = -sd.spring.k * sd.spring.x - sd.damper.c * sd.damper.vel + @named mass1 = Mass(; m = 1) + @named mass2 = Mass(; m = 1) + @named sd = SpringDamper(; k = 1000, c = 10) + function Model(u, d = 0) + eqs = [connect_sd(sd, mass1, mass2) + Dₜ(mass1.vel) ~ (sd_force(sd) + u) / mass1.m + Dₜ(mass2.vel) ~ (-sd_force(sd) + d) / mass2.m] + @named _model = ODESystem(eqs, t; observed = [y ~ mass2.pos]) + @named model = compose(_model, mass1, mass2, sd) + end + model = Model(sin(30t)) + sys = structural_simplify(model) + @test isempty(ModelingToolkit.continuous_events(sys)) end -let - function testsol(ssys, u0, p, tspan; tstops = Float64[], paramtotest = nothing, +@testset "SDE/ODESystem Discrete Callbacks" begin + function testsol( + sys, probtype, solver, u0, p, tspan; tstops = Float64[], paramtotest = nothing, kwargs...) - sprob = SDEProblem(complete(ssys), u0, tspan, p; kwargs...) - sol = solve(sprob, RI5(); tstops = tstops, abstol = 1e-10, reltol = 1e-10) - @test isapprox(sol(1.0000000001)[1] - sol(0.999999999)[1], 1.0; rtol = 1e-4) - paramtotest === nothing || (@test sol.ps[paramtotest] == 1.0) - @test isapprox(sol(4.0)[1], 2 * exp(-2.0), atol = 1e-4) + prob = probtype(complete(sys), u0, tspan, p; kwargs...) + sol = solve(prob, solver(); tstops = tstops, abstol = 1e-10, reltol = 1e-10) + @test isapprox(sol(1.0000000001)[1] - sol(0.999999999)[1], 1.0; rtol = 1e-6) + paramtotest === nothing || (@test sol.ps[paramtotest] == [0.0, 1.0]) + @test isapprox(sol(4.0)[1], 2 * exp(-2.0); rtol = 1e-6) sol end - @parameters k t1 t2 + @parameters k(t) t1 t2 @variables A(t) B(t) cond1 = (t == t1) - affect1 = [A ~ A + 1] + affect1 = [A ~ Pre(A) + 1] cb1 = cond1 => affect1 cond2 = (t == t2) affect2 = [k ~ 1.0] cb2 = cond2 => affect2 + cb2 = SymbolicDiscreteCallback(cb2, discrete_parameters = [k], iv = t) ∂ₜ = D eqs = [∂ₜ(A) ~ -k * A] - @named ssys = SDESystem(eqs, [0.0], t, [A], [k, t1, t2], - discrete_events = [cb1, cb2]) + @named osys = ODESystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb1, cb2]) + @named ssys = SDESystem(eqs, [0.0], t, [A], [k, t1, t2], discrete_events = [cb1, cb2]) u0 = [A => 1.0] p = [k => 0.0, t1 => 1.0, t2 => 2.0] tspan = (0.0, 4.0) - testsol(ssys, u0, p, tspan; tstops = [1.0, 2.0], paramtotest = k) + testsol(osys, ODEProblem, Tsit5, u0, p, tspan; tstops = [1.0, 2.0], paramtotest = k) + testsol(ssys, SDEProblem, RI5, u0, p, tspan; tstops = [1.0, 2.0], paramtotest = k) cond1a = (t == t1) - affect1a = [A ~ A + 1, B ~ A] + affect1a = [A ~ Pre(A) + 1, B ~ A] cb1a = cond1a => affect1a + @named osys1 = ODESystem(eqs, t, [A, B], [k, t1, t2], discrete_events = [cb1a, cb2]) @named ssys1 = SDESystem(eqs, [0.0], t, [A, B], [k, t1, t2], discrete_events = [cb1a, cb2]) u0′ = [A => 1.0, B => 0.0] - sol = testsol( - ssys1, u0′, p, tspan; tstops = [1.0, 2.0], check_length = false, paramtotest = k) - @test sol(1.0000001, idxs = 2) == 2.0 + sol = testsol(osys1, ODEProblem, Tsit5, u0′, p, tspan; + tstops = [1.0, 2.0], check_length = false, paramtotest = k) + @test sol(1.0000001, idxs = B) == 2.0 + + sol = testsol(ssys1, SDEProblem, RI5, u0′, p, tspan; tstops = [1.0, 2.0], + check_length = false, paramtotest = k) + @test sol(1.0000001, idxs = B) == 2.0 # same as above - but with set-time event syntax cb1‵ = [1.0] => affect1 # needs to be a Vector for the event to happen only once - cb2‵ = [2.0] => affect2 + cb2‵ = SymbolicDiscreteCallback([2.0] => affect2, discrete_parameters = [k], iv = t) + @named osys‵ = ODESystem(eqs, t, [A], [k], discrete_events = [cb1‵, cb2‵]) @named ssys‵ = SDESystem(eqs, [0.0], t, [A], [k], discrete_events = [cb1‵, cb2‵]) - testsol(ssys‵, u0, p, tspan; paramtotest = k) + testsol(osys‵, ODEProblem, Tsit5, u0, p, tspan; paramtotest = k) + testsol(ssys‵, SDEProblem, RI5, u0, p, tspan; paramtotest = k) # mixing discrete affects + @named osys3 = ODESystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb1, cb2‵]) @named ssys3 = SDESystem(eqs, [0.0], t, [A], [k, t1, t2], discrete_events = [cb1, cb2‵]) - testsol(ssys3, u0, p, tspan; tstops = [1.0], paramtotest = k) + testsol(osys3, ODEProblem, Tsit5, u0, p, tspan; tstops = [1.0], paramtotest = k) + testsol(ssys3, SDEProblem, RI5, u0, p, tspan; tstops = [1.0], paramtotest = k) # mixing with a func affect function affect!(integrator, u, p, ctx) - setp(integrator, p.k)(integrator, 1.0) + integrator.ps[p.k] = 1.0 nothing end cb2‵‵ = [2.0] => (affect!, [], [k], [k], nothing) + @named osys4 = ODESystem(eqs, t, [A], [k, t1], discrete_events = [cb1, cb2‵‵]) @named ssys4 = SDESystem(eqs, [0.0], t, [A], [k, t1], discrete_events = [cb1, cb2‵‵]) - testsol(ssys4, u0, p, tspan; tstops = [1.0], paramtotest = k) + oprob4 = ODEProblem(complete(osys4), u0, tspan, p) + testsol(osys4, ODEProblem, Tsit5, u0, p, tspan; tstops = [1.0], paramtotest = k) + testsol(ssys4, SDEProblem, RI5, u0, p, tspan; tstops = [1.0], paramtotest = k) # mixing with symbolic condition in the func affect cb2‵‵‵ = (t == t2) => (affect!, [], [k], [k], nothing) + @named osys5 = ODESystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb1, cb2‵‵‵]) @named ssys5 = SDESystem(eqs, [0.0], t, [A], [k, t1, t2], discrete_events = [cb1, cb2‵‵‵]) - testsol(ssys5, u0, p, tspan; tstops = [1.0, 2.0], paramtotest = k) + testsol(osys5, ODEProblem, Tsit5, u0, p, tspan; tstops = [1.0, 2.0]) + testsol(ssys5, SDEProblem, RI5, u0, p, tspan; tstops = [1.0, 2.0]) + @named osys6 = ODESystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb2‵‵‵, cb1]) @named ssys6 = SDESystem(eqs, [0.0], t, [A], [k, t1, t2], discrete_events = [cb2‵‵‵, cb1]) - testsol(ssys6, u0, p, tspan; tstops = [1.0, 2.0], paramtotest = k) + testsol(osys6, ODEProblem, Tsit5, u0, p, tspan; tstops = [1.0, 2.0], paramtotest = k) + testsol(ssys6, SDEProblem, RI5, u0, p, tspan; tstops = [1.0, 2.0], paramtotest = k) # mix a continuous event too cond3 = A ~ 0.1 affect3 = [k ~ 0.0] - cb3 = cond3 => affect3 + cb3 = SymbolicContinuousCallback(cond3 => affect3, discrete_parameters = [k], iv = t) + @named osys7 = ODESystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb1, cb2‵‵‵], + continuous_events = [cb3]) @named ssys7 = SDESystem(eqs, [0.0], t, [A], [k, t1, t2], discrete_events = [cb1, cb2‵‵‵], continuous_events = [cb3]) - sol = testsol(ssys7, u0, p, (0.0, 10.0); tstops = [1.0, 2.0]) + + sol = testsol(osys7, ODEProblem, Tsit5, u0, p, (0.0, 10.0); tstops = [1.0, 2.0]) + @test isapprox(sol(10.0)[1], 0.1; atol = 1e-10, rtol = 1e-10) + sol = testsol(ssys7, SDEProblem, RI5, u0, p, (0.0, 10.0); tstops = [1.0, 2.0]) @test isapprox(sol(10.0)[1], 0.1; atol = 1e-10, rtol = 1e-10) end -let rng = rng +@testset "JumpSystem Discrete Callbacks" begin function testsol(jsys, u0, p, tspan; tstops = Float64[], paramtotest = nothing, N = 40000, kwargs...) jsys = complete(jsys) @@ -751,22 +559,23 @@ let rng = rng jprob = JumpProblem(jsys, dprob, Direct(); kwargs...) sol = solve(jprob, SSAStepper(); tstops = tstops) @test (sol(1.000000000001)[1] - sol(0.99999999999)[1]) == 1 - paramtotest === nothing || (@test sol.ps[paramtotest] == 1.0) + paramtotest === nothing || (@test sol.ps[paramtotest] == [0.0, 1.0]) @test sol(40.0)[1] == 0 sol end - @parameters k t1 t2 + @parameters k(t) t1 t2 @variables A(t) B(t) + eqs = [MassActionJump(k, [A => 1], [A => -1])] cond1 = (t == t1) - affect1 = [A ~ A + 1] + affect1 = [A ~ Pre(A) + 1] cb1 = cond1 => affect1 cond2 = (t == t2) affect2 = [k ~ 1.0] cb2 = cond2 => affect2 + cb2 = SymbolicDiscreteCallback(cb2, discrete_parameters = [k], iv = t) - eqs = [MassActionJump(k, [A => 1], [A => -1])] @named jsys = JumpSystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb1, cb2]) u0 = [A => 1] p = [k => 0.0, t1 => 1.0, t2 => 2.0] @@ -774,7 +583,7 @@ let rng = rng testsol(jsys, u0, p, tspan; tstops = [1.0, 2.0], rng, paramtotest = k) cond1a = (t == t1) - affect1a = [A ~ A + 1, B ~ A] + affect1a = [A ~ Pre(A) + 1, B ~ A] cb1a = cond1a => affect1a @named jsys1 = JumpSystem(eqs, t, [A, B], [k, t1, t2], discrete_events = [cb1a, cb2]) u0′ = [A => 1, B => 0] @@ -784,7 +593,7 @@ let rng = rng # same as above - but with set-time event syntax cb1‵ = [1.0] => affect1 # needs to be a Vector for the event to happen only once - cb2‵ = [2.0] => affect2 + cb2‵ = SymbolicDiscreteCallback([2.0] => affect2, discrete_parameters = [k], iv = t) @named jsys‵ = JumpSystem(eqs, t, [A], [k], discrete_events = [cb1‵, cb2‵]) testsol(jsys‵, u0, [p[1]], tspan; rng, paramtotest = k) @@ -810,7 +619,7 @@ let rng = rng testsol(jsys6, u0, p, tspan; tstops = [1.0, 2.0], rng, paramtotest = k) end -let +@testset "Namespacing" begin function oscillator_ce(k = 1.0; name) sts = @variables x(t)=1.0 v(t)=0.0 F(t) ps = @parameters k=k Θ=0.5 @@ -829,8 +638,8 @@ let sol = solve(prob, Tsit5(), saveat = 0.1) @test typeof(oneosc_ce_simpl) == ODESystem - @test sol[oscce.x, 6] < 1.0 # test whether x(t) decreases over time - @test sol[oscce.x, 18] > 0.5 # test whether event happened + @test sol(0.5, idxs = oscce.x) < 1.0 # test whether x(t) decreases over time + @test sol(1.5, idxs = oscce.x) > 0.5 # test whether event happened end @testset "Additional SymbolicContinuousCallback options" begin @@ -1066,12 +875,12 @@ end @testset "Discrete variable timeseries" begin @variables x(t) @parameters a(t) b(t) c(t) - cb1 = [x ~ 1.0] => [a ~ -a] + cb1 = SymbolicContinuousCallback([x ~ 1.0] => [a ~ -Pre(a)], discrete_parameters = [a]) function save_affect!(integ, u, p, ctx) integ.ps[p.b] = 5.0 end cb2 = [x ~ 0.5] => (save_affect!, [], [b], [b], nothing) - cb3 = 1.0 => [c ~ t] + cb3 = SymbolicDiscreteCallback(1.0 => [c ~ t], discrete_parameters = [c]) @mtkbuild sys = ODESystem(D(x) ~ cos(t), t, [x], [a, b, c]; continuous_events = [cb1, cb2], discrete_events = [cb3]) @@ -1083,6 +892,7 @@ end @test sol[b] == [2.0, 5.0, 5.0] @test sol[c] == [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0] end + @testset "Heater" begin @variables temp(t) params = @parameters furnace_on_threshold=0.5 furnace_off_threshold=0.7 furnace_power=1.0 leakage=0.1 furnace_on::Bool=false @@ -1257,7 +1067,7 @@ end f = ModelingToolkit.FunctionalAffect( f = (i, u, p, c) -> seen = true, sts = [], pars = [], discretes = []) cb1 = ModelingToolkit.SymbolicContinuousCallback( - [x ~ 0], Equation[], initialize = [x ~ 1.5], finalize = f) + [x ~ 0], nothing, initialize = [x ~ 1.5], finalize = f) @mtkbuild sys = ODESystem(D(x) ~ -1, t, [x], []; continuous_events = [cb1]) prob = ODEProblem(sys, [x => 1.0], (0.0, 2), []) sol = solve(prob, Tsit5(); dtmax = 0.01) @@ -1270,7 +1080,7 @@ end f = ModelingToolkit.FunctionalAffect( f = (i, u, p, c) -> seen = true, sts = [], pars = [], discretes = []) cb1 = ModelingToolkit.SymbolicContinuousCallback( - [x ~ 0], Equation[], initialize = [x ~ 1.5], finalize = f) + [x ~ 0], nothing, initialize = [x ~ 1.5], finalize = f) inited = false finaled = false a = ModelingToolkit.FunctionalAffect( @@ -1278,7 +1088,7 @@ end b = ModelingToolkit.FunctionalAffect( f = (i, u, p, c) -> finaled = true, sts = [], pars = [], discretes = []) cb2 = ModelingToolkit.SymbolicContinuousCallback( - [x ~ 0.1], Equation[], initialize = a, finalize = b) + [x ~ 0.1], nothing, initialize = a, finalize = b) @mtkbuild sys = ODESystem(D(x) ~ -1, t, [x], []; continuous_events = [cb1, cb2]) prob = ODEProblem(sys, [x => 1.0], (0.0, 2), []) sol = solve(prob, Tsit5()) @@ -1344,7 +1154,7 @@ end cb = [x ~ 0.0] => [x ~ 0, y ~ 1] @mtkbuild pend = ODESystem(eqs, t; continuous_events = [cb]) prob = ODEProblem(pend, [x => 1], (0.0, 3.0), guesses = [y => x]) - @test_throws "DAE initialization failed" solve(prob, Rodas5()) + @test_broken !SciMLBase.successful_retcode(solve(prob, Rodas5())) cb = [x ~ 0.0] => [y ~ 1] @mtkbuild pend = ODESystem(eqs, t; continuous_events = [cb]) @@ -1361,7 +1171,7 @@ end @mtkmodel DECAY begin @parameters begin unrelated[1:2] = zeros(2) - k = 0.0 + k(t) = 0.0 end @variables begin x(t) = 10.0 @@ -1370,7 +1180,7 @@ end D(x) ~ -k * x end @discrete_events begin - (t == 1.0) => [k ~ 1.0] + (t == 1.0) => [k ~ 1.0], [discrete_parameters = k] end end @mtkbuild decay = DECAY() @@ -1378,7 +1188,7 @@ end @test_nowarn solve(prob, Tsit5(), tstops = [1.0]) end -@testset "Array parameter updates in ImperativeEffect" begin +@testset "Array parameter updates in ImperativeAffect" begin function weird1(max_time; name) params = @parameters begin θ(t) = 0.0 @@ -1434,3 +1244,64 @@ end sol2 = solve(ODEProblem(sys2, [], (0.0, 1.0)), Tsit5()) @test 100.0 ∈ sol2[sys2.wd2.θ] end + +@testset "Implicit affects with Pre" begin + using ModelingToolkit: UnsolvableCallbackError + @parameters g + @variables x(t) y(t) λ(t) + eqs = [D(D(x)) ~ λ * x + D(D(y)) ~ λ * y - g + x^2 + y^2 ~ 1] + c_evt = [t ~ 5.0] => [x ~ Pre(x) + 0.1] + @mtkbuild pend = ODESystem(eqs, t, continuous_events = c_evt) + prob = ODEProblem(pend, [x => -1, y => 0], (0.0, 10.0), [g => 1], guesses = [λ => 1]) + sol = solve(prob, FBDF()) + @test ≈(sol(5.000001, idxs = x) - sol(4.999999, idxs = x), 0.1, rtol = 1e-4) + @test ≈(sol(5.000001, idxs = x)^2 + sol(5.000001, idxs = y)^2, 1, rtol = 1e-4) + + # Implicit affect with Pre + c_evt = [t ~ 5.0] => [x ~ Pre(x) + y^2] + @mtkbuild pend = ODESystem(eqs, t, continuous_events = c_evt) + prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 10.0), [g => 1], guesses = [λ => 1]) + sol = solve(prob, FBDF()) + @test ≈(sol(5.000001, idxs = y)^2 + sol(4.999999, idxs = x), + sol(5.000001, idxs = x), rtol = 1e-4) + @test ≈(sol(5.000001, idxs = x)^2 + sol(5.000001, idxs = y)^2, 1, rtol = 1e-4) + + # Impossible affect errors + c_evt = [t ~ 5.0] => [x ~ Pre(x) + 2] + @mtkbuild pend = ODESystem(eqs, t, continuous_events = c_evt) + prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 10.0), [g => 1], guesses = [λ => 1]) + @test_throws UnsolvableCallbackError sol=solve(prob, FBDF()) + + # Changing both variables and parameters in the same affect. + @parameters g(t) + eqs = [D(D(x)) ~ λ * x + D(D(y)) ~ λ * y - g + x^2 + y^2 ~ 1] + c_evt = SymbolicContinuousCallback( + [t ~ 5.0], [x ~ Pre(x) + 0.1, g ~ Pre(g) + 1], discrete_parameters = [g], iv = t) + @mtkbuild pend = ODESystem(eqs, t, continuous_events = c_evt) + prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 10.0), [g => 1], guesses = [λ => 1]) + sol = solve(prob, FBDF()) + @test sol.ps[g] ≈ [1, 2] + @test ≈(sol(5.0000001, idxs = x) - sol(4.999999, idxs = x), 0.1, rtol = 1e-4) + + # Proper re-initialization after parameter change + eqs = [y ~ g^2, D(x) ~ x] + c_evt = SymbolicContinuousCallback( + [t ~ 5.0], [x ~ Pre(x) + 1, g ~ Pre(g) + 1], discrete_parameters = [g], iv = t) + @mtkbuild sys = ODESystem(eqs, t, continuous_events = c_evt) + prob = ODEProblem(sys, [x => 1.0], (0.0, 10.0), [g => 2]) + sol = solve(prob, FBDF()) + @test sol.ps[g] ≈ [2.0, 3.0] + @test ≈(sol(5.00000001, idxs = x) - sol(4.9999999, idxs = x), 1; rtol = 1e-4) + @test ≈(sol(5.00000001, idxs = y), 9, rtol = 1e-4) + + # Parameters that don't appear in affects should not be mutated. + c_evt = [t ~ 5.0] => [x ~ Pre(x) + 1] + @mtkbuild sys = ODESystem(eqs, t, continuous_events = c_evt) + prob = ODEProblem(sys, [x => 0.5], (0.0, 10.0), [g => 2], guesses = [y => 0]) + sol = solve(prob, FBDF()) + @test prob.ps[g] == sol.ps[g] +end diff --git a/test/symbolic_indexing_interface.jl b/test/symbolic_indexing_interface.jl index 8b3da5fd72..613bfc8213 100644 --- a/test/symbolic_indexing_interface.jl +++ b/test/symbolic_indexing_interface.jl @@ -1,5 +1,6 @@ using ModelingToolkit, SymbolicIndexingInterface, SciMLBase -using ModelingToolkit: t_nounits as t, D_nounits as D, ParameterIndex +using ModelingToolkit: t_nounits as t, D_nounits as D, ParameterIndex, + SymbolicContinuousCallback using SciMLStructures: Tunable @testset "ODESystem" begin @@ -230,7 +231,8 @@ end @testset "`timeseries_parameter_index` on unwrapped scalarized timeseries parameter" begin @variables x(t)[1:2] @parameters p(t)[1:2, 1:2] - ev = [x[1] ~ 2.0] => [p ~ -ones(2, 2)] + ev = SymbolicContinuousCallback( + [x[1] ~ 2.0] => [p ~ -ones(2, 2)], discrete_parameters = [p]) @mtkbuild sys = ODESystem(D(x) ~ p * x, t; continuous_events = [ev]) p = ModelingToolkit.unwrap(p) @test timeseries_parameter_index(sys, p) === ParameterTimeseriesIndex(1, (1, 1)) diff --git a/test/symbolic_parameters.jl b/test/symbolic_parameters.jl index a29090912c..f4fa21e614 100644 --- a/test/symbolic_parameters.jl +++ b/test/symbolic_parameters.jl @@ -28,7 +28,7 @@ resolved = ModelingToolkit.varmap_to_vars(Dict(), parameters(ns), prob = NonlinearProblem(complete(ns), [u => 1.0], Pair[]) @test prob.u0 == [1.0, 1.1, 0.9] -@show sol = solve(prob, NewtonRaphson()) +sol = solve(prob, NewtonRaphson()) @variables a @parameters b @@ -43,12 +43,12 @@ res = ModelingToolkit.varmap_to_vars(Dict(), parameters(top), top = complete(top) prob = NonlinearProblem(top, [unknowns(ns, u) => 1.0, a => 1.0], []) @test prob.u0 == [1.0, 0.5, 1.1, 0.9] -@show sol = solve(prob, NewtonRaphson()) +sol = solve(prob, NewtonRaphson()) # test NullParameters+defaults prob = NonlinearProblem(top, [unknowns(ns, u) => 1.0, a => 1.0]) @test prob.u0 == [1.0, 0.5, 1.1, 0.9] -@show sol = solve(prob, NewtonRaphson()) +sol = solve(prob, NewtonRaphson()) # test initial conditions and parameters at the problem level pars = @parameters(begin diff --git a/test/variable_scope.jl b/test/variable_scope.jl index bd1d3cb0cf..59647bf441 100644 --- a/test/variable_scope.jl +++ b/test/variable_scope.jl @@ -51,13 +51,11 @@ end @test renamed([:foo :bar :baz], c) == Symbol("foo₊c") @test renamed([:foo :bar :baz], d) == :d -@parameters a b c d e f +@parameters a b c d p = [a ParentScope(b) ParentScope(ParentScope(c)) - DelayParentScope(d) - DelayParentScope(e, 2) - GlobalScope(f)] + GlobalScope(d)] level0 = ODESystem(Equation[], t, [], p; name = :level0) level1 = ODESystem(Equation[], t, [], []; name = :level1) ∘ level0 @@ -69,9 +67,7 @@ ps = ModelingToolkit.getname.(parameters(level3)) @test isequal(ps[1], :level2₊level1₊level0₊a) @test isequal(ps[2], :level2₊level1₊b) @test isequal(ps[3], :level2₊c) -@test isequal(ps[4], :level2₊level0₊d) -@test isequal(ps[5], :level1₊level0₊e) -@test isequal(ps[6], :f) +@test isequal(ps[4], :d) # Issue@2252 # Tests from PR#2354 @@ -102,40 +98,36 @@ defs = ModelingToolkit.defaults(bar) @test defs[bar.p] == 2 @test isequal(defs[bar.foo.p], bar.p) -# Issue#3101 -@variables x1(t) x2(t) x3(t) x4(t) x5(t) -x2 = ParentScope(x2) -x3 = ParentScope(ParentScope(x3)) -x4 = DelayParentScope(x4) -x5 = GlobalScope(x5) -@parameters p1 p2 p3 p4 p5 -p2 = ParentScope(p2) -p3 = ParentScope(ParentScope(p3)) -p4 = DelayParentScope(p4) -p5 = GlobalScope(p5) - -@named sys1 = ODESystem([D(x1) ~ p1, D(x2) ~ p2, D(x3) ~ p3, D(x4) ~ p4, D(x5) ~ p5], t) -@test isequal(x1, only(unknowns(sys1))) -@test isequal(p1, only(parameters(sys1))) -@named sys2 = ODESystem(Equation[], t; systems = [sys1]) -@test length(unknowns(sys2)) == 2 -@test any(isequal(x2), unknowns(sys2)) -@test length(parameters(sys2)) == 2 -@test any(isequal(p2), parameters(sys2)) -@named sys3 = ODESystem(Equation[], t) -sys3 = sys3 ∘ sys2 -@test length(unknowns(sys3)) == 4 -@test any(isequal(x3), unknowns(sys3)) -@test any(isequal(ModelingToolkit.renamespace(sys1, x4)), unknowns(sys3)) -@test length(parameters(sys3)) == 4 -@test any(isequal(p3), parameters(sys3)) -@test any(isequal(ModelingToolkit.renamespace(sys1, p4)), parameters(sys3)) -sys4 = complete(sys3) -@test length(unknowns(sys3)) == 4 -@test length(parameters(sys4)) == 5 -@test any(isequal(p5), parameters(sys4)) -sys5 = structural_simplify(sys3) -@test length(unknowns(sys5)) == 5 -@test any(isequal(x5), unknowns(sys5)) -@test length(parameters(sys5)) == 5 -@test any(isequal(p5), parameters(sys5)) +@testset "Issue#3101" begin + @variables x1(t) x2(t) x3(t) x4(t) + x2 = ParentScope(x2) + x3 = ParentScope(ParentScope(x3)) + x4 = GlobalScope(x4) + @parameters p1 p2 p3 p4 + p2 = ParentScope(p2) + p3 = ParentScope(ParentScope(p3)) + p4 = GlobalScope(p4) + + @named sys1 = ODESystem([D(x1) ~ p1, D(x2) ~ p2, D(x3) ~ p3, D(x4) ~ p4], t) + @test isequal(x1, only(unknowns(sys1))) + @test isequal(p1, only(parameters(sys1))) + @named sys2 = ODESystem(Equation[], t; systems = [sys1]) + @test length(unknowns(sys2)) == 2 + @test any(isequal(x2), unknowns(sys2)) + @test length(parameters(sys2)) == 2 + @test any(isequal(p2), parameters(sys2)) + @named sys3 = ODESystem(Equation[], t) + sys3 = sys3 ∘ sys2 + @test length(unknowns(sys3)) == 3 + @test any(isequal(x3), unknowns(sys3)) + @test length(parameters(sys3)) == 3 + @test any(isequal(p3), parameters(sys3)) + sys4 = complete(sys3) + @test length(unknowns(sys4)) == 3 + @test length(parameters(sys4)) == 4 + sys5 = structural_simplify(sys3) + @test length(unknowns(sys5)) == 4 + @test any(isequal(x4), unknowns(sys5)) + @test length(parameters(sys5)) == 4 + @test any(isequal(p4), parameters(sys5)) +end