initial weighted degree implementation

CovidAlertVaccinationModel/src/ABM/output.jl

@@ -24,10 +24,10 @@ struct Recorder{ElType,ArrT1,ArrT2,ArrT3,StatAccumulator}
     avg_weighted_degree_of_vaccinators_no_EN::Vector{StatAccumulator}
     total_app_weight::ArrT2
     avg_weighted_degree_of_vaccinator_belief::Vector{StatAccumulator}
-    vac_join_counts::ArrT2
-    vac_join_11_counts::ArrT2
-    vac_join_01_counts::ArrT2
-    vac_join_00_counts::ArrT2
+    assortativity_at_end::StatAccumulator
+    vac_join_11_counts_at_end::StatAccumulator
+    vac_join_01_counts_at_end::StatAccumulator
+    vac_join_00_counts_at_end::StatAccumulator
     
 
     function Recorder(val::T,sim_length; record_degrees = false) where T
@@ -51,10 +51,6 @@ struct Recorder{ElType,ArrT1,ArrT2,ArrT3,StatAccumulator}
         final_size_by_age =  [copy(val) for i in 1:3]
         total_app_weight = [copy(val) for i in 1:sim_length]
 
-        vac_join_counts =[copy(val) for j in 1:sim_length]
-        vac_join_11_counts =[copy(val) for j in 1:sim_length]
-        vac_join_01_counts =[copy(val) for j in 1:sim_length]
-        vac_join_00_counts =[copy(val) for j in 1:sim_length]
 
         avg_weighted_degree_of_vaccinators = [Variance() for _ in 1:3]
         avg_weighted_degree_of_vaccinators_no_EN = [Variance() for _ in 1:3]
@@ -78,10 +74,10 @@ struct Recorder{ElType,ArrT1,ArrT2,ArrT3,StatAccumulator}
             avg_weighted_degree_of_vaccinators_no_EN,
             total_app_weight,
             avg_weighted_degree_of_vaccinator_belief,
-            vac_join_counts,
-            vac_join_11_counts,
-            vac_join_01_counts,
-            vac_join_00_counts
+            Variance(),
+            Variance(),
+            Variance(),
+            Variance()
         )
     end
 end

CovidAlertVaccinationModel/src/ABM/plotting.jl

@@ -33,8 +33,8 @@ function plot_model(varname,univariate_series, output_list::Vector{T},p) where T
         zoom((x_axis, data.daily_immunized_by_age.Y, "new Y vaccinations each day (zoomed in)", "no. individuals")),
         zoom((x_axis, data.daily_immunized_by_age.M, "new M vaccinations each day (zoomed in)", "no. individuals")),
         zoom((x_axis, data.daily_immunized_by_age.O, "new O vaccinations each day (zoomed in)", "no. individuals")),
-        (x_axis, data.vac_join_counts, "join counts", "no. joins"),
-        (x_axis, data.vac_join_01_counts, "join 01 counts", "no. joins"),
+        # (x_axis, data.vac_join_counts, "join counts", "no. joins"),
+        # (x_axis, data.vac_join_01_counts, "join 01 counts", "no. joins"),
     ]
     sample_data = ts_list(output_list[1])
     l = length(sample_data)

CovidAlertVaccinationModel/src/ABM/solve.jl

@@ -80,7 +80,7 @@ function vaccinate_agent!(t,agent,immunization_countdown)
     end
 end
 
-function update_vaccination_opinion_state!(t,agent,agent_demo,modelsol,vaccinate_today_flag,total_infections,π_base_vec,weighted_degree_of_i)
+function update_vaccination_opinion_state!(t,agent,agent_demo,modelsol,vaccinate_today_flag,total_infections,π_base_vec)
     @unpack infection_introduction_day,immunization_intervals, η,Γ,ζ, ω, ω_en,ξ = modelsol.params
     @unpack immunization_countdown,demographics,time_of_last_alert, soc_network,u_vac,u_next_vac,is_app_user,output_data = modelsol
     random_soc_network = sample(Random.default_rng(Threads.threadid()), soc_network.graph_list[t])
@@ -105,9 +105,9 @@ function update_vaccination_opinion_state!(t,agent,agent_demo,modelsol,vaccinate
 
                     if vaccinate_today_flag && t>infection_introduction_day
                         vaccinate_agent!(t, agent, immunization_countdown)
-                        fit!(output_data.avg_weighted_degree_of_vaccinators[Int(agent_demo)],weighted_degree_of_i)
+                        # fit!(output_data.avg_weighted_degree_of_vaccinators[Int(agent_demo)],weighted_degree_of_i)
                         if !is_app_user[agent]
-                            fit!(output_data.avg_weighted_degree_of_vaccinators_no_EN[Int(agent_demo)],weighted_degree_of_i)
+                            # fit!(output_data.avg_weighted_degree_of_vaccinators_no_EN[Int(agent_demo)],weighted_degree_of_i)
                         end
                     end
                 end
@@ -118,26 +118,47 @@ function update_vaccination_opinion_state!(t,agent,agent_demo,modelsol,vaccinate
     end
     return 0
 end
+using WeightedOnlineStats
+function assortativity(t,network::TimeDepMixingGraph,u_inf)
+    
+    total_weight = 0.0
+    x_stats = WeightedVariance()
+    y_stats = WeightedVariance()
 
-function weighted_degree_and_join_count(t,node,network::TimeDepMixingGraph,u_inf)
-    weighted_degree = 0
-    join_count_11 = 0
-    join_count_01 = 0
-    join_count_00 = 0
-
-    node_is_immunized = u_inf[node] == Immunized
     for g in network.graph_list[t]
-        for j in neighbors(g,node)
-            weight = get_weight(g,GraphEdge(node,j))
-            weighted_degree += weight
+        for e in edges(g.g)
+            v = src(e)
+            w = dst(e)
+            weight = get_weight(g,GraphEdge(v,w))
+            total_weight += weight
+            fit!(x_stats,(u_inf[v] == Immunized),weight)
+            fit!(y_stats,(u_inf[v] != Immunized),weight)
+        end
+    end
+    σ_x = std(x_stats)
+    σ_y = std(y_stats)
+    x_bar = mean(x_stats)
+    y_bar = mean(y_stats)
 
-            j_is_immunized = u_inf[j] == Immunized
-            join_count_11 += ifelse(node_is_immunized && j_is_immunized,weight,0)
-            join_count_01 += ifelse(xor(node_is_immunized, j_is_immunized),weight,0)
-            join_count_00 += ifelse(!(node_is_immunized) && !(j_is_immunized),weight,0)
+    numerator = 0.0
+    join_count_11 = 0.0
+    join_count_01 = 0.0
+    join_count_00 = 0.0
+    for g in network.graph_list[t]
+        for e in edges(g.g)
+            v = src(e)
+            w = dst(e)
+            weight = get_weight(g,GraphEdge(v,w))
+            v_inf = u_inf[v] == Immunized
+            w_inf = u_inf[w] != Immunized
+            numerator += weight*(v_inf - x_bar) * (w_inf - y_bar)
+            #ughhhh
+            join_count_11 += ifelse( u_inf[v] == Immunized &&  u_inf[w] == Immunized,weight,0)
+            join_count_01 += ifelse(xor( u_inf[v] == Immunized, u_inf[w] == Immunized),weight,0)
+            join_count_00 += ifelse(!(u_inf[v] == Immunized) && !(u_inf[w] == Immunized),weight,0)
         end
     end
-    return weighted_degree,(join_count_11,join_count_01,join_count_00)
+    return numerator/(total_weight*σ_x*σ_y),(join_count_11,join_count_01,join_count_00)
 end
 
 function sample_initial_nodes(nodes,graphs,I_0_fraction)
@@ -158,7 +179,7 @@ end
 
 function agents_step!(t,modelsol)
     @unpack EN_intervals,immunizing,infection_introduction_day,recovery_rate,β_y,β_m,β_o,α_y,α_m,α_o,π_base_y,π_base_m,π_base_o,ζ,immunization_intervals = modelsol.params
-    @unpack u_inf,u_vac,u_next_vac,u_next_inf,demographics,inf_network, is_app_user, immunization_countdown, output_data = modelsol
+    @unpack u_inf,u_vac,u_next_vac,u_next_inf,demographics,inf_network,soc_network, is_app_user, immunization_countdown, output_data = modelsol
 
     vaccinate_today_flag = 
         !immunizing ? false :
@@ -178,20 +199,20 @@ function agents_step!(t,modelsol)
     else
         π_base_vec = @SVector [π_base_y*ζ,π_base_m*ζ,π_base_o*ζ]
     end
-
-    for (agent,(agent_inf_status,agent_is_vaccinator,agent_demo)) in enumerate(zip(u_inf,u_vac,demographics))
-        weighted_degree_of_i, join_counts = output_data.record_degrees_flag ? weighted_degree_and_join_count(t,agent,inf_network,u_inf) : (0,(0,0,0))
-       
-        output_data.vac_join_counts[t] += join_counts[1] + join_counts[3]
-        output_data.vac_join_11_counts[t] += join_counts[1]
-        output_data.vac_join_01_counts[t] += join_counts[2]
-        output_data.vac_join_00_counts[t] += join_counts[3]
+    if t == modelsol.sim_length
+        weighted_assortativity, join_counts = assortativity(t,soc_network,u_inf)
+        fit!(modelsol.output_data.assortativity_at_end,weighted_assortativity)
+        fit!(modelsol.output_data.vac_join_11_counts_at_end,join_counts[1])
+        fit!(modelsol.output_data.vac_join_01_counts_at_end,join_counts[2])
+        fit!(modelsol.output_data.vac_join_00_counts_at_end,join_counts[3])
+    end        
+    for (agent,(agent_inf_status,agent_is_vaccinator,agent_demo)) in enumerate(zip(u_inf,u_vac,demographics))       
         if agent_inf_status == Susceptible 
             is_susceptible_infected!(t, modelsol, agent,agent_inf_status,agent_demo,α_vec,β_vec)
             if vaccinate_today_flag && agent_is_vaccinator && t<(infection_introduction_day-1)
-                fit!(output_data.avg_weighted_degree_of_vaccinators[Int(agent_demo)],weighted_degree_of_i)
+                # fit!(output_data.avg_weighted_degree_of_vaccinators[Int(agent_demo)],weighted_degree_of_i)
                 if !is_app_user[agent]
-                    fit!(output_data.avg_weighted_degree_of_vaccinators_no_EN[Int(agent_demo)],weighted_degree_of_i)
+                    # fit!(output_data.avg_weighted_degree_of_vaccinators_no_EN[Int(agent_demo)],weighted_degree_of_i)
                 end
                 vaccinate_agent!(t, agent, immunization_countdown)
             end
@@ -213,9 +234,9 @@ function agents_step!(t,modelsol)
             immunization_countdown[agent] -= 1
         end
         if agent_is_vaccinator
-            fit!(output_data.avg_weighted_degree_of_vaccinator_belief[t],weighted_degree_of_i)
+            # fit!(output_data.avg_weighted_degree_of_vaccinator_belief[t],weighted_degree_of_i)
         end
-        update_vaccination_opinion_state!(t, agent,agent_demo, modelsol,vaccinate_today_flag, total_infected,π_base_vec,weighted_degree_of_i)
+        update_vaccination_opinion_state!(t, agent,agent_demo, modelsol,vaccinate_today_flag, total_infected,π_base_vec)
     end
     output_data.total_vaccinators[t] = count(==(true),u_vac) 
     u_vac .= u_next_vac