diff --git a/CovidAlertVaccinationModel/src/ABM/abm.jl b/CovidAlertVaccinationModel/src/ABM/abm.jl
index 4fc434d4685a8a5484be2200df50d5be3b75be45..f63f50b8d1f0a724b6bf9b5607c6d5b5eaf0bc0a 100644
--- a/CovidAlertVaccinationModel/src/ABM/abm.jl
+++ b/CovidAlertVaccinationModel/src/ABM/abm.jl
@@ -12,8 +12,7 @@ function bench()
steps = 500
model_sol = ModelSolution(steps,get_parameters(),5000)
recording = DebugRecorder(steps)
- solve!(model_sol,recording )
- display(output.recorded_status_totals[2,:])
+ output = solve!(model_sol,recording )
end
function abm()
diff --git a/CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl b/CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl
index 1414a36d5e00fc958f1829988e6cefd48d8bee1c..7a113d106c8a03a31d4e04e7594d35e1fa3f1034 100644
--- a/CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl
+++ b/CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl
@@ -1,17 +1,39 @@
-#A type that defines a matrix of vectors, such that the sum of the ijth vector is equal to the sum of the jith vector
-struct MixingEdges{V}
+#A type that defines a matrix of vectors, such that the ijth vector contains the nodes to connect to the jith vector, ala. Chung-Lu bipartite
+struct MixingEdges
total_edges::Int
- contact_array::V
+ contact_array::Matrix{Tuple{Vector{Int},Vector{Int}}}
function MixingEdges(demographic_index_vectors,mixing_matrix)
- contacts = map(CartesianIndices(mixing_matrix)) do ind
- zeros(Int,length(demographic_index_vectors[ind[1]]))
- end
+ contact_array =[(Vector{Int}(),Vector{Int}()) for i in 1:length(demographic_index_vectors),j in 1:length(demographic_index_vectors)]
tot = 0
- for i in 1:length(mixing_matrix[:,1]), j in 1:i #diagonal
- generate_contact_vectors!(mixing_matrix[i,j],mixing_matrix[j,i],contacts[i,j],contacts[j,i] )
- tot += sum(contacts[i,j])
+ for i in 1:size(mixing_matrix)[1], j in 1:i #diagonal
+ num_degrees_ij = zeros(Int,length(demographic_index_vectors[i]))
+ num_degrees_ji = zeros(Int,length(demographic_index_vectors[j]))
+ generate_contact_vectors!(mixing_matrix[i,j],mixing_matrix[j,i],num_degrees_ij,num_degrees_ji)
+
+ m = sum(num_degrees_ij)
+ stubs_i = Vector{Int}(undef,m)
+ stubs_j = Vector{Int}(undef,m)
+ if m>0
+ # @assert m == sum(num_degrees_ji)
+ rand!(RNG,DiscreteNonParametric(demographic_index_vectors[i],num_degrees_ij./m),stubs_i)
+ rand!(RNG,DiscreteNonParametric(demographic_index_vectors[j],num_degrees_ji./m),stubs_j)
+ tot += m
+ end
+ contact_array[j,i] = (stubs_i,stubs_j)
end
- new{typeof(contacts)}(tot,contacts)
+
+ new(tot,contact_array)
+ end
+ function MixingEdges(g::SimpleGraph,demographics,demographic_index_vectors)
+ contact_array = [(Vector{Int}(),Vector{Int}()) for i in 1:length(demographic_index_vectors),j in 1:length(demographic_index_vectors)]
+ for e in edges(g)
+ i = src(e)
+ j = dst(e)
+ push!(contact_array[Int(demographics[j]),Int(demographics[i])][1],i)
+ push!(contact_array[Int(demographics[j]),Int(demographics[i])][2],j)
+ end
+ # @show ne(g)
+ new(ne(g),contact_array)
end
end
@@ -25,38 +47,9 @@ struct TimeDepMixingGraph{N,G}
)
end
end
-#modify g so that nodes specified in anodes and bnodes are connected by a bipartite graph with expected degrees given by aseq and bseq
-#implemented from Aksoy, S. G., Kolda, T. G., & Pinar, A. (2017). Measuring and modeling bipartite graphs with community structure
-#simple algorithm, lightgraphs does not allow parallel edges.
-#(also not actually that fast, this is a major bottleneck)
-function random_bipartite_graphs_fast!(g, demographic_index_vectors,mixing_contacts)
- astubs = Vector{Int}(undef,mixing_contacts.total_edges)
- bstubs = Vector{Int}(undef,mixing_contacts.total_edges)
- stubs_ptr = 1
- for i in 1:length(demographic_index_vectors), j in 1:i #diagonal
- anodes = demographic_index_vectors[i]
- bnodes = demographic_index_vectors[j]
- aseq = mixing_contacts.contact_array[i,j]
- bseq = mixing_contacts.contact_array[j,i]
- # @show sum(contacts[i,j])
- m = Int(sum(aseq))
- if m>0
- astubs_ij = @view astubs[stubs_ptr:stubs_ptr + m - 1]
- bstubs_ij = @view bstubs[stubs_ptr:stubs_ptr + m - 1]
- @assert sum(aseq) == sum(bseq) "degree sequences must have equal sum"
- rand!(RNG,DiscreteNonParametric(anodes,aseq./m),astubs_ij)
- rand!(RNG,DiscreteNonParametric(bnodes,bseq./m),bstubs_ij)
- for k in 1:m
- add_edge!(g,astubs_ij[k], bstubs_ij[k])
- end
- stubs_ptr += m
- end
- end
- return astubs,bstubs
-end
function time_dep_mixing_graphs(len,base_network,demographics,index_vectors,ws_matrix_tuple,rest_matrix_tuple)
- home_static_edges = WeightedGraph(base_network,contact_time_distributions.hh) #network with households and LTC homes
+ home_static_edges = WeightedGraph(base_network,demographics,index_vectors,contact_time_distributions.hh) #network with households and LTC homes
ws_static_edges = WeightedGraph(demographics,index_vectors,ws_matrix_tuple.daily,contact_time_distributions.ws)
ws_weekly_edges = WeightedGraph(demographics,index_vectors,ws_matrix_tuple.twice_a_week,contact_time_distributions.ws)
ws_daily_edges = WeightedGraph(demographics,index_vectors,ws_matrix_tuple.otherwise,contact_time_distributions.ws)
@@ -90,63 +83,76 @@ function remake!(time_dep_mixing_graph,demographic_index_vectors,mixing_matrix)
empty!.(weighted_graph.g.fadjlist) #empty all the vector edgelists
empty!(weighted_graph.weights_dict)
weighted_graph.g.ne = 0
- contacts = MixingEdges(demographic_index_vectors,mixing_matrix)
- astubs,bstubs = random_bipartite_graphs_fast!(weighted_graph.g,demographic_index_vectors,contacts)
- weighted_graph.astubs = astubs
- weighted_graph.bstubs = bstubs
+ weighted_graph.mixing_edges = MixingEdges(demographic_index_vectors,mixing_matrix)
+ graph_from_mixing_edges(weighted_graph.g,weighted_graph.mixing_edges)
+ weighted_graph.sample_cache = map(v-> Vector{Int}(undef,length(v[1])),weighted_graph.mixing_edges.contact_array)
+ end
+end
+
+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][1])
+ add_edge!(g,mixing_edges.contact_array[i,j][1][k],mixing_edges.contact_array[i,j][2][k])
+ end
end
end
+
#defining my own weighted graph type cause we need to be able to resample the edge weights in a very particular way
mutable struct WeightedGraph{G, V,M}
g::G
weights_dict::V
weights_distribution_matrix::M
- astubs::Vector{Int}
- bstubs::Vector{Int}
+ mixing_edges::MixingEdges
+ sample_cache::Matrix{Vector{Int}}
function WeightedGraph(demographics,demographic_index_vectors,mixing_matrix,weights_distribution_matrix)
- mixing_contacts = MixingEdges(demographic_index_vectors,mixing_matrix)
+ mixing_edges = MixingEdges(demographic_index_vectors,mixing_matrix)
g = Graph(length(demographics))
+ graph_from_mixing_edges(g,mixing_edges)
weights_dict = RobinDict{Tuple{Int,Int},UInt8}()
sampler_matrix = map(m -> Distributions.PoissonCountSampler.(m),weights_distribution_matrix)
- astubs,bstubs = random_bipartite_graphs_fast!(g,demographic_index_vectors,mixing_contacts)
+ sample_cache = map(v-> Vector{Int}(undef,length(v[1])),mixing_edges.contact_array)
+
return new{typeof(g),typeof(weights_dict),typeof(sampler_matrix)}(
g,
weights_dict,
sampler_matrix,
- astubs,
- bstubs
+ mixing_edges,
+ sample_cache
)
end
- function WeightedGraph(g::SimpleGraph,weights_distribution_matrix)
- astubs = Vector{Int}(undef,g.ne)
- bstubs = Vector{Int}(undef,g.ne)
- for (k,e) in enumerate(edges(g))
- astubs[k] = src(e)
- bstubs[k] = dst(e)
- end
+ function WeightedGraph(g::SimpleGraph,demographics,demographic_index_vectors,weights_distribution_matrix)
+ mixing_edges = MixingEdges(g, demographics, demographic_index_vectors)
weights_dict = RobinDict{Tuple{Int,Int},UInt8}()
sampler_matrix = map(m -> Distributions.PoissonCountSampler.(m),weights_distribution_matrix)
+ # @show map(t->length.(t), mixing_edges.contact_array)
+ sample_cache = map(v-> Vector{Int}(undef,length(v[1])),mixing_edges.contact_array)
return new{typeof(g),typeof(weights_dict),typeof(sampler_matrix)}(
g,
weights_dict,
sampler_matrix,
- astubs,
- bstubs
+ mixing_edges,
+ sample_cache
)
end
end
function Base.show(io::IO, g::WeightedGraph)
- print(io, "WG")
+ print(io, "WG $(ne(g.g))")
end
function sample_mixing_graph!(mixing_graph,population_demographics)
- for k in eachindex(mixing_graph.astubs)
- i = mixing_graph.astubs[k]
- j = mixing_graph.bstubs[k]
+ @unpack mixing_edges = mixing_graph
+
+ for i in 1:size(mixing_edges.contact_array)[1], j in 1:i #diagonal
demo_i = Int(population_demographics[i])
demo_j = Int(population_demographics[j])
- contact_time = rand(RNG, mixing_graph.weights_distribution_matrix[demo_i,demo_j])
- mixing_graph.weights_dict[(i,j)] = contact_time
- mixing_graph.weights_dict[(j,i)] = contact_time
+ rand!(RNG, mixing_graph.weights_distribution_matrix[demo_i,demo_j],mixing_graph.sample_cache[i,j])
+ for k in 1:length(mixing_edges.contact_array[i,j][1])
+ kth_node_i = mixing_edges.contact_array[i,j][1][k]
+ kth_node_j = mixing_edges.contact_array[i,j][2][k]
+ edge_weight_k = mixing_graph.sample_cache[i,j][k]
+ mixing_graph.weights_dict[(kth_node_i,kth_node_j)] = edge_weight_k
+ mixing_graph.weights_dict[(kth_node_j,kth_node_i)] = edge_weight_k
+ end
+
end
end
@inline function reindex!(k,csum,index_list_i,index_list_j,j_to_i_contacts,i_to_j_contacts,sample_list_i,sample_list_j)
diff --git a/CovidAlertVaccinationModel/src/ABM/output.jl b/CovidAlertVaccinationModel/src/ABM/output.jl
index 7bd34defba666a5fb640810dadfa2b3dfda81214..61748f47892381dd5e7268c1f9135e54bf6d2ef4 100644
--- a/CovidAlertVaccinationModel/src/ABM/output.jl
+++ b/CovidAlertVaccinationModel/src/ABM/output.jl
@@ -12,12 +12,12 @@ struct DebugRecorder <: AbstractRecorder
Total_Vaccinator::Vector{Int}
function DebugRecorder(sim_length)
return new(
- zeros(Int,AgentStatus.size,sim_length),
- zeros(Int,sim_length),
- zeros(Int,sim_length),
- zeros(Int,sim_length),
- zeros(Int,sim_length),
- zeros(Int,sim_length),
+ Matrix{Int}(undef,AgentStatus.size,sim_length),
+ Vector{Int}(undef,sim_length),
+ Vector{Int}(undef,sim_length),
+ Vector{Int}(undef,sim_length),
+ Vector{Int}(undef,sim_length),
+ Vector{Int}(undef,sim_length),
)
end
end
diff --git a/CovidAlertVaccinationModel/src/ABM/solve.jl b/CovidAlertVaccinationModel/src/ABM/solve.jl
index 8926e22ca2272474b9026fb684675c2d81f5c173..efcd7360c3e471a93840fc615772f9686642e264 100644
--- a/CovidAlertVaccinationModel/src/ABM/solve.jl
+++ b/CovidAlertVaccinationModel/src/ABM/solve.jl
@@ -12,6 +12,7 @@ function update_alert_durations!(t,modelsol)
covid_alert_times[i,j-1] = covid_alert_times[i,j] #shift them all back
end
for j in neighbors(mixing_graph.g,node)
+
if app_user[j]
covid_alert_times[i,end] += mixing_graph.weights_dict[(node,j)] #add the contact times for today to the back
end
@@ -26,6 +27,7 @@ function update_infection_state!(t,modelsol)
@unpack u_inf,u_vac,u_next_inf,u_next_vac,demographics,inf_network,status_totals = modelsol
function agent_transition!(node, from::AgentStatus,to::AgentStatus)
+ # @show ((from,to))
status_totals[Int(from)] -= 1
status_totals[Int(to)] += 1
u_next_inf[node] = to
@@ -44,6 +46,8 @@ function update_infection_state!(t,modelsol)
for j in neighbors(mixing_graph.g,i)
if u_inf[j] == Infected && u_next_inf[i] != Infected
if rand(RNG) < contact_weight(base_transmission_probability,mixing_graph.weights_dict[(i,j)])
+
+
agent_transition!(i, Susceptible,Infected)
end
end
@@ -60,6 +64,7 @@ function update_infection_state!(t,modelsol)
end
end
end
+ # display(u_next_inf)
end
function update_vaccination_opinion_state!(t,modelsol,total_infections)
@@ -122,7 +127,6 @@ function agents_step!(t,modelsol)
modelsol.u_vac .= modelsol.u_next_vac
modelsol.u_inf .= modelsol.u_next_inf
- # modelsol.status_totals .= modelsol.status_totals_next
end