Callbacks and Events
Callback-functions are a way of handling discontinuities in differential equations. In a nutshell, the solver checks for some "condition" (i.e. a zero crossing of some variable) and calls some "affect" if the condition is fulfilled. Within the affect function, it is safe to modify the integrator, e.g. changing some state or some parameter.
Since NetworkDynamics.jl provides nothing more than a RHS for DifferentialEquations.jl, please check their docs on event handling as a general reference. This page is introducing the general concepts, for a hands on example of a simulation with callbacks refer to the Cascading Failure example.
The ODEProblem contains a reference to exactly one copy of the flat parameter array. If you use callbacks to change those parameters (as we often do), it is advised to copy the parameter array before passing it to the ODEProblem! Also, this means you need to be careful when using the same prob for multiple subsequent solve calls, as the initial state of the prob object might have changed!
Component-based Callback functions
In practice, events often act locally, meaning they only depend and act on a specific component or type of component. NetworkDynamics.jl provides a way of defining those callbacks on a component level and automatically combine them into performant VectorContinuousCallback and DiscreteCallback for the whole network.
The main entry points are the types ContinuousComponentCallback, VectorContinuousComponentCallback and DiscreteComponentCallback. All of those objects combine a ComponentCondition with an ComponentAffect.
The "normal" ContinuousComponentCallback and DiscreteComponentCallback have a condition which returns a single value. The corresponding affect is triggered when the return value hits zero. In contrast, the "vector" version has an in-place condition which writes len outputs. When any of those outputs hits zero, the affect is triggered once with an additional argument event_signs: a length-len vector where entry i is 0 (no crossing), +1 (upcrossing) or -1 (downcrossing) for the i-th output. The affect resolves the crossing direction (and any simultaneous crossings) itself. This mirrors the affect!(integrator, event_signs) interface of the underlying VectorContinuousCallback.
There is a special type PresetTimeComponentCallback which has no explicit condition and triggers the affect at given times. This internally generates a PresetTimeCallback object from DiffEqCallbacks.jl.
Defining the Callback
To construct a condition function, you need to tell NetworkDynamics which symbols of the component you'd like to "observe" within the condition. Any named symbol of the component model works: states, parameters, inputs, outputs and observed. Within the actual condition, those values are made available through u:
condition = ComponentCondition([:x, :y, :p1]) do u, t
u[:x] == u[1] # access a state or observable :x at current time
u[:p1] == u[3] # parameters are listed like any other symbol
return some_condition(u[:x], u[:y], u[:p1])
endIn case of a VectorContinuousComponentCallback, the function signature looks slightly different:
vectorcondition = ComponentCondition([:x, :y, :p1]) do out, u, t
out[1] = some_condition(u[...])
out[2] = some_condition(u[...])
return nothing
endThe argument u will be passed as a SymbolicView object, which means it is possible to use the getindex syntax to access the desired values by name.
The affect takes the same kind of symbol list:
affect = ComponentAffect([:u, :p, :obs]) do u, ctx
t = ctx.t # extract data from context
u[:u] = u[:obs] # states and parameters are writable, observed are readable
u[:p] = 0
println("Trigger affect at t=$t")
end
vectoraffect = ComponentAffect([:u, :p]) do u, event_signs, ctx
for i in eachindex(event_signs)
event_signs[i] == 0 && continue # skip outputs that did not cross
if i == 1
u[:u] = 0 # first output crossed: change state
else
u[:p] = 0 # second output crossed: change parameter
end
# event_signs[i] is +1 for an upcrossing and -1 for a downcrossing
end
endEntries of u which are states or parameters of the component can be written to, all other entries (inputs, outputs, observed) are read only. The values in u are a snapshot taken when the affect fires: writing a state or parameter updates the integrator immediately, but an observed entry which depends on it is not recomputed within the same affect.
The affect gets passed a ctx "context" object, which is a named tuple which holds additional context like the integrator object, the component model, the index of the component model, the current time and so on. Please refer to the ComponentAffect docstring for a detailed list.
Earlier versions took two symbol lists, ComponentCondition(f, sym, psym) with f(u, p, t) and ComponentAffect(f, sym, psym) with f(u, p, ctx), where p gave access to the parameters. This form is still accepted and behaves as before.
Lastly we need to define the actual callback object using ContinuousComponentCallback/VectorContinuousComponentCallback:
ccb = ContinuousComponentCallback(condition, affect; kwargs...)
vccb = VectorContinuousComponentCallback(condition, affect, len; kwargs...)where the kwargs are passed to the underlying SciMLBase.VectorContinuousCallback to finetune the zerocrossing-detection.
Registering the Callback
Once the callback is defined, we need to "attach" it to the component, for that you can use the methods add_callback! and set_callback!:
vert = VertexModel(...)
add_callback!(vert, ccb)
add_callback!(vert, vccb)Extracting the Callback
Component-level callbacks are automatically extracted and combined when constructing an ODEProblem:
u0 = NWState(nw)
prob = ODEProblem(nw, u0, (0, 10))
sol = solve(prob, ...)For more control over callback handling—such as adding network/system-level callbacks (e.g., PeriodicCallback), temporary component callbacks, or overriding the default callbacks, the ODEProblem constructor provides keyword arguments add_comp_cb, add_nw_cb, and override_cb. See the ODEProblem(nw::Network,...) documentation for details.
When executing component callbacks, NetworkDynamics automatically checks whether states or parameters changed during the affect and calls SciMLBase.auto_dt_reset! and save_parameters! if necessary. An affect which only does bookkeeping, for example counting events in a parameter, can keep the current step size by setting ctx.dt_reset[] = false. The parameter change is saved either way.
Event Iteration
Discrete logic such as hysteresis switches, timers or sample-and-hold blocks keeps its memory in parameters which callbacks set. Such blocks are meant to act synchronously. At one event instant they should all see the same state, and one block's switch may cause another block to switch at the same instant. Plain discrete callbacks don't work like that: separate callbacks run one after another, so the result depends on their order, and a cascade within one instant is only noticed in the next step.
For this, a discrete callback can join the event iteration of the network:
cb = DiscreteComponentCallback(condition, affect; iterative=true)All iterative callbacks of the network run together as one set, after every other callback at an event instant. They run in rounds:
- All conditions are evaluated on the same state.
- The affects of all fired callbacks run. They all read one snapshot, taken before the first of them writes.
- If the network is a DAE, its algebraic states are reinitialized with the
initializealgof the solve. Then the next round starts at step 1.
The loop ends as soon as no condition fires anymore. Step size reset and parameter saving happen once for the whole instant.
Because the conditions are evaluated several times within one instant, they should be predicates on the current state alone, for example "the stored switch state contradicts the input". The affect which fixes the inconsistency then also makes its own condition false. Iterative affects should read and write only through u. Writes through ctx.integrator are invisible to the iteration, so they trigger no reinit and no further round.
The order of the callbacks at an event instant is:
- the continuous callback which found the earliest root, if any,
- the preset-time callbacks,
- the other discrete callbacks,
- the iterative set.
The callbacks passed as add_nw_cb to the ODEProblem are sorted into the same groups. Within a group, the network's own callbacks come first.
So the event iteration reacts to the jumps of all other callbacks. This holds when the callbacks come from the ODEProblem constructor. If you combine get_callbacks(nw) with your own callbacks in a CallbackSet by hand, your discrete callbacks run after the iterative set. If the set still fires after event_maxiter rounds (default 10), it keeps the state of the last round and warns, naming the components which still fire. Pass event_failure=:error to throw instead. Both are keyword arguments of get_callbacks and of the ODEProblem constructor.
Normal DiffEq Callbacks
Besides component based callbacks, it is also possible to use "normal" DiffEq callbacks together with NetworkDynamics.jl. It is far more powerful but also more cumbersome compared to the component based callback functions. To access states and parameters of specific components, we heavily rely on the Symbolic Indexing features.
using SymbolicIndexingInterface as SII
nw = Network(#= some network =#)
condition = let getvalue = SII.getsym(nw, VIndex(1:5, :some_state))
function(out, u, t, integrator)
s = NWState(integrator, u, integrator.p, t)
some_state = getvalue(s)
out .= some_condition(some_state)
end
endPlease note a few important things here:
- Symbolic indexing can be costly, and the condition function gets called very often. By using
SII.getsymwe did some of the work before the callback by creating the accessor function. When handling with "normal states" and parameters consider usingSII.variable_indexandSII.parameter_indexfor even better access patterns. trefers to the current time of the zerocrossing-detection-algorithm. This is different fromintegrator.twhich refers to the current timestep in which the zerocross-detectio takes place..
function affect!(integrator, vidx)
p = NWParameter(integrator) # get symbolically indexable parameter object
p.v[vidx, :some_vertex_parameter] = 0 # change some parameter
auto_dt_reset!(integrator)
save_parameters!(integrator)
endThe affect function is much more straight forward, as it (typically) is called far less frequent and thus less perfomance critical.
Once the condition and affect! is defined, you can use the SciMLBase.ContinuousCallback and SciMLBase.VectorContinuousCallback constructors to create the callback.
Since changes to u and p mostly introduce discontinuities in the solution, it is recommend to call auto_dt_reset! within the affect to restart integration with small steps afterwards.
An "observable" is kind of a "virtual" state, which can be reconstructed for a given time t, a given state u and a given set of parameters p
\[o = f(u(t), p(t), t)\]
To extract or plot timeseries of observed states under time variant parameters (i.e. parameters that are changed in a callback), those changes need to be recorded using the save_parameters! function whenever p is changed. When using ComponentCallback, NetworkDynamics will automatically check for changes in p and save them if necessary.