add docs and tests on Discrete/PresetTime CB

Author: hexaeder

docs/src/API.md

@@ -130,6 +130,8 @@ set_bounds!
 NetworkDynamics.ComponentCallback
 ContinousComponentCallback
 VectorContinousComponentCallback
+DiscreteComponentCallback
+PresetTimeComponentCallback
 ComponentCondition
 ComponentAffect
 SymbolicView

docs/src/callbacks.md

@@ -16,14 +16,18 @@ refer to the [Cascading Failure](@ref) example.
 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 automaticially combine them into performant
-[`VectorContinuousCallback`](@extref SciMLBase.VectorContinuousCallback) for the whole network.
+[`VectorContinuousCallback`](@extref SciMLBase.VectorContinuousCallback) and [`DiscreteCallback`](@extref SciMLBase.DiscreteCallback) for the whole network.
 
-The main entry points are the types [`ContinousComponentCallback`](@ref) and
-[`VectorContinousComponentCallback`](@ref). Both of those objects combine a [`ComponentCondition`](@ref)
-with an [`ComponentAffect`](@ref).
-The "normal" `ContinousComponentCallback` has a condition which returns a single value. The corresponding affect is triggered when the return value hits zero.
+The main entry points are the types [`ContinousComponentCallback`](@ref),
+[`VectorContinousComponentCallback`](@ref) and [`DiscreteComponentCallback`](@ref). All of those objects combine a [`ComponentCondition`](@ref) with an [`ComponentAffect`](@ref).
+
+The "normal" [`ContinousComponentCallback`](@ref) and [`DiscreteComponentCallback`](@ref) 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 with an additional argument `event_idx` which tells the effect which dimension encountered the zerocrossing.
 
+There is a special type [`PresetTimeComponentCallback`](@ref) which has no explicit condition and triggers the affect at given times.
+This internally generates a [`PresetTimeCallback`](@ref DiffEqCallbacks.PresetTimeCallback) object from `DiffEqCallbacks.jl`.
+
+
 ### Defining the Callback
 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:
 ```julia

src/NetworkDynamics.jl

@@ -21,7 +21,7 @@ using ForwardDiff: ForwardDiff
 using Printf: @sprintf
 using Random: Random
 using Static: Static, StaticInt
-using SciMLBase: VectorContinuousCallback, CallbackSet
+using SciMLBase: VectorContinuousCallback, CallbackSet, DiscreteCallback
 using DiffEqCallbacks: DiffEqCallbacks
 
 @static if VERSION ≥ v"1.11.0-0"

src/callbacks.jl

@@ -19,8 +19,9 @@ abstract type ComponentCallback end
     ComponentCondition(f::Function, sym, psym)
 
 Creates a callback condition for a [`ComponentCallback`].
-- `f`: The condition function. Must be a function of the form `out=f(u, p, t)` when used
-  for [`ContinousComponentCallback`](@ref) or `f!(out, u, p, t)` when used for
+- `f`: The condition function. Must be a function of the form `out=f(u, p, t)`
+  when used for [`ContinousComponentCallback`](@ref) or
+  [`DiscreteComponentcallback`](@ref) and `f!(out, u, p, t)` when used for
   [`VectorContinousComponentCallback`](@ref).
   - Arguments of `f`
     - `u`: The current value of the selecte `sym` states, provided as a [`SymbolicView`](@ref) object.
@@ -130,7 +131,7 @@ The `affect` will be triggered with the additional `event_idx` argument to know
 dimension the zerocrossing was detected.
 
 The `kwargs` will be forwarded to the `VectorContinuousCallback` when the component based
-callbacks are collected for the whole network using `get_callbacks`.
+callbacks are collected for the whole network using [`get_callbacks(::Network)`](@ref).
 [`DiffEq.jl docs`](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/)
 for available options.
 """
@@ -147,6 +148,21 @@ function VectorContinousComponentCallback(condition, affect, len; kwargs...)
     VectorContinousComponentCallback(condition, affect, len, NamedTuple(kwargs))
 end
 
+"""
+    DiscreteComponentCallback(condition, affect; kwargs...)
+
+Connect a [`ComponentCondition`](@ref) and a [`ComponentAffect`)[@ref] to a
+discrete callback which can be attached to a component model using
+[`add_callback!`](@ref) or [`set_callback!`](@ref).
+
+Note that the `condition` function returns a boolean value, as the discrete
+callback perform no rootfinding.
+
+The `kwargs` will be forwarded to the `DiscreteCallback` when the component based
+callbacks are collected for the whole network using [`get_callbacks(::Network)`](@ref).
+[`DiffEq.jl docs`](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/)
+for available options.
+"""
 struct DiscreteComponentCallback{C<:ComponentCondition,A<:ComponentAffect} <: ComponentCallback
     condition::C
     affect::A
@@ -156,6 +172,20 @@ function DiscreteComponentCallback(condition, affect; kwargs...)
     DiscreteComponentCallback(condition, affect, NamedTuple(kwargs))
 end
 
+"""
+    PresetTimeComponentCallback(ts, affect; kwargs...)
+
+Tirgger a [`ComponentAffect`](@ref) at given timesteps `ts` in discrete
+callback, which can be attached to a component model using
+[`add_callback!`](@ref) or [`set_callback!`](@ref).
+
+The `kwargs` will be forwarded to the [`PresetTimeCallback`](@ref DiffEqCallbacks.PresetTimeCallback)
+when the component based callbacks are collected for the whole network using
+[`get_callbacks(::Network)`](@ref).
+
+The `PresetTimeCallback` will take care of adding the timesteps to the solver, ensuring to
+exactly trigger at the correct times.
+"""
 struct PresetTimeComponentCallback{T,A} <: ComponentCallback
     ts::T
     affect::A
@@ -417,12 +447,12 @@ end
 function _batch_condition(dcw::DiscreteCallbackWrapper)
     uidxtype = dcw.component isa EIndex ? EIndex : VIndex
     pidxtype = dcw.component isa EIndex ? EPIndex : VPIndex
-    usymidxs = uidxtype(dcw.component.compidx, dcw.callback.cond.sym)
-    psymidxs = pidxtype(dcw.component.compidx, dcw.callback.cond.psym)
+    usymidxs = uidxtype(dcw.component.compidx, dcw.callback.condition.sym)
+    psymidxs = pidxtype(dcw.component.compidx, dcw.callback.condition.psym)
     ucache = DiffCache(zeros(length(usymidxs)), 12)
 
-    obsf = SII.observed(ccw.nw, usymidxs)
-    pidxs = SII.parameter_index.(Ref(ccw.nw), psymidxs)
+    obsf = SII.observed(dcw.nw, usymidxs)
+    pidxs = SII.parameter_index.(Ref(dcw.nw), psymidxs)
 
     (u, t, integrator) -> begin
         us = PreallocationTools.get_tmp(ucache, u)

test/callbacks_test.jl

@@ -76,14 +76,29 @@ end
     b = @b $batchcond($out, $u, NaN, $integrator)
     @test b.allocs == 0
 
+    # test the preste time callback
     tripfirst = PresetTimeComponentCallback(1.0, affect) # reuse the same affect
     add_callback!(nw[EIndex(5)], tripfirst)
 
-    nwcb = NetworkDynamics.get_callbacks(nw)
+    # add a useless discrete callback
+    useless_triggertime = Ref{Float64}(0.0)
+    usless_cond = ComponentCondition([:P, :₋P, :srcθ], [:limit, :K]) do u, p, t
+        t > 0.1 && iszero(useless_triggertime[])
+    end
+    usless_affect = ComponentAffect([], [:limit, :K]) do u, p, ctx
+        @info "Usless effect triggered at $(ctx.t)"
+        useless_triggertime[] = ctx.t
+    end
+    useless_cb = DiscreteComponentCallback(usless_cond, usless_affect)
+    add_callback!(nw[EIndex(1)], useless_cb)
+
+    nwcb = NetworkDynamics.get_callbacks(nw);
     s0 = NWState(nw)
     prob = ODEProblem(nw, uflat(s0), (0,6), copy(pflat(s0)), callback=nwcb)
     sol = solve(prob, Tsit5());
 
+    @test 0.1 < useless_triggertime[] <= 1.0
+
     @test tripi == [5,7,4,1,3,2]
     tref = [1, 2.247676397005474, 2.502523192233235, 3.1947647115093654, 3.3380530127462587, 3.4042696241577888]
     @test maximum(abs.(tript - tref)) < 1e-5