From 691e542fda3623822600b98f1c208f62f276efea Mon Sep 17 00:00:00 2001 From: lindnemi Date: Mon, 18 Oct 2021 10:46:47 +0200 Subject: [PATCH 01/30] refactor submodules --- Project.toml | 4 +--- src/ComponentFunctions.jl | 12 ++++------ src/{NetworkDE_mod.jl => NetworkDiffEq.jl} | 8 ------- src/NetworkDynamics.jl | 9 +------- src/NetworkStructures.jl | 22 +++++++++---------- src/Utilities.jl | 6 +---- ...etworkDE_test.jl => NetworkDiffEq_test.jl} | 12 +++++----- test/runtests.jl | 2 +- 8 files changed, 25 insertions(+), 50 deletions(-) rename src/{NetworkDE_mod.jl => NetworkDiffEq.jl} (98%) rename test/{NetworkDE_test.jl => NetworkDiffEq_test.jl} (93%) diff --git a/Project.toml b/Project.toml index 526321761..d4c884c8f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,19 +1,17 @@ name = "NetworkDynamics" uuid = "22e9dc34-2a0d-11e9-0de0-8588d035468b" authors = ["Frank Hellmann , Michael Lindner "] -version = "0.5.7" +version = "0.6.0" [deps] DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" -Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] DiffEqBase = "^6.9.4" LightGraphs = "^1.3.0" NLsolve = "^4.2.0" -Reexport = "^1.0" julia = "1.5" diff --git a/src/ComponentFunctions.jl b/src/ComponentFunctions.jl index 661982897..9c3d7f8a7 100644 --- a/src/ComponentFunctions.jl +++ b/src/ComponentFunctions.jl @@ -1,9 +1,7 @@ -""" -Together with Constructors this module forms the backbone of the core API. -It provide the basic types to construct Arrays of VertexFunction and -EdgeFunction which can be handled by network_dynamics. -""" -module ComponentFunctions +# Together with Constructors this file forms the backbone of the core API. +# It provide the basic types to construct Arrays of VertexFunction and +# EdgeFunction which can be handled by network_dynamics. + using LinearAlgebra @@ -427,5 +425,3 @@ function ODEEdge(se::StaticEdge) end end end - -end diff --git a/src/NetworkDE_mod.jl b/src/NetworkDiffEq.jl similarity index 98% rename from src/NetworkDE_mod.jl rename to src/NetworkDiffEq.jl index f69236f24..3e9aac65d 100644 --- a/src/NetworkDE_mod.jl +++ b/src/NetworkDiffEq.jl @@ -1,9 +1,3 @@ -module NetworkDE_mod - -using ..NetworkStructures -using ..ComponentFunctions -using ..Utilities - export NetworkDE @@ -174,5 +168,3 @@ end function (d::NetworkDE)(::Type{GetGS}) d.graph_structure end - -end # module diff --git a/src/NetworkDynamics.jl b/src/NetworkDynamics.jl index 6aebcb32e..c1496519e 100644 --- a/src/NetworkDynamics.jl +++ b/src/NetworkDynamics.jl @@ -1,20 +1,13 @@ module NetworkDynamics -using Reexport using DiffEqBase using LightGraphs include("Utilities.jl") -@reexport using .Utilities - include("ComponentFunctions.jl") -@reexport using .ComponentFunctions - include("NetworkStructures.jl") -@reexport using .NetworkStructures +include("NetworkDiffEq.jl") -include("NetworkDE_mod.jl") -@reexport using .NetworkDE_mod export network_dynamics diff --git a/src/NetworkStructures.jl b/src/NetworkStructures.jl index 2edaabe33..fd5a21bf0 100644 --- a/src/NetworkStructures.jl +++ b/src/NetworkStructures.jl @@ -1,19 +1,21 @@ -""" - This module contains the logic that calculate the index structures - and data access structs that Network Dynamics makes use of. +# +# This file contains the logic that calculate the index structures +# and data access structs that Network Dynamics makes use of. +# +# The key structure is the GraphData structure that allows accessing data on +# vertices and edges of the graph in an efficient manner. The neccessary indices +# are precomputed in GraphStructure. - The key structure is the GraphData structure that allows accessing data on - vertices and edges of the graph in an efficient manner. The neccessary indices - are precomputed in GraphStructure. -""" -module NetworkStructures using LightGraphs using LinearAlgebra -export GraphStruct, GraphData, EdgeData, VertexData const Idx = UnitRange{Int} +export GraphStruct, GraphData, EdgeData, VertexData + + + """ Create offsets for stacked array of dimensions dims """ @@ -401,5 +403,3 @@ end function (nds::ND_Solution)(t) nds.nd(nds.sol(t), nds.p, t, GetGD) end - -end # module diff --git a/src/Utilities.jl b/src/Utilities.jl index f4c7fd804..bfc5ee9c3 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -1,6 +1,4 @@ -module Utilities - -using NLsolve +import NLsolve: nlsolve, converged using LinearAlgebra using SparseArrays @@ -311,5 +309,3 @@ function construct_mass_matrix(mmv_array, mme_array, gs) end mass_matrix end - -end #module diff --git a/test/NetworkDE_test.jl b/test/NetworkDiffEq_test.jl similarity index 93% rename from test/NetworkDE_test.jl rename to test/NetworkDiffEq_test.jl index f5c507a65..e696b1219 100644 --- a/test/NetworkDE_test.jl +++ b/test/NetworkDiffEq_test.jl @@ -41,14 +41,14 @@ using NetworkDynamics gs = nd(GetGS) # first test array of same type a = rand(length(x)) - gd_ret = NetworkDE_mod.prep_gd(dx, a, gd, gs) + gd_ret = NetworkDynamics.prep_gd(dx, a, gd, gs) @test gd_ret === gd @test gd.gdb.v_array == a # check with type missmatch a = rand(Int, length(x)) dx = similar(a) - gd_ret = NetworkDE_mod.prep_gd(dx, a, gd, gs) + gd_ret = NetworkDynamics.prep_gd(dx, a, gd, gs) @test gd_ret !== gd @test gd_ret.gdb.v_array == a end @@ -99,14 +99,14 @@ end gs = nd(GetGS) # first test array of same type a = rand(length(x)) - gd_ret = NetworkDE_mod.prep_gd(dx, a, gd, gs) + gd_ret = NetworkDynamics.prep_gd(dx, a, gd, gs) @test gd_ret === gd @test gd.gdb.v_array == a # check with type missmatch a = rand(Int, length(x)) dx = similar(a) - gd_ret = NetworkDE_mod.prep_gd(dx, a, gd, gs) + gd_ret = NetworkDynamics.prep_gd(dx, a, gd, gs) @test gd_ret !== gd @test gd_ret.gdb.v_array == a end @@ -151,7 +151,7 @@ end gs = nd(GetGS) # first test array of same type a = rand(length(x)) - gd_ret = NetworkDE_mod.prep_gd(dx, a, gd, gs) + gd_ret = NetworkDynamics.prep_gd(dx, a, gd, gs) @test gd_ret === gd @test gd.gdb.v_array == a[1:gs.dim_v] @test gd.gdb.e_array == a[gs.dim_v+1:end] @@ -159,7 +159,7 @@ end # check with type missmatch a = rand(Int, length(x)) dx = similar(a) - gd_ret = NetworkDE_mod.prep_gd(dx, a, gd, gs) + gd_ret = NetworkDynamics.prep_gd(dx, a, gd, gs) @test gd_ret !== gd @test gd_ret.gdb.v_array == a[1:gs.dim_v] @test gd_ret.gdb.e_array == a[gs.dim_v+1:end] diff --git a/test/runtests.jl b/test/runtests.jl index 222b1babd..1394ac0b0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,7 +3,7 @@ using Test # basic tests for individual modules include("NetworkStructures_test.jl") include("ComponentFunctions_test.jl") - include("NetworkDE_test.jl") + include("NetworkDiffEq_test.jl") include("checkbounds_test.jl") # complex tests of networks From 9d85d937c6c5fb5948f4474347e204e3183cc095 Mon Sep 17 00:00:00 2001 From: lindnemi Date: Mon, 18 Oct 2021 10:46:55 +0200 Subject: [PATCH 02/30] remove unused code --- src/NetworkDynamics.jl | 10 ---- test/time_and_test_refactor.jl | 99 ---------------------------------- 2 files changed, 109 deletions(-) delete mode 100644 test/time_and_test_refactor.jl diff --git a/src/NetworkDynamics.jl b/src/NetworkDynamics.jl index c1496519e..981ef5522 100644 --- a/src/NetworkDynamics.jl +++ b/src/NetworkDynamics.jl @@ -322,14 +322,4 @@ end error("Reconstruction of ODEEdges is not implemented at the moment.") end -""" -Allow initializing StaticEdgeFunction for Power Dynamics -""" -function StaticEdgeFunction(vertices!, edges!, graph; parallel = false) - # For reasons I don't fully understand we have to qualify the call to - # the constructor of StaticEdgeFunction here. - NetworkDE_mod.StaticEdgeFunction(network_dynamics(vertices!, edges!, graph, parallel = parallel)) -end - - end # module diff --git a/test/time_and_test_refactor.jl b/test/time_and_test_refactor.jl deleted file mode 100644 index 5dd441628..000000000 --- a/test/time_and_test_refactor.jl +++ /dev/null @@ -1,99 +0,0 @@ -using Pkg -Pkg.activate(".") -using BenchmarkTools -using NetworkDynamics -using Test -using LightGraphs -using Random - - -for N = [20,200] - # Pure Julia "ground truth" - - g = watts_strogatz(N, 5, 0.) - println("\nTesting ND with $N nodes and $(ne(g)) edges:\n") - - B = incidence_matrix(g, oriented=true) - - struct kuramoto_dyn{T, T2, U} - B::T - B_trans::T2 - ω::U - N::Int - end - function (dd::kuramoto_dyn)(dx, x, p, t) - dx[1:N] .= x[N+1:2N] - dx[N+1:2N] .= dd.ω .- x[N+1:2N] .- 5. * dd.B * sin.(dd.B_trans * x[1:N]) - nothing - end - - ω = Array{Float64}(1:N) ./ N - kn! = kuramoto_dyn(B, transpose(B), ω, N) - - # NetworkDynamics - - function kuramoto_edge!(e, θ_s, θ_d, K, t) - e[1] = K * sin(θ_s[1] - θ_d[1]) - nothing - end - - function kuramoto_inertia!(dv, v, dst_edges, ω, t) - dv[1] = v[2] - dv[2] = ω - v[2] - @inbounds for e in dst_edges - dv[2] += e[1] - end - nothing - end - - inertia! = ODEVertex(f! = kuramoto_inertia!, dim = 2, sym= [:θ, :ω]); - edge! = StaticEdge(f! = kuramoto_edge!, dim = 1) - edge_dir! = StaticEdge(f! = kuramoto_edge!, dim = 1, coupling = :directed) - edge_anti! = StaticEdge(f! = kuramoto_edge!, dim = 1, coupling = :antisymmetric ) - - - - nd! = network_dynamics(inertia!, edge!, g) - nd_anti! = network_dynamics(inertia!, edge_anti!, g) - nd_dir! = network_dynamics(inertia!, edge!, SimpleDiGraph(g)) - - # Testing accuracy - - @testset "Accuracy" begin - for i in 1:20 - x0_θ = randn(N) - x0_ω = randn(N) - - x0_kn = [x0_θ; x0_ω ] - - x0_nd = Vector(vec([x0_θ x0_ω]')) - x0_nd_anti = Vector(vec([x0_θ x0_ω]')) - x0_nd_dir = Vector(vec([x0_θ x0_ω]')) - - dx_kn = zeros(2N) - dx_nd = zeros(2N) - dx_nd_anti = zeros(2N) - dx_nd_dir = zeros(2N) - - kn!(dx_kn, x0_kn, nothing, 0.) - nd!(dx_nd, x0_nd, (ω, 5.), 0.) - nd_anti!(dx_nd_anti, x0_nd_anti, (ω, 5.), 0.) - nd_dir!(dx_nd_dir, x0_nd_dir, (ω, 5.), 0.) - - @test isapprox(dx_kn, [dx_nd[1:2:2N]; dx_nd[2:2:2N]]) - @test isapprox(dx_kn, [dx_nd_anti[1:2:2N]; dx_nd_anti[2:2:2N]]) - @test isapprox(dx_kn, [dx_nd_dir[1:2:2N]; dx_nd_dir[2:2:2N]]) - - end - end - println("") - - dx0 = zeros(2N) - x0 = randn(2N) - Random.seed!(42) - display(@benchmark($nd!($dx0, $x0, ($ω, 5.), 0.))) - display(@benchmark($nd_dir!($dx0, $x0, ($ω, 5.), 0.))) - - display(@benchmark($nd_anti!($dx0, $x0, ($ω, 5.), 0.))) - -end From b8a0b0219a4e729aca9d0eb915333fa60bd05fc1 Mon Sep 17 00:00:00 2001 From: lindnemi Date: Mon, 18 Oct 2021 14:12:35 +0200 Subject: [PATCH 03/30] set up a simple test --- test/delay_test.jl | 66 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) create mode 100644 test/delay_test.jl diff --git a/test/delay_test.jl b/test/delay_test.jl new file mode 100644 index 000000000..3cc07d062 --- /dev/null +++ b/test/delay_test.jl @@ -0,0 +1,66 @@ +using Test +using NetworkDynamics +using LightGraphs +using OrdinaryDiffEq +using DelayDiffEq + +N = 10 +g = watts_strogatz(N,2,0.) + +@inline function diffusionedge!(e, v_s, v_d, p, t) + e[1] = v_s[1] - v_d[1] + nothing +end + +@inline function delayedge!(e, v_s, v_d, h_v_s, h_v_d, p, t) + e[1] = h_v_s[1] - v_d[1] + nothing +end + +@inline function diffusionvertex!(dv, v, edges, p, t) + dv[1] = 0. + for e in edges + dv[1] += e[1] + end + nothing +end + +nd_diffusion_vertex = ODEVertex(f! = diffusionvertex!, dim = 1) +nd_diffusion_edge = StaticEdge(f! = diffusionedge!, dim = 1, coupling = :antisymmetric) +nd_delay_edge = StaticDelayEdge(f! = delayedge!, dim = 1) + + +nd! = network_dynamics(nd_diffusion_vertex, nd_diffusion_edge, g) + +x0 = [ones(N÷2); -ones(N÷2)] +tspan = (0.,200.) +prob = ODEProblem(nd!, x0, tspan, nothing) +sol = solve(prob, Tsit5(), abstol=1e-6) + + +@test isapprox(sol[end], zeros(N), atol=1e-6) + + + +dnd! = network_dynamics(nd_diffusion_vertex, nd_delay_edge, g) + +x0 = [ones(N÷2); -ones(N÷2)] +tspan = (0.,200.) +# constant initial conditions on interval +h!(out, p, t) = out .= x0 +# (vertexp, edgep, τ) +p = (nothing, nothing, 1.) + +prob = DDEProblem(dnd!, x0, h!, tspan, p)#, constant_lags=[p[end]]) +sol = solve(prob, MethodOfSteps(Tsit5()), abstol=1e-6) + +isapprox(sol[end], zeros(N), atol=1e-6) + +dx0 = similar(x0) + +nd!.f + +@code_warntype dnd!.f(dx0, x0,nothing,0.) + +@allocated nd!(dx0, x0,nothing,0.) +@allocated dnd!(dx0, x0, p, 0.) From c3c82dd449cc7a5ae04418fc48816d913ee14888 Mon Sep 17 00:00:00 2001 From: lindnemi Date: Mon, 18 Oct 2021 14:21:20 +0200 Subject: [PATCH 04/30] remove internal exports --- src/NetworkDiffEq.jl | 3 --- src/NetworkStructures.jl | 4 ---- src/Utilities.jl | 3 --- test/NetworkStructures_test.jl | 1 + 4 files changed, 1 insertion(+), 10 deletions(-) diff --git a/src/NetworkDiffEq.jl b/src/NetworkDiffEq.jl index 3e9aac65d..546a51eaa 100644 --- a/src/NetworkDiffEq.jl +++ b/src/NetworkDiffEq.jl @@ -1,6 +1,3 @@ -export NetworkDE - - # In order to match the type, we need to pass both, a view that matches the type # to be constructed, and the original array we want to construct a GD on top of. @inline function prep_gd(dx::AbstractArray{T}, x::AbstractArray{T}, gd::GraphData{GDB, T, T}, gs) where {GDB, T} diff --git a/src/NetworkStructures.jl b/src/NetworkStructures.jl index fd5a21bf0..80aa3133c 100644 --- a/src/NetworkStructures.jl +++ b/src/NetworkStructures.jl @@ -12,10 +12,6 @@ using LinearAlgebra const Idx = UnitRange{Int} -export GraphStruct, GraphData, EdgeData, VertexData - - - """ Create offsets for stacked array of dimensions dims """ diff --git a/src/Utilities.jl b/src/Utilities.jl index bfc5ee9c3..5f426b491 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -2,9 +2,6 @@ import NLsolve: nlsolve, converged using LinearAlgebra using SparseArrays -export maybe_idx, p_v_idx, p_e_idx, @nd_threads, construct_mass_matrix, warn_parallel, checkbounds_p - - function warn_parallel(b::Bool) if b haskey(ENV, "JULIA_NUM_THREADS") && diff --git a/test/NetworkStructures_test.jl b/test/NetworkStructures_test.jl index d55c32751..fbdfdd0b2 100644 --- a/test/NetworkStructures_test.jl +++ b/test/NetworkStructures_test.jl @@ -1,6 +1,7 @@ using Test using LightGraphs using NetworkDynamics +import NetworkDynamics: GraphStruct, GraphData, get_vertex, get_edge, get_src_vertex, get_src_edges, get_dst_vertex, get_dst_edges, swap_v_array!, swap_e_array! @testset "Test GraphData Accessors" begin From 115f16362459807d02676871d00f9437eaa6d924 Mon Sep 17 00:00:00 2001 From: lindnemi Date: Thu, 21 Oct 2021 11:16:24 +0200 Subject: [PATCH 05/30] first attempts --- src/ComponentFunctions.jl | 6 +- src/NetworkDiffEq.jl | 23 ++++--- test/delay_test.jl | 129 ++++++++++++++++++++++++++++++-------- 3 files changed, 121 insertions(+), 37 deletions(-) diff --git a/src/ComponentFunctions.jl b/src/ComponentFunctions.jl index 9c3d7f8a7..f22461075 100644 --- a/src/ComponentFunctions.jl +++ b/src/ComponentFunctions.jl @@ -352,9 +352,9 @@ Like a static edge but with extra arguments for the history of the source and de # 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) + f! = @inline (e, v_s, v_d, h, v_s_idx, v_d_idx, p, t) -> begin + @inbounds user_f!(view(e,1:dim), v_s, v_d, h, v_s_idx, v_d_idx, p, t) + @inbounds user_f!(view(e,dim+1:2dim), v_d, v_s, h, v_d_idx, v_s_idx, p, t) nothing end elseif coupling == :antisymmetric diff --git a/src/NetworkDiffEq.jl b/src/NetworkDiffEq.jl index 546a51eaa..8f21628f7 100644 --- a/src/NetworkDiffEq.jl +++ b/src/NetworkDiffEq.jl @@ -82,13 +82,14 @@ function _inner_loop!(component::DDEVertex, indices, end function _inner_loop!(component::StaticDelayEdge, 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), get_dst_vertex(gd, i), - view(history, gs.s_e_idx[i]), - view(history, gs.d_e_idx[i]), + h, # history function + gs.s_e_idx[i], + gs.d_e_idx[i], p_e_idx(p, i), t) end @@ -142,10 +143,18 @@ function (d::NetworkDE)(dx, x, p, t) 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, p, t) + 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, d.parallel) + + component_loop!(d.unique_vertices!, d.unique_v_indices, + dx, p, t, gd, gs, h, d.parallel) return nothing end diff --git a/test/delay_test.jl b/test/delay_test.jl index 3cc07d062..afd325aa8 100644 --- a/test/delay_test.jl +++ b/test/delay_test.jl @@ -3,19 +3,17 @@ using NetworkDynamics using LightGraphs using OrdinaryDiffEq using DelayDiffEq +using BenchmarkTools N = 10 -g = watts_strogatz(N,2,0.) +g = watts_strogatz(N,4,0.) # more edges than vertices! @inline function diffusionedge!(e, v_s, v_d, p, t) e[1] = v_s[1] - v_d[1] nothing end -@inline function delayedge!(e, v_s, v_d, h_v_s, h_v_d, p, t) - e[1] = h_v_s[1] - v_d[1] - nothing -end + @inline function diffusionvertex!(dv, v, edges, p, t) dv[1] = 0. @@ -27,40 +25,117 @@ end nd_diffusion_vertex = ODEVertex(f! = diffusionvertex!, dim = 1) nd_diffusion_edge = StaticEdge(f! = diffusionedge!, dim = 1, coupling = :antisymmetric) -nd_delay_edge = StaticDelayEdge(f! = delayedge!, dim = 1) +nd! = network_dynamics(nd_diffusion_vertex, nd_diffusion_edge, g) -nd! = network_dynamics(nd_diffusion_vertex, nd_diffusion_edge, g) +@testset "ODE diffusion" begin + nd! = network_dynamics(nd_diffusion_vertex, nd_diffusion_edge, g) -x0 = [ones(N÷2); -ones(N÷2)] -tspan = (0.,200.) -prob = ODEProblem(nd!, x0, tspan, nothing) -sol = solve(prob, Tsit5(), abstol=1e-6) + x0 = [ones(N÷2); -ones(N÷2)] + tspan = (0.,5000.) + dx0 = similar(x0) + @btime $nd!.f($dx0,$x0,nothing,0.) + prob = ODEProblem(nd!, x0, tspan, nothing) + sol = solve(prob, Tsit5(), abstol=1e-6) + @btime solve($prob, Tsit5(), abstol=1e-6) + @test isapprox(sol[end], zeros(N), atol=1e-4) +end -@test isapprox(sol[end], zeros(N), atol=1e-6) +@testset "Homogeneous DDE diffusion" begin + @inline function delayedge!(e, v_s, v_d, h_v_s, h_v_d, p, t) + e[1] = h_v_s[1] - v_d[1] + nothing + end + nd_delay_edge = StaticDelayEdge(f! = delayedge!, dim = 1) + dnd! = network_dynamics(nd_diffusion_vertex, nd_delay_edge, g) + x0 = [ones(N÷2); -ones(N÷2)] + tspan = (0.,5000.) + # constant initial conditions on interval + h!(out, p, t) = out .= x0 + # (vertexp, edgep, τ) + p = (nothing, nothing, 1.) -dnd! = network_dynamics(nd_diffusion_vertex, nd_delay_edge, g) + dx0 = similar(x0) + @btime $dnd!.f($dx0,$x0,$h!,$p,0.) -x0 = [ones(N÷2); -ones(N÷2)] -tspan = (0.,200.) -# constant initial conditions on interval -h!(out, p, t) = out .= x0 -# (vertexp, edgep, τ) -p = (nothing, nothing, 1.) + prob = DDEProblem(dnd!, x0, h!, tspan, p)#, constant_lags=[p[end]]) + sol = solve(prob, MethodOfSteps(Tsit5()), abstol=1e-6) + @btime solve($prob, MethodOfSteps(Tsit5()), abstol=1e-6) + @test isapprox(sol[end], zeros(N), atol=1e-4) +end -prob = DDEProblem(dnd!, x0, h!, tspan, p)#, constant_lags=[p[end]]) -sol = solve(prob, MethodOfSteps(Tsit5()), abstol=1e-6) -isapprox(sol[end], zeros(N), atol=1e-6) -dx0 = similar(x0) +@testset "New homogeneous DDE diffusion" begin + @inline function delayedge!(e, v_s, v_d, h, v_s_idx, v_d_idx, p, t) + hist1 = h(nothing, t-1., idxs=v_s_idx[1])::Float64 + e[1] = hist1 - v_d[1] + hist2 = h(nothing, t-1., idxs=v_d_idx[1])::Float64 + e[2] = hist2 - v_s[1] + nothing + end + + nd_delay_edge = StaticDelayEdge(f! = delayedge!, dim = 2, coupling = :fiducial) + + dnd! = network_dynamics(nd_diffusion_vertex, nd_delay_edge, g) -nd!.f + x0 = [ones(N÷2); -ones(N÷2)] + tspan = (0.,5000.) + # constant initial conditions on interval + h(p, t; idxs=nothing) = x0[idxs] + # (vertexp, edgep, τ) + p = (nothing, nothing, 1.) + + dx0 = similar(x0) + @btime $dnd!.f($dx0,$x0,$h,$p,0.) + + prob = DDEProblem(dnd!, x0, h, tspan, p)#, constant_lags=[p[end]]) + sol = solve(prob, MethodOfSteps(Tsit5()), abstol=1e-6) + @btime solve($prob, MethodOfSteps(Tsit5()), abstol=1e-6) + + @test isapprox(sol[end], zeros(N), atol=1e-3) +end -@code_warntype dnd!.f(dx0, x0,nothing,0.) -@allocated nd!(dx0, x0,nothing,0.) -@allocated dnd!(dx0, x0, p, 0.) + +## Allocations + +@inline function delayedge!(e, v_s, v_d, f::F, v_s_idx, v_d_idx, p, t) where F + hist1 = f(nothing, t-1., idxs=v_s_idx[1])::Float64 + e[1] = hist1 - v_d[1] + hist2 = f(nothing, t-1., idxs=v_d_idx[1])::Float64 + e[2] = hist2 - v_s[1] + nothing +end + +nd_delay_edge = StaticDelayEdge(f! = delayedge!, dim = 2, coupling = :fiducial) + + +begin + # (vertexp, edgep, τ) + p = (nothing, nothing, 1.) + x0 = [ones(N÷2); -ones(N÷2)] + const x1 = copy(x0) + dx0 = similar(x0) + const dx1 = similar(dx0) +end +function hh(p::Nothing, t::Float64; idxs::Int64=nothing)::Float64 + return x1[idxs] +end +@btime delayedge!(dx1,dx1,dx1,hh,1:1,1:1,p,0) + + +dnd! = network_dynamics(nd_diffusion_vertex, nd_delay_edge, g) + +@btime $dnd!.f($dx1,$x1,$hh,$p,0.) +#@btime $nd!.f($dx0,$x0,$hh,$p,0.) + +tspan = (0.,5000.) +prob = DDEProblem(dnd!, x1, hh, tspan, p)#, constant_lags=[p[end]]) +sol = solve(prob, MethodOfSteps(Tsit5()), abstol=1e-6) +#@btime solve($prob, MethodOfSteps(Tsit5()), abstol=1e-6) + +@test isapprox(sol[end], zeros(N), atol=1e-3) From 7b95aa20b572683f819b6114edf75e6eac53bb9a Mon Sep 17 00:00:00 2001 From: lindnemi Date: Thu, 21 Oct 2021 11:47:04 +0200 Subject: [PATCH 06/30] adding type information everywhere --- src/NetworkDiffEq.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/NetworkDiffEq.jl b/src/NetworkDiffEq.jl index 8f21628f7..a65a7c006 100644 --- a/src/NetworkDiffEq.jl +++ b/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 @@ -82,12 +82,12 @@ function _inner_loop!(component::DDEVertex, indices, end function _inner_loop!(component::StaticDelayEdge, indices, - dx, p, t, gd, gs, h, parallel) + dx, p, t, gd, gs, h::H, parallel) where H @nd_threads parallel for i in indices component.f!(get_edge(gd, i), get_src_vertex(gd, i), get_dst_vertex(gd, i), - h, # history function + h::H, # history function gs.s_e_idx[i], gs.d_e_idx[i], p_e_idx(p, i), @@ -143,7 +143,7 @@ function (d::NetworkDE)(dx, x, p, t) return nothing end # for DDE case -function (d::NetworkDE)(dx, x, h, 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) @@ -151,10 +151,10 @@ function (d::NetworkDE)(dx, x, h, p, t) @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, d.parallel) + dx, p, t, gd, gs, h::H, d.parallel) component_loop!(d.unique_vertices!, d.unique_v_indices, - dx, p, t, gd, gs, h, d.parallel) + dx, p, t, gd, gs, h::H, d.parallel) return nothing end From 06236f4336cd0a9597f2eb0a9f86b37d64e6805c Mon Sep 17 00:00:00 2001 From: lindnemi Date: Mon, 25 Oct 2021 16:41:28 +0200 Subject: [PATCH 07/30] minor --- test/delay_test.jl | 122 ++++++++++++++++++++------------------------- 1 file changed, 53 insertions(+), 69 deletions(-) diff --git a/test/delay_test.jl b/test/delay_test.jl index afd325aa8..f305dd0f8 100644 --- a/test/delay_test.jl +++ b/test/delay_test.jl @@ -5,8 +5,9 @@ using OrdinaryDiffEq using DelayDiffEq using BenchmarkTools -N = 10 +N = 40 g = watts_strogatz(N,4,0.) # more edges than vertices! +g = complete_graph(N) @inline function diffusionedge!(e, v_s, v_d, p, t) e[1] = v_s[1] - v_d[1] @@ -32,7 +33,7 @@ nd! = network_dynamics(nd_diffusion_vertex, nd_diffusion_edge, g) nd! = network_dynamics(nd_diffusion_vertex, nd_diffusion_edge, g) x0 = [ones(N÷2); -ones(N÷2)] - tspan = (0.,5000.) + tspan = (0.,500.) dx0 = similar(x0) @btime $nd!.f($dx0,$x0,nothing,0.) prob = ODEProblem(nd!, x0, tspan, nothing) @@ -42,100 +43,83 @@ nd! = network_dynamics(nd_diffusion_vertex, nd_diffusion_edge, g) @test isapprox(sol[end], zeros(N), atol=1e-4) end -@testset "Homogeneous DDE diffusion" begin - @inline function delayedge!(e, v_s, v_d, h_v_s, h_v_d, p, t) - e[1] = h_v_s[1] - v_d[1] +# @testset "Homogeneous DDE diffusion" begin +# @inline function delayedge!(e, v_s, v_d, h_v_s, h_v_d, p, t) +# e[1] = h_v_s[1] - v_d[1] +# nothing +# end +# nd_delay_edge = StaticDelayEdge(f! = delayedge!, dim = 1) +# +# dnd! = network_dynamics(nd_diffusion_vertex, nd_delay_edge, g) +# +# x0 = [ones(N÷2); -ones(N÷2)] +# tspan = (0.,500.) +# # constant initial conditions on interval +# h!(out, p, t) = out .= x0 +# # (vertexp, edgep, τ) +# p = (nothing, nothing, 1.) +# +# dx0 = similar(x0) +# @btime $dnd!.f($dx0,$x0,$h!,$p,0.) +# +# prob = DDEProblem(dnd!, x0, h!, tspan, p)#, constant_lags=[p[end]]) +# sol = solve(prob, MethodOfSteps(Tsit5()), abstol=1e-6) +# @btime solve($prob, MethodOfSteps(Tsit5()), abstol=1e-6) +# @test isapprox(sol[end], zeros(N), atol=1e-4) +# end + + + +@testset "New homogeneous DDE diffusion" begin + @inline function delayedge!(e, v_s, v_d, h, v_s_idx, v_d_idx, p, t) + hist1 = h(p, t - 1., idxs=v_s_idx[1]) + e[1] = hist1 - v_d[1] nothing end - nd_delay_edge = StaticDelayEdge(f! = delayedge!, dim = 1) + + nd_delay_edge = StaticDelayEdge(f! = delayedge!, dim = 1, coupling = :undirected) dnd! = network_dynamics(nd_diffusion_vertex, nd_delay_edge, g) x0 = [ones(N÷2); -ones(N÷2)] - tspan = (0.,5000.) + tspan = (0.,500.) # constant initial conditions on interval - h!(out, p, t) = out .= x0 + h(p, t; idxs=nothing) = x0[idxs] # (vertexp, edgep, τ) - p = (nothing, nothing, 1.) dx0 = similar(x0) - @btime $dnd!.f($dx0,$x0,$h!,$p,0.) + @btime $dnd!.f($dx0,$x0,$h,nothing,0.) + prob = DDEProblem(dnd!, x0, h, tspan, constant_lags=[1.]) + sol = solve(prob, MethodOfSteps(Tsit5())) + @btime solve($prob, MethodOfSteps(Tsit5())) - prob = DDEProblem(dnd!, x0, h!, tspan, p)#, constant_lags=[p[end]]) - sol = solve(prob, MethodOfSteps(Tsit5()), abstol=1e-6) - @btime solve($prob, MethodOfSteps(Tsit5()), abstol=1e-6) @test isapprox(sol[end], zeros(N), atol=1e-4) end - - -@testset "New homogeneous DDE diffusion" begin +@testset "New heterogeneous DDE diffusion" begin @inline function delayedge!(e, v_s, v_d, h, v_s_idx, v_d_idx, p, t) - hist1 = h(nothing, t-1., idxs=v_s_idx[1])::Float64 + hist1 = h(p, t - p, idxs=v_s_idx[1]) e[1] = hist1 - v_d[1] - hist2 = h(nothing, t-1., idxs=v_d_idx[1])::Float64 - e[2] = hist2 - v_s[1] nothing end - nd_delay_edge = StaticDelayEdge(f! = delayedge!, dim = 2, coupling = :fiducial) + nd_delay_edge = StaticDelayEdge(f! = delayedge!, dim = 1, coupling = :undirected) dnd! = network_dynamics(nd_diffusion_vertex, nd_delay_edge, g) x0 = [ones(N÷2); -ones(N÷2)] - tspan = (0.,5000.) + tspan = (0.,500.) # constant initial conditions on interval - h(p, t; idxs=nothing) = x0[idxs] - # (vertexp, edgep, τ) - p = (nothing, nothing, 1.) + h(p, t; idxs=nothing) = 0. + # (vertexp, edgep) + p = (nothing, collect(1:ne(g)) ./ ne(g)) dx0 = similar(x0) - @btime $dnd!.f($dx0,$x0,$h,$p,0.) + @btime $dnd!.f($dx0,$x0,$h,1.,0.) - prob = DDEProblem(dnd!, x0, h, tspan, p)#, constant_lags=[p[end]]) - sol = solve(prob, MethodOfSteps(Tsit5()), abstol=1e-6) - @btime solve($prob, MethodOfSteps(Tsit5()), abstol=1e-6) + prob = DDEProblem(dnd!, x0, h, tspan, p) # too many lags, hence no constant_lags=[1.]) + sol = solve(prob, MethodOfSteps(Tsit5()),save_everystep=false) + @btime solve($prob, MethodOfSteps(Tsit5()),save_everystep=false) @test isapprox(sol[end], zeros(N), atol=1e-3) end - - - -## Allocations - -@inline function delayedge!(e, v_s, v_d, f::F, v_s_idx, v_d_idx, p, t) where F - hist1 = f(nothing, t-1., idxs=v_s_idx[1])::Float64 - e[1] = hist1 - v_d[1] - hist2 = f(nothing, t-1., idxs=v_d_idx[1])::Float64 - e[2] = hist2 - v_s[1] - nothing -end - -nd_delay_edge = StaticDelayEdge(f! = delayedge!, dim = 2, coupling = :fiducial) - - -begin - # (vertexp, edgep, τ) - p = (nothing, nothing, 1.) - x0 = [ones(N÷2); -ones(N÷2)] - const x1 = copy(x0) - dx0 = similar(x0) - const dx1 = similar(dx0) -end -function hh(p::Nothing, t::Float64; idxs::Int64=nothing)::Float64 - return x1[idxs] -end -@btime delayedge!(dx1,dx1,dx1,hh,1:1,1:1,p,0) - - -dnd! = network_dynamics(nd_diffusion_vertex, nd_delay_edge, g) - -@btime $dnd!.f($dx1,$x1,$hh,$p,0.) -#@btime $nd!.f($dx0,$x0,$hh,$p,0.) - -tspan = (0.,5000.) -prob = DDEProblem(dnd!, x1, hh, tspan, p)#, constant_lags=[p[end]]) -sol = solve(prob, MethodOfSteps(Tsit5()), abstol=1e-6) -#@btime solve($prob, MethodOfSteps(Tsit5()), abstol=1e-6) - -@test isapprox(sol[end], zeros(N), atol=1e-3) From 872cf3d43821ca6afcf25a6e47dbb1d4192d6621 Mon Sep 17 00:00:00 2001 From: lindnemi Date: Mon, 25 Oct 2021 16:44:25 +0200 Subject: [PATCH 08/30] Static Edge promotion is not used anymore --- src/ComponentFunctions.jl | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/src/ComponentFunctions.jl b/src/ComponentFunctions.jl index f22461075..bde75c659 100644 --- a/src/ComponentFunctions.jl +++ b/src/ComponentFunctions.jl @@ -378,17 +378,6 @@ Like a static edge but with extra arguments for the history of the source and de 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 From c0ac04c6882806a7e4d1e3f4faf473c9e22f806f Mon Sep 17 00:00:00 2001 From: lindnemi Date: Mon, 25 Oct 2021 17:00:29 +0200 Subject: [PATCH 09/30] remove history from NetworkDE --- src/NetworkDiffEq.jl | 23 +++++++++++------------ src/NetworkDynamics.jl | 11 +++-------- 2 files changed, 14 insertions(+), 20 deletions(-) diff --git a/src/NetworkDiffEq.jl b/src/NetworkDiffEq.jl index a65a7c006..75de4beec 100644 --- a/src/NetworkDiffEq.jl +++ b/src/NetworkDiffEq.jl @@ -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,12 +69,13 @@ 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 component.f!(view(dx, gs.v_idx[i]), get_vertex(gd, i), get_dst_edges(gd, i), - view(history, gs.v_idx[i]), + h::H, + gs.v_idx[i], p_v_idx(p, i), t) end @@ -97,7 +98,7 @@ function _inner_loop!(component::StaticDelayEdge, indices, 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), @@ -112,7 +113,7 @@ 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 @@ -120,7 +121,6 @@ end 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 parallel::Bool # enables multithreading for the core loop end @@ -133,13 +133,12 @@ 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 @@ -166,7 +165,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 diff --git a/src/NetworkDynamics.jl b/src/NetworkDynamics.jl index 981ef5522..179fa1d25 100644 --- a/src/NetworkDynamics.jl +++ b/src/NetworkDynamics.jl @@ -149,7 +149,6 @@ end function _network_dynamics(vertices!::Union{Array{T, 1}, T}, edges!::Union{Array{U, 1}, U}, graph; - initial_history=nothing, x_prototype=zeros(1), parallel=false) where {T <: VertexFunction, U <: EdgeFunction} @@ -171,15 +170,11 @@ function _network_dynamics(vertices!::Union{Array{T, 1}, 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 hasDelay ? DDEFunction(nd!; mass_matrix = mass_matrix, syms=symbols_v) : @@ -188,7 +183,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) @@ -214,7 +209,7 @@ function _network_dynamics(vertices!::Union{Vector{T}, T}, edges!::Union{Vector{ 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) From a8abe62ae1d4b909e22211d4ee845ebca6a6c59a Mon Sep 17 00:00:00 2001 From: lindnemi Date: Mon, 25 Oct 2021 17:00:37 +0200 Subject: [PATCH 10/30] tuple syntax not necessary anymore --- src/Utilities.jl | 24 ------------------------ 1 file changed, 24 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index 5f426b491..770442421 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -73,18 +73,6 @@ end ## non-allocating but code duplication -@inline Base.@propagate_inbounds function p_v_idx(p::Tuple{T1,T2,T3}, i) where {T1 <: AbstractArray, T2, T3} - @view p[1][:, i] -end - -@inline Base.@propagate_inbounds function p_v_idx(p::Tuple{T1,T2,T3}, i) where {T1 <: AbstractArray{T4, 1} where T4, T3, T2} - p[1][i] -end - -@inline Base.@propagate_inbounds function p_v_idx(p::Tuple{T1,T2, T3}, i) where {T1, T2, T3} - p[1] -end - @inline Base.@propagate_inbounds function p_v_idx(p::Tuple{T1,T2}, i) where {T1 <: AbstractArray, T2} @view p[1][:, i] end @@ -106,18 +94,6 @@ end ## non-allocating but code duplication -@inline Base.@propagate_inbounds function p_e_idx(p::Tuple{T1,T2,T3}, i) where {T1 <: AbstractArray, T2, T3} - @view p[2][:, i] -end - -@inline Base.@propagate_inbounds function p_e_idx(p::Tuple{T1,T2,T3}, i) where {T1 <: AbstractArray{T4, 1} where T4, T3, T2} - p[2][i] -end - -@inline Base.@propagate_inbounds function p_e_idx(p::Tuple{T1,T2, T3}, i) where {T1, T2, T3} - p[2] -end - @inline Base.@propagate_inbounds function p_e_idx(p::Tuple{T1,T2}, i) where {T1, T2 <: AbstractArray} @view p[2][:, i] end From f03df36fea25e32523dfdb2488ce285b6f644120 Mon Sep 17 00:00:00 2001 From: lindnemi Date: Mon, 25 Oct 2021 17:59:21 +0200 Subject: [PATCH 11/30] fixing the tests --- src/ComponentFunctions.jl | 8 ++--- test/ComponentFunctions_test.jl | 16 +++++---- test/NetworkDiffEq_test.jl | 10 +++--- test/delay_test.jl | 62 +++++++++++++++++++-------------- test/runtests.jl | 1 + 5 files changed, 54 insertions(+), 43 deletions(-) diff --git a/src/ComponentFunctions.jl b/src/ComponentFunctions.jl index bde75c659..237dd3dcb 100644 --- a/src/ComponentFunctions.jl +++ b/src/ComponentFunctions.jl @@ -358,16 +358,16 @@ Like a static edge but with extra arguments for the history of the source and de 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) + f! = @inline (e, v_s, v_d, h, v_s_idx, v_d_idx, p, t) -> begin + @inbounds user_f!(view(e,1:dim), v_s, v_d, h, v_s_idx, v_d_idx, 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) + f! = @inline (e, v_s, v_d, h, v_s_idx, v_d_idx, p, t) -> begin + @inbounds user_f!(view(e,1:dim), v_s, v_d, h, v_s_idx, v_d_idx, p, t) @inbounds for i in 1:dim e[dim + i] = e[i] end diff --git a/test/ComponentFunctions_test.jl b/test/ComponentFunctions_test.jl index f751d9ce8..94bb8db23 100644 --- a/test/ComponentFunctions_test.jl +++ b/test/ComponentFunctions_test.jl @@ -52,8 +52,10 @@ using LinearAlgebra end @testset "StaticDelayEdge constructor" begin - f! = (e, v_s, v_d, h_v_s, h_v_d, p, t) -> begin - e .= h_v_s .- h_v_d + f! = (e, v_s, v_d, h, v_s_idx, v_d_idx, p, t) -> begin + h_v_s = h(p,t - 1., idxs = v_s_idx[1]) + h_v_d = h(p,t - 1., idxs = v_d_idx[1]) + e[1] = h_v_s - h_v_d nothing end @@ -84,11 +86,11 @@ end edir = zeros(2) efid = zeros(2) - fundir.f!(eundir, nothing, nothing, x, y, nothing, nothing) - fanti.f!(eanti, nothing, nothing, x, y, nothing, nothing) - fsym.f!(esym, nothing, nothing, x, y, nothing, nothing) - fdir.f!(edir, nothing, nothing, x, y, nothing, nothing) - ffid.f!(efid, nothing, nothing, x, y, nothing, nothing) + fundir.f!(eundir, nothing, nothing, (p,t; idxs = nothing) -> 1., [1.], [1.], nothing, 0.) + fanti.f!(eanti, nothing, nothing, (p,t; idxs = nothing) -> 1., [1.], [1.], nothing, 0.) + fsym.f!(esym, nothing, nothing, (p,t; idxs = nothing) -> 1., [1.], [1.], nothing, 0.) + fdir.f!(edir, nothing, nothing, (p,t; idxs = nothing) -> 1., [1.], [1.], nothing, 0.) + ffid.f!(efid, nothing, nothing, (p,t; idxs = nothing) -> 1., [1.], [1.], nothing, 0.) @test eundir == eanti @test eanti[1:2] == -eanti[3:4] diff --git a/test/NetworkDiffEq_test.jl b/test/NetworkDiffEq_test.jl index e696b1219..576a4ddcc 100644 --- a/test/NetworkDiffEq_test.jl +++ b/test/NetworkDiffEq_test.jl @@ -61,8 +61,8 @@ end nothing end - @inline function diffusionvertex!(dv, v, edges, h_v, p, t) - dv .= -h_v + @inline function diffusionvertex!(dv, v, edges, h, v_idx, p, t) + dv[1] = -h(p, t - 1., idxs=v_idx[1]) sum_coupling!(dv, edges) nothing end @@ -79,11 +79,11 @@ end x = rand(N) dx = similar(x) - p = (nothing, nothing, 1.0) - h(out, p, t) = (out .= x) + p = nothing + h(p, t; idxs=1) = x[idxs] # first run to compile - nd(dx, x, h, p, 0.0) # signature found in nd_DDE_Static.jl line 50 + nd(dx, x, h, p, 0.0) # function barrier, otherwise there are still allocations function alloc(nd, dx, x, h) diff --git a/test/delay_test.jl b/test/delay_test.jl index f305dd0f8..5e39b3bcf 100644 --- a/test/delay_test.jl +++ b/test/delay_test.jl @@ -3,11 +3,10 @@ using NetworkDynamics using LightGraphs using OrdinaryDiffEq using DelayDiffEq -using BenchmarkTools +#using BenchmarkTools -N = 40 -g = watts_strogatz(N,4,0.) # more edges than vertices! -g = complete_graph(N) +N = 10 +g = complete_graph(N) # more edges than vertices! @inline function diffusionedge!(e, v_s, v_d, p, t) e[1] = v_s[1] - v_d[1] @@ -25,7 +24,7 @@ end end nd_diffusion_vertex = ODEVertex(f! = diffusionvertex!, dim = 1) -nd_diffusion_edge = StaticEdge(f! = diffusionedge!, dim = 1, coupling = :antisymmetric) +nd_diffusion_edge = StaticEdge(f! = diffusionedge!, dim = 1, coupling = :undirected) nd! = network_dynamics(nd_diffusion_vertex, nd_diffusion_edge, g) @@ -33,17 +32,18 @@ nd! = network_dynamics(nd_diffusion_vertex, nd_diffusion_edge, g) nd! = network_dynamics(nd_diffusion_vertex, nd_diffusion_edge, g) x0 = [ones(N÷2); -ones(N÷2)] - tspan = (0.,500.) + tspan = (0.,100.) + dx0 = similar(x0) - @btime $nd!.f($dx0,$x0,nothing,0.) + #@btime $nd!.f($dx0,$x0,nothing,0.) prob = ODEProblem(nd!, x0, tspan, nothing) - sol = solve(prob, Tsit5(), abstol=1e-6) + sol = solve(prob, Tsit5()) - @btime solve($prob, Tsit5(), abstol=1e-6) - @test isapprox(sol[end], zeros(N), atol=1e-4) + #@btime solve($prob, Tsit5(), abstol=1e-6) + @test isapprox(sol[end], zeros(N), atol=1e-6) end -# @testset "Homogeneous DDE diffusion" begin +# @testset "Homogeneous DDE diffusion old" begin # @inline function delayedge!(e, v_s, v_d, h_v_s, h_v_d, p, t) # e[1] = h_v_s[1] - v_d[1] # nothing @@ -53,24 +53,25 @@ end # dnd! = network_dynamics(nd_diffusion_vertex, nd_delay_edge, g) # # x0 = [ones(N÷2); -ones(N÷2)] -# tspan = (0.,500.) +# tspan = (0.,100.) + # # constant initial conditions on interval # h!(out, p, t) = out .= x0 # # (vertexp, edgep, τ) # p = (nothing, nothing, 1.) # # dx0 = similar(x0) -# @btime $dnd!.f($dx0,$x0,$h!,$p,0.) +# # @btime $dnd!.f($dx0,$x0,$h!,$p,0.) # # prob = DDEProblem(dnd!, x0, h!, tspan, p)#, constant_lags=[p[end]]) # sol = solve(prob, MethodOfSteps(Tsit5()), abstol=1e-6) -# @btime solve($prob, MethodOfSteps(Tsit5()), abstol=1e-6) +# # @btime solve($prob, MethodOfSteps(Tsit5()), abstol=1e-6) # @test isapprox(sol[end], zeros(N), atol=1e-4) # end -@testset "New homogeneous DDE diffusion" begin +@testset "Homogeneous DDE diffusion" begin @inline function delayedge!(e, v_s, v_d, h, v_s_idx, v_d_idx, p, t) hist1 = h(p, t - 1., idxs=v_s_idx[1]) e[1] = hist1 - v_d[1] @@ -80,23 +81,24 @@ end nd_delay_edge = StaticDelayEdge(f! = delayedge!, dim = 1, coupling = :undirected) dnd! = network_dynamics(nd_diffusion_vertex, nd_delay_edge, g) - + dnd! = ODEFunction(dnd!.f) x0 = [ones(N÷2); -ones(N÷2)] - tspan = (0.,500.) + tspan = (0.,100.) + # constant initial conditions on interval h(p, t; idxs=nothing) = x0[idxs] # (vertexp, edgep, τ) dx0 = similar(x0) - @btime $dnd!.f($dx0,$x0,$h,nothing,0.) + #@btime $dnd!.f($dx0,$x0,$h,nothing,0.) prob = DDEProblem(dnd!, x0, h, tspan, constant_lags=[1.]) sol = solve(prob, MethodOfSteps(Tsit5())) - @btime solve($prob, MethodOfSteps(Tsit5())) + #@btime solve($prob, MethodOfSteps(Tsit5())) - @test isapprox(sol[end], zeros(N), atol=1e-4) + @test isapprox(sol[end], zeros(N), atol=1e-6) end -@testset "New heterogeneous DDE diffusion" begin +@testset "Heterogeneous DDE diffusion" begin @inline function delayedge!(e, v_s, v_d, h, v_s_idx, v_d_idx, p, t) hist1 = h(p, t - p, idxs=v_s_idx[1]) e[1] = hist1 - v_d[1] @@ -108,18 +110,24 @@ end dnd! = network_dynamics(nd_diffusion_vertex, nd_delay_edge, g) x0 = [ones(N÷2); -ones(N÷2)] - tspan = (0.,500.) + tspan = (0.,100.) + # constant initial conditions on interval h(p, t; idxs=nothing) = 0. # (vertexp, edgep) p = (nothing, collect(1:ne(g)) ./ ne(g)) dx0 = similar(x0) - @btime $dnd!.f($dx0,$x0,$h,1.,0.) + #@btime $dnd!.f($dx0,$x0,$h,1.,0.) + + prob = DDEProblem(dnd!, x0, h, tspan, p, constant_lags = p[2]) + sol = solve(prob, MethodOfSteps(Tsit5())) - prob = DDEProblem(dnd!, x0, h, tspan, p) # too many lags, hence no constant_lags=[1.]) - sol = solve(prob, MethodOfSteps(Tsit5()),save_everystep=false) - @btime solve($prob, MethodOfSteps(Tsit5()),save_everystep=false) + # for too many lags, declaring lags is very slow + #prob = DDEProblem(dnd!, x0, h, tspan, p) + #@btime solve($prob, MethodOfSteps(Tsit5()),save_everystep=false) - @test isapprox(sol[end], zeros(N), atol=1e-3) + # With inhomogenities we lose a lot of precision + # increase reltol? + @test isapprox(sol[end], zeros(N), atol=1e-2) end diff --git a/test/runtests.jl b/test/runtests.jl index 1394ac0b0..55eed74d1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,4 +11,5 @@ using Test include("inhomogeneous_test.jl") include("autodiff_test.jl") include("massmatrix_test.jl") + include("delay_test.jl") end From 0bdb6d450f7de0a1b84c5c6b14a4f875f2cfe75e Mon Sep 17 00:00:00 2001 From: lindnemi Date: Mon, 25 Oct 2021 17:59:50 +0200 Subject: [PATCH 12/30] simplify tests (minor) --- test/checkbounds_test.jl | 2 -- test/inhomogeneous_test.jl | 15 +-------------- 2 files changed, 1 insertion(+), 16 deletions(-) diff --git a/test/checkbounds_test.jl b/test/checkbounds_test.jl index 02ac4775d..41ad7d9e6 100644 --- a/test/checkbounds_test.jl +++ b/test/checkbounds_test.jl @@ -7,8 +7,6 @@ g = erdos_renyi(10,10) ov = ODEVertex(f! = (dv, v, edges, p, t) -> @inbounds( dv[1] = v[1]), dim = 1) se = StaticEdge(f! = (e, v_s, v_d, p, t) -> nothing, dim = 1, coupling = :undirected) -sde = StaticDelayEdge(se) -oe = ODEEdge(se) nd! = network_dynamics(ov,se,g) diff --git a/test/inhomogeneous_test.jl b/test/inhomogeneous_test.jl index cd1e45a73..5a7b828a0 100644 --- a/test/inhomogeneous_test.jl +++ b/test/inhomogeneous_test.jl @@ -46,22 +46,9 @@ x0[1] = pi prob_st_ver = ODEProblem(diff_network_st_ver, x0, (0.,500.)) -@time sol_st_ver = solve(prob_st_ver, Rodas4()) +sol_st_ver = solve(prob_st_ver, Rodas4()) -function test(vertex_list, edge_list, g) - nd = network_dynamics(vertex_list, edge_list, g) - - x0 = rand(nv(g)) - x0[1] = pi - - - prob_st_ver = ODEProblem(nd, x0, (0.,500.)) - @time sol_st_ver = solve(prob_st_ver, Rodas4()) - @time sol_st_ver = solve(prob_st_ver, Rodas4(), saveat=[250,500]) -end -test(vertex_list, edge_list, barabasi_albert(10,5)) - println("These dynamics should flow to π, at t=500. they are there up to $(maximum(abs.(sol_st_ver(500.) .- pi)))") @test maximum(abs.(sol_st_ver(500.) .- pi)) < 10^-7 From 2ef532580c5febed8e9f63defd65923a9e118192 Mon Sep 17 00:00:00 2001 From: lindnemi Date: Thu, 18 Nov 2021 12:21:36 +0100 Subject: [PATCH 13/30] new example with delay --- examples/Project.toml | 3 +- examples/kura_delay.jl | 94 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 96 insertions(+), 1 deletion(-) create mode 100644 examples/kura_delay.jl diff --git a/examples/Project.toml b/examples/Project.toml index 195456f2d..5169fa884 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -5,6 +5,7 @@ authors = ["Frank Hellmann ", "Michael Lindner Date: Mon, 22 Nov 2021 11:31:21 +0100 Subject: [PATCH 14/30] wrap `h`, change function signature --- src/ComponentFunctions.jl | 31 ++++++++----------------------- src/NetworkDiffEq.jl | 13 ++++++++++--- test/ComponentFunctions_test.jl | 32 ++++++++++++++++---------------- test/delay_test.jl | 8 ++++---- 4 files changed, 38 insertions(+), 46 deletions(-) diff --git a/src/ComponentFunctions.jl b/src/ComponentFunctions.jl index 237dd3dcb..1fa39e53f 100644 --- a/src/ComponentFunctions.jl +++ b/src/ComponentFunctions.jl @@ -347,30 +347,15 @@ Like a static edge but with extra arguments for the history of the source and de 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. - f! = @inline (e, v_s, v_d, h, v_s_idx, v_d_idx, p, t) -> begin - @inbounds user_f!(view(e,1:dim), v_s, v_d, h, v_s_idx, v_d_idx, p, t) - @inbounds user_f!(view(e,dim+1:2dim), v_d, v_s, h, v_d_idx, v_s_idx, p, t) - nothing - end - elseif coupling == :antisymmetric - f! = @inline (e, v_s, v_d, h, v_s_idx, v_d_idx, p, t) -> begin - @inbounds user_f!(view(e,1:dim), v_s, v_d, h, v_s_idx, v_d_idx, 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_idx, v_d_idx, p, t) -> begin - @inbounds user_f!(view(e,1:dim), v_s, v_d, h, v_s_idx, v_d_idx, p, t) - @inbounds for i in 1:dim - e[dim + i] = e[i] - end + # 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 end diff --git a/src/NetworkDiffEq.jl b/src/NetworkDiffEq.jl index 75de4beec..6c373bfe1 100644 --- a/src/NetworkDiffEq.jl +++ b/src/NetworkDiffEq.jl @@ -70,7 +70,11 @@ end function _inner_loop!(component::DDEVertex, indices, 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), @@ -85,12 +89,15 @@ end function _inner_loop!(component::StaticDelayEdge, indices, 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), - h::H, # history function - gs.s_e_idx[i], - gs.d_e_idx[i], + h_v_s, + h_v_d, p_e_idx(p, i), t) end diff --git a/test/ComponentFunctions_test.jl b/test/ComponentFunctions_test.jl index 94bb8db23..52270075a 100644 --- a/test/ComponentFunctions_test.jl +++ b/test/ComponentFunctions_test.jl @@ -52,10 +52,10 @@ using LinearAlgebra end @testset "StaticDelayEdge constructor" begin - f! = (e, v_s, v_d, h, v_s_idx, v_d_idx, p, t) -> begin - h_v_s = h(p,t - 1., idxs = v_s_idx[1]) - h_v_d = h(p,t - 1., idxs = v_d_idx[1]) - e[1] = h_v_s - h_v_d + f! = (e, v_s, v_d, h_v_s, h_v_d, p, t) -> begin + hist_v_s = h_v_s(t - 1., idxs = 1) + hist_v_d = h_v_d(t - 1., idxs = 1) + e[1] = hist_v_s - hist_v_d nothing end @@ -72,29 +72,29 @@ end @test StaticDelayEdge(f!, 1, :undirected, [:e]) isa StaticDelayEdge fundir = StaticDelayEdge(f! = f!, dim = 2, coupling = :undirected) - fanti = StaticDelayEdge(f! = f!, dim = 2, coupling = :antisymmetric) + #fanti = StaticDelayEdge(f! = f!, dim = 2, coupling = :antisymmetric) fdir = StaticDelayEdge(f! = f!, dim = 2, coupling = :directed) - fsym = StaticDelayEdge(f! = f!, dim = 2, coupling = :symmetric) + #fsym = StaticDelayEdge(f! = f!, dim = 2, coupling = :symmetric) ffid = StaticDelayEdge(f! = f!, dim = 2, coupling = :fiducial) x = rand(2) y = rand(2) eundir = zeros(4) - eanti = zeros(4) - esym = zeros(4) + #eanti = zeros(4) + #esym = zeros(4) edir = zeros(2) efid = zeros(2) - fundir.f!(eundir, nothing, nothing, (p,t; idxs = nothing) -> 1., [1.], [1.], nothing, 0.) - fanti.f!(eanti, nothing, nothing, (p,t; idxs = nothing) -> 1., [1.], [1.], nothing, 0.) - fsym.f!(esym, nothing, nothing, (p,t; idxs = nothing) -> 1., [1.], [1.], nothing, 0.) - fdir.f!(edir, nothing, nothing, (p,t; idxs = nothing) -> 1., [1.], [1.], nothing, 0.) - ffid.f!(efid, nothing, nothing, (p,t; idxs = nothing) -> 1., [1.], [1.], nothing, 0.) + fundir.f!(eundir, nothing, nothing, (t; idxs) -> 1., (t; idxs) -> 1., nothing, 0.) + #fanti.f!(eanti, nothing, nothing, (p,t; idxs = nothing) -> 1., [1.], [1.], nothing, 0.) + #fsym.f!(esym, nothing, nothing, (p,t; idxs = nothing) -> 1., [1.], [1.], nothing, 0.) + fdir.f!(edir, nothing, nothing, (t; idxs) -> 1., (t; idxs) -> 1., nothing, 0.) + ffid.f!(efid, nothing, nothing, (t; idxs) -> 1., (t; idxs) -> 1., nothing, 0.) - @test eundir == eanti - @test eanti[1:2] == -eanti[3:4] - @test esym[1:2] == esym[3:4] + #@test eundir == eanti + #@test eanti[1:2] == -eanti[3:4] + #@test esym[1:2] == esym[3:4] @test eundir[1:2] == edir @test edir == efid end diff --git a/test/delay_test.jl b/test/delay_test.jl index 5e39b3bcf..aa4c1581b 100644 --- a/test/delay_test.jl +++ b/test/delay_test.jl @@ -72,8 +72,8 @@ end @testset "Homogeneous DDE diffusion" begin - @inline function delayedge!(e, v_s, v_d, h, v_s_idx, v_d_idx, p, t) - hist1 = h(p, t - 1., idxs=v_s_idx[1]) + @inline function delayedge!(e, v_s, v_d, h_v_s, h_v_d, p, t) + hist1 = h_v_s(t - 1., idxs=1) e[1] = hist1 - v_d[1] nothing end @@ -99,8 +99,8 @@ end end @testset "Heterogeneous DDE diffusion" begin - @inline function delayedge!(e, v_s, v_d, h, v_s_idx, v_d_idx, p, t) - hist1 = h(p, t - p, idxs=v_s_idx[1]) + @inline function delayedge!(e, v_s, v_d, h_v_s, h_v_d, p, t) + hist1 = h_v_s(t - p, idxs=1) e[1] = hist1 - v_d[1] nothing end From a2cb157e171be3c09b2da4b21604aec25a2e80c2 Mon Sep 17 00:00:00 2001 From: lindnemi Date: Tue, 21 Jun 2022 12:02:19 +0200 Subject: [PATCH 15/30] adjustments to new conventions --- src/ComponentFunctions.jl | 12 ++++++------ src/NetworkDiffEq.jl | 4 ++-- test/ComponentFunctions_test.jl | 10 +--------- test/delay_test.jl | 17 ++++++++--------- 4 files changed, 17 insertions(+), 26 deletions(-) diff --git a/src/ComponentFunctions.jl b/src/ComponentFunctions.jl index 9850b3b5b..ec2c90f8b 100644 --- a/src/ComponentFunctions.jl +++ b/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. @@ -332,16 +332,16 @@ Base.@kwdef struct StaticDelayEdge{T} <: EdgeFunction 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) + 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. - 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) + 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 end diff --git a/src/NetworkDiffEq.jl b/src/NetworkDiffEq.jl index 0be95c76f..d122e9ea7 100644 --- a/src/NetworkDiffEq.jl +++ b/src/NetworkDiffEq.jl @@ -75,7 +75,7 @@ function _inner_loop!(component::DDEVertex, 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]), + component.f(view(dx, gs.v_idx[i]), get_vertex(gd, i), get_dst_edges(gd, i), h::H, @@ -93,7 +93,7 @@ function _inner_loop!(component::StaticDelayEdge, indices, 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), + component.f(get_edge(gd, i), get_src_vertex(gd, i), get_dst_vertex(gd, i), h_v_s, diff --git a/test/ComponentFunctions_test.jl b/test/ComponentFunctions_test.jl index 4caed604a..a5f6ef184 100644 --- a/test/ComponentFunctions_test.jl +++ b/test/ComponentFunctions_test.jl @@ -114,29 +114,21 @@ end @test StaticDelayEdge(f, 1, :undirected, [:e]) isa StaticDelayEdge fundir = StaticDelayEdge(; f=f, dim=2, coupling=:undirected) - #fanti = StaticDelayEdge(; f=f, dim=2, coupling=:antisymmetric) fdir = StaticDelayEdge(; f=f, dim=2, coupling=:directed) - #fsym = StaticDelayEdge(; f=f, dim=2, coupling=:symmetric) ffid = StaticDelayEdge(; f=f, dim=2, coupling=:fiducial) x = rand(2) y = rand(2) eundir = zeros(4) - #eanti = zeros(4) - #esym = zeros(4) edir = zeros(2) efid = zeros(2) fundir.f(eundir, nothing, nothing, (t; idxs) -> 1., (t; idxs) -> 1., nothing, 0.) - #fanti.f!(eanti, nothing, nothing, (p,t; idxs = nothing) -> 1., [1.], [1.], nothing, 0.) - #fsym.f!(esym, nothing, nothing, (p,t; idxs = nothing) -> 1., [1.], [1.], nothing, 0.) fdir.f(edir, nothing, nothing, (t; idxs) -> 1., (t; idxs) -> 1., nothing, 0.) ffid.f(efid, nothing, nothing, (t; idxs) -> 1., (t; idxs) -> 1., nothing, 0.) - #@test eundir == eanti - #@test eanti[1:2] == -eanti[3:4] - #@test esym[1:2] == esym[3:4] + @test eundir[1:2] == edir @test edir == efid end diff --git a/test/delay_test.jl b/test/delay_test.jl index aa4c1581b..5f40f53a2 100644 --- a/test/delay_test.jl +++ b/test/delay_test.jl @@ -1,6 +1,6 @@ using Test using NetworkDynamics -using LightGraphs +using Graphs using OrdinaryDiffEq using DelayDiffEq #using BenchmarkTools @@ -23,8 +23,8 @@ end nothing end -nd_diffusion_vertex = ODEVertex(f! = diffusionvertex!, dim = 1) -nd_diffusion_edge = StaticEdge(f! = diffusionedge!, dim = 1, coupling = :undirected) +nd_diffusion_vertex = ODEVertex(f = diffusionvertex!, dim = 1) +nd_diffusion_edge = StaticEdge(f = diffusionedge!, dim = 1, coupling = :undirected) nd! = network_dynamics(nd_diffusion_vertex, nd_diffusion_edge, g) @@ -48,7 +48,7 @@ end # e[1] = h_v_s[1] - v_d[1] # nothing # end -# nd_delay_edge = StaticDelayEdge(f! = delayedge!, dim = 1) +# nd_delay_edge = StaticDelayEdge(f = delayedge!, dim = 1) # # dnd! = network_dynamics(nd_diffusion_vertex, nd_delay_edge, g) # @@ -78,7 +78,7 @@ end nothing end - nd_delay_edge = StaticDelayEdge(f! = delayedge!, dim = 1, coupling = :undirected) + nd_delay_edge = StaticDelayEdge(f = delayedge!, dim = 1, coupling = :undirected) dnd! = network_dynamics(nd_diffusion_vertex, nd_delay_edge, g) dnd! = ODEFunction(dnd!.f) @@ -95,7 +95,7 @@ end sol = solve(prob, MethodOfSteps(Tsit5())) #@btime solve($prob, MethodOfSteps(Tsit5())) - @test isapprox(sol[end], zeros(N), atol=1e-6) + @test isapprox(sol[end], zeros(N), atol=1e-5) end @testset "Heterogeneous DDE diffusion" begin @@ -105,7 +105,7 @@ end nothing end - nd_delay_edge = StaticDelayEdge(f! = delayedge!, dim = 1, coupling = :undirected) + nd_delay_edge = StaticDelayEdge(f = delayedge!, dim = 1, coupling = :undirected) dnd! = network_dynamics(nd_diffusion_vertex, nd_delay_edge, g) @@ -127,7 +127,6 @@ end #prob = DDEProblem(dnd!, x0, h, tspan, p) #@btime solve($prob, MethodOfSteps(Tsit5()),save_everystep=false) - # With inhomogenities we lose a lot of precision - # increase reltol? + # Not sure if the failure to reach 0 is due to loss of precision or due to the problem formulation @test isapprox(sol[end], zeros(N), atol=1e-2) end From 25effce77abba17852c18fb2b79ad8dd12ce85a5 Mon Sep 17 00:00:00 2001 From: lindnemi Date: Tue, 21 Jun 2022 12:02:29 +0200 Subject: [PATCH 16/30] bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index b19b544a7..e14fd3683 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "NetworkDynamics" uuid = "22e9dc34-2a0d-11e9-0de0-8588d035468b" authors = ["Frank Hellmann , Michael Lindner "] -version = "0.7.1" +version = "0.8" [deps] Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" From ad40cd4f7950b3dd684f23885c749e48438e319e Mon Sep 17 00:00:00 2001 From: lindnemi Date: Tue, 21 Jun 2022 13:24:11 +0200 Subject: [PATCH 17/30] test on 1.8 --- .github/workflows/tests.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 5585471e9..b7a700410 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -12,7 +12,7 @@ jobs: fail-fast: false matrix: version: - - '1.6' + - '1.8' - '1' os: - ubuntu-latest From f9681b1f362d2f2cab8032d6dcd42cfe1ef2bb0f Mon Sep 17 00:00:00 2001 From: lindnemi Date: Wed, 22 Jun 2022 14:15:07 +0200 Subject: [PATCH 18/30] update docs with DDE --- docs/Project.toml | 8 +- docs/localmake.jl | 6 +- docs/make.jl | 2 +- docs/src/BasicConstructors.md | 6 +- docs/src/getting_started_with_DDEs.md | 194 -------------------------- docs/src/kuramoto_delay.md | 126 +++++++++++++++++ docs/src/parameters.md | 17 +-- 7 files changed, 142 insertions(+), 217 deletions(-) delete mode 100644 docs/src/getting_started_with_DDEs.md create mode 100644 docs/src/kuramoto_delay.md diff --git a/docs/Project.toml b/docs/Project.toml index 110fad9f6..6a3556842 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,13 +1,13 @@ [deps] +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" DelayDiffEq = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" +LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" NetworkDynamics = "22e9dc34-2a0d-11e9-0de0-8588d035468b" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622" StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" - -[compat] -Documenter = "0.26" diff --git a/docs/localmake.jl b/docs/localmake.jl index a6789dbcd..206a36e86 100644 --- a/docs/localmake.jl +++ b/docs/localmake.jl @@ -3,7 +3,6 @@ julia localmake.jl python3 -m http.server --bind localhost =# - using Pkg Pkg.activate(@__DIR__) Pkg.develop(PackageSpec(; path=dirname(@__DIR__))) # adds the package this script is called from @@ -22,6 +21,7 @@ makedocs(; sitename="NetworkDynamics", "accessing_edge_variables.md", "Tutorials" => ["Getting started" => "getting_started_with_network_dynamics.md", "Directed and weighted graphs" => "directed_and_weighted_graphs.md", - "Delay differential equations" => "getting_started_with_DDEs.md", + "Delay differential equations" => "kura_delay.md", "Heterogeneous systems" => "heterogeneous_system.md", - "Stochastic differential equations" => "SDEVertex.md"]]) + "Stochastic differential equations" => "SDEVertex.md" + ]]) diff --git a/docs/make.jl b/docs/make.jl index 70ede2cb4..df3bed647 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -16,7 +16,7 @@ makedocs(; sitename = "NetworkDynamics", "Directed and weighted graphs" => "directed_and_weighted_graphs.md", "Heterogeneous systems" => "heterogeneous_system.md", "Stochastic differential equations" => "SDEVertex.md", - "Delay differential equations" => "getting_started_with_DDEs.md"]]) + "Delay differential equations" => "kuramoto_delay.md"]]) deploydocs(; repo="github.com/PIK-ICoNe/NetworkDynamics.jl.git", devbranch="main") diff --git a/docs/src/BasicConstructors.md b/docs/src/BasicConstructors.md index 74f75811e..615bc96bd 100644 --- a/docs/src/BasicConstructors.md +++ b/docs/src/BasicConstructors.md @@ -80,7 +80,7 @@ The function then defaults to using the identity as mass matrix and `[:v for i i ### [`DDEVertex`](@ref) -If a vertex has local dynamics described by a delay differential equation (DDE) the local dynamics need to have the signature `vertexfunction!(dv, v, edges, h_v, p, t)`, where `h_v` are the history values of `v`. Then the `VertexFunction` is constructed as +If a vertex has local dynamics described by a delay differential equation (DDE) the local dynamics need to have the signature `vertexfunction!(dv, v, edges, h_v, p, t)`, where `h_v` is an automatically generated wrapper of the global history function, providing history values for`v`. Then the `VertexFunction` is constructed as ```julia DDEVertex(vertexfunction!, dim, mass_matrix, sym) @@ -150,12 +150,12 @@ In this case the function defaults to using the identity matrix as mass matrix a ### [`StaticDelayEdge`](@ref) -This constructor is used when edge variables depend on past values of the vertex variables. In this case the `edgefunction!` has to accept two additional arguments `h_v_s` and `h_v_d` that hold the history of `v_s` and `v_d`. *Static* means that the edge depends only on the dynamics of the vertices the edge is connected to and not on an internal derivative of the edge variables itself. +This constructor is used when edge variables depend on past values of the vertex variables. In this case the `edgefunction!` has to accept two additional arguments `h_v_s` and `h_v_d` that hold the (automatically generated) localised history functions of `v_s` and `v_d`. *Static* means that the edge depends only on the dynamics of the vertices the edge is connected to and not on an internal derivative of the edge variables itself. As an example for such system, we show a diffusive coupling with delay: ```julia -edgefunction! = (e, v_s, v_d, h_v_s, h_v_d, p, t) -> e .= 0.1 * (h_v_s .- v_d) +edgefunction! = (e, v_s, v_d, h_v_s, h_v_d, p, t) -> e .= 0.1 * (h_v_s(t - p) .- v_d) ``` The `EdgeFunction` object is constructed as diff --git a/docs/src/getting_started_with_DDEs.md b/docs/src/getting_started_with_DDEs.md deleted file mode 100644 index 5308053f4..000000000 --- a/docs/src/getting_started_with_DDEs.md +++ /dev/null @@ -1,194 +0,0 @@ -# Time-delayed network dynamics - -#### Topics covered in this tutorial include: - - - modeling time-delays with Delay Differential Equations (DDEs) - - handling parameters for DDEs - - time-delayed vertex variables - - time-delayed coupling - -## Delayed network diffusion - -We revisit the introductory example on [Network diffusion](@ref getting_started), this time with time-delayed vertex dynamics. In this section our goal is not to model real physical phenomena but rather to learn how time-delays can be modeled in `NetworkDynamics.jl`. - -Let $g$ be a graph with $N$ nodes and adjacency matrix $A$. Let $v = (v_1, \dots, v_n)$ be a vector of (abstract) temperatures or concentrations at each node $i = 1, \dots, N$. The dynamics of state $v_i$ at time $t$ are determined by its past value $v_i(t-\tau)$, where $\tau$ is a constant time lag and its difference with its neighbors with coupling strength $\sigma$. - -```math -\begin{aligned} -\dot v_i(t) = - v_i(t-\tau) - \sigma * \sum_{i=1}^N A_{ij} (v_i(t) - v_j(t)) -\end{aligned} -``` - -## Modeling diffusion with NetworkDynamics.jl - -The sum on the right hand side plays the role of a (discrete) gradient. If the temperature at node $i$ is higher than at its neighboring node $j$, it will decrease along that edge. - -The coupling function is the same as in the [previous example](@ref getting_started) - -```@example DDEVertex -function diffusionedge!(e, v_s, v_d, p, t) - e .= 0.1 * (v_s .- v_d) - nothing -end -nothing # hide -``` - -However, the internal vertex dynamics are now determined by the time-delayed equation $\dot v_i(t) = - v_i(t-\tau)$ and are described in the vertex function with help of the history array $h_v$ containing the past values of the vertex variables. - -```@example DDEVertex -function delaydiffusionvertex!(dv, v, edges, h_v, p, t) - dv .= -h_v - sum_coupling!(dv, edges) - nothing -end -nothing # hide -``` - -The Watts-Strogatz algorithm constructs a regular ring network with $N$ nodes connected to $k$ neighbors and a probability $p$ of rewiring links. Since $p=0$ no rewiring happens and `g` is a simple ring network. - -```@example DDEVertex -using Graphs - -### Defining a graph - -N = 10 # number of nodes -k = 4 # average degree -g = watts_strogatz(N, k, 0.0) # ring-network -nothing # hide -``` - -While the `EdgeFunction` is constructed via `StaticEdge` just as before, the vertex functions is wrapped in `DDEVertex` to account for the time-lag in the internal dynamics. `network_dynamics` then returns a `DDEFunction` compatible with the solvers of DifferentialEquations.jl. Combining those into a `DDEFunction` via `network_dynamics` is then straightforward. - -```@example DDEVertex -using NetworkDynamics - -nd_diffusion_vertex = DDEVertex(; f=delaydiffusionvertex!, dim=1) -nd_diffusion_edge = StaticEdge(; f=diffusionedge!, dim=1) - -nd = network_dynamics(nd_diffusion_vertex, nd_diffusion_edge, g) -nothing # hide -``` - -# Simulation - -When simulating a time-delayed system the initial conditions have to be specified on the whole interval $[-\tau, 0 ]$. This is usually done by specifying initial conditions at `t=0` and a history function `h` for extrapolating these values to time $t= - \tau$. For simplicity `h` is often chosen constant. When solving a DDE with DifferentialEquations.jl `h` is later overloaded to interpolate between time-steps, for details [see the docs](https://diffeq.sciml.ai/stable/tutorials/dde_example/). - -```@example DDEVertex - -const x0 = randn(N) # random initial conditions -tspan = (0.0, 10.0) # time interval - -# history function keeps the initial conditions constant and is in-place to save allocations -h(out, p, t) = (out .= x0) -nothing # hide -``` - -We still have to specify the *constant* time-lag $\tau$. At the moment this is only possible by using the ND-specific [tuple syntax](@ref parameters). For DDEs these tuples have to contain a third argument specifying the delay time. - -```@example DDEVertex -# parameters p = (vertex parameters, edge parameters, delay time) -p = (nothing, nothing, 1.0) -nothing # hide -``` - -We solve the problem on the time interval $[0, 10]$ with the `Tsit5()` algorithm, recommended for solving non-stiff problems. The `MethodOfSteps(..)` translates an `OrdinaryDiffEq.jl` solver into an efficient method for delay differential equations. `DDEProblem` addtional accepts the keyword `constant_lags` that can be useful in some situation, see [their docs](https://diffeq.sciml.ai/stable/tutorials/dde_example/) for details. - -```@example DDEVertex -using DelayDiffEq, Plots - -dde_prob = DDEProblem(nd, x0, h, tspan, p) -sol = solve(dde_prob, MethodOfSteps(Tsit5())) -plot(sol; vars=syms_containing(nd, "v"), fmt=:png) -``` - -## Bonus: Two independent diffusions - -In this extension of the first example, there are two independent diffusions on the same network with variables $x$ and $\phi$ - hence the dimension is set to `dim=2`. - -```@example DDEVertex -nd_diffusion_vertex_2 = DDEVertex(; f=delaydiffusionvertex!, dim=2, sym=[:x, :ϕ]) -nd_diffusion_edge_2 = StaticEdge(; f=diffusionedge!, dim=2) -nd_2 = network_dynamics(nd_diffusion_vertex_2, nd_diffusion_edge_2, g) -nothing # hide -``` - -The initial conditions are sampled from (squared) normal distributinos such that the first $N$ values correspond to variable `x` and the values with indices from $N+1$ to $2N$ belong to variable `ϕ`, where $x \sim \mathcal{N}(0,1)$; $ϕ \sim \mathcal{N}(0,1)^2$. - -```@example DDEVertex -const x0_2 = Vector{Float64}(vec([randn(N) .- 10 randn(N) .^ 2]')) # x ~ \mathcal{N}(0,1); ϕ ~ \mathcal{N}(0,1)^2 - -h_2(out, p, t) = (out .= x0_2) - -p = (nothing, nothing, 1.0) # p = (vertexparameters, edgeparameters, delaytime) -nothing # hide -``` - -Now we can define the `DDEProblem` and solve it. - -```@example DDEVertex -dde_prob_2 = DDEProblem(nd_2, x0_2, h_2, tspan, p) -sol_2 = solve(dde_prob_2, MethodOfSteps(Tsit5())); -plot(sol_2; legend=false, fmt=:png) -``` - -## Kuramoto model with delay - -A common modification of the [Kuramoto model](https://en.wikipedia.org/wiki/Kuramoto_model) is to include a time-lag in the coupling function. In neuroscience such a coupling is used to account for transmission delays along synapses connecting different neurons. - -Unlike in the diffusion example, the edges depend on past values of the vertex variables . For this reason the edgefunction has the history arrays of the destination and source vertices as arguments. - -```@example DDEVertex -function kuramoto_delay_edge!(e, v_s, v_d, h_v_s, h_v_d, p, t) - e[1] = p * sin(v_s[1] - h_v_d[1]) - nothing -end -nothing # hide -``` - -The dynamics of vertices in the Kuramoto model are determined by a constant rotation frequency `p`. Contributions of the edges are summed up according to their destinations. - -```@example DDEVertex -function kuramoto_vertex!(dv, v, edges, p, t) - dv[1] = p - for e in edges - dv .+= e - end - nothing -end -nothing # hide -``` - -In this case the vertex function does not depend on any past values and is simply constructed as an `ODEVertex`. Since the edges depend on the history of the vertex variables their calling signature has changed and they are handed to the delay-aware constructor `StaticDelayEdge`. They are called *Static* since don't have internal dynamical variables. - -```@example DDEVertex -kdedge! = StaticDelayEdge(; f=kuramoto_delay_edge!, dim=1) -kdvertex! = ODEVertex(; f=kuramoto_vertex!, dim=1) - -nd! = network_dynamics(kdvertex!, kdedge!, g) -nothing # hide -``` - -The remaining steps for simulating the system are the same as above. - -```@example DDEVertex -const x0_3 = randn(N) # random initial conditions -h_3(out, p, t) = (out .= x0_3) -# p = (vertexparameters, edgeparameters, delaytime) -ω = randn(N) -ω .-= sum(ω) / N # center at 0 -p = (ω, 2.0, 1.0) -tspan = (0.0, 10.0) - -dde_prob = DDEProblem(nd!, x0_3, h_3, tspan, p) - -sol = solve(dde_prob, MethodOfSteps(Tsit5())); - -### Plotting - -plot(sol; vars=syms_containing(nd, "v"), legend=false, fmt=:png) - -``` - -## Delay differential equations with algebraic constraints - -The vertex dynamics may in principle contain algebraic equations in mass matrix form. For an experimental test case have a look at the last section of this [example on github](https://github.com/pik-icone/NetworkDynamics.jl/blob/master/examples/getting_started_with_DDEs.jl). diff --git a/docs/src/kuramoto_delay.md b/docs/src/kuramoto_delay.md new file mode 100644 index 000000000..fb27b5bf9 --- /dev/null +++ b/docs/src/kuramoto_delay.md @@ -0,0 +1,126 @@ +# Network with time delays - Kuramoto model with delayed coupling + +````@example kuramoto_delay +using DelayDiffEq +using Random +using Graphs +using NetworkDynamics +using Distributions +using BenchmarkTools +```` + +A common modification of the [Kuramoto model](https://en.wikipedia.org/wiki/Kuramoto_model) is to include time-lags in the coupling function. In neuroscience this may be used to account for transmission delays along synapses connecting different neurons. + +In this tutorial we solve the system + +```math +\begin{aligned} +\dot \theta_i(t) = \omega_i + \sum_{j=1}^N \kappa_{ji} (\theta_j(t - \tau_{ji}) - \theta_i(t)) +\end{aligned} +``` +where $\kappa_{ji}$ denotes the coupling strength of the edge from $j$ to $i$, $\tau_{ji}$ its time delay and $\theta_j(t - \tau_{ji})$ the past state of its source vertex + +For a general introduction to solving delay differential equations in Julia see this [tutorial]( https://diffeq.sciml.ai/stable/tutorials/dde_example/). + +To implement this in NetworkDynamics.jl a `StaticDelayEdge` has to be defined. Such a coupling function receives two additional arguments `h_v_s` and `h_v_d`. +These are wrappers of the global history function that directly compute the history values of the local variables at the +source and destination vertex respectively. The delay time should be +passed like any other parameter. The `idxs` keyword argument of the history function can be used to access specific local variables. + +````@example kuramoto_delay +function delay_coupling(e, θ_s, θ_d, h_θ_s, h_θ_d, p, t) + τ, κ = p + hist1 = h_θ_s(t - τ; idxs=1) + e[1] = κ * sin(hist1 - θ_d[1]) + return nothing +end + +edge = StaticDelayEdge(; f=delay_coupling, dim=1, coupling=:directed); +nothing #hide +```` + +The vertex dynamics are simple ODE vertices. + +````@example kuramoto_delay +function kuramoto_vertex(dθ, θ, edges, p, t) + ω = p + dθ[1] = ω + for e in edges + dθ[1] += e[1] + end + return nothing +end + +vertex = ODEVertex(; f=kuramoto_vertex, dim=1); +nothing #hide +```` + +For this example we use a complete graph. Bear in mind however that the data structures of Network Dynamics are best suited for sparse problems and might introduce some additional overhead for dense graphs. + +````@example kuramoto_delay +N = 10 +g = SimpleDiGraph(complete_graph(N)) +nd = network_dynamics(vertex, edge, g) + +Random.seed!(1) + +ω = rand(nv(g)) # internal frequency +τ = rand(ne(g)) # time-lag per edge +κ = rand(ne(g)) # coupling strength per edge +p = (ω, [τ'; κ']) + +θ₀ = rand(Uniform(0, 2π), N); # initial conditions +nothing #hide +```` + +Define random initial history function + +````@example kuramoto_delay +const past = rand(Uniform(0, 2π), N) +h(p, t; idxs=nothing) = typeof(idxs) <: Number ? past[idxs] : past; +nothing #hide +```` + +When constructing the DDEProblem lags may be specified. +The solver will then track discontinuities arising from the evaluation of +the history function and step to each of the discontinuities. +Since each multiplication combination of the lags may be connected to a +discontinuity, this may be slow if many different lags are specified. + +````@example kuramoto_delay +prob = DDEProblem(nd, θ₀, h, (0.0, 1.0), p; constant_lags=τ) +@btime solve(prob, MethodOfSteps(Tsit5()); abstol=1e-10, reltol=1e-5); # ~50000 steps because of discontinuities +nothing #hide +```` + +We recommend to solve with lowered absolute and relative error tolerances, since the Kuramoto system is highly multistable and the simulations may else produce very different results. + + +The discontinuities arise from the initial history function and quickly get smoothed out (i.e. reduced in order) when the integration time is larger than the maximum lag. If the asymptotic behaviour is more interesting than the correct solution for a specific initial condition, it is possible to trade accuracy for computational speed by leaving the `constant_lags` undeclared. + +````@example kuramoto_delay +fast_prob = DDEProblem(nd, θ₀, h, (0.0, 1.0), p) +@btime solve(fast_prob, MethodOfSteps(BS3()); abstol=1e-10, reltol=1e-5); # ~200 steps +nothing #hide +```` + +The `MethodOfSteps` algortihm extends an ODE solver to DDEs. For an overview of available solvers consult the manual of DifferentialEquations.jl. For example, for stiff systems, such as this one, it might be beneficial to use a stiff solver such as `TRBDF2`. + +````@example kuramoto_delay +@btime solve(fast_prob, MethodOfSteps(TRBDF2()); abstol=1e-10, reltol=1e-5); +nothing #hide +```` + +Some further helpful comments for dealing within initial discontinuities in DDEs may be found in the [manual](https://jitcdde.readthedocs.io/en/stable/#dealing-with-initial-discontinuities) of the Python software JiTCDDE + +````@example kuramoto_delay +using Plots +fast_prob = DDEProblem(nd, θ₀, h, (0.0, 10.0), p) +sol = solve(fast_prob, MethodOfSteps(BS3()); abstol=1e-10, reltol=1e-5, saveat=0.01) +plot(sol; fmt=:png) +```` + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + diff --git a/docs/src/parameters.md b/docs/src/parameters.md index 43ecd6893..25d648794 100644 --- a/docs/src/parameters.md +++ b/docs/src/parameters.md @@ -12,7 +12,7 @@ nd! = network_dynamics(vertices!, edges!, graph) * When `p = (p_v, p_e)` is a Tuple of two values, then the first value will be passed to all vertices and the second to all edges. * If `p = (p_v_arr, p_e_arr)` is a Tuple of two Arrays with lengths corresponding to the number of nodes and number of edges respectively, then the edges or nodes receive only the parameter with the corresponding index. * If all nodes and/or edges have no internal parameters the value `nothing` may be passed. Using `nothing` instead of dummy parameters is usually faster, since then less data are copied. - * If you are working with delay differential equations the parameter tuple should have a third entry that specifies the delay time, `p = (p_v, p_e, delay_time)`. + Another option for specifying heterogeneous parameters is to make each `VertexFunction` a callable struct with the parameters hardcoded as fields. This approach is used in [PowerDynamics.jl](https://github.com/JuliaEnergy/PowerDynamics.jl). However it provides considerably less flexibility and interoperability with other packages. @@ -20,22 +20,15 @@ For its greater speed and flexibility in modeling we recommend to use the tuple -## Compatability with DiffEqFlux +## Compatability with specific packages -Most of the sensitivity algorithms that DiffEqFlux makes use of assume that the parameters `p` are a subtype of `AbstractArray`. Therefore they are not compatible with the tuple syntax.`TrackerAdjoint` does not have these limitations but works best on out-of-place problems. Unfortunately, `network_dynamics` returns an ODEProblem that is in-place (mutating its inputs) leading to slow performance with Tracker. +Some other packages from the Julia ecosystem, e.g. DiffEqFlux, assume that the parameters `p` are a subtype of `AbstractArray`. Therefore they are not fully compatible with the tuple syntax. -However, wrapping the function in such a way that it accepts arrays of parameters that are later pasted into the tuple syntax sidesteps these issues and enables the use of adjoint methods. Depending on the use case such a wrapper might look like this: +However, wrapping the function in such a way that it accepts arrays of parameters that are later pasted into the tuple syntax sidesteps potential issues. Depending on the use case such a wrapper might look like this: ```julia function nd_wrapper!(dx, x, p, t) nd!(dx, x, (p, nothing), t) end ``` - -`nd_wrapper!` should now work with `BacksolveAdjoint` and `InterpolatingAdjoint`. At the moment we recommend this way for combining NetworDynamics and DiffEqFlux. - -Forward mode (ForwardDiff.jl) and source-to-source (Zygote.jl) automatic differentiation is not fully-supported yet. For more detailed discussion see [this issue](https://github.com/pik-icone/NetworkDynamics.jl/issues/34). - -Further resources: - -* [DiffEqSensitivity algorithms](https://diffeq.sciml.ai/stable/analysis/sensitivity/) +At the moment we recommend this way for combining NetworkDynamics and DiffEqFlux. \ No newline at end of file From e53caa40346fba053c7220245684e48aef0eacab Mon Sep 17 00:00:00 2001 From: lindnemi Date: Wed, 22 Jun 2022 14:15:40 +0200 Subject: [PATCH 19/30] update dde test --- test/delay_test.jl | 96 +++++++++++++++++++++++----------------------- 1 file changed, 47 insertions(+), 49 deletions(-) diff --git a/test/delay_test.jl b/test/delay_test.jl index 5f40f53a2..f9f5c30e7 100644 --- a/test/delay_test.jl +++ b/test/delay_test.jl @@ -3,7 +3,6 @@ using NetworkDynamics using Graphs using OrdinaryDiffEq using DelayDiffEq -#using BenchmarkTools N = 10 g = complete_graph(N) # more edges than vertices! @@ -16,111 +15,110 @@ end @inline function diffusionvertex!(dv, v, edges, p, t) - dv[1] = 0. + dv[1] = 0.0 for e in edges dv[1] += e[1] end nothing end -nd_diffusion_vertex = ODEVertex(f = diffusionvertex!, dim = 1) -nd_diffusion_edge = StaticEdge(f = diffusionedge!, dim = 1, coupling = :undirected) +nd_diffusion_vertex = ODEVertex(; f=diffusionvertex!, dim=1) +nd_diffusion_edge = StaticEdge(; f=diffusionedge!, dim=1, coupling=:undirected) nd! = network_dynamics(nd_diffusion_vertex, nd_diffusion_edge, g) @testset "ODE diffusion" begin nd! = network_dynamics(nd_diffusion_vertex, nd_diffusion_edge, g) - x0 = [ones(N÷2); -ones(N÷2)] - tspan = (0.,100.) + x0 = [ones(N ÷ 2); -ones(N ÷ 2)] + tspan = (0.0, 100.0) dx0 = similar(x0) - #@btime $nd!.f($dx0,$x0,nothing,0.) prob = ODEProblem(nd!, x0, tspan, nothing) sol = solve(prob, Tsit5()) #@btime solve($prob, Tsit5(), abstol=1e-6) - @test isapprox(sol[end], zeros(N), atol=1e-6) + @test isapprox(sol[end], zeros(N); atol=1e-6) end -# @testset "Homogeneous DDE diffusion old" begin -# @inline function delayedge!(e, v_s, v_d, h_v_s, h_v_d, p, t) -# e[1] = h_v_s[1] - v_d[1] -# nothing -# end -# nd_delay_edge = StaticDelayEdge(f = delayedge!, dim = 1) -# -# dnd! = network_dynamics(nd_diffusion_vertex, nd_delay_edge, g) -# -# x0 = [ones(N÷2); -ones(N÷2)] -# tspan = (0.,100.) - -# # constant initial conditions on interval -# h!(out, p, t) = out .= x0 -# # (vertexp, edgep, τ) -# p = (nothing, nothing, 1.) -# -# dx0 = similar(x0) -# # @btime $dnd!.f($dx0,$x0,$h!,$p,0.) -# -# prob = DDEProblem(dnd!, x0, h!, tspan, p)#, constant_lags=[p[end]]) -# sol = solve(prob, MethodOfSteps(Tsit5()), abstol=1e-6) -# # @btime solve($prob, MethodOfSteps(Tsit5()), abstol=1e-6) -# @test isapprox(sol[end], zeros(N), atol=1e-4) -# end - @testset "Homogeneous DDE diffusion" begin @inline function delayedge!(e, v_s, v_d, h_v_s, h_v_d, p, t) - hist1 = h_v_s(t - 1., idxs=1) + hist1 = h_v_s(t - 1.0; idxs=1) e[1] = hist1 - v_d[1] nothing end - nd_delay_edge = StaticDelayEdge(f = delayedge!, dim = 1, coupling = :undirected) + nd_delay_edge = StaticDelayEdge(; f=delayedge!, dim=1, coupling=:undirected) dnd! = network_dynamics(nd_diffusion_vertex, nd_delay_edge, g) - dnd! = ODEFunction(dnd!.f) - x0 = [ones(N÷2); -ones(N÷2)] - tspan = (0.,100.) + x0 = [ones(N ÷ 2); -ones(N ÷ 2)] + tspan = (0.0, 100.0) # constant initial conditions on interval h(p, t; idxs=nothing) = x0[idxs] - # (vertexp, edgep, τ) dx0 = similar(x0) #@btime $dnd!.f($dx0,$x0,$h,nothing,0.) - prob = DDEProblem(dnd!, x0, h, tspan, constant_lags=[1.]) + prob = DDEProblem(dnd!, x0, h, tspan; constant_lags=[1.0]) + sol = solve(prob, MethodOfSteps(Tsit5())) + #@btime solve($prob, MethodOfSteps(Tsit5())) + + @test isapprox(sol[end], zeros(N); atol=1e-5) +end + + +@testset "Vertex DDE diffusion" begin + + @inline function delayvertex!(dv, v, edges, h_v, p, t) + dv[1] = h_v(t - 1.0; idxs=1) + for e in edges + dv[1] += e[1] + end + nothing + end + nd_delay_vertex = DDEVertex(; f=delayvertex!, dim=1) + + dnd! = network_dynamics(nd_delay_vertex, nd_diffusion_edge, g) + dnd! = ODEFunction(dnd!.f) + x0 = [ones(N ÷ 2); -ones(N ÷ 2)] + tspan = (0.0, 100.0) + + # constant initial conditions on interval + h(p, t; idxs=nothing) = x0[idxs] + dx0 = similar(x0) + #@btime $dnd!.f($dx0,$x0,$h,nothing,0.) + prob = DDEProblem(dnd!, x0, h, tspan; constant_lags=[1.0]) sol = solve(prob, MethodOfSteps(Tsit5())) #@btime solve($prob, MethodOfSteps(Tsit5())) - @test isapprox(sol[end], zeros(N), atol=1e-5) + @test isapprox(sol[end], zeros(N); atol=1e-5) end @testset "Heterogeneous DDE diffusion" begin @inline function delayedge!(e, v_s, v_d, h_v_s, h_v_d, p, t) - hist1 = h_v_s(t - p, idxs=1) + hist1 = h_v_s(t - p; idxs=1) e[1] = hist1 - v_d[1] nothing end - nd_delay_edge = StaticDelayEdge(f = delayedge!, dim = 1, coupling = :undirected) + nd_delay_edge = StaticDelayEdge(; f=delayedge!, dim=1, coupling=:undirected) dnd! = network_dynamics(nd_diffusion_vertex, nd_delay_edge, g) - x0 = [ones(N÷2); -ones(N÷2)] - tspan = (0.,100.) + x0 = [ones(N ÷ 2); -ones(N ÷ 2)] + tspan = (0.0, 100.0) # constant initial conditions on interval - h(p, t; idxs=nothing) = 0. + h(p, t; idxs=nothing) = 0.0 # (vertexp, edgep) p = (nothing, collect(1:ne(g)) ./ ne(g)) dx0 = similar(x0) #@btime $dnd!.f($dx0,$x0,$h,1.,0.) - prob = DDEProblem(dnd!, x0, h, tspan, p, constant_lags = p[2]) + prob = DDEProblem(dnd!, x0, h, tspan, p; constant_lags=p[2]) sol = solve(prob, MethodOfSteps(Tsit5())) # for too many lags, declaring lags is very slow @@ -128,5 +126,5 @@ end #@btime solve($prob, MethodOfSteps(Tsit5()),save_everystep=false) # Not sure if the failure to reach 0 is due to loss of precision or due to the problem formulation - @test isapprox(sol[end], zeros(N), atol=1e-2) + @test isapprox(sol[end], zeros(N); atol=1e-2) end From 9f1caf83869b993c77509a408634af7a677ff7e0 Mon Sep 17 00:00:00 2001 From: lindnemi Date: Wed, 22 Jun 2022 14:16:00 +0200 Subject: [PATCH 20/30] Literate compatible example --- examples/kuramoto_delay.jl | 96 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 96 insertions(+) create mode 100644 examples/kuramoto_delay.jl diff --git a/examples/kuramoto_delay.jl b/examples/kuramoto_delay.jl new file mode 100644 index 000000000..c07f87dd0 --- /dev/null +++ b/examples/kuramoto_delay.jl @@ -0,0 +1,96 @@ +# # Network with time delays - Kuramoto model with delayed coupling + +using DelayDiffEq +using Random +using Graphs +using NetworkDynamics +using Distributions +using BenchmarkTools + +# A common modification of the [Kuramoto model](https://en.wikipedia.org/wiki/Kuramoto_model) is to include time-lags in the coupling function. In neuroscience this may be used to account for transmission delays along synapses connecting different neurons. +# +# In this tutorial we solve the system + +# ```math +# \begin{aligned} +# \dot \theta_i(t) = \omega_i + \sum_{j=1}^N \kappa_{ji} (\theta_j(t - \tau_{ji}) - \theta_i(t)) +# \end{aligned} +# ``` +# where $\kappa_{ji}$ denotes the coupling strength of the edge from $j$ to $i$, $\tau_{ji}$ its time delay and $\theta_j(t - \tau_{ji})$ the past state of its source vertex +# +# For a general introduction to solving delay differential equations in Julia see this [tutorial]( https://diffeq.sciml.ai/stable/tutorials/dde_example/). +# +# To implement this in NetworkDynamics.jl a `StaticDelayEdge` has to be defined. Such a coupling function receives two additional arguments `h_v_s` and `h_v_d`. +# These are wrappers of the global history function that directly compute the history values of the local variables at the +# source and destination vertex respectively. The delay time should be +# passed like any other parameter. The `idxs` keyword argument of the history function can be used to access specific local variables. + +function delay_coupling(e, θ_s, θ_d, h_θ_s, h_θ_d, p, t) + τ, κ = p + hist1 = h_θ_s(t - τ; idxs=1) + e[1] = κ * sin(hist1 - θ_d[1]) + return nothing +end + +edge = StaticDelayEdge(; f=delay_coupling, dim=1, coupling=:directed); + +# The vertex dynamics are simple ODE vertices. + +function kuramoto_vertex(dθ, θ, edges, p, t) + ω = p + dθ[1] = ω + for e in edges + dθ[1] += e[1] + end + return nothing +end + +vertex = ODEVertex(; f=kuramoto_vertex, dim=1); + +# For this example we use a complete graph. Bear in mind however that the data structures of Network Dynamics are best suited for sparse problems and might introduce some additional overhead for dense graphs. + +N = 10 +g = SimpleDiGraph(complete_graph(N)) +nd = network_dynamics(vertex, edge, g) + +Random.seed!(1) + +ω = rand(nv(g)) # internal frequency +τ = rand(ne(g)) # time-lag per edge +κ = rand(ne(g)) # coupling strength per edge +p = (ω, [τ'; κ']) + +θ₀ = rand(Uniform(0, 2π), N); # initial conditions + +# Define random initial history function + +const past = rand(Uniform(0, 2π), N) +h(p, t; idxs=nothing) = typeof(idxs) <: Number ? past[idxs] : past; + +# When constructing the DDEProblem lags may be specified. +# The solver will then track discontinuities arising from the evaluation of +# the history function and step to each of the discontinuities. +# Since each multiplication combination of the lags may be connected to a +# discontinuity, this may be slow if many different lags are specified. +prob = DDEProblem(nd, θ₀, h, (0.0, 1.0), p; constant_lags=τ) +@btime solve(prob, MethodOfSteps(BS3()); abstol=1e-10, reltol=1e-5); # ~50000 steps because of discontinuities + +# We recommend to solve with lowered absolute and relative error tolerances, since the Kuramoto system is highly multistable and the simulations may else produce very different results. +# +# +# The discontinuities arise from the initial history function and quickly get smoothed out (i.e. reduced in order) when the integration time is larger than the maximum lag. If the asymptotic behaviour is more interesting than the correct solution for a specific initial condition, it is possible to trade accuracy for computational speed by leaving the `constant_lags` undeclared. + +fast_prob = DDEProblem(nd, θ₀, h, (0.0, 1.0), p) +@btime solve(fast_prob, MethodOfSteps(BS3()); abstol=1e-10, reltol=1e-5); # ~200 steps + +# The `MethodOfSteps` algortihm extends an ODE solver to DDEs. For an overview of available solvers consult the manual of DifferentialEquations.jl. For example, for stiff systems, such as this one, it might be beneficial to use a stiff solver such as `TRBDF2`. + +@btime solve(fast_prob, MethodOfSteps(TRBDF2()); abstol=1e-10, reltol=1e-5); + +# Some further helpful comments for dealing within initial discontinuities in DDEs may be found in the [manual](https://jitcdde.readthedocs.io/en/stable/#dealing-with-initial-discontinuities) of the Python software JiTCDDE + +using Plots +fast_prob = DDEProblem(nd, θ₀, h, (0.0, 10.0), p) +sol = solve(fast_prob, MethodOfSteps(BS3()); abstol=1e-10, reltol=1e-5, saveat=0.01) +plot(sol; fmt=:png) + From 5c5a58259fce54f6f4256a2eaa0ac94c26abef41 Mon Sep 17 00:00:00 2001 From: lindnemi Date: Wed, 22 Jun 2022 14:18:04 +0200 Subject: [PATCH 21/30] rename --- examples/kura_delay.jl | 94 ------------------------------------------ 1 file changed, 94 deletions(-) delete mode 100644 examples/kura_delay.jl diff --git a/examples/kura_delay.jl b/examples/kura_delay.jl deleted file mode 100644 index c3a5ff7e9..000000000 --- a/examples/kura_delay.jl +++ /dev/null @@ -1,94 +0,0 @@ -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) From aff297e9bc5584a66c1fbc6e168cc70d20b8d012 Mon Sep 17 00:00:00 2001 From: lindnemi Date: Wed, 22 Jun 2022 14:18:16 +0200 Subject: [PATCH 22/30] remove compat --- examples/Project.toml | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/examples/Project.toml b/examples/Project.toml index b36b26de5..2f4bae0b2 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -16,7 +16,4 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" - -[compat] -NetworkDynamics = "≥ 0.5" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" \ No newline at end of file From eeadf9b7d42449038e83df228fa4fc3de67e364d Mon Sep 17 00:00:00 2001 From: lindnemi Date: Wed, 22 Jun 2022 14:18:52 +0200 Subject: [PATCH 23/30] fix arguments --- src/NetworkDiffEq.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/NetworkDiffEq.jl b/src/NetworkDiffEq.jl index d122e9ea7..9c90aa2d1 100644 --- a/src/NetworkDiffEq.jl +++ b/src/NetworkDiffEq.jl @@ -73,13 +73,13 @@ function _inner_loop!(component::DDEVertex, indices, @nd_threads parallel for i in indices # Wrappers for the history function correct for global p and global idx + # should the default argument be idxs=eachindex(gs.v_idx[i]) and should we use views? 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), - h::H, - gs.v_idx[i], + h_v, p_v_idx(p, i), t) end From 1f800067c6a4f56136baaa3d5ed89caac1b5ccc3 Mon Sep 17 00:00:00 2001 From: Micha Date: Wed, 22 Jun 2022 14:19:55 +0200 Subject: [PATCH 24/30] Update tests.yml --- .github/workflows/tests.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index b7a700410..5585471e9 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -12,7 +12,7 @@ jobs: fail-fast: false matrix: version: - - '1.8' + - '1.6' - '1' os: - ubuntu-latest From c5993ef332a68516e80905a9ce0566069176c7f1 Mon Sep 17 00:00:00 2001 From: Micha Date: Wed, 22 Jun 2022 14:23:21 +0200 Subject: [PATCH 25/30] Update src/NetworkDiffEq.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Hans Würfel --- src/NetworkDiffEq.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/NetworkDiffEq.jl b/src/NetworkDiffEq.jl index 9c90aa2d1..0aa02a986 100644 --- a/src/NetworkDiffEq.jl +++ b/src/NetworkDiffEq.jl @@ -34,7 +34,7 @@ end for j in 1:length(unique_components) # Function barrier _inner_loop!(unique_components[j], unique_c_indices[j], - dx, p, t, gd, gs, h::H, parallel) + dx, p, t, gd, gs, h, parallel) end return nothing end From ba0ff1137817ae1af0739c26d7327f31274a1b8c Mon Sep 17 00:00:00 2001 From: Micha Date: Wed, 22 Jun 2022 14:24:55 +0200 Subject: [PATCH 26/30] Update NetworkDiffEq.jl --- src/NetworkDiffEq.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/NetworkDiffEq.jl b/src/NetworkDiffEq.jl index 0aa02a986..aa520631d 100644 --- a/src/NetworkDiffEq.jl +++ b/src/NetworkDiffEq.jl @@ -119,8 +119,7 @@ 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} +@Base.kwdef struct NetworkDE{G,GDB,elV,elE,TUV,TUE} unique_vertices!::TUV unique_v_indices::Vector{Vector{Int}} unique_edges!::TUE From 18efdc47f729f107f341538a49f4e3c43d5a7a85 Mon Sep 17 00:00:00 2001 From: Micha Date: Wed, 22 Jun 2022 17:43:41 +0200 Subject: [PATCH 27/30] Improve docstrings --- src/ComponentFunctions.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/ComponentFunctions.jl b/src/ComponentFunctions.jl index ec2c90f8b..852319bd6 100644 --- a/src/ComponentFunctions.jl +++ b/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, i.e. the edge gets passed a vector of `2dim` edge states, where the first `dim` elements will be presented to the dst and the second `dim` elements will be presented to src, For more details see the documentation. """ @@ -285,7 +285,6 @@ 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. - `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. @@ -302,7 +301,7 @@ end """ - StaticDilayEdge(; f, dim, coupling, sym) + StaticDelayEdge(; f, dim, coupling, sym) Like a static edge but with extra arguments for the history of the source and destination vertices. This is NOT a DDEEdge. """ @@ -329,9 +328,9 @@ Base.@kwdef struct StaticDelayEdge{T} <: EdgeFunction return new{T}(user_f, dim, coupling, sym) elseif coupling == :fiducial - 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.") + 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.") 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.") From f77f25e9f1541d4b0524f7ae02b72b40c7d9f804 Mon Sep 17 00:00:00 2001 From: lindnemi Date: Wed, 22 Jun 2022 17:38:20 +0200 Subject: [PATCH 28/30] fix test --- test/NetworkDiffEq_test.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/NetworkDiffEq_test.jl b/test/NetworkDiffEq_test.jl index 37f681d7c..6d494baa3 100644 --- a/test/NetworkDiffEq_test.jl +++ b/test/NetworkDiffEq_test.jl @@ -61,8 +61,8 @@ end nothing end - @inline function diffusionvertex!(dv, v, edges, h, v_idx, p, t) - dv[1] = -h(p, t - 1., idxs=v_idx[1]) + @inline function diffusionvertex!(dv, v, edges, h_v, p, t) + dv[1] = -h_v(t - 1., idxs=1) sum_coupling!(dv, edges) nothing end From 8691c06686239a5952696f167b65f148c7113632 Mon Sep 17 00:00:00 2001 From: Micha Date: Wed, 22 Jun 2022 17:47:42 +0200 Subject: [PATCH 29/30] Apply suggestions from code review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Hans Würfel --- src/Utilities.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Utilities.jl b/src/Utilities.jl index f398ee9b7..701201a1e 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -84,7 +84,7 @@ end ## non-allocating but code duplication -@inline Base.@propagate_inbounds function p_v_idx(p::Tuple{T1,T2}, i) where {T1<:AbstractArray,T2} +Base.@propagate_inbounds function p_v_idx(p::Tuple{T1,T2}, i) where {T1<:AbstractArray,T2} @view p[1][:, i] end @@ -105,7 +105,7 @@ end ## non-allocating but code duplication -@inline Base.@propagate_inbounds function p_e_idx(p::Tuple{T1,T2}, i) where {T1,T2<:AbstractArray} +Base.@propagate_inbounds function p_e_idx(p::Tuple{T1,T2}, i) where {T1,T2<:AbstractArray} @view p[2][:, i] end From 066d8fa30b61a889f8aed03a47c96d7718f23c5f Mon Sep 17 00:00:00 2001 From: lindnemi Date: Wed, 22 Jun 2022 17:58:03 +0200 Subject: [PATCH 30/30] hopefully lower docs build time --- docs/Project.toml | 1 - docs/src/kuramoto_delay.md | 12 +++++++----- examples/kuramoto_delay.jl | 12 +++++++----- 3 files changed, 14 insertions(+), 11 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 6a3556842..177cae6d2 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,5 +1,4 @@ [deps] -BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" DelayDiffEq = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" diff --git a/docs/src/kuramoto_delay.md b/docs/src/kuramoto_delay.md index fb27b5bf9..e9ec84d26 100644 --- a/docs/src/kuramoto_delay.md +++ b/docs/src/kuramoto_delay.md @@ -6,7 +6,6 @@ using Random using Graphs using NetworkDynamics using Distributions -using BenchmarkTools ```` A common modification of the [Kuramoto model](https://en.wikipedia.org/wiki/Kuramoto_model) is to include time-lags in the coupling function. In neuroscience this may be used to account for transmission delays along synapses connecting different neurons. @@ -58,7 +57,7 @@ nothing #hide For this example we use a complete graph. Bear in mind however that the data structures of Network Dynamics are best suited for sparse problems and might introduce some additional overhead for dense graphs. ````@example kuramoto_delay -N = 10 +N = 6 g = SimpleDiGraph(complete_graph(N)) nd = network_dynamics(vertex, edge, g) @@ -89,7 +88,8 @@ discontinuity, this may be slow if many different lags are specified. ````@example kuramoto_delay prob = DDEProblem(nd, θ₀, h, (0.0, 1.0), p; constant_lags=τ) -@btime solve(prob, MethodOfSteps(Tsit5()); abstol=1e-10, reltol=1e-5); # ~50000 steps because of discontinuities +@time solve(prob, MethodOfSteps(BS3()); abstol=1e-10, reltol=1e-5); # ~50000 steps because of discontinuities +@time solve(prob, MethodOfSteps(BS3()); abstol=1e-10, reltol=1e-5); nothing #hide ```` @@ -100,14 +100,16 @@ The discontinuities arise from the initial history function and quickly get smoo ````@example kuramoto_delay fast_prob = DDEProblem(nd, θ₀, h, (0.0, 1.0), p) -@btime solve(fast_prob, MethodOfSteps(BS3()); abstol=1e-10, reltol=1e-5); # ~200 steps +@time solve(fast_prob, MethodOfSteps(BS3()); abstol=1e-10, reltol=1e-5); # ~200 steps +@time solve(fast_prob, MethodOfSteps(BS3()); abstol=1e-10, reltol=1e-5); nothing #hide ```` The `MethodOfSteps` algortihm extends an ODE solver to DDEs. For an overview of available solvers consult the manual of DifferentialEquations.jl. For example, for stiff systems, such as this one, it might be beneficial to use a stiff solver such as `TRBDF2`. ````@example kuramoto_delay -@btime solve(fast_prob, MethodOfSteps(TRBDF2()); abstol=1e-10, reltol=1e-5); +@time solve(fast_prob, MethodOfSteps(TRBDF2()); abstol=1e-10, reltol=1e-5); +@time solve(fast_prob, MethodOfSteps(TRBDF2()); abstol=1e-10, reltol=1e-5); nothing #hide ```` diff --git a/examples/kuramoto_delay.jl b/examples/kuramoto_delay.jl index c07f87dd0..816847468 100644 --- a/examples/kuramoto_delay.jl +++ b/examples/kuramoto_delay.jl @@ -5,7 +5,6 @@ using Random using Graphs using NetworkDynamics using Distributions -using BenchmarkTools # A common modification of the [Kuramoto model](https://en.wikipedia.org/wiki/Kuramoto_model) is to include time-lags in the coupling function. In neuroscience this may be used to account for transmission delays along synapses connecting different neurons. # @@ -49,7 +48,7 @@ vertex = ODEVertex(; f=kuramoto_vertex, dim=1); # For this example we use a complete graph. Bear in mind however that the data structures of Network Dynamics are best suited for sparse problems and might introduce some additional overhead for dense graphs. -N = 10 +N = 6 g = SimpleDiGraph(complete_graph(N)) nd = network_dynamics(vertex, edge, g) @@ -73,7 +72,8 @@ h(p, t; idxs=nothing) = typeof(idxs) <: Number ? past[idxs] : past; # Since each multiplication combination of the lags may be connected to a # discontinuity, this may be slow if many different lags are specified. prob = DDEProblem(nd, θ₀, h, (0.0, 1.0), p; constant_lags=τ) -@btime solve(prob, MethodOfSteps(BS3()); abstol=1e-10, reltol=1e-5); # ~50000 steps because of discontinuities +@time solve(prob, MethodOfSteps(BS3()); abstol=1e-10, reltol=1e-5); # ~50000 steps because of discontinuities +@time solve(prob, MethodOfSteps(BS3()); abstol=1e-10, reltol=1e-5); # We recommend to solve with lowered absolute and relative error tolerances, since the Kuramoto system is highly multistable and the simulations may else produce very different results. # @@ -81,11 +81,13 @@ prob = DDEProblem(nd, θ₀, h, (0.0, 1.0), p; constant_lags=τ) # The discontinuities arise from the initial history function and quickly get smoothed out (i.e. reduced in order) when the integration time is larger than the maximum lag. If the asymptotic behaviour is more interesting than the correct solution for a specific initial condition, it is possible to trade accuracy for computational speed by leaving the `constant_lags` undeclared. fast_prob = DDEProblem(nd, θ₀, h, (0.0, 1.0), p) -@btime solve(fast_prob, MethodOfSteps(BS3()); abstol=1e-10, reltol=1e-5); # ~200 steps +@time solve(fast_prob, MethodOfSteps(BS3()); abstol=1e-10, reltol=1e-5); # ~200 steps +@time solve(fast_prob, MethodOfSteps(BS3()); abstol=1e-10, reltol=1e-5); # The `MethodOfSteps` algortihm extends an ODE solver to DDEs. For an overview of available solvers consult the manual of DifferentialEquations.jl. For example, for stiff systems, such as this one, it might be beneficial to use a stiff solver such as `TRBDF2`. -@btime solve(fast_prob, MethodOfSteps(TRBDF2()); abstol=1e-10, reltol=1e-5); +@time solve(fast_prob, MethodOfSteps(TRBDF2()); abstol=1e-10, reltol=1e-5); +@time solve(fast_prob, MethodOfSteps(TRBDF2()); abstol=1e-10, reltol=1e-5); # Some further helpful comments for dealing within initial discontinuities in DDEs may be found in the [manual](https://jitcdde.readthedocs.io/en/stable/#dealing-with-initial-discontinuities) of the Python software JiTCDDE