Heterogeneous lags for networks with delay

.github/workflows/tests.yml

@@ -12,7 +12,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.6'
+          - '1.8'
           - '1'
         os:
           - ubuntu-latest

Project.toml

@@ -1,7 +1,7 @@
 name = "NetworkDynamics"
 uuid = "22e9dc34-2a0d-11e9-0de0-8588d035468b"
 authors = ["Frank Hellmann <[email protected]>, Michael Lindner <[email protected]>"]
-version = "0.7.1"
+version = "0.8"
 
 [deps]
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"

examples/Project.toml

@@ -5,6 +5,7 @@ authors = ["Frank Hellmann <[email protected]>", "Michael Lindner <mlindne
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 DelayDiffEq = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GraphPlot = "a2cc645c-3eea-5389-862e-a155d0052231"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -18,4 +19,4 @@ StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-NetworkDynamics = "0.5"
+NetworkDynamics = "≥ 0.5"

examples/kura_delay.jl

@@ -0,0 +1,94 @@
+using DelayDiffEq
+using Random
+using LightGraphs
+using NetworkDynamics
+using Distributions
+using BenchmarkTools
+
+N = 6
+g = SimpleDiGraph(complete_graph(N))
+
+
+function kuramoto!(dθ, θ, edges, p, t)
+    ω = p
+    dθ[1] = ω
+    for e in edges
+        dθ[1] += e[1]
+    end
+    return nothing
+end
+
+function delay_coupling(e, θ_s, θ_d, uh, θ_s_idxs, θ_d_idxs, p, global_p, t)
+    hist1 = uh(t - p[1]; idxs=1)
+    e[1] = p[2] * sin(hist1 - θ_d[1])
+    return nothing
+end
+
+vertex = ODEVertex(f! = kuramoto!, dim = 1)
+edge = StaticDelayEdge(f! = delay_coupling, dim=1, coupling=:directed)
+nd = network_dynamics(vertex,edge,g)
+
+Random.seed!(1)
+
+ω = rand(nv(g))
+τ = rand(ne(g))
+minimum(τ)
+κ = rand(ne(g))
+p = (ω, [τ'; κ'])
+
+θ₀ = rand(Uniform(0, 2π), N)
+
+const past = rand(Uniform(0, 2π), N)
+h(p, t; idxs=nothing) = past[idxs]
+
+x0 = ones(N)
+dx0 = similar(x0)
+@allocated nd.f(dx0, x0,h,p,0) # Works as expected
+@allocated nd.f(dx0, x0,h,p,0) # Works as expected
+@btime $nd.f($dx0, $x0,$h,$p,0) # Works as expected
+
+prob = DDEProblem(nd, θ₀, h, (0.0, 1.0), p, constant_lags=τ)
+
+@btime solve(prob, MethodOfSteps(BS3()), saveat=0.01)
+
+#### Other system
+
+using DelayDiffEq
+using Distributions
+using Random
+
+function f!(dθ, θ, h::H, p, t) where H
+    ω, A = p
+    n = length(θ)
+    lags = reshape(lag, n,n)
+    @inbounds for j in 1:n
+        coupling = 0.0
+        @inbounds for i in 1:n
+            coupling += foo(dθ, θ, h::H, p, t, i, j)
+        end
+        dθ[j] = ω[j] + coupling
+    end
+    nothing
+end
+
+function foo(dθ, θ, h::H, p, t, i, j) where H
+    ω, A = p
+    n = length(θ)
+    lags = reshape(lag, n,n)
+    hist = h(p, t-lags[i,j]; idxs=i)
+    return A[i,j]*sin(hist - θ[j])
+end
+
+n = 10
+Random.seed!(1)
+ω = rand(n)
+A = rand(n,n)
+const lag = rand(n*n)
+θ₀ = rand(Uniform(0, 2π), n)
+p = (ω, A)
+const past = rand(Uniform(0, 2π), n)
+h(p, t; idxs=nothing) = typeof(idxs) <: Number ? past[idxs] : past
+
+prob = DDEProblem(f!, θ₀, h, (0.0, 1.0), p, constant_lags=lag)
+
+@btime solve(prob, MethodOfSteps(BS3()), saveat=0.01, reltol=0, abstol=1e-5)

src/ComponentFunctions.jl

@@ -48,7 +48,7 @@ the source and destination of the described edge.
 
   - `dim` is the number of independent variables in the edge equations and
   - `sym` is an array of symbols for these variables.
-  - `coupling` is a Symbol describing if the EdgeFunction is intended for a directed graph (`:directed`) or for an undirected graph (`{:undirected, :symmetric, :antisymmetric, :fiducial}`). `:directed` is intended for directed graphs. `:undirected` is the default option and is only compatible with SimpleGraph. in this case f! should specify the coupling from a source vertex to a destination vertex. `:symmetric` and `:antisymmetric` trigger performance optimizations, if `f!` has that symmetry property. `:fiducial` lets the user specify both the coupling from src to dst, as well as the coupling from dst to src and is intended for advanced users.
+  - `coupling` is a Symbol describing if the EdgeFunction is intended for a directed graph (`:directed`) or for an undirected graph (`{:undirected, :symmetric, :antisymmetric, :fiducial}`). `:directed` is intended for directed graphs. `:undirected` is the default option and is only compatible with SimpleGraph. in this case f should specify the coupling from a source vertex to a destination vertex. `:symmetric` and `:antisymmetric` trigger performance optimizations, if `f` has that symmetry property. `:fiducial` lets the user specify both the coupling from src to dst, as well as the coupling from dst to src and is intended for advanced users.
 
 For more details see the documentation.
 """
@@ -285,7 +285,7 @@ Here `dv`, `v`, `p` and `t` are the usual ODE arguments, while
 
   - `dim` is the number of independent variables in the edge equations and
   - `sym` is an array of symbols for these variables.
-  - `coupling` is a Symbol describing if the EdgeFunction is intended for a directed graph (`:directed`) or for an undirected graph (`{:undirected, :fiducial}`). `:directed` is intended for directed graphs. `:undirected` is the default option and is only compatible with SimpleGraph. in this case f! should specify the coupling from a source vertex to a destination vertex.  `:fiducial` lets the user specify both the coupling from src to dst, as well as the coupling from dst to src and is intended for advanced users.
+  - `coupling` is a Symbol describing if the EdgeFunction is intended for a directed graph (`:directed`) or for an undirected graph (`{:undirected, :fiducial}`). `:directed` is intended for directed graphs. `:undirected` is the default option and is only compatible with SimpleGraph. in this case f should specify the coupling from a source vertex to a destination vertex.  `:fiducial` lets the user specify both the coupling from src to dst, as well as the coupling from dst to src and is intended for advanced users.
   - `mass_matrix` is an optional argument that defaults to the identity
     matrix `I`. If a mass matrix M is given the system `M * de = f` will be
     solved.
@@ -329,52 +329,26 @@ Base.@kwdef struct StaticDelayEdge{T} <: EdgeFunction
             return new{T}(user_f, dim, coupling, sym)
 
         elseif coupling == :fiducial
-            dim % 2 != 0 && error("Fiducial edges are required to have even dim. ",
-                                  "The first dim args are used for src -> dst ",
-                                  "the second for dst -> src coupling.")
+            dim % 2 == 0 ? nothing : error("Fiducial edges are required to have even dim.
+                                            The first dim args are used for src -> dst,
+                                            the second for dst -> src coupling.")
             return new{T}(user_f, dim, coupling, sym)
-
+        elseif coupling ∈ (:antisymmetric, :symmetric)
+            error("$coupling coupling not implemented for edges with delay. If you need it please open an issue on GitHub.")
         elseif coupling == :undirected
-            # This might cause unexpected behaviour if source and destination vertex don't
-            # have the same internal arguments.
-            # Make sure to explicitly define the edge is :fiducial in that case.
+        # This might cause unexpected behaviour if source and destination vertex don't
+        # have the same internal arguments.
+        # Make sure to explicitly define the edge is :fiducial in that case.
             f = @inline (e, v_s, v_d, h_v_s, h_v_d, p, t) -> begin
-                @inbounds user_f(view(e, 1:dim), v_s, v_d, h_v_s, h_v_d, p, t)
-                @inbounds user_f(view(e, dim+1:2dim), v_d, v_s, h_v_d, h_v_s, p, t)
-                nothing
-            end
-        elseif coupling == :antisymmetric
-            f = @inline (e, v_s, v_d, h_v_s, h_v_d, p, t) -> begin
-                @inbounds user_f(view(e, 1:dim), v_s, v_d, h_v_s, h_v_d, p, t)
-                @inbounds for i in 1:dim
-                    e[dim+i] = -1.0 * e[i]
-                end
-                nothing
-            end
-        elseif coupling == :symmetric
-            f = @inline (e, v_s, v_d, h_v_s, h_v_d, p, t) -> begin
-                @inbounds user_f(view(e, 1:dim), v_s, v_d, h_v_s, h_v_d, p, t)
-                @inbounds for i in 1:dim
-                    e[dim+i] = e[i]
-                end
+                @inbounds user_f(view(e,1:dim), v_s, v_d, h_v_s, h_v_d, p, t)
+                @inbounds user_f(view(e,dim+1:2dim), v_d, v_s, h_v_d, h_v_s, p, t)
                 nothing
             end
         end
         return new{typeof(f)}(f, 2dim, coupling, repeat(sym, 2))
     end
 end
 
-function StaticDelayEdge(se::StaticEdge)
-    let _f = se.f, dim = se.dim, coupling = se.coupling, sym = se.sym
-        f = @inline (e, v_s, v_d, h_v_s, h_v_d, p, t) -> _f(e, v_s, v_d, p, t)
-        if coupling ∈ (:undefined, :fiducial, :directed)
-            return StaticDelayEdge(f, dim, coupling, sym)
-        else
-            return StaticDelayEdge(f, dim, :fiducial, sym)
-        end
-    end
-end
-
 # Promotion rules
 
 

src/NetworkDiffEq.jl

@@ -30,11 +30,11 @@ end
 
 
 @inline function component_loop!(unique_components, unique_c_indices,
-                                 dx, p, t, gd, gs, history, parallel)
+    dx, p, t, gd, gs, h::H, parallel) where H
     for j in 1:length(unique_components)
         # Function barrier
         _inner_loop!(unique_components[j], unique_c_indices[j],
-                     dx, p, t, gd, gs, history, parallel)
+                     dx, p, t, gd, gs, h::H, parallel)
     end
     return nothing
 end
@@ -43,7 +43,7 @@ end
 # inner loops for ODE + Static case
 
 function _inner_loop!(component::ODEVertex, indices,
-                      dx, p, t, gd, gs, history, parallel)
+                      dx, p, t, gd, gs, h, parallel)
     @nd_threads parallel for i in indices
         component.f(view(dx, gs.v_idx[i]),
                     get_vertex(gd, i),
@@ -55,7 +55,7 @@ function _inner_loop!(component::ODEVertex, indices,
 end
 
 function _inner_loop!(component::StaticEdge, indices,
-                      dx, p, t, gd, gs, history, parallel)
+                      dx, p, t, gd, gs, h, parallel)
     @nd_threads parallel for i in indices
         component.f(get_edge(gd, i),
                     get_src_vertex(gd, i),
@@ -69,34 +69,43 @@ end
 # inner loops for DDE + Static Delay
 
 function _inner_loop!(component::DDEVertex, indices,
-                      dx, p, t, gd, gs, history, parallel)
+                      dx, p, t, gd, gs, h::H, parallel) where H
+
     @nd_threads parallel for i in indices
+        # Wrappers for the history function correct for global p and global idx
+        h_v = @inline((t; idxs) -> h(p,t;idxs=gs.v_idx[i][idxs]))
+
         component.f(view(dx, gs.v_idx[i]),
-                    get_vertex(gd, i),
-                    get_dst_edges(gd, i),
-                    view(history, gs.v_idx[i]),
-                    p_v_idx(p, i),
-                    t)
+                  get_vertex(gd, i),
+                  get_dst_edges(gd, i),
+                  h::H,
+                  gs.v_idx[i],
+                  p_v_idx(p, i),
+                  t)
     end
     return nothing
 end
 
 function _inner_loop!(component::StaticDelayEdge, indices,
-                      dx, p, t, gd, gs, history, parallel)
+                      dx, p, t, gd, gs, h::H, parallel) where H
     @nd_threads parallel for i in indices
+        # Wrappers for the history function correct for global p and global idx
+        h_v_s = @inline((t; idxs) -> h(p,t;idxs=gs.s_e_idx[i][idxs]))
+        h_v_d = @inline((t; idxs) -> h(p,t;idxs=gs.d_e_idx[i][idxs]))
+
         component.f(get_edge(gd, i),
-                    get_src_vertex(gd, i),
-                    get_dst_vertex(gd, i),
-                    view(history, gs.s_e_idx[i]),
-                    view(history, gs.d_e_idx[i]),
-                    p_e_idx(p, i),
-                    t)
+                     get_src_vertex(gd, i),
+                     get_dst_vertex(gd, i),
+                     h_v_s,
+                     h_v_d,
+                     p_e_idx(p, i),
+                     t)
     end
     return nothing
 end
 
 function _inner_loop!(component::ODEEdge, indices,
-                      dx, p, t, gd, gs, history, parallel)
+                      dx, p, t, gd, gs, h, parallel)
     @nd_threads parallel for i in indices
         component.f(view(dx, gs.e_idx[i] .+ gs.dim_v),
                     get_edge(gd, i),
@@ -111,15 +120,14 @@ end
 # struct for both cases
 
 #@Base.kwdef struct NetworkDE{G, GDB, elV, elE, TUV, TUE, Th<:AbstractArray{elV}}
-Base.@kwdef struct NetworkDE{G,GDB,elV,elE,TUV,TUE,Th<:Union{AbstractArray,Nothing}}
+@Base.kwdef struct NetworkDE{G, GDB, elV, elE, TUV, TUE}
     unique_vertices!::TUV
     unique_v_indices::Vector{Vector{Int}}
     unique_edges!::TUE
     unique_e_indices::Vector{Vector{Int}}
     graph::G #redundant?
     graph_structure::GraphStruct
-    graph_data::GraphData{GDB,elV,elE}
-    history::Th # for ODE + Static case: nothing, for DDE + Static Delay case: Th
+    graph_data::GraphData{GDB, elV, elE}
     parallel::Bool # enables multithreading for the core loop
 end
 
@@ -132,20 +140,27 @@ function (d::NetworkDE)(dx, x, p, t)
 
     @assert size(dx) == size(x) "Sizes of dx and x do not match"
 
-    # Pass nothing, because here we have only Static Edges or Static Delay Edges (if we
-    # include the ODE ODE case than we have to pass dx)
+    # Pass nothing instead of the history function
     component_loop!(d.unique_edges!, d.unique_e_indices,
-                    dx, p, t, gd, gs, d.history, d.parallel)
+                    dx, p, t, gd, gs, nothing, d.parallel)
 
     component_loop!(d.unique_vertices!, d.unique_v_indices,
-                    dx, p, t, gd, gs, d.history, d.parallel)
+                    dx, p, t, gd, gs, nothing, d.parallel)
     return nothing
 end
 # for DDE case
-function (d::NetworkDE)(dx, x, h!, p, t)
-    # History computation happens beforehand and is cached in d.history
-    h!(d.history, p, t - p[end])
-    d(dx, x, p, t)
+function (d::NetworkDE)(dx, x, h::H, p, t) where H
+    gs = d.graph_structure
+    checkbounds_p(p, gs.num_v, gs.num_e)
+    gd = prep_gd(dx, x, d.graph_data, gs)
+
+    @assert size(dx) == size(x) "Sizes of dx and x do not match"
+
+    component_loop!(d.unique_edges!, d.unique_e_indices,
+                    dx, p, t, gd, gs, h::H, d.parallel)
+
+    component_loop!(d.unique_vertices!, d.unique_v_indices,
+                    dx, p, t, gd, gs, h::H, d.parallel)
     return nothing
 end
 
@@ -157,7 +172,7 @@ function (d::NetworkDE)(x, p, t, ::Type{GetGD})
     if size(x) == (gs.dim_v,)
         checkbounds_p(p, gs.num_v, gs.num_e)
         component_loop!(d.unique_edges!, d.unique_e_indices,
-                        nothing, p, t, gd, gs, d.history, d.parallel)
+                        nothing, p, t, gd, gs, nothing, d.parallel)
     end
     gd
 end

src/NetworkDynamics.jl

@@ -175,16 +175,11 @@ function _network_dynamics(vertices!::Union{Vector{T},T},
     hasDelay = any(v -> v isa DDEVertex, unique_vertices!) ||
                any(e -> e isa StaticDelayEdge, unique_edges!)
 
-    if initial_history === nothing && hasDelay
-        initial_history = ones(sum(v_dims))
-    end
-
     graph_stucture = GraphStruct(graph, v_dims, e_dims, symbols_v, symbols_e)
 
     graph_data = GraphData(v_array, e_array, graph_stucture)
 
-    nd! = NetworkDE(unique_vertices!, unique_v_indices, unique_edges!, unique_e_indices, graph, graph_stucture,
-                    graph_data, initial_history, parallel)
+    nd! = NetworkDE(unique_vertices!, unique_v_indices, unique_edges!, unique_e_indices, graph, graph_stucture, graph_data, parallel)
     mass_matrix = construct_mass_matrix(mmv_array, graph_stucture)
 
     return if hasDelay
@@ -196,8 +191,7 @@ end
 
 ## ODEEdge
 
-function _network_dynamics(vertices!::Union{Vector{T},T}, edges!::Union{Vector{U},U}, graph; initial_history=nothing,
-                           x_prototype=zeros(1), parallel=false) where {T<:ODEVertex,U<:ODEEdge}
+function _network_dynamics(vertices!::Union{Vector{T},T}, edges!::Union{Vector{U},U}, graph; x_prototype=zeros(1), parallel=false) where {T<:ODEVertex, U<:ODEEdge}
 
     warn_parallel(parallel)
 
@@ -223,8 +217,7 @@ function _network_dynamics(vertices!::Union{Vector{T},T}, edges!::Union{Vector{U
 
     graph_data = GraphData(v_array, e_array, graph_stucture)
 
-    nd! = NetworkDE(unique_vertices!, unique_v_indices, unique_edges!, unique_e_indices, graph, graph_stucture,
-                    graph_data, initial_history, parallel)
+    nd! = NetworkDE(unique_vertices!, unique_v_indices, unique_edges!, unique_e_indices, graph, graph_stucture, graph_data, parallel)
 
     mass_matrix = construct_mass_matrix(mmv_array, mme_array, graph_stucture)