Update benchmark_w_matlab.jl

examples/benchmark_w_matlab.jl

@@ -1,11 +1,10 @@
-
 ## Comparing the performance of the Julia and MATLAB implementations
 
 # We can compare the performance of the Julia and MATLAB implementations
 # by running the same model for the same number of epochs and measuring
 # the time taken.
 
-using BeforeIT, StatsPlots, Statistics, ThreadPinning
+using BeforeIT, CairoMakie, Statistics, ThreadPinning
 
 pinthreads(:cores)
 
@@ -29,78 +28,98 @@ T = 12
 @time run(parameters, initial_conditions, T; multi_threading = true);
 
 # Time taken by the MATLAB code and the Generated C code with MATLAB Coder
-# (6 threads for the parallel version), computed independently on an AMD Ryzen 5 5600H
-matlab_times = [4.399592, 4.398576, 4.352314, 4.385039, 4.389989]
-matlab_times ./= T
-matlab_times = [matlab_times, [2000.189556, 1990.485305, 2002.295329, 1994.215843, 1993.651854]]
-matlab_time = mean.(matlab_times)
-matlab_time_std = std.(matlab_times)
-
-c_times = [0.952, 0.940, 0.951, 0.942, 0.938]
-c_times ./= T
-c_times = [c_times, [574.692, 576.725, 572.382, 573.921, 575.949]]
-c_time = mean.(c_times)
-c_time_std = std.(c_times)
-
-c_times_multi_thread = [0.305, 0.324, 0.330, 0.334, 0.323]
-c_times_multi_thread ./= T
-c_times_multi_thread = [c_times_multi_thread, [162.603, 160.197, 161.610, 160.759, 161.174]]
-c_time_multi_thread = mean.(c_times_multi_thread)
-c_time_multi_thread_std = std.(c_times_multi_thread)
-
-# Time taken by the Julia code (same platform as in the the matlab benchmarks),
+# (4 threads for the parallel version), computed independently on an 
+# AMD Ryzen 5 5600H
+matlab_times_small = [4.399592, 4.398576, 4.352314, 4.385039, 4.389989] ./ T
+matlab_times_big = [2000.189556, 1990.485305, 2002.295329, 1994.215843, 1993.651854]
+matlab_mtime_small = mean(matlab_times_small)
+matlab_stime_small = std(matlab_times_small)
+matlab_mtime_big = mean(matlab_times_big)
+matlab_stime_big = std(matlab_times_big)
+
+c_times_small = [0.952, 0.940, 0.951, 0.942, 0.938] ./ T
+c_times_big = [574.692, 576.725, 572.382, 573.921, 575.949]
+c_mtime_small = mean(c_times_small)
+c_stime_small = std(c_times_small)
+c_mtime_big = mean(c_times_big)
+c_stime_big = std(c_times_big)
+
+c_times_small_multi = [0.305, 0.324, 0.330, 0.334, 0.323] ./ T
+c_times_big_multi = [209.603, 210.197, 209.610, 211.759, 209.174]
+c_mtime_small_multi = mean(c_times_small_multi)
+c_stime_small_multi = std(c_times_small_multi)
+c_mtime_big_multi = mean(c_times_big_multi)
+c_stime_big_multi = std(c_times_big_multi)
+
+# timings from "High-performance computing implementations of agent-based
+# economic models for realizing 1: 1 scale simulations of large
+# economies", A. Gill et al. (2021)
+hpc_mtime_big_multi = 22.92
+hpc_mtime_big = hpc_mtime_big_multi*4
+
+# Time taken by the Julia code (same platform as in the matlab benchmarks),
 # computed as the average of 5 runs
 n_runs = 5
 
-julia_times_1_thread = zeros(n_runs)   
+julia_times_small = zeros(n_runs)   
 for i in 1:n_runs
-    julia_times_1_thread[i] = @elapsed run(parameters, initial_conditions, T; multi_threading = false);
+    julia_times_small[i] = @elapsed run(parameters, initial_conditions, T; multi_threading = false);
 end
-julia_times_1_thread ./= T
-julia_times_1_thread = [julia_times_1_thread, [56.593792, 56.297404, 54.682004, 55.079674, 55.192496]]
-julia_time_1_thread = mean.(julia_times_1_thread)
-julia_time_1_thread_std = std.(julia_times_1_thread)
-
-julia_times_multi_thread = zeros(n_runs)
+julia_times_small ./= T
+julia_times_big = [56.593792, 56.297404, 54.682004, 55.079674, 55.192496]
+julia_mtime_small = mean(julia_times_small)
+julia_stime_small = std(julia_times_small)
+julia_mtime_big = mean(julia_times_big)
+julia_stime_big = std(julia_times_big)
+
+julia_times_small_multi = zeros(n_runs)
 for i in 1:5
-    julia_times_multi_thread[i] =  @elapsed run(parameters, initial_conditions, T; multi_threading = true);
+    julia_times_small_multi[i] =  @elapsed run(parameters, initial_conditions, T; multi_threading = true);
 end
-julia_times_multi_thread ./= T
-julia_times_multi_thread = [julia_times_multi_thread, [30.775353, 29.533283, 30.594063, 30.469412, 30.521846]]
-julia_time_multi_thread = mean.(julia_times_multi_thread)
-julia_time_multi_thread_std = std.(julia_times_multi_thread)
-
-# Get the number of threads used
-n_threads = Threads.nthreads()
-
-theme(:default, bg = :white)
-
-means = reduce(hcat, [matlab_time, c_time, c_time_multi_thread, 
-           julia_time_1_thread, julia_time_multi_thread])'
-
-stds = reduce(hcat, [matlab_time_std, c_time_std, c_time_multi_thread_std, 
-               julia_time_1_thread_std, julia_time_multi_thread_std])'
-
-ix1, ix2 = [1], [2]
-means1, means2 = means[:, ix1], means[:, ix2]
-
-scaler = maximum(means1) / maximum(means2)
-
-means[:, ix2] .*= scaler
-stds[:, ix2] .*= scaler
-
-labels = ["MATLAB", "Gen. C, 1 thread", "Gen. C, 6 threads", "BeforeIT.jl, 1 thread", "BeforeIT.jl, 6 threads"]
-
-p3 = groupedbar(means, yerr = stds, ylabel="mean time for one epoch (s)", 
-                labels=["\$8\\cdot10^3\\textrm{\\:agents}\$" "\$8\\cdot10^6\\textrm{\\:agents}\$"],
-                xticks=(1:size(means,1), labels), foreground_color_legend=nothing, xtickfont = font(5),
-                ytickfont = font(8), guidefont = font(8), grid=false, dpi=1200);
-
-plot!(twinx(),  ylims=ylims(p3)./scaler, ylabel="",
-      xtickfont = font(5), ytickfont = font(8), guidefont = font(8))
-
-# Save the image
-savefig("benchmark_w_matlab.png")
-
-# The Julia implementation is faster than the MATLAB implementation, and the multi-threaded version is
-# faster than the single-threaded version.
+julia_times_small_multi ./= T
+julia_times_big_multi = [35.775353, 34.933283, 35.594063, 35.469412, 35.521846]
+julia_mtime_small_multi = mean(julia_times_small_multi)
+julia_stime_small_multi = std(julia_times_small_multi)
+julia_mtime_big_multi = mean(julia_times_big_multi)
+julia_stime_big_multi = std(julia_times_big_multi)
+
+
+labels = ["MATLAB", "Gen. C, 1 core", "Gen. C, 4 cores", "HPC, 1 core*", "HPC, 4 cores*", 
+          "BeforeIT.jl, 1 core", "BeforeIT.jl, 4 cores"]
+
+# Create the layout
+fig = Figure(size = (800, 400));
+
+ax1 = Axis(fig[1, 1], ylabel="time for one epoch (s)", title="Model with 8 thousands agents")
+ax2 = Axis(fig[1, 2], title="Model with 8 millions agents")
+
+times_small = [matlab_mtime_small, c_mtime_small, c_mtime_small_multi, julia_mtime_small, julia_mtime_small_multi]
+barplot!(ax1,
+    1:5,
+    times_small,
+    bar_labels = :y,
+);
+ylims!(ax1, 0, 1.15*maximum(times_small))
+
+ax1.yticklabelspace = 25.0
+ax1.xticks = (1:5, [labels[1:3]..., labels[6:end]...])
+ax1.xticklabelrotation = 45.0
+ax1.xgridvisible = false
+
+times_big = [matlab_mtime_big, c_mtime_big, c_mtime_big_multi, hpc_mtime_big, hpc_mtime_big_multi, 
+             julia_mtime_big, julia_mtime_big_multi]
+
+ylims!(ax2, 0, 1.15*maximum(times_big))
+barplot!(ax2,
+    1:7,
+    round.(times_big, digits=1),
+    bar_labels = :y
+);
+ax2.xticks = (1:7, labels)
+ax2.xticklabelrotation = 45.0
+ax2.xgridvisible = false
+
+# Save or display the layout
+display(fig)
+
+save("benchmark_w_matlab.png", fig)