Commit cb05a139 authored by Peter Jentsch's avatar Peter Jentsch
Browse files

resim intervalsmodel so output format is actually usable, covidalert usage in ABM

parent db98ecd1
......@@ -94,9 +94,9 @@ version = "0.12.7"
[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "4fecfd5485d3c5de4003e19f00c6898cccd40667"
git-tree-sha1 = "ac4132ad78082518ec2037ae5770b6e796f7f956"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.26.0"
version = "3.27.0"
[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
......@@ -255,16 +255,16 @@ uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89"
version = "3.3.3+0"
[[GR]]
deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "LinearAlgebra", "Pkg", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"]
git-tree-sha1 = "f74b42150042d11a1d94badf27c6125867c7e9bd"
deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"]
git-tree-sha1 = "82a03d9c331d181bf33c99c9af04202cc4533d48"
uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
version = "0.57.1"
version = "0.57.3"
[[GR_jll]]
deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt5Base_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "578527027f2d6c29a8de3e2eb6887d8850ef755c"
git-tree-sha1 = "90acee5c38f4933342fa9a3bbc483119d20e7033"
uuid = "d2c73de3-f751-5644-a686-071e5b155ba9"
version = "0.57.1+0"
version = "0.57.2+0"
[[GeometryBasics]]
deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
......@@ -793,9 +793,9 @@ uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics"]
git-tree-sha1 = "a83fa3021ac4c5a918582ec4721bc0cf70b495a9"
git-tree-sha1 = "4bc58880426274277a066de306ef19ecc22a6863"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.4"
version = "0.33.5"
[[StatsFuns]]
deps = ["Rmath", "SpecialFunctions"]
......@@ -805,9 +805,9 @@ version = "0.9.7"
[[StructArrays]]
deps = ["Adapt", "DataAPI", "Tables"]
git-tree-sha1 = "26ea43b4be7e919a2390c3c0f824e7eb4fc19a0a"
git-tree-sha1 = "44b3afd37b17422a62aea25f04c1f7e09ce6b07f"
uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.5.0"
version = "0.5.1"
[[StructTypes]]
deps = ["Dates", "UUIDs"]
......@@ -1029,7 +1029,7 @@ uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
version = "1.4.0+3"
[[ZeroWeightedDistributions]]
deps = ["Distributions", "Random"]
deps = ["Distributions", "Random", "StatsBase"]
path = "../ZeroWeightedDistributions"
uuid = "24733ad3-391a-4e41-8839-c7177de7dea4"
version = "0.1.0"
......
......@@ -7,6 +7,7 @@ version = "0.1.0"
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
......@@ -21,6 +22,8 @@ Pipe = "b98c9c47-44ae-5843-9183-064241ee97a0"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
RandomNumbers = "e6cf234a-135c-5ec9-84dd-332b85af5143"
Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d"
UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
......
......@@ -14,7 +14,7 @@ using NamedTupleTools
# using CUDA
using NetworkLayout:Stress
using NetworkLayout:SFDP
using ZeroWeightedDistributions
const PACKAGE_FOLDER = dirname(dirname(pathof(CovidAlertVaccinationModel)))
......@@ -41,13 +41,12 @@ function main()
agent_model = AgentModel(5000)
steps = 100
# return agent_model
sol,graphs = solve!(get_parameters(),steps,agent_model,vaccinate_uniformly!);
b1 = @benchmark solve!($get_parameters(),$steps,$agent_model,$vaccinate_uniformly!)
# display(b1)
# println("done")
# sol,graphs = solve_record!(u_0,get_parameters(),steps,agent_model,vaccinate_uniformly!);
sol,graphs = solve!(get_parameters(),steps,agent_model);
b1 = @benchmark solve!($get_parameters(),$steps,$agent_model) seconds = 20
display(b1)
println("done")
# sol,graphs = solve_record!(u_0,get_parameters(),steps,agent_model);
# return aggregate_timeseries(sol)
# plot_model_spatial_gif(agent_model.base_network,gddddddddraphs,sol)
end
end
......@@ -14,6 +14,8 @@ end
struct AgentModel{GraphType,L1,L2}
demographics::Vector{AgentDemographic}
demographic_index_vectors::Vector{Vector{Int}}
app_user::Vector{Bool}
app_user_index::Vector{Int}
base_network::GraphType
ws_matrix_list::L1
rest_matrix_list::L2
......@@ -25,7 +27,9 @@ struct AgentModel{GraphType,L1,L2}
workschool_contacts_mean_adjusted = map_symmetrize(workschool_mixing)
rest_contacts_mean_adjusted = map_symmetrize(rest_mixing)
# vaccinator_indices = sample(RNG,1:length(pop_list),length(pop_list) // 2)
app_user = rand(RNG,length(pop_list)) .< 0.1 #replace with parameter
app_user_index = findall(==(true),app_user)
return new{
typeof(base_network),
typeof(workschool_contacts_mean_adjusted),
......@@ -33,6 +37,8 @@ struct AgentModel{GraphType,L1,L2}
}(
pop_list,
index_vectors,
app_user,
app_user_index,
base_network,
workschool_contacts_mean_adjusted,
rest_contacts_mean_adjusted
......
function reindex!(k,csum,index_list_i,index_list_j,j_to_i_contacts,i_to_j_contacts,sample_list_i,sample_list_j)
@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)
i_index = index_list_i[k]
j_index = index_list_j[k]
csum += j_to_i_contacts[j_index] - i_to_j_contacts[i_index]
i_to_j_contacts[i_index] = sample_list_i[k]
csum += j_to_i_contacts[j_index] - i_to_j_contacts[i_index] + sample_list_i[k] - sample_list_j[k]
i_to_j_contacts[i_index] = sample_list_i[k]
j_to_i_contacts[j_index] = sample_list_j[k]
csum += i_to_j_contacts[i_index] - j_to_i_contacts[j_index]
return csum
end
function generate_contact_vectors!(ij_dist,ji_dist,i_to_j_contacts, j_to_i_contacts)
using StaticArrays
function generate_contact_vectors!(ij_dist,ji_dist,i_to_j_contacts::Vector{T}, j_to_i_contacts::Vector{T}) where T
rand!(RNG,ij_dist,i_to_j_contacts)
rand!(RNG,ji_dist,j_to_i_contacts)
l_i = length(i_to_j_contacts)
......@@ -19,17 +16,17 @@ function generate_contact_vectors!(ij_dist,ji_dist,i_to_j_contacts, j_to_i_conta
csum = sum(i_to_j_contacts) - sum(j_to_i_contacts)
inner_iter = 10
index_list_i = zeros(Int,inner_iter)
index_list_i = @MVector zeros(Int,inner_iter)
index_list_j = similar(index_list_i)
sample_list_i = zeros(eltype(i_to_j_contacts),inner_iter)
sample_list_i = @MVector zeros(T,inner_iter)
sample_list_j = similar(sample_list_i)
while csum != 0
index_list_i = sample!(RNG,1:l_i,index_list_i)
index_list_j = sample!(RNG,1:l_j,index_list_j)
sample_list_i = rand!(RNG,ij_dist,sample_list_i)
sample_list_j = rand!(RNG,ji_dist,sample_list_j)
for i = 1:inner_iter
sample!(RNG,1:l_i,index_list_i)
sample!(RNG,1:l_j,index_list_j)
rand!(RNG,ij_dist,sample_list_i)
rand!(RNG,ji_dist,sample_list_j)
@inbounds @simd for i = 1:inner_iter
if csum != 0
csum = reindex!(i,csum,index_list_i,index_list_j,j_to_i_contacts,i_to_j_contacts,sample_list_i,sample_list_j)
end
......
......@@ -79,4 +79,11 @@ function load_mixing_matrices()
end
end
return map(t -> from_mean.(Geometric{Float64},t),workschool_mixing), map(t -> from_mean.(Geometric{Float64}, t),rest_mixing)
end
\ No newline at end of file
end
# using Serialization
# using IntervalsModel
# function load_contact_time_distributions()
# dat = deserialize(joinpath(PACKAGE_FOLDER,"data/intervals_model_particle/hh.dat"))
# return dat
# end
\ No newline at end of file
......@@ -6,7 +6,7 @@ struct MixingContacts{V}
contact_array::V
function MixingContacts(index_vectors,mixing_matrix)
contacts = map(CartesianIndices(mixing_matrix)) do ind
zeros(length(index_vectors[ind[1]]))
zeros(Int,length(index_vectors[ind[1]]))
end
tot = 0
for i in 1:length(mixing_matrix[:,1]), j in 1:i #diagonal
......@@ -33,7 +33,7 @@ function random_bipartite_graph_fast_CL!(g::SimpleGraph,anodes,bnodes,aseq,bseq,
return g
end
using DataStructures
struct MixingGraph{G, V, M}
g::G
weights::V
......@@ -44,7 +44,7 @@ struct MixingGraph{G, V, M}
g = Graph(length(population_list))
weights = Dict{Tuple{Int,Int},UInt8}()
weights = RobinDict{Tuple{Int,Int},UInt8}()
for i in 1:length(index_vectors), j in 1:i #diagonal
random_bipartite_graph_fast_CL!(g,index_vectors[i],index_vectors[j],contacts.contact_array[i,j],contacts.contact_array[j,i],weights)
end
......@@ -57,7 +57,7 @@ struct MixingGraph{G, V, M}
)
end
function MixingGraph(g::SimpleGraph)
weights = Dict{Tuple{Int,Int},UInt8}()
weights = RobinDict{Tuple{Int,Int},UInt8}()
covid_alert_time = zeros(nv(g))
return new{typeof(g),typeof(weights),typeof(covid_alert_time)}(
g,
......@@ -69,12 +69,10 @@ end
function MixingGraph(g::SimpleGraph,population_list,contact_time_distribution_matrix)
mg = MixingGraph(g)
sample_mixing_graph!(mg, population_list, contact_time_distribution_matrix)
return mg
end
function MixingGraph(population_list,index_vectors,mixing_matrix,contact_time_distribution_matrix)
mg = MixingGraph(population_list,index_vectors,mixing_matrix)
sample_mixing_graph!(mg, population_list, contact_time_distribution_matrix)
return mg
end
......
......@@ -25,6 +25,13 @@ function from_mean_workschool(_,μ)
return Poisson(μ)
end
const workschool_mixing, rest_mixing = load_mixing_matrices()
const contact_time_distribution_matrix = [Geometric() for i in 1:(AgentDemographic.size-1), j in 1:(AgentDemographic.size-1)]
function alpha_matrix(alphas)
M = zeros(length(alphas),length(alphas))
for i in 1:length(alphas), j in 1:length(alphas)
......@@ -41,8 +48,27 @@ const unemployment_matrix = alpha_matrix(
)
)
const workschool_mixing, rest_mixing = load_mixing_matrices()
function make_workschool_mixing_matrix()
#all geometric with means given
ws_mixing = map(t->from_mean(t...),[
(Geometric{Float64}, 4.104848) (Geometric{Float64},2.568782) (Geometric{Float64},0.017729)
(Geometric{Float64}, 0.975688) (Geometric{Float64},5.057572) (Geometric{Float64},0.021307)
(Geometric{Float64},0.001937) (Geometric{Float64},0.00722) (Geometric{Float64}, 0.022134)
])
const contact_time_distribution_matrix = [Geometric() for i in 1:(AgentDemographic.size-1), j in 1:(AgentDemographic.size-1)]
#symmetrize WS mixing with respect to the population proportions in the data
ws_mixing_w_unemployment_symmetrized = symmetrize_means(get_household_data_proportions(),ws_mixing)
#define a function that adjusts the means of W according to the unemployment_matrix
ws_adjust_mean(W) = (5/7) .* (1 .- unemployment_matrix) ./ ( mean.(W) + (5/7) .* (1 .- unemployment_matrix))
#create a zero weighted distribution where the zero-weight is given by the unemployment_matrix, and the non-zero weight is given by ws_mixing, symmetrized, with means adjust upwards
ws_mixing_w_unemployment_symmetrized_weekday_adjusted = ZWDist.(unemployment_matrix,Geometric.(ws_adjust_mean(ws_mixing_w_unemployment_symmetrized)))
return ws_mixing_w_unemployment_symmetrized_weekday_adjusted
end
const initial_workschool_mixing_matrix = make_workschool_mixing_matrix()
@views function ws_sample(age)
return rand.(initial_workschool_mixing_matrix[age,:]) * (rand(RNG) < (5/7))
end
function get_u_0(nodes)
u_0 = fill(Susceptible,nodes)
#vaccinator with 50% probability
is_vaccinator = rand(nodes) .> 0.5
u_0 = (
status = fill(Susceptible,nodes),
is_vaccinator = is_vaccinator
)
init_indices = rand(RNG, 1:nodes, 5)
u_0[init_indices] .= Infected #start with two infected
u_0.status[init_indices] .= Infected #start with five infected
return u_0
end
......@@ -9,31 +16,54 @@ function contact_weight(p, contact_time)
return 1 - (1-p)^contact_time
end
function vaccinate_uniformly!(num_vaccines,u_next,u,graph,population_list,index_vectors) #vaccination_algorithm
vaccination_indices = rand(RNG,1:length(u),num_vaccines)
u_next[vaccination_indices] .= Immunized
function vac_switch_probability()
return 0.1
end
function hotspotting!(num_vaccines,u_next,u,graph,population_list,index_vectors) #vaccination_algorithm
# vaccination_indices = rand(RNG,1:length(u),num_vaccines)
# u_next[vaccination_indices] .= Immunized
end
function agents_step!(t,u_next,u,population_list,network_list,params,index_vectors,vaccination_algorithm!)
@unpack p, recovery_rate, vaccines_per_day,vaccination_start_day = params
u_next .= u
@inbounds for i in 1:length(u)
agent_status = u[i]
function agents_step!(t,u_next,u,population_list,network_list,params,covid_alert_times,app_user_index)
@unpack p, recovery_rate, immunization_loss_prob = params
u_next.status .= u.status
u_next.is_vaccinator .= u.is_vaccinator
for (i,node) in enumerate(app_user_index), mixing_graph in network_list
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)
covid_alert_times[i,end] += mixing_graph.weights[(node,j)] #add the contact times for today to the back
end
end
for i in 1:length(population_list)
if rand(RNG)<vac_switch_probability()
u_next.is_vaccinator[i] = !u.is_vaccinator[i]
end
end
for i in 1:length(population_list)
agent_status = u.status[i]
is_vaccinator = u.is_vaccinator[i]
agent_demo = population_list[i]
if agent_status == Susceptible
for mixing_graph in network_list
for j in neighbors(mixing_graph.g,i)
if u[j] == Infected && rand(RNG) < contact_weight(p,mixing_graph.weights[(i,j)])
u_next[i] = Infected
if is_vaccinator
u_next.status[i] = Immunized
else
for mixing_graph in network_list
for j in neighbors(mixing_graph.g,i)
if u.status[j] == Infected && rand(RNG) < contact_weight(p,mixing_graph.weights[(i,j)])
u_next.status[i] = Infected
end
end
end
end
elseif agent_status == Infected
if rand(RNG) < recovery_rate
u_next[i] = Recovered
u_next.status[i] = Recovered
end
elseif agent_status == Immunized
if rand(RNG) < immunization_loss_prob
u_next.status[i] = Susceptible
end
end
end
......@@ -41,12 +71,14 @@ end
function solve!(params,steps,agent_model,vaccination_algorithm!; record = false)
function solve!(params,steps,agent_model; record = false)
population_list = agent_model.demographics #list of demographic status for each agent
popsize = length(population_list)
index_vectors = agent_model.demographic_index_vectors #lists of indices of each agent with a given demographic
ws_contacts = agent_model.ws_matrix_list
rest_contacts = agent_model.rest_matrix_list
app_user_index = agent_model.app_user_index
p = params.p
u_0 = get_u_0(length(agent_model.demographics))
......@@ -61,12 +93,14 @@ function solve!(params,steps,agent_model,vaccination_algorithm!; record = false)
solution = vcat([u_0], [fill(Susceptible,length(u_0)) for i in 1:steps])
solution = vcat([u_0], [deepcopy(u_0) for i in 1:steps])
# @show solution
network_lists = [
[home_static_edges,rest_static_edges] for i in 1:steps
]
covid_alert_times = zeros(length(app_user_index),14) #two weeks worth of values
for (t,l) in enumerate(network_lists)
day_of_week = mod(t,7)
if !(day_of_week == 3 || day_of_week == 4) #simulation begins on thursday I guess
......@@ -82,8 +116,11 @@ function solve!(params,steps,agent_model,vaccination_algorithm!; record = false)
for t in 1:steps
network_list = network_lists[t]
#advance agent states based on the new network, vaccination process given by vaccination_algorithm!, which is just a function defined as above
agents_step!(t,solution[t+1],solution[t],population_list,network_list,params,index_vectors,vaccination_algorithm!)
for network in network_list
sample_mixing_graph!(network,population_list, contact_time_distribution_matrix) #get new contact weights
end
#advance agent states based on the new network
agents_step!(t,solution[t+1],solution[t],population_list,network_list,params,covid_alert_times,app_user_index)
end
return solution, network_lists
end
......@@ -92,7 +129,8 @@ function get_parameters()
p = 0.05,
recovery_rate = 0.1,
vaccines_per_day = 0.0, #percentage of population vaccinated per day
vaccination_start_day = 60
vaccination_start_day = 60,
immunization_loss_prob = 0.01
)
return params
end
......
This diff is collapsed.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment