Use a dynamic sampler to speed-up weighted sampling

.github/workflows/documentation.yml

@@ -28,7 +28,7 @@ jobs:
       - uses: julia-actions/setup-julia@v2
         with:
           version: '1.10'
-      - uses: julia-actions/cache@v1
+      - uses: julia-actions/cache@v2
       - name: Install dependencies
         run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy

Project.toml

@@ -6,8 +6,10 @@ version = "0.1.2"
 [deps]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+DynamicSampling = "2083aeaf-6258-5d07-89fc-32cf5060c837"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
 MutableNamedTuples = "af6c499f-54b4-48cc-bbd2-094bba7533c7"
@@ -20,8 +22,10 @@ StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
 [compat]
 Dates = "1"
 Distributions = "0.25"
+DynamicSampling = "0.4"
 FileIO = "1.16"
 JLD2 = "0.4"
+LazyArrays = "2"
 LinearAlgebra = "1"
 MAT = "0.10"
 MutableNamedTuples = "0.1.3"

src/BeforeIT.jl

@@ -1,9 +1,12 @@
 module BeforeIT
 
+import Base: length
+
+using DynamicSampling
+using LazyArrays
 using LinearAlgebra
 using Random
 using StatsBase
-import Base: length
 
 # definition of agents
 include("model_init/agents.jl")

src/markets/search_and_matching.jl

@@ -1,5 +1,4 @@
 
-
 """
     search_and_matching!(model, multi_threading::Bool = false)
 
@@ -54,7 +53,6 @@ function search_and_matching!(model, multi_threading = false)
         S_fg = copy(S_f)
         S_fg_ = copy(S_f_)
 
-
         perform_firms_market!(
             g,
             firms,
@@ -144,9 +142,6 @@ function search_and_matching!(model, multi_threading = false)
 
 end
 
-
-
-
 function update_aggregate_variables!(
     agg,
     w_act,
@@ -238,20 +233,20 @@ function initialize_variables_retail_market(firms, rotw, prop, agg, w_act, w_ina
     # ... Initialize all the variables ...
 
     # change some variables according to arguments of matlab function
-    b_HH_g = agg.P_bar_g .* prop.products.b_HH_g / sum(agg.P_bar_g .* prop.products.b_HH_g) #prop.products.b_HH_g
-    b_CFH_g = agg.P_bar_g .* prop.products.b_CFH_g / sum(agg.P_bar_g .* prop.products.b_CFH_g)   #prop.products.b_CFH_g
-    c_G_g = agg.P_bar_g .* prop.products.c_G_g / sum(agg.P_bar_g .* prop.products.c_G_g)   #prop.products.c_G_g
-    c_E_g = agg.P_bar_g .* prop.products.c_E_g / sum(agg.P_bar_g .* prop.products.c_E_g)   #prop.products.c_E_g
+    b_HH_g = agg.P_bar_g .* prop.products.b_HH_g / sum(agg.P_bar_g .* prop.products.b_HH_g)    #prop.products.b_HH_g
+    b_CFH_g = agg.P_bar_g .* prop.products.b_CFH_g / sum(agg.P_bar_g .* prop.products.b_CFH_g) #prop.products.b_CFH_g
+    c_G_g = agg.P_bar_g .* prop.products.c_G_g / sum(agg.P_bar_g .* prop.products.c_G_g)       #prop.products.c_G_g
+    c_E_g = agg.P_bar_g .* prop.products.c_E_g / sum(agg.P_bar_g .* prop.products.c_E_g)       #prop.products.c_E_g
 
     G = size(agg.P_bar_g, 1)
 
     # retrieve some general lengths from existing arrays
-    I = size(firms.P_i, 1)            # number of firms
+    I = size(firms.P_i, 1)        # number of firms
     H_W = length(w_act)           # number of active households
     H_inact = length(w_inact)     # number of inactive households
-    H = H_W + H_inact + I + 1      # number of households
-    L = size(rotw.C_d_l, 1)        # number of export partners
-    J = size(gov.C_d_j, 1)       # number of government entities
+    H = H_W + H_inact + I + 1     # number of households
+    L = size(rotw.C_d_l, 1)       # number of export partners
+    J = size(gov.C_d_j, 1)        # number of government entities
 
 
     # define a global C_d_h and I_d_h
@@ -305,14 +300,14 @@ function initialize_variables_firms_market(firms, rotw, prop)
     b_CF_g = prop.products.b_CF_g
 
     G = length(prop.products.b_HH_g) # number of goods
-    I = length(firms)            # number of firms
+    I = length(firms)                # number of firms
 
     # join internal and foreign firms arrays
-    P_f = [firms.P_i; rotw.P_m]                    # price array (firms + foreign firms)) 
-    S_f = [firms.Y_i + firms.S_i; rotw.Y_m]        # size array (firms + foreign firms)
+    P_f = [firms.P_i; rotw.P_m]                     # price array (firms + foreign firms)) 
+    S_f = [firms.Y_i + firms.S_i; rotw.Y_m]         # size array (firms + foreign firms)
     S_i_ = firms.K_i .* firms.kappa_i .- firms.Y_i  # (from matlab inputs)
     S_f_ = [S_i_; ones(size(rotw.Y_m)) .* Inf]      # Join S_i_ with an array of Infs of size(Y_m)
-    G_f = [firms.G_i; collect(1:G)]                # enlarge vector of final goods with foreign firms
+    G_f = [firms.G_i; collect(1:G)]                 # enlarge vector of final goods with foreign firms
 
     I_i_g = zeros(G, I)         # output
     P_CF_i_g = zeros(G, I)
@@ -322,6 +317,9 @@ function initialize_variables_firms_market(firms, rotw, prop)
     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
 
+"""
+Perform the firms market exchange process
+"""
 function perform_firms_market!(
     g,
     firms,
@@ -339,8 +337,6 @@ function perform_firms_market!(
     S_fg,
     S_fg_,
 )
-    # ... Perform the firms market exchange process ...
-
     ##############################
     ######## FIRMS MARKET ########
     ##############################
@@ -351,70 +347,54 @@ function perform_firms_market!(
     # firms that have demand for good "g" participate as buyers
     I_g = findall(DM_d_ig .> 0)
 
-    # remove firms that have no stock of good "g"
-    #F_g[S_fg[F_g] .<= 0] .= []
-    to_delete = findall(@view(S_fg[F_g]) .<= 0)
-    deleteat!(F_g, to_delete)
+    # keep firms that have positive stock of good "g"
+    filter!(i -> S_fg[i] > 0, F_g)
 
     # continue exchanges until either demand or supply terminates
 
-    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]))))
+    # weights according to size and price
+    F_g_sampler = create_weighted_sampler(P_f, S_f, F_g)
 
-        # size probability of being selected
-        pr_size_f = @view(S_f[F_g]) ./ sum(@view(S_f[F_g]))
-
-        # total probabilities of being selected
-        pr_cum_f_ = (pr_price_f + pr_size_f) ./ sum(pr_price_f + pr_size_f)
+    while !isempty(I_g) && !isempty(F_g_sampler)
 
         # select buyers at random
         shuffle!(I_g)
-        for j in eachindex(I_g)
-            i = I_g[j]
-
+        for i in I_g
             # select a random firm according to the probabilities
-            e = wsample(1:length(F_g), pr_cum_f_)
+            e = rand(F_g_sampler)
             f = F_g[e]
 
             # selected firm has sufficient stock
             if S_fg[f] > DM_d_ig[i]
                 S_fg[f] -= DM_d_ig[i]
-                DM_nominal_ig[i] += DM_d_ig[i] .* @view(P_f[f])
+                DM_nominal_ig[i] += DM_d_ig[i] * P_f[f]
                 DM_d_ig[i] = 0
             else
                 DM_d_ig[i] -= S_fg[f]
-                DM_nominal_ig[i] += @view(S_fg[f]) .* @view(P_f[f])
+                DM_nominal_ig[i] += S_fg[f] .* P_f[f]
                 S_fg[f] = 0
-                deleteat!(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_size_f = @view(S_f[F_g]) ./ sum(@view(S_f[F_g]))
-                pr_cum_f_ = (pr_price_f + pr_size_f) ./ sum(pr_price_f + pr_size_f)
+                delete!(F_g_sampler, e)
+                isempty(F_g_sampler) && break
             end
         end
-        I_g = findall(DM_d_ig .> 0)
+        I_g = I_g[@view(DM_d_ig[I_g]) .> 0]
     end
 
     if !isempty(I_g)
         DM_d_ig_ = copy(DM_d_ig)
+
         I_g = findall(DM_d_ig_ .> 0)
         F_g = findall(G_f .== g)
+        filter!(i -> S_fg_[i] > 0 && S_f[i] > 0, F_g)
 
-        to_delete = findall((@view(S_fg_[F_g]) .<= 0.0) .|| (@view(S_f[F_g]) .<= 0.0))
-        deleteat!(F_g, to_delete)
+        # weights according to size and price
+        F_g_sampler = create_weighted_sampler(P_f, S_f, F_g)
 
-        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_size_f = @view(S_f[F_g]) ./ sum(@view(S_f[F_g]))
-            pr_cum_f_ = (pr_price_f + pr_size_f) ./ sum(pr_price_f + pr_size_f)
+        while !isempty(I_g) && !isempty(F_g_sampler)
 
             shuffle!(I_g)
-            for j in eachindex(I_g)
-                i = I_g[j]
-
-                e = wsample(1:length(F_g), pr_cum_f_)
+            for i in I_g
+                e = rand(F_g_sampler)
                 f = F_g[e]
 
                 if S_fg_[f] > DM_d_ig_[i]
@@ -425,29 +405,30 @@ function perform_firms_market!(
                     DM_d_ig_[i] -= S_fg_[f]
                     S_fg[f] -= S_fg_[f]
                     S_fg_[f] = 0
-                    deleteat!(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_size_f = @view(S_f[F_g]) ./ sum(@view(S_f[F_g]))
-                    pr_cum_f_ = (pr_price_f + pr_size_f) ./ sum(pr_price_f + pr_size_f)
+                    delete!(F_g_sampler, e)
+                    isempty(F_g_sampler) && break
                 end
             end
-            I_g = findall(DM_d_ig_ .> 0)
+            I_g = I_g[@view(DM_d_ig_[I_g]) .> 0]
         end
     end
 
-    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
+    F_g = F_g[allinds(F_g_sampler)]
 
-    DM_i_g[g, :] .= a
+    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, :] .= b
+    @~ DM_i_g[g, :] .= a
+    @~ 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)
+    @~ P_bar_i_g[g, :] .= pos.(DM_nominal_ig .* a ./ c)
+    @~ P_CF_i_g[g, :] .= pos.(DM_nominal_ig .* b ./ c)
 end
 
+"""
+Perform the retail market exchange process
+"""
 function perform_retail_market!(
     g,
     agg,
@@ -480,11 +461,10 @@ function perform_retail_market!(
     S_f,
     G_f,
 )
-    # ... Perform the retail market exchange process ...
-
     ###############################
     ######## RETAIL MARKET ########
     ###############################
+
     C_d_hg = [
         b_HH_g[g] .* C_d_h .+ b_CFH_g[g] .* I_d_h
         c_E_g[g] .* rotw.C_d_l
@@ -493,19 +473,16 @@ function perform_retail_market!(
     C_real_hg = zeros(size(C_d_hg))
     H_g = findall(C_d_hg .> 0.0)
 
-    to_delete = findall(@view(S_fg[F_g]) .<= 0)
-    deleteat!(F_g, to_delete)
+    filter!(i -> S_fg[i] > 0, F_g)
 
-    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_size_f = @view(S_f[F_g]) ./ sum(@view(S_f[F_g]))
-        pr_cum_f_ = (pr_price_f + pr_size_f) ./ sum(pr_price_f + pr_size_f)
+    # weights according to size and price
+    F_g_sampler = create_weighted_sampler(P_f, S_f, F_g)
 
-        shuffle!(H_g)
-        for j in eachindex(H_g)
-            h = H_g[j]
+    while !isempty(H_g) && !isempty(F_g_sampler)
 
-            e = wsample(1:length(F_g), pr_cum_f_) # SLOW
+        shuffle!(H_g)
+        for h in H_g
+            e = rand(F_g_sampler)
             f = F_g[e]
 
             if S_fg[f] > C_d_hg[h] / P_f[f]
@@ -516,30 +493,28 @@ function perform_retail_market!(
                 C_d_hg[h] -= S_fg[f] * P_f[f]
                 C_real_hg[h] += S_fg[f]
                 S_fg[f] = 0
-                deleteat!(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_size_f = @view(S_f[F_g]) ./ sum(@view(S_f[F_g]))
-                pr_cum_f_ = (pr_price_f + pr_size_f) ./ sum(pr_price_f + pr_size_f)
+                delete!(F_g_sampler, e)
+                isempty(F_g_sampler) && break
             end
         end
-        H_g = findall(C_d_hg .> 0)
+        H_g = H_g[@view(C_d_hg[H_g]) .> 0]
     end
 
     if !isempty(H_g)
         C_d_hg_ = copy(C_d_hg)
+
         H_g = findall(C_d_hg_ .> 0)
         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_size_f = @view(S_f[F_g]) ./ sum(@view(S_f[F_g]))
-            pr_cum_f_ = (pr_price_f + pr_size_f) ./ sum(pr_price_f + pr_size_f)
-
-            H_g = shuffle(H_g)
-            for j in eachindex(H_g)
-                h = H_g[j]
-                e = wsample(1:length(F_g), pr_cum_f_)
+        filter!(i -> S_fg_[i] > 0 && S_f[i] > 0, F_g)
+
+        # weights according to size and price
+        F_g_sampler = create_weighted_sampler(P_f, S_f, F_g)
+
+        while !isempty(H_g) && !isempty(F_g_sampler)
+
+            shuffle!(H_g)
+            for h in H_g
+                e = rand(F_g_sampler)
                 f = F_g[e]
 
                 if S_fg_[f] > C_d_hg_[h] / P_f[f]
@@ -550,34 +525,55 @@ function perform_retail_market!(
                     C_d_hg_[h] -= S_fg_[f] * P_f[f]
                     S_fg[f] -= S_fg_[f]
                     S_fg_[f] = 0
-                    deleteat!(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_size_f = @view(S_f[F_g]) ./ sum(@view(S_f[F_g]))
-                    pr_cum_f_ = (pr_price_f + pr_size_f) ./ sum(pr_price_f + pr_size_f)
+                    delete!(F_g_sampler, e)
+                    isempty(F_g_sampler) && break
                 end
             end
-            H_g = findall(C_d_hg_ .> 0)
+            H_g = H_g[@view(C_d_hg_[H_g]) .> 0]
         end
     end
 
+    F_g = F_g[allinds(F_g_sampler)]
+
     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]))
+    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])
+    @~ 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
-    I_h_g[g, :] .= d
+    @~ 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)]))
+    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)]))
 
     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
+
+function compute_price_size_weights(P_f, S_f, F_g)
+    # price probability of being selected
+    pr_price_f_v = @~ exp.(-2 .* @view(P_f[F_g]))
+    pr_price_f_sum = sum(pr_price_f_v)
+    pr_price_f = @~ pos.(pr_price_f_v ./ pr_price_f_sum)
+    # size probability of being selected
+    pr_size_f_v = @view(S_f[F_g])
+    pr_size_f_sum = sum(pr_size_f_v)
+    pr_size_f = @~ pr_size_f_v ./ pr_size_f_sum
+    # total weight of being selected
+    w_cum_f_ = @~ pr_price_f .+ pr_size_f
+    return w_cum_f_
+end
+
+function create_weighted_sampler(P_f, S_f, F_g)
+    sampler = DynamicSampler()
+    isempty(F_g) && return sampler
+    w_cum_f_ = compute_price_size_weights(P_f, S_f, F_g)
+    append!(sampler, 1:length(F_g), w_cum_f_)
+    return sampler
+end

src/markets/search_and_matching_labour.jl

@@ -25,6 +25,7 @@ function search_and_matching_labour(firms::AbstractFirms, model)
     O_h = model.w_act.O_h
 
     V_i = N_d_i .- N_i
+
     # get employed workers in random order
     H_E = findall(O_h .> 0)
     shuffle!(H_E)
@@ -46,25 +47,26 @@ function search_and_matching_labour(firms::AbstractFirms, model)
 
     # find unemployed workers and positive vacancies
     H_U = findall(O_h .== 0)
+    shuffle!(H_U)
     I_V = findall(V_i .> 0)
+    shuffle!(I_V)
 
     # while there are no more vacancies or unemployed workers
     while !isempty(H_U) && !isempty(I_V)
-        shuffle!(I_V)
-
         for f in eachindex(I_V)
-            i = I_V[f]   # select random vacancy
-            e = rand(1:length(H_U))
-            h = H_U[e]   # select random unemployed worker
-            O_h[h] = i   # employ worker
+            # select random vacancy
+            i = I_V[f]
+            # select random unemployed worker
+            h = H_U[1]
+            # employ worker
+            O_h[h] = i
             N_i[i] += 1
             V_i[i] -= 1
-            deleteat!(H_U, e)
-            if isempty(H_U)
-                break
-            end
+            popfirst!(H_U)
+            isempty(H_U) && break
         end
-        I_V = findall(V_i .> 0)
+        I_V = I_V[@view(V_i[I_V]) .> 0]
+        shuffle!(I_V)
     end
 
     return N_i, O_h

src/model_init/init_properties.jl

@@ -62,4 +62,4 @@ function init_properties(parameters::Dict{String, Any}, T::Integer; typeInt::Dat
     properties = recursive_namedtuple(properties)
 
     return properties
-end
+end

src/utils/positive.jl

@@ -29,12 +29,8 @@ function pos(vector::AbstractArray)
     return r
 end
 
-function pos(number::T) where {T <: Number}
-    if isnan(number) || number < zero(T)
-        return zero(T)
-    else
-        return number
-    end
+@inline function pos(x::T) where {T <: Number}
+    return isnan(x) || x < zero(T) ? zero(T) : x
 end
 
 """
@@ -81,12 +77,8 @@ function neg(vector::AbstractArray)
     return r
 end
 
-function neg(number::T) where {T <: Number}
-    if isnan(number) || number > zero(T)
-        return zero(T)
-    else
-        return number
-    end
+@inline function neg(x::T) where {T <: Number}
+    return isnan(x) || x > zero(T) ? zero(T) : x
 end
 
 """

test/Project.toml

@@ -3,6 +3,7 @@
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+DynamicSampling = "2083aeaf-6258-5d07-89fc-32cf5060c837"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"

test/deterministic/make_model_deterministic.jl

@@ -2,14 +2,26 @@ using BeforeIT, MAT
 using Test
 
 import BeforeIT: randpl, epsilon
+import DynamicSampling: DynamicSampler, IndexInfo, getlevel, allinds
 import Random: shuffle!, rand, randn
 import StatsBase: wsample
 using Distributions
 
+function allinds(s::DynamicSampler)
+    return sort!(reduce(vcat, s.level_buckets))
+end
+
 function randn()
     return 0.0
 end
 
+function rand(s::DynamicSampler)
+    idx = minimum(minimum.(s.level_buckets; init=typemax(Int)))
+    weight = s.weights[idx]
+    level = ceil(Int, log2(weight)) - first(s.level_inds) + 1
+    idx_in_level = findfirst(x -> x == idx, s.level_buckets[level])
+    return idx
+end
 function rand(n::UnitRange)
     return 1
 end

test/markets/search_and_matching.jl

@@ -4,35 +4,36 @@ using Random
 @testset "search and matching" begin
     Random.seed!(1)
 
-    dir = @__DIR__
-
     parameters = BeforeIT.AUSTRIA2010Q1.parameters
     initial_conditions = BeforeIT.AUSTRIA2010Q1.initial_conditions
 
     T = 1
     model = BeforeIT.init_model(parameters, initial_conditions, T;)
 
-
-    gov = model.gov # government
-    cb = model.cb # central bank
-    rotw = model.rotw # rest of the world
-    firms = model.firms # firms
-    bank = model.bank # bank
-    w_act = model.w_act # active workers
+    gov = model.gov         # government
+    cb = model.cb           # central bank
+    rotw = model.rotw       # rest of the world
+    firms = model.firms     # firms
+    bank = model.bank       # bank
+    w_act = model.w_act     # active workers
     w_inact = model.w_inact # inactive workers
-    agg = model.agg # aggregate variables
+    agg = model.agg         # aggregates
+    prop = model.prop       # model properties
+
+    BeforeIT.finance_insolvent_firms!(firms, bank, model)
 
-    prop = model.prop # model properties
+    agg.Y_e, agg.gamma_e, agg.pi_e = BeforeIT.growth_inflation_expectations(model)
 
+    agg.epsilon_Y_EA, agg.epsilon_E, agg.epsilon_I = BeforeIT.epsilon(prop.C)
 
-    gamma_e = 0.01 # set expected growth in euro area
-    pi_e = 0.001   # set expected inflation in euro area
+    rotw.Y_EA, rotw.gamma_EA, rotw.pi_EA = BeforeIT.growth_inflation_EA(rotw, model)
 
-    agg.gamma_e = gamma_e
-    agg.pi_e = pi_e
+    cb.r_bar = BeforeIT.central_bank_rate(cb, model)
+
+    bank.r = BeforeIT.bank_rate(bank, model)
 
     Q_s_i, I_d_i, DM_d_i, N_d_i, Pi_e_i, DL_d_i, K_e_i, L_e_i, P_i =
-        BeforeIT.firms_expectations_and_decisions(model.firms, model)
+        BeforeIT.firms_expectations_and_decisions(firms, model)
 
     firms.Q_s_i .= Q_s_i
     firms.I_d_i .= I_d_i
@@ -44,8 +45,20 @@ using Random
     firms.K_e_i .= K_e_i
     firms.L_e_i .= L_e_i
 
-    Pi_e_k = BeforeIT.bank_expected_profits(bank, model)
-    bank.Pi_e_k = Pi_e_k
+    firms.DL_i .= BeforeIT.search_and_matching_credit(firms, model)
+
+    N_i, Oh = BeforeIT.search_and_matching_labour(firms, model)
+    firms.N_i .= N_i
+    w_act.O_h .= Oh
+
+    firms.w_i .= BeforeIT.firms_wages(firms)
+    firms.Y_i .= BeforeIT.firms_production(firms)
+
+    BeforeIT.update_workers_wages!(w_act, firms.w_i)
+
+    gov.sb_other, gov.sb_inact = BeforeIT.gov_social_benefits(gov, model)
+
+    bank.Pi_e_k = BeforeIT.bank_expected_profits(bank, model)
 
     C_d_h, I_d_h = BeforeIT.households_budget_act(w_act, model)
     w_act.C_d_h .= C_d_h
@@ -62,13 +75,6 @@ using Random
     gov.C_G = C_G
     gov.C_d_j .= C_d_j
 
-
-    epsilon_E = 0.28
-    epsilon_I = 0.36
-
-    agg.epsilon_E = epsilon_E
-    agg.epsilon_I = epsilon_I
-
     C_E, Y_I, C_d_l, Y_m, P_m = BeforeIT.rotw_import_export(rotw, model)
     rotw.C_E = C_E
     rotw.Y_I = Y_I
@@ -78,27 +84,18 @@ using Random
 
     BeforeIT.search_and_matching!(model, false)
 
-    rtol = 0.0001
-
-    # NOTE: the expected numbers come out of the original implementation, 
-    # and only hold for the serial code (without multithreading)
-    @test isapprox(mean(w_act.C_h), 4.148850396106796, rtol = rtol)
-    @test isapprox(mean(w_inact.C_h), 2.205381003981018, rtol = rtol)
-    @test isapprox(mean(firms.C_h), 9.060799641122962, rtol = rtol)
-    @test isapprox(bank.C_h, 2931.5395701704915, rtol = rtol)
-
-    @test isapprox(mean(w_act.I_h), 0.34186063655926524, rtol = rtol)
-    @test isapprox(mean(w_inact.I_h), 0.18217582636296747, rtol = rtol)
-    @test isapprox(mean(firms.I_h), 0.7442975169996757, rtol = rtol)
-    @test isapprox(bank.I_h, 233.5381841004737, rtol = rtol)
-
-    @test isapprox(gov.C_j, 14686.094833493271, rtol = rtol)
-    @test isapprox(rotw.C_l, 44241.742486622454, rtol = rtol)
-
-    @test isapprox(mean(firms.I_i), 20.671016463479898, rtol = rtol)
-    @test isapprox(mean(firms.DM_i), 110.18635469222951, rtol = rtol)
-    @test isapprox(mean(firms.P_bar_i), 1.0010000000000023, rtol = rtol)
-    @test isapprox(mean(firms.P_CF_i), 1.0010000000000023, rtol = rtol)
-    @test isapprox(mean(firms.Q_d_i), 216.70740037345882, rtol = rtol)
-    @test isapprox(mean(rotw.Q_d_m), 719.2385742449192, rtol = rtol)
+    # NOTE: as a test we use the expected values and standard deviations of the
+    #       original implementation, with tolerance = 3*(standard deviation)
+    @test isapprox(mean([bank.I_h, w_act.I_h..., w_inact.I_h..., firms.I_h...]), 
+                   0.32975, atol = 3*0.0025351)
+    @test isapprox(mean([bank.C_h, w_act.C_h..., w_inact.C_h..., firms.C_h...]), 
+                   3.973, atol = 3*0.029366)
+    @test isapprox(mean(firms.I_i), 20.5075, atol = 3*0.12763)
+    @test isapprox(mean(firms.DM_i), 109.3163, atol = 3*0.68033)
+    @test isapprox(mean(firms.P_bar_i), 1.0031, atol = 3*0.0044726)
+    @test isapprox(mean(firms.P_CF_i), 1.0031, atol = 3*0.0044726)
+    @test isapprox(gov.C_j, 14752.2413, atol = 3*126.7441)
+    @test isapprox(rotw.C_l, 34188.1258, atol = 3*666.275)
+    @test isapprox(mean(firms.Q_d_i), 216.2474, atol = 3*1.2275)
+    @test isapprox(mean(rotw.Q_d_m), 535.7522, atol = 3*9.6082)
 end