new vaccination mechanics

CovidAlertVaccinationModel/abm_timeseries.jl

@@ -5,19 +5,21 @@ const samples = 1
 const parameters = (
     sim_length = 600,
     num_households = 5000,
-    I_0_fraction = 0.002,
+    I_0_fraction = 0.0,
     base_transmission_probability = 0.001,
     recovery_rate = 1/7,
     immunization_loss_prob = 0.0055, #mean time of 6 months
-    π_base = -0.3,
+    π_base = -1.0,
     η = 0.0,
-    κ = 0.8,
+    κ = 0.0,
     ω = 0.00005,
-    ρ = [0.0,0.0,0.0],
+    ρ_y = 0.0,
+    ρ_m = 0.0,
+    ρ_o = 0.0,
     ω_en = 0.00,
     ρ_en =  [0.0,0.0,0.0],
     γ = 0.0,
-    β = 10.0,
+    β = 5.0,
     notification_parameter = 0.001,
     vaccinator_prob = 0.2,
     app_user_fraction = 0.2,

CovidAlertVaccinationModel/src/ABM/solve.jl

@@ -89,95 +89,94 @@ Base.@propagate_inbounds @views function update_infection_state!(t,modelsol)
         end
     end
 end
-Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,modelsol,total_infections)
-    @unpack π_base, η,γ, κ, ω, ρ_y,ρ_m,ρ_o, ω_en,ρ_en,γ,β = modelsol.params
-    @unpack demographics,time_of_last_alert, nodes, soc_network,u_vac,u_next_vac,app_user,app_user_list = modelsol
-    app_user_pointer = 0
-    for i in 1:nodes
-        vac_payoff = 0
-        soc_nbrs_vac = @MArray [0,0,0]
-        soc_nbrs_nonvac = 0
-        num_soc_nbrs = 0
-        for sc_g in soc_network.graph_list[t]
-            soc_nbrs = neighbors(sc_g.g,i)
-            num_soc_nbrs += length(soc_nbrs)
-            for nbr in soc_nbrs
-                if u_vac[nbr]
-                    soc_nbrs_vac[Int(demographics[nbr])] += 1
-                else
-                    soc_nbrs_nonvac += 1
-                end
-            end
-        end
-        ρ = @SVector [ρ_y,ρ_m,ρ_o]
-
-        vac_payoff += π_base + dot(ρ,soc_nbrs_vac) + total_infections*ω +
-            ifelse(num_soc_nbrs> 0, κ * ((sum(soc_nbrs_vac) - soc_nbrs_nonvac)/num_soc_nbrs),0)
+# Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,modelsol,total_infections)
+#     @unpack π_base, η,γ, κ, ω, ρ_y,ρ_m,ρ_o, ω_en,ρ_en,γ,β = modelsol.params
+#     @unpack demographics,time_of_last_alert, nodes, soc_network,u_vac,u_next_vac,app_user,app_user_list = modelsol
+#     app_user_pointer = 0
+#     for i in 1:nodes
+#         vac_payoff = 0
+#         soc_nbrs_vac = @MArray [0,0,0]
+#         soc_nbrs_nonvac = 0
+#         num_soc_nbrs = 0
+#         for sc_g in soc_network.graph_list[t]
+#             soc_nbrs = neighbors(sc_g.g,i)
+#             num_soc_nbrs += length(soc_nbrs)
+#             for nbr in soc_nbrs
+#                 if u_vac[nbr]
+#                     soc_nbrs_vac[Int(demographics[nbr])] += 1
+#                 else
+#                     soc_nbrs_nonvac += 1
+#                 end
+#             end
+#         end
+#         ρ = @SVector [ρ_y,ρ_m,ρ_o]
+
+#         vac_payoff += π_base + dot(ρ,soc_nbrs_vac) + total_infections*ω +
+#             ifelse(num_soc_nbrs> 0, κ * ((sum(soc_nbrs_vac) - soc_nbrs_nonvac)/num_soc_nbrs),0)
         
-        if app_user[i] && time_of_last_alert[app_user_list[i]]>=0
-            vac_payoff += γ^(-1*(t - time_of_last_alert[app_user_list[i]]))* (η + dot(ρ_en,soc_nbrs_vac) + total_infections*ω_en) 
-        end
+#         if app_user[i] && time_of_last_alert[app_user_list[i]]>=0
+#             vac_payoff += γ^(-1*(t - time_of_last_alert[app_user_list[i]]))* (η + dot(ρ_en,soc_nbrs_vac) + total_infections*ω_en) 
+#         end
         
-        if u_vac[i]
-            if rand(Random.default_rng(Threads.threadid())) < 1 - Φ(vac_payoff,β)
-                u_next_vac[i] = !u_vac[i]
-            else
-                u_next_vac[i] = u_vac[i]
-            end
-        else
-            if rand(Random.default_rng(Threads.threadid())) < Φ(vac_payoff,β)
-                u_next_vac[i] = !u_vac[i]
-            else
-                u_next_vac[i] = u_vac[i]
-            end
-        end
-    end
+#         if u_vac[i]
+#             if rand(Random.default_rng(Threads.threadid())) < 1 - Φ(vac_payoff,β)
+#                 u_next_vac[i] = !u_vac[i]
+#             else
+#                 u_next_vac[i] = u_vac[i]
+#             end
+#         else
+#             if rand(Random.default_rng(Threads.threadid())) < Φ(vac_payoff,β)
+#                 u_next_vac[i] = !u_vac[i]
+#             else
+#                 u_next_vac[i] = u_vac[i]
+#             end
+#         end
+#     end
+
+#     modelsol.daily_vaccinators = count(==(true),u_vac) #could maybe make this more efficient
+# end
 
-    modelsol.daily_vaccinators = count(==(true),u_vac) #could maybe make this more efficient
-end
 
-
-Base.@propagate_inbounds @views function update_vaccination_opinion_state_2!(t,modelsol,total_infections)
+Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,modelsol,total_infections)
     @unpack π_base, η,γ, κ, ω, ρ_y,ρ_m,ρ_o, ω_en,ρ_en,γ,β = modelsol.params
     @unpack demographics,time_of_last_alert, nodes, soc_network,u_vac,u_next_vac,app_user,app_user_list = modelsol
     app_user_pointer = 0
     for i in 1:nodes
         vac_payoff = 0
-        soc_nbrs = @MArray [0,0,0]
+        soc_nbrs_count = @MArray [0,0,0]
         num_soc_nbrs = 0
         for sc_g in soc_network.graph_list[t]
             soc_nbrs = neighbors(sc_g.g,i)
             num_soc_nbrs += length(soc_nbrs)
             for nbr in soc_nbrs
-                soc_nbrs[Int(demographics[nbr])] += 1
+                soc_nbrs_count[Int(demographics[nbr])] += 1
             end
         end
         ρ = @SVector [ρ_y,ρ_m,ρ_o]
-        # random_neighbour = 
+        random_soc_network = sample(Random.default_rng(Threads.threadid()), soc_network.graph_list[t])
+        # display(random_soc_network)
 
+        if !isempty(neighbors(random_soc_network.g,i))
+            random_neighbour = sample(Random.default_rng(Threads.threadid()), neighbors(random_soc_network.g,i))
+
+            if u_vac[random_neighbour] == u_vac[i]
+                vac_payoff += π_base + dot(ρ,soc_nbrs_count) + total_infections*ω             
+                if app_user[i] && time_of_last_alert[app_user_list[i]]>=0
+                    vac_payoff += γ^(-1*(t - time_of_last_alert[app_user_list[i]]))* (η + dot(ρ_en,soc_nbrs_count) + total_infections*ω_en) 
+                end
         
-            vac_payoff += π_base + dot(ρ,soc_nbrs_vac) + total_infections*ω +
-                ifelse(num_soc_nbrs> 0, κ * ((sum(soc_nbrs_vac) - soc_nbrs_nonvac)/num_soc_nbrs),0)
-            
-            if app_user[i] && time_of_last_alert[app_user_list[i]]>=0
-                vac_payoff += γ^(-1*(t - time_of_last_alert[app_user_list[i]]))* (η + dot(ρ_en,soc_nbrs_vac) + total_infections*ω_en) 
-            end
-        
-        if u_vac[i]
-            if rand(Random.default_rng(Threads.threadid())) < 1 - Φ(vac_payoff,β)
-                u_next_vac[i] = !u_vac[i]
-            else
-                u_next_vac[i] = u_vac[i]
-            end
-        else
-            if rand(Random.default_rng(Threads.threadid())) < Φ(vac_payoff,β)
-                u_next_vac[i] = !u_vac[i]
-            else
-                u_next_vac[i] = u_vac[i]
+                if u_vac[i]
+                    if rand(Random.default_rng(Threads.threadid())) < 1 - Φ(vac_payoff,β)
+                        u_next_vac[i] = false
+                    end
+                else
+                    if rand(Random.default_rng(Threads.threadid())) < Φ(vac_payoff,β)
+                        u_next_vac[i] = true
+                    end
+                end
             end
         end
     end
-
     modelsol.daily_vaccinators = count(==(true),u_vac) #could maybe make this more efficient
 end
 

CovidAlertVaccinationModel/src/data.jl

@@ -175,10 +175,10 @@ function filter_priors(location_key)
     end
 
     priors_dict = filter(row-> row["location"] == location_key, eachrow(priors)) |>
-                    df -> map(r -> (r[:index_out],r[:index_in]) => convert(Vector,r[string.(distribution_support)]),df) |> #this is pretty inflexible to changing support unfortunately
+                    df -> map(r -> (r[:index_out],r[:index_in]) =>Vector{Float64}(r[string.(distribution_support)]),df) |> #this is pretty inflexible to changing support unfortunately
                     Dict
 
-    # display(priors_dict)
+    display(priors_dict)
 
     priors_probs = Dict(map(v -> v => DiscreteNonParametric(distribution_support,priors_dict[v]), collect(keys(priors_dict))))