diff --git a/CovidAlertVaccinationModel/src/ABM/abm.jl b/CovidAlertVaccinationModel/src/ABM/abm.jl index 084fd4d9b50bec979248cd63be3c53e967553a1a..563b084c61ea848ca058a99101025307595f01e2 100644 --- a/CovidAlertVaccinationModel/src/ABM/abm.jl +++ b/CovidAlertVaccinationModel/src/ABM/abm.jl @@ -10,8 +10,8 @@ default(framestyle = :box) function bench() steps = 100 model_sol = ModelSolution(steps,get_parameters(),5000); - recording = recorder(steps) - output = solve!(model_sol,DebugRecorder) + recording = DebugRecorder(steps) + 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 c038e83329d3754a8e7f239ef4ecd647637ae913..30e0fea605dc521f80b4cd1218efe0882aa9dd67 100644 --- a/CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl +++ b/CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl @@ -31,11 +31,71 @@ function random_bipartite_graph_fast_CL!(g::SimpleGraph,anodes,bnodes,aseq,bseq, return g end +struct TimeDepMixingGraph{N,G} + resampled_graphs::NTuple{N,G} + graph_list::Vector{Vector{G}} + function TimeDepMixingGraph(len,resampled_graphs::NTuple{N,G},base_graph_list::Vector{G}) where {G,N} + return new{N,G}( + resampled_graphs, + [copy(base_graph_list) for i in 1:len] + ) + end +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 + 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) + rest_static_edges = WeightedGraph(demographics,index_vectors,rest_matrix_tuple.daily,contact_time_distributions.rest) + rest_weekly_edges = WeightedGraph(demographics,index_vectors,rest_matrix_tuple.twice_a_week,contact_time_distributions.rest) + rest_daily_edges = WeightedGraph(demographics,index_vectors,rest_matrix_tuple.otherwise,contact_time_distributions.rest) + inf_network_list = [home_static_edges,rest_static_edges] + soc_network_list = [home_static_edges,rest_static_edges,ws_static_edges] + + infected_mixing_graph = TimeDepMixingGraph(len,(ws_daily_edges,rest_daily_edges),inf_network_list) + soc_mixing_graph = TimeDepMixingGraph(len,(ws_daily_edges,rest_daily_edges),soc_network_list) + # display(infected_mixing_graph.graph_list) + for (t,l) in enumerate(infected_mixing_graph.graph_list) + day_of_week = mod(t,7) + if !(day_of_week == 3 || day_of_week == 4) #simulation begins on thursday I guess + push!(l, ws_static_edges) + end + if rand(RNG)<5/7 + push!(l, ws_weekly_edges) + push!(l, rest_weekly_edges) + end + push!(l,ws_daily_edges) + push!(l,rest_daily_edges) + end + + return infected_mixing_graph,soc_mixing_graph +end + +function remake!(time_dep_mixing_graph,demographic_index_vectors,mixing_matrix) + for weighted_graph in time_dep_mixing_graph.resampled_graphs + empty!.(weighted_graph.g.fadjlist) #empty all the vector edgelists + empty!(weighted_graph.weights_dict) + weighted_graph.g.ne = 0 + + contacts = MixingContacts(demographic_index_vectors,mixing_matrix) + + for i in 1:length(demographic_index_vectors), j in 1:i #diagonal + random_bipartite_graph_fast_CL!( + weighted_graph.g, + demographic_index_vectors[i], + demographic_index_vectors[j], + contacts.contact_array[i,j], + contacts.contact_array[j,i], + weighted_graph.weights_dict + ) + end + end +end #defining my own weighted graph type cause we need to be able to resample the edge weights in a particular way struct WeightedGraph{G, V,M} g::G - weights::V + weights_dict::V weights_distribution_matrix::M function WeightedGraph(demographics,demographic_index_vectors,mixing_matrix,weights_distribution_matrix) @@ -43,29 +103,31 @@ struct WeightedGraph{G, V,M} g = Graph(length(demographics)) - weights = RobinDict{Tuple{Int,Int},UInt8}() + weights_dict = RobinDict{Tuple{Int,Int},UInt8}() for i in 1:length(demographic_index_vectors), j in 1:i #diagonal - random_bipartite_graph_fast_CL!(g,demographic_index_vectors[i],demographic_index_vectors[j],contacts.contact_array[i,j],contacts.contact_array[j,i],weights) + random_bipartite_graph_fast_CL!(g,demographic_index_vectors[i],demographic_index_vectors[j],contacts.contact_array[i,j],contacts.contact_array[j,i],weights_dict) end covid_alert_time = zeros(nv(g)) - return new{typeof(g),typeof(weights),typeof(weights_distribution_matrix)}( + return new{typeof(g),typeof(weights_dict),typeof(weights_distribution_matrix)}( g, - weights, + weights_dict, weights_distribution_matrix ) end function WeightedGraph(g::SimpleGraph,weights_distribution_matrix) - weights = RobinDict{Tuple{Int,Int},UInt8}() - return new{typeof(g),typeof(weights),typeof(weights_distribution_matrix)}( + weights_dict = RobinDict{Tuple{Int,Int},UInt8}() + return new{typeof(g),typeof(weights_dict),typeof(weights_distribution_matrix)}( g, - weights, + weights_dict, weights_distribution_matrix ) end end - +function Base.show(io::IO, g::WeightedGraph) + print(io, "WG") +end function sample_mixing_graph!(mixing_graph,population_demographics) for (k,e) in enumerate(edges(mixing_graph.g)) i = src(e) @@ -74,8 +136,8 @@ function sample_mixing_graph!(mixing_graph,population_demographics) demo_j = Int(population_demographics[j]) contact_time = rand(RNG, mixing_graph.weights_distribution_matrix[demo_i,demo_j]) - mixing_graph.weights[(i,j)] = contact_time - mixing_graph.weights[(j,i)] = contact_time + mixing_graph.weights_dict[(i,j)] = contact_time + mixing_graph.weights_dict[(j,i)] = contact_time 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/model_setup.jl b/CovidAlertVaccinationModel/src/ABM/model_setup.jl index 6d623e82ac61196f3e620148b7b7513459732cc7..55ad33c8ff0ef06c34dbf1804c985f95dcc1359b 100644 --- a/CovidAlertVaccinationModel/src/ABM/model_setup.jl +++ b/CovidAlertVaccinationModel/src/ABM/model_setup.jl @@ -30,7 +30,7 @@ function get_u_0(nodes,I_0_fraction,vaccinator_prob) return status,is_vaccinator end -struct ModelSolution{T,G} +struct ModelSolution{T,InfNet,SocNet,WSMixingDist,RestMixingDist} sim_length::Int nodes::Int params::T @@ -40,8 +40,8 @@ struct ModelSolution{T,G} u_vac::Vector{Bool} covid_alert_times::Array{Int,2} time_of_last_alert::Vector{Int} - inf_network_lists::Vector{Vector{G}} - soc_networks::Vector{G} + inf_network::InfNet + soc_network::SocNet index_vectors::Vector{Vector{Int}} demographics::Vector{AgentDemographic} app_user::Vector{Bool} @@ -49,62 +49,32 @@ struct ModelSolution{T,G} app_user_index::Vector{Int} status_totals::Vector{Int} status_totals_next::Vector{Int} + ws_matrix_tuple::WSMixingDist + rest_matrix_tuple::RestMixingDist function ModelSolution(sim_length,params::T,num_households) where T demographics,base_network,index_vectors = generate_population(num_households) pop_sizes = length.(index_vectors) map_symmetrize(m_tuple) = map(md -> symmetrize_means(pop_sizes,md), m_tuple) - ws_matrix_list = map_symmetrize(workschool_mixing) - rest_matrix_list = map_symmetrize(rest_mixing) + ws_matrix_tuple = map_symmetrize(workschool_mixing) + rest_matrix_tuple = map_symmetrize(rest_mixing) app_user_list = zeros(length(demographics)) is_app_user = rand(RNG,length(demographics)) .< params.app_user_fraction app_user_index = findall(==(true),is_app_user) app_user_list[is_app_user] .= collect(1:length(app_user_index)) - # display(app) nodes = length(demographics) u_0_inf,u_0_vac = get_u_0(nodes,params.I_0_fraction,params.vaccinator_prob) - home_static_edges = WeightedGraph(base_network,contact_time_distributions.hh) #network with households and LTC homes - ws_static_edges = WeightedGraph(demographics,index_vectors,ws_matrix_list.daily,contact_time_distributions.ws) - ws_weekly_edges = WeightedGraph(demographics,index_vectors,ws_matrix_list.twice_a_week,contact_time_distributions.ws) - ws_daily_edges_vector = [WeightedGraph(demographics,index_vectors,ws_matrix_list.otherwise,contact_time_distributions.ws) for i in 1:sim_length] - - rest_static_edges = WeightedGraph(demographics,index_vectors,rest_matrix_list.daily,contact_time_distributions.rest) - rest_weekly_edges = WeightedGraph(demographics,index_vectors,rest_matrix_list.twice_a_week,contact_time_distributions.rest) - rest_daily_edges_vector = [WeightedGraph(demographics,index_vectors,rest_matrix_list.otherwise,contact_time_distributions.rest) for i in 1:sim_length] - - inf_network_lists = [ - [home_static_edges,rest_static_edges] for i in 1:sim_length - ] - soc_network_list = [home_static_edges,rest_static_edges,ws_static_edges] - - for (t,l) in enumerate(inf_network_lists) - day_of_week = mod(t,7) - if !(day_of_week == 3 || day_of_week == 4) #simulation begins on thursday I guess - push!(l, ws_static_edges) - end - if rand(RNG)<5/7 - push!(l, ws_weekly_edges) - push!(l, rest_weekly_edges) - end - push!(l,ws_daily_edges_vector[t]) - push!(l,rest_daily_edges_vector[t]) - end - + infected_mixing_graph,soc_mixing_graph = time_dep_mixing_graphs(sim_length,base_network,demographics,index_vectors,ws_matrix_tuple,rest_matrix_tuple) + covid_alert_times = zeros(Int,length(app_user_index),14) #two weeks worth of values - time_of_last_alert = fill(-1,length(app_user_index)) #two weeks worth of values - - status_totals = zeros(Int, AgentStatus.size) + time_of_last_alert = fill(-1,length(app_user_index)) #time of last alert is negative if no alert has been recieved - status_totals[1] = count(==(Susceptible), u_0_inf) - status_totals[2] = count(==(Infected), u_0_inf) - status_totals[3] = count(==(Recovered), u_0_inf) - status_totals[4] = count(==(Immunized), u_0_inf) + status_totals = [count(==(i), u_0_inf) for i in 1:AgentStatus.size] - - return new{T,typeof(home_static_edges)}( + return new{T,typeof(infected_mixing_graph),typeof(soc_mixing_graph),typeof(ws_matrix_tuple),typeof(rest_matrix_tuple)}( sim_length, nodes, params, @@ -114,15 +84,17 @@ struct ModelSolution{T,G} copy(u_0_vac), covid_alert_times, time_of_last_alert, - inf_network_lists, - soc_network_list, + infected_mixing_graph, + soc_mixing_graph, index_vectors, demographics, is_app_user, app_user_list, app_user_index, status_totals, - copy(status_totals) + copy(status_totals), + ws_matrix_tuple, + rest_matrix_tuple ) end end diff --git a/CovidAlertVaccinationModel/src/ABM/solve.jl b/CovidAlertVaccinationModel/src/ABM/solve.jl index 817aa6a9c4ed9f1f7edb84d143e2409436975345..fc8684de1af695e54171381de3ca475d3d434712 100644 --- a/CovidAlertVaccinationModel/src/ABM/solve.jl +++ b/CovidAlertVaccinationModel/src/ABM/solve.jl @@ -4,16 +4,16 @@ function contact_weight(p, contact_time) end function update_alert_durations!(t,modelsol) @unpack notification_parameter = modelsol.params - @unpack time_of_last_alert, app_user_index,inf_network_lists,covid_alert_times,app_user = modelsol + @unpack time_of_last_alert, app_user_index,inf_network,covid_alert_times,app_user = modelsol - for (i,node) in enumerate(modelsol.app_user_index), mixing_graph in modelsol.inf_network_lists[t] + for (i,node) in enumerate(modelsol.app_user_index), mixing_graph in modelsol.inf_network.graph_list[t] for j in 2:14 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[(node,j)] #add the contact times for today to the back + covid_alert_times[i,end] += mixing_graph.weights_dict[(node,j)] #add the contact times for today to the back end end if rand(RNG) < 1 - (1- notification_parameter)^sum(covid_alert_times[i,:]) @@ -23,7 +23,7 @@ function update_alert_durations!(t,modelsol) end function update_infection_state!(t,modelsol) @unpack base_transmission_probability,immunization_loss_prob,recovery_rate = modelsol.params - @unpack u_inf,u_vac,u_next_inf,u_next_vac,demographics,inf_network_lists,status_totals,status_totals_next = modelsol + @unpack u_inf,u_vac,u_next_inf,u_next_vac,demographics,inf_network,status_totals,status_totals_next = modelsol function agent_transition!(node, from::AgentStatus,to::AgentStatus) status_totals_next[Int(from)] -= 1 @@ -41,9 +41,9 @@ function update_infection_state!(t,modelsol) if is_vaccinator agent_transition!(i, Susceptible,Immunized) else - for mixing_graph in inf_network_lists[t] + for mixing_graph in inf_network.graph_list[t] for j in neighbors(mixing_graph.g,i) - if u_inf[j] == Infected && rand(RNG) < contact_weight(base_transmission_probability,mixing_graph.weights[(i,j)]) + if u_inf[j] == Infected && rand(RNG) < contact_weight(base_transmission_probability,mixing_graph.weights_dict[(i,j)]) agent_transition!(i, Susceptible,Infected) end end @@ -63,7 +63,7 @@ end function update_vaccination_opinion_state!(t,modelsol,total_infections) @unpack Ï€_base, η,γ, κ, ω, Ï, ω_en,Ï_en,γ,β = modelsol.params - @unpack demographics,time_of_last_alert, nodes, soc_networks,u_vac,u_next_vac,app_user,app_user_list = modelsol + @unpack demographics,time_of_last_alert, nodes, soc_network,u_vac,u_next_vac,app_user,app_user_list = modelsol app_user_pointer = 0 @@ -72,8 +72,8 @@ function update_vaccination_opinion_state!(t,modelsol,total_infections) soc_nbrs_vac = [0,0,0] soc_nbrs_nonvac = 0 num_soc_nbrs = 0 - for soc_network in soc_networks - soc_nbrs = neighbors(soc_network.g,i) + for sc_g in soc_network.graph_list[t] + soc_nbrs = neighbors(sc_g.g,i) num_soc_nbrs += length(soc_nbrs) for nbr in soc_nbrs if u_vac[nbr] @@ -108,7 +108,10 @@ end function agents_step!(t,modelsol) - for network in modelsol.inf_network_lists[t] + + remake!(modelsol.inf_network,modelsol.index_vectors,modelsol.ws_matrix_tuple.daily) + remake!(modelsol.inf_network,modelsol.index_vectors,modelsol.rest_matrix_tuple.daily) + for network in modelsol.inf_network.graph_list[t] sample_mixing_graph!(network,modelsol.demographics) #get new contact weights end diff --git a/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl b/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl index 6ef2ce693d8f2443d7add25c39352c0b7a0e54e1..4cc5eb67090e496582ec07a3b9adda0046a778d8 100644 --- a/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl +++ b/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl @@ -33,9 +33,9 @@ const RNG = Xoroshiro128Star(1) const color_palette = palette(:seaborn_pastel) #color theme for the plots include("utils.jl") include("data.jl") -include("ABM/agents.jl") include("ABM/mixing_distributions.jl") include("ABM/mixing_graphs.jl") +include("ABM/agents.jl") include("ABM/plotting.jl") include("ABM/model_setup.jl") include("ABM/output.jl")