From a6f781d45cbbc885edf9167695682124cd9655bb Mon Sep 17 00:00:00 2001 From: Ben Chung Date: Thu, 1 Aug 2024 20:36:50 -0700 Subject: [PATCH 01/47] First pass at MutatingFunctionalAffect --- Project.toml | 1 + src/ModelingToolkit.jl | 1 + src/systems/callbacks.jl | 144 +++++++++++++++++++++++++++++-- src/systems/diffeqs/odesystem.jl | 27 +++++- test/symbolic_events.jl | 49 +++++++++++ 5 files changed, 215 insertions(+), 7 deletions(-) diff --git a/Project.toml b/Project.toml index 808d68af06..729966fde0 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" +ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index 2f57bb1765..f5262a1526 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -54,6 +54,7 @@ using Reexport using RecursiveArrayTools import Graphs: SimpleDiGraph, add_edge!, incidence_matrix import BlockArrays: BlockedArray, Block, blocksize, blocksizes +import ComponentArrays using RuntimeGeneratedFunctions using RuntimeGeneratedFunctions: drop_expr diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index 86cab57634..f7d4baa4cb 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -60,8 +60,6 @@ function Base.hash(a::FunctionalAffect, s::UInt) hash(a.ctx, s) end -has_functional_affect(cb) = affects(cb) isa FunctionalAffect - namespace_affect(affect, s) = namespace_equation(affect, s) function namespace_affect(affect::FunctionalAffect, s) FunctionalAffect(func(affect), @@ -73,6 +71,67 @@ function namespace_affect(affect::FunctionalAffect, s) context(affect)) end +""" +`MutatingFunctionalAffect` differs from `FunctionalAffect` in two key ways: +* First, insetad of the `u` vector passed to `f` being a vector of indices into `integ.u` it's instead the result of evaluating `obs` at the current state, named as specified in `obs_syms`. This allows affects to easily access observed states and decouples affect inputs from the system structure. +* Second, it abstracts the assignment back to system states away. Instead of writing `integ.u[u.myvar] = [whatever]`, you instead declare in `mod_params` that you want to modify `myvar` and then either (out of place) return a named tuple with `myvar` or (in place) modify the associated element in the ComponentArray that's given. +Initially, we only support "flat" states in `modified`; these states will be marked as irreducible in the overarching system and they will simply be bulk assigned at mutation. In the future, this will be extended to perform a nonlinear solve to further decouple the affect from the system structure. +""" +@kwdef struct MutatingFunctionalAffect + f::Any + obs::Vector + obs_syms::Vector{Symbol} + modified::Vector + mod_syms::Vector{Symbol} + ctx::Any +end + +MutatingFunctionalAffect(f::Function; + observed::NamedTuple = NamedTuple{()}(()), + modified::NamedTuple = NamedTuple{()}(()), + ctx=nothing) = MutatingFunctionalAffect(f, collect(values(observed)), collect(keys(observed)), collect(values(modified)), collect(keys(modified)), ctx) +MutatingFunctionalAffect(f::Function, observed::NamedTuple; modified::NamedTuple = NamedTuple{()}(()), ctx=nothing) = + MutatingFunctionalAffect(f, observed=observed, modified=modified, ctx=ctx) +MutatingFunctionalAffect(f::Function, observed::NamedTuple, modified::NamedTuple; ctx=nothing) = + MutatingFunctionalAffect(f, observed=observed, modified=modified, ctx=ctx) +MutatingFunctionalAffect(f::Function, observed::NamedTuple, modified::NamedTuple, ctx) = + MutatingFunctionalAffect(f, observed=observed, modified=modified, ctx=ctx) + +func(f::MutatingFunctionalAffect) = f.f +context(a::MutatingFunctionalAffect) = a.ctx +observed(a::MutatingFunctionalAffect) = a.obs +observed_syms(a::MutatingFunctionalAffect) = a.obs_syms +discretes(a::MutatingFunctionalAffect) = filter(ModelingToolkit.isparameter, a.modified) +modified(a::MutatingFunctionalAffect) = a.modified +modified_syms(a::MutatingFunctionalAffect) = a.mod_syms + +function Base.:(==)(a1::MutatingFunctionalAffect, a2::MutatingFunctionalAffect) + isequal(a1.f, a2.f) && isequal(a1.obs, a2.obs) && isequal(a1.modified, a2.modified) && + isequal(a1.obs_syms, a2.obs_syms) && isequal(a1.mod_syms, a2.mod_syms)&& isequal(a1.ctx, a2.ctx) +end + +function Base.hash(a::MutatingFunctionalAffect, s::UInt) + s = hash(a.f, s) + s = hash(a.obs, s) + s = hash(a.obs_syms, s) + s = hash(a.modified, s) + s = hash(a.mod_syms, s) + hash(a.ctx, s) +end + +function namespace_affect(affect::MutatingFunctionalAffect, s) + MutatingFunctionalAffect(func(affect), + renamespace.((s,), observed(affect)), + observed_syms(affect), + renamespace.((s,), modified(affect)), + modified_syms(affect), + context(affect)) +end + +function has_functional_affect(cb) + (affects(cb) isa FunctionalAffect || affects(cb) isa MutatingFunctionalAffect) +end + #################################### continuous events ##################################### const NULL_AFFECT = Equation[] @@ -109,8 +168,8 @@ Affects (i.e. `affect` and `affect_neg`) can be specified as either: """ struct SymbolicContinuousCallback eqs::Vector{Equation} - affect::Union{Vector{Equation}, FunctionalAffect} - affect_neg::Union{Vector{Equation}, FunctionalAffect, Nothing} + affect::Union{Vector{Equation}, FunctionalAffect, MutatingFunctionalAffect} + affect_neg::Union{Vector{Equation}, FunctionalAffect, MutatingFunctionalAffect, Nothing} rootfind::SciMLBase.RootfindOpt function SymbolicContinuousCallback(; eqs::Vector{Equation}, affect = NULL_AFFECT, affect_neg = affect, rootfind = SciMLBase.LeftRootFind) @@ -250,6 +309,7 @@ scalarize_affects(affects) = scalarize(affects) scalarize_affects(affects::Tuple) = FunctionalAffect(affects...) scalarize_affects(affects::NamedTuple) = FunctionalAffect(; affects...) scalarize_affects(affects::FunctionalAffect) = affects +scalarize_affects(affects::MutatingFunctionalAffect) = affects SymbolicDiscreteCallback(p::Pair) = SymbolicDiscreteCallback(p[1], p[2]) SymbolicDiscreteCallback(cb::SymbolicDiscreteCallback) = cb # passthrough @@ -257,7 +317,7 @@ 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 + if db.affects isa FunctionalAffect || db.affects isa MutatingFunctionalAffect # TODO println(io, " ", db.affects) else @@ -749,6 +809,80 @@ function compile_user_affect(affect::FunctionalAffect, cb, sys, dvs, ps; kwargs. end end +invalid_variables(sys, expr) = filter(x -> !any(isequal(x), all_symbols(sys)), vars(expr)) +function unassignable_variables(sys, expr) + assignable_syms = vcat(unknowns(sys), parameters(sys)) + return filter(x -> !any(isequal(x), assignable_syms), vars(expr)) +end + +function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwargs...) + #= + Implementation sketch: + generate observed function (oop), should save to a component array under obs_syms + do the same stuff as the normal FA for pars_syms + call the affect method - test if it's OOP or IP using applicable + unpack and apply the resulting values + =# + obs_exprs = observed(affect) + for oexpr in obs_exprs + invalid_vars = invalid_variables(sys, oexpr) + if length(invalid_vars) > 0 + error("Observed equation $(oexpr) in affect refers to missing variable(s) $(invalid_vars); the variables may not have been added (e.g. if a component is missing).") + end + end + obs_syms = observed_syms(affect) + obs_size = size.(obs_exprs) # we will generate a work buffer of a ComponentArray that maps obs_syms to arrays of size obs_size + + mod_exprs = modified(affect) + for mexpr in mod_exprs + if !is_observed(sys, mexpr) && parameter_index(sys, mexpr) === nothing + error("Expression $mexpr cannot be assigned to; currently only unknowns and parameters may be updated by an affect.") + end + invalid_vars = unassignable_variables(sys, mexpr) + if length(invalid_vars) > 0 + error("Observed equation $(mexpr) in affect refers to missing variable(s) $(invalid_vars); the variables may not have been added (e.g. if a component is missing) or they may have been reduced away.") + end + end + mod_syms = modified_syms(affect) + _, mod_og_val_fun = build_explicit_observed_function(sys, mod_exprs; return_inplace=true) + + # sanity checks done! now build the data and update function for observed values + mkzero(sz) = if sz === () 0.0 else zeros(sz) end + _, obs_fun = build_explicit_observed_function(sys, reduce(vcat, Symbolics.scalarize.(obs_exprs); init = []); return_inplace=true) + obs_component_array = ComponentArrays.ComponentArray(NamedTuple{(obs_syms..., )}(mkzero.(obs_size))) + + # okay so now to generate the stuff to assign it back into the system + # note that we reorder the componentarray to make the views coherent wrt the base array + mod_pairs = mod_exprs .=> mod_syms + mod_param_pairs = filter(v -> is_parameter(sys, v[1]), mod_pairs) + mod_unk_pairs = filter(v -> !is_parameter(sys, v[1]), mod_pairs) + _, mod_og_val_fun = build_explicit_observed_function(sys, reduce(vcat, [first.(mod_param_pairs); first.(mod_unk_pairs)]; init = []); return_inplace=true) + upd_params_fun = setu(sys, reduce(vcat, Symbolics.scalarize.(first.(mod_param_pairs)); init = [])) + upd_unk_fun = setu(sys, reduce(vcat, Symbolics.scalarize.(first.(mod_unk_pairs)); init = [])) + + upd_component_array = ComponentArrays.ComponentArray(NamedTuple{([last.(mod_param_pairs); last.(mod_unk_pairs)]...,)}( + [collect(mkzero(size(e)) for e in first.(mod_param_pairs)); + collect(mkzero(size(e)) for e in first.(mod_unk_pairs))])) + upd_params_view = view(upd_component_array, last.(mod_param_pairs)) + upd_unks_view = view(upd_component_array, last.(mod_unk_pairs)) + let user_affect = func(affect), ctx = context(affect) + function (integ) + # update the to-be-mutated values; this ensures that if you do a no-op then nothing happens + mod_og_val_fun(upd_component_array, integ.u, integ.p..., integ.t) + + # update the observed values + obs_fun(obs_component_array, integ.u, integ.p..., integ.t) + + # let the user do their thing + user_affect(upd_component_array, obs_component_array, integ, ctx) + + # write the new values back to the integrator + upd_params_fun(integ, upd_params_view) + upd_unk_fun(integ, upd_unks_view) + end + end +end + function compile_affect(affect::FunctionalAffect, cb, sys, dvs, ps; kwargs...) compile_user_affect(affect, cb, sys, dvs, ps; kwargs...) end diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index ec2c5f8157..daa4321ed0 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -404,8 +404,31 @@ ODESystem(eq::Equation, args...; kwargs...) = ODESystem([eq], args...; kwargs... """ $(SIGNATURES) -Build the observed function assuming the observed equations are all explicit, -i.e. there are no cycles. +Generates a function that computes the observed value(s) `ts` in the system `sys` assuming that there are no cycles in the equations. + +The return value will be either: +* a single function if the input is a scalar or if the input is a Vector but `return_inplace` is false +* the out of place and in-place functions `(ip, oop)` if `return_inplace` is true and the input is a `Vector` + +The function(s) will be: +* `RuntimeGeneratedFunction`s by default, +* A Julia `Expr` if `expression` is true, +* A directly evaluated Julia function in the module `eval_module` if `eval_expression` is true + +The signatures will be of the form `g(...)` with arguments: +* `output` for in-place functions +* `unknowns` if `params_only` is `false` +* `inputs` if `inputs` is an array of symbolic inputs that should be available in `ts` +* `p...` unconditionally; note that in the case of `MTKParameters` more than one parameters argument may be present, so it must be splatted +* `t` if the system is time-dependent; for example `NonlinearSystem` will not have `t` +For example, a function `g(op, unknowns, p, inputs, t)` will be the in-place function generated if `return_inplace` is true, `ts` is a vector, an array of inputs `inputs` is given, and `params_only` is false for a time-dependent system. + +Options not otherwise specified are: +* `output_type = Array` the type of the array generated by the out-of-place vector-valued function +* `checkbounds = true` checks bounds if true when destructuring parameters +* `op = Operator` sets the recursion terminator for the walk done by `vars` to identify the variables that appear in `ts`. See the documentation for `vars` for more detail. +* `throw = true` if true, throw an error when generating a function for `ts` that reference variables that do not exist +* `drop_expr` is deprecated. """ function build_explicit_observed_function(sys, ts; inputs = nothing, diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index e1d12814ef..351f8111ee 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -887,3 +887,52 @@ 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 + eqs = [ + D(temp) ~ furnace_on * furnace_power - temp^2 * leakage + ] + + furnace_off = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_off_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on)) do x, o, i, c + x.furnace_on = false + end) + furnace_enable = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_on_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on)) do x, o, i, c + x.furnace_on = true + end) + + @named sys = ODESystem(eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) + ss = structural_simplify(sys) + prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + sol = solve(prob, Tsit5(); dtmax=0.01) + @test all(sol[temp][sol.t .> 1.0] .<= 0.79) && all(sol[temp][sol.t .> 1.0] .>= 0.49) +end + +@testset "Quadrature" begin + @variables theta(t) omega(t) + params = @parameters qA=0 qB=0 + eqs = [ + D(theta) ~ omega + omega ~ sin(0.5*t) + ] + qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(1000 * theta) ~ 0], + ModelingToolkit.MutatingFunctionalAffect(modified=(; qA)) do x, o, i, c + x.qA = 1 + end, + affect_neg = ModelingToolkit.MutatingFunctionalAffect(modified=(; qA)) do x, o, i, c + x.qA = 0 + end) + qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(1000 * theta + π/2) ~ 0], + ModelingToolkit.MutatingFunctionalAffect(modified=(; qB)) do x, o, i, c + x.qB = 1 + end, + affect_neg = ModelingToolkit.MutatingFunctionalAffect(modified=(; qB)) do x, o, i, c + x.qB = 0 + end) + @named sys = ODESystem(eqs, t, [theta, omega], params; continuous_events = [qAevt, qBevt]) + ss = structural_simplify(sys) + prob = ODEProblem(ss, [theta => 0.0], (0.0, 1.0)) + sol = solve(prob, Tsit5(); dtmax=0.01) +end From f151e429cd9edb33b55ec019be069af12cb9d108 Mon Sep 17 00:00:00 2001 From: Ben Chung Date: Fri, 2 Aug 2024 16:19:51 -0700 Subject: [PATCH 02/47] Clarify documentation for SCC --- src/systems/callbacks.jl | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index f7d4baa4cb..a57d5c006d 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -145,17 +145,26 @@ By default `affect_neg = affect`; to only get rising edges specify `affect_neg = Assume without loss of generality that the equation is of the form `c(u,p,t) ~ 0`; we denote the integrator state as `i.u`. For compactness, we define `prev_sign = sign(c(u[t-1], p[t-1], t-1))` and `cur_sign = sign(c(u[t], p[t], t))`. A condition edge will be detected and the callback will be invoked iff `prev_sign * cur_sign <= 0`. +The positive edge `affect` will be triggered iff an edge is detected and if `prev_sign < 0`; similarly, `affect_neg` will be +triggered iff an edge is detected and `prev_sign > 0`. + Inter-sample condition activation is not guaranteed; for example if we use the dirac delta function as `c` to insert a sharp discontinuity between integrator steps (which in this example would not normally be identified by adaptivity) then the condition is not guaranteed to be triggered. Once detected the integrator will "wind back" through a root-finding process to identify the point when the condition became active; the method used -is specified by `rootfind` from [`SciMLBase.RootfindOpt`](@ref). Multiple callbacks in the same system with different `rootfind` operations will be resolved -into separate VectorContinuousCallbacks in the enumeration order of `SciMLBase.RootfindOpt`, which may cause some callbacks to not fire if several become -active at the same instant. See the `SciMLBase` documentation for more information on the semantic rules. - -The positive edge `affect` will be triggered iff an edge is detected and if `prev_sign < 0`; similarly, `affect_neg` will be -triggered iff an edge is detected `prev_sign > 0`. +is specified by `rootfind` from [`SciMLBase.RootfindOpt`](@ref). If we denote the time when the condition becomes active at tc, +the value in the integrator after windback will be: +* `u[tc-epsilon], p[tc-epsilon], tc` if `LeftRootFind` is used, +* `u[tc+epsilon], p[tc+epsilon], tc` if `RightRootFind` is used, +* or `u[t], p[t], t` if `NoRootFind` is used. +For example, if we want to detect when an unknown variable `x` satisfies `x > 0` using the condition `x ~ 0` on a positive edge (that is, `D(x) > 0`), +then left root finding will get us `x=-epsilon`, right root finding `x=epsilon` and no root finding whatever the next step of the integrator was after +it passed through 0. + +Multiple callbacks in the same system with different `rootfind` operations will be grouped +by their `rootfind` value into separate VectorContinuousCallbacks in the enumeration order of `SciMLBase.RootfindOpt`. This may cause some callbacks to not fire if several become +active at the same instant. See the `SciMLBase` documentation for more information on the semantic rules. Affects (i.e. `affect` and `affect_neg`) can be specified as either: * A list of equations that should be applied when the callback is triggered (e.g. `x ~ 3, y ~ 7`) which must be of the form `unknown ~ observed value` where each `unknown` appears only once. Equations will be applied in the order that they appear in the vector; parameters and state updates will become immediately visible to following equations. From 49d48b833645a17f18f7581009f2a15a9c708161 Mon Sep 17 00:00:00 2001 From: Ben Chung Date: Fri, 2 Aug 2024 16:20:05 -0700 Subject: [PATCH 03/47] MutatingFunctionalAffect test cases --- test/symbolic_events.jl | 160 +++++++++++++++++++++++++++++++++++++--- 1 file changed, 149 insertions(+), 11 deletions(-) diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index 351f8111ee..ca9f0ad9c3 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -227,6 +227,117 @@ affect_neg = [x ~ 1] @test e[].affect == affect end +@testset "MutatingFunctionalAffect constructors" begin + fmfa(o, x, i, c) = nothing + m = ModelingToolkit.MutatingFunctionalAffect(fmfa) + @test m isa ModelingToolkit.MutatingFunctionalAffect + @test m.f == fmfa + @test m.obs == [] + @test m.obs_syms == [] + @test m.modified == [] + @test m.mod_syms == [] + @test m.ctx === nothing + + m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (;)) + @test m isa ModelingToolkit.MutatingFunctionalAffect + @test m.f == fmfa + @test m.obs == [] + @test m.obs_syms == [] + @test m.modified == [] + @test m.mod_syms == [] + @test m.ctx === nothing + + m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (; x)) + @test m isa ModelingToolkit.MutatingFunctionalAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:x] + @test m.modified == [] + @test m.mod_syms == [] + @test m.ctx === nothing + + m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (; y=x)) + @test m isa ModelingToolkit.MutatingFunctionalAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:y] + @test m.modified == [] + @test m.mod_syms == [] + @test m.ctx === nothing + + m = ModelingToolkit.MutatingFunctionalAffect(fmfa; observed=(; y=x)) + @test m isa ModelingToolkit.MutatingFunctionalAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:y] + @test m.modified == [] + @test m.mod_syms == [] + @test m.ctx === nothing + + m = ModelingToolkit.MutatingFunctionalAffect(fmfa; modified=(; x)) + @test m isa ModelingToolkit.MutatingFunctionalAffect + @test m.f == fmfa + @test isequal(m.obs, []) + @test m.obs_syms == [] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:x] + @test m.ctx === nothing + + m = ModelingToolkit.MutatingFunctionalAffect(fmfa; modified=(; y=x)) + @test m isa ModelingToolkit.MutatingFunctionalAffect + @test m.f == fmfa + @test isequal(m.obs, []) + @test m.obs_syms == [] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:y] + @test m.ctx === nothing + + m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (; x), (; x)) + @test m isa ModelingToolkit.MutatingFunctionalAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:x] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:x] + @test m.ctx === nothing + + m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (; y=x), (; y=x)) + @test m isa ModelingToolkit.MutatingFunctionalAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:y] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:y] + @test m.ctx === nothing + + m = ModelingToolkit.MutatingFunctionalAffect(fmfa; modified=(; y=x), observed=(; y=x)) + @test m isa ModelingToolkit.MutatingFunctionalAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:y] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:y] + @test m.ctx === nothing + + m = ModelingToolkit.MutatingFunctionalAffect(fmfa; modified=(; y=x), observed=(; y=x), ctx=3) + @test m isa ModelingToolkit.MutatingFunctionalAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:y] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:y] + @test m.ctx === 3 + + m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (; x), (; x), 3) + @test m isa ModelingToolkit.MutatingFunctionalAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:x] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:x] + @test m.ctx === 3 +end + ## @named sys = ODESystem(eqs, t, continuous_events = [x ~ 1]) @@ -912,27 +1023,54 @@ end @testset "Quadrature" begin @variables theta(t) omega(t) - params = @parameters qA=0 qB=0 + params = @parameters qA=0 qB=0 hA=0 hB=0 cnt=0 eqs = [ D(theta) ~ omega - omega ~ sin(0.5*t) + omega ~ 1.0 ] - qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(1000 * theta) ~ 0], - ModelingToolkit.MutatingFunctionalAffect(modified=(; qA)) do x, o, i, c + function decoder(oldA, oldB, newA, newB) + state = (oldA, oldB, newA, newB) + if state == (0, 0, 1, 0) || state == (1, 0, 1, 1) || state == (1, 1, 0, 1) || state == (0, 1, 0, 0) + return 1 + elseif state == (0, 0, 0, 1) || state == (0, 1, 1, 1) || state == (1, 1, 1, 0) || state == (1, 0, 0, 0) + return -1 + elseif state == (0, 0, 0, 0) || state == (0, 1, 0, 1) || state == (1, 0, 1, 0) || state == (1, 1, 1, 1) + return 0 + else + return 0 # err is interpreted as no movement + end + end + # todo: warn about dups + # todo: warn if a variable appears in both observed and modified + qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0], + ModelingToolkit.MutatingFunctionalAffect((; qB), (; qA, hA, hB, cnt)) do x, o, i, c + x.hA = x.qA + x.hB = o.qB x.qA = 1 + x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) end, - affect_neg = ModelingToolkit.MutatingFunctionalAffect(modified=(; qA)) do x, o, i, c + affect_neg = ModelingToolkit.MutatingFunctionalAffect((; qB), (; qA, hA, hB, cnt)) do x, o, i, c + x.hA = x.qA + x.hB = o.qB x.qA = 0 - end) - qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(1000 * theta + π/2) ~ 0], - ModelingToolkit.MutatingFunctionalAffect(modified=(; qB)) do x, o, i, c + x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) + end; rootfind=SciMLBase.RightRootFind) + qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π/2) ~ 0], + ModelingToolkit.MutatingFunctionalAffect((; qA), (; qB, hA, hB, cnt)) do x, o, i, c + x.hA = o.qA + x.hB = x.qB x.qB = 1 + x.cnt += decoder(x.hA, x.hB, o.qA, x.qB) end, - affect_neg = ModelingToolkit.MutatingFunctionalAffect(modified=(; qB)) do x, o, i, c + affect_neg = ModelingToolkit.MutatingFunctionalAffect((; qA), (; qB, hA, hB, cnt)) do x, o, i, c + x.hA = o.qA + x.hB = x.qB x.qB = 0 - end) + x.cnt += decoder(x.hA, x.hB, o.qA, x.qB) + end; rootfind=SciMLBase.RightRootFind) @named sys = ODESystem(eqs, t, [theta, omega], params; continuous_events = [qAevt, qBevt]) ss = structural_simplify(sys) - prob = ODEProblem(ss, [theta => 0.0], (0.0, 1.0)) + prob = ODEProblem(ss, [theta => 0.0], (0.0, pi)) sol = solve(prob, Tsit5(); dtmax=0.01) + @test sol[cnt] == 197 # we get 2 pulses per phase cycle (cos 0 crossing) and we go to 100 cycles; we miss a few due to the initial state end From b96cd2eb273d640bf39d3a7f3c01a5e994cb7cc9 Mon Sep 17 00:00:00 2001 From: Ben Chung Date: Fri, 2 Aug 2024 16:42:55 -0700 Subject: [PATCH 04/47] More sanity checking --- src/systems/callbacks.jl | 28 ++++++++++++++++++++++++--- test/symbolic_events.jl | 42 ++++++++++++++++++++++++++++++++++++++-- 2 files changed, 65 insertions(+), 5 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index a57d5c006d..612d80536c 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -818,10 +818,10 @@ function compile_user_affect(affect::FunctionalAffect, cb, sys, dvs, ps; kwargs. end end -invalid_variables(sys, expr) = filter(x -> !any(isequal(x), all_symbols(sys)), vars(expr)) +invalid_variables(sys, expr) = filter(x -> !any(isequal(x), all_symbols(sys)), reduce(vcat, vars(expr); init=[])) function unassignable_variables(sys, expr) assignable_syms = vcat(unknowns(sys), parameters(sys)) - return filter(x -> !any(isequal(x), assignable_syms), vars(expr)) + return filter(x -> !any(isequal(x), assignable_syms), reduce(vcat, vars(expr); init=[])) end function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwargs...) @@ -832,6 +832,21 @@ function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwa call the affect method - test if it's OOP or IP using applicable unpack and apply the resulting values =# + function check_dups(syms, exprs) # = (syms_dedup, exprs_dedup) + seen = Set{Symbol}() + syms_dedup = []; exprs_dedup = [] + for (sym, exp) in Iterators.zip(syms, exprs) + if !in(sym, seen) + push!(syms_dedup, sym) + push!(exprs_dedup, exp) + push!(seen, sym) + else + @warn "Expression $(expr) is aliased as $sym, which has already been used. The first definition will be used." + end + end + return (syms_dedup, exprs_dedup) + end + obs_exprs = observed(affect) for oexpr in obs_exprs invalid_vars = invalid_variables(sys, oexpr) @@ -840,6 +855,7 @@ function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwa end end obs_syms = observed_syms(affect) + obs_syms, obs_exprs = check_dups(obs_syms, obs_exprs) obs_size = size.(obs_exprs) # we will generate a work buffer of a ComponentArray that maps obs_syms to arrays of size obs_size mod_exprs = modified(affect) @@ -849,12 +865,18 @@ function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwa end invalid_vars = unassignable_variables(sys, mexpr) if length(invalid_vars) > 0 - error("Observed equation $(mexpr) in affect refers to missing variable(s) $(invalid_vars); the variables may not have been added (e.g. if a component is missing) or they may have been reduced away.") + error("Modified equation $(mexpr) in affect refers to missing variable(s) $(invalid_vars); the variables may not have been added (e.g. if a component is missing) or they may have been reduced away.") end end mod_syms = modified_syms(affect) + mod_syms, mod_exprs = check_dups(mod_syms, mod_exprs) _, mod_og_val_fun = build_explicit_observed_function(sys, mod_exprs; return_inplace=true) + overlapping_syms = intersect(mod_syms, obs_syms) + if length(overlapping_syms) > 0 + @warn "The symbols $overlapping_syms are declared as both observed and modified; this is a code smell because it becomes easy to confuse them and assign/not assign a value." + end + # sanity checks done! now build the data and update function for observed values mkzero(sz) = if sz === () 0.0 else zeros(sz) end _, obs_fun = build_explicit_observed_function(sys, reduce(vcat, Symbolics.scalarize.(obs_exprs); init = []); return_inplace=true) diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index ca9f0ad9c3..d24b590970 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -1021,6 +1021,46 @@ end @test all(sol[temp][sol.t .> 1.0] .<= 0.79) && all(sol[temp][sol.t .> 1.0] .>= 0.49) end +@testset "MutatingFunctionalAffect errors and warnings" 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 + eqs = [ + D(temp) ~ furnace_on * furnace_power - temp^2 * leakage + ] + + furnace_off = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_off_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on), observed=(; furnace_on)) do x, o, i, c + x.furnace_on = false + end) + @named sys = ODESystem(eqs, t, [temp], params; continuous_events = [furnace_off]) + ss = structural_simplify(sys) + @test_logs (:warn, "The symbols Any[:furnace_on] are declared as both observed and modified; this is a code smell because it becomes easy to confuse them and assign/not assign a value.") prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + + @variables tempsq(t) # trivially eliminated + eqs = [ + tempsq ~ temp^2 + D(temp) ~ furnace_on * furnace_power - temp^2 * leakage + ] + + furnace_off = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_off_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on, tempsq), observed=(; furnace_on)) do x, o, i, c + x.furnace_on = false + end) + @named sys = ODESystem(eqs, t, [temp, tempsq], params; continuous_events = [furnace_off]) + ss = structural_simplify(sys) + @test_throws "refers to missing variable(s)" prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + + + @parameters not_actually_here + furnace_off = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_off_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on), observed=(; furnace_on, not_actually_here)) do x, o, i, c + x.furnace_on = false + end) + @named sys = ODESystem(eqs, t, [temp, tempsq], params; continuous_events = [furnace_off]) + ss = structural_simplify(sys) + @test_throws "refers to missing variable(s)" prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) +end + @testset "Quadrature" begin @variables theta(t) omega(t) params = @parameters qA=0 qB=0 hA=0 hB=0 cnt=0 @@ -1040,8 +1080,6 @@ end return 0 # err is interpreted as no movement end end - # todo: warn about dups - # todo: warn if a variable appears in both observed and modified qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0], ModelingToolkit.MutatingFunctionalAffect((; qB), (; qA, hA, hB, cnt)) do x, o, i, c x.hA = x.qA From eec24cfb95546fb5b8a2415a8cd3a2be0d60a6a8 Mon Sep 17 00:00:00 2001 From: Ben Chung Date: Fri, 2 Aug 2024 16:56:56 -0700 Subject: [PATCH 05/47] Document MutatingFunctionalAffect --- src/systems/callbacks.jl | 28 ++++++++++++++++++++++++---- 1 file changed, 24 insertions(+), 4 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index 612d80536c..fefff94442 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -72,10 +72,29 @@ function namespace_affect(affect::FunctionalAffect, s) end """ -`MutatingFunctionalAffect` differs from `FunctionalAffect` in two key ways: -* First, insetad of the `u` vector passed to `f` being a vector of indices into `integ.u` it's instead the result of evaluating `obs` at the current state, named as specified in `obs_syms`. This allows affects to easily access observed states and decouples affect inputs from the system structure. -* Second, it abstracts the assignment back to system states away. Instead of writing `integ.u[u.myvar] = [whatever]`, you instead declare in `mod_params` that you want to modify `myvar` and then either (out of place) return a named tuple with `myvar` or (in place) modify the associated element in the ComponentArray that's given. -Initially, we only support "flat" states in `modified`; these states will be marked as irreducible in the overarching system and they will simply be bulk assigned at mutation. In the future, this will be extended to perform a nonlinear solve to further decouple the affect from the system structure. + MutatingFunctionalAffect(f::Function; observed::NamedTuple, modified::NamedTuple, ctx) + +`MutatingFunctionalAffect` is a helper for writing affect functions that will compute observed values and +ensure that modified values are correctly written back into the system. The affect function `f` needs to have +one of three signatures: +* `f(observed::ComponentArray)` if the function only reads observed values back from the system, +* `f(observed::ComponentArray, modified::ComponentArray)` if the function also writes values (unknowns or parameters) into the system, +* `f(observed::ComponentArray, modified::ComponentArray, ctx)` if the function needs the user-defined context, +* `f(observed::ComponentArray, modified::ComponentArray, ctx, integrator)` if the function needs the low-level integrator. + +The function `f` will be called with `observed` and `modified` `ComponentArray`s that are derived from their respective `NamedTuple` definitions. +Each `NamedTuple` should map an expression to a symbol; for example if we pass `observed=(; x = a + b)` this will alias the result of executing `a+b` in the system as `x` +so the value of `a + b` will be accessible as `observed.x` in `f`. `modified` currently restricts symbolic expressions to only bare variables, so only tuples of the form +`(; x = y)` or `(; x)` (which aliases `x` as itself) are allowed. + +Both `observed` and `modified` will be automatically populated with the current values of their corresponding expressions on function entry. +The values in `modified` will be written back to the system after `f` returns. For example, if we want to update the value of `x` to be the result of `x + y` we could write + + MutatingFunctionalAffect(observed=(; x_plus_y = x + y), modified=(; x)) do o, m + m.x = o.x_plus_y + end + +The affect function updates the value at `x` in `modified` to be the result of evaluating `x + y` as passed in the observed values. """ @kwdef struct MutatingFunctionalAffect f::Any @@ -174,6 +193,7 @@ Affects (i.e. `affect` and `affect_neg`) can be specified as either: + `read_parameters` is a vector of the parameters that are *used* by `f!`. Their indices are passed to `f` in `p` similarly to the indices of `unknowns` passed in `u`. + `modified_parameters` is a vector of the parameters that are *modified* by `f!`. Note that a parameter will not appear in `p` if it only appears in `modified_parameters`; it must appear in both `parameters` and `modified_parameters` if it is used in the affect definition. + `ctx` is a user-defined context object passed to `f!` when invoked. This value is aliased for each problem. +* A [`MutatingFunctionalAffect`](@ref); refer to its documentation for details. """ struct SymbolicContinuousCallback eqs::Vector{Equation} From fd0125d36715693c0b7e2c5d5fe130e402c137ae Mon Sep 17 00:00:00 2001 From: Ben Chung Date: Fri, 2 Aug 2024 17:15:38 -0700 Subject: [PATCH 06/47] Flip modified and observed order; write docstring --- src/systems/callbacks.jl | 34 +++++++++++++------ test/symbolic_events.jl | 73 +++++++++++++++++++++++++++++++--------- 2 files changed, 80 insertions(+), 27 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index fefff94442..715b0e9a91 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -72,15 +72,16 @@ function namespace_affect(affect::FunctionalAffect, s) end """ - MutatingFunctionalAffect(f::Function; observed::NamedTuple, modified::NamedTuple, ctx) + MutatingFunctionalAffect(f::Function; modified::NamedTuple, observed::NamedTuple, ctx) `MutatingFunctionalAffect` is a helper for writing affect functions that will compute observed values and ensure that modified values are correctly written back into the system. The affect function `f` needs to have -one of three signatures: -* `f(observed::ComponentArray)` if the function only reads observed values back from the system, -* `f(observed::ComponentArray, modified::ComponentArray)` if the function also writes values (unknowns or parameters) into the system, -* `f(observed::ComponentArray, modified::ComponentArray, ctx)` if the function needs the user-defined context, -* `f(observed::ComponentArray, modified::ComponentArray, ctx, integrator)` if the function needs the low-level integrator. +one of four signatures: +* `f(modified::ComponentArray)` if the function only writes values (unknowns or parameters) to the system, +* `f(modified::ComponentArray, observed::ComponentArray)` if the function also reads observed values from the system, +* `f(modified::ComponentArray, observed::ComponentArray, ctx)` if the function needs the user-defined context, +* `f(modified::ComponentArray, observed::ComponentArray, ctx, integrator)` if the function needs the low-level integrator. +These will be checked in reverse order (that is, the four-argument version first, than the 3, etc). The function `f` will be called with `observed` and `modified` `ComponentArray`s that are derived from their respective `NamedTuple` definitions. Each `NamedTuple` should map an expression to a symbol; for example if we pass `observed=(; x = a + b)` this will alias the result of executing `a+b` in the system as `x` @@ -90,7 +91,7 @@ so the value of `a + b` will be accessible as `observed.x` in `f`. `modified` cu Both `observed` and `modified` will be automatically populated with the current values of their corresponding expressions on function entry. The values in `modified` will be written back to the system after `f` returns. For example, if we want to update the value of `x` to be the result of `x + y` we could write - MutatingFunctionalAffect(observed=(; x_plus_y = x + y), modified=(; x)) do o, m + MutatingFunctionalAffect(observed=(; x_plus_y = x + y), modified=(; x)) do m, o m.x = o.x_plus_y end @@ -109,11 +110,11 @@ MutatingFunctionalAffect(f::Function; observed::NamedTuple = NamedTuple{()}(()), modified::NamedTuple = NamedTuple{()}(()), ctx=nothing) = MutatingFunctionalAffect(f, collect(values(observed)), collect(keys(observed)), collect(values(modified)), collect(keys(modified)), ctx) -MutatingFunctionalAffect(f::Function, observed::NamedTuple; modified::NamedTuple = NamedTuple{()}(()), ctx=nothing) = +MutatingFunctionalAffect(f::Function, modified::NamedTuple; observed::NamedTuple = NamedTuple{()}(()), ctx=nothing) = MutatingFunctionalAffect(f, observed=observed, modified=modified, ctx=ctx) -MutatingFunctionalAffect(f::Function, observed::NamedTuple, modified::NamedTuple; ctx=nothing) = +MutatingFunctionalAffect(f::Function, modified::NamedTuple, observed::NamedTuple; ctx=nothing) = MutatingFunctionalAffect(f, observed=observed, modified=modified, ctx=ctx) -MutatingFunctionalAffect(f::Function, observed::NamedTuple, modified::NamedTuple, ctx) = +MutatingFunctionalAffect(f::Function, modified::NamedTuple, observed::NamedTuple, ctx) = MutatingFunctionalAffect(f, observed=observed, modified=modified, ctx=ctx) func(f::MutatingFunctionalAffect) = f.f @@ -925,7 +926,18 @@ function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwa obs_fun(obs_component_array, integ.u, integ.p..., integ.t) # let the user do their thing - user_affect(upd_component_array, obs_component_array, integ, ctx) + if applicable(user_affect, upd_component_array, obs_component_array, ctx, integ) + user_affect(upd_component_array, obs_component_array, ctx, integ) + elseif applicable(user_affect, upd_component_array, obs_component_array, ctx) + user_affect(upd_component_array, obs_component_array, ctx) + elseif applicable(user_affect, upd_component_array, obs_component_array) + user_affect(upd_component_array, obs_component_array) + elseif applicable(user_affect, upd_component_array) + user_affect(upd_component_array) + else + @error "User affect function $user_affect needs to implement one of the supported MutatingFunctionalAffect callback forms; see the MutatingFunctionalAffect docstring for more details" + user_affect(upd_component_array, obs_component_array, integ, ctx) # this WILL error but it'll give a more sensible message + end # write the new values back to the integrator upd_params_fun(integ, upd_params_view) diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index d24b590970..779f41471a 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -250,19 +250,19 @@ end m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (; x)) @test m isa ModelingToolkit.MutatingFunctionalAffect @test m.f == fmfa - @test isequal(m.obs, [x]) - @test m.obs_syms == [:x] - @test m.modified == [] - @test m.mod_syms == [] + @test isequal(m.obs, []) + @test m.obs_syms == [] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:x] @test m.ctx === nothing m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (; y=x)) @test m isa ModelingToolkit.MutatingFunctionalAffect @test m.f == fmfa - @test isequal(m.obs, [x]) - @test m.obs_syms == [:y] - @test m.modified == [] - @test m.mod_syms == [] + @test isequal(m.obs, []) + @test m.obs_syms == [] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:y] @test m.ctx === nothing m = ModelingToolkit.MutatingFunctionalAffect(fmfa; observed=(; y=x)) @@ -1013,7 +1013,48 @@ end ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on)) do x, o, i, c x.furnace_on = true end) - + @named sys = ODESystem(eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) + ss = structural_simplify(sys) + prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + sol = solve(prob, Tsit5(); dtmax=0.01) + @test all(sol[temp][sol.t .> 1.0] .<= 0.79) && all(sol[temp][sol.t .> 1.0] .>= 0.49) + + furnace_off = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_off_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on)) do x, o, i + x.furnace_on = false + end) + furnace_enable = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_on_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on)) do x, o, i + x.furnace_on = true + end) + @named sys = ODESystem(eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) + ss = structural_simplify(sys) + prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + sol = solve(prob, Tsit5(); dtmax=0.01) + @test all(sol[temp][sol.t .> 1.0] .<= 0.79) && all(sol[temp][sol.t .> 1.0] .>= 0.49) + + furnace_off = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_off_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on)) do x, o + x.furnace_on = false + end) + furnace_enable = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_on_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on)) do x, o + x.furnace_on = true + end) + @named sys = ODESystem(eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) + ss = structural_simplify(sys) + prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + sol = solve(prob, Tsit5(); dtmax=0.01) + @test all(sol[temp][sol.t .> 1.0] .<= 0.79) && all(sol[temp][sol.t .> 1.0] .>= 0.49) + + furnace_off = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_off_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on)) do x + x.furnace_on = false + end) + furnace_enable = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_on_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on)) do x + x.furnace_on = true + end) @named sys = ODESystem(eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) ss = structural_simplify(sys) prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) @@ -1029,7 +1070,7 @@ end ] furnace_off = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_off_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on), observed=(; furnace_on)) do x, o, i, c + ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on), observed=(; furnace_on)) do x, o, c, i x.furnace_on = false end) @named sys = ODESystem(eqs, t, [temp], params; continuous_events = [furnace_off]) @@ -1043,7 +1084,7 @@ end ] furnace_off = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_off_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on, tempsq), observed=(; furnace_on)) do x, o, i, c + ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on, tempsq), observed=(; furnace_on)) do x, o, c, i x.furnace_on = false end) @named sys = ODESystem(eqs, t, [temp, tempsq], params; continuous_events = [furnace_off]) @@ -1053,7 +1094,7 @@ end @parameters not_actually_here furnace_off = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_off_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on), observed=(; furnace_on, not_actually_here)) do x, o, i, c + ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on), observed=(; furnace_on, not_actually_here)) do x, o, c, i x.furnace_on = false end) @named sys = ODESystem(eqs, t, [temp, tempsq], params; continuous_events = [furnace_off]) @@ -1081,26 +1122,26 @@ end end end qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0], - ModelingToolkit.MutatingFunctionalAffect((; qB), (; qA, hA, hB, cnt)) do x, o, i, c + ModelingToolkit.MutatingFunctionalAffect((; qA, hA, hB, cnt), (; qB)) do x, o, i, c x.hA = x.qA x.hB = o.qB x.qA = 1 x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) end, - affect_neg = ModelingToolkit.MutatingFunctionalAffect((; qB), (; qA, hA, hB, cnt)) do x, o, i, c + affect_neg = ModelingToolkit.MutatingFunctionalAffect((; qA, hA, hB, cnt), (; qB)) do x, o, c, i x.hA = x.qA x.hB = o.qB x.qA = 0 x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) end; rootfind=SciMLBase.RightRootFind) qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π/2) ~ 0], - ModelingToolkit.MutatingFunctionalAffect((; qA), (; qB, hA, hB, cnt)) do x, o, i, c + ModelingToolkit.MutatingFunctionalAffect((; qB, hA, hB, cnt), (; qA)) do x, o, i, c x.hA = o.qA x.hB = x.qB x.qB = 1 x.cnt += decoder(x.hA, x.hB, o.qA, x.qB) end, - affect_neg = ModelingToolkit.MutatingFunctionalAffect((; qA), (; qB, hA, hB, cnt)) do x, o, i, c + affect_neg = ModelingToolkit.MutatingFunctionalAffect((; qB, hA, hB, cnt), (; qA)) do x, o, c, i x.hA = o.qA x.hB = x.qB x.qB = 0 From 9948de076b491a114b395f8a9168eaaebe0959b1 Mon Sep 17 00:00:00 2001 From: Ben Chung Date: Fri, 2 Aug 2024 17:17:12 -0700 Subject: [PATCH 07/47] Run formatter --- src/systems/callbacks.jl | 87 ++++++++++++------- test/symbolic_events.jl | 175 ++++++++++++++++++++++----------------- 2 files changed, 158 insertions(+), 104 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index 715b0e9a91..5c6f9b2801 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -106,16 +106,25 @@ The affect function updates the value at `x` in `modified` to be the result of e ctx::Any end -MutatingFunctionalAffect(f::Function; - observed::NamedTuple = NamedTuple{()}(()), - modified::NamedTuple = NamedTuple{()}(()), - ctx=nothing) = MutatingFunctionalAffect(f, collect(values(observed)), collect(keys(observed)), collect(values(modified)), collect(keys(modified)), ctx) -MutatingFunctionalAffect(f::Function, modified::NamedTuple; observed::NamedTuple = NamedTuple{()}(()), ctx=nothing) = - MutatingFunctionalAffect(f, observed=observed, modified=modified, ctx=ctx) -MutatingFunctionalAffect(f::Function, modified::NamedTuple, observed::NamedTuple; ctx=nothing) = - MutatingFunctionalAffect(f, observed=observed, modified=modified, ctx=ctx) -MutatingFunctionalAffect(f::Function, modified::NamedTuple, observed::NamedTuple, ctx) = - MutatingFunctionalAffect(f, observed=observed, modified=modified, ctx=ctx) +function MutatingFunctionalAffect(f::Function; + observed::NamedTuple = NamedTuple{()}(()), + modified::NamedTuple = NamedTuple{()}(()), + ctx = nothing) + MutatingFunctionalAffect(f, collect(values(observed)), collect(keys(observed)), + collect(values(modified)), collect(keys(modified)), ctx) +end +function MutatingFunctionalAffect(f::Function, modified::NamedTuple; + observed::NamedTuple = NamedTuple{()}(()), ctx = nothing) + MutatingFunctionalAffect(f, observed = observed, modified = modified, ctx = ctx) +end +function MutatingFunctionalAffect( + f::Function, modified::NamedTuple, observed::NamedTuple; ctx = nothing) + MutatingFunctionalAffect(f, observed = observed, modified = modified, ctx = ctx) +end +function MutatingFunctionalAffect( + f::Function, modified::NamedTuple, observed::NamedTuple, ctx) + MutatingFunctionalAffect(f, observed = observed, modified = modified, ctx = ctx) +end func(f::MutatingFunctionalAffect) = f.f context(a::MutatingFunctionalAffect) = a.ctx @@ -126,8 +135,9 @@ modified(a::MutatingFunctionalAffect) = a.modified modified_syms(a::MutatingFunctionalAffect) = a.mod_syms function Base.:(==)(a1::MutatingFunctionalAffect, a2::MutatingFunctionalAffect) - isequal(a1.f, a2.f) && isequal(a1.obs, a2.obs) && isequal(a1.modified, a2.modified) && - isequal(a1.obs_syms, a2.obs_syms) && isequal(a1.mod_syms, a2.mod_syms)&& isequal(a1.ctx, a2.ctx) + isequal(a1.f, a2.f) && isequal(a1.obs, a2.obs) && isequal(a1.modified, a2.modified) && + isequal(a1.obs_syms, a2.obs_syms) && isequal(a1.mod_syms, a2.mod_syms) && + isequal(a1.ctx, a2.ctx) end function Base.hash(a::MutatingFunctionalAffect, s::UInt) @@ -839,10 +849,13 @@ function compile_user_affect(affect::FunctionalAffect, cb, sys, dvs, ps; kwargs. end end -invalid_variables(sys, expr) = filter(x -> !any(isequal(x), all_symbols(sys)), reduce(vcat, vars(expr); init=[])) -function unassignable_variables(sys, expr) +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 = vcat(unknowns(sys), parameters(sys)) - return filter(x -> !any(isequal(x), assignable_syms), reduce(vcat, vars(expr); init=[])) + return filter( + x -> !any(isequal(x), assignable_syms), reduce(vcat, vars(expr); init = [])) end function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwargs...) @@ -855,7 +868,8 @@ function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwa =# function check_dups(syms, exprs) # = (syms_dedup, exprs_dedup) seen = Set{Symbol}() - syms_dedup = []; exprs_dedup = [] + syms_dedup = [] + exprs_dedup = [] for (sym, exp) in Iterators.zip(syms, exprs) if !in(sym, seen) push!(syms_dedup, sym) @@ -869,7 +883,7 @@ function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwa end obs_exprs = observed(affect) - for oexpr in obs_exprs + for oexpr in obs_exprs invalid_vars = invalid_variables(sys, oexpr) if length(invalid_vars) > 0 error("Observed equation $(oexpr) in affect refers to missing variable(s) $(invalid_vars); the variables may not have been added (e.g. if a component is missing).") @@ -880,7 +894,7 @@ function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwa obs_size = size.(obs_exprs) # we will generate a work buffer of a ComponentArray that maps obs_syms to arrays of size obs_size mod_exprs = modified(affect) - for mexpr in mod_exprs + for mexpr in mod_exprs if !is_observed(sys, mexpr) && parameter_index(sys, mexpr) === nothing error("Expression $mexpr cannot be assigned to; currently only unknowns and parameters may be updated by an affect.") end @@ -891,7 +905,8 @@ function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwa end mod_syms = modified_syms(affect) mod_syms, mod_exprs = check_dups(mod_syms, mod_exprs) - _, mod_og_val_fun = build_explicit_observed_function(sys, mod_exprs; return_inplace=true) + _, mod_og_val_fun = build_explicit_observed_function( + sys, mod_exprs; return_inplace = true) overlapping_syms = intersect(mod_syms, obs_syms) if length(overlapping_syms) > 0 @@ -899,21 +914,33 @@ function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwa end # sanity checks done! now build the data and update function for observed values - mkzero(sz) = if sz === () 0.0 else zeros(sz) end - _, obs_fun = build_explicit_observed_function(sys, reduce(vcat, Symbolics.scalarize.(obs_exprs); init = []); return_inplace=true) - obs_component_array = ComponentArrays.ComponentArray(NamedTuple{(obs_syms..., )}(mkzero.(obs_size))) + mkzero(sz) = + if sz === () + 0.0 + else + zeros(sz) + end + _, obs_fun = build_explicit_observed_function( + sys, reduce(vcat, Symbolics.scalarize.(obs_exprs); init = []); + return_inplace = true) + obs_component_array = ComponentArrays.ComponentArray(NamedTuple{(obs_syms...,)}(mkzero.(obs_size))) # okay so now to generate the stuff to assign it back into the system # note that we reorder the componentarray to make the views coherent wrt the base array mod_pairs = mod_exprs .=> mod_syms mod_param_pairs = filter(v -> is_parameter(sys, v[1]), mod_pairs) mod_unk_pairs = filter(v -> !is_parameter(sys, v[1]), mod_pairs) - _, mod_og_val_fun = build_explicit_observed_function(sys, reduce(vcat, [first.(mod_param_pairs); first.(mod_unk_pairs)]; init = []); return_inplace=true) - upd_params_fun = setu(sys, reduce(vcat, Symbolics.scalarize.(first.(mod_param_pairs)); init = [])) - upd_unk_fun = setu(sys, reduce(vcat, Symbolics.scalarize.(first.(mod_unk_pairs)); init = [])) - - upd_component_array = ComponentArrays.ComponentArray(NamedTuple{([last.(mod_param_pairs); last.(mod_unk_pairs)]...,)}( - [collect(mkzero(size(e)) for e in first.(mod_param_pairs)); + _, mod_og_val_fun = build_explicit_observed_function( + sys, reduce(vcat, [first.(mod_param_pairs); first.(mod_unk_pairs)]; init = []); + return_inplace = true) + upd_params_fun = setu( + sys, reduce(vcat, Symbolics.scalarize.(first.(mod_param_pairs)); init = [])) + upd_unk_fun = setu( + sys, reduce(vcat, Symbolics.scalarize.(first.(mod_unk_pairs)); init = [])) + + upd_component_array = ComponentArrays.ComponentArray(NamedTuple{([last.(mod_param_pairs); + last.(mod_unk_pairs)]...,)}( + [collect(mkzero(size(e)) for e in first.(mod_param_pairs)); collect(mkzero(size(e)) for e in first.(mod_unk_pairs))])) upd_params_view = view(upd_component_array, last.(mod_param_pairs)) upd_unks_view = view(upd_component_array, last.(mod_unk_pairs)) @@ -921,7 +948,7 @@ function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwa function (integ) # update the to-be-mutated values; this ensures that if you do a no-op then nothing happens mod_og_val_fun(upd_component_array, integ.u, integ.p..., integ.t) - + # update the observed values obs_fun(obs_component_array, integ.u, integ.p..., integ.t) @@ -934,7 +961,7 @@ function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwa user_affect(upd_component_array, obs_component_array) elseif applicable(user_affect, upd_component_array) user_affect(upd_component_array) - else + else @error "User affect function $user_affect needs to implement one of the supported MutatingFunctionalAffect callback forms; see the MutatingFunctionalAffect docstring for more details" user_affect(upd_component_array, obs_component_array, integ, ctx) # this WILL error but it'll give a more sensible message end diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index 779f41471a..7aba69dc7f 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -237,7 +237,7 @@ end @test m.modified == [] @test m.mod_syms == [] @test m.ctx === nothing - + m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (;)) @test m isa ModelingToolkit.MutatingFunctionalAffect @test m.f == fmfa @@ -246,7 +246,7 @@ end @test m.modified == [] @test m.mod_syms == [] @test m.ctx === nothing - + m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (; x)) @test m isa ModelingToolkit.MutatingFunctionalAffect @test m.f == fmfa @@ -255,8 +255,8 @@ end @test isequal(m.modified, [x]) @test m.mod_syms == [:x] @test m.ctx === nothing - - m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (; y=x)) + + m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (; y = x)) @test m isa ModelingToolkit.MutatingFunctionalAffect @test m.f == fmfa @test isequal(m.obs, []) @@ -264,8 +264,8 @@ end @test isequal(m.modified, [x]) @test m.mod_syms == [:y] @test m.ctx === nothing - - m = ModelingToolkit.MutatingFunctionalAffect(fmfa; observed=(; y=x)) + + m = ModelingToolkit.MutatingFunctionalAffect(fmfa; observed = (; y = x)) @test m isa ModelingToolkit.MutatingFunctionalAffect @test m.f == fmfa @test isequal(m.obs, [x]) @@ -273,8 +273,8 @@ end @test m.modified == [] @test m.mod_syms == [] @test m.ctx === nothing - - m = ModelingToolkit.MutatingFunctionalAffect(fmfa; modified=(; x)) + + m = ModelingToolkit.MutatingFunctionalAffect(fmfa; modified = (; x)) @test m isa ModelingToolkit.MutatingFunctionalAffect @test m.f == fmfa @test isequal(m.obs, []) @@ -283,7 +283,7 @@ end @test m.mod_syms == [:x] @test m.ctx === nothing - m = ModelingToolkit.MutatingFunctionalAffect(fmfa; modified=(; y=x)) + m = ModelingToolkit.MutatingFunctionalAffect(fmfa; modified = (; y = x)) @test m isa ModelingToolkit.MutatingFunctionalAffect @test m.f == fmfa @test isequal(m.obs, []) @@ -291,7 +291,7 @@ end @test isequal(m.modified, [x]) @test m.mod_syms == [:y] @test m.ctx === nothing - + m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (; x), (; x)) @test m isa ModelingToolkit.MutatingFunctionalAffect @test m.f == fmfa @@ -300,8 +300,8 @@ end @test isequal(m.modified, [x]) @test m.mod_syms == [:x] @test m.ctx === nothing - - m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (; y=x), (; y=x)) + + m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (; y = x), (; y = x)) @test m isa ModelingToolkit.MutatingFunctionalAffect @test m.f == fmfa @test isequal(m.obs, [x]) @@ -309,8 +309,9 @@ end @test isequal(m.modified, [x]) @test m.mod_syms == [:y] @test m.ctx === nothing - - m = ModelingToolkit.MutatingFunctionalAffect(fmfa; modified=(; y=x), observed=(; y=x)) + + m = ModelingToolkit.MutatingFunctionalAffect( + fmfa; modified = (; y = x), observed = (; y = x)) @test m isa ModelingToolkit.MutatingFunctionalAffect @test m.f == fmfa @test isequal(m.obs, [x]) @@ -318,8 +319,9 @@ end @test isequal(m.modified, [x]) @test m.mod_syms == [:y] @test m.ctx === nothing - - m = ModelingToolkit.MutatingFunctionalAffect(fmfa; modified=(; y=x), observed=(; y=x), ctx=3) + + m = ModelingToolkit.MutatingFunctionalAffect( + fmfa; modified = (; y = x), observed = (; y = x), ctx = 3) @test m isa ModelingToolkit.MutatingFunctionalAffect @test m.f == fmfa @test isequal(m.obs, [x]) @@ -327,7 +329,7 @@ end @test isequal(m.modified, [x]) @test m.mod_syms == [:y] @test m.ctx === 3 - + m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (; x), (; x), 3) @test m isa ModelingToolkit.MutatingFunctionalAffect @test m.f == fmfa @@ -1005,151 +1007,176 @@ end D(temp) ~ furnace_on * furnace_power - temp^2 * leakage ] - furnace_off = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_off_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on)) do x, o, i, c + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o, i, c x.furnace_on = false end) - furnace_enable = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_on_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on)) do x, o, i, c + furnace_enable = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_on_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o, i, c x.furnace_on = true end) - @named sys = ODESystem(eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) + @named sys = ODESystem( + eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) ss = structural_simplify(sys) prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) - sol = solve(prob, Tsit5(); dtmax=0.01) + sol = solve(prob, Tsit5(); dtmax = 0.01) @test all(sol[temp][sol.t .> 1.0] .<= 0.79) && all(sol[temp][sol.t .> 1.0] .>= 0.49) - furnace_off = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_off_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on)) do x, o, i + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o, i x.furnace_on = false end) - furnace_enable = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_on_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on)) do x, o, i + furnace_enable = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_on_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o, i x.furnace_on = true end) - @named sys = ODESystem(eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) + @named sys = ODESystem( + eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) ss = structural_simplify(sys) prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) - sol = solve(prob, Tsit5(); dtmax=0.01) + sol = solve(prob, Tsit5(); dtmax = 0.01) @test all(sol[temp][sol.t .> 1.0] .<= 0.79) && all(sol[temp][sol.t .> 1.0] .>= 0.49) - furnace_off = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_off_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on)) do x, o + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o x.furnace_on = false end) - furnace_enable = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_on_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on)) do x, o + furnace_enable = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_on_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o x.furnace_on = true end) - @named sys = ODESystem(eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) + @named sys = ODESystem( + eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) ss = structural_simplify(sys) prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) - sol = solve(prob, Tsit5(); dtmax=0.01) + sol = solve(prob, Tsit5(); dtmax = 0.01) @test all(sol[temp][sol.t .> 1.0] .<= 0.79) && all(sol[temp][sol.t .> 1.0] .>= 0.49) - furnace_off = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_off_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on)) do x + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x x.furnace_on = false end) - furnace_enable = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_on_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on)) do x + furnace_enable = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_on_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x x.furnace_on = true end) - @named sys = ODESystem(eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) + @named sys = ODESystem( + eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) ss = structural_simplify(sys) prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) - sol = solve(prob, Tsit5(); dtmax=0.01) + sol = solve(prob, Tsit5(); dtmax = 0.01) @test all(sol[temp][sol.t .> 1.0] .<= 0.79) && all(sol[temp][sol.t .> 1.0] .>= 0.49) end -@testset "MutatingFunctionalAffect errors and warnings" begin +@testset "MutatingFunctionalAffect errors and warnings" 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 eqs = [ D(temp) ~ furnace_on * furnace_power - temp^2 * leakage ] - furnace_off = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_off_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on), observed=(; furnace_on)) do x, o, c, i + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.MutatingFunctionalAffect( + modified = (; furnace_on), observed = (; furnace_on)) do x, o, c, i x.furnace_on = false end) @named sys = ODESystem(eqs, t, [temp], params; continuous_events = [furnace_off]) ss = structural_simplify(sys) - @test_logs (:warn, "The symbols Any[:furnace_on] are declared as both observed and modified; this is a code smell because it becomes easy to confuse them and assign/not assign a value.") prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + @test_logs (:warn, + "The symbols Any[:furnace_on] are declared as both observed and modified; this is a code smell because it becomes easy to confuse them and assign/not assign a value.") prob=ODEProblem( + ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) @variables tempsq(t) # trivially eliminated - eqs = [ - tempsq ~ temp^2 - D(temp) ~ furnace_on * furnace_power - temp^2 * leakage - ] + eqs = [tempsq ~ temp^2 + D(temp) ~ furnace_on * furnace_power - temp^2 * leakage] - furnace_off = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_off_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on, tempsq), observed=(; furnace_on)) do x, o, c, i + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.MutatingFunctionalAffect( + modified = (; furnace_on, tempsq), observed = (; furnace_on)) do x, o, c, i x.furnace_on = false end) - @named sys = ODESystem(eqs, t, [temp, tempsq], params; continuous_events = [furnace_off]) + @named sys = ODESystem( + eqs, t, [temp, tempsq], params; continuous_events = [furnace_off]) ss = structural_simplify(sys) - @test_throws "refers to missing variable(s)" prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + @test_throws "refers to missing variable(s)" prob=ODEProblem( + ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) - @parameters not_actually_here - furnace_off = ModelingToolkit.SymbolicContinuousCallback([temp ~ furnace_off_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified=(; furnace_on), observed=(; furnace_on, not_actually_here)) do x, o, c, i + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on), + observed = (; furnace_on, not_actually_here)) do x, o, c, i x.furnace_on = false end) - @named sys = ODESystem(eqs, t, [temp, tempsq], params; continuous_events = [furnace_off]) + @named sys = ODESystem( + eqs, t, [temp, tempsq], params; continuous_events = [furnace_off]) ss = structural_simplify(sys) - @test_throws "refers to missing variable(s)" prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + @test_throws "refers to missing variable(s)" prob=ODEProblem( + ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) end -@testset "Quadrature" begin +@testset "Quadrature" begin @variables theta(t) omega(t) params = @parameters qA=0 qB=0 hA=0 hB=0 cnt=0 - eqs = [ - D(theta) ~ omega - omega ~ 1.0 - ] + eqs = [D(theta) ~ omega + omega ~ 1.0] function decoder(oldA, oldB, newA, newB) state = (oldA, oldB, newA, newB) - if state == (0, 0, 1, 0) || state == (1, 0, 1, 1) || state == (1, 1, 0, 1) || state == (0, 1, 0, 0) + if state == (0, 0, 1, 0) || state == (1, 0, 1, 1) || state == (1, 1, 0, 1) || + state == (0, 1, 0, 0) return 1 - elseif state == (0, 0, 0, 1) || state == (0, 1, 1, 1) || state == (1, 1, 1, 0) || state == (1, 0, 0, 0) + elseif state == (0, 0, 0, 1) || state == (0, 1, 1, 1) || state == (1, 1, 1, 0) || + state == (1, 0, 0, 0) return -1 - elseif state == (0, 0, 0, 0) || state == (0, 1, 0, 1) || state == (1, 0, 1, 0) || state == (1, 1, 1, 1) + elseif state == (0, 0, 0, 0) || state == (0, 1, 0, 1) || state == (1, 0, 1, 0) || + state == (1, 1, 1, 1) return 0 else return 0 # err is interpreted as no movement end end - qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0], + qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0], ModelingToolkit.MutatingFunctionalAffect((; qA, hA, hB, cnt), (; qB)) do x, o, i, c x.hA = x.qA x.hB = o.qB x.qA = 1 x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) end, - affect_neg = ModelingToolkit.MutatingFunctionalAffect((; qA, hA, hB, cnt), (; qB)) do x, o, c, i + affect_neg = ModelingToolkit.MutatingFunctionalAffect( + (; qA, hA, hB, cnt), (; qB)) do x, o, c, i x.hA = x.qA x.hB = o.qB x.qA = 0 x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) - end; rootfind=SciMLBase.RightRootFind) - qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π/2) ~ 0], + end; rootfind = SciMLBase.RightRootFind) + qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π / 2) ~ 0], ModelingToolkit.MutatingFunctionalAffect((; qB, hA, hB, cnt), (; qA)) do x, o, i, c x.hA = o.qA x.hB = x.qB x.qB = 1 x.cnt += decoder(x.hA, x.hB, o.qA, x.qB) end, - affect_neg = ModelingToolkit.MutatingFunctionalAffect((; qB, hA, hB, cnt), (; qA)) do x, o, c, i + affect_neg = ModelingToolkit.MutatingFunctionalAffect( + (; qB, hA, hB, cnt), (; qA)) do x, o, c, i x.hA = o.qA x.hB = x.qB x.qB = 0 x.cnt += decoder(x.hA, x.hB, o.qA, x.qB) - end; rootfind=SciMLBase.RightRootFind) - @named sys = ODESystem(eqs, t, [theta, omega], params; continuous_events = [qAevt, qBevt]) + end; rootfind = SciMLBase.RightRootFind) + @named sys = ODESystem( + eqs, t, [theta, omega], params; continuous_events = [qAevt, qBevt]) ss = structural_simplify(sys) prob = ODEProblem(ss, [theta => 0.0], (0.0, pi)) - sol = solve(prob, Tsit5(); dtmax=0.01) + sol = solve(prob, Tsit5(); dtmax = 0.01) @test sol[cnt] == 197 # we get 2 pulses per phase cycle (cos 0 crossing) and we go to 100 cycles; we miss a few due to the initial state end From 61cf6762c43e595448ed6f37f37bd889db663b10 Mon Sep 17 00:00:00 2001 From: Ben Chung Date: Fri, 2 Aug 2024 18:44:18 -0700 Subject: [PATCH 08/47] FIx SCC reconstruction --- src/systems/callbacks.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index 5c6f9b2801..8ea7d9f506 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -284,13 +284,15 @@ 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(af::MutatingFunctionalAffect, s) = namespace_affect(af, s) namespace_affects(::Nothing, s) = nothing 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)) + SymbolicContinuousCallback(; + eqs = namespace_equation.(equations(cb), (s,)), + affect = namespace_affects(affects(cb), s), + affect_neg = namespace_affects(affect_negs(cb), s), + rootfind = cb.rootfind) end """ From 8edef14b63637a4132e5139540d3ec89fa9d39fd Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Mon, 19 Aug 2024 07:49:33 -0700 Subject: [PATCH 09/47] Implement initialize and finalize affects for symbolic callbacks --- src/systems/callbacks.jl | 206 +++++++++++++++++++++++++++++---------- 1 file changed, 154 insertions(+), 52 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index 8ea7d9f506..d14fd50691 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -104,28 +104,38 @@ The affect function updates the value at `x` in `modified` to be the result of e modified::Vector mod_syms::Vector{Symbol} ctx::Any + skip_checks::Bool end function MutatingFunctionalAffect(f::Function; observed::NamedTuple = NamedTuple{()}(()), modified::NamedTuple = NamedTuple{()}(()), - ctx = nothing) - MutatingFunctionalAffect(f, collect(values(observed)), collect(keys(observed)), - collect(values(modified)), collect(keys(modified)), ctx) + ctx = nothing, + skip_checks = false) + MutatingFunctionalAffect(f, + collect(values(observed)), collect(keys(observed)), + collect(values(modified)), collect(keys(modified)), + ctx, skip_checks) end function MutatingFunctionalAffect(f::Function, modified::NamedTuple; - observed::NamedTuple = NamedTuple{()}(()), ctx = nothing) - MutatingFunctionalAffect(f, observed = observed, modified = modified, ctx = ctx) + observed::NamedTuple = NamedTuple{()}(()), ctx = nothing, skip_checks=false) + MutatingFunctionalAffect(f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) end function MutatingFunctionalAffect( - f::Function, modified::NamedTuple, observed::NamedTuple; ctx = nothing) - MutatingFunctionalAffect(f, observed = observed, modified = modified, ctx = ctx) + f::Function, modified::NamedTuple, observed::NamedTuple; ctx = nothing, skip_checks=false) + MutatingFunctionalAffect(f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) end function MutatingFunctionalAffect( - f::Function, modified::NamedTuple, observed::NamedTuple, ctx) - MutatingFunctionalAffect(f, observed = observed, modified = modified, ctx = ctx) + f::Function, modified::NamedTuple, observed::NamedTuple, ctx; skip_checks=false) + MutatingFunctionalAffect(f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) end +function Base.show(io::IO, mfa::MutatingFunctionalAffect) + obs_vals = join(map((ob,nm) -> "$ob => $nm", mfa.obs, mfa.obs_syms), ", ") + mod_vals = join(map((md,nm) -> "$md => $nm", mfa.modified, mfa.mod_syms), ", ") + affect = mfa.f + print(io, "MutatingFunctionalAffect(observed: [$obs_vals], modified: [$mod_vals], affect:$affect)") +end func(f::MutatingFunctionalAffect) = f.f context(a::MutatingFunctionalAffect) = a.ctx observed(a::MutatingFunctionalAffect) = a.obs @@ -208,12 +218,19 @@ Affects (i.e. `affect` and `affect_neg`) can be specified as either: """ struct SymbolicContinuousCallback eqs::Vector{Equation} + initialize::Union{Vector{Equation}, FunctionalAffect, MutatingFunctionalAffect} + finalize::Union{Vector{Equation}, FunctionalAffect, MutatingFunctionalAffect} affect::Union{Vector{Equation}, FunctionalAffect, MutatingFunctionalAffect} affect_neg::Union{Vector{Equation}, FunctionalAffect, MutatingFunctionalAffect, Nothing} rootfind::SciMLBase.RootfindOpt - function SymbolicContinuousCallback(; eqs::Vector{Equation}, affect = NULL_AFFECT, - affect_neg = affect, rootfind = SciMLBase.LeftRootFind) - new(eqs, make_affect(affect), make_affect(affect_neg), rootfind) + function SymbolicContinuousCallback(; + eqs::Vector{Equation}, + affect = NULL_AFFECT, + affect_neg = affect, + rootfind = SciMLBase.LeftRootFind, + initialize=NULL_AFFECT, + finalize=NULL_AFFECT) + new(eqs, initialize, finalize, make_affect(affect), make_affect(affect_neg), rootfind) end # Default affect to nothing end make_affect(affect) = affect @@ -221,18 +238,81 @@ 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.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) 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(cb.affect, s) s = foldr(hash, cb.eqs, init = s) - s = cb.affect isa AbstractVector ? foldr(hash, cb.affect, init = s) : hash(cb.affect, s) - s = cb.affect_neg isa AbstractVector ? foldr(hash, cb.affect_neg, init = s) : - hash(cb.affect_neg, 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) hash(cb.rootfind, s) end + +function Base.show(io::IO, cb::SymbolicContinuousCallback) + indent = get(io, :indent, 0) + iio = IOContext(io, :indent => indent+1) + print(io, "SymbolicContinuousCallback(") + print(iio, "Equations:") + show(iio, equations(cb)) + print(iio, "; ") + if affects(cb) != NULL_AFFECT + print(iio, "Affect:") + show(iio, affects(cb)) + print(iio, ", ") + end + if affect_negs(cb) != NULL_AFFECT + print(iio, "Negative-edge affect:") + show(iio, affect_negs(cb)) + print(iio, ", ") + end + if initialize_affects(cb) != NULL_AFFECT + print(iio, "Initialization affect:") + show(iio, initialize_affects(cb)) + print(iio, ", ") + end + if finalize_affects(cb) != NULL_AFFECT + print(iio, "Finalization affect:") + show(iio, finalize_affects(cb)) + end + print(iio, ")") +end + +function Base.show(io::IO, mime::MIME"text/plain", cb::SymbolicContinuousCallback) + indent = get(io, :indent, 0) + iio = IOContext(io, :indent => indent+1) + println(io, "SymbolicContinuousCallback:") + println(iio, "Equations:") + show(iio, mime, equations(cb)) + print(iio, "\n") + if affects(cb) != NULL_AFFECT + println(iio, "Affect:") + show(iio, mime, affects(cb)) + print(iio, "\n") + end + if affect_negs(cb) != NULL_AFFECT + println(iio, "Negative-edge affect:") + show(iio, mime, affect_negs(cb)) + print(iio, "\n") + end + if initialize_affects(cb) != NULL_AFFECT + println(iio, "Initialization affect:") + show(iio, mime, initialize_affects(cb)) + print(iio, "\n") + end + if finalize_affects(cb) != NULL_AFFECT + 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}) @@ -246,14 +326,14 @@ 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; - affect_neg = affect, rootfind = SciMLBase.LeftRootFind) + affect_neg = affect, rootfind = SciMLBase.LeftRootFind, initialize = NULL_AFFECT, finalize = NULL_AFFECT) SymbolicContinuousCallback( - eqs = [eqs], affect = affect, affect_neg = affect_neg, rootfind = rootfind) + eqs = [eqs], affect = affect, affect_neg = affect_neg, rootfind = rootfind, initialize=initialize, finalize=finalize) end function SymbolicContinuousCallback(eqs::Vector{Equation}, affect = NULL_AFFECT; - affect_neg = affect, rootfind = SciMLBase.LeftRootFind) + affect_neg = affect, rootfind = SciMLBase.LeftRootFind, initialize = NULL_AFFECT, finalize = NULL_AFFECT) SymbolicContinuousCallback( - eqs = eqs, affect = affect, affect_neg = affect_neg, rootfind = rootfind) + eqs = eqs, affect = affect, affect_neg = affect_neg, rootfind = rootfind, initialize=initialize, finalize=finalize) end SymbolicContinuousCallbacks(cb::SymbolicContinuousCallback) = [cb] @@ -282,6 +362,16 @@ function affect_negs(cbs::Vector{SymbolicContinuousCallback}) mapreduce(affect_negs, vcat, cbs, init = Equation[]) 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.initialize +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(af::MutatingFunctionalAffect, s) = namespace_affect(af, s) @@ -292,6 +382,8 @@ function namespace_callback(cb::SymbolicContinuousCallback, s)::SymbolicContinuo 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 @@ -681,8 +773,9 @@ function generate_single_rootfinding_callback( initfn = SciMLBase.INITIALIZE_DEFAULT end return ContinuousCallback( - cond, affect_function.affect, affect_function.affect_neg, - rootfind = cb.rootfind, initialize = initfn) + cond, affect_function.affect, affect_function.affect_neg, rootfind = cb.rootfind, + initialize = isnothing(affect_function.initialize) ? SciMLBase.INITIALIZE_DEFAULT : (c, u, t, i) -> affect_function.initialize(i), + finalize = isnothing(affect_function.finalize) ? SciMLBase.FINALIZE_DEFAULT : (c, u, t, i) -> affect_function.finalize(i)) end function generate_vector_rootfinding_callback( @@ -702,13 +795,12 @@ function generate_vector_rootfinding_callback( _, rf_ip = generate_custom_function( sys, rhss, dvs, ps; expression = Val{false}, kwargs...) - affect_functions = @NamedTuple{affect::Function, affect_neg::Union{Function, Nothing}}[compile_affect_fn( - cb, - sys, - dvs, - ps, - kwargs) - for cb in cbs] + 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 @@ -734,25 +826,27 @@ function generate_vector_rootfinding_callback( affect_neg(integ) end end - if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing - save_idxs = mapreduce( - cb -> get(ic.callback_to_clocks, cb, Int[]), vcat, cbs; init = Int[]) - initfn = if isempty(save_idxs) - SciMLBase.INITIALIZE_DEFAULT + function handle_optional_setup_fn(funs, default) + if all(isnothing, funs) + return default else - let save_idxs = save_idxs - function (cb, u, t, integrator) - for idx in save_idxs - SciMLBase.save_discretes!(integrator, idx) + 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 - else - initfn = SciMLBase.INITIALIZE_DEFAULT end + initialize = handle_optional_setup_fn(map(fn -> fn.initialize, affect_functions), SciMLBase.INITIALIZE_DEFAULT) + 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 = initfn) + cond, affect, affect_neg, length(eqs), rootfind = rootfind, initialize = initialize, finalize = finalize) end """ @@ -762,15 +856,23 @@ function compile_affect_fn(cb, sys::AbstractODESystem, 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) + if isnothing(aff) + return nothing + else + affspr = compile_affect(aff, cb, sys, dvs, ps; expression = Val{true}, kwargs...) + @show affspr + return compile_affect(aff, cb, sys, dvs, ps; expression = Val{false}, kwargs...) + end + end if eq_neg_aff === eq_aff affect_neg = affect - elseif isnothing(eq_neg_aff) - affect_neg = nothing else - affect_neg = compile_affect( - eq_neg_aff, cb, sys, dvs, ps; expression = Val{false}, kwargs...) + affect_neg = compile_optional_affect(eq_neg_aff) end - (affect = affect, affect_neg = affect_neg) + initialize = compile_optional_affect(initialize_affects(cb)) + finalize = compile_optional_affect(finalize_affects(cb)) + (affect = affect, affect_neg = affect_neg, initialize = initialize, finalize = finalize) end function generate_rootfinding_callback(cbs, sys::AbstractODESystem, dvs = unknowns(sys), @@ -877,7 +979,7 @@ function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwa push!(syms_dedup, sym) push!(exprs_dedup, exp) push!(seen, sym) - else + elseif !affect.skip_checks @warn "Expression $(expr) is aliased as $sym, which has already been used. The first definition will be used." end end @@ -887,7 +989,7 @@ function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwa obs_exprs = observed(affect) for oexpr in obs_exprs invalid_vars = invalid_variables(sys, oexpr) - if length(invalid_vars) > 0 + if length(invalid_vars) > 0 && !affect.skip_checks error("Observed equation $(oexpr) in affect refers to missing variable(s) $(invalid_vars); the variables may not have been added (e.g. if a component is missing).") end end @@ -897,11 +999,11 @@ function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwa mod_exprs = modified(affect) for mexpr in mod_exprs - if !is_observed(sys, mexpr) && parameter_index(sys, mexpr) === nothing - error("Expression $mexpr cannot be assigned to; currently only unknowns and parameters may be updated by an affect.") + if !is_observed(sys, mexpr) && parameter_index(sys, mexpr) === nothing && !affect.skip_checks + @warn ("Expression $mexpr cannot be assigned to; currently only unknowns and parameters may be updated by an affect.") end invalid_vars = unassignable_variables(sys, mexpr) - if length(invalid_vars) > 0 + if length(invalid_vars) > 0 && !affect.skip_checks error("Modified equation $(mexpr) in affect refers to missing variable(s) $(invalid_vars); the variables may not have been added (e.g. if a component is missing) or they may have been reduced away.") end end @@ -911,7 +1013,7 @@ function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwa sys, mod_exprs; return_inplace = true) overlapping_syms = intersect(mod_syms, obs_syms) - if length(overlapping_syms) > 0 + if length(overlapping_syms) > 0 && !affect.skip_checks @warn "The symbols $overlapping_syms are declared as both observed and modified; this is a code smell because it becomes easy to confuse them and assign/not assign a value." end From 3c7afd8c40b526dd6f139aa2d119095ab6debeb7 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Mon, 19 Aug 2024 08:08:13 -0700 Subject: [PATCH 10/47] Test some simple initialization affects --- test/symbolic_events.jl | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index 7aba69dc7f..b0f230d75e 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -1074,6 +1074,26 @@ end prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) sol = solve(prob, Tsit5(); dtmax = 0.01) @test all(sol[temp][sol.t .> 1.0] .<= 0.79) && all(sol[temp][sol.t .> 1.0] .>= 0.49) + + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x + x.furnace_on = false + end; initialize = ModelingToolkit.MutatingFunctionalAffect(modified = (; temp)) do x + x.temp = 0.2 + end) + furnace_enable = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_on_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o, c, i + x.furnace_on = true + end) + @named sys = ODESystem( + eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) + ss = structural_simplify(sys) + prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + sol = solve(prob, Tsit5(); dtmax = 0.01) + @test all(sol[temp][sol.t .> 1.0] .<= 0.79) && all(sol[temp][sol.t .> 1.0] .>= 0.49) + @test all(sol[temp][sol.t .!= 0.0] .<= 0.79) && all(sol[temp][sol.t .!= 0.0] .>= 0.2) end @testset "MutatingFunctionalAffect errors and warnings" begin From 95956f363eaa8bbc5b492452946c6c2a21e35bbd Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Thu, 22 Aug 2024 07:17:50 -0700 Subject: [PATCH 11/47] Properly pass skip_checks through --- src/systems/callbacks.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index d14fd50691..ca0a64b452 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -165,7 +165,8 @@ function namespace_affect(affect::MutatingFunctionalAffect, s) observed_syms(affect), renamespace.((s,), modified(affect)), modified_syms(affect), - context(affect)) + context(affect), + affect.skip_checks) end function has_functional_affect(cb) From 4f928ae57e33b2d8bb5e52faed229941b5e589af Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Thu, 22 Aug 2024 07:35:38 -0700 Subject: [PATCH 12/47] Support time-indexed parameters --- src/systems/callbacks.jl | 16 ++++++++++++++-- src/systems/index_cache.jl | 2 +- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index ca0a64b452..279d31483e 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -963,7 +963,7 @@ function unassignable_variables(sys, expr) x -> !any(isequal(x), assignable_syms), reduce(vcat, vars(expr); init = [])) end -function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwargs...) +function compile_user_affect(affect::MutatingFunctionalAffect, cb, sys, dvs, ps; kwargs...) #= Implementation sketch: generate observed function (oop), should save to a component array under obs_syms @@ -1049,6 +1049,13 @@ function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwa collect(mkzero(size(e)) for e in first.(mod_unk_pairs))])) upd_params_view = view(upd_component_array, last.(mod_param_pairs)) upd_unks_view = view(upd_component_array, last.(mod_unk_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) # update the to-be-mutated values; this ensures that if you do a no-op then nothing happens @@ -1074,11 +1081,16 @@ function compile_user_affect(affect::MutatingFunctionalAffect, sys, dvs, ps; kwa # write the new values back to the integrator upd_params_fun(integ, upd_params_view) upd_unk_fun(integ, upd_unks_view) + + + for idx in save_idxs + SciMLBase.save_discretes!(integ, idx) + end end end end -function compile_affect(affect::FunctionalAffect, cb, sys, dvs, ps; kwargs...) +function compile_affect(affect::Union{FunctionalAffect, MutatingFunctionalAffect}, cb, sys, dvs, ps; kwargs...) compile_user_affect(affect, cb, sys, dvs, ps; kwargs...) end diff --git a/src/systems/index_cache.jl b/src/systems/index_cache.jl index 00f7837407..55d819990b 100644 --- a/src/systems/index_cache.jl +++ b/src/systems/index_cache.jl @@ -126,7 +126,7 @@ function IndexCache(sys::AbstractSystem) for affect in affs if affect isa Equation is_parameter(sys, affect.lhs) && push!(discs, affect.lhs) - elseif affect isa FunctionalAffect + elseif affect isa FunctionalAffect || affect isa MutatingFunctionalAffect union!(discs, unwrap.(discretes(affect))) else error("Unhandled affect type $(typeof(affect))") From 3eca9b96567008469ee6963745a6a4f2ea64ed50 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Tue, 27 Aug 2024 09:28:17 -0700 Subject: [PATCH 13/47] Fix bugs relating to array arguments to callbacks --- src/systems/callbacks.jl | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index 279d31483e..884d1d87f3 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -958,9 +958,10 @@ 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 = vcat(unknowns(sys), parameters(sys)) + assignable_syms = reduce(vcat, Symbolics.scalarize.(vcat(unknowns(sys), parameters(sys))); init=[]) + written = reduce(vcat, Symbolics.scalarize.(vars(expr)); init = []) return filter( - x -> !any(isequal(x), assignable_syms), reduce(vcat, vars(expr); init = [])) + x -> !any(isequal(x), assignable_syms), written) end function compile_user_affect(affect::MutatingFunctionalAffect, cb, sys, dvs, ps; kwargs...) @@ -1000,7 +1001,10 @@ function compile_user_affect(affect::MutatingFunctionalAffect, cb, sys, dvs, ps; mod_exprs = modified(affect) for mexpr in mod_exprs - if !is_observed(sys, mexpr) && parameter_index(sys, mexpr) === nothing && !affect.skip_checks + if affect.skip_checks + continue + end + if !is_variable(sys, mexpr) && parameter_index(sys, mexpr) === nothing && !affect.skip_checks @warn ("Expression $mexpr cannot be assigned to; currently only unknowns and parameters may be updated by an affect.") end invalid_vars = unassignable_variables(sys, mexpr) @@ -1036,7 +1040,7 @@ function compile_user_affect(affect::MutatingFunctionalAffect, cb, sys, dvs, ps; mod_param_pairs = filter(v -> is_parameter(sys, v[1]), mod_pairs) mod_unk_pairs = filter(v -> !is_parameter(sys, v[1]), mod_pairs) _, mod_og_val_fun = build_explicit_observed_function( - sys, reduce(vcat, [first.(mod_param_pairs); first.(mod_unk_pairs)]; init = []); + sys, reduce(vcat, Symbolics.scalarize.([first.(mod_param_pairs); first.(mod_unk_pairs)]); init = []); return_inplace = true) upd_params_fun = setu( sys, reduce(vcat, Symbolics.scalarize.(first.(mod_param_pairs)); init = [])) From c41e7d4677077f69024c0f4f0a8e9ef6f83b2d38 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Tue, 27 Aug 2024 18:27:22 -0700 Subject: [PATCH 14/47] Remove debug logging --- src/systems/callbacks.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index 884d1d87f3..a44b5093c8 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -861,8 +861,6 @@ function compile_affect_fn(cb, sys::AbstractODESystem, dvs, ps, kwargs) if isnothing(aff) return nothing else - affspr = compile_affect(aff, cb, sys, dvs, ps; expression = Val{true}, kwargs...) - @show affspr return compile_affect(aff, cb, sys, dvs, ps; expression = Val{false}, kwargs...) end end From 95fa1ee554c0cf449c3c2ac9abdcff85546c80d4 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Mon, 9 Sep 2024 20:06:45 -0700 Subject: [PATCH 15/47] Fix the namespace operation used while namespacing MutatingFunctionalAffects --- src/systems/callbacks.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index a44b5093c8..a7834d918e 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -161,7 +161,7 @@ end function namespace_affect(affect::MutatingFunctionalAffect, s) MutatingFunctionalAffect(func(affect), - renamespace.((s,), observed(affect)), + namespace_expr.(observed(affect), (s,)), observed_syms(affect), renamespace.((s,), modified(affect)), modified_syms(affect), From f57215ab42f5c4c3747744827b8793403c2e6761 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Tue, 10 Sep 2024 12:45:14 -0700 Subject: [PATCH 16/47] Add support for the initializealg argument in SciMLBase callbacks --- src/systems/callbacks.jl | 41 ++++++++++++++++++++++++++++++---------- test/symbolic_events.jl | 6 +++--- 2 files changed, 34 insertions(+), 13 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index a7834d918e..d2657d2ced 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -216,6 +216,11 @@ Affects (i.e. `affect` and `affect_neg`) can be specified as either: + `modified_parameters` is a vector of the parameters that are *modified* by `f!`. Note that a parameter will not appear in `p` if it only appears in `modified_parameters`; it must appear in both `parameters` and `modified_parameters` if it is used in the affect definition. + `ctx` is a user-defined context object passed to `f!` when invoked. This value is aliased for each problem. * A [`MutatingFunctionalAffect`](@ref); refer to its documentation for details. + +Callbacks that impact a DAE are applied, then the DAE is reinitialized using `reinitializealg` (which defaults to `SciMLBase.CheckInit`). +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. """ struct SymbolicContinuousCallback eqs::Vector{Equation} @@ -224,14 +229,16 @@ struct SymbolicContinuousCallback affect::Union{Vector{Equation}, FunctionalAffect, MutatingFunctionalAffect} affect_neg::Union{Vector{Equation}, FunctionalAffect, MutatingFunctionalAffect, Nothing} rootfind::SciMLBase.RootfindOpt + reinitializealg::SciMLBase.DAEInitializationAlgorithm function SymbolicContinuousCallback(; eqs::Vector{Equation}, affect = NULL_AFFECT, affect_neg = affect, rootfind = SciMLBase.LeftRootFind, initialize=NULL_AFFECT, - finalize=NULL_AFFECT) - new(eqs, initialize, finalize, make_affect(affect), make_affect(affect_neg), rootfind) + finalize=NULL_AFFECT, + reinitializealg=SciMLBase.CheckInit()) + new(eqs, initialize, finalize, make_affect(affect), make_affect(affect_neg), rootfind, reinitializealg) end # Default affect to nothing end make_affect(affect) = affect @@ -373,6 +380,10 @@ function finalize_affects(cbs::Vector{SymbolicContinuousCallback}) mapreduce(finalize_affects, vcat, cbs, init = Equation[]) end +reinitialization_alg(cb::SymbolicContinuousCallback) = cb.reinitializealg +reinitialization_algs(cbs::Vector{SymbolicContinuousCallback}) = + mapreduce(reinitialization_alg, vcat, cbs, init = SciMLBase.DAEInitializationAlgorithm[]) + 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(af::MutatingFunctionalAffect, s) = namespace_affect(af, s) @@ -419,11 +430,12 @@ struct SymbolicDiscreteCallback # TODO: Iterative condition::Any affects::Any + reinitializealg::SciMLBase.DAEInitializationAlgorithm - function SymbolicDiscreteCallback(condition, affects = NULL_AFFECT) + function SymbolicDiscreteCallback(condition, affects = NULL_AFFECT, reinitializealg=SciMLBase.CheckInit()) c = scalarize_condition(condition) a = scalarize_affects(affects) - new(c, a) + new(c, a, reinitializealg) end # Default affect to nothing end @@ -481,6 +493,10 @@ function affects(cbs::Vector{SymbolicDiscreteCallback}) reduce(vcat, affects(cb) for cb in cbs; init = []) end +reinitialization_alg(cb::SymbolicDiscreteCallback) = cb.reinitializealg +reinitialization_algs(cbs::Vector{SymbolicDiscreteCallback}) = + mapreduce(reinitialization_alg, vcat, cbs, init = SciMLBase.DAEInitializationAlgorithm[]) + function namespace_callback(cb::SymbolicDiscreteCallback, s)::SymbolicDiscreteCallback af = affects(cb) af = af isa AbstractVector ? namespace_affect.(af, Ref(s)) : namespace_affect(af, s) @@ -776,12 +792,13 @@ function generate_single_rootfinding_callback( return ContinuousCallback( cond, affect_function.affect, affect_function.affect_neg, rootfind = cb.rootfind, initialize = isnothing(affect_function.initialize) ? SciMLBase.INITIALIZE_DEFAULT : (c, u, t, i) -> affect_function.initialize(i), - finalize = isnothing(affect_function.finalize) ? SciMLBase.FINALIZE_DEFAULT : (c, u, t, i) -> affect_function.finalize(i)) + finalize = isnothing(affect_function.finalize) ? SciMLBase.FINALIZE_DEFAULT : (c, u, t, i) -> affect_function.finalize(i), + initializealg = reinitialization_alg(cb)) end function generate_vector_rootfinding_callback( cbs, sys::AbstractODESystem, dvs = unknowns(sys), - ps = parameters(sys); rootfind = SciMLBase.RightRootFind, kwargs...) + ps = parameters(sys); rootfind = SciMLBase.RightRootFind, reinitialization = SciMLBase.CheckInit(), kwargs...) eqs = map(cb -> flatten_equations(cb.eqs), cbs) num_eqs = length.(eqs) # fuse equations to create VectorContinuousCallback @@ -847,7 +864,7 @@ function generate_vector_rootfinding_callback( initialize = handle_optional_setup_fn(map(fn -> fn.initialize, affect_functions), SciMLBase.INITIALIZE_DEFAULT) 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) + cond, affect, affect_neg, length(eqs), rootfind = rootfind, initialize = initialize, finalize = finalize, initializealg = reinitialization) end """ @@ -893,10 +910,14 @@ function generate_rootfinding_callback(cbs, sys::AbstractODESystem, dvs = unknow # group the cbs by what rootfind op they use # groupby would be very useful here, but alas cb_classes = Dict{ - @NamedTuple{rootfind::SciMLBase.RootfindOpt}, Vector{SymbolicContinuousCallback}}() + @NamedTuple{ + rootfind::SciMLBase.RootfindOpt, + reinitialization::SciMLBase.DAEInitializationAlgorithm}, Vector{SymbolicContinuousCallback}}() for cb in cbs push!( - get!(() -> SymbolicContinuousCallback[], cb_classes, (rootfind = cb.rootfind,)), + get!(() -> SymbolicContinuousCallback[], cb_classes, ( + rootfind = cb.rootfind, + reinitialization = reinitialization_alg(cb))), cb) end @@ -904,7 +925,7 @@ function generate_rootfinding_callback(cbs, sys::AbstractODESystem, dvs = unknow 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, kwargs...) + cbs_in_class, sys, dvs, ps; rootfind = equiv_class.rootfind, reinitialization=equiv_class.reinitialization, kwargs...) end if length(compiled_callbacks) == 1 return compiled_callbacks[] diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index b0f230d75e..dd60b99003 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -996,8 +996,8 @@ end @test sort(canonicalize(Discrete(), prob.p)[1]) == [0.0, 1.0, 2.0] sol = solve(prob, Tsit5()) - @test sol[a] == [1.0, -1.0] - @test sol[b] == [2.0, 5.0, 5.0] + @test sol[a] == [-1.0] + @test sol[b] == [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 @@ -1198,5 +1198,5 @@ end ss = structural_simplify(sys) prob = ODEProblem(ss, [theta => 0.0], (0.0, pi)) sol = solve(prob, Tsit5(); dtmax = 0.01) - @test sol[cnt] == 197 # we get 2 pulses per phase cycle (cos 0 crossing) and we go to 100 cycles; we miss a few due to the initial state + @test getp(sol, cnt)(sol) == 197 # we get 2 pulses per phase cycle (cos 0 crossing) and we go to 100 cycles; we miss a few due to the initial state end From d79d49de426ece41ab43df618c7960d88611b6c5 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Tue, 10 Sep 2024 16:23:58 -0700 Subject: [PATCH 17/47] Switch MutatingFunctionalAffect from using ComponentArrays to using NamedTuples for heterotyped operation support. --- src/systems/callbacks.jl | 78 +++++++++++++++--------------- src/systems/diffeqs/odesystem.jl | 10 ++-- test/symbolic_events.jl | 81 ++++++++++++++++++++------------ 3 files changed, 94 insertions(+), 75 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index d2657d2ced..34acb5b2ae 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -77,25 +77,29 @@ end `MutatingFunctionalAffect` is a helper for writing affect functions that will compute observed values and ensure that modified values are correctly written back into the system. The affect function `f` needs to have one of four signatures: -* `f(modified::ComponentArray)` if the function only writes values (unknowns or parameters) to the system, -* `f(modified::ComponentArray, observed::ComponentArray)` if the function also reads observed values from the system, -* `f(modified::ComponentArray, observed::ComponentArray, ctx)` if the function needs the user-defined context, -* `f(modified::ComponentArray, observed::ComponentArray, ctx, integrator)` if the function needs the low-level integrator. +* `f(modified::NamedTuple)::NamedTuple` if the function only writes values (unknowns or parameters) to the system, +* `f(modified::NamedTuple, observed::NamedTuple)::NamedTuple` if the function also reads observed values from the system, +* `f(modified::NamedTuple, observed::NamedTuple, ctx)::NamedTuple` if the function needs the user-defined context, +* `f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple` if the function needs the low-level integrator. These will be checked in reverse order (that is, the four-argument version first, than the 3, etc). -The function `f` will be called with `observed` and `modified` `ComponentArray`s that are derived from their respective `NamedTuple` definitions. -Each `NamedTuple` should map an expression to a symbol; for example if we pass `observed=(; x = a + b)` this will alias the result of executing `a+b` in the system as `x` +The function `f` will be called with `observed` and `modified` `NamedTuple`s that are derived from their respective `NamedTuple` definitions. +Each declaration`NamedTuple` should map an expression to a symbol; for example if we pass `observed=(; x = a + b)` this will alias the result of executing `a+b` in the system as `x` so the value of `a + b` will be accessible as `observed.x` in `f`. `modified` currently restricts symbolic expressions to only bare variables, so only tuples of the form `(; x = y)` or `(; x)` (which aliases `x` as itself) are allowed. -Both `observed` and `modified` will be automatically populated with the current values of their corresponding expressions on function entry. -The values in `modified` will be written back to the system after `f` returns. For example, if we want to update the value of `x` to be the result of `x + y` we could write +The argument NamedTuples (for instance `(;x=y)`) will be populated with the declared values on function entry; if we require `(;x=y)` in `observed` and `y=2`, for example, +then the NamedTuple `(;x=2)` will be passed as `observed` to the affect function `f`. + +The NamedTuple returned from `f` includes the values to be written back to the system after `f` returns. For example, if we want to update the value of `x` to be the result of `x + y` we could write MutatingFunctionalAffect(observed=(; x_plus_y = x + y), modified=(; x)) do m, o - m.x = o.x_plus_y + @set! m.x = o.x_plus_y end -The affect function updates the value at `x` in `modified` to be the result of evaluating `x + y` as passed in the observed values. +Where we use Setfield to copy the tuple `m` with a new value for `x`, then return the modified value of `m`. All values updated by the tuple must have names originally declared in +`modified`; a runtime error will be produced if a value is written that does not appear in `modified`. The user can dynamically decide not to write a value back by not including it +in the returned tuple, in which case the associated field will not be updated. """ @kwdef struct MutatingFunctionalAffect f::Any @@ -983,6 +987,18 @@ function unassignable_variables(sys, expr) 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 compile_user_affect(affect::MutatingFunctionalAffect, cb, sys, dvs, ps; kwargs...) #= Implementation sketch: @@ -1016,7 +1032,6 @@ function compile_user_affect(affect::MutatingFunctionalAffect, cb, sys, dvs, ps; end obs_syms = observed_syms(affect) obs_syms, obs_exprs = check_dups(obs_syms, obs_exprs) - obs_size = size.(obs_exprs) # we will generate a work buffer of a ComponentArray that maps obs_syms to arrays of size obs_size mod_exprs = modified(affect) for mexpr in mod_exprs @@ -1033,8 +1048,6 @@ function compile_user_affect(affect::MutatingFunctionalAffect, cb, sys, dvs, ps; end mod_syms = modified_syms(affect) mod_syms, mod_exprs = check_dups(mod_syms, mod_exprs) - _, mod_og_val_fun = build_explicit_observed_function( - sys, mod_exprs; return_inplace = true) overlapping_syms = intersect(mod_syms, obs_syms) if length(overlapping_syms) > 0 && !affect.skip_checks @@ -1048,31 +1061,20 @@ function compile_user_affect(affect::MutatingFunctionalAffect, cb, sys, dvs, ps; else zeros(sz) end - _, obs_fun = build_explicit_observed_function( + obs_fun = build_explicit_observed_function( sys, reduce(vcat, Symbolics.scalarize.(obs_exprs); init = []); - return_inplace = true) - obs_component_array = ComponentArrays.ComponentArray(NamedTuple{(obs_syms...,)}(mkzero.(obs_size))) + array_type = :tuple) + obs_sym_tuple = (obs_syms...,) # okay so now to generate the stuff to assign it back into the system - # note that we reorder the componentarray to make the views coherent wrt the base array mod_pairs = mod_exprs .=> mod_syms - mod_param_pairs = filter(v -> is_parameter(sys, v[1]), mod_pairs) - mod_unk_pairs = filter(v -> !is_parameter(sys, v[1]), mod_pairs) - _, mod_og_val_fun = build_explicit_observed_function( - sys, reduce(vcat, Symbolics.scalarize.([first.(mod_param_pairs); first.(mod_unk_pairs)]); init = []); - return_inplace = true) - upd_params_fun = setu( - sys, reduce(vcat, Symbolics.scalarize.(first.(mod_param_pairs)); init = [])) - upd_unk_fun = setu( - sys, reduce(vcat, Symbolics.scalarize.(first.(mod_unk_pairs)); init = [])) - - upd_component_array = ComponentArrays.ComponentArray(NamedTuple{([last.(mod_param_pairs); - last.(mod_unk_pairs)]...,)}( - [collect(mkzero(size(e)) for e in first.(mod_param_pairs)); - collect(mkzero(size(e)) for e in first.(mod_unk_pairs))])) - upd_params_view = view(upd_component_array, last.(mod_param_pairs)) - upd_unks_view = view(upd_component_array, last.(mod_unk_pairs)) + mod_names = (mod_syms..., ) + mod_og_val_fun = build_explicit_observed_function( + sys, reduce(vcat, Symbolics.scalarize.(first.(mod_pairs)); init = []); + array_type = :tuple) + 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 @@ -1082,13 +1084,13 @@ function compile_user_affect(affect::MutatingFunctionalAffect, cb, sys, dvs, ps; let user_affect = func(affect), ctx = context(affect) function (integ) # update the to-be-mutated values; this ensures that if you do a no-op then nothing happens - mod_og_val_fun(upd_component_array, integ.u, integ.p..., integ.t) + upd_component_array = NamedTuple{mod_names}(mod_og_val_fun(integ.u, integ.p..., integ.t)) # update the observed values - obs_fun(obs_component_array, integ.u, integ.p..., integ.t) + obs_component_array = NamedTuple{obs_sym_tuple}(obs_fun(integ.u, integ.p..., integ.t)) # let the user do their thing - if applicable(user_affect, upd_component_array, obs_component_array, ctx, integ) + modvals = if applicable(user_affect, upd_component_array, obs_component_array, ctx, integ) user_affect(upd_component_array, obs_component_array, ctx, integ) elseif applicable(user_affect, upd_component_array, obs_component_array, ctx) user_affect(upd_component_array, obs_component_array, ctx) @@ -1102,9 +1104,7 @@ function compile_user_affect(affect::MutatingFunctionalAffect, cb, sys, dvs, ps; end # write the new values back to the integrator - upd_params_fun(integ, upd_params_view) - upd_unk_fun(integ, upd_unks_view) - + _generated_writeback(integ, upd_funs, modvals) for idx in save_idxs SciMLBase.save_discretes!(integ, idx) diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index daa4321ed0..e99911eec6 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -429,6 +429,7 @@ Options not otherwise specified are: * `op = Operator` sets the recursion terminator for the walk done by `vars` to identify the variables that appear in `ts`. See the documentation for `vars` for more detail. * `throw = true` if true, throw an error when generating a function for `ts` that reference variables that do not exist * `drop_expr` is deprecated. +* `array_type`; only used if the output is an array (that is, `!isscalar(ts)`). If `:array`, then it will generate an array, if `:tuple` then it will generate a tuple. """ function build_explicit_observed_function(sys, ts; inputs = nothing, @@ -442,7 +443,8 @@ function build_explicit_observed_function(sys, ts; return_inplace = false, param_only = false, op = Operator, - throw = true) + throw = true, + array_type=:array) if (isscalar = symbolic_type(ts) !== NotSymbolic()) ts = [ts] end @@ -587,12 +589,10 @@ function build_explicit_observed_function(sys, ts; oop_mtkp_wrapper = mtkparams_wrapper end + output_expr = isscalar ? ts[1] : (array_type == :array ? MakeArray(ts, output_type) : MakeTuple(ts)) # Need to keep old method of building the function since it uses `output_type`, # which can't be provided to `build_function` - oop_fn = Func(args, [], - pre(Let(obsexprs, - isscalar ? ts[1] : MakeArray(ts, output_type), - false))) |> array_wrapper[1] |> oop_mtkp_wrapper |> toexpr + oop_fn = Func(args, [], pre(Let(obsexprs, output_expr, false))) |> array_wrapper[1] |> oop_mtkp_wrapper |> toexpr oop_fn = expression ? oop_fn : eval_or_rgf(oop_fn; eval_expression, eval_module) if !isscalar diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index dd60b99003..bc455ec06e 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -8,6 +8,7 @@ using ModelingToolkit: SymbolicContinuousCallback, using StableRNGs import SciMLBase using SymbolicIndexingInterface +using Setfield rng = StableRNG(12345) @variables x(t) = 0 @@ -1010,12 +1011,12 @@ end furnace_off = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_off_threshold], ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o, i, c - x.furnace_on = false + @set! x.furnace_on = false end) furnace_enable = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_on_threshold], ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o, i, c - x.furnace_on = true + @set! x.furnace_on = true end) @named sys = ODESystem( eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) @@ -1027,12 +1028,12 @@ end furnace_off = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_off_threshold], ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o, i - x.furnace_on = false + @set! x.furnace_on = false end) furnace_enable = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_on_threshold], ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o, i - x.furnace_on = true + @set! x.furnace_on = true end) @named sys = ODESystem( eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) @@ -1044,12 +1045,12 @@ end furnace_off = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_off_threshold], ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o - x.furnace_on = false + @set! x.furnace_on = false end) furnace_enable = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_on_threshold], ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o - x.furnace_on = true + @set! x.furnace_on = true end) @named sys = ODESystem( eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) @@ -1061,12 +1062,12 @@ end furnace_off = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_off_threshold], ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x - x.furnace_on = false + @set! x.furnace_on = false end) furnace_enable = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_on_threshold], ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x - x.furnace_on = true + @set! x.furnace_on = true end) @named sys = ODESystem( eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) @@ -1078,14 +1079,14 @@ end furnace_off = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_off_threshold], ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x - x.furnace_on = false + @set! x.furnace_on = false end; initialize = ModelingToolkit.MutatingFunctionalAffect(modified = (; temp)) do x - x.temp = 0.2 + @set! x.temp = 0.2 end) furnace_enable = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_on_threshold], ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o, c, i - x.furnace_on = true + @set! x.furnace_on = true end) @named sys = ODESystem( eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) @@ -1107,7 +1108,7 @@ end [temp ~ furnace_off_threshold], ModelingToolkit.MutatingFunctionalAffect( modified = (; furnace_on), observed = (; furnace_on)) do x, o, c, i - x.furnace_on = false + @set! x.furnace_on = false end) @named sys = ODESystem(eqs, t, [temp], params; continuous_events = [furnace_off]) ss = structural_simplify(sys) @@ -1123,7 +1124,7 @@ end [temp ~ furnace_off_threshold], ModelingToolkit.MutatingFunctionalAffect( modified = (; furnace_on, tempsq), observed = (; furnace_on)) do x, o, c, i - x.furnace_on = false + @set! x.furnace_on = false end) @named sys = ODESystem( eqs, t, [temp, tempsq], params; continuous_events = [furnace_off]) @@ -1136,18 +1137,32 @@ end [temp ~ furnace_off_threshold], ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on), observed = (; furnace_on, not_actually_here)) do x, o, c, i - x.furnace_on = false + @set! x.furnace_on = false end) @named sys = ODESystem( eqs, t, [temp, tempsq], params; continuous_events = [furnace_off]) ss = structural_simplify(sys) @test_throws "refers to missing variable(s)" prob=ODEProblem( ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + + + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on), + observed = (; furnace_on)) do x, o, c, i + return (;fictional2 = false) + end) + @named sys = ODESystem( + eqs, t, [temp, tempsq], params; continuous_events = [furnace_off]) + ss = structural_simplify(sys) + prob=ODEProblem( + ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + @test_throws "Tried to write back to" solve(prob, Tsit5()) end @testset "Quadrature" begin @variables theta(t) omega(t) - params = @parameters qA=0 qB=0 hA=0 hB=0 cnt=0 + params = @parameters qA=0 qB=0 hA=0 hB=0 cnt::Int=0 eqs = [D(theta) ~ omega omega ~ 1.0] function decoder(oldA, oldB, newA, newB) @@ -1167,31 +1182,35 @@ end end qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0], ModelingToolkit.MutatingFunctionalAffect((; qA, hA, hB, cnt), (; qB)) do x, o, i, c - x.hA = x.qA - x.hB = o.qB - x.qA = 1 - x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) + @set! x.hA = x.qA + @set! x.hB = o.qB + @set! x.qA = 1 + @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) + x end, affect_neg = ModelingToolkit.MutatingFunctionalAffect( (; qA, hA, hB, cnt), (; qB)) do x, o, c, i - x.hA = x.qA - x.hB = o.qB - x.qA = 0 - x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) + @set! x.hA = x.qA + @set! x.hB = o.qB + @set! x.qA = 0 + @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) + x end; rootfind = SciMLBase.RightRootFind) qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π / 2) ~ 0], ModelingToolkit.MutatingFunctionalAffect((; qB, hA, hB, cnt), (; qA)) do x, o, i, c - x.hA = o.qA - x.hB = x.qB - x.qB = 1 - x.cnt += decoder(x.hA, x.hB, o.qA, x.qB) + @set! x.hA = o.qA + @set! x.hB = x.qB + @set! x.qB = 1 + @set! x.cnt += decoder(x.hA, x.hB, o.qA, x.qB) + x end, affect_neg = ModelingToolkit.MutatingFunctionalAffect( (; qB, hA, hB, cnt), (; qA)) do x, o, c, i - x.hA = o.qA - x.hB = x.qB - x.qB = 0 - x.cnt += decoder(x.hA, x.hB, o.qA, x.qB) + @set! x.hA = o.qA + @set! x.hB = x.qB + @set! x.qB = 0 + @set! x.cnt += decoder(x.hA, x.hB, o.qA, x.qB) + x end; rootfind = SciMLBase.RightRootFind) @named sys = ODESystem( eqs, t, [theta, omega], params; continuous_events = [qAevt, qBevt]) From c940e5ed27eeed1cafdb3d2a408a245b31950c3e Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Fri, 13 Sep 2024 09:10:20 -0700 Subject: [PATCH 18/47] Fix support for array forms in the NamedTuple version of MutatingFunctionalAffect --- src/systems/callbacks.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index 34acb5b2ae..ad41807ab8 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -1062,7 +1062,7 @@ function compile_user_affect(affect::MutatingFunctionalAffect, cb, sys, dvs, ps; zeros(sz) end obs_fun = build_explicit_observed_function( - sys, reduce(vcat, Symbolics.scalarize.(obs_exprs); init = []); + sys, Symbolics.scalarize.(obs_exprs); array_type = :tuple) obs_sym_tuple = (obs_syms...,) @@ -1070,7 +1070,7 @@ function compile_user_affect(affect::MutatingFunctionalAffect, cb, sys, dvs, ps; mod_pairs = mod_exprs .=> mod_syms mod_names = (mod_syms..., ) mod_og_val_fun = build_explicit_observed_function( - sys, reduce(vcat, Symbolics.scalarize.(first.(mod_pairs)); init = []); + sys, Symbolics.scalarize.(first.(mod_pairs)); array_type = :tuple) upd_funs = NamedTuple{mod_names}((setu.((sys,), first.(mod_pairs))...,)) @@ -1084,7 +1084,8 @@ function compile_user_affect(affect::MutatingFunctionalAffect, cb, sys, dvs, ps; let user_affect = func(affect), ctx = context(affect) function (integ) # update the to-be-mutated values; this ensures that if you do a no-op then nothing happens - upd_component_array = NamedTuple{mod_names}(mod_og_val_fun(integ.u, integ.p..., integ.t)) + modvals = mod_og_val_fun(integ.u, integ.p..., integ.t) + upd_component_array = NamedTuple{mod_names}(modvals) # update the observed values obs_component_array = NamedTuple{obs_sym_tuple}(obs_fun(integ.u, integ.p..., integ.t)) From 2206425cae06f57c2985607c1c18d740bb7d439c Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Mon, 23 Sep 2024 08:18:16 -0700 Subject: [PATCH 19/47] Remove ComponentArrays dep, cleanup handling of skip_checks --- Project.toml | 1 - src/systems/callbacks.jl | 29 +++++++++++++++-------------- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/Project.toml b/Project.toml index 729966fde0..808d68af06 100644 --- a/Project.toml +++ b/Project.toml @@ -9,7 +9,6 @@ ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" -ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index ad41807ab8..cf12b9078b 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -1024,26 +1024,27 @@ function compile_user_affect(affect::MutatingFunctionalAffect, cb, sys, dvs, ps; end obs_exprs = observed(affect) - for oexpr in obs_exprs - invalid_vars = invalid_variables(sys, oexpr) - if length(invalid_vars) > 0 && !affect.skip_checks - error("Observed equation $(oexpr) in affect refers to missing variable(s) $(invalid_vars); the variables may not have been added (e.g. if a component is missing).") + if !affect.skip_checks + for oexpr in obs_exprs + invalid_vars = invalid_variables(sys, oexpr) + if length(invalid_vars) > 0 + error("Observed equation $(oexpr) in affect refers to missing variable(s) $(invalid_vars); the variables may not have been added (e.g. if a component is missing).") + end end end obs_syms = observed_syms(affect) obs_syms, obs_exprs = check_dups(obs_syms, obs_exprs) mod_exprs = modified(affect) - for mexpr in mod_exprs - if affect.skip_checks - continue - end - if !is_variable(sys, mexpr) && parameter_index(sys, mexpr) === nothing && !affect.skip_checks - @warn ("Expression $mexpr cannot be assigned to; currently only unknowns and parameters may be updated by an affect.") - end - invalid_vars = unassignable_variables(sys, mexpr) - if length(invalid_vars) > 0 && !affect.skip_checks - error("Modified equation $(mexpr) in affect refers to missing variable(s) $(invalid_vars); the variables may not have been added (e.g. if a component is missing) or they may have been reduced away.") + if !affect.skip_checks + for mexpr in mod_exprs + if !is_variable(sys, mexpr) && parameter_index(sys, mexpr) === nothing + @warn ("Expression $mexpr cannot be assigned to; currently only unknowns and parameters may be updated by an affect.") + end + invalid_vars = unassignable_variables(sys, mexpr) + if length(invalid_vars) > 0 + error("Modified equation $(mexpr) in affect refers to missing variable(s) $(invalid_vars); the variables may not have been added (e.g. if a component is missing) or they may have been reduced away.") + end end end mod_syms = modified_syms(affect) From 98dcd4ee1ff90a572c273977f780c4cd23f4ce0d Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Mon, 23 Sep 2024 09:59:43 -0700 Subject: [PATCH 20/47] Improve detection of writeback values --- src/systems/callbacks.jl | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index cf12b9078b..23b7ccac2c 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -999,6 +999,18 @@ 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_user_affect(affect::MutatingFunctionalAffect, cb, sys, dvs, ps; kwargs...) #= Implementation sketch: @@ -1038,7 +1050,7 @@ function compile_user_affect(affect::MutatingFunctionalAffect, cb, sys, dvs, ps; mod_exprs = modified(affect) if !affect.skip_checks for mexpr in mod_exprs - if !is_variable(sys, mexpr) && parameter_index(sys, mexpr) === nothing + if !check_assignable(sys, mexpr) @warn ("Expression $mexpr cannot be assigned to; currently only unknowns and parameters may be updated by an affect.") end invalid_vars = unassignable_variables(sys, mexpr) From 3fd4462d31d0c10b77fae3a0c6678d36f9fe0046 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Mon, 23 Sep 2024 10:06:16 -0700 Subject: [PATCH 21/47] Remove ComponentArrays dep --- src/ModelingToolkit.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index f5262a1526..2f57bb1765 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -54,7 +54,6 @@ using Reexport using RecursiveArrayTools import Graphs: SimpleDiGraph, add_edge!, incidence_matrix import BlockArrays: BlockedArray, Block, blocksize, blocksizes -import ComponentArrays using RuntimeGeneratedFunctions using RuntimeGeneratedFunctions: drop_expr From e6ce6ab65fe9e62b5c34f8adc546ceafa91f475f Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Thu, 17 Oct 2024 12:00:44 -0700 Subject: [PATCH 22/47] Rename MutatingFunctionalAffect to ImperativeAffect --- src/systems/callbacks.jl | 74 +++++++++++++++---------------- src/systems/index_cache.jl | 2 +- test/symbolic_events.jl | 90 +++++++++++++++++++------------------- 3 files changed, 83 insertions(+), 83 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index 23b7ccac2c..73f2a7a6cf 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -72,9 +72,9 @@ function namespace_affect(affect::FunctionalAffect, s) end """ - MutatingFunctionalAffect(f::Function; modified::NamedTuple, observed::NamedTuple, ctx) + ImperativeAffect(f::Function; modified::NamedTuple, observed::NamedTuple, ctx) -`MutatingFunctionalAffect` is a helper for writing affect functions that will compute observed values and +`ImperativeAffect` is a helper for writing affect functions that will compute observed values and ensure that modified values are correctly written back into the system. The affect function `f` needs to have one of four signatures: * `f(modified::NamedTuple)::NamedTuple` if the function only writes values (unknowns or parameters) to the system, @@ -93,7 +93,7 @@ then the NamedTuple `(;x=2)` will be passed as `observed` to the affect function The NamedTuple returned from `f` includes the values to be written back to the system after `f` returns. For example, if we want to update the value of `x` to be the result of `x + y` we could write - MutatingFunctionalAffect(observed=(; x_plus_y = x + y), modified=(; x)) do m, o + ImperativeAffect(observed=(; x_plus_y = x + y), modified=(; x)) do m, o @set! m.x = o.x_plus_y end @@ -101,7 +101,7 @@ Where we use Setfield to copy the tuple `m` with a new value for `x`, then retur `modified`; a runtime error will be produced if a value is written that does not appear in `modified`. The user can dynamically decide not to write a value back by not including it in the returned tuple, in which case the associated field will not be updated. """ -@kwdef struct MutatingFunctionalAffect +@kwdef struct ImperativeAffect f::Any obs::Vector obs_syms::Vector{Symbol} @@ -111,50 +111,50 @@ in the returned tuple, in which case the associated field will not be updated. skip_checks::Bool end -function MutatingFunctionalAffect(f::Function; +function ImperativeAffect(f::Function; observed::NamedTuple = NamedTuple{()}(()), modified::NamedTuple = NamedTuple{()}(()), ctx = nothing, skip_checks = false) - MutatingFunctionalAffect(f, + ImperativeAffect(f, collect(values(observed)), collect(keys(observed)), collect(values(modified)), collect(keys(modified)), ctx, skip_checks) end -function MutatingFunctionalAffect(f::Function, modified::NamedTuple; +function ImperativeAffect(f::Function, modified::NamedTuple; observed::NamedTuple = NamedTuple{()}(()), ctx = nothing, skip_checks=false) - MutatingFunctionalAffect(f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) + ImperativeAffect(f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) end -function MutatingFunctionalAffect( +function ImperativeAffect( f::Function, modified::NamedTuple, observed::NamedTuple; ctx = nothing, skip_checks=false) - MutatingFunctionalAffect(f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) + ImperativeAffect(f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) end -function MutatingFunctionalAffect( +function ImperativeAffect( f::Function, modified::NamedTuple, observed::NamedTuple, ctx; skip_checks=false) - MutatingFunctionalAffect(f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) + ImperativeAffect(f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) end -function Base.show(io::IO, mfa::MutatingFunctionalAffect) +function Base.show(io::IO, mfa::ImperativeAffect) obs_vals = join(map((ob,nm) -> "$ob => $nm", mfa.obs, mfa.obs_syms), ", ") mod_vals = join(map((md,nm) -> "$md => $nm", mfa.modified, mfa.mod_syms), ", ") affect = mfa.f - print(io, "MutatingFunctionalAffect(observed: [$obs_vals], modified: [$mod_vals], affect:$affect)") + print(io, "ImperativeAffect(observed: [$obs_vals], modified: [$mod_vals], affect:$affect)") end -func(f::MutatingFunctionalAffect) = f.f -context(a::MutatingFunctionalAffect) = a.ctx -observed(a::MutatingFunctionalAffect) = a.obs -observed_syms(a::MutatingFunctionalAffect) = a.obs_syms -discretes(a::MutatingFunctionalAffect) = filter(ModelingToolkit.isparameter, a.modified) -modified(a::MutatingFunctionalAffect) = a.modified -modified_syms(a::MutatingFunctionalAffect) = a.mod_syms +func(f::ImperativeAffect) = f.f +context(a::ImperativeAffect) = a.ctx +observed(a::ImperativeAffect) = a.obs +observed_syms(a::ImperativeAffect) = a.obs_syms +discretes(a::ImperativeAffect) = filter(ModelingToolkit.isparameter, a.modified) +modified(a::ImperativeAffect) = a.modified +modified_syms(a::ImperativeAffect) = a.mod_syms -function Base.:(==)(a1::MutatingFunctionalAffect, a2::MutatingFunctionalAffect) +function Base.:(==)(a1::ImperativeAffect, a2::ImperativeAffect) isequal(a1.f, a2.f) && isequal(a1.obs, a2.obs) && isequal(a1.modified, a2.modified) && isequal(a1.obs_syms, a2.obs_syms) && isequal(a1.mod_syms, a2.mod_syms) && isequal(a1.ctx, a2.ctx) end -function Base.hash(a::MutatingFunctionalAffect, s::UInt) +function Base.hash(a::ImperativeAffect, s::UInt) s = hash(a.f, s) s = hash(a.obs, s) s = hash(a.obs_syms, s) @@ -163,8 +163,8 @@ function Base.hash(a::MutatingFunctionalAffect, s::UInt) hash(a.ctx, s) end -function namespace_affect(affect::MutatingFunctionalAffect, s) - MutatingFunctionalAffect(func(affect), +function namespace_affect(affect::ImperativeAffect, s) + ImperativeAffect(func(affect), namespace_expr.(observed(affect), (s,)), observed_syms(affect), renamespace.((s,), modified(affect)), @@ -174,7 +174,7 @@ function namespace_affect(affect::MutatingFunctionalAffect, s) end function has_functional_affect(cb) - (affects(cb) isa FunctionalAffect || affects(cb) isa MutatingFunctionalAffect) + (affects(cb) isa FunctionalAffect || affects(cb) isa ImperativeAffect) end #################################### continuous events ##################################### @@ -219,7 +219,7 @@ Affects (i.e. `affect` and `affect_neg`) can be specified as either: + `read_parameters` is a vector of the parameters that are *used* by `f!`. Their indices are passed to `f` in `p` similarly to the indices of `unknowns` passed in `u`. + `modified_parameters` is a vector of the parameters that are *modified* by `f!`. Note that a parameter will not appear in `p` if it only appears in `modified_parameters`; it must appear in both `parameters` and `modified_parameters` if it is used in the affect definition. + `ctx` is a user-defined context object passed to `f!` when invoked. This value is aliased for each problem. -* A [`MutatingFunctionalAffect`](@ref); refer to its documentation for details. +* A [`ImperativeAffect`](@ref); refer to its documentation for details. Callbacks that impact a DAE are applied, then the DAE is reinitialized using `reinitializealg` (which defaults to `SciMLBase.CheckInit`). This reinitialization algorithm ensures that the DAE is satisfied after the callback runs. The default value of `CheckInit` will simply validate @@ -228,10 +228,10 @@ initialization. """ struct SymbolicContinuousCallback eqs::Vector{Equation} - initialize::Union{Vector{Equation}, FunctionalAffect, MutatingFunctionalAffect} - finalize::Union{Vector{Equation}, FunctionalAffect, MutatingFunctionalAffect} - affect::Union{Vector{Equation}, FunctionalAffect, MutatingFunctionalAffect} - affect_neg::Union{Vector{Equation}, FunctionalAffect, MutatingFunctionalAffect, Nothing} + 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 reinitializealg::SciMLBase.DAEInitializationAlgorithm function SymbolicContinuousCallback(; @@ -390,7 +390,7 @@ reinitialization_algs(cbs::Vector{SymbolicContinuousCallback}) = 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(af::MutatingFunctionalAffect, s) = namespace_affect(af, s) +namespace_affects(af::ImperativeAffect, s) = namespace_affect(af, s) namespace_affects(::Nothing, s) = nothing function namespace_callback(cb::SymbolicContinuousCallback, s)::SymbolicContinuousCallback @@ -460,7 +460,7 @@ scalarize_affects(affects) = scalarize(affects) scalarize_affects(affects::Tuple) = FunctionalAffect(affects...) scalarize_affects(affects::NamedTuple) = FunctionalAffect(; affects...) scalarize_affects(affects::FunctionalAffect) = affects -scalarize_affects(affects::MutatingFunctionalAffect) = affects +scalarize_affects(affects::ImperativeAffect) = affects SymbolicDiscreteCallback(p::Pair) = SymbolicDiscreteCallback(p[1], p[2]) SymbolicDiscreteCallback(cb::SymbolicDiscreteCallback) = cb # passthrough @@ -468,7 +468,7 @@ 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 MutatingFunctionalAffect + if db.affects isa FunctionalAffect || db.affects isa ImperativeAffect # TODO println(io, " ", db.affects) else @@ -1011,7 +1011,7 @@ function check_assignable(sys, sym) end end -function compile_user_affect(affect::MutatingFunctionalAffect, cb, sys, dvs, ps; kwargs...) +function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs...) #= Implementation sketch: generate observed function (oop), should save to a component array under obs_syms @@ -1113,7 +1113,7 @@ function compile_user_affect(affect::MutatingFunctionalAffect, cb, sys, dvs, ps; elseif applicable(user_affect, upd_component_array) user_affect(upd_component_array) else - @error "User affect function $user_affect needs to implement one of the supported MutatingFunctionalAffect callback forms; see the MutatingFunctionalAffect docstring for more details" + @error "User affect function $user_affect needs to implement one of the supported ImperativeAffect callback forms; see the ImperativeAffect docstring for more details" user_affect(upd_component_array, obs_component_array, integ, ctx) # this WILL error but it'll give a more sensible message end @@ -1127,7 +1127,7 @@ function compile_user_affect(affect::MutatingFunctionalAffect, cb, sys, dvs, ps; end end -function compile_affect(affect::Union{FunctionalAffect, MutatingFunctionalAffect}, cb, sys, dvs, ps; kwargs...) +function compile_affect(affect::Union{FunctionalAffect, ImperativeAffect}, cb, sys, dvs, ps; kwargs...) compile_user_affect(affect, cb, sys, dvs, ps; kwargs...) end diff --git a/src/systems/index_cache.jl b/src/systems/index_cache.jl index 55d819990b..bd295055fa 100644 --- a/src/systems/index_cache.jl +++ b/src/systems/index_cache.jl @@ -126,7 +126,7 @@ function IndexCache(sys::AbstractSystem) for affect in affs if affect isa Equation is_parameter(sys, affect.lhs) && push!(discs, affect.lhs) - elseif affect isa FunctionalAffect || affect isa MutatingFunctionalAffect + elseif affect isa FunctionalAffect || affect isa ImperativeAffect union!(discs, unwrap.(discretes(affect))) else error("Unhandled affect type $(typeof(affect))") diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index bc455ec06e..b301b65650 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -228,10 +228,10 @@ affect_neg = [x ~ 1] @test e[].affect == affect end -@testset "MutatingFunctionalAffect constructors" begin +@testset "ImperativeAffect constructors" begin fmfa(o, x, i, c) = nothing - m = ModelingToolkit.MutatingFunctionalAffect(fmfa) - @test m isa ModelingToolkit.MutatingFunctionalAffect + m = ModelingToolkit.ImperativeAffect(fmfa) + @test m isa ModelingToolkit.ImperativeAffect @test m.f == fmfa @test m.obs == [] @test m.obs_syms == [] @@ -239,8 +239,8 @@ end @test m.mod_syms == [] @test m.ctx === nothing - m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (;)) - @test m isa ModelingToolkit.MutatingFunctionalAffect + m = ModelingToolkit.ImperativeAffect(fmfa, (;)) + @test m isa ModelingToolkit.ImperativeAffect @test m.f == fmfa @test m.obs == [] @test m.obs_syms == [] @@ -248,8 +248,8 @@ end @test m.mod_syms == [] @test m.ctx === nothing - m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (; x)) - @test m isa ModelingToolkit.MutatingFunctionalAffect + m = ModelingToolkit.ImperativeAffect(fmfa, (; x)) + @test m isa ModelingToolkit.ImperativeAffect @test m.f == fmfa @test isequal(m.obs, []) @test m.obs_syms == [] @@ -257,8 +257,8 @@ end @test m.mod_syms == [:x] @test m.ctx === nothing - m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (; y = x)) - @test m isa ModelingToolkit.MutatingFunctionalAffect + m = ModelingToolkit.ImperativeAffect(fmfa, (; y = x)) + @test m isa ModelingToolkit.ImperativeAffect @test m.f == fmfa @test isequal(m.obs, []) @test m.obs_syms == [] @@ -266,8 +266,8 @@ end @test m.mod_syms == [:y] @test m.ctx === nothing - m = ModelingToolkit.MutatingFunctionalAffect(fmfa; observed = (; y = x)) - @test m isa ModelingToolkit.MutatingFunctionalAffect + m = ModelingToolkit.ImperativeAffect(fmfa; observed = (; y = x)) + @test m isa ModelingToolkit.ImperativeAffect @test m.f == fmfa @test isequal(m.obs, [x]) @test m.obs_syms == [:y] @@ -275,8 +275,8 @@ end @test m.mod_syms == [] @test m.ctx === nothing - m = ModelingToolkit.MutatingFunctionalAffect(fmfa; modified = (; x)) - @test m isa ModelingToolkit.MutatingFunctionalAffect + m = ModelingToolkit.ImperativeAffect(fmfa; modified = (; x)) + @test m isa ModelingToolkit.ImperativeAffect @test m.f == fmfa @test isequal(m.obs, []) @test m.obs_syms == [] @@ -284,8 +284,8 @@ end @test m.mod_syms == [:x] @test m.ctx === nothing - m = ModelingToolkit.MutatingFunctionalAffect(fmfa; modified = (; y = x)) - @test m isa ModelingToolkit.MutatingFunctionalAffect + m = ModelingToolkit.ImperativeAffect(fmfa; modified = (; y = x)) + @test m isa ModelingToolkit.ImperativeAffect @test m.f == fmfa @test isequal(m.obs, []) @test m.obs_syms == [] @@ -293,8 +293,8 @@ end @test m.mod_syms == [:y] @test m.ctx === nothing - m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (; x), (; x)) - @test m isa ModelingToolkit.MutatingFunctionalAffect + m = ModelingToolkit.ImperativeAffect(fmfa, (; x), (; x)) + @test m isa ModelingToolkit.ImperativeAffect @test m.f == fmfa @test isequal(m.obs, [x]) @test m.obs_syms == [:x] @@ -302,8 +302,8 @@ end @test m.mod_syms == [:x] @test m.ctx === nothing - m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (; y = x), (; y = x)) - @test m isa ModelingToolkit.MutatingFunctionalAffect + m = ModelingToolkit.ImperativeAffect(fmfa, (; y = x), (; y = x)) + @test m isa ModelingToolkit.ImperativeAffect @test m.f == fmfa @test isequal(m.obs, [x]) @test m.obs_syms == [:y] @@ -311,9 +311,9 @@ end @test m.mod_syms == [:y] @test m.ctx === nothing - m = ModelingToolkit.MutatingFunctionalAffect( + m = ModelingToolkit.ImperativeAffect( fmfa; modified = (; y = x), observed = (; y = x)) - @test m isa ModelingToolkit.MutatingFunctionalAffect + @test m isa ModelingToolkit.ImperativeAffect @test m.f == fmfa @test isequal(m.obs, [x]) @test m.obs_syms == [:y] @@ -321,9 +321,9 @@ end @test m.mod_syms == [:y] @test m.ctx === nothing - m = ModelingToolkit.MutatingFunctionalAffect( + m = ModelingToolkit.ImperativeAffect( fmfa; modified = (; y = x), observed = (; y = x), ctx = 3) - @test m isa ModelingToolkit.MutatingFunctionalAffect + @test m isa ModelingToolkit.ImperativeAffect @test m.f == fmfa @test isequal(m.obs, [x]) @test m.obs_syms == [:y] @@ -331,8 +331,8 @@ end @test m.mod_syms == [:y] @test m.ctx === 3 - m = ModelingToolkit.MutatingFunctionalAffect(fmfa, (; x), (; x), 3) - @test m isa ModelingToolkit.MutatingFunctionalAffect + m = ModelingToolkit.ImperativeAffect(fmfa, (; x), (; x), 3) + @test m isa ModelingToolkit.ImperativeAffect @test m.f == fmfa @test isequal(m.obs, [x]) @test m.obs_syms == [:x] @@ -1010,12 +1010,12 @@ end furnace_off = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_off_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o, i, c + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, i, c @set! x.furnace_on = false end) furnace_enable = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_on_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o, i, c + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, i, c @set! x.furnace_on = true end) @named sys = ODESystem( @@ -1027,12 +1027,12 @@ end furnace_off = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_off_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o, i + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, i @set! x.furnace_on = false end) furnace_enable = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_on_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o, i + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, i @set! x.furnace_on = true end) @named sys = ODESystem( @@ -1044,12 +1044,12 @@ end furnace_off = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_off_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o @set! x.furnace_on = false end) furnace_enable = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_on_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o @set! x.furnace_on = true end) @named sys = ODESystem( @@ -1061,12 +1061,12 @@ end furnace_off = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_off_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x @set! x.furnace_on = false end) furnace_enable = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_on_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x @set! x.furnace_on = true end) @named sys = ODESystem( @@ -1078,14 +1078,14 @@ end furnace_off = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_off_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x @set! x.furnace_on = false - end; initialize = ModelingToolkit.MutatingFunctionalAffect(modified = (; temp)) do x + end; initialize = ModelingToolkit.ImperativeAffect(modified = (; temp)) do x @set! x.temp = 0.2 end) furnace_enable = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_on_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on)) do x, o, c, i + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, c, i @set! x.furnace_on = true end) @named sys = ODESystem( @@ -1097,7 +1097,7 @@ end @test all(sol[temp][sol.t .!= 0.0] .<= 0.79) && all(sol[temp][sol.t .!= 0.0] .>= 0.2) end -@testset "MutatingFunctionalAffect errors and warnings" begin +@testset "ImperativeAffect errors and warnings" 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 eqs = [ @@ -1106,7 +1106,7 @@ end furnace_off = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_off_threshold], - ModelingToolkit.MutatingFunctionalAffect( + ModelingToolkit.ImperativeAffect( modified = (; furnace_on), observed = (; furnace_on)) do x, o, c, i @set! x.furnace_on = false end) @@ -1122,7 +1122,7 @@ end furnace_off = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_off_threshold], - ModelingToolkit.MutatingFunctionalAffect( + ModelingToolkit.ImperativeAffect( modified = (; furnace_on, tempsq), observed = (; furnace_on)) do x, o, c, i @set! x.furnace_on = false end) @@ -1135,7 +1135,7 @@ end @parameters not_actually_here furnace_off = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_off_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on), + ModelingToolkit.ImperativeAffect(modified = (; furnace_on), observed = (; furnace_on, not_actually_here)) do x, o, c, i @set! x.furnace_on = false end) @@ -1148,7 +1148,7 @@ end furnace_off = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_off_threshold], - ModelingToolkit.MutatingFunctionalAffect(modified = (; furnace_on), + ModelingToolkit.ImperativeAffect(modified = (; furnace_on), observed = (; furnace_on)) do x, o, c, i return (;fictional2 = false) end) @@ -1181,14 +1181,14 @@ end end end qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0], - ModelingToolkit.MutatingFunctionalAffect((; qA, hA, hB, cnt), (; qB)) do x, o, i, c + ModelingToolkit.ImperativeAffect((; qA, hA, hB, cnt), (; qB)) do x, o, i, c @set! x.hA = x.qA @set! x.hB = o.qB @set! x.qA = 1 @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) x end, - affect_neg = ModelingToolkit.MutatingFunctionalAffect( + affect_neg = ModelingToolkit.ImperativeAffect( (; qA, hA, hB, cnt), (; qB)) do x, o, c, i @set! x.hA = x.qA @set! x.hB = o.qB @@ -1197,14 +1197,14 @@ end x end; rootfind = SciMLBase.RightRootFind) qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π / 2) ~ 0], - ModelingToolkit.MutatingFunctionalAffect((; qB, hA, hB, cnt), (; qA)) do x, o, i, c + ModelingToolkit.ImperativeAffect((; qB, hA, hB, cnt), (; qA)) do x, o, i, c @set! x.hA = o.qA @set! x.hB = x.qB @set! x.qB = 1 @set! x.cnt += decoder(x.hA, x.hB, o.qA, x.qB) x end, - affect_neg = ModelingToolkit.MutatingFunctionalAffect( + affect_neg = ModelingToolkit.ImperativeAffect( (; qB, hA, hB, cnt), (; qA)) do x, o, c, i @set! x.hA = o.qA @set! x.hB = x.qB From 54bb95af05f6212375ba3a9e2717147b45246ad9 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Thu, 17 Oct 2024 17:02:04 -0700 Subject: [PATCH 23/47] Fix tests --- src/systems/callbacks.jl | 19 +++++++++++-------- test/symbolic_events.jl | 14 ++++++++++++++ 2 files changed, 25 insertions(+), 8 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index 73f2a7a6cf..96758ccf04 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -781,21 +781,24 @@ function generate_single_rootfinding_callback( 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 else - initfn = SciMLBase.INITIALIZE_DEFAULT + initfn = user_initfun end + return ContinuousCallback( cond, affect_function.affect, affect_function.affect_neg, rootfind = cb.rootfind, - initialize = isnothing(affect_function.initialize) ? SciMLBase.INITIALIZE_DEFAULT : (c, u, t, i) -> affect_function.initialize(i), + initialize = initfn, finalize = isnothing(affect_function.finalize) ? SciMLBase.FINALIZE_DEFAULT : (c, u, t, i) -> affect_function.finalize(i), initializealg = reinitialization_alg(cb)) end @@ -878,8 +881,8 @@ function compile_affect_fn(cb, sys::AbstractODESystem, 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) - if isnothing(aff) + 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...) @@ -890,8 +893,8 @@ function compile_affect_fn(cb, sys::AbstractODESystem, dvs, ps, kwargs) else affect_neg = compile_optional_affect(eq_neg_aff) end - initialize = compile_optional_affect(initialize_affects(cb)) - finalize = compile_optional_affect(finalize_affects(cb)) + initialize = compile_optional_affect(initialize_affects(cb), NULL_AFFECT) + finalize = compile_optional_affect(finalize_affects(cb), NULL_AFFECT) (affect = affect, affect_neg = affect_neg, initialize = initialize, finalize = finalize) end @@ -1097,11 +1100,11 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs. let user_affect = func(affect), ctx = context(affect) 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) + modvals = mod_og_val_fun(integ.u, integ.p, integ.t) upd_component_array = NamedTuple{mod_names}(modvals) # update the observed values - obs_component_array = NamedTuple{obs_sym_tuple}(obs_fun(integ.u, integ.p..., integ.t)) + obs_component_array = NamedTuple{obs_sym_tuple}(obs_fun(integ.u, integ.p, integ.t)) # let the user do their thing modvals = if applicable(user_affect, upd_component_array, obs_component_array, ctx, integ) diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index b301b65650..c2c26aae7f 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -1219,3 +1219,17 @@ end sol = solve(prob, Tsit5(); dtmax = 0.01) @test getp(sol, cnt)(sol) == 197 # we get 2 pulses per phase cycle (cos 0 crossing) and we go to 100 cycles; we miss a few due to the initial state end + + + +import RuntimeGeneratedFunctions +function (f::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{argnames, cache_tag, context_tag, id})(args::Vararg{Any, N}) where {N, argnames, cache_tag, context_tag, id} + try + RuntimeGeneratedFunctions.generated_callfunc(f, args...) + catch e + @error "Caught error in RuntimeGeneratedFunction; source code follows" + func_expr = Expr(:->, Expr(:tuple, argnames...), RuntimeGeneratedFunctions._lookup_body(cache_tag, id)) + @show func_expr + rethrow(e) + end +end From 065490971ad5bb91718f0cd2a61463060b962794 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Tue, 22 Oct 2024 18:36:18 -0700 Subject: [PATCH 24/47] Document ImperativeEffect and the SymbolicContinousCallback changes --- docs/Project.toml | 1 + docs/src/basics/Events.md | 198 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 199 insertions(+) diff --git a/docs/Project.toml b/docs/Project.toml index 078df7d696..15f1a6a7f0 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -14,6 +14,7 @@ OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226" +Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index f425fdce5b..9d3ba30780 100644 --- a/docs/src/basics/Events.md +++ b/docs/src/basics/Events.md @@ -378,3 +378,201 @@ sol.ps[c] # sol[c] will error, since `c` is not a timeseries value ``` It can be seen that the timeseries for `c` is not saved. + + +## [(Experimental) Imperative affects](@id imp_affects) +The `ImperativeAffect` can be used as an alternative to the aforementioned functional affect form. Note +that `ImperativeAffect` is still experimental; to emphasize this, we do not export it and it should be +included as `ModelingToolkit.ImperativeAffect`. It abstracts over how values are written back to the +system, simplifying the definitions and (in the future) allowing assignments back to observed values +by solving the nonlinear reinitialization problem afterwards. + +We will use two examples to describe `ImperativeAffect`: a simple heater and a quadrature encoder. +These examples will also demonstrate advanced usage of `ModelingToolkit.SymbolicContinousCallback`, +the low-level interface that the aforementioned tuple form converts into and allows control over the +exact SciMLCallbacks event that is generated for a continous event. + +### [Heater](@id heater_events) +Bang-bang control of a heater connected to a leaky plant requires hysteresis in order to prevent control oscillation. + +```@example events +@variables temp(t) +params = @parameters furnace_on_threshold=0.5 furnace_off_threshold=0.7 furnace_power=1.0 leakage=0.1 furnace_on(t)::Bool=false +eqs = [ + D(temp) ~ furnace_on * furnace_power - temp^2 * leakage +] +``` +Our plant is simple. We have a heater that's turned on and off by the clocked parameter `furnace_on` +which adds `furnace_power` forcing to the system when enabled. We then leak heat porportional to `leakage` +as a function of the square of the current temperature. + +We need a controller with hysteresis to conol the plant. We wish the furnace to turn on when the temperature +is below `furnace_on_threshold` and off when above `furnace_off_threshold`, while maintaining its current state +in between. To do this, we create two continous callbacks: +```@example events +using Setfield +furnace_disable = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, i, c + @set! x.furnace_on = false + end) +furnace_enable = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_on_threshold], + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, i, c + @set! x.furnace_on = true + end) +``` +We're using the explicit form of `SymbolicContinuousCallback` here, though +so far we aren't using anything that's not possible with the implicit interface. +You can also write +```julia +[temp ~ furnace_off_threshold] => ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, i, c + @set! x.furnace_on = false +end +``` +and it would work the same. + +The `ImperativeAffect` is the larger change in this example. `ImperativeAffect` has the constructor signature +```julia + ImperativeAffect(f::Function; modified::NamedTuple, observed::NamedTuple, ctx) +``` +that accepts the function to call, a named tuple of both the names of and symbolic values representing +values in the system to be modified, a named tuple of the values that are merely observed (that is, used from +the system but not modified), and a context that's passed to the affect function. + +In our example, each event merely changes whether the furnace is on or off. Accordingly, we pass a `modified` tuple +`(; furnace_on)` (creating a `NamedTuple` equivalent to `(furnace_on = furnace_on)`). `ImperativeAffect` will then +evaluate this before calling our function to fill out all of the numerical values, then apply them back to the system +once our affect function returns. Furthermore, it will check that it is possible to do this assignment. + +The function given to `ImperativeAffect` needs to have one of four signatures, checked in this order: +* `f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple` if the function needs the low-level integrator, +* `f(modified::NamedTuple, observed::NamedTuple, ctx)::NamedTuple` if the function needs the user-defined context, +* `f(modified::NamedTuple, observed::NamedTuple)::NamedTuple` if the function also reads observed values from the system, +* `f(modified::NamedTuple)::NamedTuple` if the function only writes values (unknowns or parameters) to the system. +The `do` block in the example implicitly constructs said function inline. For exposition, we use the full version (e.g. `x, o, i, c`) but this could be simplified to merely `x`. + +The function `f` will be called with `observed` and `modified` `NamedTuple`s that are derived from their respective `NamedTuple` definitions. +In our example, if `furnace_on` is `false`, then the value of the `x` that's passed in as `modified` will be `(furnace_on = false)`. +The modified values should be passed out in the same format: to set `furnace_on` to `true` we need to return a tuple `(furnace_on = true)`. +We use Setfield to do this in the example, recreating the result tuple before returning it. + +Accordingly, we can now interpret the `ImperativeAffect` definitions to mean that when `temp = furnace_off_threshold` we +will write `furnace_on = false` back to the system, and when `temp = furnace_on_threshold` we will write `furnace_on = true` back +to the system. + +```@example events +@named sys = ODESystem( + eqs, t, [temp], params; continuous_events = [furnace_disable, furnace_enable]) +ss = structural_simplify(sys) +prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 10.0)) +sol = solve(prob, Tsit5()) +plot(sol) +hline!([sol.ps[furnace_off_threshold], sol.ps[furnace_on_threshold]], l = (:black, 1), primary = false) +``` + +Here we see exactly the desired hysteresis. The heater starts on until the temperature hits +`furnace_off_threshold`. The temperature then bleeds away until `furnace_on_threshold` at which +point the furnace turns on again until `furnace_off_threshold` and so on and so forth. The controller +is effectively regulating the temperature of the plant. + +### [Quadrature Encoder](@id quadrature) +For a more complex application we'll look at modeling a quadrature encoder attached to a shaft spinning at a constant speed. +Traditionally, a quadrature encoder is built out of a code wheel that interrupts the sensors at constant intervals and two sensors slightly out of phase with one another. +A state machine can take the pattern of pulses produced by the two sensors and determine the number of steps that the shaft has spun. The state machine takes the new value +from each sensor and the old values and decodes them into the direction that the wheel has spun in this step. + +```@example events + @variables theta(t) omega(t) + params = @parameters qA=0 qB=0 hA=0 hB=0 cnt::Int=0 + eqs = [D(theta) ~ omega + omega ~ 1.0] +``` +Our continous-time system is extremely simple. We have two states, `theta` for the angle of the shaft +and `omega` for the rate at which it's spinning. We then have parameters for the state machine `qA, qB, hA, hB` +and a step count `cnt`. + +We'll then implement the decoder as a simple Julia function. +```@example events + function decoder(oldA, oldB, newA, newB) + state = (oldA, oldB, newA, newB) + if state == (0, 0, 1, 0) || state == (1, 0, 1, 1) || state == (1, 1, 0, 1) || + state == (0, 1, 0, 0) + return 1 + elseif state == (0, 0, 0, 1) || state == (0, 1, 1, 1) || state == (1, 1, 1, 0) || + state == (1, 0, 0, 0) + return -1 + elseif state == (0, 0, 0, 0) || state == (0, 1, 0, 1) || state == (1, 0, 1, 0) || + state == (1, 1, 1, 1) + return 0 + else + return 0 # err is interpreted as no movement + end + end +``` +Based on the current and old state, this function will return 1 if the wheel spun in the positive direction, +-1 if in the negative, and 0 otherwise. + +The encoder state advances when the occlusion begins or ends. We model the +code wheel as simply detecting when `cos(100*theta)` is 0; if we're at a positive +edge of the 0 crossing, then we interpret that as occlusion (so the discrete `qA` goes to 1). Otherwise, if `cos` is +going negative, we interpret that as lack of occlusion (so the discrete goes to 0). The decoder function is +then invoked to update the count with this new information. + +We can implement this in one of two ways: using edge sign detection or right root finding. For exposition, we +will implement each sensor differently. + +For sensor A, we're using the edge detction method. By providing a different affect to `SymbolicContinuousCallback`'s +`affect_neg` argument, we can specify different behaviour for the negative crossing vs. the positive crossing of the root. +In our encoder, we interpret this as occlusion or nonocclusion of the sensor, update the internal state, and tick the decoder. +```@example events + qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0], + ModelingToolkit.ImperativeAffect((; qA, hA, hB, cnt), (; qB)) do x, o, i, c + @set! x.hA = x.qA + @set! x.hB = o.qB + @set! x.qA = 1 + @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) + x + end, + affect_neg = ModelingToolkit.ImperativeAffect( + (; qA, hA, hB, cnt), (; qB)) do x, o, c, i + @set! x.hA = x.qA + @set! x.hB = o.qB + @set! x.qA = 0 + @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) + x + end) +``` + +The other way we can implement a sensor is by changing the root find. +Normally, we use left root finding; the affect will be invoked instantaneously before +the root is crossed. This makes it trickier to figure out what the new state is. +Instead, we can use right root finding: + +```@example events + qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π / 2) ~ 0], + ModelingToolkit.ImperativeAffect((; qB, hA, hB, cnt), (; qA, theta)) do x, o, i, c + @set! x.hA = o.qA + @set! x.hB = x.qB + @set! x.qB = clamp(sign(cos(100 * o.theta - π / 2)), 0.0, 1.0) + @set! x.cnt += decoder(x.hA, x.hB, o.qA, x.qB) + x + end; rootfind = SciMLBase.RightRootFind) +``` +Here, sensor B is located `π / 2` behind sensor A in angular space, so we're adjusting our +trigger function accordingly. We here ask for right root finding on the callback, so we know +that the value of said function will have the "new" sign rather than the old one. Thus, we can +determine the new state of the sensor from the sign of the indicator function evaluated at the +affect activation point, with -1 mapped to 0. + +We can now simulate the encoder. +```@example events + @named sys = ODESystem( + eqs, t, [theta, omega], params; continuous_events = [qAevt, qBevt]) + ss = structural_simplify(sys) + prob = ODEProblem(ss, [theta => 0.0], (0.0, pi)) + sol = solve(prob, Tsit5(); dtmax = 0.01) + sol.ps[cnt] +``` +`cos(100*theta)` will have 200 crossings in the half rotation we've gone through, so the encoder would notionally count 200 steps. +Our encoder counts 198 steps (it loses one step to initialization and one step due to the final state falling squarely on an edge). \ No newline at end of file From aecd59bfea213d5e0650057135afe7eb6262b988 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Tue, 22 Oct 2024 18:40:34 -0700 Subject: [PATCH 25/47] Formatter --- src/systems/callbacks.jl | 164 +++++++++++++++++++++++---------------- test/symbolic_events.jl | 20 ++--- 2 files changed, 107 insertions(+), 77 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index 96758ccf04..e7198817b1 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -116,29 +116,33 @@ function ImperativeAffect(f::Function; modified::NamedTuple = NamedTuple{()}(()), ctx = nothing, skip_checks = false) - ImperativeAffect(f, + ImperativeAffect(f, collect(values(observed)), collect(keys(observed)), - collect(values(modified)), collect(keys(modified)), + collect(values(modified)), collect(keys(modified)), ctx, skip_checks) end function ImperativeAffect(f::Function, modified::NamedTuple; - observed::NamedTuple = NamedTuple{()}(()), ctx = nothing, skip_checks=false) - ImperativeAffect(f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) + observed::NamedTuple = NamedTuple{()}(()), ctx = nothing, skip_checks = false) + ImperativeAffect( + f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) end function ImperativeAffect( - f::Function, modified::NamedTuple, observed::NamedTuple; ctx = nothing, skip_checks=false) - ImperativeAffect(f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) + f::Function, modified::NamedTuple, observed::NamedTuple; ctx = nothing, skip_checks = false) + ImperativeAffect( + f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) end function ImperativeAffect( - f::Function, modified::NamedTuple, observed::NamedTuple, ctx; skip_checks=false) - ImperativeAffect(f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) + f::Function, modified::NamedTuple, observed::NamedTuple, ctx; skip_checks = false) + ImperativeAffect( + f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) end -function Base.show(io::IO, mfa::ImperativeAffect) - obs_vals = join(map((ob,nm) -> "$ob => $nm", mfa.obs, mfa.obs_syms), ", ") - mod_vals = join(map((md,nm) -> "$md => $nm", mfa.modified, mfa.mod_syms), ", ") +function Base.show(io::IO, mfa::ImperativeAffect) + obs_vals = join(map((ob, nm) -> "$ob => $nm", mfa.obs, mfa.obs_syms), ", ") + mod_vals = join(map((md, nm) -> "$md => $nm", mfa.modified, mfa.mod_syms), ", ") affect = mfa.f - print(io, "ImperativeAffect(observed: [$obs_vals], modified: [$mod_vals], affect:$affect)") + print(io, + "ImperativeAffect(observed: [$obs_vals], modified: [$mod_vals], affect:$affect)") end func(f::ImperativeAffect) = f.f context(a::ImperativeAffect) = a.ctx @@ -234,15 +238,16 @@ struct SymbolicContinuousCallback affect_neg::Union{Vector{Equation}, FunctionalAffect, ImperativeAffect, Nothing} rootfind::SciMLBase.RootfindOpt reinitializealg::SciMLBase.DAEInitializationAlgorithm - function SymbolicContinuousCallback(; - eqs::Vector{Equation}, - affect = NULL_AFFECT, - affect_neg = affect, - rootfind = SciMLBase.LeftRootFind, - initialize=NULL_AFFECT, - finalize=NULL_AFFECT, - reinitializealg=SciMLBase.CheckInit()) - new(eqs, initialize, finalize, make_affect(affect), make_affect(affect_neg), rootfind, reinitializealg) + function SymbolicContinuousCallback(; + eqs::Vector{Equation}, + affect = NULL_AFFECT, + affect_neg = affect, + rootfind = SciMLBase.LeftRootFind, + initialize = NULL_AFFECT, + finalize = NULL_AFFECT, + reinitializealg = SciMLBase.CheckInit()) + new(eqs, initialize, finalize, make_affect(affect), + make_affect(affect_neg), rootfind, reinitializealg) end # Default affect to nothing end make_affect(affect) = affect @@ -250,8 +255,8 @@ 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.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) end Base.isempty(cb::SymbolicContinuousCallback) = isempty(cb.eqs) @@ -266,10 +271,9 @@ function Base.hash(cb::SymbolicContinuousCallback, s::UInt) hash(cb.rootfind, s) end - function Base.show(io::IO, cb::SymbolicContinuousCallback) indent = get(io, :indent, 0) - iio = IOContext(io, :indent => indent+1) + iio = IOContext(io, :indent => indent + 1) print(io, "SymbolicContinuousCallback(") print(iio, "Equations:") show(iio, equations(cb)) @@ -298,7 +302,7 @@ end function Base.show(io::IO, mime::MIME"text/plain", cb::SymbolicContinuousCallback) indent = get(io, :indent, 0) - iio = IOContext(io, :indent => indent+1) + iio = IOContext(io, :indent => indent + 1) println(io, "SymbolicContinuousCallback:") println(iio, "Equations:") show(iio, mime, equations(cb)) @@ -338,14 +342,18 @@ 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; - affect_neg = affect, rootfind = SciMLBase.LeftRootFind, initialize = NULL_AFFECT, finalize = NULL_AFFECT) + affect_neg = affect, rootfind = SciMLBase.LeftRootFind, + initialize = NULL_AFFECT, finalize = NULL_AFFECT) SymbolicContinuousCallback( - eqs = [eqs], affect = affect, affect_neg = affect_neg, rootfind = rootfind, initialize=initialize, finalize=finalize) + eqs = [eqs], affect = affect, affect_neg = affect_neg, rootfind = rootfind, + initialize = initialize, finalize = finalize) end function SymbolicContinuousCallback(eqs::Vector{Equation}, affect = NULL_AFFECT; - affect_neg = affect, rootfind = SciMLBase.LeftRootFind, initialize = NULL_AFFECT, finalize = NULL_AFFECT) + affect_neg = affect, rootfind = SciMLBase.LeftRootFind, + initialize = NULL_AFFECT, finalize = NULL_AFFECT) SymbolicContinuousCallback( - eqs = eqs, affect = affect, affect_neg = affect_neg, rootfind = rootfind, initialize=initialize, finalize=finalize) + eqs = eqs, affect = affect, affect_neg = affect_neg, rootfind = rootfind, + initialize = initialize, finalize = finalize) end SymbolicContinuousCallbacks(cb::SymbolicContinuousCallback) = [cb] @@ -385,8 +393,10 @@ function finalize_affects(cbs::Vector{SymbolicContinuousCallback}) end reinitialization_alg(cb::SymbolicContinuousCallback) = cb.reinitializealg -reinitialization_algs(cbs::Vector{SymbolicContinuousCallback}) = - mapreduce(reinitialization_alg, vcat, cbs, init = SciMLBase.DAEInitializationAlgorithm[]) +function reinitialization_algs(cbs::Vector{SymbolicContinuousCallback}) + mapreduce( + reinitialization_alg, vcat, cbs, init = SciMLBase.DAEInitializationAlgorithm[]) +end namespace_affects(af::Vector, s) = Equation[namespace_affect(a, s) for a in af] namespace_affects(af::FunctionalAffect, s) = namespace_affect(af, s) @@ -436,7 +446,8 @@ struct SymbolicDiscreteCallback affects::Any reinitializealg::SciMLBase.DAEInitializationAlgorithm - function SymbolicDiscreteCallback(condition, affects = NULL_AFFECT, reinitializealg=SciMLBase.CheckInit()) + function SymbolicDiscreteCallback( + condition, affects = NULL_AFFECT, reinitializealg = SciMLBase.CheckInit()) c = scalarize_condition(condition) a = scalarize_affects(affects) new(c, a, reinitializealg) @@ -498,8 +509,10 @@ function affects(cbs::Vector{SymbolicDiscreteCallback}) end reinitialization_alg(cb::SymbolicDiscreteCallback) = cb.reinitializealg -reinitialization_algs(cbs::Vector{SymbolicDiscreteCallback}) = - mapreduce(reinitialization_alg, vcat, cbs, init = SciMLBase.DAEInitializationAlgorithm[]) +function reinitialization_algs(cbs::Vector{SymbolicDiscreteCallback}) + mapreduce( + reinitialization_alg, vcat, cbs, init = SciMLBase.DAEInitializationAlgorithm[]) +end function namespace_callback(cb::SymbolicDiscreteCallback, s)::SymbolicDiscreteCallback af = affects(cb) @@ -781,7 +794,8 @@ function generate_single_rootfinding_callback( end end - user_initfun = isnothing(affect_function.initialize) ? SciMLBase.INITIALIZE_DEFAULT : (c, u, t, i) -> affect_function.initialize(i) + 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 @@ -795,17 +809,19 @@ function generate_single_rootfinding_callback( else initfn = user_initfun end - + 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), + 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)) end function generate_vector_rootfinding_callback( cbs, sys::AbstractODESystem, dvs = unknowns(sys), - ps = parameters(sys); rootfind = SciMLBase.RightRootFind, reinitialization = SciMLBase.CheckInit(), kwargs...) + ps = parameters(sys); rootfind = SciMLBase.RightRootFind, + reinitialization = SciMLBase.CheckInit(), kwargs...) eqs = map(cb -> flatten_equations(cb.eqs), cbs) num_eqs = length.(eqs) # fuse equations to create VectorContinuousCallback @@ -821,11 +837,12 @@ function generate_vector_rootfinding_callback( sys, rhss, dvs, ps; expression = Val{false}, kwargs...) affect_functions = @NamedTuple{ - affect::Function, - affect_neg::Union{Function, Nothing}, - initialize::Union{Function, Nothing}, + 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] + 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 @@ -861,17 +878,20 @@ function generate_vector_rootfinding_callback( if isnothing(func) continue else - func(integ) + func(integ) end end end end end end - initialize = handle_optional_setup_fn(map(fn -> fn.initialize, affect_functions), SciMLBase.INITIALIZE_DEFAULT) - finalize = handle_optional_setup_fn(map(fn -> fn.finalize, affect_functions), SciMLBase.FINALIZE_DEFAULT) + initialize = handle_optional_setup_fn( + map(fn -> fn.initialize, affect_functions), SciMLBase.INITIALIZE_DEFAULT) + 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) + cond, affect, affect_neg, length(eqs), rootfind = rootfind, initialize = initialize, + finalize = finalize, initializealg = reinitialization) end """ @@ -881,8 +901,8 @@ function compile_affect_fn(cb, sys::AbstractODESystem, 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 + 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...) @@ -918,13 +938,14 @@ function generate_rootfinding_callback(cbs, sys::AbstractODESystem, dvs = unknow # groupby would be very useful here, but alas cb_classes = Dict{ @NamedTuple{ - rootfind::SciMLBase.RootfindOpt, + rootfind::SciMLBase.RootfindOpt, reinitialization::SciMLBase.DAEInitializationAlgorithm}, Vector{SymbolicContinuousCallback}}() for cb in cbs push!( - get!(() -> SymbolicContinuousCallback[], cb_classes, ( - rootfind = cb.rootfind, - reinitialization = reinitialization_alg(cb))), + get!(() -> SymbolicContinuousCallback[], cb_classes, + ( + rootfind = cb.rootfind, + reinitialization = reinitialization_alg(cb))), cb) end @@ -932,7 +953,8 @@ function generate_rootfinding_callback(cbs, sys::AbstractODESystem, dvs = unknow 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...) + cbs_in_class, sys, dvs, ps; rootfind = equiv_class.rootfind, + reinitialization = equiv_class.reinitialization, kwargs...) end if length(compiled_callbacks) == 1 return compiled_callbacks[] @@ -984,29 +1006,34 @@ 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))); init=[]) + assignable_syms = reduce( + vcat, Symbolics.scalarize.(vcat(unknowns(sys), parameters(sys))); 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} +@generated function _generated_writeback(integ, setters::NamedTuple{NS1, <:Tuple}, + values::NamedTuple{NS2, <:Tuple}) where {NS1, NS2} setter_exprs = [] - for name in NS2 + 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) + 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)) + 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 @@ -1084,13 +1111,13 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs. # okay so now to generate the stuff to assign it back into the system mod_pairs = mod_exprs .=> mod_syms - mod_names = (mod_syms..., ) + mod_names = (mod_syms...,) mod_og_val_fun = build_explicit_observed_function( sys, Symbolics.scalarize.(first.(mod_pairs)); array_type = :tuple) 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 @@ -1104,10 +1131,12 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs. upd_component_array = NamedTuple{mod_names}(modvals) # update the observed values - obs_component_array = NamedTuple{obs_sym_tuple}(obs_fun(integ.u, integ.p, integ.t)) + obs_component_array = NamedTuple{obs_sym_tuple}(obs_fun( + integ.u, integ.p, integ.t)) # let the user do their thing - modvals = if applicable(user_affect, upd_component_array, obs_component_array, ctx, integ) + modvals = if applicable( + user_affect, upd_component_array, obs_component_array, ctx, integ) user_affect(upd_component_array, obs_component_array, ctx, integ) elseif applicable(user_affect, upd_component_array, obs_component_array, ctx) user_affect(upd_component_array, obs_component_array, ctx) @@ -1122,7 +1151,7 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs. # write the new values back to the integrator _generated_writeback(integ, upd_funs, modvals) - + for idx in save_idxs SciMLBase.save_discretes!(integ, idx) end @@ -1130,7 +1159,8 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs. end end -function compile_affect(affect::Union{FunctionalAffect, ImperativeAffect}, cb, sys, dvs, ps; kwargs...) +function compile_affect( + affect::Union{FunctionalAffect, ImperativeAffect}, cb, sys, dvs, ps; kwargs...) compile_user_affect(affect, cb, sys, dvs, ps; kwargs...) end diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index c2c26aae7f..f4f97fb2b9 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -1001,7 +1001,7 @@ end @test sol[b] == [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 +@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 eqs = [ @@ -1080,7 +1080,7 @@ end [temp ~ furnace_off_threshold], ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x @set! x.furnace_on = false - end; initialize = ModelingToolkit.ImperativeAffect(modified = (; temp)) do x + end; initialize = ModelingToolkit.ImperativeAffect(modified = (; temp)) do x @set! x.temp = 0.2 end) furnace_enable = ModelingToolkit.SymbolicContinuousCallback( @@ -1145,17 +1145,16 @@ end @test_throws "refers to missing variable(s)" prob=ODEProblem( ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) - furnace_off = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_off_threshold], ModelingToolkit.ImperativeAffect(modified = (; furnace_on), observed = (; furnace_on)) do x, o, c, i - return (;fictional2 = false) + return (; fictional2 = false) end) @named sys = ODESystem( eqs, t, [temp, tempsq], params; continuous_events = [furnace_off]) ss = structural_simplify(sys) - prob=ODEProblem( + prob = ODEProblem( ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) @test_throws "Tried to write back to" solve(prob, Tsit5()) end @@ -1220,15 +1219,16 @@ end @test getp(sol, cnt)(sol) == 197 # we get 2 pulses per phase cycle (cos 0 crossing) and we go to 100 cycles; we miss a few due to the initial state end - - import RuntimeGeneratedFunctions -function (f::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{argnames, cache_tag, context_tag, id})(args::Vararg{Any, N}) where {N, argnames, cache_tag, context_tag, id} +function (f::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{ + argnames, cache_tag, context_tag, + id})(args::Vararg{Any, N}) where {N, argnames, cache_tag, context_tag, id} try RuntimeGeneratedFunctions.generated_callfunc(f, args...) - catch e + catch e @error "Caught error in RuntimeGeneratedFunction; source code follows" - func_expr = Expr(:->, Expr(:tuple, argnames...), RuntimeGeneratedFunctions._lookup_body(cache_tag, id)) + func_expr = Expr(:->, Expr(:tuple, argnames...), + RuntimeGeneratedFunctions._lookup_body(cache_tag, id)) @show func_expr rethrow(e) end From 89954e46437a7a6fd3536d699d2e3d1846020a48 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Tue, 22 Oct 2024 18:43:31 -0700 Subject: [PATCH 26/47] Formatter (2) --- docs/src/basics/Events.md | 161 +++++++++++++++++-------------- src/systems/diffeqs/odesystem.jl | 8 +- 2 files changed, 95 insertions(+), 74 deletions(-) diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index 9d3ba30780..4a59149308 100644 --- a/docs/src/basics/Events.md +++ b/docs/src/basics/Events.md @@ -379,36 +379,39 @@ sol.ps[c] # sol[c] will error, since `c` is not a timeseries value It can be seen that the timeseries for `c` is not saved. - ## [(Experimental) Imperative affects](@id imp_affects) + The `ImperativeAffect` can be used as an alternative to the aforementioned functional affect form. Note that `ImperativeAffect` is still experimental; to emphasize this, we do not export it and it should be -included as `ModelingToolkit.ImperativeAffect`. It abstracts over how values are written back to the +included as `ModelingToolkit.ImperativeAffect`. It abstracts over how values are written back to the system, simplifying the definitions and (in the future) allowing assignments back to observed values by solving the nonlinear reinitialization problem afterwards. -We will use two examples to describe `ImperativeAffect`: a simple heater and a quadrature encoder. +We will use two examples to describe `ImperativeAffect`: a simple heater and a quadrature encoder. These examples will also demonstrate advanced usage of `ModelingToolkit.SymbolicContinousCallback`, the low-level interface that the aforementioned tuple form converts into and allows control over the exact SciMLCallbacks event that is generated for a continous event. ### [Heater](@id heater_events) + Bang-bang control of a heater connected to a leaky plant requires hysteresis in order to prevent control oscillation. -```@example events +```@example events @variables temp(t) params = @parameters furnace_on_threshold=0.5 furnace_off_threshold=0.7 furnace_power=1.0 leakage=0.1 furnace_on(t)::Bool=false eqs = [ D(temp) ~ furnace_on * furnace_power - temp^2 * leakage ] ``` + Our plant is simple. We have a heater that's turned on and off by the clocked parameter `furnace_on` which adds `furnace_power` forcing to the system when enabled. We then leak heat porportional to `leakage` -as a function of the square of the current temperature. +as a function of the square of the current temperature. We need a controller with hysteresis to conol the plant. We wish the furnace to turn on when the temperature is below `furnace_on_threshold` and off when above `furnace_off_threshold`, while maintaining its current state in between. To do this, we create two continous callbacks: + ```@example events using Setfield furnace_disable = ModelingToolkit.SymbolicContinuousCallback( @@ -422,42 +425,49 @@ furnace_enable = ModelingToolkit.SymbolicContinuousCallback( @set! x.furnace_on = true end) ``` + We're using the explicit form of `SymbolicContinuousCallback` here, though so far we aren't using anything that's not possible with the implicit interface. -You can also write +You can also write + ```julia -[temp ~ furnace_off_threshold] => ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, i, c +[temp ~ furnace_off_threshold] => ModelingToolkit.ImperativeAffect(modified = (; + furnace_on)) do x, o, i, c @set! x.furnace_on = false end ``` + and it would work the same. The `ImperativeAffect` is the larger change in this example. `ImperativeAffect` has the constructor signature + ```julia - ImperativeAffect(f::Function; modified::NamedTuple, observed::NamedTuple, ctx) +ImperativeAffect(f::Function; modified::NamedTuple, observed::NamedTuple, ctx) ``` + that accepts the function to call, a named tuple of both the names of and symbolic values representing values in the system to be modified, a named tuple of the values that are merely observed (that is, used from the system but not modified), and a context that's passed to the affect function. In our example, each event merely changes whether the furnace is on or off. Accordingly, we pass a `modified` tuple -`(; furnace_on)` (creating a `NamedTuple` equivalent to `(furnace_on = furnace_on)`). `ImperativeAffect` will then +`(; furnace_on)` (creating a `NamedTuple` equivalent to `(furnace_on = furnace_on)`). `ImperativeAffect` will then evaluate this before calling our function to fill out all of the numerical values, then apply them back to the system once our affect function returns. Furthermore, it will check that it is possible to do this assignment. The function given to `ImperativeAffect` needs to have one of four signatures, checked in this order: -* `f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple` if the function needs the low-level integrator, -* `f(modified::NamedTuple, observed::NamedTuple, ctx)::NamedTuple` if the function needs the user-defined context, -* `f(modified::NamedTuple, observed::NamedTuple)::NamedTuple` if the function also reads observed values from the system, -* `f(modified::NamedTuple)::NamedTuple` if the function only writes values (unknowns or parameters) to the system. -The `do` block in the example implicitly constructs said function inline. For exposition, we use the full version (e.g. `x, o, i, c`) but this could be simplified to merely `x`. + + - `f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple` if the function needs the low-level integrator, + - `f(modified::NamedTuple, observed::NamedTuple, ctx)::NamedTuple` if the function needs the user-defined context, + - `f(modified::NamedTuple, observed::NamedTuple)::NamedTuple` if the function also reads observed values from the system, + - `f(modified::NamedTuple)::NamedTuple` if the function only writes values (unknowns or parameters) to the system. + The `do` block in the example implicitly constructs said function inline. For exposition, we use the full version (e.g. `x, o, i, c`) but this could be simplified to merely `x`. The function `f` will be called with `observed` and `modified` `NamedTuple`s that are derived from their respective `NamedTuple` definitions. -In our example, if `furnace_on` is `false`, then the value of the `x` that's passed in as `modified` will be `(furnace_on = false)`. +In our example, if `furnace_on` is `false`, then the value of the `x` that's passed in as `modified` will be `(furnace_on = false)`. The modified values should be passed out in the same format: to set `furnace_on` to `true` we need to return a tuple `(furnace_on = true)`. We use Setfield to do this in the example, recreating the result tuple before returning it. -Accordingly, we can now interpret the `ImperativeAffect` definitions to mean that when `temp = furnace_off_threshold` we +Accordingly, we can now interpret the `ImperativeAffect` definitions to mean that when `temp = furnace_off_threshold` we will write `furnace_on = false` back to the system, and when `temp = furnace_on_threshold` we will write `furnace_on = true` back to the system. @@ -468,7 +478,8 @@ ss = structural_simplify(sys) prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 10.0)) sol = solve(prob, Tsit5()) plot(sol) -hline!([sol.ps[furnace_off_threshold], sol.ps[furnace_on_threshold]], l = (:black, 1), primary = false) +hline!([sol.ps[furnace_off_threshold], sol.ps[furnace_on_threshold]], + l = (:black, 1), primary = false) ``` Here we see exactly the desired hysteresis. The heater starts on until the temperature hits @@ -477,71 +488,76 @@ point the furnace turns on again until `furnace_off_threshold` and so on and so is effectively regulating the temperature of the plant. ### [Quadrature Encoder](@id quadrature) + For a more complex application we'll look at modeling a quadrature encoder attached to a shaft spinning at a constant speed. Traditionally, a quadrature encoder is built out of a code wheel that interrupts the sensors at constant intervals and two sensors slightly out of phase with one another. A state machine can take the pattern of pulses produced by the two sensors and determine the number of steps that the shaft has spun. The state machine takes the new value from each sensor and the old values and decodes them into the direction that the wheel has spun in this step. ```@example events - @variables theta(t) omega(t) - params = @parameters qA=0 qB=0 hA=0 hB=0 cnt::Int=0 - eqs = [D(theta) ~ omega - omega ~ 1.0] +@variables theta(t) omega(t) +params = @parameters qA=0 qB=0 hA=0 hB=0 cnt::Int=0 +eqs = [D(theta) ~ omega + omega ~ 1.0] ``` + Our continous-time system is extremely simple. We have two states, `theta` for the angle of the shaft -and `omega` for the rate at which it's spinning. We then have parameters for the state machine `qA, qB, hA, hB` +and `omega` for the rate at which it's spinning. We then have parameters for the state machine `qA, qB, hA, hB` and a step count `cnt`. We'll then implement the decoder as a simple Julia function. + ```@example events - function decoder(oldA, oldB, newA, newB) - state = (oldA, oldB, newA, newB) - if state == (0, 0, 1, 0) || state == (1, 0, 1, 1) || state == (1, 1, 0, 1) || - state == (0, 1, 0, 0) - return 1 - elseif state == (0, 0, 0, 1) || state == (0, 1, 1, 1) || state == (1, 1, 1, 0) || - state == (1, 0, 0, 0) - return -1 - elseif state == (0, 0, 0, 0) || state == (0, 1, 0, 1) || state == (1, 0, 1, 0) || - state == (1, 1, 1, 1) - return 0 - else - return 0 # err is interpreted as no movement - end +function decoder(oldA, oldB, newA, newB) + state = (oldA, oldB, newA, newB) + if state == (0, 0, 1, 0) || state == (1, 0, 1, 1) || state == (1, 1, 0, 1) || + state == (0, 1, 0, 0) + return 1 + elseif state == (0, 0, 0, 1) || state == (0, 1, 1, 1) || state == (1, 1, 1, 0) || + state == (1, 0, 0, 0) + return -1 + elseif state == (0, 0, 0, 0) || state == (0, 1, 0, 1) || state == (1, 0, 1, 0) || + state == (1, 1, 1, 1) + return 0 + else + return 0 # err is interpreted as no movement end +end ``` + Based on the current and old state, this function will return 1 if the wheel spun in the positive direction, -1 if in the negative, and 0 otherwise. -The encoder state advances when the occlusion begins or ends. We model the +The encoder state advances when the occlusion begins or ends. We model the code wheel as simply detecting when `cos(100*theta)` is 0; if we're at a positive edge of the 0 crossing, then we interpret that as occlusion (so the discrete `qA` goes to 1). Otherwise, if `cos` is going negative, we interpret that as lack of occlusion (so the discrete goes to 0). The decoder function is then invoked to update the count with this new information. We can implement this in one of two ways: using edge sign detection or right root finding. For exposition, we -will implement each sensor differently. +will implement each sensor differently. For sensor A, we're using the edge detction method. By providing a different affect to `SymbolicContinuousCallback`'s `affect_neg` argument, we can specify different behaviour for the negative crossing vs. the positive crossing of the root. In our encoder, we interpret this as occlusion or nonocclusion of the sensor, update the internal state, and tick the decoder. + ```@example events - qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0], - ModelingToolkit.ImperativeAffect((; qA, hA, hB, cnt), (; qB)) do x, o, i, c - @set! x.hA = x.qA - @set! x.hB = o.qB - @set! x.qA = 1 - @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) - x - end, - affect_neg = ModelingToolkit.ImperativeAffect( - (; qA, hA, hB, cnt), (; qB)) do x, o, c, i - @set! x.hA = x.qA - @set! x.hB = o.qB - @set! x.qA = 0 - @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) - x - end) +qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0], + ModelingToolkit.ImperativeAffect((; qA, hA, hB, cnt), (; qB)) do x, o, i, c + @set! x.hA = x.qA + @set! x.hB = o.qB + @set! x.qA = 1 + @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) + x + end, + affect_neg = ModelingToolkit.ImperativeAffect( + (; qA, hA, hB, cnt), (; qB)) do x, o, c, i + @set! x.hA = x.qA + @set! x.hB = o.qB + @set! x.qA = 0 + @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) + x + end) ``` The other way we can implement a sensor is by changing the root find. @@ -550,29 +566,32 @@ the root is crossed. This makes it trickier to figure out what the new state is. Instead, we can use right root finding: ```@example events - qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π / 2) ~ 0], - ModelingToolkit.ImperativeAffect((; qB, hA, hB, cnt), (; qA, theta)) do x, o, i, c - @set! x.hA = o.qA - @set! x.hB = x.qB - @set! x.qB = clamp(sign(cos(100 * o.theta - π / 2)), 0.0, 1.0) - @set! x.cnt += decoder(x.hA, x.hB, o.qA, x.qB) - x - end; rootfind = SciMLBase.RightRootFind) +qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π / 2) ~ 0], + ModelingToolkit.ImperativeAffect((; qB, hA, hB, cnt), (; qA, theta)) do x, o, i, c + @set! x.hA = o.qA + @set! x.hB = x.qB + @set! x.qB = clamp(sign(cos(100 * o.theta - π / 2)), 0.0, 1.0) + @set! x.cnt += decoder(x.hA, x.hB, o.qA, x.qB) + x + end; rootfind = SciMLBase.RightRootFind) ``` -Here, sensor B is located `π / 2` behind sensor A in angular space, so we're adjusting our + +Here, sensor B is located `π / 2` behind sensor A in angular space, so we're adjusting our trigger function accordingly. We here ask for right root finding on the callback, so we know -that the value of said function will have the "new" sign rather than the old one. Thus, we can +that the value of said function will have the "new" sign rather than the old one. Thus, we can determine the new state of the sensor from the sign of the indicator function evaluated at the affect activation point, with -1 mapped to 0. We can now simulate the encoder. + ```@example events - @named sys = ODESystem( - eqs, t, [theta, omega], params; continuous_events = [qAevt, qBevt]) - ss = structural_simplify(sys) - prob = ODEProblem(ss, [theta => 0.0], (0.0, pi)) - sol = solve(prob, Tsit5(); dtmax = 0.01) - sol.ps[cnt] +@named sys = ODESystem( + eqs, t, [theta, omega], params; continuous_events = [qAevt, qBevt]) +ss = structural_simplify(sys) +prob = ODEProblem(ss, [theta => 0.0], (0.0, pi)) +sol = solve(prob, Tsit5(); dtmax = 0.01) +sol.ps[cnt] ``` + `cos(100*theta)` will have 200 crossings in the half rotation we've gone through, so the encoder would notionally count 200 steps. -Our encoder counts 198 steps (it loses one step to initialization and one step due to the final state falling squarely on an edge). \ No newline at end of file +Our encoder counts 198 steps (it loses one step to initialization and one step due to the final state falling squarely on an edge). diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index e99911eec6..1cc8273b4d 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -444,7 +444,7 @@ function build_explicit_observed_function(sys, ts; param_only = false, op = Operator, throw = true, - array_type=:array) + array_type = :array) if (isscalar = symbolic_type(ts) !== NotSymbolic()) ts = [ts] end @@ -589,10 +589,12 @@ function build_explicit_observed_function(sys, ts; oop_mtkp_wrapper = mtkparams_wrapper end - output_expr = isscalar ? ts[1] : (array_type == :array ? MakeArray(ts, output_type) : MakeTuple(ts)) + output_expr = isscalar ? ts[1] : + (array_type == :array ? MakeArray(ts, output_type) : MakeTuple(ts)) # Need to keep old method of building the function since it uses `output_type`, # which can't be provided to `build_function` - oop_fn = Func(args, [], pre(Let(obsexprs, output_expr, false))) |> array_wrapper[1] |> oop_mtkp_wrapper |> toexpr + oop_fn = Func(args, [], pre(Let(obsexprs, output_expr, false))) |> array_wrapper[1] |> + oop_mtkp_wrapper |> toexpr oop_fn = expression ? oop_fn : eval_or_rgf(oop_fn; eval_expression, eval_module) if !isscalar From 711fb8c654be38ce7377b85f3025ee7934cd5de8 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Tue, 22 Oct 2024 18:44:45 -0700 Subject: [PATCH 27/47] Spelling --- docs/src/basics/Events.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index 4a59149308..8088e0f52a 100644 --- a/docs/src/basics/Events.md +++ b/docs/src/basics/Events.md @@ -388,9 +388,9 @@ system, simplifying the definitions and (in the future) allowing assignments bac by solving the nonlinear reinitialization problem afterwards. We will use two examples to describe `ImperativeAffect`: a simple heater and a quadrature encoder. -These examples will also demonstrate advanced usage of `ModelingToolkit.SymbolicContinousCallback`, +These examples will also demonstrate advanced usage of `ModelingToolkit.SymbolicContinuousCallback`, the low-level interface that the aforementioned tuple form converts into and allows control over the -exact SciMLCallbacks event that is generated for a continous event. +exact SciMLCallbacks event that is generated for a continuous event. ### [Heater](@id heater_events) @@ -405,7 +405,7 @@ eqs = [ ``` Our plant is simple. We have a heater that's turned on and off by the clocked parameter `furnace_on` -which adds `furnace_power` forcing to the system when enabled. We then leak heat porportional to `leakage` +which adds `furnace_power` forcing to the system when enabled. We then leak heat proportional to `leakage` as a function of the square of the current temperature. We need a controller with hysteresis to conol the plant. We wish the furnace to turn on when the temperature @@ -537,7 +537,7 @@ then invoked to update the count with this new information. We can implement this in one of two ways: using edge sign detection or right root finding. For exposition, we will implement each sensor differently. -For sensor A, we're using the edge detction method. By providing a different affect to `SymbolicContinuousCallback`'s +For sensor A, we're using the edge detection method. By providing a different affect to `SymbolicContinuousCallback`'s `affect_neg` argument, we can specify different behaviour for the negative crossing vs. the positive crossing of the root. In our encoder, we interpret this as occlusion or nonocclusion of the sensor, update the internal state, and tick the decoder. From a8ea3698e71afe91dd9040a6fb80bab22c7d6c95 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Tue, 22 Oct 2024 18:46:54 -0700 Subject: [PATCH 28/47] Spelling (2) --- docs/src/basics/Events.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index 8088e0f52a..e583af8d9e 100644 --- a/docs/src/basics/Events.md +++ b/docs/src/basics/Events.md @@ -410,7 +410,7 @@ as a function of the square of the current temperature. We need a controller with hysteresis to conol the plant. We wish the furnace to turn on when the temperature is below `furnace_on_threshold` and off when above `furnace_off_threshold`, while maintaining its current state -in between. To do this, we create two continous callbacks: +in between. To do this, we create two continuous callbacks: ```@example events using Setfield @@ -501,7 +501,7 @@ eqs = [D(theta) ~ omega omega ~ 1.0] ``` -Our continous-time system is extremely simple. We have two states, `theta` for the angle of the shaft +Our continuous-time system is extremely simple. We have two states, `theta` for the angle of the shaft and `omega` for the rate at which it's spinning. We then have parameters for the state machine `qA, qB, hA, hB` and a step count `cnt`. From 6a304e7d7f8960d7f2cc63d1be93d8f41f6db824 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Tue, 22 Oct 2024 18:48:15 -0700 Subject: [PATCH 29/47] Remove debug RGF shim --- test/symbolic_events.jl | 17 +---------------- 1 file changed, 1 insertion(+), 16 deletions(-) diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index f4f97fb2b9..96a57ec40e 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -1217,19 +1217,4 @@ end prob = ODEProblem(ss, [theta => 0.0], (0.0, pi)) sol = solve(prob, Tsit5(); dtmax = 0.01) @test getp(sol, cnt)(sol) == 197 # we get 2 pulses per phase cycle (cos 0 crossing) and we go to 100 cycles; we miss a few due to the initial state -end - -import RuntimeGeneratedFunctions -function (f::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{ - argnames, cache_tag, context_tag, - id})(args::Vararg{Any, N}) where {N, argnames, cache_tag, context_tag, id} - try - RuntimeGeneratedFunctions.generated_callfunc(f, args...) - catch e - @error "Caught error in RuntimeGeneratedFunction; source code follows" - func_expr = Expr(:->, Expr(:tuple, argnames...), - RuntimeGeneratedFunctions._lookup_body(cache_tag, id)) - @show func_expr - rethrow(e) - end -end +end \ No newline at end of file From 35ec1c59b31ab3d0c6a421482d7e48459a00111c Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Tue, 22 Oct 2024 18:50:32 -0700 Subject: [PATCH 30/47] Formatter (3) --- test/symbolic_events.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index 96a57ec40e..c021f99eea 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -1217,4 +1217,4 @@ end prob = ODEProblem(ss, [theta => 0.0], (0.0, pi)) sol = solve(prob, Tsit5(); dtmax = 0.01) @test getp(sol, cnt)(sol) == 197 # we get 2 pulses per phase cycle (cos 0 crossing) and we go to 100 cycles; we miss a few due to the initial state -end \ No newline at end of file +end From 4c3d6fcc1eaedbec33cfd3f979a958110ff4eace Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 23 Oct 2024 06:21:05 +0000 Subject: [PATCH 31/47] Update docs/src/basics/Events.md Co-authored-by: Fredrik Bagge Carlson --- docs/src/basics/Events.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index e583af8d9e..41fa96d619 100644 --- a/docs/src/basics/Events.md +++ b/docs/src/basics/Events.md @@ -394,7 +394,7 @@ exact SciMLCallbacks event that is generated for a continuous event. ### [Heater](@id heater_events) -Bang-bang control of a heater connected to a leaky plant requires hysteresis in order to prevent control oscillation. +Bang-bang control of a heater connected to a leaky plant requires hysteresis in order to prevent rapid control oscillation. ```@example events @variables temp(t) From cdddeb9f62ffcd12357dc96648167fe8514993be Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 23 Oct 2024 06:21:34 +0000 Subject: [PATCH 32/47] Update docs/src/basics/Events.md Co-authored-by: Fredrik Bagge Carlson --- docs/src/basics/Events.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index 41fa96d619..add870b7f3 100644 --- a/docs/src/basics/Events.md +++ b/docs/src/basics/Events.md @@ -408,7 +408,7 @@ Our plant is simple. We have a heater that's turned on and off by the clocked pa which adds `furnace_power` forcing to the system when enabled. We then leak heat proportional to `leakage` as a function of the square of the current temperature. -We need a controller with hysteresis to conol the plant. We wish the furnace to turn on when the temperature +We need a controller with hysteresis to control the plant. We wish the furnace to turn on when the temperature is below `furnace_on_threshold` and off when above `furnace_off_threshold`, while maintaining its current state in between. To do this, we create two continuous callbacks: From 2ac2d230664916fdd0151cb6257adbf9dd7b2ccb Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 23 Oct 2024 06:21:54 +0000 Subject: [PATCH 33/47] Update docs/src/basics/Events.md Co-authored-by: Fredrik Bagge Carlson --- docs/src/basics/Events.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index add870b7f3..0d9149e23c 100644 --- a/docs/src/basics/Events.md +++ b/docs/src/basics/Events.md @@ -561,7 +561,7 @@ qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0], ``` The other way we can implement a sensor is by changing the root find. -Normally, we use left root finding; the affect will be invoked instantaneously before +Normally, we use left root finding; the affect will be invoked instantaneously _before_ the root is crossed. This makes it trickier to figure out what the new state is. Instead, we can use right root finding: From 9bf734dcbc03123b72fb179a47370d9aa484869f Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Wed, 23 Oct 2024 15:28:56 -0700 Subject: [PATCH 34/47] Simplify callback interface, fix references --- docs/src/basics/Events.md | 22 +++++++-------- src/systems/callbacks.jl | 31 ++++++-------------- test/symbolic_events.jl | 59 +++------------------------------------ 3 files changed, 23 insertions(+), 89 deletions(-) diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index 0d9149e23c..de1c1000ce 100644 --- a/docs/src/basics/Events.md +++ b/docs/src/basics/Events.md @@ -416,12 +416,12 @@ in between. To do this, we create two continuous callbacks: using Setfield furnace_disable = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_off_threshold], - ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, i, c + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, c, i @set! x.furnace_on = false end) furnace_enable = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_on_threshold], - ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, i, c + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, c, i @set! x.furnace_on = true end) ``` @@ -454,18 +454,16 @@ In our example, each event merely changes whether the furnace is on or off. Acco evaluate this before calling our function to fill out all of the numerical values, then apply them back to the system once our affect function returns. Furthermore, it will check that it is possible to do this assignment. -The function given to `ImperativeAffect` needs to have one of four signatures, checked in this order: - - - `f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple` if the function needs the low-level integrator, - - `f(modified::NamedTuple, observed::NamedTuple, ctx)::NamedTuple` if the function needs the user-defined context, - - `f(modified::NamedTuple, observed::NamedTuple)::NamedTuple` if the function also reads observed values from the system, - - `f(modified::NamedTuple)::NamedTuple` if the function only writes values (unknowns or parameters) to the system. - The `do` block in the example implicitly constructs said function inline. For exposition, we use the full version (e.g. `x, o, i, c`) but this could be simplified to merely `x`. +The function given to `ImperativeAffect` needs to have the signature: +```julia + f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple +``` The function `f` will be called with `observed` and `modified` `NamedTuple`s that are derived from their respective `NamedTuple` definitions. In our example, if `furnace_on` is `false`, then the value of the `x` that's passed in as `modified` will be `(furnace_on = false)`. The modified values should be passed out in the same format: to set `furnace_on` to `true` we need to return a tuple `(furnace_on = true)`. -We use Setfield to do this in the example, recreating the result tuple before returning it. +The examples does this with Setfield, recreating the result tuple before returning it; the returned tuple may optionally be missing values as +well, in which case those values will not be written back to the problem. Accordingly, we can now interpret the `ImperativeAffect` definitions to mean that when `temp = furnace_off_threshold` we will write `furnace_on = false` back to the system, and when `temp = furnace_on_threshold` we will write `furnace_on = true` back @@ -543,7 +541,7 @@ In our encoder, we interpret this as occlusion or nonocclusion of the sensor, up ```@example events qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0], - ModelingToolkit.ImperativeAffect((; qA, hA, hB, cnt), (; qB)) do x, o, i, c + ModelingToolkit.ImperativeAffect((; qA, hA, hB, cnt), (; qB)) do x, o, c, i @set! x.hA = x.qA @set! x.hB = o.qB @set! x.qA = 1 @@ -567,7 +565,7 @@ Instead, we can use right root finding: ```@example events qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π / 2) ~ 0], - ModelingToolkit.ImperativeAffect((; qB, hA, hB, cnt), (; qA, theta)) do x, o, i, c + ModelingToolkit.ImperativeAffect((; qB, hA, hB, cnt), (; qA, theta)) do x, o, c, i @set! x.hA = o.qA @set! x.hB = x.qB @set! x.qB = clamp(sign(cos(100 * o.theta - π / 2)), 0.0, 1.0) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index e7198817b1..618a1f3a83 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -76,12 +76,11 @@ end `ImperativeAffect` is a helper for writing affect functions that will compute observed values and ensure that modified values are correctly written back into the system. The affect function `f` needs to have -one of four signatures: -* `f(modified::NamedTuple)::NamedTuple` if the function only writes values (unknowns or parameters) to the system, -* `f(modified::NamedTuple, observed::NamedTuple)::NamedTuple` if the function also reads observed values from the system, -* `f(modified::NamedTuple, observed::NamedTuple, ctx)::NamedTuple` if the function needs the user-defined context, -* `f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple` if the function needs the low-level integrator. -These will be checked in reverse order (that is, the four-argument version first, than the 3, etc). +the signature + +``` + f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple +``` The function `f` will be called with `observed` and `modified` `NamedTuple`s that are derived from their respective `NamedTuple` definitions. Each declaration`NamedTuple` should map an expression to a symbol; for example if we pass `observed=(; x = a + b)` this will alias the result of executing `a+b` in the system as `x` @@ -1046,7 +1045,7 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs. Implementation sketch: generate observed function (oop), should save to a component array under obs_syms do the same stuff as the normal FA for pars_syms - call the affect method - test if it's OOP or IP using applicable + call the affect method unpack and apply the resulting values =# function check_dups(syms, exprs) # = (syms_dedup, exprs_dedup) @@ -1135,22 +1134,10 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs. integ.u, integ.p, integ.t)) # let the user do their thing - modvals = if applicable( - user_affect, upd_component_array, obs_component_array, ctx, integ) - user_affect(upd_component_array, obs_component_array, ctx, integ) - elseif applicable(user_affect, upd_component_array, obs_component_array, ctx) - user_affect(upd_component_array, obs_component_array, ctx) - elseif applicable(user_affect, upd_component_array, obs_component_array) - user_affect(upd_component_array, obs_component_array) - elseif applicable(user_affect, upd_component_array) - user_affect(upd_component_array) - else - @error "User affect function $user_affect needs to implement one of the supported ImperativeAffect callback forms; see the ImperativeAffect docstring for more details" - user_affect(upd_component_array, obs_component_array, integ, ctx) # this WILL error but it'll give a more sensible message - end - + upd_vals = user_affect(upd_component_array, obs_component_array, ctx, integ) + # write the new values back to the integrator - _generated_writeback(integ, upd_funs, modvals) + _generated_writeback(integ, upd_funs, upd_vals) for idx in save_idxs SciMLBase.save_discretes!(integ, idx) diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index c021f99eea..3c3e9168e7 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -1027,60 +1027,9 @@ end furnace_off = ModelingToolkit.SymbolicContinuousCallback( [temp ~ furnace_off_threshold], - ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, i - @set! x.furnace_on = false - end) - furnace_enable = ModelingToolkit.SymbolicContinuousCallback( - [temp ~ furnace_on_threshold], - ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, i - @set! x.furnace_on = true - end) - @named sys = ODESystem( - eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) - ss = structural_simplify(sys) - prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) - sol = solve(prob, Tsit5(); dtmax = 0.01) - @test all(sol[temp][sol.t .> 1.0] .<= 0.79) && all(sol[temp][sol.t .> 1.0] .>= 0.49) - - furnace_off = ModelingToolkit.SymbolicContinuousCallback( - [temp ~ furnace_off_threshold], - ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o - @set! x.furnace_on = false - end) - furnace_enable = ModelingToolkit.SymbolicContinuousCallback( - [temp ~ furnace_on_threshold], - ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o - @set! x.furnace_on = true - end) - @named sys = ODESystem( - eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) - ss = structural_simplify(sys) - prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) - sol = solve(prob, Tsit5(); dtmax = 0.01) - @test all(sol[temp][sol.t .> 1.0] .<= 0.79) && all(sol[temp][sol.t .> 1.0] .>= 0.49) - - furnace_off = ModelingToolkit.SymbolicContinuousCallback( - [temp ~ furnace_off_threshold], - ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x - @set! x.furnace_on = false - end) - furnace_enable = ModelingToolkit.SymbolicContinuousCallback( - [temp ~ furnace_on_threshold], - ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x - @set! x.furnace_on = true - end) - @named sys = ODESystem( - eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) - ss = structural_simplify(sys) - prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) - sol = solve(prob, Tsit5(); dtmax = 0.01) - @test all(sol[temp][sol.t .> 1.0] .<= 0.79) && all(sol[temp][sol.t .> 1.0] .>= 0.49) - - furnace_off = ModelingToolkit.SymbolicContinuousCallback( - [temp ~ furnace_off_threshold], - ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, c, i @set! x.furnace_on = false - end; initialize = ModelingToolkit.ImperativeAffect(modified = (; temp)) do x + end; initialize = ModelingToolkit.ImperativeAffect(modified = (; temp)) do x, o, c, i @set! x.temp = 0.2 end) furnace_enable = ModelingToolkit.SymbolicContinuousCallback( @@ -1180,7 +1129,7 @@ end end end qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0], - ModelingToolkit.ImperativeAffect((; qA, hA, hB, cnt), (; qB)) do x, o, i, c + ModelingToolkit.ImperativeAffect((; qA, hA, hB, cnt), (; qB)) do x, o, c, i @set! x.hA = x.qA @set! x.hB = o.qB @set! x.qA = 1 @@ -1196,7 +1145,7 @@ end x end; rootfind = SciMLBase.RightRootFind) qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π / 2) ~ 0], - ModelingToolkit.ImperativeAffect((; qB, hA, hB, cnt), (; qA)) do x, o, i, c + ModelingToolkit.ImperativeAffect((; qB, hA, hB, cnt), (; qA)) do x, o, c, i @set! x.hA = o.qA @set! x.hB = x.qB @set! x.qB = 1 From 1a3f7d4e4916c12d26a4659334a859a9c3d5ca8b Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Wed, 23 Oct 2024 15:33:20 -0700 Subject: [PATCH 35/47] Make array_type an actual type, though only in a limited sense --- src/systems/callbacks.jl | 4 ++-- src/systems/diffeqs/odesystem.jl | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index 618a1f3a83..c34681a20c 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -1105,7 +1105,7 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs. end obs_fun = build_explicit_observed_function( sys, Symbolics.scalarize.(obs_exprs); - array_type = :tuple) + array_type = Tuple) obs_sym_tuple = (obs_syms...,) # okay so now to generate the stuff to assign it back into the system @@ -1113,7 +1113,7 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs. mod_names = (mod_syms...,) mod_og_val_fun = build_explicit_observed_function( sys, Symbolics.scalarize.(first.(mod_pairs)); - array_type = :tuple) + array_type = Tuple) upd_funs = NamedTuple{mod_names}((setu.((sys,), first.(mod_pairs))...,)) diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index 1cc8273b4d..374cbab5f8 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -429,7 +429,7 @@ Options not otherwise specified are: * `op = Operator` sets the recursion terminator for the walk done by `vars` to identify the variables that appear in `ts`. See the documentation for `vars` for more detail. * `throw = true` if true, throw an error when generating a function for `ts` that reference variables that do not exist * `drop_expr` is deprecated. -* `array_type`; only used if the output is an array (that is, `!isscalar(ts)`). If `:array`, then it will generate an array, if `:tuple` then it will generate a tuple. +* `array_type`; only used if the output is an array (that is, `!isscalar(ts)`). If it is `Vector`, then it will generate an array, if `Tuple` then it will generate a tuple. """ function build_explicit_observed_function(sys, ts; inputs = nothing, @@ -444,7 +444,7 @@ function build_explicit_observed_function(sys, ts; param_only = false, op = Operator, throw = true, - array_type = :array) + array_type = Vector) if (isscalar = symbolic_type(ts) !== NotSymbolic()) ts = [ts] end @@ -590,7 +590,7 @@ function build_explicit_observed_function(sys, ts; end output_expr = isscalar ? ts[1] : - (array_type == :array ? MakeArray(ts, output_type) : MakeTuple(ts)) + (array_type <: Vector ? MakeArray(ts, output_type) : MakeTuple(ts)) # Need to keep old method of building the function since it uses `output_type`, # which can't be provided to `build_function` oop_fn = Func(args, [], pre(Let(obsexprs, output_expr, false))) |> array_wrapper[1] |> From 6c7e4b84ae40e0bf1cbe1c61f6e4c41202ee449c Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Wed, 23 Oct 2024 15:34:59 -0700 Subject: [PATCH 36/47] Clean up language indocs --- docs/src/basics/Events.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index de1c1000ce..902be98fd4 100644 --- a/docs/src/basics/Events.md +++ b/docs/src/basics/Events.md @@ -404,7 +404,7 @@ eqs = [ ] ``` -Our plant is simple. We have a heater that's turned on and off by the clocked parameter `furnace_on` +Our plant is simple. We have a heater that's turned on and off by the time-indexed parameter `furnace_on` which adds `furnace_power` forcing to the system when enabled. We then leak heat proportional to `leakage` as a function of the square of the current temperature. @@ -499,9 +499,9 @@ eqs = [D(theta) ~ omega omega ~ 1.0] ``` -Our continuous-time system is extremely simple. We have two states, `theta` for the angle of the shaft +Our continuous-time system is extremely simple. We have two unknown variables `theta` for the angle of the shaft and `omega` for the rate at which it's spinning. We then have parameters for the state machine `qA, qB, hA, hB` -and a step count `cnt`. +(corresponding to the current quadrature of the A/B sensors and the historical ones) and a step count `cnt`. We'll then implement the decoder as a simple Julia function. From 06ecd2ee541d94bcf92c39543583f9bc0646286d Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Wed, 23 Oct 2024 21:24:46 -0700 Subject: [PATCH 37/47] Format --- docs/src/basics/Events.md | 7 ++++--- src/systems/callbacks.jl | 2 +- test/symbolic_events.jl | 3 ++- 3 files changed, 7 insertions(+), 5 deletions(-) diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index 902be98fd4..4c9803ef57 100644 --- a/docs/src/basics/Events.md +++ b/docs/src/basics/Events.md @@ -456,13 +456,14 @@ once our affect function returns. Furthermore, it will check that it is possible The function given to `ImperativeAffect` needs to have the signature: -```julia - f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple +```julia +f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple ``` + The function `f` will be called with `observed` and `modified` `NamedTuple`s that are derived from their respective `NamedTuple` definitions. In our example, if `furnace_on` is `false`, then the value of the `x` that's passed in as `modified` will be `(furnace_on = false)`. The modified values should be passed out in the same format: to set `furnace_on` to `true` we need to return a tuple `(furnace_on = true)`. -The examples does this with Setfield, recreating the result tuple before returning it; the returned tuple may optionally be missing values as +The examples does this with Setfield, recreating the result tuple before returning it; the returned tuple may optionally be missing values as well, in which case those values will not be written back to the problem. Accordingly, we can now interpret the `ImperativeAffect` definitions to mean that when `temp = furnace_off_threshold` we diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index c34681a20c..1061139a8a 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -1135,7 +1135,7 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs. # let the user do their thing upd_vals = user_affect(upd_component_array, obs_component_array, ctx, integ) - + # write the new values back to the integrator _generated_writeback(integ, upd_funs, upd_vals) diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index 3c3e9168e7..77961cb01f 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -1029,7 +1029,8 @@ end [temp ~ furnace_off_threshold], ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, c, i @set! x.furnace_on = false - end; initialize = ModelingToolkit.ImperativeAffect(modified = (; temp)) do x, o, c, i + end; initialize = ModelingToolkit.ImperativeAffect(modified = (; + temp)) do x, o, c, i @set! x.temp = 0.2 end) furnace_enable = ModelingToolkit.SymbolicContinuousCallback( From 7f7f65cc3e9dde8fec99037f4df27221088ed818 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Fri, 25 Oct 2024 10:17:26 -0700 Subject: [PATCH 38/47] Clear up some of the documentation language --- docs/src/basics/Events.md | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index 4c9803ef57..90b1b4036d 100644 --- a/docs/src/basics/Events.md +++ b/docs/src/basics/Events.md @@ -383,14 +383,13 @@ It can be seen that the timeseries for `c` is not saved. The `ImperativeAffect` can be used as an alternative to the aforementioned functional affect form. Note that `ImperativeAffect` is still experimental; to emphasize this, we do not export it and it should be -included as `ModelingToolkit.ImperativeAffect`. It abstracts over how values are written back to the -system, simplifying the definitions and (in the future) allowing assignments back to observed values -by solving the nonlinear reinitialization problem afterwards. +included as `ModelingToolkit.ImperativeAffect`. `ImperativeAffect` aims to simplify the manipulation of +system state. We will use two examples to describe `ImperativeAffect`: a simple heater and a quadrature encoder. These examples will also demonstrate advanced usage of `ModelingToolkit.SymbolicContinuousCallback`, -the low-level interface that the aforementioned tuple form converts into and allows control over the -exact SciMLCallbacks event that is generated for a continuous event. +the low-level interface of the tuple form converts into that allows control over the SciMLBase-level +event that is generated for a continuous event. ### [Heater](@id heater_events) From e750219fcb2d746b5b3bf2c52bcf22c19c48955c Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Fri, 25 Oct 2024 10:49:42 -0700 Subject: [PATCH 39/47] Format --- docs/src/basics/Events.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index 90b1b4036d..23e1e6d7d1 100644 --- a/docs/src/basics/Events.md +++ b/docs/src/basics/Events.md @@ -383,12 +383,12 @@ It can be seen that the timeseries for `c` is not saved. The `ImperativeAffect` can be used as an alternative to the aforementioned functional affect form. Note that `ImperativeAffect` is still experimental; to emphasize this, we do not export it and it should be -included as `ModelingToolkit.ImperativeAffect`. `ImperativeAffect` aims to simplify the manipulation of +included as `ModelingToolkit.ImperativeAffect`. `ImperativeAffect` aims to simplify the manipulation of system state. We will use two examples to describe `ImperativeAffect`: a simple heater and a quadrature encoder. These examples will also demonstrate advanced usage of `ModelingToolkit.SymbolicContinuousCallback`, -the low-level interface of the tuple form converts into that allows control over the SciMLBase-level +the low-level interface of the tuple form converts into that allows control over the SciMLBase-level event that is generated for a continuous event. ### [Heater](@id heater_events) From 3bc5d236053cebf4b08315f64379c3666d5c9e84 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Tue, 12 Nov 2024 16:29:47 -0800 Subject: [PATCH 40/47] Fix merge issues --- src/systems/callbacks.jl | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index c2c895626c..fab573cffb 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -247,8 +247,6 @@ struct SymbolicContinuousCallback initialize = NULL_AFFECT, finalize = NULL_AFFECT, rootfind = SciMLBase.LeftRootFind, - initialize = NULL_AFFECT, - finalize = NULL_AFFECT, reinitializealg = SciMLBase.CheckInit()) new(eqs, initialize, finalize, make_affect(affect), make_affect(affect_neg), rootfind, reinitializealg) @@ -387,16 +385,6 @@ function affect_negs(cbs::Vector{SymbolicContinuousCallback}) mapreduce(affect_negs, vcat, cbs, init = Equation[]) 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.initialize -function finalize_affects(cbs::Vector{SymbolicContinuousCallback}) - mapreduce(finalize_affects, vcat, cbs, init = Equation[]) -end - reinitialization_alg(cb::SymbolicContinuousCallback) = cb.reinitializealg function reinitialization_algs(cbs::Vector{SymbolicContinuousCallback}) mapreduce( From e69ef194e685fb94016541255030b75764acbbb5 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Tue, 12 Nov 2024 23:45:49 -0800 Subject: [PATCH 41/47] Clean up usage of custom init --- src/systems/callbacks.jl | 17 +++++------------ 1 file changed, 5 insertions(+), 12 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index fab573cffb..9c8187f6b7 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -922,18 +922,11 @@ function generate_vector_rootfinding_callback( (cb, fn) -> begin if (save_idxs = get(ic.callback_to_clocks, cb, nothing)) !== nothing let save_idxs = save_idxs - if !isnothing(fn.initialize) - (i) -> begin - fn.initialize(i) - for idx in save_idxs - SciMLBase.save_discretes!(i, idx) - end - end - else - (i) -> begin - for idx in save_idxs - SciMLBase.save_discretes!(i, idx) - end + 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 From 3c6b37ed3852d3590afe9b29f9690b0824f84a09 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Tue, 12 Nov 2024 23:46:52 -0800 Subject: [PATCH 42/47] Simplify setup function construction --- src/systems/callbacks.jl | 27 ++++++++++++--------------- 1 file changed, 12 insertions(+), 15 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index 9c8187f6b7..b9423b8f07 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -918,24 +918,21 @@ function generate_vector_rootfinding_callback( initialize = nothing if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing initialize = handle_optional_setup_fn( - map( - (cb, fn) -> begin - 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 + 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 - else - fn.initialize end - end, - cbs, - affect_functions), + else + fn.initialize + end + end, SciMLBase.INITIALIZE_DEFAULT) else From 16d3f5c7a57bffd698a6c69da3ea89e0089b2d54 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Mon, 9 Dec 2024 02:52:15 -0800 Subject: [PATCH 43/47] Refactor ImperativeAffect into its own file --- src/ModelingToolkit.jl | 1 + src/systems/callbacks.jl | 221 +------------------------------ src/systems/diffeqs/odesystem.jl | 2 - src/systems/imperative_affect.jl | 220 ++++++++++++++++++++++++++++++ 4 files changed, 224 insertions(+), 220 deletions(-) create mode 100644 src/systems/imperative_affect.jl diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index de8e69c41f..ccdaa15e99 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -145,6 +145,7 @@ include("systems/parameter_buffer.jl") include("systems/abstractsystem.jl") include("systems/model_parsing.jl") include("systems/connectors.jl") +include("systems/imperative_affect.jl") include("systems/callbacks.jl") include("systems/problem_utils.jl") diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index 0cb7f16c9f..4533534029 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -73,111 +73,6 @@ function namespace_affect(affect::FunctionalAffect, s) context(affect)) end -""" - ImperativeAffect(f::Function; modified::NamedTuple, observed::NamedTuple, ctx) - -`ImperativeAffect` is a helper for writing affect functions that will compute observed values and -ensure that modified values are correctly written back into the system. The affect function `f` needs to have -the signature - -``` - f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple -``` - -The function `f` will be called with `observed` and `modified` `NamedTuple`s that are derived from their respective `NamedTuple` definitions. -Each declaration`NamedTuple` should map an expression to a symbol; for example if we pass `observed=(; x = a + b)` this will alias the result of executing `a+b` in the system as `x` -so the value of `a + b` will be accessible as `observed.x` in `f`. `modified` currently restricts symbolic expressions to only bare variables, so only tuples of the form -`(; x = y)` or `(; x)` (which aliases `x` as itself) are allowed. - -The argument NamedTuples (for instance `(;x=y)`) will be populated with the declared values on function entry; if we require `(;x=y)` in `observed` and `y=2`, for example, -then the NamedTuple `(;x=2)` will be passed as `observed` to the affect function `f`. - -The NamedTuple returned from `f` includes the values to be written back to the system after `f` returns. For example, if we want to update the value of `x` to be the result of `x + y` we could write - - ImperativeAffect(observed=(; x_plus_y = x + y), modified=(; x)) do m, o - @set! m.x = o.x_plus_y - end - -Where we use Setfield to copy the tuple `m` with a new value for `x`, then return the modified value of `m`. All values updated by the tuple must have names originally declared in -`modified`; a runtime error will be produced if a value is written that does not appear in `modified`. The user can dynamically decide not to write a value back by not including it -in the returned tuple, in which case the associated field will not be updated. -""" -@kwdef struct ImperativeAffect - f::Any - obs::Vector - obs_syms::Vector{Symbol} - modified::Vector - mod_syms::Vector{Symbol} - ctx::Any - skip_checks::Bool -end - -function ImperativeAffect(f::Function; - observed::NamedTuple = NamedTuple{()}(()), - modified::NamedTuple = NamedTuple{()}(()), - ctx = nothing, - skip_checks = false) - ImperativeAffect(f, - collect(values(observed)), collect(keys(observed)), - collect(values(modified)), collect(keys(modified)), - ctx, skip_checks) -end -function ImperativeAffect(f::Function, modified::NamedTuple; - observed::NamedTuple = NamedTuple{()}(()), ctx = nothing, skip_checks = false) - ImperativeAffect( - f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) -end -function ImperativeAffect( - f::Function, modified::NamedTuple, observed::NamedTuple; ctx = nothing, skip_checks = false) - ImperativeAffect( - f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) -end -function ImperativeAffect( - f::Function, modified::NamedTuple, observed::NamedTuple, ctx; skip_checks = false) - ImperativeAffect( - f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) -end - -function Base.show(io::IO, mfa::ImperativeAffect) - obs_vals = join(map((ob, nm) -> "$ob => $nm", mfa.obs, mfa.obs_syms), ", ") - mod_vals = join(map((md, nm) -> "$md => $nm", mfa.modified, mfa.mod_syms), ", ") - affect = mfa.f - print(io, - "ImperativeAffect(observed: [$obs_vals], modified: [$mod_vals], affect:$affect)") -end -func(f::ImperativeAffect) = f.f -context(a::ImperativeAffect) = a.ctx -observed(a::ImperativeAffect) = a.obs -observed_syms(a::ImperativeAffect) = a.obs_syms -discretes(a::ImperativeAffect) = filter(ModelingToolkit.isparameter, a.modified) -modified(a::ImperativeAffect) = a.modified -modified_syms(a::ImperativeAffect) = a.mod_syms - -function Base.:(==)(a1::ImperativeAffect, a2::ImperativeAffect) - isequal(a1.f, a2.f) && isequal(a1.obs, a2.obs) && isequal(a1.modified, a2.modified) && - isequal(a1.obs_syms, a2.obs_syms) && isequal(a1.mod_syms, a2.mod_syms) && - isequal(a1.ctx, a2.ctx) -end - -function Base.hash(a::ImperativeAffect, s::UInt) - s = hash(a.f, s) - s = hash(a.obs, s) - s = hash(a.obs_syms, s) - s = hash(a.modified, s) - s = hash(a.mod_syms, s) - hash(a.ctx, s) -end - -function namespace_affect(affect::ImperativeAffect, s) - ImperativeAffect(func(affect), - namespace_expr.(observed(affect), (s,)), - observed_syms(affect), - renamespace.((s,), modified(affect)), - modified_syms(affect), - context(affect), - affect.skip_checks) -end - function has_functional_affect(cb) (affects(cb) isa FunctionalAffect || affects(cb) isa ImperativeAffect) end @@ -203,13 +98,13 @@ sharp discontinuity between integrator steps (which in this example would not no guaranteed to be triggered. Once detected the integrator will "wind back" through a root-finding process to identify the point when the condition became active; the method used -is specified by `rootfind` from [`SciMLBase.RootfindOpt`](@ref). If we denote the time when the condition becomes active at tc, +is specified by `rootfind` from [`SciMLBase.RootfindOpt`](@ref). If we denote the time when the condition becomes active as `tc``, the value in the integrator after windback will be: * `u[tc-epsilon], p[tc-epsilon], tc` if `LeftRootFind` is used, * `u[tc+epsilon], p[tc+epsilon], tc` if `RightRootFind` is used, * or `u[t], p[t], t` if `NoRootFind` is used. For example, if we want to detect when an unknown variable `x` satisfies `x > 0` using the condition `x ~ 0` on a positive edge (that is, `D(x) > 0`), -then left root finding will get us `x=-epsilon`, right root finding `x=epsilon` and no root finding whatever the next step of the integrator was after +then left root finding will get us `x=-epsilon`, right root finding `x=epsilon` and no root finding will produce whatever the next step of the integrator was after it passed through 0. Multiple callbacks in the same system with different `rootfind` operations will be grouped @@ -405,7 +300,6 @@ 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(af::ImperativeAffect, s) = namespace_affect(af, s) namespace_affects(::Nothing, s) = nothing function namespace_callback(cb::SymbolicContinuousCallback, s)::SymbolicContinuousCallback @@ -480,7 +374,6 @@ scalarize_affects(affects) = scalarize(affects) scalarize_affects(affects::Tuple) = FunctionalAffect(affects...) scalarize_affects(affects::NamedTuple) = FunctionalAffect(; affects...) scalarize_affects(affects::FunctionalAffect) = affects -scalarize_affects(affects::ImperativeAffect) = affects SymbolicDiscreteCallback(p::Pair) = SymbolicDiscreteCallback(p[1], p[2]) SymbolicDiscreteCallback(cb::SymbolicDiscreteCallback) = cb # passthrough @@ -1099,117 +992,9 @@ function check_assignable(sys, sym) end end -function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs...) - #= - Implementation sketch: - generate observed function (oop), should save to a component array under obs_syms - do the same stuff as the normal FA for pars_syms - call the affect method - unpack and apply the resulting values - =# - function check_dups(syms, exprs) # = (syms_dedup, exprs_dedup) - seen = Set{Symbol}() - syms_dedup = [] - exprs_dedup = [] - for (sym, exp) in Iterators.zip(syms, exprs) - if !in(sym, seen) - push!(syms_dedup, sym) - push!(exprs_dedup, exp) - push!(seen, sym) - elseif !affect.skip_checks - @warn "Expression $(expr) is aliased as $sym, which has already been used. The first definition will be used." - end - end - return (syms_dedup, exprs_dedup) - end - - obs_exprs = observed(affect) - if !affect.skip_checks - for oexpr in obs_exprs - invalid_vars = invalid_variables(sys, oexpr) - if length(invalid_vars) > 0 - error("Observed equation $(oexpr) in affect refers to missing variable(s) $(invalid_vars); the variables may not have been added (e.g. if a component is missing).") - end - end - end - obs_syms = observed_syms(affect) - obs_syms, obs_exprs = check_dups(obs_syms, obs_exprs) - - mod_exprs = modified(affect) - if !affect.skip_checks - for mexpr in mod_exprs - if !check_assignable(sys, mexpr) - @warn ("Expression $mexpr cannot be assigned to; currently only unknowns and parameters may be updated by an affect.") - end - invalid_vars = unassignable_variables(sys, mexpr) - if length(invalid_vars) > 0 - error("Modified equation $(mexpr) in affect refers to missing variable(s) $(invalid_vars); the variables may not have been added (e.g. if a component is missing) or they may have been reduced away.") - end - end - end - mod_syms = modified_syms(affect) - mod_syms, mod_exprs = check_dups(mod_syms, mod_exprs) - - overlapping_syms = intersect(mod_syms, obs_syms) - if length(overlapping_syms) > 0 && !affect.skip_checks - @warn "The symbols $overlapping_syms are declared as both observed and modified; this is a code smell because it becomes easy to confuse them and assign/not assign a value." - end - - # sanity checks done! now build the data and update function for observed values - mkzero(sz) = - if sz === () - 0.0 - else - zeros(sz) - end - obs_fun = build_explicit_observed_function( - sys, Symbolics.scalarize.(obs_exprs); - array_type = Tuple) - obs_sym_tuple = (obs_syms...,) - - # okay so now to generate the stuff to assign it back into the system - mod_pairs = mod_exprs .=> mod_syms - mod_names = (mod_syms...,) - mod_og_val_fun = build_explicit_observed_function( - sys, Symbolics.scalarize.(first.(mod_pairs)); - array_type = Tuple) - - 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) - # 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) - - # update the observed values - obs_component_array = NamedTuple{obs_sym_tuple}(obs_fun( - integ.u, integ.p, integ.t)) - - # let the user do their thing - upd_vals = user_affect(upd_component_array, obs_component_array, ctx, integ) - - # 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 - -function compile_affect( - affect::Union{FunctionalAffect, ImperativeAffect}, cb, sys, dvs, ps; kwargs...) +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 diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index 4040b7a646..2b0bd8c8d7 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -629,8 +629,6 @@ function build_explicit_observed_function(sys, ts; oop_mtkp_wrapper = mtkparams_wrapper end - output_expr = isscalar ? ts[1] : - (array_type <: Vector ? MakeArray(ts, output_type) : MakeTuple(ts)) # Need to keep old method of building the function since it uses `output_type`, # which can't be provided to `build_function` return_value = if isscalar diff --git a/src/systems/imperative_affect.jl b/src/systems/imperative_affect.jl new file mode 100644 index 0000000000..f20ad941da --- /dev/null +++ b/src/systems/imperative_affect.jl @@ -0,0 +1,220 @@ + +""" + ImperativeAffect(f::Function; modified::NamedTuple, observed::NamedTuple, ctx) + +`ImperativeAffect` is a helper for writing affect functions that will compute observed values and +ensure that modified values are correctly written back into the system. The affect function `f` needs to have +the signature + +``` + f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple +``` + +The function `f` will be called with `observed` and `modified` `NamedTuple`s that are derived from their respective `NamedTuple` definitions. +Each declaration`NamedTuple` should map an expression to a symbol; for example if we pass `observed=(; x = a + b)` this will alias the result of executing `a+b` in the system as `x` +so the value of `a + b` will be accessible as `observed.x` in `f`. `modified` currently restricts symbolic expressions to only bare variables, so only tuples of the form +`(; x = y)` or `(; x)` (which aliases `x` as itself) are allowed. + +The argument NamedTuples (for instance `(;x=y)`) will be populated with the declared values on function entry; if we require `(;x=y)` in `observed` and `y=2`, for example, +then the NamedTuple `(;x=2)` will be passed as `observed` to the affect function `f`. + +The NamedTuple returned from `f` includes the values to be written back to the system after `f` returns. For example, if we want to update the value of `x` to be the result of `x + y` we could write + + ImperativeAffect(observed=(; x_plus_y = x + y), modified=(; x)) do m, o + @set! m.x = o.x_plus_y + end + +Where we use Setfield to copy the tuple `m` with a new value for `x`, then return the modified value of `m`. All values updated by the tuple must have names originally declared in +`modified`; a runtime error will be produced if a value is written that does not appear in `modified`. The user can dynamically decide not to write a value back by not including it +in the returned tuple, in which case the associated field will not be updated. +""" +@kwdef struct ImperativeAffect + f::Any + obs::Vector + obs_syms::Vector{Symbol} + modified::Vector + mod_syms::Vector{Symbol} + ctx::Any + skip_checks::Bool +end + +function ImperativeAffect(f::Function; + observed::NamedTuple = NamedTuple{()}(()), + modified::NamedTuple = NamedTuple{()}(()), + ctx = nothing, + skip_checks = false) + ImperativeAffect(f, + collect(values(observed)), collect(keys(observed)), + collect(values(modified)), collect(keys(modified)), + ctx, skip_checks) +end +function ImperativeAffect(f::Function, modified::NamedTuple; + observed::NamedTuple = NamedTuple{()}(()), ctx = nothing, skip_checks = false) + ImperativeAffect( + f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) +end +function ImperativeAffect( + f::Function, modified::NamedTuple, observed::NamedTuple; ctx = nothing, skip_checks = false) + ImperativeAffect( + f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) +end +function ImperativeAffect( + f::Function, modified::NamedTuple, observed::NamedTuple, ctx; skip_checks = false) + ImperativeAffect( + f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) +end + +function Base.show(io::IO, mfa::ImperativeAffect) + obs_vals = join(map((ob, nm) -> "$ob => $nm", mfa.obs, mfa.obs_syms), ", ") + mod_vals = join(map((md, nm) -> "$md => $nm", mfa.modified, mfa.mod_syms), ", ") + affect = mfa.f + print(io, + "ImperativeAffect(observed: [$obs_vals], modified: [$mod_vals], affect:$affect)") +end +func(f::ImperativeAffect) = f.f +context(a::ImperativeAffect) = a.ctx +observed(a::ImperativeAffect) = a.obs +observed_syms(a::ImperativeAffect) = a.obs_syms +discretes(a::ImperativeAffect) = filter(ModelingToolkit.isparameter, a.modified) +modified(a::ImperativeAffect) = a.modified +modified_syms(a::ImperativeAffect) = a.mod_syms + +function Base.:(==)(a1::ImperativeAffect, a2::ImperativeAffect) + isequal(a1.f, a2.f) && isequal(a1.obs, a2.obs) && isequal(a1.modified, a2.modified) && + isequal(a1.obs_syms, a2.obs_syms) && isequal(a1.mod_syms, a2.mod_syms) && + isequal(a1.ctx, a2.ctx) +end + +function Base.hash(a::ImperativeAffect, s::UInt) + s = hash(a.f, s) + s = hash(a.obs, s) + s = hash(a.obs_syms, s) + s = hash(a.modified, s) + s = hash(a.mod_syms, s) + hash(a.ctx, s) +end + + +namespace_affects(af::ImperativeAffect, s) = namespace_affect(af, s) +function namespace_affect(affect::ImperativeAffect, s) + ImperativeAffect(func(affect), + namespace_expr.(observed(affect), (s,)), + observed_syms(affect), + renamespace.((s,), modified(affect)), + modified_syms(affect), + context(affect), + affect.skip_checks) +end + +function compile_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs...) + compile_user_affect(affect, cb, sys, dvs, ps; kwargs...) +end + +function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs...) + #= + Implementation sketch: + generate observed function (oop), should save to a component array under obs_syms + do the same stuff as the normal FA for pars_syms + call the affect method + unpack and apply the resulting values + =# + function check_dups(syms, exprs) # = (syms_dedup, exprs_dedup) + seen = Set{Symbol}() + syms_dedup = [] + exprs_dedup = [] + for (sym, exp) in Iterators.zip(syms, exprs) + if !in(sym, seen) + push!(syms_dedup, sym) + push!(exprs_dedup, exp) + push!(seen, sym) + elseif !affect.skip_checks + @warn "Expression $(expr) is aliased as $sym, which has already been used. The first definition will be used." + end + end + return (syms_dedup, exprs_dedup) + end + + obs_exprs = observed(affect) + if !affect.skip_checks + for oexpr in obs_exprs + invalid_vars = invalid_variables(sys, oexpr) + if length(invalid_vars) > 0 + error("Observed equation $(oexpr) in affect refers to missing variable(s) $(invalid_vars); the variables may not have been added (e.g. if a component is missing).") + end + end + end + obs_syms = observed_syms(affect) + obs_syms, obs_exprs = check_dups(obs_syms, obs_exprs) + + mod_exprs = modified(affect) + if !affect.skip_checks + for mexpr in mod_exprs + if !check_assignable(sys, mexpr) + @warn ("Expression $mexpr cannot be assigned to; currently only unknowns and parameters may be updated by an affect.") + end + invalid_vars = unassignable_variables(sys, mexpr) + if length(invalid_vars) > 0 + error("Modified equation $(mexpr) in affect refers to missing variable(s) $(invalid_vars); the variables may not have been added (e.g. if a component is missing) or they may have been reduced away.") + end + end + end + mod_syms = modified_syms(affect) + mod_syms, mod_exprs = check_dups(mod_syms, mod_exprs) + + overlapping_syms = intersect(mod_syms, obs_syms) + if length(overlapping_syms) > 0 && !affect.skip_checks + @warn "The symbols $overlapping_syms are declared as both observed and modified; this is a code smell because it becomes easy to confuse them and assign/not assign a value." + end + + # sanity checks done! now build the data and update function for observed values + mkzero(sz) = + if sz === () + 0.0 + else + zeros(sz) + end + obs_fun = build_explicit_observed_function( + sys, Symbolics.scalarize.(obs_exprs); + array_type = Tuple) + obs_sym_tuple = (obs_syms...,) + + # okay so now to generate the stuff to assign it back into the system + mod_pairs = mod_exprs .=> mod_syms + mod_names = (mod_syms...,) + mod_og_val_fun = build_explicit_observed_function( + sys, Symbolics.scalarize.(first.(mod_pairs)); + array_type = Tuple) + + 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) + # 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) + + # update the observed values + obs_component_array = NamedTuple{obs_sym_tuple}(obs_fun( + integ.u, integ.p, integ.t)) + + # let the user do their thing + upd_vals = user_affect(upd_component_array, obs_component_array, ctx, integ) + + # 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 + + +scalarize_affects(affects::ImperativeAffect) = affects From 1c78dee591ac5623a240084005545dbca04e2941 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Mon, 9 Dec 2024 03:16:24 -0800 Subject: [PATCH 44/47] Fix tests & update API usage --- src/systems/callbacks.jl | 2 +- src/systems/imperative_affect.jl | 4 ++-- test/symbolic_events.jl | 8 ++++---- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index 4533534029..57db5e097c 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -818,7 +818,7 @@ function generate_vector_rootfinding_callback( let save_idxs = save_idxs custom_init = fn.initialize (i) -> begin - isnothing(custom_init) && custom_init(i) + !isnothing(custom_init) && custom_init(i) for idx in save_idxs SciMLBase.save_discretes!(i, idx) end diff --git a/src/systems/imperative_affect.jl b/src/systems/imperative_affect.jl index f20ad941da..19f7f7590d 100644 --- a/src/systems/imperative_affect.jl +++ b/src/systems/imperative_affect.jl @@ -175,7 +175,7 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs. end obs_fun = build_explicit_observed_function( sys, Symbolics.scalarize.(obs_exprs); - array_type = Tuple) + mkarray = (es,_) -> MakeTuple(es)) obs_sym_tuple = (obs_syms...,) # okay so now to generate the stuff to assign it back into the system @@ -183,7 +183,7 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs. mod_names = (mod_syms...,) mod_og_val_fun = build_explicit_observed_function( sys, Symbolics.scalarize.(first.(mod_pairs)); - array_type = Tuple) + mkarray = (es,_) -> MakeTuple(es)) upd_funs = NamedTuple{mod_names}((setu.((sys,), first.(mod_pairs))...,)) diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index 717d2438ac..858a30f4fd 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -936,7 +936,7 @@ end @named trigsys = ODESystem(eqs, t; continuous_events = [evt1, evt2]) trigsys_ss = structural_simplify(trigsys) prob = ODEProblem(trigsys_ss, [], (0.0, 2π)) - sol = solve(prob, Tsit5()) + sol = solve(prob, Tsit5(); dtmax=0.01) required_crossings_c1 = [π / 2, 3 * π / 2] required_crossings_c2 = [π / 6, π / 2, 5 * π / 6, 7 * π / 6, 3 * π / 2, 11 * π / 6] @test maximum(abs.(first.(cr1) .- required_crossings_c1)) < 1e-4 @@ -1079,8 +1079,8 @@ end @test sort(canonicalize(Discrete(), prob.p)[1]) == [0.0, 1.0, 2.0] sol = solve(prob, Tsit5()) - @test sol[a] == [-1.0] - @test sol[b] == [5.0, 5.0] + @test sol[a] == [1.0,-1.0] + @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 @@ -1248,7 +1248,7 @@ end ss = structural_simplify(sys) prob = ODEProblem(ss, [theta => 0.0], (0.0, pi)) sol = solve(prob, Tsit5(); dtmax = 0.01) - @test getp(sol, cnt)(sol) == 197 # we get 2 pulses per phase cycle (cos 0 crossing) and we go to 100 cycles; we miss a few due to the initial state + @test getp(sol, cnt)(sol) == 198 # we get 2 pulses per phase cycle (cos 0 crossing) and we go to 100 cycles; we miss a few due to the initial state end @testset "Initialization" begin From aa556d612a12f2cc297b7af7282a959aa8674775 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Mon, 9 Dec 2024 03:59:31 -0800 Subject: [PATCH 45/47] Adjust quadrature test forward a little to avoid numerical issues --- test/symbolic_events.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index 858a30f4fd..1138b7c96f 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -1246,7 +1246,7 @@ end @named sys = ODESystem( eqs, t, [theta, omega], params; continuous_events = [qAevt, qBevt]) ss = structural_simplify(sys) - prob = ODEProblem(ss, [theta => 0.0], (0.0, pi)) + prob = ODEProblem(ss, [theta => 1e-5], (0.0, pi)) sol = solve(prob, Tsit5(); dtmax = 0.01) @test getp(sol, cnt)(sol) == 198 # we get 2 pulses per phase cycle (cos 0 crossing) and we go to 100 cycles; we miss a few due to the initial state end From c169b9e97629c3031bc9f0e82f3c0899e6e17f71 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Mon, 9 Dec 2024 04:01:00 -0800 Subject: [PATCH 46/47] Formatter --- src/systems/imperative_affect.jl | 6 ++---- test/symbolic_events.jl | 6 +++--- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/src/systems/imperative_affect.jl b/src/systems/imperative_affect.jl index 19f7f7590d..2f489913ad 100644 --- a/src/systems/imperative_affect.jl +++ b/src/systems/imperative_affect.jl @@ -94,7 +94,6 @@ 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) ImperativeAffect(func(affect), @@ -175,7 +174,7 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs. end obs_fun = build_explicit_observed_function( sys, Symbolics.scalarize.(obs_exprs); - mkarray = (es,_) -> MakeTuple(es)) + mkarray = (es, _) -> MakeTuple(es)) obs_sym_tuple = (obs_syms...,) # okay so now to generate the stuff to assign it back into the system @@ -183,7 +182,7 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs. mod_names = (mod_syms...,) mod_og_val_fun = build_explicit_observed_function( sys, Symbolics.scalarize.(first.(mod_pairs)); - mkarray = (es,_) -> MakeTuple(es)) + mkarray = (es, _) -> MakeTuple(es)) upd_funs = NamedTuple{mod_names}((setu.((sys,), first.(mod_pairs))...,)) @@ -216,5 +215,4 @@ function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs. end end - scalarize_affects(affects::ImperativeAffect) = affects diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index 1138b7c96f..ccf0a17a40 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -936,7 +936,7 @@ end @named trigsys = ODESystem(eqs, t; continuous_events = [evt1, evt2]) trigsys_ss = structural_simplify(trigsys) prob = ODEProblem(trigsys_ss, [], (0.0, 2π)) - sol = solve(prob, Tsit5(); dtmax=0.01) + sol = solve(prob, Tsit5(); dtmax = 0.01) required_crossings_c1 = [π / 2, 3 * π / 2] required_crossings_c2 = [π / 6, π / 2, 5 * π / 6, 7 * π / 6, 3 * π / 2, 11 * π / 6] @test maximum(abs.(first.(cr1) .- required_crossings_c1)) < 1e-4 @@ -1079,8 +1079,8 @@ end @test sort(canonicalize(Discrete(), prob.p)[1]) == [0.0, 1.0, 2.0] sol = solve(prob, Tsit5()) - @test sol[a] == [1.0,-1.0] - @test sol[b] == [2.0,5.0, 5.0] + @test sol[a] == [1.0, -1.0] + @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 From b24e5256cd9d3e15b25e9b9468ee2996af6febdc Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Tue, 10 Dec 2024 17:23:24 -0800 Subject: [PATCH 47/47] Update src/systems/callbacks.jl Co-authored-by: Aayush Sabharwal --- src/systems/callbacks.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index 57db5e097c..eaf31a9c5f 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -98,7 +98,7 @@ sharp discontinuity between integrator steps (which in this example would not no guaranteed to be triggered. Once detected the integrator will "wind back" through a root-finding process to identify the point when the condition became active; the method used -is specified by `rootfind` from [`SciMLBase.RootfindOpt`](@ref). If we denote the time when the condition becomes active as `tc``, +is specified by `rootfind` from [`SciMLBase.RootfindOpt`](@ref). If we denote the time when the condition becomes active as `tc`, the value in the integrator after windback will be: * `u[tc-epsilon], p[tc-epsilon], tc` if `LeftRootFind` is used, * `u[tc+epsilon], p[tc+epsilon], tc` if `RightRootFind` is used,