Refactor network data

Kuramoto-oscillators.jl

@@ -1,25 +1,27 @@
 include("src/NetworkDynamics.jl")
-using .NetworkDynamics
+import .NetworkDynamics
+ND = NetworkDynamics
 using LightGraphs
 using LinearAlgebra
 using DifferentialEquations
 
-g = barabasi_albert(500,5)
+g = barabasi_albert(50,5)
 
 #= vertex! is basically dv = sum(e_d) - sum(e_s), so basically simple diffusion with the addition
 of staticedge! and odeedge! below. =#
 
-function kuramoto_edge!(e,v_s,v_d,p,t)
-    e[1] = sin(v_s[1] - v_d[1])
+@inline function kuramoto_edge!(e,v_s,v_d,p,t)
+    e[1] = sin(v_s[1] - v_d[1]) * (1. + t * p.T_inv)
     # e[2] = sin(v_d[1] - v_s[1])
     nothing
 end
 
 struct kuramoto_parameters
     ω
+    T_inv
 end
 
-function kuramoto_vertex!(dv, v, e_s, e_d, p, t)
+@inline function kuramoto_vertex!(dv, v, e_s, e_d, p, t)
     # Note that e_s and e_d might be empty, the code needs to be able to deal
     # with this situation.
     dv .= p.ω
@@ -32,26 +34,41 @@ function kuramoto_vertex!(dv, v, e_s, e_d, p, t)
     nothing
 end
 
-odevertex = ODEVertex(f! = kuramoto_vertex!, dim = 1)
-staticedge = StaticEdge(f! = kuramoto_edge!, dim = 1)
+odevertex = ND.ODEVertex(f! = kuramoto_vertex!, dim = 1)
+staticedge = ND.StaticEdge(f! = kuramoto_edge!, dim = 1)
 
 vertexes = [odevertex for v in vertices(g)]
 edgices = [staticedge for e in edges(g)]
 
-parameters = [kuramoto_parameters(10. * randn()) for v in vertices(g)]
-append!(parameters, [kuramoto_parameters(0.) for v in edges(g)])
+parameters = [kuramoto_parameters(10. * randn(), 0.) for v in vertices(g)]
+append!(parameters, [kuramoto_parameters(0., 1. /30.) for v in edges(g)])
 
-kuramoto_network! = network_dynamics(vertexes,edgices,g)
+kuramoto_network! = ND.network_dynamics(vertexes,edgices,g)
+#
+#
+# using ForwardDiff
+# T_dual = ForwardDiff.Dual{nothing,Float64, ForwardDiff.pickchunksize(length(kuramoto_network!.f.e_int))}
+# kuramoto_network! = network_dynamics(vertexes,edgices,g; T=T_dual)
 
 x0 = randn(nv(g))
 dx = similar(x0)
 
 kuramoto_network!(dx, x0, parameters, 0.)
 
-prob = ODEProblem(kuramoto_network!, x0, (0.,5.), parameters)
+prob = ODEProblem(kuramoto_network!, x0, (0.,150.), parameters)
 
-sol = solve(prob)
+sol = solve(prob, Rodas4(autodiff=false))
+sol2 = solve(prob)
 
 using Plots
 
-plot(sol, vars=1:10:500)
+plot(sol, tspan=(147.,150.))
+plot(sol.t[end-200:end],mod2pi.(sol'[end-200:end,:]),legend=false)
+plot(sol2)
+
+using Profile
+using ProfileView
+
+@profile sol = solve(prob)
+
+ProfileView.view()

Project.toml

@@ -4,12 +4,17 @@ authors = ["Frank Hellmann <hellmann@pik-potsdam.de>, Alexander Penner <alex301p
 version = "0.1.0"
 
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+ProfileView = "c46f51b8-102a-5cf2-8d2c-8597cb0e0da7"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

REQUIRE

@@ -5,3 +5,4 @@ LinearAlgebra
 NLsolve
 Parameters
 SparseArrays
+Reexport

getting-started0.jl

@@ -0,0 +1,67 @@
+using Pkg
+Pkg.activate(".")
+using Revise
+
+struct nd_dummy{T}
+    e::AbstractArray{T}
+end
+
+function nd_dummy()
+    nd_dummy{Float64}(zeros(4))
+end
+
+function (ndd::nd_dummy{T})(x::AbstractArray{T}) where T
+    ndd.e[1] + x[1]
+end
+
+function (ndd::nd_dummy{U})(x::AbstractArray{T}) where T where U
+    println("Types not compatible")
+end
+
+
+include("src/NetworkDynamics.jl")
+using .NetworkDynamics
+using LightGraphs
+using LinearAlgebra
+using DifferentialEquations
+
+g = barabasi_albert(10,5)
+
+#= vertex! is basically dv = sum(e_d) - sum(e_s), so basically simple diffusion with the addition
+of staticedge! and odeedge! below. =#
+
+function vertex!(dv, v, e_s, e_d, p, t)
+    # Note that e_s and e_d might be empty, the code needs to be able to deal
+    # with this situation.
+    dv .= 0
+    for e in e_s
+        dv .-= e
+    end
+    for e in e_d
+        dv .+= e
+    end
+    nothing
+end
+
+odeedge! = (dl,l,v_s,v_d,p,t) -> dl .= 1000*(v_s - v_d - l)
+staticedge! = (l,v_s,v_d,p,t) -> l .= v_s - v_d
+
+# We construct the Vertices and Edges with dimension 2.
+
+odevertex = ODEVertex(vertex!,2,[1 0; 0 1],[:v,:v])
+odeedge = ODEEdge(odeedge! ,2)
+staticedge = StaticEdge(staticedge!, 2)
+
+vertexes = [odevertex for v in vertices(g)]
+edgices = [odeedge for e in edges(g)]
+
+nd! = network_dynamics(vertexes,edgices,g)
+
+x0 = rand(nd!.dim_nd)
+
+test_prob = ODEProblem(test,x0,h0,(0.,5.))
+test_sol = solve(test_prob)
+
+using Plots
+
+plot(test_sol, vars = 21:40)

src/NetworkDynamics.jl

@@ -1,34 +1,24 @@
 module NetworkDynamics
 
+using Reexport
+
+include("Functions.jl")
+@reexport using .NDFunctions
+
 include("nd_ODE_ODE.jl")
-using .nd_ODE_ODE_mod
-export nd_ODE_ODE
+@reexport using .nd_ODE_ODE_mod
 
 include("nd_ODE_Static.jl")
-using .nd_ODE_Static_mod
-export nd_ODE_Static
-export StaticEdgeFunction
+@reexport using .nd_ODE_Static_mod
 
 include("nd_DDE_DDE.jl")
-using .nd_DDE_DDE_mod
-export nd_DDE_DDE
-
-include("Functions.jl")
-using .NDFunctions
-export StaticVertex
-export StaticEdge
-export ODEVertex
-export ODEEdge
-export VertexFunction
-export EdgeFunction
-export DDEVertex
-export DDEEdge
+@reexport using .nd_DDE_DDE_mod
 
 include("Utilities.jl")
-using .Utilities
-export RootRhs
-export find_valid_ic
-export syms_containing
+@reexport using .Utilities
+
+include("NetworkStructures.jl")
+@reexport using .NetworkStructures
 
 export network_dynamics
 
@@ -40,93 +30,56 @@ using DifferentialEquations
 #= network_dynamics: The Main Constructor of the Package. It takes Arrays of Vertex- and Edgefunction + a graph and
 spits out an ODEFunction or DDEFunction. Others still need to be implemented. =#
 
-function network_dynamics(vertices!::Array{VertexFunction}, edges!::Array{EdgeFunction}, graph)
-    @assert length(vertices!) = length(vertices(graph))
-    for i in 1:length(vertices!)
-        if typeof(vertices![i]) == DDEVertex
-            for i in 1:length(vertices!)
-                vertices![i] = DDEVertex(vertices![i])
-            end
-            break
-        end
-    end
-    for i in 1:length(edges!)
-        if typeof(edges![i]) == DDEEdge
-            for i in 1:length(edges!)
-                edges![i] = DDEEdge(edges![i])
-            end
-            break
-        end
-    end
-    network_dynamics(vertices!,edges!,graph)
-end
-
 function network_dynamics(vertices!::Array{ODEVertex,1}, edges!::Array{StaticEdge,1}, graph)
+    @assert length(vertices!) == length(vertices(graph))
+    @assert length(edges!) == length(edges(graph))
 
-    vertex_functions = [v.f! for v in vertices!]
     dim_v = [v.dim for v in vertices!]
-    edge_functions = [e.f! for e in edges!]
     dim_e = [e.dim for e in edges!]
     dim_nd = sum(dim_v)
 
-    symbols = [Symbol(vertices![i].sym[j],"_",i) for i in 1:length(vertices!) for j in 1:dim_v[i]]
+    e_int = zeros(sum(dim_e))
+
+    graph_stucture = create_graph_structure(graph, dim_v, dim_e, e_int)
 
-    nd! = nd_ODE_Static(vertex_functions, edge_functions, graph, dim_v, dim_e)
+    nd! = nd_ODE_Static(vertices!, edges!, graph, graph_stucture)
 
     # Construct mass matrix
-    mm_array = [v.mass_matrix for v in vertices!]
-    if all([mm == I for mm in mm_array])
-        mass_matrix = I
-    else
-        mass_matrix = sparse(1.0I,dim_nd,dim_nd)
-        for i in 1:length(vertex_functions)
-            if vertices![i].mass_matrix != I
-                mass_matrix[nd!.v_idx[i],nd!.v_idx[i]] .= vertices![i].mass_matrix
-            end
-        end
-    end
+    mass_matrix = construct_mass_matrix([v.mass_matrix for v in vertices!], dim_nd, graph_stucture)
+
+    symbols = [Symbol(vertices![i].sym[j],"_",i) for i in 1:length(vertices!) for j in 1:dim_v[i]]
 
     ODEFunction(nd!,mass_matrix = mass_matrix,syms=symbols)
 end
 
+
 function network_dynamics(vertices!::Array{ODEVertex}, edges!::Array{ODEEdge}, graph)
+    @assert length(vertices!) == length(vertices(graph))
+    @assert length(edges!) == length(edges(graph))
 
-    vertex_functions = [v.f! for v in vertices!]
     dim_v = [v.dim for v in vertices!]
-    edge_functions = [e.f! for e in edges!]
     dim_e = [e.dim for e in edges!]
-    dim_nd = sum(dim_e) + sum(dim_v)
+    dim_nd = sum(dim_v)
 
-    vsymbols = [Symbol(vertices![i].sym[j],"_",i) for i in 1:length(vertices!) for j in 1:dim_v[i]]
-    esymbols = [Symbol(edges![i].sym[j],"_",i) for i in 1:length(edges!) for j in 1:dim_e[i]]
-    symbols = append!(vsymbols,esymbols)
+    e_int = zeros(sum(dim_e))
 
-    nd! = nd_ODE_ODE(vertex_functions, edge_functions, graph, dim_v, dim_e)
+    graph_stucture = create_graph_structure(graph, dim_v, dim_e, e_int)
+
+    nd! = nd_ODE_ODE(vertices!, edges!, graph, graph_stucture)
 
     # Construct mass matrix
-    mmv_array = [v.mass_matrix for v in vertices!]
-    mme_array = [e.mass_matrix for e in edges!]
-    if all([mm == I for mm in mmv_array]) && all([mm == I for mm in mme_array])
-        mass_matrix = I
-    else
-        mass_matrix = sparse(1.0I,dim_nd,dim_nd)
-        for i in 1:length(vertex_functions)
-            if vertices![i].mass_matrix != I
-                mass_matrix[nd!.v_idx[i],nd!.v_idx[i]] .= vertices![i].mass_matrix
-            end
-        end
-        for i in 1:length(edge_functions)
-            if edges![i].mass_matrix != I
-                mass_matrix[nd!.e_x_idx[i],nd!.e_x_idx[i]] = edges![i].mass_matrix
-            end
-        end
-    end
+    mass_matrix = construct_mass_matrix([v.mass_matrix for v in vertices!], [e.mass_matrix for e in edges!], dim_nd, graph_stucture)
+
+    vsymbols = [Symbol(vertices![i].sym[j],"_",i) for i in 1:length(vertices!) for j in 1:dim_v[i]]
+    esymbols = [Symbol(edges![i].sym[j],"_",i) for i in 1:length(edges!) for j in 1:dim_e[i]]
+    symbols = vcat(vsymbols,esymbols)
 
-    ODEFunction(nd!,mass_matrix = mass_matrix,syms = symbols)
+    ODEFunction(nd!,mass_matrix = mass_matrix,syms=symbols)
 end
 
 function network_dynamics(vertices!::Array{DDEVertex}, edges!::Array{DDEEdge}, graph)
-
+    @assert 1==0
+    # This isn't really ready for testing/production yet.
     vertex_functions = [v.f! for v in vertices!]
     dim_v = [v.dim for v in vertices!]
     edge_functions = [e.f! for e in edges!]
@@ -162,7 +115,6 @@ function network_dynamics(vertices!::Array{DDEVertex}, edges!::Array{DDEEdge}, g
 
     DDEFunction(nd!,mass_matrix = mass_matrix, syms = symbols)
 end
-end
 
 
 function network_dynamic(vertices!,  edges!, graph)
@@ -182,3 +134,44 @@ function network_dynamic(vertices!,  edges!, graph)
     end
     network_dynamic(va!,  ea!, graph)
 end
+
+function network_dynamics(vertices!::Array{VertexFunction}, edges!::Array{EdgeFunction}, graph)
+    @assert length(vertices!) == length(vertices(graph))
+    @assert length(edges!) == length(edges(graph))
+
+    contains_delays = false
+    contains_stochastic = false
+    contains_dyn_edge = false
+
+    for v in vertices!
+        if typeof(v) == DDEVertex
+            contains_delays = true
+        end
+        # if typeof(v) == SDEVertex
+        #     contains_stochastic = true
+        # end
+    end
+    for e in edges!
+        if typeof(e) == DDEEdge
+            contains_delays = true
+        end
+        if typeof(e) == ODEEdge
+            contains_dyn_edge = true
+        end
+        # if typeof(e) == SDEEdge
+        #     contains_stochastic = true
+        # end
+    end
+
+    # ToDo... more logic about what to construct when... A lot of this should
+    # be covered by the casts.
+    if contains_delay && contains_stochastic
+        println("Stochasticity and delay are not supported together.")
+        return nothing
+    elseif contains_delay
+        return network_dynamics(Array{DDEVertex}(vertices!),Array{DDEEdge}(edges!),graph)
+    end
+    nothing
+end
+
+end # module