checking correctness

CovidAlertVaccinationModel/src/ABM/abm.jl

@@ -12,7 +12,7 @@ Runs the model with default parameters.
 """
 function bench()
     p = get_parameters()
-    out = abm(p)
+    out = abm(p;record_degrees = false)
     return out
 end
 

CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl

@@ -33,9 +33,9 @@ Resample all the weights in `mixing_graph`
 """
 function sample_mixing_edges!(wg,sampler_matrix,demographics)
     for k in 1:ne(wg)
-        i = wg.srcs[k]
-        j = wg.sinks[k]
-        wg.weights[k].x = rand(Random.default_rng(Threads.threadid()),sampler_matrix[Int(demographics[j]),Int(demographics[i])])
+        i = wg.edges[k][1]
+        j = wg.edges[k][2]
+        wg.edges[k][3].x = rand(Random.default_rng(Threads.threadid()),sampler_matrix[Int(demographics[j]),Int(demographics[i])])
     end
 end
 
@@ -127,10 +127,9 @@ function remake_all!(t,time_dep_mixing_graph,index_vectors,demographics)
 end
 
 struct Neighborhood
-    neighbourhood::Vector{Int}
-    weights::Vector{Base.RefValue{UInt8}}
+    neighbourhood::Vector{Tuple{Int,Base.RefValue{UInt8}}}
     function Neighborhood()
-        return new(Vector{Int}(),Vector{Base.RefValue{UInt8}}())
+        return new(Vector{Tuple{Int,Base.RefValue{UInt8}}}())
     end
 end
 
@@ -157,19 +156,15 @@ Matrix of distributions determining the edge weights
 
 """
 struct WeightedGraph{M1,M2} 
-    srcs::Vector{Int}
-    sinks::Vector{Int}
-    weights::Vector{Base.RefValue{UInt8}}
+    edges::Vector{Tuple{Int,Int,Base.RefValue{UInt8}}}
     fadjlist::Vector{Neighborhood}
     mixing_matrix::M1
     sampler_matrix::M2
     function WeightedGraph(demographics::AbstractVector,index_vectors,mixing_matrix::M1, sampler_matrix::M2) where {M1,M2}
-        srcs = Vector{Int}()
-        sinks = Vector{Int}()
-        weights = Vector{Base.RefValue{UInt8}}()
+        edges = Vector{Tuple{Int,Int,Base.RefValue{UInt8}}}()
         fadjlist = [Neighborhood() for _ in 1:length(demographics)]
         wg = new{M1,M2}(
-            srcs,sinks,weights,fadjlist,
+            edges,fadjlist,
             mixing_matrix,
             sampler_matrix,
         )
@@ -177,12 +172,10 @@ struct WeightedGraph{M1,M2}
         return wg
     end
     function WeightedGraph(nv,sampler_matrix::M2) where M2
-        srcs = Vector{Int}()
-        sinks = Vector{Int}()
-        weights = Vector{Base.RefValue{UInt8}}()
+        edges = Vector{Tuple{Int,Int,Base.RefValue{UInt8}}}()
         fadjlist = [Neighborhood() for _ in 1:nv]
         wg = new{Nothing,M2}(
-            srcs,sinks,weights,fadjlist,
+            edges,fadjlist,
             nothing,
             sampler_matrix,
         )
@@ -193,14 +186,11 @@ function empty!(wg::WeightedGraph)
     for nbhd in wg.fadjlist    
         empty!(nbhd)
     end
-    empty!(wg.srcs)
-    empty!(wg.sinks)
-    empty!(wg.weights)
+    empty!(wg.edges)
 end
 
 function empty!(nbhd::Neighborhood)
     empty!(nbhd.neighbourhood) #empty all the vector edgelists
-    empty!(nbhd.weights)
 end
 
 function assemble_graph!(g,index_vectors,mixing_matrix,sampler_matrix)
@@ -244,7 +234,6 @@ function fast_chung_lu(pop_i,mixing_dist)
     return stubs_i,stubs_j
 end
 
-neighbors(g::WeightedGraph,i) = g.fadjlist[i].neighbourhood
 function Base.show(io::IO, g::WeightedGraph) 
     print(io, "WG $(ne(g)))")
 end 
@@ -258,21 +247,17 @@ function add_edge!(g::WeightedGraph{T}, a,b,weight) where {T}
     (b in a_nbhd.neighbourhood) && return false  # edge already in graph
     weight_boxed = Ref(weight)
 
-    push!(a_nbhd.neighbourhood, b)
-    push!(a_nbhd.weights,weight_boxed)
-    push!(g.srcs,a)
-    push!(g.sinks,b)
-    push!(g.weights,weight_boxed)
+    push!(a_nbhd.neighbourhood, ( b, weight_boxed))
+    push!(g.edges,(a,b,weight_boxed))
 
     a == b && return true  # selfloop
     
     b_nbhd = g.fadjlist[b]
-    push!(b_nbhd.neighbourhood, a)
-    push!(b_nbhd.weights,weight_boxed)
+    push!(b_nbhd.neighbourhood,( a, weight_boxed))
     
     return true  # edge successfully added
 end
 is_directed(g::WeightedGraph) = false
-ne(g::WeightedGraph) = length(g.srcs)
+ne(g::WeightedGraph) = length(g.edges)
 vertices(g::WeightedGraph) = 1:length(g.fadjlist)
-neighbors_and_weights(g::WeightedGraph,i) = zip(g.fadjlist[i].neighbourhood,g.fadjlist[i].weights)
\ No newline at end of file
+neighbors_and_weights(g::WeightedGraph,i) = g.fadjlist[i].neighbourhood
\ No newline at end of file

CovidAlertVaccinationModel/src/ABM/solve.jl

@@ -110,8 +110,8 @@ Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,mod
         (π_base_y,π_base_m,π_base_o) : 
         (π_base_y*ζ,π_base_m*ζ,π_base_o*ζ)
         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,i))
+        if !isempty(neighbors_and_weights(random_soc_network,i))
+                random_neighbour,_ = sample(Random.default_rng(Threads.threadid()), neighbors_and_weights(random_soc_network,i))
                 app_vac_payoff = 0.0    
                 if app_user[i] && time_of_last_alert[app_user_list[i]]>=0
                     output_data.daily_total_notified_agents[t] += 1
@@ -140,7 +140,7 @@ end
 function weighted_degree(t,node,network::TimeDepMixingGraph)
     weighted_degree = 0
     for g in network.graph_list[t]
-        for (j,w) in neighbors_and_weights(mixing_graph,node)
+        for (j,w) in neighbors_and_weights(g,node)
             weighted_degree += w[]
         end
     end
@@ -192,6 +192,7 @@ function solve!(modelsol)
 
         modelsol.u_vac .= modelsol.u_next_vac
         modelsol.u_inf .= modelsol.u_next_inf
+        display(mean.(modelsol.output_data.avg_weighted_degree))
     end
 end
 

CovidAlertVaccinationModel/test/ABM/mixing_test.jl

 
 using CovidAlertVaccinationModel:abm,generate_contact_vectors!
-using CovidAlertVaccinationModel:ModelSolution,AgentDemographic,mean,get_weight
-using CovidAlertVaccinationModel:AgentStatus,get_u_0,get_parameters,solve!, WeightedGraph, contact_time_distributions,GraphEdge
+using CovidAlertVaccinationModel:ModelSolution,AgentDemographic,mean,neighbors_and_weights, vertices
+using CovidAlertVaccinationModel:AgentStatus,get_u_0,get_parameters,solve!, WeightedGraph, contact_time_distributions
 using LightGraphs
 using Random
 using Plots
@@ -19,7 +19,7 @@ const model_sizes = [1000,5000]
         m = ModelSolution(100,get_parameters(),sz)
         ws_dist = m.ws_matrix_tuple
         r_dist = m.rest_matrix_tuple 
-        index_vec =m.index_vectors
+        index_vec = m.index_vectors
         @testset "workschool" for i in dem_cat, j in dem_cat
             for t in 1:length(ws_dist)
                 @test mean(ws_dist[t][i,j])*length(index_vec[i]) == mean(ws_dist[t][j,i])*length(index_vec[j])
@@ -33,54 +33,7 @@ const model_sizes = [1000,5000]
     end
 end
 
-function approximate_mixing_matricies_weights_dict(p,sol)
-    mean_mixing_degree = zeros(3,3)
-    mean_mixing_weighted_degree = zeros(3,3)
-        for g_list in sol.inf_network.graph_list
-            for (k,g) in enumerate(g_list)
-                for e in keys(g.weights_dict)
-                    demo_i = sol.demographics[e.a]
-                    demo_j = sol.demographics[e.b]
-                    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))
-                end 
-        end
-    end
-    mean_mixing_degree ./= (p.sim_length .* length.(sol.index_vectors) * 2)
-    mean_mixing_weighted_degree ./= (p.sim_length .* length.(sol.index_vectors) * 2)
-    return mean_mixing_degree, mean_mixing_weighted_degree
-end
-
-function approximate_mixing_matricies_graph(p,sol)
-    mean_mixing_degree = zeros(3,3)
-    mean_mixing_weighted_degree = zeros(3,3)
-        for g_list in sol.inf_network.graph_list
-            for g in g_list
-                for e in edges(g.g)  
-                        demo_i = sol.demographics[src(e)]
-                        demo_j = sol.demographics[dst(e)]
-                        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_weighted_degree[Int(demo_i), Int(demo_j)] += get_weight(g,GraphEdge(src(e),dst(e)))
-                        mean_mixing_weighted_degree[Int(demo_j), Int(demo_i)] += get_weight(g,GraphEdge(src(e),dst(e)))
-                end 
-            end
-        end
-        mean_mixing_degree ./= (p.sim_length .* length.(sol.index_vectors) * 2)
-        mean_mixing_weighted_degree ./= (p.sim_length .* length.(sol.index_vectors) * 2)
-    return mean_mixing_degree, mean_mixing_weighted_degree
-end
-@testset "weights dict and graph match" begin
-    p = get_parameters()
-    sol = abm(p)
-    graph_deg_matrix,graph_wdeg_matrix = approximate_mixing_matricies_graph(p,sol)
-    dict_deg_matrix,dict_wdeg_matrix = approximate_mixing_matricies_weights_dict(p,sol)
 
-    @test all(graph_deg_matrix .≈ dict_deg_matrix)
-    @test all(graph_wdeg_matrix .≈ dict_wdeg_matrix)
-end
 
 using OnlineStats
 
@@ -131,10 +84,10 @@ end
     for (expected_dist_mean,dist) in zip(dist_means,(ws_matrix_tuple...,rest_matrix_tuple...))
         g = WeightedGraph(demographics,index_vectors,dist,contact_time_distributions.ws)
         mixing_dist = [Variance() for _ in 1:3, _ in 1:3]
-        for v in vertices(g.g)  
+        for v in vertices(g)  
             demo_v = demographics[v]
             degs = zeros(3)
-            for w in LightGraphs.neighbors(g.g,v)
+            for (w,_) in neighbors_and_weights(g,v)
                 demo_w = demographics[w]
                 degs[Int(demo_w)] += 1
             end