|
1 | | -# Callbacks |
| 1 | +# Callbacks and Events (@id Callbacks) |
2 | 2 |
|
3 | | -tbd |
| 3 | +Callback-functions are a way of handling discontinuities in differential equations. |
| 4 | +In a nutshell, the solver checks for some "condition" (i.e. a zero crossing of some variable) |
| 5 | +and calls some "affect" if the condition is fulfilled. |
| 6 | +In the affect it is safe to modify the integrator, e.g. changing some state or some parameter. |
| 7 | + |
| 8 | +Since `NetworkDynamics.jl` provides nothing more than a RHS for DifferentialEquations.jl, please check |
| 9 | +[their docs on event handling](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/) |
| 10 | +as a general reference. |
| 11 | + |
| 12 | +This page at introducing the general concepts, for a hands on example of a simulation with callbacks |
| 13 | +refer to the [Cascading Failure](@ref) example. |
| 14 | + |
| 15 | + |
| 16 | +## Component-based Callback functions |
| 17 | +In practice, events often act locally, meaning they only depend and act on a |
| 18 | +specific component or type of component. `NetworkDynamics.jl` provides a way of |
| 19 | +defining those callbacks on a component level and automaticially combine them into performant |
| 20 | +[`VectorContinuousCallback`](@extref SciMLBase.VectorContinuousCallback) for the whole network. |
| 21 | + |
| 22 | +The main entry points are the types [`ContinousComponentCallback`](@ref) and |
| 23 | +[`VectorContinousComponentCallback`](@ref). Both of those objects combine a [`ComponentCondition`](@ref) |
| 24 | +with an [`ComponentAffect`](@ref). |
| 25 | +The "normal" `ContinousComponentCallback` has a condition which returns a single value. The corresponding affect is triggered when the return value hits zero. |
| 26 | +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 with an additional argument `event_idx` which tells the effect which dimension encountered the zerocrossing. |
| 27 | + |
| 28 | +### Defining the Callback |
| 29 | +To construct a condition function, you need to tell network dynamics which states and parameters you'd like to "observe" within the condition. Within the actual condition, those states will be made available: |
| 30 | +```julia |
| 31 | +condition = ComponentCond([:x, :y], [:p1, :p2]) do u, p, t |
| 32 | + u[:x] == u[1] # access a state or observable :x at current time |
| 33 | + p[:p2] == p[2] # access a parameter at current time |
| 34 | + return some_condition(u[:x], u[:y], ...) |
| 35 | +end |
| 36 | +``` |
| 37 | +In case of a `VectorContinousComponentCallback`, the function signature looks slightly different: |
| 38 | +```julia |
| 39 | +vectorcondition = ComponentCond([:x, :y], [:p1, :p2]) do out, u, p, t |
| 40 | + out[1] = some_condition(u[...], p[...]) |
| 41 | + out[2] = some_condition(u[...], p[...]) |
| 42 | + return nothing |
| 43 | +end |
| 44 | +``` |
| 45 | +Note that the `syms` argument (here `[:x, :y]`) can be used to reference **any** |
| 46 | +named state of the component model, this includes "ordinary" states, observed, |
| 47 | +inputs and outputs. |
| 48 | +The arguments `u` and `p` will be passed as [`SymbolicView`](@ref) objects, which mean |
| 49 | +it is possible to use the getindex syntax to acces the desired states by name. |
| 50 | + |
| 51 | +The affect takes a similar form: |
| 52 | +```julia |
| 53 | +affect = ComponentAffect([:u], [:p]) do u, p, ctx |
| 54 | + t = ctx.t # extract data from context |
| 55 | + obs = NWState(ctx.integrator)[VIndex(ctx.vidx, :obs)] # extract some observed state from context |
| 56 | + println("Trigger affect at t=$t") |
| 57 | +end |
| 58 | +vectoraffect = ComponentAffect([:u], [:p]) do u, p, event_idx, ctx |
| 59 | + if event_idx == 1 |
| 60 | + u[:u] = 0 # change state |
| 61 | + else |
| 62 | + u[:p] = 0 # change parameter |
| 63 | + end |
| 64 | +end |
| 65 | +``` |
| 66 | +Notably, the `syms` (here `:u`) can *exclusivly* refer to "ordinary" states, since they are now writable. |
| 67 | +However 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 refere to the [`ComponentAffect`](@ref) docstring for a detailed list. |
| 68 | + |
| 69 | +Lastly we need to define the actuall callback object using [`ContinousComponentCallback`](@ref)/[`VectorContinousComponentCallback`](@ref): |
| 70 | +```julia |
| 71 | +ccb = ContinousComponentCallback(condition, affect; kwargs...) |
| 72 | +vccb = VectorContinousComponentCallback(condition, affect; kwargs...) |
| 73 | +``` |
| 74 | +where the `kwargs` are passed to the underlying [`SciMLBase.VectorContinuousCallback`](@extref) to finetune the zerocrossing-detection. |
| 75 | + |
| 76 | + |
| 77 | +### Registering the Callback |
| 78 | +Once the callback is defined, we need to "attach" it to the component, for that you can use the methods [`add_callback!`](@ref) and [`set_callback!`](@ref): |
| 79 | +```julia |
| 80 | +vert = VertexModel(...) |
| 81 | +add_callback!(vert, ccb) |
| 82 | +add_callback!(vert, vccb) |
| 83 | +``` |
| 84 | + |
| 85 | + |
| 86 | +### Extracting the Callback |
| 87 | +In order to use the callback during simulation, we need to generate a [`SciMLBase.CallbackSet`](@extref) which contains the conditions and affects of all the component based callbacks in the network. For that we use [`get_callbacks(::Network)`](@ref `get_callbacks(::NetworkDynamics.Network)`): |
| 88 | +```julia |
| 89 | +u0 = NWState(u0) |
| 90 | +cbs = get_callbacks(nw) |
| 91 | +prob = ODEProblem(nw, uflat(u0), (0,10), pflat(u0); callback=cbs) |
| 92 | +sol = solve(prob, ...) |
| 93 | +``` |
| 94 | + |
| 95 | +When combining the component based callbacks to a single callback, NetworkDynamics will check whether states and or parameters changed during the affect and automaticially call [`SciMLBase.auto_dt_reset!`](@extref) and [`save_parameters!`](@ref) if necessary. |
| 96 | + |
| 97 | + |
| 98 | +## Normal DiffEq Callbacks |
| 99 | +Besides component based callbacks, it is also possible to use "normal" DiffEq |
| 100 | +callbacks together with `NetworkDynamics.jl`. |
| 101 | +It is far more powerful but also more cumbersome compared to the component based callback functions. |
| 102 | +To access states and parameters of specific components, we havily rely on the [Symbolic Indexing](@ref) features. |
| 103 | + |
| 104 | +```julia |
| 105 | +using SymbolicIndexingInterface as SII |
| 106 | +nw = Network(#= some network =#) |
| 107 | + |
| 108 | +condition = let getvalue = SII.getu(nw, VIndex(1:5, :some_state)) |
| 109 | + function(out, u, t, integrator) |
| 110 | + s = NWState(integrator, u, integrator.p, t) |
| 111 | + some_state = getvalue(s) |
| 112 | + out .= some_condition(some_state) |
| 113 | + end |
| 114 | +end |
| 115 | +``` |
| 116 | +Please not a few important things here: |
| 117 | + - Symbolic indexing can be costly, and the condition function gets called very |
| 118 | + often. By using [`SII.getu`](@extref `SymbolicIndexingInterface.getu`) we did |
| 119 | + some of the work *before* the callback by creating the accessor function. |
| 120 | + When handling with "normal states" and parameters consider using |
| 121 | + [`SII.variable_index`](@ref `SymbolicIndexingInterface.variable_index`) and |
| 122 | + [`SII.parameter_index`](@ref `SymbolicIndexingInterface.parameter_index`) for |
| 123 | + even better access patterns. |
| 124 | + - `t` refers to the current time of the zerocrossing-detection-algorithm. This is different from `integrator.t` which refers to the current timestep in which the zerocross-detectio takes place.. |
| 125 | + |
| 126 | +```julia |
| 127 | +function affect!(integrator, vidx) |
| 128 | + p = NWParameter(integrator) # get symbolicially indexable parameter object |
| 129 | + p.v[vidx, :some_vertex_parameter] = 0 # change some parameter |
| 130 | + auto_dt_reset!(integrator) |
| 131 | + save_parameters!(integrator) |
| 132 | +end |
| 133 | +``` |
| 134 | +The affect function is much more straight forward, as it (typically) is called far less frequent and thus less perfomance critical. |
| 135 | + |
| 136 | +Once the `condition` and `affect!` is defined, you can use the [`SciMLBase.ContinuousCallback`](@extref) and [`SciMLBase.VectorContinuousCallback`](@extref) constructors to create the callback. |
| 137 | + |
| 138 | +!!! note "Introducing discontinuities with adaptive timestepping" |
| 139 | + Since changes to `u` and `p` mostly introduce discontinuities in the |
| 140 | + solution, it is recommend to call [`auto_dt_reset!`](@extref |
| 141 | + `SciMLBase.auto_dt_reset!`) within the affect to restart integration with |
| 142 | + small steps afterwards. |
| 143 | + |
| 144 | +!!! note "Changing Parameters and Observables" |
| 145 | + An "observable" is kind of a "virtual" state, which can be reconstructed for |
| 146 | + a given time `t`, a given state `u` and a given set of parameters `p` |
| 147 | + ```math |
| 148 | + o = f(u(t), p(t), t) |
| 149 | + ``` |
| 150 | + To extract or plot timeseries of observed states under *time variant |
| 151 | + parameters* (i.e. parameters that are changed in a callback), those changes |
| 152 | + need to be recorded using the [`save_parameters!`](@ref) function whenever `p` is changed. |
| 153 | + When using [ComponentCallback](@ref), NetworkDynamics will automaticially check for changes in `p` and save them if necessary. |
0 commit comments