mixing graph bug fix

CovidAlertVaccinationModel/abm_timeseries.jl

 using CovidAlertVaccinationModel
 using OnlineStats
 using Plots
-const samples = 5
+const samples = 1
 
 const vaccination_data = [0.0,0.043,0.385,0.424,0.115,0.03,0.005] #by month starting in august
 const ymo_vac = [0.255,0.278,0.602]
@@ -32,4 +32,5 @@ function solve_and_plot_parameters()
     return out
 end
 
+
 out = solve_and_plot_parameters()

CovidAlertVaccinationModel/prem_mixing_approx.jl

@@ -4,23 +4,24 @@ function approximate_mixing_matricies_1(sol)
     mean_mixing_degree = zeros(3,3)
     mean_mixing_weighted_degree = zeros(3,3)
     display(p.O_distribution_shift)
-    te = 0
+    te_missed = zeros(8)
+    te = zeros(8)
+        edgelist = []
         for g_list in sol.inf_network.graph_list
-            for g in g_list
+            for (k,g) in enumerate(g_list)
                 # display(g.mixing_edges.weights_dict)
                 for e in keys(g.mixing_edges.weights_dict)
                     demo_i = sol.demographics[e.a]
                     demo_j = sol.demographics[e.b]
                     if(has_edge(g.g,e.a,e.b))
-                        mean_mixing_degree[Int(demo_i), Int(demo_j)] += 1#get_weight(g,GraphEdge(node_i,node_j)) /2
-                        mean_mixing_degree[Int(demo_j), Int(demo_i)] += 1#get_weight(g,GraphEdge(node_i,node_j)) /2
+                        mean_mixing_degree[Int(demo_i), Int(demo_j)] += 1
+                        mean_mixing_degree[Int(demo_j), Int(demo_i)] += 1
                         mean_mixing_weighted_degree[Int(demo_i), Int(demo_j)] += get_weight(g,GraphEdge(e.a,e.b))
                         mean_mixing_weighted_degree[Int(demo_j), Int(demo_i)] += get_weight(g,GraphEdge(e.a,e.b))
-                        
-                    else
-                        if(demo_i == demo_j)
-                            te += 1
-                        end
+                        te[k] += 1
+                    else 
+                        push!(edgelist,e)
+                        te_missed[k] += 1
                     end
                 end 
              end
@@ -29,12 +30,13 @@ function approximate_mixing_matricies_1(sol)
     mean_mixing_weighted_degree ./= (500 .* length.(sol.index_vectors) * 2)
     display(mean_mixing_degree)
     display(mean_mixing_weighted_degree)
-    display(te)
+    display(te_missed)
     return sol
 end
 function approximate_mixing_matricies_2(sol)
 mean_mixing_degree = zeros(3,3)
 mean_mixing_weighted_degree = zeros(3,3)
+te_missed = zeros(8)
     display(p.O_distribution_shift)
         for g_list in sol.inf_network.graph_list
             for g in g_list
@@ -51,8 +53,6 @@ mean_mixing_weighted_degree = zeros(3,3)
     end
     mean_mixing_degree ./= (500 .* length.(sol.index_vectors) * 2)
     mean_mixing_weighted_degree ./= (500 .* length.(sol.index_vectors) * 2)
-    display(mean_mixing_degree)
-    display(mean_mixing_weighted_degree)
     return sol
 end
 

CovidAlertVaccinationModel/src/ABM/abm.jl

@@ -40,7 +40,7 @@ end
 Run the model with given parameter tuple and output recorder. See `get_parameters` for list of parameters. See `output.jl` for the list of recorders. Currently just 'DebugRecorder`.
 """
 function abm(parameters, recorder)
-    model_sol = ModelSolution(parameters.sim_length,parameters,5000)
+    model_sol = ModelSolution(parameters.sim_length,parameters,parameters.num_households)
     output = solve!(model_sol,recorder )
 
     # total_weighted_degree = map(modelsol.index_vectors) do age_group_indices

CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl

@@ -33,11 +33,10 @@ Resample all the weights in `mixing_graph`
 """
 function sample_mixing_graph!(mixing_graph)
     mixing_edges = mixing_graph.mixing_edges
-    for i in 1:size(mixing_edges.contact_array)[1], j in 1:i  #diagonal     
+    for i in 1:3, j in 1:i  #diagonal     
         rand!(Random.default_rng(Threads.threadid()), mixing_edges.sampler_matrix[j,i],mixing_edges.sample_cache[j,i])
         for k in 1:length(mixing_edges.contact_array[j,i])
             edge_weight_k = mixing_edges.sample_cache[j,i][k]
-            
             set!(mixing_edges.weights_dict, mixing_edges.contact_array[j,i][k], edge_weight_k)
         end
     end
@@ -87,7 +86,7 @@ This function constructs the `MixingEdges` type for rest and WS graphs. Calls `g
 function create_mixing_edges(demographic_index_vectors,mixing_matrix,weights_distribution_matrix)
     contact_array =[Vector{GraphEdge}() for i in 1:length(demographic_index_vectors),j in 1:length(demographic_index_vectors)]
     tot = 0
-    for i in 1:size(mixing_matrix)[1], j in 1:i  #diagonal
+    for i in 1:3, j in 1:i  #diagonal
         num_degrees_ij = similar(demographic_index_vectors[i])
         num_degrees_ji = similar(demographic_index_vectors[j])
         generate_contact_vectors!(mixing_matrix[i,j],mixing_matrix[j,i],num_degrees_ij,num_degrees_ji)
@@ -200,13 +199,10 @@ end
 Add the edges defined by MixingEdges to the actual graph G. Another big bottleneck since adjancency lists don't add edges super efficiently, and there are a ton of them. 
 """
 function graph_from_mixing_edges(g,mixing_edges)
-    for i in 1:size(mixing_edges.contact_array)[1], j in 1:i  #diagonal
-        for k in 1:length(mixing_edges.contact_array[i,j])
-            e = mixing_edges.contact_array[i,j][k]
-            success = add_edge!(g,e.a,e.b)
-            # if !success
-            #     display((has_edge(g,e.a,e.b)))
-            # end
+    for i in 1:3, j in 1:i  #diagonal
+        for k in eachindex(mixing_edges.contact_array[j,i])
+            e = mixing_edges.contact_array[j,i][k]
+            add_edge!(g,e.a,e.b)
         end
     end
 end

CovidAlertVaccinationModel/src/ABM/model_setup.jl

 function get_parameters()#(0.0000,0.00048,0.0005,0.16,-1.30,-1.24,-0.8,0.35,0.35,0.35,0.2)
     params = (
-        sim_length = 1300,
+        sim_length = 500,
         num_households = 5000,
         I_0_fraction = 0.005,
-        β_y =0.00051,
-        β_m = 0.00047,
-        β_o = 0.07,
+        β_y = 0.00042,
+        β_m = 0.00042,
+        β_o = 0.1,
         α_y = 0.35,
         α_m = 0.35,
         α_o = 0.35,
         recovery_rate = 1/5,
-        π_base_y = -1.35,
-        π_base_m = -1.35,
-        π_base_o = -0.9,
+        π_base_y = -1.0,
+        π_base_m = -1.0,
+        π_base_o = -1.0,
         η = 0.0,
         κ = 0.0,
         ω = 0.00,
@@ -25,14 +25,14 @@ function get_parameters()#(0.0000,0.00048,0.0005,0.16,-1.30,-1.24,-0.8,0.35,0.35
         notification_threshold = 20,
         immunizing = true,
         immunization_delay = 14,
-        immunization_begin_day = 1060, 
-        infection_introduction_day = 1180,
-        O_distribution_shift = 0.5,
+        immunization_begin_day = 60, 
+        infection_introduction_day = 180,
+        O_distribution_shift = 0.0,
     )
     return params
 end
 
-function get_u_0(nodes,I_0_fraction,vaccinator_prob)
+function get_u_0(nodes,vaccinator_prob)
     is_vaccinator = rand(Random.default_rng(Threads.threadid()),nodes) .< vaccinator_prob
     status = fill(Susceptible,nodes)
     return status,is_vaccinator
@@ -105,7 +105,7 @@ mutable struct ModelSolution{T,InfNet,SocNet,WSMixingDist,RestMixingDist}
         app_user_index = findall(==(true),is_app_user)
         app_user_list[is_app_user] .= collect(1:length(app_user_index))
         
-        u_0_inf,u_0_vac = get_u_0(nodes,params.I_0_fraction,params.vaccinator_prob)
+        u_0_inf,u_0_vac = get_u_0(nodes,params.vaccinator_prob)
 
         infected_mixing_graph,soc_mixing_graph = time_dep_mixing_graphs(sim_length,base_network,demographics,index_vectors,ws_matrix_tuple,rest_matrix_tuple)
       
@@ -114,6 +114,12 @@ mutable struct ModelSolution{T,InfNet,SocNet,WSMixingDist,RestMixingDist}
 
         status_totals = [count(==(AgentStatus(i)), u_0_inf) for i in 1:AgentStatus.size]
         immunization_countdown = fill(-1, nodes) #immunization countdown is negative if not counting down
+
+
+        for network in infected_mixing_graph.graph_list[begin] #this also resamples the soc network weights since they point to the same objects, but those are never used
+            sample_mixing_graph!(network) #get new contact weights
+        end
+        
         return new{T,typeof(infected_mixing_graph),typeof(soc_mixing_graph),typeof(ws_matrix_tuple),typeof(rest_matrix_tuple)}(
             sim_length,
             nodes,

CovidAlertVaccinationModel/src/ABM/parameter_optimization.jl

@@ -15,7 +15,7 @@ function solve_w_parameters(default_p, p_names, new_p_list)
 end
 
 function fit_distribution_parameters(p_tuple)
-    p_names = (:ω,:β_y,:β_m,:β_o,:π_base_y,:π_base_m,:π_base_o,:α_y,:α_m,:α_o,:O_distribution_shift)
+    p_names = (:ω,:β_y,:β_m,:β_o,:π_base_y,:π_base_m,:π_base_o,:α_y,:α_m,:α_o)
     priors = Factored(
         Uniform(0.0,0.01),
         Uniform(0.0,0.5),  
@@ -27,7 +27,6 @@ function fit_distribution_parameters(p_tuple)
         TriangularDist(23.0, 44.0, 34.0),
         TriangularDist(22.0, 54.0, 40.0),
         TriangularDist(9.0, 59.0, 39.0),   
-        Uniform(0.0,1.0),   
     )
     #simulation begins in july
     #60 days for opinion dynamics to stabilize, then immunization begins in september,
@@ -54,18 +53,18 @@ function fit_distribution_parameters(p_tuple)
         target_preinf_vac = ymo_vac .* sum(vaccination_data[1:4]) .* length.(model.index_vectors)
         target_postinf_vac = ymo_vac .* sum(vaccination_data[5:end]) .* length.(model.index_vectors)
     
-        display((final_size,target_final_size))
-        display((total_preinf_vaccination,target_preinf_vac))
-        display((total_postinf_vaccination,target_postinf_vac))
-        display(sum.(eachrow(ymo_vaccination_ts)) ./length.(model.index_vectors))
+        # display((final_size,target_final_size))
+        # display((total_preinf_vaccination,target_preinf_vac))
+        # display((total_postinf_vaccination,target_postinf_vac))
+        # display(sum.(eachrow(ymo_vaccination_ts)) ./length.(model.index_vectors))
 
-        return sum((total_preinf_vaccination .- target_preinf_vac).^2)  
-        + sum((total_postinf_vaccination .- total_postinf_vaccination).^2)  
-        + 5*sum((final_size .- target_final_size).^2)
+        # return sum((total_preinf_vaccination .- target_preinf_vac).^2)  
+        # + sum((total_postinf_vaccination .- total_postinf_vaccination).^2)  
+        return sum((final_size .- target_final_size).^2)
     end 
-       display(cost((0.0000,0.00051,0.00046,0.18,-1.30,-1.24,-0.8,0.35,0.35,0.35,0.2)))
+    #    display(cost((0.0000,0.00051,0.00046,0.18,-1.30,-1.24,-0.8,0.35,0.35,0.35,0.2)))
     #
-    # out = smc(priors,cost; verbose = true, nparticles = 600, parallel = true)
+    out = smc(priors,cost; verbose = true, nparticles = 600, parallel = true)
     return NamedTuple{p_names}(ntuple(i -> out.P[i].particles,length(p_names)))
 end
 function fit_parameters(default_parameters)

CovidAlertVaccinationModel/src/ABM/solve.jl

@@ -8,7 +8,7 @@ function contact_weight(β, contact_time)
 end
 
 function Φ(payoff,ξ)
-    return 1 / (exp(-1*ξ*payoff))
+    return 1 / (1 + exp(-1*ξ*payoff))
 end
 
 Base.@propagate_inbounds @views function update_alert_durations!(t,modelsol) # Base.@propagate_inbounds 
@@ -103,33 +103,37 @@ Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,mod
     app_user_pointer = 0
     for i in 1:nodes
         π_base = @SVector [π_base_y,π_base_m,π_base_o]
-        vac_payoff = 0
-        num_soc_nbrs = 0
         random_soc_network = sample(Random.default_rng(Threads.threadid()), soc_network.graph_list[t])
 
         if !isempty(neighbors(random_soc_network,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[Int(demographics[i])] + 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]])) * (η + total_infections*ω_en) 
-                end
-        
-                if u_vac[i]
-                    if rand(Random.default_rng(Threads.threadid())) < 1 - Φ(vac_payoff,ξ)
-                        u_next_vac[i] = false
+            for _ = 1:10
+                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[Int(demographics[i])] + 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]])) * (η + total_infections*ω_en) 
                     end
-                else
-                    if rand(Random.default_rng(Threads.threadid())) < Φ(vac_payoff,ξ)
-                        u_next_vac[i] = true
+            
+                    if u_vac[i]
+                        # display(1 - Φ(vac_payoff,ξ))
+                        if rand(Random.default_rng(Threads.threadid())) < 1 - Φ(vac_payoff,ξ)
+                            u_next_vac[i] = false
+                        end
+                    else
+                        # display( Φ(vac_payoff,ξ))
+                        if rand(Random.default_rng(Threads.threadid())) < Φ(vac_payoff,ξ)
+                            u_next_vac[i] = true
+                        end
                     end
+                    break
                 end
+
             end
         end
     end
 
-    modelsol.daily_vaccinators = count(==(true),u_vac) #could maybe make this more efficient
-
+    modelsol.daily_vaccinators = count(==(true),u_vac) 
 end
 
 function weighted_degree(node,network::TimeDepMixingGraph)
@@ -146,10 +150,6 @@ end
 
 function agents_step!(t,modelsol)
 
-    remake!(modelsol.inf_network,modelsol.index_vectors)
-    for network in modelsol.inf_network.graph_list[t] #this also resamples the soc network weights since they point to the same objects, but those are never used
-        sample_mixing_graph!(network) #get new contact weights
-    end
     if t == modelsol.params.infection_introduction_day
         init_indices = rand(Random.default_rng(Threads.threadid()), 1:modelsol.nodes, round(Int,modelsol.nodes*modelsol.params.I_0_fraction))
         modelsol.u_inf[init_indices] .= Infected
@@ -163,6 +163,11 @@ function agents_step!(t,modelsol)
 
     modelsol.u_vac .= modelsol.u_next_vac
     modelsol.u_inf .= modelsol.u_next_inf
+
+    remake!(modelsol.inf_network,modelsol.index_vectors)
+    for network in modelsol.inf_network.graph_list[t] #this also resamples the soc network weights since they point to the same objects, but those are never used
+        sample_mixing_graph!(network) #get new contact weights
+    end
 end