diff --git a/examples/benchmark_w_matlab.jl b/examples/benchmark_w_matlab.jl index 9004e28..367bfc6 100644 --- a/examples/benchmark_w_matlab.jl +++ b/examples/benchmark_w_matlab.jl @@ -7,11 +7,7 @@ using BeforeIT, Plots, Statistics -parameters = BeforeIT.AUSTRIA2010Q1.parameters -initial_conditions = BeforeIT.AUSTRIA2010Q1.initial_conditions -T = 12*2 - -function run(; multi_threading = false) +function run(parameters, initial_conditions, T; multi_threading = false) model = BeforeIT.init_model(parameters, initial_conditions, T) data = BeforeIT.init_data(model); @@ -22,28 +18,42 @@ function run(; multi_threading = false) return model, data end +parameters = BeforeIT.AUSTRIA2010Q1.parameters +initial_conditions = BeforeIT.AUSTRIA2010Q1.initial_conditions +T = 12 + # we run the code to compile it first -@time run(); -@time run(;multi_threading = true); +@time run(parameters, initial_conditions, T; multi_threading = false); +@time run(parameters, initial_conditions, T; multi_threading = true); -# time taken by the MATLAB code, computed independently on an Apple M1 chip -matlab_times = [3.1919, 3.2454, 3.1501, 3.1074, 3.1551] +# 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_time = mean(matlab_times) matlab_time_std = std(matlab_times) -# time taken by the Julia code, computed as the average of 5 runs +c_times = [0.952, 0.940, 0.951, 0.942, 0.938] +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_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), +# computed as the average of 5 runs n_runs = 5 julia_times_1_thread = zeros(n_runs) for i in 1:n_runs - julia_times_1_thread[i] = @elapsed run(); + julia_times_1_thread[i] = @elapsed run(parameters, initial_conditions, T; multi_threading = false); end 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) for i in 1:5 - julia_times_multi_thread[i] = @elapsed run(multi_threading = true); + julia_times_multi_thread[i] = @elapsed run(parameters, initial_conditions, T; multi_threading = true); end julia_time_multi_thread = mean(julia_times_multi_thread) julia_time_multi_thread_std = std(julia_times_multi_thread) @@ -55,14 +65,23 @@ theme(:default, bg = :white) # bar chart of the time taken vs the time taken by the MATLAB code, also plot the stds as error bars # make a white background with no grid -bar(["MATLAB", "Julia, 1 thread", "Julia, $n_threads threads"], [matlab_time, julia_time_1_thread, julia_time_multi_thread], -yerr = [matlab_time_std, julia_time_1_thread_std, julia_time_multi_thread_std], -legend = false, dpi=300, size=(400, 300), grid = false, ylabel = "Time for one simulation (s)") +bar( + ["MATLAB", "Gen. C, 1 thread", "Gen. C, 6 threads", "Julia, 1 thread", "Julia, $n_threads threads"], + [matlab_time, c_time, c_time_multi_thread, julia_time_1_thread, julia_time_multi_thread], + yerr = [matlab_time_std, c_time_std, c_time_multi_thread_std, julia_time_1_thread_std, julia_time_multi_thread_std], + legend = false, + dpi = 300, + size = (400, 300), + grid = false, + ylabel = "Time for one simulation (s)", + xtickfont = font(4), + ytickfont = font(6), + guidefont = font(6) +) + # the Julia implementation is faster than the MATLAB implementation, and the multi-threaded version is # faster than the single-threaded version. -# increase - # save the image savefig("benchmark_w_matlab.png") diff --git a/src/markets/search_and_matching.jl b/src/markets/search_and_matching.jl index 5cabe0d..84023bd 100644 --- a/src/markets/search_and_matching.jl +++ b/src/markets/search_and_matching.jl @@ -320,7 +320,6 @@ function initialize_variables_firms_market(firms, rotw, prop) P_bar_i_g = zeros(G, I) return a_sg, b_CF_g, P_f, S_f, S_f_, G_f, I_i_g, DM_i_g, P_bar_i_g, P_CF_i_g - end function perform_firms_market!( @@ -362,7 +361,7 @@ function perform_firms_market!( while length(I_g) != 0 && length(F_g) != 0 # price probability of being selected - pr_price_f = pos(exp.(-2 .* @view(P_f[F_g])) ./ sum(exp.(-2 .* @view(P_f[F_g])))) + pr_price_f = pos!(exp.(-2 .* @view(P_f[F_g])) ./ sum(exp.(-2 .* @view(P_f[F_g])))) # size probability of being selected pr_size_f = @view(S_f[F_g]) ./ sum(@view(S_f[F_g])) @@ -390,7 +389,7 @@ function perform_firms_market!( S_fg[f] = 0 swap_pop!(F_g, e) isempty(F_g) && break - pr_price_f = pos(exp.(-2 .* @view(P_f[F_g])) ./ sum(exp.(-2 .* @view(P_f[F_g])))) + pr_price_f = pos!(exp.(-2 .* @view(P_f[F_g])) ./ sum(exp.(-2 .* @view(P_f[F_g])))) pr_size_f = @view(S_f[F_g]) ./ sum(@view(S_f[F_g])) pr_cum_f_ = pr_price_f + pr_size_f end @@ -407,7 +406,7 @@ function perform_firms_market!( deleteat!(F_g, to_delete) while !isempty(I_g) && !isempty(F_g) - pr_price_f = pos(exp.(-2 .* @view(P_f[F_g])) ./ sum(exp.(-2 .* @view(P_f[F_g])))) + pr_price_f = pos!(exp.(-2 .* @view(P_f[F_g])) ./ sum(exp.(-2 .* @view(P_f[F_g])))) pr_size_f = @view(S_f[F_g]) ./ sum(@view(S_f[F_g])) pr_cum_f_ = pr_price_f + pr_size_f @@ -428,7 +427,7 @@ function perform_firms_market!( S_fg_[f] = 0 swap_pop!(F_g, e) isempty(F_g) && break - pr_price_f = pos(exp.(-2 .* @view(P_f[F_g])) ./ sum(exp.(-2 .* @view(P_f[F_g])))) + pr_price_f = pos!(exp.(-2 .* @view(P_f[F_g])) ./ sum(exp.(-2 .* @view(P_f[F_g])))) pr_size_f = @view(S_f[F_g]) ./ sum(@view(S_f[F_g])) pr_cum_f_ = pr_price_f + pr_size_f end @@ -437,20 +436,16 @@ function perform_firms_market!( end end - DM_i_g[g, :] .= @view(a_sg[g, firms.G_i]) .* firms.DM_d_i .- pos(DM_d_ig .- b_CF_g[g] .* firms.I_d_i) + a = @view(a_sg[g, firms.G_i]) .* firms.DM_d_i .- pos!(DM_d_ig .- b_CF_g[g] .* firms.I_d_i) + b = pos!(b_CF_g[g] .* firms.I_d_i .- DM_d_ig) + c = @view(a_sg[g, firms.G_i]) .* firms.DM_d_i .+ b_CF_g[g] .* firms.I_d_i .- DM_d_ig - I_i_g[g, :] .= pos(b_CF_g[g] .* firms.I_d_i .- DM_d_ig) + DM_i_g[g, :] .= a - P_bar_i_g[g, :] .= pos( - DM_nominal_ig .* (@view(a_sg[g, firms.G_i]) .* firms.DM_d_i .- pos(DM_d_ig .- b_CF_g[g] .* firms.I_d_i)) ./ - (@view(a_sg[g, firms.G_i]) .* firms.DM_d_i .+ b_CF_g[g] .* firms.I_d_i .- DM_d_ig), - ) - - P_CF_i_g[g, :] .= pos( - DM_nominal_ig .* pos(b_CF_g[g] .* firms.I_d_i .- DM_d_ig) ./ - (@view(a_sg[g, firms.G_i]) .* firms.DM_d_i .+ b_CF_g[g] .* firms.I_d_i .- DM_d_ig), - ) + I_i_g[g, :] .= b + P_bar_i_g[g, :] .= pos!(DM_nominal_ig .* a ./ c) + P_CF_i_g[g, :] .= pos!(DM_nominal_ig .* b ./ c) end function perform_retail_market!( @@ -501,9 +496,8 @@ function perform_retail_market!( to_delete = findall(@view(S_fg[F_g]) .<= 0) deleteat!(F_g, to_delete) - while !isempty(H_g) && !isempty(F_g) - pr_price_f = pos(exp.(-2 .* @view(P_f[F_g])) ./ sum(exp.(-2 .* @view(P_f[F_g])))) + pr_price_f = pos!(exp.(-2 .* @view(P_f[F_g])) ./ sum(exp.(-2 .* @view(P_f[F_g])))) pr_size_f = @view(S_f[F_g]) ./ sum(@view(S_f[F_g])) pr_cum_f_ = pr_price_f + pr_size_f @@ -524,7 +518,7 @@ function perform_retail_market!( S_fg[f] = 0 swap_pop!(F_g, e) isempty(F_g) && break - pr_price_f = pos(exp.(-2 .* @view(P_f[F_g])) ./ sum(exp.(-2 .* @view(P_f[F_g])))) + pr_price_f = pos!(exp.(-2 .* @view(P_f[F_g])) ./ sum(exp.(-2 .* @view(P_f[F_g])))) pr_size_f = @view(S_f[F_g]) ./ sum(@view(S_f[F_g])) pr_cum_f_ = pr_price_f + pr_size_f end @@ -538,7 +532,7 @@ function perform_retail_market!( F_g = findall(G_f .== g) F_g = F_g[(@view(S_fg_[F_g]) .> 0) .& (@view(S_f[F_g]) .> 0)] while !isempty(H_g) && !isempty(F_g) - pr_price_f = pos(exp.(-2 .* @view(P_f[F_g])) ./ sum(exp.(-2 .* @view(P_f[F_g])))) + pr_price_f = pos!(exp.(-2 .* @view(P_f[F_g])) ./ sum(exp.(-2 .* @view(P_f[F_g])))) pr_size_f = @view(S_f[F_g]) ./ sum(@view(S_f[F_g])) pr_cum_f_ = pr_price_f + pr_size_f @@ -558,8 +552,7 @@ function perform_retail_market!( S_fg_[f] = 0 swap_pop!(F_g, e) isempty(F_g) && break - pr_price_f = max.(0, exp.(-2 .* @view(P_f[F_g])) ./ sum(exp.(-2 .* @view(P_f[F_g])))) - pr_price_f[isnan.(pr_price_f)] .= 0.0 + pr_price_f = pos!(exp.(-2 .* @view(P_f[F_g])) ./ sum(exp.(-2 .* @view(P_f[F_g])))) pr_size_f = @view(S_f[F_g]) ./ sum(@view(S_f[F_g])) pr_cum_f_ = pr_price_f + pr_size_f end @@ -568,27 +561,23 @@ function perform_retail_market!( end end + a = @view(C_real_hg[1:H]) + b = C_d_h .* b_HH_g[g] .- pos!(@view(C_d_hg[1:H]) .- b_CFH_g[g] .* I_d_h) + c = C_d_h .* b_HH_g[g] .+ b_CFH_g[g] .* I_d_h .- @view(C_d_hg[1:H]) + d = pos!(b_CFH_g[g] .* I_d_h .- @view(C_d_hg[1:H])) + Q_d_i_g[g, :] .= @view(S_f[1:I]) .- @view(S_fg[1:I]) Q_d_m_g[g, :] .= @view(S_f[(I + 1):end]) .- @view(S_fg[(I + 1):end]) - C_h_g[g, :] .= b_HH_g[g] .* C_d_h .- pos(@view(C_d_hg[1:H]) .- b_CFH_g[g] .* I_d_h) - I_h_g[g, :] .= pos(b_CFH_g[g] .* I_d_h .- @view(C_d_hg[1:H])) + C_h_g[g, :] .= b + I_h_g[g, :] .= d C_j_g[g] = sum(c_G_g[g] .* gov.C_d_j) - sum(@view(C_d_hg[(H + L + 1):(H + L + J)])) C_l_g[g] = sum(c_E_g[g] .* rotw.C_d_l) - sum(@view(C_d_hg[(H + 1):(H + L)])) - a = sum(@view(C_real_hg[1:H])) - b = sum(C_d_h .* b_HH_g[g] .- pos(@view(C_d_hg[1:H]) .- b_CFH_g[g] .* I_d_h)) - c = sum((C_d_h .* b_HH_g[g] .+ b_CFH_g[g] .* I_d_h .- @view(C_d_hg[1:H]))) - - P_bar_h_g[g] = pos(a * b / c) - - P_bar_CF_h_g[g] = pos( - sum(@view(C_real_hg[1:H])) * sum(pos(b_CFH_g[g] .* I_d_h .- @view(C_d_hg[1:H]))) / - sum((C_d_h .* b_HH_g[g] .+ b_CFH_g[g] .* I_d_h .- @view(C_d_hg[1:H]))), - ) + P_bar_h_g[g] = pos(sum(a) * sum(b) / sum(c)) + P_bar_CF_h_g[g] = pos(sum(a) * sum(d) / sum(c)) P_j_g[g] = sum(@view(C_real_hg[(H + L + 1):(H + L + J)])) P_l_g[g] = sum(@view(C_real_hg[(H + 1):(H + L)])) - end diff --git a/src/utils/positive.jl b/src/utils/positive.jl index 9f9ff66..ee047e3 100644 --- a/src/utils/positive.jl +++ b/src/utils/positive.jl @@ -37,6 +37,19 @@ function pos(number::T) where {T <: Number} end end +""" + pos!(vector) -> vector + +In-place version of `pos`. Mimicks max(0, vector) in Matlab. +Returns the updated vector. +""" +function pos!(A) + @simd for i in eachindex(A) + A[i] = ifelse(isnan(A[i]), zero(eltype(A)), max(zero(eltype(A)), A[i])) + end + return A +end + """ neg(vector) -> vector @@ -76,6 +89,19 @@ function neg(number::T) where {T <: Number} end end +""" + neg!(vector) -> vector + +In-place version of `neg`. Mimicks min(0, vector) in Matlab. +Returns the updated vector. +""" +function neg!(A) + @simd for i in eachindex(A) + A[i] = ifelse(isnan(A[i]), zero(eltype(A)), min(zero(eltype(A)), A[i])) + end + return A +end + # like in the original code function matlab_round(x) return Base.round(x, RoundNearestTiesUp) diff --git a/test/calibration_script.jl b/test/calibration_script.jl index ddd4f6f..8c93d36 100644 --- a/test/calibration_script.jl +++ b/test/calibration_script.jl @@ -1,14 +1,12 @@ -using BeforeIT, MAT -using FileIO, Dates + +using MAT, FileIO, Dates dir = @__DIR__ parameters_mat = matread(joinpath(dir, "matlab_code/italy_calibration/parameters/2010Q1.mat")) initial_conditions_mat = matread(joinpath(dir, "matlab_code/italy_calibration/initial_conditions/2010Q1.mat")) -# parameters = BeforeIT.ITALY2010Q1.parameters -# initial_conditions = BeforeIT.ITALY2010Q1.initial_conditions - +using BeforeIT calibration_data = BeforeIT.ITALY_CALIBRATION.calibration figaro = BeforeIT.ITALY_CALIBRATION.figaro @@ -23,13 +21,12 @@ estimation_date = DateTime(1996, 12, 31) # Calibrate on a specific quarter calibration_date = DateTime(2010, 03, 31)#-Dates.Month(3) parameters, initial_conditions = BeforeIT.get_params_and_initial_conditions( - calibration_date, - calibration_data, - figaro, - data, - ea, - max_calibration_date, - estimation_date; + (calibration = calibration_data, + figaro = figaro, + data = data, + ea = ea, + max_calibration_date = max_calibration_date, + estimation_date = estimation_date), calibration_date; scale = 0.001, )