Further steps for Delay differential equations

[lindnemi]

If we want our package to be used by other communities like neurodynamics, we should support DDEs. Any ideas on how to do that are welcome!

[FHell]

DDE Documentation of DiffEq:

https://docs.sciml.ai/stable/types/dde_types/

The DDEFunction of DiffEq:

https://docs.sciml.ai/stable/features/performance_overloads/#DDEFunction-1

A simple start to have some support could look like that:

We will consider systems with only one delay tau. I will consider only delays in the vertex dynamics and static links. We need to provide a function with this signature: f(du, u, h, p, t) where h is a function. We can do this quite simply like this:

function (d::nd_DDE_Static)(dx, x, h, p, t)
    gd = prep_gd(dx, x, d.graph_data, d.graph_structure)
    h_tau = h(p, t - d.tau)

    @nd_threads d.parallel for i in 1:d.graph_structure.num_e
        maybe_idx(d.edges!, i).f!(gd.e[i], gd.v_s_e[i], gd.v_d_e[i], p_e_idx(p, i), t)
    end

    @nd_threads d.parallel for i in 1:d.graph_structure.num_v
        maybe_idx(d.vertices!,i).f!(
            view(dx,d.graph_structure.v_idx[i]),
            gd.v[i], 
            view(h_tau, d.graph_structure.v_idx[i]),
            gd.e_s_v[i], gd.e_d_v[i],
            p_v_idx(p, i),  t)
    end

    nothing
end

So the only thing we need to add is a field tau in the callable struct nd_DDE_Static, and the machinery for a DDEVertex.

Disadvantages of this approach other than being highly restrictive: allocates a new h_tau (could be replaced by an AD aware cache), evaluates the whole history vector even if we don't need it.

[lindnemi]

To enable varying tau, it could be either a global const or passed as a parameter at index p[end].

To avoid allocations h could be in-place.

Further optimization by indexing seem to be more complicated.

[FHell]

To enable varying tau, it could be either a global const or passed as a parameter at index p[end].

If it's part of the network right hand side

To avoid allocations h could be in-place.

Yes, that's what I meant with: "h_tau (could be replaced by an AD aware cache)"

The relevant bit of DiffEq documentation is here:

https://docs.sciml.ai/stable/basics/faq/#I-get-Dual-number-errors-when-I-solve-my-ODE-with-Rosenbrock-or-SDIRK-methods-1

[FHell]

Except that the example code doesn't currently work... might be best to figure out why and fill an issue at OrdinaryDiffEq to get this sorted first.

[lindnemi]

I set up a separate branch to develop DDEs. So far i implemented a DDEVertex with a new call signature that is compatible with DDE solvers. It's working and looking good so far.

Here is a plot of a kuramoto model with delay on a small network.

Example script can be found here: https://github.com/FHell/NetworkDynamics.jl/blob/DDEVertex/examples/getting-started-with-DDEs.jl

With a few more changes this branch can be merged into master:

• [x] Improve the interface for heterogeneous vertices / nodes. At the moment, if there is a DDEVertex all edges have to take additional arguments h_v_s and h_v_d for the history of source and destination vertices and all vertex function have to take h_v. All vertices have to be specified as DDEVertex, even if they don't use the history. We need promotion rules for these cases.

• [ ] DDEEdges (though this is not high priority)

• [x] For performance: Pre-allocate the views into the history array

• [ ] ODEEdges that are aware of vertex delays (ODEDelayEdge)

With what I learned from implementing DDEVertices I would now favour an approach for multilayer structure (#37 ) where every layer is represented by a function that wraps the inner loop according to the type of edge (layer_ODE, layer_Static, layer_DDE). This will break our code into smaller pieces, make it more readable and potentially allow better optimizations. For a list of Vertex Functions of different types (Static, ODE, DDE) that each may have different signatures we could also try a vararges approach.

[lindnemi]

Promotion rules added: needs more testing https://github.com/FHell/NetworkDynamics.jl/commit/5ba6a6c537a28574705f9883ce4189e35d4f9eb2

[lindnemi]

putting the more fancy things aside the branch should be almost ready for merging. we just need to add some tests and documentation

[lindnemi]

Our DDE handling should support hetereogeneous delays. We might implement it in the following way by reforming our current approach:

Define a new type

struct VertexHistoryFunction{h}
    h # global history function
    ghb::GraphDataBuffer
    idxs::Idxs
end

This type can be called and returns a view into the GraphHistoryBuffer ghb (which is essentially another GraphDataBuffer just for history values)

function (vhf::VertexHistoryFunction)(p,t)
    vhf.h(vhf.ghb, p, t, idxs=vhf.idxs)
    view(vhf.ghb, vhf.idxs)
end

This VertexHistoryFunctions can be passed as arguments to StaticDelayEdge to do something like

function sdedge!(e, v_s, v_d, h_v_s, h_v_d, τ, t)
    sin(v_d .- h_v_s(τ, t))
end

And it should be able to do everything allocation free since views dont allocate since Julia 1.5 and the vertexhistoryfunction can be preallocated when constructing the DDEFunction.

[lindnemi]

The history function h is not known when network_dynamics constructs the NetworkDE.

At the moment i see essentially two options to move on.

• Have three arrays v_history, v_s_history, v_d_history. Each ComponentFunction gets views into these arrays. History values are computed before the core loop by looping over h!(buffer[idxs], p,t-tau; idxs=idxs) for approriate idxs and tau. This will allow each edge and each vertex to access history values at exactly one, possibly heterogeneous lag time.

• OR: Pass the history function h on towards the ComponentFunctions. The Components will then additionally need to know their indices in the buffer array. To avoid allocations these indices need to be passed as views, which means components get write access into the internal indexing structure. ComponentFunctions need to call h(p,t-tau, idxs=v_s_idx[1]). No buffer is necessary. Multiple lags per component are possible.

[lindnemi]

The most reasonable way at the moment seems to be to pass h to each component. I set up an attempt in the branch hetero_delay. Even though directly calling h within a single component functions causes no allocations, allocations appear within the full NetworkDE. The reason seems to be, that the type of h can then not be inferred, cf. https://github.com/PIK-ICoNe/NetworkDynamics.jl/blob/115f16362459807d02676871d00f9437eaa6d924/test/delay_test.jl#L104

@hexaeder could you take a look?

[lindnemi]

adding type information everywhere fixed it. thanks for the pointer!