PlantBiophysics.jl benchmark

The main objective of this notebook is to compare the computational times of PlantBiophysics.jl against the plantecophys R package and the LeafGasExchange.jl Julia package from the Cropbox.jl framework. The comparison follows three steps:

• create an N-large basis of random conditions.

• benchmark the computational time of the three packages via similar functions (i.e. photosynthesis-stomatal conductance-energy balance coupled model for C3 leaves): energy_balance, photosynEB and simulate with ModelC3MD.

• compare the results with plots and statistics.

This notebook does not perform the benchmark by itself for the obvious reason that it takes forever to run, and because there is an overhead cause by Pluto (benchmark and reactive are not a good mix). Instead, it shows the outputs of a script from this repository that implements the code shown here. If you want to perform the benchmark by yourself, you can run this script from the command line. The versions used for the above dependencies are available in the Project.toml of the repository.

Importing the dependencies:

Note

Make sure to have R installed on your computer first.

Loading the Julia packages:

begin
    using CairoMakie
	using BenchmarkTools
	using PlantBiophysics
    using Cropbox
    using LeafGasExchange
    using RCall
end

Parameters

Benchmark parameters

You'll find below the main parameters of the benchmark. In a few words, each package runs a simulation for N different time-steps microbenchmark_steps times repeated microbenchmark_evals times. We make N different simulations because the simulation duration can vary depending on the inputs due to iterative computations in the code, i.e. different initial conditions can make the algorithms converge more or less rapidly.

Random input simulation dataset

We create possible ranges for input parameters. These ranges where chosen so all of the three packages don't return errors during computation (plantecophys has issues with low temperatures).

• Ta: air temperature ($°C$)

• Wind: wind speed ($m.s^{-1}$)

• P: ambient pressure ($kPa$)

• Rh: relative humidity (between 0 and 1)

• Ca: air CO₂ concentration ($ppm$)

• Jmax: potential rate of electron transport ($\mu mol_{CO2}.m^{-2}.s^{-1}$)

• Vmax: maximum rate of Rubisco activity ($\mu mol_{CO2}.m^{-2}.s^{-1}$)

• Rd: mitochondrial respiration in the light at reference temperature ($\mu mol_{CO2}.m^{-2}.s^{-1}$)

• TPU: triose phosphate utilization-limited photosynthesis rate ($\mu mol_{CO2}.m^{-2}.s^{-1}$)

• Rs: short-wave net radiation ($W.m^{-1}$)

• skyF: Sun-visible fraction of the leaf (between 0 and 1)

• d: characteristic length ($m$)

• g0: residual stomatal conductance ($mol_{CO2}.m^{-2}.s^{-1}$)

• g1: slope of the stomatal conductance relationship.

We then sample N conditions from the given ranges:

Benchmarking

plantecophys

Preparing R to make the benchmark:

R"""
if(!require("plantecophys")){
    install.packages("plantecophys", repos = "https://cloud.r-project.org")
}
if(!require("microbenchmark")){
    install.packages("microbenchmark", repos = "https://cloud.r-project.org")
}
"""

# Make variables available to the R session
@rput set N microbenchmark_steps

Making the benchmark:

R"""
# Define the function call in a function that takes a list as input to limit DataFrame overhead
function_EB <- function(input) {
    PhotosynEB(
        Tair = input$Tair, VPD = input$VPD, Wind = input$Wind,
        Wleaf = input$Wleaf,Ca = input$Ca,  StomatalRatio = 1,
        LeafAbs = input$LeafAbs, gsmodel = "BBOpti", g0 = input$g0, g1 = input$g1,
        alpha = 0.24, theta = 0.7, Jmax = input$Jmax,
        Vcmax = input$Vcmax, TPU = input$TPU, Rd = input$Rd,
        RH = input$RH, PPFD=input$PPFD, Patm = input$Patm
    )
}

time_PE = c()
for(i in seq_len(N)){
    # Put the inputs into a vector to limit dataframe overhead:
    input = list(
        Tair = set$T[i], VPD = set$vpd[i], Wind = set$Wind[i], Wleaf = set$d[i],
        Ca = set$Ca[i], LeafAbs = set$sky_fraction[i], g0 = set$g0[i], g1 = set$g1[i],
        Jmax = set$JMaxRef[i], Vcmax = set$VcMaxRef[i], TPU = set$TPURef[i],
        Rd = set$RdRef[i], RH = set$Rh[i]*100, PPFD=set$PPFD[i],Patm = set$P[i]
    )

    m = microbenchmark(function_EB(input), times = microbenchmark_steps)

    time_PE = append(time_PE,m$time * 10e-9) # transform in seconds
}
"""

@rget time_PE

LeafGasExchange.jl

Note that we benchmark LeafGasExchange.jl with the nounit flag to compute a fair comparison with PlantBiophysics.jl in case computing units takes time (it shouldn't much).

time_LG = []
n_lg = fill(0, N)
for i = 1:N
    config =
        :Weather => (
            PFD=set.PPFD[i],
            CO2=set.Ca[i],
            RH=set.Rh[i] * 100,
            T_air=set.T[i],
            wind=set.Wind[i],
            P_air=set.P[i],
            g0=set.g0[i],
            g1=set.g1[i],
            Vcmax=set.VcMaxRef[i],
            Jmax=set.JMaxRef[i],
            Rd=set.RdRef[i],
            TPU=set.TPURef[i],
        )
    b_LG =
        @benchmark simulate($ModelC3MD; config=$config) evals = microbenchmark_evals samples =
            microbenchmark_steps
    append!(time_LG, b_LG.times .* 1e-9) # transform in seconds
    n_lg[i] = 1
end

PlantBiophysics.jl

Benchmarking PlantBiophysics.jl:

constants = Constants()
time_PB = []
for i = 1:N
    leaf = ModelList(
        energy_balance=Monteith(),
        photosynthesis=Fvcb(
            VcMaxRef=set.VcMaxRef[i],
            JMaxRef=set.JMaxRef[i],
            RdRef=set.RdRef[i],
            TPURef=set.TPURef[i],
        ),
        stomatal_conductance=Medlyn(set.g0[i], set.g1[i]),
        status=(
            Rₛ=set.Rs[i],
            sky_fraction=set.sky_fraction[i],
            PPFD=set.PPFD[i],
            d=set.d[i],
        ),
    )
    deps = PlantSimEngine.dep(leaf)
    meteo = Atmosphere(T=set.T[i], Wind=set.Wind[i], P=set.P[i], Rh=set.Rh[i], Cₐ=set.Ca[i])
    st = PlantMeteo.row_struct(leaf.status[1])
    b_PB = @benchmark run!($leaf, $deps, $st, $meteo, $constants, nothing) evals =
        microbenchmark_evals samples = microbenchmark_steps
    append!(time_PB, b_PB.times .* 1e-9) # transform in seconds
end

Comparison

Statistics

We compute here basic statistics, i.e. mean, median, min, max, standard deviation.

statsPB = basic_stat(time_PB)
statsPE = basic_stat(time_PE)
statsLG = basic_stat(time_LG)

factorPE = mean(time_PE) / mean(time_PB)
factorLG = mean(time_LG) / mean(time_PB)

# Write overall timings:
df = DataFrame(
	[getfield(j, i) for i in fieldnames(StatResults), j in [statsPB, statsPE, statsLG]],
	["PlantBiophysics", "plantecophys", "LeafGasExchange"]
)
insertcols!(df, 1, :Stat => [fieldnames(StatResults)...])
CSV.write("benchmark.csv", df)

# Write timing for each sample:
CSV.write("benchmark_full.csv",
	DataFrame(
		"package" => vcat(
			[
				repeat([i.first], length(i.second)) for i in [
					"PlantBiophysics" => time_PB,
					"plantecophys" => time_PE,
					"LeafGasExchange" => time_LG
				]
			]...
		),
		"sample_time" => vcat(time_PB, time_PE, time_LG)
	)
)

Histogram plotting

We here display the computational time histogram of each package on the same scale in order to compare them: PlantBiophysics.jl (a), plantecophys (b) and LeafGasExchange.jl (c). The y-axis represents the density (i.e. reaching 0.3 means that 30% of the computed times are in this bar). Orange zone represents the interval [mean - standard deviation; mean + standard deviation]. Red dashed line represents the mean. Note that the x-axis is logarithmic.

    fig = plot_benchmark_Makie(statsPB, statsPE, statsLG, time_PB, time_PE, time_LG)
    save("benchmark_each_time_steps.png", fig, px_per_unit=3)

References

Base.show

Base.show(io::IO, m::StatResults)
Base.show(io::IO, ::MIME"text/plain", m::StatResults)

Add a show method for our StatResults type.