diff --git a/.gitignore b/.gitignore index 82ac67d5..e2f1c9fd 100644 --- a/.gitignore +++ b/.gitignore @@ -7,3 +7,4 @@ Manifest.toml scratch.jl docs/build bugs +design_experiments diff --git a/Project.toml b/Project.toml index c16775fd..5bfa8db2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "NetworkDynamics" uuid = "22e9dc34-2a0d-11e9-0de0-8588d035468b" -authors = ["Frank Hellmann , Alexander Penner "] -version = "0.2.0" +authors = ["Frank Hellmann , Michael Lindner , Alexander Penner "] +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" diff --git a/benchmarking/Project.toml b/benchmarking/Project.toml index c6f11dfa..7aeb6ce6 100644 --- a/benchmarking/Project.toml +++ b/benchmarking/Project.toml @@ -1,9 +1,11 @@ name = "NetworkDynamics-Benchmarking" -author = ["Frank Hellmann "] +authors = ["Frank Hellmann ", + "Michael Lindner "] [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" diff --git a/benchmarking/multithreading_benchmarks.jl b/benchmarking/multithreading_benchmarks.jl new file mode 100644 index 00000000..eb4b10e6 --- /dev/null +++ b/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")) + +=# diff --git a/docs/localmake.jl b/docs/localmake.jl index 35895fa0..ccc7c6ba 100644 --- a/docs/localmake.jl +++ b/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", diff --git a/docs/make.jl b/docs/make.jl index b77c18c1..8b1cd948 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -13,6 +13,7 @@ makedocs( pages = [ "General" => "index.md", "BasicConstructors.md", + "Multithreading.md", "Library.md", "Tutorials" => [ "Getting started" => "getting_started_with_network_dynamics.md", @@ -20,15 +21,8 @@ makedocs( ] ]) -# 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 ) diff --git a/docs/src/BasicConstructors.md b/docs/src/BasicConstructors.md index bce430af..7d988413 100644 --- a/docs/src/BasicConstructors.md +++ b/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: diff --git a/docs/src/Multithreading.md b/docs/src/Multithreading.md new file mode 100644 index 00000000..18efc3e1 --- /dev/null +++ b/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). diff --git a/docs/src/directed_and_weighted_graphs.md b/docs/src/directed_and_weighted_graphs.md index b719eb34..e61d4f15 100644 --- a/docs/src/directed_and_weighted_graphs.md +++ b/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 diff --git a/docs/src/getting_started_with_network_dynamics.md b/docs/src/getting_started_with_network_dynamics.md index b8186b89..7e4640eb 100644 --- a/docs/src/getting_started_with_network_dynamics.md +++ b/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. diff --git a/docs/src/index.md b/docs/src/index.md index 3e04786f..8112a145 100644 --- a/docs/src/index.md +++ b/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. :) diff --git a/src/NetworkDynamics.jl b/src/NetworkDynamics.jl index 2666d1e9..d3f54cfb 100644 --- a/src/NetworkDynamics.jl +++ b/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,9 +176,9 @@ 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 @@ -159,10 +186,10 @@ 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 diff --git a/src/Utilities.jl b/src/Utilities.jl index a0d6e07a..a5750283 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -3,7 +3,51 @@ module Utilities using NLsolve using LinearAlgebra -export maybe_idx, p_v_idx, p_e_idx +export maybe_idx, p_v_idx, p_e_idx, @nd_threads + + +""" +nd_threads: Allows control over threading by the 1st argument, +a boolean (likely from the network object). +This allows you to control for multithreading at runtime without +code duplication. +""" + +macro nd_threads(trigger,args...) + na = length(args) + if na != 1 + throw(ArgumentError("wrong number of arguments in @nd_threads")) + end + ex = args[1] + if ex.head == :for + # Need to escape the bounds of the loop, but not the whole expression, + # in order to work with the base Threads.@thread macro. + # That macro does it's own escaping, but also it's own expression + # checking (the head of the first expression must be :for. So we will + # escape the bounds here, keeping the expression formed in a suitable + # way for @threads + # + # This is a brute-force drill down an AST of a for loop to the limits + # of said loop. See: dump(quote for i in e.a:e.b print(i) end end) + ex.args[1].args[1] = esc(ex.args[1].args[1]) + ex.args[1].args[2].args[2] = esc(ex.args[1].args[2].args[2]) + ex.args[1].args[2].args[3] = esc(ex.args[1].args[2].args[3]) + + # Same for the body of the loop + ex.args[2] = esc(ex.args[2]) + + return quote + @inbounds if $(esc(trigger)) + Threads.@threads $ex + else + $(ex) + end + end + else + throw(ArgumentError("unrecognized argument to @nd_threads")) + end +end + """ Utility function that drops the indexing operation when the argument is not diff --git a/src/nd_ODE_ODE.jl b/src/nd_ODE_ODE.jl index b48adae3..455243d8 100644 --- a/src/nd_ODE_ODE.jl +++ b/src/nd_ODE_ODE.jl @@ -38,23 +38,21 @@ end graph::G graph_structure::GraphStruct graph_data::GraphData{T, T} + parallel::Bool end function (d::nd_ODE_ODE)(dx, x, p, t) gd = prep_gd(view(dx, 1:2), view(x, 1:2), x, d.graph_data, d.graph_structure) - @inbounds begin - - for i in 1:d.graph_structure.num_e - maybe_idx(d.edges!, i).f!(view(dx,d.graph_structure.e_idx[i] .+ d.graph_structure.dim_v), gd.e[i], gd.v_s_e[i], gd.v_d_e[i], p_e_idx(p, i), t) - end - - for i in 1:d.graph_structure.num_v - maybe_idx(d.vertices!,i).f!(view(dx,d.graph_structure.v_idx[i]), gd.v[i], gd.e_s_v[i], gd.e_d_v[i], p_v_idx(p, i), t) + @nd_threads d.parallel for i in 1:d.graph_structure.num_e + maybe_idx(d.edges!, i).f!(view(dx,d.graph_structure.e_idx[i] .+ d.graph_structure.dim_v), gd.e[i], gd.v_s_e[i], gd.v_d_e[i], p_e_idx(p, i), t) end + @nd_threads d.parallel for i in 1:d.graph_structure.num_v + maybe_idx(d.vertices!,i).f!(view(dx,d.graph_structure.v_idx[i]), gd.v[i], gd.e_s_v[i], gd.e_d_v[i], p_v_idx(p, i), t) end + nothing end diff --git a/src/nd_ODE_Static.jl b/src/nd_ODE_Static.jl index 8550444a..9e9beb81 100644 --- a/src/nd_ODE_Static.jl +++ b/src/nd_ODE_Static.jl @@ -1,6 +1,5 @@ module nd_ODE_Static_mod - using ..NetworkStructures using ..NDFunctions using ..Utilities @@ -22,32 +21,28 @@ end end + @Base.kwdef struct nd_ODE_Static{G, Tv, Te, T1, T2} vertices!::T1 edges!::T2 graph::G graph_structure::GraphStruct graph_data::GraphData{Tv, Te} + parallel::Bool # enables multithreading for the core loop end function (d::nd_ODE_Static)(dx, x, p, t) gd = prep_gd(dx, x, d.graph_data, d.graph_structure) - @inbounds begin - - for i in 1:d.graph_structure.num_e - # maybe_idx(d.edges!,i).f!(gd.e[i], gd.v_s_e[i], gd.v_d_e[i], maybe_idx(p, i+d.graph_structure.num_v), t) + @nd_threads d.parallel for i in 1:d.graph_structure.num_e maybe_idx(d.edges!, i).f!(gd.e[i], gd.v_s_e[i], gd.v_d_e[i], p_e_idx(p, i), t) end - for i in 1:d.graph_structure.num_v - # maybe_idx(d.vertices!,i).f!(view(dx,d.graph_structure.v_idx[i]), gd.v[i], gd.e_s_v[i], gd.e_d_v[i], maybe_idx(p, i), t) + @nd_threads d.parallel for i in 1:d.graph_structure.num_v maybe_idx(d.vertices!,i).f!(view(dx,d.graph_structure.v_idx[i]), gd.v[i], gd.e_s_v[i], gd.e_d_v[i], p_v_idx(p, i), t) end - end - nothing end @@ -55,14 +50,11 @@ end function (d::nd_ODE_Static)(x, p, t, ::Type{GetGD}) gd = prep_gd(x, x, d.graph_data, d.graph_structure) - @inbounds begin - for i in 1:d.graph_structure.num_e + @nd_threads d.parallel for i in 1:d.graph_structure.num_e maybe_idx(d.edges!,i).f!(gd.e[i], gd.v_s_e[i], gd.v_d_e[i], p_e_idx(p, i), t) end - end - gd end @@ -71,6 +63,7 @@ function (d::nd_ODE_Static)(::Type{GetGS}) d.graph_structure end + # For compatibility with PowerDynamics struct StaticEdgeFunction diff --git a/test/diffusion_test.jl b/test/diffusion_test.jl index 880fb39d..474e06b4 100644 --- a/test/diffusion_test.jl +++ b/test/diffusion_test.jl @@ -142,11 +142,3 @@ println("Maximum difference to analytic solution with fake ode edge ND: $max_ode syms_v = syms_containing(diff_network_ode, "v") idx_v = idx_containing(diff_network_ode, "v") @test all( diff_network_ode.syms[idx_v] .== syms_v ) - -sef = StaticEdgeFunction(diff_network_st.f) -t = 1. -e_s, e_d = sef(sol_st(t), nothing, t) -e_int_ode = sol_ode(t)[N+1:N+ne(g)] -e_int_st = diff_network_st.f.graph_data.e_array - -println("Maximum difference of the static and dynamic edge variables at t=$t: $(abs.(e_int_st .- e_int_ode) |> maximum)")