[v10] refactor: use ImplicitDiscreteSystem to implement affects in callback systems

.github/workflows/Documentation.yml

@@ -4,6 +4,7 @@ on:
   push:
     branches:
       - master
+      - v10
     tags: '*'
   pull_request:
 

.github/workflows/FormatCheck.yml

@@ -4,6 +4,7 @@ on:
   push:
     branches:
       - 'master'
+      - v10
     tags: '*'
   pull_request:
 

.github/workflows/Tests.yml

@@ -5,6 +5,7 @@ on:
     branches:
       - master
       - 'release-'
+      - v10
     paths-ignore:
       - 'docs/**'
   push:

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

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Yingbo Ma <[email protected]>", "Chris Rackauckas <[email protected]> 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"]

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"

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"

docs/src/basics/Composition.md

@@ -135,42 +135,34 @@ 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
 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

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,16 +361,17 @@ 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`.
 
 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())