Merge pull request #220 from JuliaDynamics/hw/update_benchmarking

update benchmark script

Author: hexaeder

benchmark/benchmark_cases.jl

@@ -0,0 +1,105 @@
+"""
+    BenchmarkCase
+
+Represents a single benchmark case for NetworkDynamics.
+
+Fields:
+- name: String identifier for the benchmark
+- constructor: Function that creates the Network for a given N
+- Ns: Vector of network sizes to benchmark
+- description: Documentation string
+"""
+struct BenchmarkCase
+    name::String
+    constructor::Function
+    Ns::Vector{Int}
+    description::String
+end
+
+# Standard benchmark cases
+BENCHMARK_CASES = [
+    BenchmarkCase(
+        "diffusion_static_edge",
+        N -> begin
+            g = watts_strogatz(N, N ÷ 2, 0.0; rng=StableRNG(1))
+            (g, diffusion_vertex(), diffusion_edge())
+        end,
+        # [10,20],
+        [100, 300, 1000, 3000],
+        "Diffusion network with static edges"
+    ),
+
+    BenchmarkCase(
+        "diffusion_ode_edge",
+        N -> begin
+            g = watts_strogatz(N, N ÷ 2, 0.0; rng=StableRNG(1))
+            (g, diffusion_vertex(), diffusion_dedge())
+        end,
+        # [10,20],
+        [100, 300, 1000, 3000],
+        "Diffusion network with ODE edges"
+    ),
+
+    BenchmarkCase(
+        "kuramoto_homogeneous",
+        N -> begin
+            g = watts_strogatz(N, 3, 0.8; rng=StableRNG(1))
+            (g, kuramoto_vertex_2d(), static_kuramoto_edge())
+        end,
+        # [10,20],
+        [100, 1_000, 10_000, 100_000],
+        "Homogeneous Kuramoto oscillators"
+    ),
+
+    BenchmarkCase(
+        "kuramoto_heterogeneous",
+        N -> begin
+            g = watts_strogatz(N, 3, 0.8; rng=StableRNG(1))
+            rng = StableRNG(1)
+            vtypes = [kuramoto_vertex_1d(), kuramoto_vertex_2d()]
+            vertices = vtypes[shuffle(rng, vcat([1 for _ in 1:N÷2], [2 for _ in 1:N÷2]))]
+            (g, vertices, static_kuramoto_edge())
+        end,
+        # [10,20],
+        [100, 1_000, 10_000, 100_000],
+        "Heterogeneous Kuramoto oscillators"
+    ),
+
+    BenchmarkCase(
+        "powergrid",
+        N -> begin
+            g = watts_strogatz(N, 3, 0.8; rng=StableRNG(1))
+            rng = StableRNG(1)
+            vtypes = [pqnode(), generator()]
+            vertices = vtypes[shuffle(rng, vcat([1 for _ in 1:N÷2], [2 for _ in 1:N÷2]))]
+            (g, vertices, piline())
+        end,
+        # [10,20],
+        [100, 1_000, 10_000, 100_000],
+        "Power grid network with PQ nodes and generators"
+    ),
+
+    BenchmarkCase(
+        "powergrid_inhomogeneous_pq",
+        N -> begin
+            g = watts_strogatz(N, 3, 0.8; rng=StableRNG(1))
+            rng = StableRNG(1)
+
+            # Create fixed P/Q values for half the nodes
+            pq_count = N ÷ 2
+            P_vals = -1*rand(rng, pq_count)
+            Q_vals = -1*rand(rng, pq_count)
+
+            # Create vertex models with baked-in P/Q values
+            pq_vertices = [pqnode_inhomogeneous(P_vals[i], Q_vals[i]) for i in 1:pq_count]
+            gen_vertices = [generator() for _ in 1:(N-pq_count)]
+
+            # Shuffle the vertex types
+            vertices = shuffle(rng, vcat(pq_vertices, gen_vertices))
+
+            (g, vertices, piline())
+        end,
+        [100, 1_000, 10_000, 100_000],
+        "Power grid with heterogeneous PQ nodes having fixed P/Q values"
+    )
+]

benchmark/benchmark_compat.jl

@@ -1,88 +1,14 @@
 #=
-Helper functions to create compatible states between NetworkDynamics.jl and NetworkDynamics.jl v0.8
+Helper functions to create compatible states between newer and older versions of NetworkDynamics.jl
+for cross version benchmarking.
 =#
 
-if pkgversion(NetworkDynamics) < v"0.9.0"
-    @info "Define compatibility functions"
-    struct Network end
-    struct SequentialExecution{P} end
-    struct KAExecution{P} end
-    struct PolyesterExecution{P} end
-    struct ThreadedExecution{P} end
-    struct KAAggregator; f; end
-    struct SequentialAggregator; f; end
-    struct PolyesterAggregator; f; end
-    struct ThreadedAggregator; f; end
-    function Network(g, v, e; execution = SequentialExecution{true}(), aggregator=SequentialAggregator(+))
-        if execution isa PolyesterExecution{true} && aggregator isa PolyesterAggregator
-            network_dynamics(v, e, g; parallel=true)
-        elseif execution isa SequentialExecution{true} && aggregator isa SequentialAggregator
-            network_dynamics(v, e, g; parallel=false)
-        else
-            error("execution type not supported")
-        end
-    end
-elseif pkgversion(NetworkDynamics) < v"0.9.1"
-    isdynamic = NetworkDynamics.isdynamic
-else
-    isdynamic(_) = true
-end
-
-function _syms_old_order(nd::Network)
-    syms = []
-    for (i,cf) in enumerate(nd.im.vertexm)
-        isdynamic(cf) || continue
-        append!(syms, collect(VIndex(i, 1:dim(cf))))
-    end
-    for (i,cf) in enumerate(nd.im.edgem)
-        isdynamic(cf) || continue
-        append!(syms, collect(EIndex(i, 1:dim(cf))))
-    end
-    syms
-end
-function _psyms_old_order(nd::Network)
-    syms = []
-    for (i,cf) in enumerate(nd.im.vertexm)
-        append!(syms, collect(VPIndex(i, 1:pdim(cf))))
-    end
-    for (i,cf) in enumerate(nd.im.edgem)
-        append!(syms, collect(EPIndex(i, 1:pdim(cf))))
-    end
-    syms
-end
 function randx0(nd::Network)
     rng = StableRNG(1)
-    _x0 = rand(rng, dim(nd))
-
-    s = NWState(nd; utype=typeof(_x0))
-    for (i,sym) in enumerate(_syms_old_order(nd::Network))
-        s[sym] = _x0[i]
-    end
-    s[:]
-end
-randx0(nd::ODEFunction) = rand(StableRNG(1), length(nd.syms))
-function legacy_order(nd::Network, _dx)
-    s = NWState(nd, _dx)
-    dx = similar(_dx)
-    for (i,sym) in enumerate(_syms_old_order(nd::Network))
-        dx[i] = s[sym]
-    end
-    dx
+    rand(rng, dim(nd))
 end
-legacy_order(nd::ODEFunction, dx) = dx
 
 function randp(nd::Network)
     rng = StableRNG(1)
-    _p = rand(rng, pdim(nd))
-    p = NWParameter(nd; ptype=typeof(_p))
-    for (i, sym) in enumerate(_psyms_old_order(nd))
-        p[sym] = _p[i]
-    end
-    p[:]
-end
-function randp(nd::ODEFunction)
-    rng = StableRNG(1)
-    g = nd.f.graph
-    _p = rand(rng, nv(g) + ne(g))
-    p = (_p[1:nv(g)], _p[nv(g)+1:end])
+    rand(rng, pdim(nd))
 end

benchmark/benchmark_models.jl

@@ -109,3 +109,15 @@ Base.@propagate_inbounds function gen!(dv, v, isum, p, T)
 end
 generator() = VertexModel(; f=gen!, sym=[:u_r, :u_i, :θ, :ω], g=StateMask(1:2),
             psym=[:H, :P, :D, :Ω, :E_f, :T_d_dash, :T_q_dash, :X_d_dash, :X_q_dash, :X_d, :X_q])
+
+function pqnode_inhomogeneous(P, Q)
+    Base.@propagate_inbounds function pq!(du, u, isum, _, t)
+        ic = Complex(isum[1], isum[2])
+        S = Complex(P, Q)
+        uc = S/conj(ic)
+        du[1] = real(uc)
+        du[2] = imag(uc)
+        nothing
+    end
+    VertexModel(; f=pq!, sym=[:u_r, :u_i], g=StateMask(1:2), psym=[], mass_matrix=0)
+end

benchmark/benchmark_models_v0.8.jl

@@ -244,10 +244,11 @@ function plot_over_N(target, baseline=nothing)
         comp = target
     end
 
-    fig = Makie.Figure(size=(2000,2000))
+    fig = Makie.Figure(size=(2000,2500))
 
     bmkeys = [
         "powergrid",
+        "powergrid_inhomogeneous_pq",
         "diffusion_static_edge",
         "diffusion_ode_edge",
         "kuramoto_homogeneous",