feat: implement new callback semantics

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

@@ -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"
@@ -139,7 +142,7 @@ RecursiveArrayTools = "3.26"
 Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.5.9"
 SCCNonlinearSolve = "1.0.0"
-SciMLBase = "2.84"
+SciMLBase = "2.89.1"
 SciMLStructures = "1.7"
 Serialization = "1"
 Setfield = "0.7, 0.8, 1"

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())

docs/src/basics/MTKLanguage.md

@@ -203,14 +203,15 @@ 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
 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
 ```

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
 

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
 

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)

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 = SciMLBase.NoInit())
 
         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)