bug squashed I think!

CovidAlertVaccinationModel/src/ABM/abm.jl

@@ -12,7 +12,7 @@ Runs the model with default parameters.
 """
 function bench()
     p = get_parameters()
-    out = abm(p;record_degrees = false)
+    out = abm(p)
     return out
 end
 
@@ -36,3 +36,15 @@ function abm(parameters; record_degrees = false)
     solve!(model_sol)
     return model_sol
 end
+# 3-element Vector{Float64}:
+#  248.40732146371218
+#  302.0494462054946
+#   79.01415331590017
+# 3-element Vector{Float64}:
+#  248.49975738800097
+#  302.13518249351813
+#   79.01634127655387
+# 3-element Vector{Float64}:
+#  248.59457168494825
+#  302.2272494288624
+#   79.01827786039279
\ No newline at end of file

CovidAlertVaccinationModel/src/ABM/agents.jl

@@ -33,6 +33,8 @@ Also returns a vector of AgentDemographic representing each agent, defined by th
     total_household_pop = sum(sum.(households_composition))
     population_list = Vector{AgentDemographic}(undef,total_household_pop)
     wg = WeightedGraph(total_household_pop,sampler_matrix)
+    # test_graph_1 = SimpleGraph(total_household_pop)
+    # test_graph_2 = SimpleGraph(total_household_pop)
     vertex_pointer = 1
     for household in households_composition
         num_vertices = sum(household)
@@ -45,10 +47,18 @@ Also returns a vector of AgentDemographic representing each agent, defined by th
             if v != w
                 weight::UInt8 = rand(Random.default_rng(Threads.threadid()),sampler_matrix[Int(population_list[v]),Int(population_list[w])])
                 add_edge!(wg,v,w,weight)
+                # add_edge!(test_graph_1,v,w)
             end
         end
+        # for v in vertex_pointer:(vertex_pointer + num_vertices - 1),  w in vertex_pointer:(vertex_pointer + num_vertices - 1)
+        #     if v != w
+        #         weight::UInt8 = rand(Random.default_rng(Threads.threadid()),sampler_matrix[Int(population_list[v]),Int(population_list[w])])
+        #         add_edge!(test_graph_2,v,w)
+        #     end
+        # end
         vertex_pointer+=num_vertices
     end
+    # display((test_graph.ne,ne(wg)))
     return wg,population_list
 end
 

CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl

-using Dictionaries
 """
 A type that stores a pair of nodes representing an edge in the graph. 
 
 We need to define a custom type for these so we can define a hash function on graph edges, in order to more efficiently use them in hashmaps (dictionaries)
 """
-struct GraphEdge
+struct HalfWeightedEdge
     a::Int
-    b::Int
+    weight::Base.RefValue{UInt8}
 end
 
-"""
-Define a hash on GraphEdge such that `hash(a,b) = hash(b,a)` (hash is commutative).
-
-This is helpful because then we only need to store (a,b) in the graph edges weights dictionary, rather than both (a,b) and (b,a).
-"""
-function Base.hash(e::GraphEdge)
-    return hash(minmax(e.a,e.b))
-end
 
 """
 Define symmetric edge equality, matches the hash function.
 """
-function Base.isequal(e1::GraphEdge,e2::GraphEdge)
-    return isequal(minmax(e1.a,e1.b),minmax(e2.a,e2.b))
+function Base.:(==)(e1::HalfWeightedEdge,e2::HalfWeightedEdge)
+    return isequal(e1.a,e2.a)
 end
 
 """
@@ -126,13 +117,6 @@ function remake_all!(t,time_dep_mixing_graph,index_vectors,demographics)
     # display_degree(time_dep_mixing_graph.resampled_graphs[1])
 end
 
-struct Neighborhood
-    neighbourhood::Vector{Tuple{Int,Base.RefValue{UInt8}}}
-    function Neighborhood()
-        return new(Vector{Tuple{Int,Base.RefValue{UInt8}}}())
-    end
-end
-
 """
 Weighted graph type. Stores the graph in `g`, and the weights and edges in `mixing_edges`. 
 Fields
@@ -157,12 +141,12 @@ Matrix of distributions determining the edge weights
 """
 struct WeightedGraph{M1,M2} 
     edges::Vector{Tuple{Int,Int,Base.RefValue{UInt8}}}
-    fadjlist::Vector{Neighborhood}
+    fadjlist::Vector{Vector{HalfWeightedEdge}}
     mixing_matrix::M1
     sampler_matrix::M2
     function WeightedGraph(demographics::AbstractVector,index_vectors,mixing_matrix::M1, sampler_matrix::M2) where {M1,M2}
         edges = Vector{Tuple{Int,Int,Base.RefValue{UInt8}}}()
-        fadjlist = [Neighborhood() for _ in 1:length(demographics)]
+        fadjlist = [Vector{HalfWeightedEdge}() for _ in 1:length(demographics)]
         wg = new{M1,M2}(
             edges,fadjlist,
             mixing_matrix,
@@ -173,7 +157,7 @@ struct WeightedGraph{M1,M2}
     end
     function WeightedGraph(nv,sampler_matrix::M2) where M2
         edges = Vector{Tuple{Int,Int,Base.RefValue{UInt8}}}()
-        fadjlist = [Neighborhood() for _ in 1:nv]
+        fadjlist = [Vector{HalfWeightedEdge}() for _ in 1:nv]
         wg = new{Nothing,M2}(
             edges,fadjlist,
             nothing,
@@ -189,10 +173,6 @@ function empty!(wg::WeightedGraph)
     empty!(wg.edges)
 end
 
-function empty!(nbhd::Neighborhood)
-    empty!(nbhd.neighbourhood) #empty all the vector edgelists
-end
-
 function assemble_graph!(g,index_vectors,mixing_matrix,sampler_matrix)
     for i in 1:3, j in 1:i #diagonal
         if i != j 
@@ -244,20 +224,23 @@ function add_edge!(g::WeightedGraph{T}, a,b,weight) where {T}
     verts = vertices(g)
     (a in verts && b in verts) || return false  # edge out of bounds
     a_nbhd = g.fadjlist[a]
-    (b in a_nbhd.neighbourhood) && return false  # edge already in graph
     weight_boxed = Ref(weight)
+    b_edge = HalfWeightedEdge(b, weight_boxed)
+
+    (b_edge in a_nbhd) && return false  # edge already in graph
+
 
-    push!(a_nbhd.neighbourhood, ( b, weight_boxed))
+    push!(a_nbhd, b_edge)
     push!(g.edges,(a,b,weight_boxed))
 
     a == b && return true  # selfloop
     
     b_nbhd = g.fadjlist[b]
-    push!(b_nbhd.neighbourhood,( a, weight_boxed))
+    push!(b_nbhd,HalfWeightedEdge( a, weight_boxed))
     
     return true  # edge successfully added
 end
 is_directed(g::WeightedGraph) = false
 ne(g::WeightedGraph) = length(g.edges)
 vertices(g::WeightedGraph) = 1:length(g.fadjlist)
-neighbors_and_weights(g::WeightedGraph,i) = g.fadjlist[i].neighbourhood
\ No newline at end of file
+neighbors_and_weights(g::WeightedGraph,i) = g.fadjlist[i]
\ No newline at end of file

CovidAlertVaccinationModel/src/ABM/solve.jl

@@ -21,9 +21,9 @@ Base.@propagate_inbounds @views function update_alert_durations!(t,modelsol) # B
         end
         total_weight_i = 0
         for mixing_graph in inf_network.graph_list[t]
-            for (j,w) in neighbors_and_weights(mixing_graph,node)
+            for he in neighbors_and_weights(mixing_graph,node)
                 if app_user[j]
-                    total_weight_i += w[]
+                    total_weight_i += he.weight[]
                 end
             end
         end
@@ -61,11 +61,12 @@ Base.@propagate_inbounds @views function update_infection_state!(t,modelsol; rec
         end
         if  agent_status == Susceptible || agent_status == Immunized 
             for mixing_graph in inf_network.graph_list[t]
-                for (j,w) in neighbors_and_weights(mixing_graph,i)   
+                for he in neighbors_and_weights(mixing_graph,i)
+                    j = he.a   
                     if u_inf[j] == Infected && u_next_inf[i] != Infected
                         β_vec = (β_y,β_m,β_o)
                         α_vec = (α_y,α_m,α_o)
-                        if rand(Random.default_rng(Threads.threadid())) < contact_weight(β_vec[Int(agent_demo)],w[])
+                        if rand(Random.default_rng(Threads.threadid())) < contact_weight(β_vec[Int(agent_demo)],he.weight[])
                             if agent_status == Immunized && rand(Random.default_rng(Threads.threadid())) < 1- α_vec[Int(agent_demo)]
                                 agent_transition!(i, Immunized,Infected)
                                 output_data.daily_cases_by_age[Int(agent_demo),t]+=1
@@ -111,7 +112,7 @@ Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,mod
         (π_base_y*ζ,π_base_m*ζ,π_base_o*ζ)
         random_soc_network = sample(Random.default_rng(Threads.threadid()), soc_network.graph_list[t])
         if !isempty(neighbors_and_weights(random_soc_network,i))
-                random_neighbour,_ = sample(Random.default_rng(Threads.threadid()), neighbors_and_weights(random_soc_network,i))
+                random_neighbour = sample(Random.default_rng(Threads.threadid()), neighbors_and_weights(random_soc_network,i)).a
                 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,8 +141,8 @@ 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(g,node)
-            weighted_degree += w[]
+        for he in neighbors_and_weights(g,node)
+            weighted_degree += he.weight[]
         end
     end
     return weighted_degree
@@ -151,8 +152,8 @@ function sample_initial_nodes(nodes,graphs,I_0_fraction)
     weighted_degrees = zeros(nodes)
     for v in 1:nodes
         for g in graphs
-            for (j,w) in neighbors_and_weights(g,v)
-                weighted_degrees[v] += w[]
+            for he in neighbors_and_weights(g,v)
+                weighted_degrees[v] += he.weight[]
             end
         end
     end
@@ -192,7 +193,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))
+        # display(mean.(modelsol.output_data.avg_weighted_degree))
     end
 end