GitHub

Package speed benchmarking

In this notebook, we will benchmark the speed of the package. Considering our framework can technically create any DAE-based model, we will limit this benchmark to the core functional modules of PlantModules.jl. Additionally, we consider two very straightforward plant structures for simplicity's sake.

Setup

using Pkg; Pkg.activate("../..")
using PlutoUI; TableOfContents()
using PlantModules
using ModelingToolkit, OrdinaryDiffEq, Plots
using PlantGraphs

System definition

Structure

For plant structure, we will consider plants of a single node type Segment connected to an environment consisting of a single soil and air node. In order to inspect the influence of the level of branching on solution time, we consider two types of structure:

  • A completely linear plant, where the base is connected to the soil and all other nodes are connected to the air. We can imagine this structure to represent something like a conifer, with each plant segment containing a part of the stem and all connected branches.

  • A plant of which each segment branches into two other segments, where the base is connected to the soil and all leaf nodes are connected to the air. We can imagine this to represent a herbaceous plant of some kind.

We will define the structure of the plants using rewrite rules so we can easily change the size of the system, using PlantGraphs.jl from the VirtualPlantLab.jl ecosystem.

Base.@kwdef struct Segment <: PlantGraphs.Node
    D::Vector{Float64} = PlantModules.default_values[:D]
end
Segment
struct Soil <: PlantGraphs.Node end
struct Air <: PlantGraphs.Node end

Linear

# all segments split into two linearly connected segments
linear_rule = Rule(
    Segment,
    rhs = seg -> Segment() + Segment()
)
Rule replacing nodes of type Segment without context capturing.

# get the structure of the linear plant (without and with connections to environment) for a given number of rewrite steps
function get_structure_linear(rewrite_steps)
    plant = Graph(axiom = Segment(), rules = linear_rule)
    for i in 1:rewrite_steps
        rewrite!(plant)
    end

    plantstructure = PlantStructure(
        [plant, Soil(), Air()],
        [
            # connect first segment to soil
            (1, 2) => (getnodes(plant)[1], :Soil), 
            # connect all others to the air
            (1, 3) => (seg, air) -> seg != getnodes(plant)[1] 
        ]
    )
    
    return plant, plantstructure
end
get_structure_linear (generic function with 1 method)

Branching

# all leaf nodes branch out into two more segments, where the cross area of the parent segment equals the summed cross areas of the child segments (as per Da Vinci's rule)
branching_rule = Rule(
    Segment,
    lhs = seg -> !has_children(seg),
    rhs = seg -> Segment(data(seg).D) + 
        (
            Segment(data(seg).D .* [1/sqrt(2), 1]), 
            Segment(data(seg).D .* [1/sqrt(2), 1])
        )
)
Rule replacing nodes of type Segment without context capturing.

# get the structure of the branching plant (without and with connections to environment) for a given number of rewrite steps
function get_structure_branching(rewrite_steps)
    plant = Graph(axiom = Segment(), rules = branching_rule)
    for i in 1:rewrite_steps
        rewrite!(plant)
    end

    plantstructure = PlantStructure(
        [plant, Soil(), Air()],
        [
            (1, 2) => (getnodes(plant)[1], :Soil),
            (1, 3) => (seg, air) -> is_leaf(seg)
        ]
    )
    
    return plant, plantstructure
end
get_structure_branching (generic function with 1 method)

Coupling

All node types are assigned default hydraulic functionality.

module_coupling = Dict(
    :Segment => [hydraulic_module, constant_carbon_module, K_module],
    :Soil => [environmental_module, Ψ_soil_module, constant_K_module],
    :Air => [environmental_module, Ψ_air_module, constant_K_module],
);
connecting_modules = Dict(
    (:Soil, :Segment) => constant_hydraulic_connection,
    (:Segment, :Segment) => hydraulic_connection,
    (:Segment, :Air) => daynight_hydraulic_connection
);
const plantcoupling = PlantCoupling(; module_coupling, connecting_modules); # mark global variables as having a constant type for increased efficiency when using them in functions

Parameters

We change the following parameter values from their default values:

  • The specific hydraulic conductivity K_s of the plant segments is increased to prevent larger systems from having unrealistically small hydraulic conductivities, considering the cross area of their segments is not increased in our model.

  • The maximum water content W_max of the soil is reduced to match our relatively small plant structures.

  • The hydraulic conductivity K of the air is set to a smaller, more realistic value, as the default value is intended to represent a well-conducting node.

  • The relative water content W_r of the air is set to a smaller value for similar reasons.

module_defaults = Dict(
    :Segment => Dict(:K_s => 500.0),
    :Soil => Dict(:W_max => 1e4),
    :Air => Dict(:K => 1e-2, :W_r => 0.7)
    )
Dict{Symbol, Dict{Symbol, Float64}} with 3 entries:
  :Soil    => Dict(:W_max=>10000.0)
  :Air     => Dict(:W_r=>0.7, :K=>0.01)
  :Segment => Dict(:K_s=>500.0)
plantparams = PlantParameters(; module_defaults);

Sanity check

Before getting to the actual benchmarks, we'll do a quick sanity check to verify our system definition returns sufficiently realistic simulations.

function test_simulation(plantstructure; tspan = (0.0, 7*24.0))
    system = generate_system(plantstructure, plantcoupling, plantparams)
    prob = ODEProblem(system, [], tspan, sparse = true)
    sol = solve(prob, FBDF())

    return sol
end
test_simulation (generic function with 1 method)
rewrite_steps_test = 5
5

Linear

test_plant_linear, test_structure_linear = get_structure_linear(rewrite_steps_test)
(Dynamic graph with 32 nodes of types Main.var"workspace#33".Segment and 1 rewriting rules.
, PlantStructure{Int64}([1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  25, 26, 27, 28, 29, 30, 31, 32, 33, 34], PlantModules.PMEdge{Int64}[PlantModules.PMEdge{Int64}(1, 2), PlantModules.PMEdge{Int64}(1, 33), PlantModules.PMEdge{Int64}(2, 3), PlantModules.PMEdge{Int64}(2, 1), PlantModules.PMEdge{Int64}(2, 34), PlantModules.PMEdge{Int64}(3, 4), PlantModules.PMEdge{Int64}(3, 2), PlantModules.PMEdge{Int64}(3, 34), PlantModules.PMEdge{Int64}(4, 5), PlantModules.PMEdge{Int64}(4, 3)  …  PlantModules.PMEdge{Int64}(34, 23), PlantModules.PMEdge{Int64}(34, 24), PlantModules.PMEdge{Int64}(34, 25), PlantModules.PMEdge{Int64}(34, 26), PlantModules.PMEdge{Int64}(34, 27), PlantModules.PMEdge{Int64}(34, 28), PlantModules.PMEdge{Int64}(34, 29), PlantModules.PMEdge{Int64}(34, 30), PlantModules.PMEdge{Int64}(34, 31), PlantModules.PMEdge{Int64}(34, 32)], Dict{Int64, PlantModules.PMVertex{Int64}}(5 => PlantModules.PMVertex{Int64}(5, :Segment, Dict(:D => [0.5, 5.0])), 16 => PlantModules.PMVertex{Int64}(16, :Segment, Dict(:D => [0.5, 5.0])), 20 => PlantModules.PMVertex{Int64}(20, :Segment, Dict(:D => [0.5, 5.0])), 12 => PlantModules.PMVertex{Int64}(12, :Segment, Dict(:D => [0.5, 5.0])), 24 => PlantModules.PMVertex{Int64}(24, :Segment, Dict(:D => [0.5, 5.0])), 28 => PlantModules.PMVertex{Int64}(28, :Segment, Dict(:D => [0.5, 5.0])), 8 => PlantModules.PMVertex{Int64}(8, :Segment, Dict(:D => [0.5, 5.0])), 17 => PlantModules.PMVertex{Int64}(17, :Segment, Dict(:D => [0.5, 5.0])), 30 => PlantModules.PMVertex{Int64}(30, :Segment, Dict(:D => [0.5, 5.0])), 1 => PlantModules.PMVertex{Int64}(1, :Segment, Dict(:D => [0.5, 5.0]))…), Dict(5 => [4, 6, 34], 16 => [15, 17, 34], 20 => [19, 21, 34], 12 => [11, 13, 34], 24 => [23, 25, 34], 28 => [27, 29, 34], 8 => [7, 9, 34], 17 => [16, 18, 34], 30 => [29, 31, 34], 1 => [2, 33]…), Dict(5 => (1, 437), 16 => (1, 448), 20 => (1, 452), 12 => (1, 444), 24 => (1, 456), 28 => (1, 460), 8 => (1, 440), 17 => (1, 449), 30 => (1, 462), 1 => (1, 433)…)))
test_sol_linear = test_simulation(test_structure_linear);

Branching

test_plant_branching, test_structure_branching = get_structure_branching(rewrite_steps_test)
(Dynamic graph with 63 nodes of types Main.var"workspace#33".Segment and 1 rewriting rules.
, PlantStructure{Int64}([1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  56, 57, 58, 59, 60, 61, 62, 63, 64, 65], PlantModules.PMEdge{Int64}[PlantModules.PMEdge{Int64}(1, 2), PlantModules.PMEdge{Int64}(1, 3), PlantModules.PMEdge{Int64}(1, 64), PlantModules.PMEdge{Int64}(2, 4), PlantModules.PMEdge{Int64}(2, 5), PlantModules.PMEdge{Int64}(2, 1), PlantModules.PMEdge{Int64}(3, 6), PlantModules.PMEdge{Int64}(3, 7), PlantModules.PMEdge{Int64}(3, 1), PlantModules.PMEdge{Int64}(4, 8)  …  PlantModules.PMEdge{Int64}(65, 50), PlantModules.PMEdge{Int64}(65, 51), PlantModules.PMEdge{Int64}(65, 53), PlantModules.PMEdge{Int64}(65, 54), PlantModules.PMEdge{Int64}(65, 56), PlantModules.PMEdge{Int64}(65, 57), PlantModules.PMEdge{Int64}(65, 59), PlantModules.PMEdge{Int64}(65, 60), PlantModules.PMEdge{Int64}(65, 62), PlantModules.PMEdge{Int64}(65, 63)], Dict{Int64, PlantModules.PMVertex{Int64}}(5 => PlantModules.PMVertex{Int64}(5, :Segment, Dict(:D => [0.24999999999999994, 5.0])), 56 => PlantModules.PMVertex{Int64}(56, :Segment, Dict(:D => [0.0883883476483184, 5.0])), 35 => PlantModules.PMVertex{Int64}(35, :Segment, Dict(:D => [0.0883883476483184, 5.0])), 55 => PlantModules.PMVertex{Int64}(55, :Segment, Dict(:D => [0.12499999999999996, 5.0])), 60 => PlantModules.PMVertex{Int64}(60, :Segment, Dict(:D => [0.0883883476483184, 5.0])), 30 => PlantModules.PMVertex{Int64}(30, :Segment, Dict(:D => [0.0883883476483184, 5.0])), 32 => PlantModules.PMVertex{Int64}(32, :Segment, Dict(:D => [0.0883883476483184, 5.0])), 6 => PlantModules.PMVertex{Int64}(6, :Segment, Dict(:D => [0.24999999999999994, 5.0])), 45 => PlantModules.PMVertex{Int64}(45, :Segment, Dict(:D => [0.0883883476483184, 5.0])), 64 => PlantModules.PMVertex{Int64}(64, :Soil, Dict{Any, Any}())…), Dict(5 => [2, 10, 11], 56 => [55, 65], 55 => [14, 56, 57], 35 => [34, 65], 60 => [58, 65], 30 => [28, 65], 32 => [31, 65], 6 => [3, 12, 13], 45 => [43, 65], 64 => [1]…), Dict(5 => (1, 478), 56 => (1, 551), 16 => (1, 511), 20 => (1, 515), 35 => (1, 530), 55 => (1, 550), 60 => (1, 555), 30 => (1, 525), 19 => (1, 514), 32 => (1, 527)…)))
test_sol_branching = test_simulation(test_structure_branching);
begin
    plotgraph(test_sol_branching, test_structure_branching, 
              varname = :P, structmod = :Segment, 
              label = "Segment turgor pressure", 
              title = "Turgor pressure over time",
              xlabel = "Time (h)", ylabel = "P (MPa)",
              size = (800, 400), margins = 5*Plots.mm,
              linewidth = 2
    )
    hline!(
        [PlantModules.default_values[:Γ]], label = "Yield turgor pressure", lw = 2
    )
end

Benchmarking

Now let's get to the actual benchmarking. We will time the system generation, problem generation and problem solving for differing numbers of rewrite steps to evaluate how computation time scales with system size.

function get_stats(plantstructure; tspan = (0.0, 7*24.0))
    system_stats = @timed generate_system(
        plantstructure, plantcoupling, plantparams
    )
    prob_stats = @timed ODEProblem(
        system_stats.value, [], tspan, sparse = true
    )
    sol_stats = @timed solve(prob_stats.value, FBDF())

    num_nodes = length(getnodes(plantstructure))
    num_variables = length(unknowns(system_stats.value))

    return num_nodes, num_variables, system_stats, prob_stats, sol_stats
end
get_stats (generic function with 1 method)
Warning

Note that all benchmarks are based on only a single run, which makes them somewhat inconsistent compared to taking the median of a sample of runs. The reason for this is that we also benchmark function compilation time, which only triggers the first time a function method is used, making it non-trivial to get a sample of.

benchmark_linear (generic function with 1 method)
benchmark_branching (generic function with 1 method)
stackedbar (generic function with 1 method)
plot_size (generic function with 1 method)
plot_time (generic function with 1 method)
rewrite_steps_set = [2, 4, 6] # note not to include `rewrite_steps_test`, as the functions have already compiled for this system size
3-element Vector{Int64}:
 2
 4
 6

Linear structure

stats_linear = [benchmark_linear(rewrite_steps) for rewrite_steps in rewrite_steps_set];

Branching structure

stats_branching = [benchmark_branching(rewrite_steps) for rewrite_steps in rewrite_steps_set];
p_time_branching = plot_time(rewrite_steps_set, stats_branching)