Merge pull request #21 from FHell/multithreading

Multithreading

.gitignore

@@ -7,3 +7,4 @@ Manifest.toml
 scratch.jl
 docs/build
 bugs
+design_experiments

Project.toml

@@ -1,7 +1,7 @@
 name = "NetworkDynamics"
 uuid = "22e9dc34-2a0d-11e9-0de0-8588d035468b"
-authors = ["Frank Hellmann <[email protected]>, Alexander Penner <[email protected]>"]
-version = "0.2.0"
+authors = ["Frank Hellmann <[email protected]>, Michael Lindner <[email protected]>, Alexander Penner <[email protected]>"]
+version = "0.3.0"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
@@ -12,8 +12,8 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-julia = "1"
 DiffEqBase = "^6.9.4"
 LightGraphs = "^1.3.0"
 NLsolve = "^4.2.0"
 Reexport = "^0.2.0"
+julia = "1"

benchmarking/Project.toml

@@ -1,9 +1,11 @@
 name = "NetworkDynamics-Benchmarking"
-author = ["Frank Hellmann <[email protected]>"]
+authors = ["Frank Hellmann <[email protected]>",
+           "Michael Lindner <[email protected]>"]
 
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Debugger = "31a5f54b-26ea-5ae9-a837-f05ce5417438"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
@@ -14,5 +16,6 @@ Profile = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
 ProfileView = "c46f51b8-102a-5cf2-8d2c-8597cb0e0da7"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

benchmarking/multithreading_benchmarks.jl

@@ -0,0 +1,101 @@
+# call from shell: env JULIA_NUM_THREADS=1 julia mutltithreading_benchmarks.jl
+
+using NetworkDynamics
+using LightGraphs
+using BenchmarkTools
+
+### Defining a graph
+
+N = 100 # number of nodes
+k = 20  # average degree
+g = barabasi_albert(N, k) # a little more exciting than a bare random graph
+
+
+### Functions for edges and vertices
+
+@inline Base.@propagate_inbounds function diffusionedge!(e, v_s, v_d, p, t)
+    # usually e, v_s, v_d are arrays, hence we use the broadcasting operator .
+    e .= p * (v_s .- v_d)
+    nothing
+end
+
+@inline Base.@propagate_inbounds function diffusionvertex!(dv, v, e_s, e_d, p, t)
+    # usually v, e_s, e_d are arrays, hence we use the broadcasting operator .
+    dv .= p
+    # edges for which v is the source
+    for e in e_s
+        dv .-= e
+    end
+    # edges for which v is the destination
+    for e in e_d
+        dv .+= e
+    end
+    nothing
+end
+
+### Constructing the network dynamics
+
+nd_diffusion_vertex = ODEVertex(f! = diffusionvertex!, dim = 1)
+nd_diffusion_edge = StaticEdge(f! = diffusionedge!, dim = 1)
+
+
+### Benchmarking
+println("Number of Threads: ", haskey(ENV, "JULIA_NUM_THREADS") ? ENV["JULIA_NUM_THREADS"] : "1")
+println("\nBenchmarking ODE_Static...")
+p  = (collect(1:nv(g))./nv(g) .- 0.5, .5 .* ones(ne(g)))
+println("\nsingle-threaded...")
+nd = network_dynamics(nd_diffusion_vertex, nd_diffusion_edge, g, parallel=false)
+display(@benchmark nd(randn(N), randn(N), p , 0,))
+println("\nparallel...")
+nd = network_dynamics(nd_diffusion_vertex, nd_diffusion_edge, g, parallel=true)
+display(@benchmark nd(randn(N), randn(N), p , 0,))
+#=
+### Simulation
+using OrdinaryDiffEq
+
+x0 = randn(N)
+
+ode_prob = ODEProblem(nd, x0, (0., 4.), p)
+sol = solve(ode_prob, Tsit5());
+
+### Plotting
+using Plots
+plot(sol, vars = syms_containing(nd, "v"))
+=#
+### Redefinition
+
+@inline Base.@propagate_inbounds function diffusion_dedge!(de, e, v_s, v_d, p, t)
+    de .= 100. * (p * sin.(v_s - v_d) - e)
+    nothing
+end
+
+nd_diffusion_dedge = ODEEdge(f! = diffusion_dedge!, dim = 1)
+
+
+
+p  = (collect(1:nv(g))./nv(g) .- 0.5, 5 .* ones(ne(g)))
+
+### Benchmarking
+println("\nBenchmarking ODE_ODE...")
+println("\nsingle-threaded...")
+nd_ode = network_dynamics(nd_diffusion_vertex, nd_diffusion_dedge, g, parallel=false)
+display(@benchmark nd_ode(randn(nv(g) + ne(g)), randn(nv(g) + ne(g)), p , 0,))
+println("\nparallel...")
+nd_ode = network_dynamics(nd_diffusion_vertex, nd_diffusion_dedge, g, parallel=true)
+display(@benchmark nd_ode(randn(nv(g) + ne(g)), randn(nv(g) + ne(g)), p , 0,))
+
+
+println("")
+#=
+### Simulation
+
+x0 = randn(nv(g) + ne(g))
+
+ode_ode_prob = ODEProblem(nd_ode, x0, (0., 2.), p)
+sol_ode = solve(ode_ode_prob, Tsit5());
+
+### Plotting
+
+plot(sol_ode, vars = syms_containing(nd, "v"))
+
+=#

docs/localmake.jl

@@ -11,6 +11,7 @@ makedocs(
     pages = [
     "General" => "index.md",
     "BasicConstructors.md",
+	"Multithreading.md",
     "Library.md",
     "Tutorials" => [
 		"Getting started" => "getting_started_with_network_dynamics.md",

docs/make.jl

@@ -13,22 +13,16 @@ makedocs(
     pages = [
     "General" => "index.md",
     "BasicConstructors.md",
+	"Multithreading.md",
     "Library.md",
 	"Tutorials" => [
 		"Getting started" => "getting_started_with_network_dynamics.md",
 		"Directed and weighted graphs" => "directed_and_weighted_graphs.md",
     ]
     ])
 
-# Documenter automatically deploys documentation to gh-pages.
-# However at the moment there is a bug when GITHUB_TOKEN is used as
-# authentication method by the GithubActions worker: The html pages build is
-# not triggered automatically and while the docs are updated in gh-pages branch
-# they are not visible on github.io. A workaround is to manually trigger the
-# html build by commiting anything to the gh-pages branch. Alternatively the
-# SSH_DEPLOY_KEY can be used as authentication method for the worker.
+# To trigger the Github pages build, one initally manual commit to gh-pages is
+# necessary
 deploydocs(
     repo = "github.com/FHell/NetworkDynamics.jl.git",
-#   no longer necessary
-#   deploy_config = Documenter.GitHubActions(), # this should work, but it's strange
 )

docs/src/BasicConstructors.md

@@ -150,6 +150,22 @@ nd = network_dynamics(vertices!::Array{VertexFunction},
 nd(dx, x, p, t)
 ```
 
+If all vertices, respectively edges share the same `VertexFunction` or `EdgeFunction`, than `network_dynamics` can be conveniently called with these functions as arguments.
+
+```julia
+nd = network_dynamics(vertexfunction!::VertexFunction,
+                      edgefunction!::EdgeFunction, g)
+
+```
+
+The optional keyword argument `parallel` is `false` by default and can be set to `true` if a multi-threaded `ODEFunction` should be used. This may significantly improve performance on multi-core machines, for more details see section [Multi-Threading](@ref).
+
+```julia
+nd = network_dynamics(vertexfunction!::VertexFunction,
+                      edgefunction!::EdgeFunction, g, parallel=true)
+
+```
+
 ### Example
 
  Let's look at an example. First we define our graph as well as the differential systems connected to its vertices and edges:

docs/src/Multithreading.md

@@ -0,0 +1,21 @@
+# Multi-Threading
+
+Since version `0.3.0` multi-threading via the `Threads.@threads` macro is possible. This allows julia to integrate different nodes and edges in different threads, and leads to significant performance gains on parallel architectures. To enable multi-threading call `network_dynamics` with the keyword argument `parallel=true`.
+
+```julia
+network_dynamics(vertices!, edges!, graph, parallel=true)
+```
+
+In order for this to take effect, multiple threads have to be available. This is achieved by setting the environment variable `JULIA_NUM_THREADS` **before** starting Julia.  To start Julia from a bash shell and with 4 threads use:
+```
+$ env JULIA_NUM_THREADS=4 julia
+```
+
+If you are using `Juno` for the `Atom` text editor `JULIA_NUM_THREADS` is set to the number of physical cores of your processor by default. This is also the number of threads we recommend to use.
+
+!!! note
+    The thread handling causes a small overhead in the order of *20 μs* per call to the ODE function which might impair performance on small networks (5 to 20 nodes) or on single core machines. In theses cases `network_dynamics` can be called without any additional arguments, since `parallel` defaults to `false`.
+
+
+
+For more information on setting environment varibales see the [Julia documentation](https://docs.julialang.org/en/v1/manual/environment-variables/index.html#JULIA_NUM_THREADS-1).

docs/src/directed_and_weighted_graphs.md

@@ -1,5 +1,7 @@
 # Neurodynamic model of synchronization in the human brain
 
+ An `IJulia` [notebook](https://github.com/FHell/NetworkDynamics.jl/tree/master/examples) corresponding to this tutorial is available on GitHub.
+
 #### Topics covered in this tutorial include:
  * constructing a directed, weighted graph from data
  * some useful macros

docs/src/getting_started_with_network_dynamics.md

@@ -1,6 +1,6 @@
 # [Network diffusion](@id getting_started)
 
-This introductory example explains the use of the basic types and constructors in NetworkDynamics.jl by modeling a simple diffusion on an undirected network.
+This introductory example explains the use of the basic types and constructors in NetworkDynamics.jl by modeling a simple diffusion on an undirected network. A corresponding `IJulia` [notebook](https://github.com/FHell/NetworkDynamics.jl/tree/master/examples) is available on GitHub.
 
 ## Theoretical background
 
@@ -63,7 +63,7 @@ The reason for this syntax is found in the `LightGraphs.jl` package on which `Ne
 
 LightGraphs.jl implements edges as pairs of node indices `i` and `j`. [Pairs](https://docs.julialang.org/en/v1/base/collections/#Base.Pair) are basic julia data types consisting of two fixed elements, definded by writing `i => j`. In directed graphs these pairs additionally represent the direction in which the edge is pointing (from the first two the second element). In undirected graphs every edge is represent by only a single pair `i => j` if index `i` is smaller than index `j` and `j => i` otherwise. Hence, even for undirected graphs every edge has an abstract direction, specified by the pair of indices of the attached nodes.
 
-A LightGraphs.jl user, who is only interested in undirected graphs, usually does not have to deal with this abstract directionality. However since NetworkDynamics.jl is interfacing directly to the underlying graph objects, we have to keep in mind that since edges are represented by pairs every edge has an *abstract direction* and thus a source and a destination.
+A LightGraphs.jl user, who is only interested in undirected graphs, usually does not have to deal with this abstract directionality. However since NetworkDynamics.jl is interfacing directly to the underlying graph objects, we have to keep in mind that every edge has an *abstract direction* and thus a source and a destination.
 
 In the diffusion example the coupling terms have to be modified accordingly. Assume node $i$ is connected to node $j$ by an undirected edge.
 

docs/src/index.md

@@ -28,7 +28,3 @@ DifferentialEquations.jl calling syntax.
 nd = network_dynamics(vertices!::Array{VertexFunction},  edges!::Array{EdgeFunction}, g)
 nd(dx, x, p, t)
 ```
-
-
-
-This page is still in development. :)

src/NetworkDynamics.jl

@@ -35,7 +35,7 @@ constructors.
 """
 function collect_ve_info(vertices!, edges!, graph)
     if vertices! isa Array
-        @assert length(vertices!) == length(vertices(graph))
+        @assert length(vertices!) == nv(graph)
         v_dims = [v.dim for v in vertices!]
         symbols_v = [Symbol(vertices![i].sym[j],"_",i) for i in 1:length(vertices!) for j in 1:v_dims[i]]
         mmv_array = [v.mass_matrix for v in vertices!]
@@ -46,7 +46,7 @@ function collect_ve_info(vertices!, edges!, graph)
     end
 
     if edges! isa Array
-        @assert length(edges!) == length(edges(graph))
+        @assert length(edges!) == ne(graph)
         e_dims = [e.dim for e in edges!]
         symbols_e = [Symbol(edges![i].sym[j],"_",i) for i in 1:length(edges!) for j in 1:e_dims[i]]
         if eltype(edges!) <: StaticEdge
@@ -68,14 +68,29 @@ function collect_ve_info(vertices!, edges!, graph)
 end
 
 """
-    network_dynamics(vertices!, edges!, g)
+    network_dynamics(vertices!, edges!, g; parallel = false)
 
 Assembles the the dynamical equations of the network problem into an `ODEFunction`
 compatible with the `DifferentialEquations.jl` solvers. Takes as arguments an array
 of VertexFunctions **`vertices!`**, an array of EdgeFunctions **`edges!`** and a
-`LightGraph.jl` object **`g`**.
+`LightGraph.jl` object **`g`**. The optional argument `parallel` is a boolean
+value that denotes if the central loop should be executed in parallel.
 """
-function network_dynamics(vertices!::Union{Array{T, 1}, T}, edges!::Union{Array{U, 1}, U}, graph; x_prototype=zeros(1)) where {T <: ODEVertex, U <: StaticEdge}
+function network_dynamics(vertices!::Union{Array{T, 1}, T}, edges!::Union{Array{U, 1}, U},
+                          graph; x_prototype=zeros(1), parallel=false) where {T <: ODEVertex, U <: StaticEdge}
+    if parallel
+        haskey(ENV, "JULIA_NUM_THREADS") &&
+        parse(Int, ENV["JULIA_NUM_THREADS"]) > 1 ? nothing :
+        print("Warning: You are using multi-threading with only one thread ",
+        "available to Julia. Consider re-starting Julia with the environment ",
+        "variable JULIA_NUM_THREADS set to the number of physical cores of your CPU.")
+    else
+        haskey(ENV, "JULIA_NUM_THREADS") && parse(Int, ENV["JULIA_NUM_THREADS"])  > 1 ?
+        print("Your instance of Julia has more than one thread available for ",
+        "executing code. Consider calling network_dynamics with the keyword ",
+        "parallel=true.") : nothing
+    end
+
     v_dims, e_dims, symbols_v, symbols_e, mmv_array, mme_array = collect_ve_info(vertices!, edges!, graph)
 
     v_array = similar(x_prototype, sum(v_dims))
@@ -87,15 +102,27 @@ function network_dynamics(vertices!::Union{Array{T, 1}, T}, edges!::Union{Array{
 
     graph_data = GraphData(v_array, e_array, graph_stucture)
 
-    nd! = nd_ODE_Static(vertices!, edges!, graph, graph_stucture, graph_data)
-
+    nd! = nd_ODE_Static(vertices!, edges!, graph, graph_stucture, graph_data, parallel)
     mass_matrix = construct_mass_matrix(mmv_array, graph_stucture)
 
     ODEFunction(nd!; mass_matrix = mass_matrix, syms=symbols)
 end
 
 
-function network_dynamics(vertices!::Union{Array{T, 1}, T}, edges!::Union{Array{U, 1}, U}, graph; x_prototype=zeros(1)) where {T <: ODEVertex, U <: ODEEdge}
+function network_dynamics(vertices!::Union{Array{T, 1}, T}, edges!::Union{Array{U, 1}, U}, graph; x_prototype=zeros(1), parallel=false) where {T <: ODEVertex, U <: ODEEdge}
+    if parallel
+        haskey(ENV, "JULIA_NUM_THREADS") &&
+        parse(Int, ENV["JULIA_NUM_THREADS"]) > 1 ? nothing :
+        println("Warning: You are using multi-threading with only one thread ",
+        "available to Julia. Consider re-starting Julia with the environment ",
+        "variable JULIA_NUM_THREADS set to the number of physical cores of your CPU.")
+    else
+        haskey(ENV, "JULIA_NUM_THREADS") && parse(Int, ENV["JULIA_NUM_THREADS"])  > 1 ?
+        println("Info: Your instance of Julia has more than one thread available for ",
+        "executing code. Consider calling network_dynamics with the keyword ",
+        "parallel=true.") : nothing
+    end
+
     v_dims, e_dims, symbols_v, symbols_e, mmv_array, mme_array = collect_ve_info(vertices!, edges!, graph)
 
     x_array = similar(x_prototype, sum(v_dims) + sum(e_dims))
@@ -109,14 +136,14 @@ function network_dynamics(vertices!::Union{Array{T, 1}, T}, edges!::Union{Array{
 
     graph_data = GraphData(v_array, e_array, graph_stucture)
 
-    nd! = nd_ODE_ODE(vertices!, edges!, graph, graph_stucture, graph_data)
+    nd! = nd_ODE_ODE(vertices!, edges!, graph, graph_stucture, graph_data, parallel)
 
     mass_matrix = construct_mass_matrix(mmv_array, mme_array, graph_stucture)
 
     ODEFunction(nd!; mass_matrix = mass_matrix, syms=symbols)
 end
 
-function network_dynamics(vertices!,  edges!, graph)
+function network_dynamics(vertices!,  edges!, graph; parallel=false)
     try
         Array{VertexFunction}(vertices!)
     catch err
@@ -133,12 +160,12 @@ function network_dynamics(vertices!,  edges!, graph)
     end
     va! = Array{VertexFunction}(vertices!)
     ea! = Array{EdgeFunction}(edges!)
-    network_dynamics(va!,  ea!, graph)
+    network_dynamics(va!,  ea!, graph, parallel = parallel)
 end
 
-function network_dynamics(vertices!::Array{VertexFunction}, edges!::Array{EdgeFunction}, graph)
-    @assert length(vertices!) == length(vertices(graph))
-    @assert length(edges!) == length(edges(graph))
+function network_dynamics(vertices!::Array{VertexFunction}, edges!::Array{EdgeFunction}, graph; parallel=false)
+    @assert length(vertices!) == nv(graph)
+    @assert length(edges!) == ne(graph)
 
     contains_dyn_edge = false
 
@@ -149,20 +176,20 @@ function network_dynamics(vertices!::Array{VertexFunction}, edges!::Array{EdgeFu
     end
 
     if contains_dyn_edge
-        return network_dynamics(Array{ODEVertex}(vertices!),Array{ODEEdge}(edges!), graph)
+        return network_dynamics(Array{ODEVertex}(vertices!),Array{ODEEdge}(edges!), graph, parallel = parallel)
     else
-        return network_dynamics(Array{ODEVertex}(vertices!),Array{StaticEdge}(edges!), graph)
+        return network_dynamics(Array{ODEVertex}(vertices!),Array{StaticEdge}(edges!), graph, parallel = parallel)
     end
     nothing
 end
 
 """
 Allow initializing StaticEdgeFunction for Power Dynamics
 """
-function StaticEdgeFunction(vertices!, edges!, graph)
+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.
-    nd_ODE_Static_mod.StaticEdgeFunction(network_dynamics(vertices!, edges!, graph))
+    nd_ODE_Static_mod.StaticEdgeFunction(network_dynamics(vertices!, edges!, graph, parallel = parallel))
 end