Feature request: support for events & callbacks

[jleugeri]

Hi!

During JuliaCon I mentioned this already in Slack, so I'm now just posting an Issue to make sure you don't forget about my feature request 😄

It would be awesome to have an interface in NetworkDynamics.jl for triggering events based on the state of individual components (nodes or edges), but right now that seems to be a little bit complicated and un-documented

For example for my work, I need to be able to trigger an event when the state of one node crosses above a critical value, which should then affect the state of a different node or edge.
A clean way to achieve that might be to only allow interactions via edges, but then there should be a way to specify a callback that is triggered based on the state of (one of) the edge's connected nodes and can affect the states of (one of) the connected nodes. Or maybe there should be an EventEdge type?

Since DifferentialEquations.jl already offers a specialized VectorContinuousCallback, it'd probably be best if NetworkDynamics.jl provided some convenient way to define individual call-backs on a per-node / per-edge basis, but then aggregated all of the individual continuous callbacks behind the scenes into one vector-valued callback. (Probably similar for other callback types, but I haven't given that much thought.)

For starters, I'd also be happy with a dirty hack, as this feature is the one thing currently keeping me from using NetworkDynamics.jl. My use-case is in neuro-science, but one of you mentioned a use-case for modeling failures in a power-grid, and I imagine there are many more, so I'm sure this would be a useful feature to have!

[antonplietzsch]

Thanks for your request! We have not forgotten about it. 😉

One of our students has re-implemented a simple example for overload line failures in a five node power network as described in Schäfer et al. (2018) using NetworkDynamics. The basic idea of this example is to set down the capacity parameter of an edge to zero if the flow on it exceeds a critical value. This is achieved by using a VectorContinuousCallback. It is possible to interact with the state of the system by calling the network rhs with the signature (u, p, t, GetGD) as shown in the example. This will return a GraphData object, that allows for reading and manipulating the underlying data. It should be also straight forward to affect other nodes/edges than the one that triggered the event.

I'm not entirely sure if this approach solves your problem. We are currently working on a more simple example and a tutorial for event handling with NetworkDynamics and eventually will also build an interface to make it more convenient for the user.

using NetworkDynamics
using LightGraphs
using DifferentialEquations
using Plots

struct SwingVertex
    P::Float64 # power of node
    I::Float64 # inertia of node
    γ::Float64 # damping of node
end
function (sv::SwingVertex)(dv, v, e_s, e_d, p, t)
    # v[1] -> δ, dv[1] = dγ = ω = v[2]
    # v[2] -> ω, dv[2] -> dω
    dv[1] = v[2]
    dv[2] = sv.P - sv.γ*v[2] + flow_sum(e_s, e_d)
    dv[2] = dv[2]/sv.I
    nothing
end

function flow_sum(e_s, e_d)
    sum = 0.0
    for e in e_d
        sum -= e[1]
    end
    for e in e_s
        sum += e[1]
    end
    return sum
end

mutable struct SwingEdge
    from::Int8
    to::Int8
    K::Float64
    max_flow::Float64
end
function (se::SwingEdge)(e, v_s, v_d, p, t)
    e[1] = se.K * sin(v_d[1] - v_s[1])
end

"""
Function returns callable object, which can be used as 'affect!' of a
callback. The affect! will set the coupling of edge idx to zero.
"""
function kill_line_affect(idx)
    function affect!(integrator)
        @info "Line $idx killed at t=$(integrator.t)."
        integrator.f.f.edges![idx].f!.K = 0
        u_modified!(integrator, false)
    end
    return affect!
end


"""
Function returns callable object, which can be used as 'affect!' of a
callback. The affect! will set the coupling of edge (source, destination)
in network to zero.
"""
function kill_line_affect(network, source, destination)
    gs = network.f.graph_structure
    s_ind = findall(x->x==source, gs.s_e)
    d_ind = findall(x->x==destination, gs.d_e)
    ind = intersect(s_ind, d_ind)
    @assert length(ind) == 1
    return kill_line_affect(ind[1])
end


"""
Function returns a VectorContinoursCallback which compares the max_flow of all
edges in the network an shuts them down once the flow is to high.
"""
function watch_line_limit_callback(network)
    num_e = network.f.graph_structure.num_e
    edges = network.f.edges!

    function condition(out, u, t, integrator)
        # get current edge values from integrator
        graph_data = integrator.f.f(u, integrator.p, integrator.t, GetGD)
        edge_values = graph_data.e_array
        for i in 1:num_e
            out[i] = edges[i].f!.max_flow - abs(edge_values[i])
        end
    end

    function affect!(integrator, idx)
        kill_line_affect(idx)(integrator)
    end
    return VectorContinuousCallback(condition, affect!, num_e)
end

function plot_flow(saved_edge_data, network)
    t = saved_edge_data.t
    vals = saved_edge_data.saveval

    data = zeros(length(t), length(vals[1]))
    for i in 1:length(vals)
        # devide by 1.63 to normalize
        data[i, :] = abs.(vals[i])/1.63
    end

    gs = network.f.graph_structure
    sym = Array{String}(undef, (1,gs.num_e))
    for i in 1:gs.num_e
        sym[1,i] = string(gs.s_e[i]) * "->" * string(gs.d_e[i])
    end

    plot(t, data, label=sym, size=(800,600))
end

# Definition of nodes and edges according to schäfer18
I = 1.0
γ = 0.1
verticies = [
    SwingVertex(-1.0, I, γ), #1
    SwingVertex( 1.5, I, γ), #2
    SwingVertex(-1.0, I, γ), #3
    SwingVertex(-1.0, I, γ), #4
    SwingVertex( 1.5, I, γ), #5
]

swingedges = [
    SwingEdge(1, 2, 1.63, 0.6*1.63), #1
    SwingEdge(1, 3, 1.63, 0.6*1.63), #2
    SwingEdge(1, 5, 1.63, 0.6*1.63), #3
    SwingEdge(2, 3, 1.63, 0.6*1.63), #4
    SwingEdge(2, 4, 1.63, 0.6*1.63), #5
    SwingEdge(3, 4, 1.63, 0.6*1.63), #6
    SwingEdge(4, 5, 1.63, 0.6*1.63)  #7
]

# BEWARE: ordering of edges ≠ ordering of add_edge calls!
# In this case I've ordered them manualy for this to match.
g = SimpleGraph(length(verticies))
for e in swingedges
    add_edge!(g, e.from, e.to)
end

# Define nodes/edges and network
odeverts = [ODEVertex(f! = vert, dim=2, sym=[:d, :w]) for vert in verticies]
staticedges = [StaticEdge(f! = edge, dim=1, sym=[:F]) for edge in swingedges]
swing_network! = network_dynamics(odeverts, staticedges, g)

# x0 determined from static solution
# x0 = [0.0 for i in 1:2*nv(g)]
x0 = [-0.12692637482862684, -1.3649456633810975e-6, 0.14641121510104085, 4.2191082676726005e-7, -0.24376507587890778, 1.567589744768255e-6, -0.12692637482862684, -1.3649456633810975e-6, 0.35120661043511864, 7.403907552948938e-7]

# define callback to save edge values (for plotting)
function save_edges(u, t, integrator)
    graph_data = integrator.f.f(u, integrator.p, integrator.t, GetGD)
    return copy(graph_data.e_array)
end
saved_edgevalues = SavedValues(Float64, Array{Float64, 1})
save_callback = SavingCallback(save_edges, saved_edgevalues)

# define problem
prob = ODEProblem(swing_network!, x0, (0., 6))

# define kill_first callback, at t=1.0s remove edge(2,4)
kill_first = PresetTimeCallback( 1.0, kill_line_affect(swing_network!, 2, 4))

# define load_watch_callback
load_watch_callback = watch_line_limit_callback(swing_network!)

@info "Start solver"
sol = solve(prob, Tsit5(), callback=CallbackSet(save_callback, kill_first, load_watch_callback), dtmax=0.01)

# plot(sol, vars = idx_containing(swing_network!, :w))
# plot(sol, vars = idx_containing(swing_network!, :d))
plot_flow(saved_edgevalues, swing_network!)

[jleugeri]

Thanks for the reply and code! I'll try to get something basic implemented this way - if I run into major roadblocks or have some useful input, I'll let you know.
In either case, I'd personally still leave the issue open until there's an actual "user-friendly" interface for this sort of stuff, but I also won't mind if you prefer to close it instead.

[KalelR]

Hey, has there been an update on this issue? I have a similar use case for synaptic coupling between neurons, in which I need to detect an event, save the time of that event, and use that time to send pulses to coupled units. Specifically, when a neuron j spikes (one of its variables increases above a threshold) at time t_s, it sends a current to its neighboring neurons that is a function of the time difference t - t_s. So the coupling term also needs to access the spike times of each neuron.

I would like to implement this with NetworkDynamics (nice package btw!), but I would first like to have an intuition if using NetworkDynamics for this will be efficient, or whether I should write a more specific code for this.

Thanks a lot!

[lindnemi]

Hey @KalelR, thank you for your interest in this package :)

I am not an expert on spiking neural networks but in general it sounds possible. In general this would need to be implemented with the DifferentialEquations.jl callback interface. One could add a further internal array of spiketimes for each neuron for a subtype of NetworkDiffEq and an interface for conveniently building callbacks with access to these arrays.

There was a BoF session at this years JuliaCon 2022 on spiking neural nets. There seems to be a fairly large community of researchers who want to use this, here is a thread on discourse (slightly outdated): https://discourse.julialang.org/t/spiking-neural-networks/63270/17

One of the packages mentioned in this thread which is still actively developed would be: https://github.com/wsphillips/Conductor.jl

[lindnemi]

That said, if you want to build this into ND.jl i would gladly merge PRs and help with integrating it into the core of the package but i don't have the resources (or background on Spiking NNs) to do it myself at the moment

[KalelR]

Hi @lindnemi, thanks for the answer! Indeed, I had a similar idea to your approach and it seems to work well. There's currently a small bug with the Callback interface when neurons are synchronized. I'll first work on fixing it and then at some point it would be nice to work in ND.jl. I quite like the flexibility of the package! Once I have time I can do this but no hurry for now :)