WIP Add discrete callbacks

[isaacsas]

This is a followup to #1661 that adds discrete callbacks.

TODO

• Simple test.
• Rebase on master when Cleanup callbacks  #1661 is merged.
• Needs wrap_code for scalar functions to enable DiscreteCallbacks in ModelingToolkit JuliaSymbolics/Symbolics.jl#632 (could possibly support older Symbolics versions by wrapping the condition function in a second function as done for continuous callbacks...)
• More ODE tests
• Add other system types.

[baggepinnen]

Any overlap with

• Periodic Callbacks with general affect function (PoC) #1649
  ?

[isaacsas]

I don't really understand the approach in that PR. This is more like the previous continuous events PR, fully translating the symbolic callback into the standard DiscreteCallback form of DifferentialEquations.jl.

[isaacsas]

I do think it would be better to build future callbacks off #1661 which collates the various callback functionality into one location, making it easier to see what is needed...

[dodoplus]

Any overlap with

• Periodic Callbacks with general affect function (PoC) #1649
  ?

The point of that was less about the callbacks themselves and more about a more general affect! functions. PeriodicCallbak was more of a placeholder (the condition part is really simple in this case).

I want to simulate a system with a complex control mechanism. As a motivating example, think MPPT inverter. The logic can be quite complex and cannot be defined by a set of equations in general.
As another example, think about "memory" effect (e.g. hysteresis, battery charging)

An affect! has the following signature:
affect!(u, p, t, state)

where t is time, u and p are the dependent and independent (parameters) variables, respectively and state is some state provided when constructing the system.

To add a periodic callback you do:
ODESystem(eqs, t, ... , periodic_events=[(dt, affect!, state),...])

which I don't like that much.

within the affect! function you can access the parameters/states of the ODESystem you defined the event for by name:

function affect!(u, p, t, state)
    pars.dir = some_logic(state, pars.dir, u.frate) # update the flow direction based on some external logic
    if Int(t) % 2 == 1 # let's assume dt = 1.0
        p[:G] *= 2 # double the gravitational constant. Note that p[:G] === p.G
    end
end

I'm not sure this is the right interface for a generalized affect!, but I think we need more flexibility than is afforded by the current callbacks.

Feel free to merge the PIs

DD

[isaacsas]

Allowing more general non-symbolic affect functions, that can in someway use the symbolic names in being written, would definitely be very helpful. It is a general MT issue to make it easier to interface compiled code for the DifferentialEquations.jl interface with MT (i.e. going between indices and the MT variable/parameter names).

This PR is more about getting a very simple notation for common discrete callback cases that can be described symbolically.

I've had several people ask how to change a parameter at a set time in Catalyst the past several months, and SBMLToolkit would also like to be able to use symbolic equations for discrete events, so it would be nice to easily support those use cases for discrete callbacks via letting users just write a few simple symbolic expressions.

[baggepinnen]

I want to simulate a system with a complex control mechanism. As a motivating example, think MPPT inverter. The logic can be quite complex and cannot be defined by a set of equations in general.

I agree with this sentiment, we absolutely have to enable this sooner or later. Now might thus be a good time to think hard about the interface, in case the interface we have with the symbolic events have to change in a breaking way. Maybe all events can be consolidated into a sufficiently general interface (split up into continuous/discrete just like in DiffEq of course).

[isaacsas]

I don't really see how this is an either/or situation. Having a simple interface for callbacks that can be easily defined symbolically seems very useful. In the chemical modeling context what multiple users have asked me for is very basic -- just having the ability to easily change a parameter during a simulation at a set time. This is something that is trivial to define symbolically for users.

That in no way precludes adding future support for more complicated events that involve compiled functions. They can subsequently be handled by dispatches or just different keywords... At the end of the day all the DifferentialEquations.jl callbacks that get generated will be merged into one CallbackSet so there is no reason we can't support a variety of ways to input them in ModelingToolkit.

I don't think we should push simple one line events that are adequately described symbolically to have to be written as a compiled function.

[dodoplus]

I agree with this sentiment, we absolutely have to enable this sooner or later. Now might thus be a good time to think hard about the interface, in case the interface we have with the symbolic events have to change in a breaking way. Maybe all events can be consolidated into a sufficiently general interface (split up into continuous/discrete just like in DiffEq of course).

Events are defined by their condition and affect!.
For continuous events, an equation based condition makes a lot of sense. However, as DiffEq shows, specialized conditions are useful for some discrete events.
It seems like this PI keeps the equation-based condition for discrete events, but I think this should be generalized to encapsulate event types that DiffEq already identifies: PeriodicCallback (as I've started to do), PresetTimeCallback, and IterativeCallback.

One can simply pass a Real (periodic, iterative) or a Vector (preset-time).

I think all of these event types should support affect equations: in the case of a iterative events, this would be done by setting a special variable (say, τ) in one of the equations (if it is not on the lhs of an equation, the event periodic). Stopping an iterative event is a bit of a challenge in this case, but returning Inf seems reasonable.

Above (and in the PI) I've described what I think is a reasonable interface for a more general affect!, but I'm open for any improvements. The difference between iterative and periodic affect! functions could be their return values: returning nothing repeats the last period, while returning a value will change the period (i.e. its an iterable). Returning Inf will not call the callback again.

I've described above how I think variables can be (symbolically) accessed in the affect!. I'll add that an affect! should only have access to variables in its own namespace (i.e. defined in the same system as the callback or in one of its subsystems). It should not have access to variables higher in the hierarchy. Variables of a subsystem should be accessed hierarchically, i.e. u.subsystem.var.

What do you think?

[isaacsas]

Just to clarify; this doesn't use an Equation based condition for discrete events. It uses a symbolic expression that should evaluate to true or false once compiled, consistent with the current DiscreteCallback interface in DifferentialEquations.jl. A user can specify any valid Symbolics.jl expression here and it should work.

[dodoplus]

Just to clarify; this doesn't use an Equation based condition for discrete events. It uses a symbolic expression that should evaluate to true or false once compiled, consistent with the current DiscreteCallback interface in DifferentialEquations.jl. A user can specify any valid Symbolics.jl expression here and it should work.

Right - I see that in the code.
This is different than ContinuousCallback, though, where you do have a Vector{Equation}.
I assume this is needed to deal with inequalities?
I think ContinuousCallback can also benefit from this. E.g. trigger an event when x==0 but only if y>1

[isaacsas]

This needed as discrete callbacks in DiffEq assume the condition returns true or false, which could include an inequality. Continuous callbacks in DiffEq expect a number and look for that number being zero. Continuous callbacks here could be modified though to not take equations, but just conditions to check if they are zero.

[dodoplus]

What about we add support for periodic and preset-time events (I think that iterative will require generalized affect) as follows:

struct SymbolicDiscreteCallback
    condition
    affects::Vector{Equation}

    function SymbolicDiscreteCallback(condition, affects = NULL_AFFECT)
        c = scalarize_condition(condition) # changed
        a = scalarize(affects)
        new(c, a)
    end # Default affect to nothing
end

function discrete_events(sys::AbstractSystem)
    obs = get_discrete_events(sys)
    systems = get_systems(sys)
    cbs = [obs;
           reduce(vcat,
                  (map(o -> namespace_callback(o, s), discrete_events(s)) for s in systems), # changed
                  init = SymbolicDiscreteCallback[])]
    cbs
end

function generate_discrete_callbacks(sys::AbstractSystem, dvs = states(sys),
                                    ps = parameters(sys); kwargs...)
    has_discrete_events(sys) || return nothing
    symcbs = discrete_events(sys)
    isempty(symcbs) && return nothing

    dbs = map(symcbs) do cb
        generate_discrete_callback(cb, sys, dvs, ps; kwargs...) # changed
    end

    dbs
end

# Added

is_timed_condition(cb) = false
is_timed_condition(::R) where {R<:Real} = true
is_timed_condition(::V) where {V<:AbstractVector} = eltype(V) <: Real
is_timed_condition(cb::SymbolicDiscreteCallback) = is_timed_condition(condition(cb))

scalarize_condition(condition) = is_timed_condition(condition) ? condition : value(scalarize(condition))
namespace_condition(condition, s) = is_timed_condition(condition) ? condition : namespace_expr(condition, s)

function namespace_callback(cb::SymbolicDiscreteCallback, s)::SymbolicDiscreteCallback
    SymbolicDiscreteCallback(namespace_condition(condition(cb), s),
                             namespace_affect.(affects(cb), Ref(s)))
end

function generate_timed_callback(cb, sys, dvs, ps; kwargs...)
    cond = condition(cb)
    as = compile_affect(affects(cb), sys, dvs, ps; expression = Val{false},
                        kwargs...)
    if cond isa AbstractVector
        # Preset Time
        return PresetTimeCallback(cond, as)
    else
        # Periodic
        return PeriodicCallback(as, cond)
    end
end

function generate_discrete_callback(cb, sys, dvs, ps; kwargs...)
    if is_timed_condition(cb)
        return generate_timed_callback(cb, sys, dvs, ps, kwargs...)
    else
        c = compile_condition(cb, sys, dvs, ps; expression=Val{false}, kwargs...)
        as = compile_affect(affects(cb), sys, dvs, ps; expression = Val{false},
                        kwargs...)
        return DiscreteCallback(c, as)
    end
end

A periodic or a preset-time callback can then be specified by passing a period (Real) or times (Vector{Real}) as condition.

[isaacsas]

Wouldn't one want to allow the use of symbolic parameter(s) to specify event times? Then we would suddenly have a confusing intersection of whether a user's parameter is supposed to be a logical condition or an event time.

I would argue we should just parallel the design in DiffEqCallbacks, which I find works quite nicely. This makes it easier on users -- they are simply creating symbolic versions of the existing DiffEqCallbacks.

So we could just have a SymbolicPresetTimeCallback, which wraps building a SymbolicDiscreteCallback (analogous to DiffEqCallbacks), or just directly wraps PresetTimeCallback. Perhaps we should change the design in ModelingToolkit to just take generic AbstractSymbolicCallbacks via a single events or callback kwarg, maybe as a CallbackSet, and then process and merge them. Then a library of symbolic callbacks could be built, paralleling DiffEqCallbacks.

Many DifferentialEquations.jl callbacks can be used with symbolic expressions, so we could potentially even drop the "SymbolicContinuousCallback" and "SymbolicDiscreteCallback" and just let users create DiffEqCallbacks that contain symbolic expressions (processing them via dispatch).

[dodoplus]

Wouldn't one want to allow the use of symbolic parameter(s) to specify event times? Then we would suddenly have a confusing intersection of whether a user's parameter is supposed to be a logical condition or an event time.

(PresetTime|Periodic)Callback accept a single parameter for timing the events. These events are fixed in advance. I don't think they gain from making them symbolic. The 3rd specialized discrete callback – IterativeCallback – on the other hand accepts a time_choice function, which could be replaced with a symbolic condition.

However, I don't think we need a symbolic IterativeCallback as it's just SymbolicDiscreteCallback.

Perhaps we should change the design in ModelingToolkit to just take generic AbstractSymbolicCallbacks via a single events or callback kwarg, maybe as a CallbackSet, and then process and merge them. Then a library of symbolic callbacks could be built, paralleling DiffEqCallbacks.

The current design accepts SymbolicContinuousCallback as a Pair. A CallbackSet interface will require the user to construct the *Callbacks E.g. something like:

ODESystem(eqs, ..., events = [DiscreteEvent([x], [d ~ c]), PresetTimeEvent([2, 3, 5, 8], [x ~ -x])])

I don't like this. Of course, you could distinguish between callback types based on the type of the condition - which is the approach I chose above. Then you have:

ODESystem(eqs, ..., discrete_events = [x  => d ~ c, [2, 3, 5, 8] => x ~ -x])

[isaacsas]

@YingboMa do you mind letting me know if this looks ok to you? If so I'll add support for other systems types and add more tests, but I don't want to invest more time unless we are good with having a discrete event approach that parallels the current continuous symbolic event approach. (I know about the formatting issues and will fix when this is completed.)

[isaacsas]

@YingboMa just a bump about this. @TorkelE and I were thinking of using this as part of the Catalyst Juliacon tutorial in a couple weeks, if it is merged, given the amount of questions we get about these types of callbacks in Catalyst. Would you mind taking a look and letting me know your thoughts (i.e. whether I should finish it off, or if this is not an approach we want in ModelingToolkit)? Thanks!