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 30580d961b..3738ad4f65 100644 --- a/Project.toml +++ b/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" diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index 23e1e6d7d1..4901910bfa 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 9f69458528..690f2cbdb7 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 @@ -157,7 +158,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") @@ -167,19 +167,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") @@ -302,6 +303,7 @@ export initialization_equations, guesses, defaults, parameter_dependencies, hier 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/parameters.jl b/src/parameters.jl index 91121b7cbb..8c0f9f1b00 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}) = Num(tovar(value(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/callbacks.jl b/src/systems/callbacks.jl index 07809bf611..6805a3f1a9 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,114 @@ 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 + system::ImplicitDiscreteSystem + unknowns::Vector + parameters::Vector + 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,467 +198,393 @@ 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) + any(a -> (a isa FunctionalAffect || a isa ImperativeAffect), + [affect, affect_neg, initialize, finalize]) ? + reinitializealg = SciMLBase.CheckInit() : + reinitializealg = SciMLBase.NoInit() + 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 +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 + +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 + isempty(alg_eqs) && warn_no_algebraic && + @warn "No algebraic equations were found for the callback defined by $(join(affect, ", ")). If the system has no algebraic equations, this can be disregarded. Otherwise pass in `alg_eqs` to the SymbolicContinuousCallback constructor." + 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)." + 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 = false) + # 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 -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) + +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::SymbolicContinuousCallback) +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 appaers 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 + @show cbs + [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 +Base.isempty(cb::AbstractCallback) = isempty(cb.conditions) + +#################################### +####### Compilation functions ###### +#################################### """ - compile_condition(cb::SymbolicDiscreteCallback, sys, dvs, ps; expression, kwargs...) + compile_condition(cb::AbstractCallback, sys, dvs, ps; expression, kwargs...) -Returns a function `condition(u,t,integrator)` returning the `condition(cb)`. +Returns a function `condition(u,t,integrator)`, condition(out,u,t,integrator)` returning the `condition(cb)`. Notes @@ -596,527 +592,379 @@ Notes If set to `Val{false}` a `RuntimeGeneratedFunction` will be returned. - `kwargs` are passed through to `Symbolics.build_function`. """ -function compile_condition(cb::SymbolicDiscreteCallback, sys, dvs, ps; - expression = Val{true}, eval_expression = false, eval_module = @__MODULE__, kwargs...) +function compile_condition( + cbs::Union{AbstractCallback, Vector{<:AbstractCallback}}, sys, dvs, ps; + expression = Val{false}, 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)) 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; expression, + kwargs...) -Notes + if expression == Val{true} + fs = eval_or_rgf.(fs; eval_expression, eval_module) + end + f_oop, f_iip = is_discrete(cbs) ? (fs, nothing) : fs # no iip function for discrete condition. - - `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(integ.sol.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 = function (integ, idx) + affects[eq2affect[idx]](integ) end - if length(compiled_callbacks) == 1 - return compiled_callbacks[] - else - return CallbackSet(compiled_callbacks...) + affect_neg = function (integ, idx) + f = affect_negs[eq2affect[idx]] + isnothing(f) && return + f(integ) end + initialize = wrap_vector_optional_affect(inits, SciMLBase.INITIALIZE_DEFAULT) + finalize = wrap_vector_optional_affect(finals, SciMLBase.FINALIZE_DEFAULT) + + return VectorContinuousCallback( + trigger, affect, affect_neg, length(eqs); initialize, finalize, + rootfind = cbs[1].rootfind, initializealg = cbs[1].reinitializealg) end -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)) +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) && get_index_cache(sys) !== nothing - p_inds = [if (pind = parameter_index(sys, sym)) === nothing - sym - else - pind - end - for sym in parameters(affect)] + 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 - ps_ind = Dict(reverse(en) for en in enumerate(ps)) - p_inds = map(sym -> get(ps_ind, sym, sym), parameters(affect)) + (cb.affect === cb.affect_neg) ? affect : + compile_affect(cb.affect_neg, cb, sys; default = EMPTY_AFFECT, kwargs...) 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 - - 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 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 + + - `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( + 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 + if f === SciMLBase.INITIALIZE_DEFAULT + (c, u, t, i) -> begin + isnothing(f) || 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 + +""" +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 + end 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...) +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 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 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 + 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 + + function implicit_affect!(integ) + pmap = Pair[] + for pre_p in parameters(affsys) + p = unPre(pre_p) + pval = isparameter(p) ? integ.ps[p] : integ[p] + push!(pmap, pre_p => pval) + end + u0 = Pair[] + for u in unknowns(affsys) + uval = isparameter(aff_map[u]) ? integ.ps[aff_map[u]] : integ[u] + push!(u0, u => uval) + end + affprob = ImplicitDiscreteProblem(affsys, u0, (integ.t, integ.t), pmap; + build_initializeprob = false, check_length = false) + @show pmap + @show u0 + affsol = init(affprob, IDSolve()) + @show affsol + (check_error(affsol) === ReturnCode.InitialFailure) && + throw(UnsolvableCallbackError(all_equations(aff))) + for u in dvs_to_update + integ[u] = affsol[sys_map[u]] + end + for p in ps_to_update + integ.ps[p] = affsol[sys_map[p]] 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) - 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 - - dbs = map(symcbs) do cb - generate_discrete_callback(cb, sys, dvs, ps; kwargs...) - end +struct UnsolvableCallbackError + eqs::Vector{Equation} +end - dbs +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 +972,96 @@ 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 - end - if has_discrete_events(sys) && !isempty(discrete_events(sys)) - discrete_cb = generate_discrete_callbacks(sys; kwargs...) - else - discrete_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 + return get_discrete_events(sys) +end + +""" + continuous_events(sys::AbstractSystem)::Vector{SymbolicContinuousCallback} - cb = merge_cb(contin_cb, callback) - (discrete_cb === nothing) ? cb : CallbackSet(contin_cb, discrete_cb...) +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 diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index cbce569a9b..9b72549e65 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -317,8 +317,13 @@ 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) + + alg_eqs = filter(eq -> eq.lhs isa Union{Symbolic, Number} && !is_diff_equation(eq), + deqs) + cont_callbacks = to_cb_vector(continuous_events; CB_TYPE = SymbolicContinuousCallback, + iv = iv, alg_eqs = alg_eqs, warn_no_algebraic = false) + disc_callbacks = to_cb_vector(discrete_events; CB_TYPE = SymbolicDiscreteCallback, + iv = iv, alg_eqs = alg_eqs, warn_no_algebraic = false) if is_dde === nothing is_dde = _check_if_dde(deqs, iv′, systems) diff --git a/src/systems/diffeqs/sdesystem.jl b/src/systems/diffeqs/sdesystem.jl index 3fa1302630..5b636c39e1 100644 --- a/src/systems/diffeqs/sdesystem.jl +++ b/src/systems/diffeqs/sdesystem.jl @@ -269,8 +269,14 @@ 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) + + alg_eqs = filter(eq -> eq.lhs isa Union{Symbolic, Number} && !is_diff_equation(eq), + deqs) + cont_callbacks = to_cb_vector(continuous_events; CB_TYPE = SymbolicContinuousCallback, + iv = iv, alg_eqs = alg_eqs, warn_no_algebraic = false) + disc_callbacks = to_cb_vector(discrete_events; CB_TYPE = SymbolicDiscreteCallback, + iv = iv, alg_eqs = alg_eqs, warn_no_algebraic = false) + 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 5f7c986659..1d0832bb36 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..eb3a094ef5 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) @@ -297,11 +299,11 @@ function shift_u0map_forward(sys::ImplicitDiscreteSystem, u0map, defs) 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)).") + @warn "Initial condition given in term of current state of the unknown. If `build_initializeprob = false`, this may be overriden by the implicit discrete solver." 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 +382,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 +402,7 @@ function SciMLBase.ImplicitDiscreteFunction{iip, specialize}( sys = sys, observed = observedfun, analytic = analytic, + resid_prototype = resid_prototype, kwargs...) end @@ -441,3 +450,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 47a784c00b..d71bcc60a9 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 57a3aee7df..7ecd6e6202 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) @@ -89,10 +75,6 @@ struct JumpSystem{U <: ArrayPartition} <: AbstractTimeDependentSystem Type of the system. """ connector_type::Any - """ - A `Vector{SymbolicContinuousCallback}` that model events. - The integrator will use root finding to guarantee that it steps at each zero crossing. - """ continuous_events::Vector{SymbolicContinuousCallback} """ A `Vector{SymbolicDiscreteCallback}` that models events. Symbolic @@ -230,8 +212,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 +266,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 +281,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 +289,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 +303,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 +311,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 +552,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 195b02118e..cbd0787e69 100644 --- a/src/systems/model_parsing.jl +++ b/src/systems/model_parsing.jl @@ -126,7 +126,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)) + name, description = $description, systems, gui_metadata = $gui_metadata, + continuous_events = [$(c_evts...)], discrete_events = [$(d_evts...)], + defaults)) if length(ext) == 0 push!(exprs.args, :(var"#___sys___" = $sys)) @@ -137,16 +139,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 @@ -1134,8 +1126,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 @@ -1143,12 +1143,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_icon!(body::String, dict, icon, mod) icon_dir = get(ENV, "MTK_ICONS_DIR", joinpath(DEPOT_PATH[1], "mtk_icons")) dict[:icon] = icon[] = if isfile(body) diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index 0750585905..8b96463e3c 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -380,9 +380,7 @@ function better_varmap_to_vars(varmap::AbstractDict, vars::Vector; vals = promote_to_concrete(vals; tofloat = tofloat, use_union = false) end - if isempty(vals) - return nothing - elseif container_type <: Tuple + if container_type <: Tuple return (vals...,) else return SymbolicUtils.Code.create_array(container_type, eltype(vals), Val{1}(), diff --git a/src/systems/systems.jl b/src/systems/systems.jl index 0f8633f31f..4d7216d081 100644 --- a/src/systems/systems.jl +++ b/src/systems/systems.jl @@ -41,10 +41,10 @@ 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. - """) + #error(""" + # Encountered algebraic equations when simplifying discrete system. Please construct \ + # an ImplicitDiscreteSystem instead. + #""") end for pass in additional_passes newsys = pass(newsys) @@ -155,9 +155,9 @@ function __structural_simplify(sys::AbstractSystem, io = nothing; simplify = fal 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/utils.jl b/src/utils.jl index 1884a91c19..6ca95341b1 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) 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/discrete_system.jl b/test/discrete_system.jl index 0d215052d8..7fa8349c8e 100644 --- a/test/discrete_system.jl +++ b/test/discrete_system.jl @@ -252,7 +252,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/extensions/ad.jl b/test/extensions/ad.jl index adaf6117c6..845d1ad818 100644 --- a/test/extensions/ad.jl +++ b/test/extensions/ad.jl @@ -59,7 +59,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/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 3004044d61..6e699d1838 100644 --- a/test/funcaffect.jl +++ b/test/funcaffect.jl @@ -24,8 +24,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 = [], @@ -46,7 +45,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 4ade3481cb..aac0566ef5 100644 --- a/test/initializationsystem.jl +++ b/test/initializationsystem.jl @@ -1285,9 +1285,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/jumpsystem.jl b/test/jumpsystem.jl index 9568990e73..82299adb8c 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() @@ -389,9 +390,9 @@ let p4 = DelayParentScope(p4) p5 = GlobalScope(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]) @named js = JumpSystem([j1, j2, j3, j4], t, [x1, x2, x3, x4, x5], [p1, p2, p3, p4, p5]) @@ -429,8 +430,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 +472,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 +540,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) diff --git a/test/model_parsing.jl b/test/model_parsing.jl index e8464707de..0cdf294a1f 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 22201b1988..9fefe40af2 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 4b76da6e9d..20d70bbfda 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 @@ -1028,24 +1029,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 @@ -1192,10 +1196,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 @@ -1533,8 +1539,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 @@ -1594,15 +1604,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)) 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/symbolic_events.jl b/test/symbolic_events.jl index 804432408b..adcb24dc76 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 @@ -980,7 +789,7 @@ end @test sign.(cos.(required_crossings_c1 .- 1e-6)) == sign.(last.(cr1)) @test sign.(cos.(3 * (required_crossings_c2 .+ 1e-6))) == sign.(last.(cr2)) end - +# @testset "Discrete event reinitialization (#3142)" begin @connector LiquidPort begin p(t)::Float64, [description = "Set pressure in bar", @@ -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