Skip to content
GitLab
Commit 91c68770 authored by Peter Jentsch's avatar Peter Jentsch
Browse files

new output plots, alert mechanic bug fixed and mechanic changed to give alerts... 

new output plots, alert mechanic bug fixed and mechanic changed to give alerts if coin flips exceed notification_threshold
parent 8bff20c7
Hide whitespace changes
Inline Side-by-side
CovidAlertVaccinationModel/Manifest.toml
View file @ 91c68770
......@@ -171,9 +171,9 @@ version = "1.6.0"
[[DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "56ff5833e5b755d2db654479993e949e73606b64"
git-tree-sha1 = "78028b8c1c5147515e0a367e6bbe679fc86ff085"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.0.0"
version = "1.0.1"
[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
......@@ -502,9 +502,9 @@ version = "1.6.0"
[[Latexify]]
deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
git-tree-sha1 = "1925f6838df247e7853f3f9727dd8a52a78f60f4"
git-tree-sha1 = "f77a16cb3804f4a74f57e5272a6a4a9a628577cb"
uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
version = "0.15.2"
version = "0.15.5"
[[Lazy]]
deps = ["MacroTools"]
......@@ -602,14 +602,14 @@ uuid = "093fc24a-ae57-5d10-9952-331d41423f4d"
version = "1.3.5"
[[LinearAlgebra]]
deps = ["Libdl"]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
[[LogExpFunctions]]
deps = ["DocStringExtensions"]
git-tree-sha1 = "49c5c32deda5999d15378b64ee10f2e87831ab25"
deps = ["DocStringExtensions", "LinearAlgebra"]
git-tree-sha1 = "ed26854d7c2c867d143f0e07c198fc9e8b721d10"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.2.2"
version = "0.2.3"
[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
......@@ -718,6 +718,22 @@ git-tree-sha1 = "a42c0f138b9ebe8b58eba2271c5053773bde52d0"
uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
version = "1.3.4+2"
[[OnlineStats]]
deps = ["AbstractTrees", "Dates", "LinearAlgebra", "OnlineStatsBase", "OrderedCollections", "Random", "RecipesBase", "Statistics", "StatsBase", "SweepOperator"]
git-tree-sha1 = "a3d0e72ffde4f596cc148c4c6555e310dcd5b2ba"
uuid = "a15396b6-48d5-5d58-9928-6d29437db91e"
version = "1.5.8"
[[OnlineStatsBase]]
deps = ["AbstractTrees", "Dates", "LinearAlgebra", "OrderedCollections", "Statistics", "StatsBase"]
git-tree-sha1 = "ed076aedbd3cb85731730268d60a7acfa4ac92f3"
uuid = "925886fa-5bf2-5e8e-b522-a9147a512338"
version = "1.4.4"
[[OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
[[OpenSSL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "71bbbc616a1d710879f5a1021bcba65ffba6ce58"
......@@ -726,9 +742,9 @@ version = "1.1.1+6"
[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "9db77584158d0ab52307f8c04f8e7c08ca76b5b3"
git-tree-sha1 = "b9b8b8ed236998f91143938a760c2112dceeb2b4"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.3+4"
version = "0.5.4+0"
[[Opus_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
......@@ -998,9 +1014,9 @@ uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics"]
git-tree-sha1 = "4bc58880426274277a066de306ef19ecc22a6863"
git-tree-sha1 = "4d8ca45223d7a28839e775d73a6f6b6b2ac64fd1"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.5"
version = "0.33.6"
[[StatsFuns]]
deps = ["LogExpFunctions", "Rmath", "SpecialFunctions"]
......@@ -1021,9 +1037,19 @@ uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.5.1"
[[SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays", "SuiteSparse_jll"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
[[SuiteSparse_jll]]
deps = ["Artifacts", "Libdl", "Pkg", "libblastrampoline_jll"]
uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c"
[[SweepOperator]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "20da2784e79fc0bfdd70592d1b47d7a6034e82d1"
uuid = "7522ee7d-7047-56d0-94d9-4bc626e7058d"
version = "0.3.0"
[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
......@@ -1103,9 +1129,9 @@ uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
[[VectorizationBase]]
deps = ["ArrayInterface", "Hwloc", "IfElse", "Libdl", "LinearAlgebra", "Static"]
git-tree-sha1 = "4f9b7b3e40da418518e0282e7397fd0ca17a7527"
git-tree-sha1 = "5d66756d2d5c538e9ab0c39efcb4512e3a8c3361"
uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
version = "0.19.28"
version = "0.19.31"
[[VersionParsing]]
git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f"
......@@ -1296,6 +1322,10 @@ git-tree-sha1 = "acc685bcf777b2202a904cdcb49ad34c2fa1880c"
uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
version = "0.14.0+4"
[[libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll", "Pkg"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
[[libfdk_aac_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7a5780a0d9c6864184b3a2eeeb833a0c871f00ab"
......
CovidAlertVaccinationModel/Project.toml
View file @ 91c68770
......@@ -23,6 +23,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
NamedTupleTools = "d9ec5142-1e00-5aa0-9d6a-321866360f50"
NetworkLayout = "46757867-2c16-5918-afeb-47bfcb05e46a"
OnlineStats = "a15396b6-48d5-5d58-9928-6d29437db91e"
Pandas = "eadc2687-ae89-51f9-a5d9-86b5a6373a9c"
Pipe = "b98c9c47-44ae-5843-9183-064241ee97a0"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
......
CovidAlertVaccinationModel/abm_test.jl
View file @ 91c68770
......@@ -2,45 +2,47 @@ using CovidAlertVaccinationModel
using OnlineStats
using ThreadsX
using Plots
const samples = 5
const samples = 10
##Univariate tests
const len = 4 #number of points to evaluate
gr()
const default_parameters = (
sim_length = 400,
sim_length = 600,
num_households = 5000,
I_0_fraction = 0.001,
base_transmission_probability = 0.005,
base_transmission_probability = 0.0015,
recovery_rate = 0.1,
immunization_loss_prob = 0.0055, #mean time of 6 months
π_base = -2.0,
π_base = -0.25,
η = 0.0,
κ = 0.0,
ω = 0.0,
ω = 0.0005,
ρ = [0.0,0.0,0.0],
ω_en = 0.02,
ω_en = 0.00,
ρ_en = [0.0,0.0,0.0],
γ = 0.0,
β = 5.0,
notification_parameter = 0.1,
notification_parameter = 0.001,
vaccinator_prob = 0.2,
app_user_fraction = 0.4,
notification_threshold = 2,
)
const univarate_test_list = (
(:I_0_fraction, range(0.0, 0.05; length = len)),
(:base_transmission_probability, range(0.001, 0.02; length = len)),
(:recovery_rate, range(0.05, 0.5; length = len)),
(:base_transmission_probability, range(0.0001, 0.002; length = len)),
(:recovery_rate, range(0.1, 0.5; length = len)),
(:immunization_loss_prob, range(0.00, 0.05; length = len)),
(:π_base, range(-0.1, 0.1; length = len)),
(:η, range(0.0, 1.0; length = len)),
(:π_base, range(-1.0, -0.0; length = len)),
(:η, range(0.0, 0.1; length = len)),
(:κ, range(0.0, 0.1; length = len)),
(:ω, range(0.0, 0.025; length = len)),
(:ω, range(0.0, 0.001; length = len)),
(:ω_en, range(0.0, 0.5; length = len)),
(:γ, range(0.0, 0.5; length = len)),
(:notification_parameter, range(0.1, 1.0; length = len)),
(:app_user_fraction, range(0.05, 0.5; length = len)),
(:notification_parameter, range(0.00, 0.05; length = len)),
(:app_user_fraction, range(0.05, 0.25; length = len)),
(:notification_threshold, (1:len)),
)
const univariate_path = "CovidAlertVaccinationModel/plots/univariate/"
......
CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl
View file @ 91c68770
......@@ -35,11 +35,13 @@ Resample all the weights in `mixing_graph`
"""
function sample_mixing_graph!(mixing_graph)
mixing_edges = mixing_graph.mixing_edges
for i in 1:size(mixing_edges.contact_array)[1], j in 1:i #diagonal
rand!(RNG, mixing_edges.sampler_matrix[i,j],mixing_edges.sample_cache[i,j])
for k in 1:length(mixing_edges.contact_array[i,j])
edge_weight_k = mixing_edges.sample_cache[i,j][k]
set!(mixing_edges.weights_dict, mixing_edges.contact_array[i,j][k], edge_weight_k)
# display(length.(mixing_edges.contact_array))
# display(length.(mixing_edges.sample_cache))
for i in 1:size(mixing_edges.contact_array)[1], j in 1:i #diagonal
rand!(RNG, mixing_edges.sampler_matrix[j,i],mixing_edges.sample_cache[j,i])
for k in 1:length(mixing_edges.contact_array[j,i])
edge_weight_k = mixing_edges.sample_cache[j,i][k]
set!(mixing_edges.weights_dict, mixing_edges.contact_array[j,i][k], edge_weight_k)
end
end
end
......@@ -116,7 +118,7 @@ function create_mixing_edges(g::SimpleGraph,demographics,demographic_index_vecto
for e in edges(g)
i = src(e)
j = dst(e)
push!(contact_array[Int(demographics[j]),Int(demographics[i])], GraphEdge(i,j))
push!(contact_array[Int(demographics[i]),Int(demographics[j])], GraphEdge(i,j))
end
sample_cache = map(v-> Vector{Int}(undef,length(v)),contact_array)
return MixingEdges(ne(g),contact_array,weights_distribution_matrix)
......
CovidAlertVaccinationModel/src/ABM/model_setup.jl
View file @ 91c68770
......@@ -3,29 +3,30 @@ function get_parameters()
params = (
sim_length = 500,
num_households = 5000,
I_0_fraction = 0.005,
base_transmission_probability = 0.01,
I_0_fraction = 0.001,
base_transmission_probability = 0.005,
recovery_rate = 0.1,
immunization_loss_prob = 0.0055,
immunization_loss_prob = 0.0055, #mean time of 6 months
π_base = -1.0,
η = 0.0,
η = 0.01,
κ = 0.0,
ω = 0.0,
ω = 0.00,
ρ = [0.0,0.0,0.0],
ω_en = 0.01,
ω_en = 0.001,
ρ_en = [0.0,0.0,0.0],
γ = 0.0,
β = 5.0,
notification_parameter = 0.2,
notification_parameter = 0.001,
vaccinator_prob = 0.2,
app_user_fraction = 0.4,
app_user_fraction = 0.5,
notification_threshold = 8,
)
return params
end
function get_u_0(nodes,I_0_fraction,vaccinator_prob)
#vaccinator with 50% probability
is_vaccinator = rand(nodes) .> vaccinator_prob
is_vaccinator = rand(nodes) .< vaccinator_prob
status = fill(Susceptible,nodes)
init_indices = rand(RNG, 1:nodes, round(Int,nodes*I_0_fraction))
status[init_indices] .= Infected #start with five infected
......@@ -40,7 +41,7 @@ struct ModelSolution{T,InfNet,SocNet,WSMixingDist,RestMixingDist}
u_next_vac::Vector{Bool}
u_inf::Vector{AgentStatus}
u_vac::Vector{Bool}
covid_alert_times::Array{Int,2}
covid_alert_notifications::Array{Bool,2}
time_of_last_alert::Vector{Int}
inf_network::InfNet
soc_network::SocNet
......@@ -71,7 +72,7 @@ struct ModelSolution{T,InfNet,SocNet,WSMixingDist,RestMixingDist}
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,14,length(app_user_index)) #two weeks worth of values
covid_alert_notifications = zeros(Bool,14,length(app_user_index)) #two weeks worth of values
time_of_last_alert = fill(-1,length(app_user_index)) #time of last alert is negative if no alert has been recieved
status_totals = [count(==(AgentStatus(i)), u_0_inf) for i in 1:AgentStatus.size]
......@@ -84,7 +85,7 @@ struct ModelSolution{T,InfNet,SocNet,WSMixingDist,RestMixingDist}
u_0_vac,
copy(u_0_inf),
copy(u_0_vac),
covid_alert_times,
covid_alert_notifications,
time_of_last_alert,
infected_mixing_graph,
soc_mixing_graph,
......
CovidAlertVaccinationModel/src/ABM/output.jl
View file @ 91c68770
......@@ -57,9 +57,13 @@ end
function record!(t,modelsol, recorder::DebugRecorder)
recorder.total_vaccinators[t] = count(==(true),modelsol.u_vac)
recorder.recorded_status_totals[:,t] .= modelsol.status_totals
mean_alert_time = mean(t .- modelsol.time_of_last_alert)
recorder.mean_time_since_last_notification[t] = round(Int,mean_alert_time)
alerts = filter(>(0),modelsol.time_of_last_alert)
if !isempty(alerts)
mean_alert_time = mean(t .- alerts)
recorder.mean_time_since_last_notification[t] = round(Int,mean_alert_time)
else
recorder.mean_time_since_last_notification[t] = 0
end
end
function record!(t,modelsol, recorder::Nothing)
......@@ -92,44 +96,50 @@ function plot_model(varname,univariate_series, output_list::Vector{T}) where T<:
plts = [plot(),plot(),plot(),plot(),plot(),plot()]
for (i,(p,data)) in enumerate(zip(univariate_series, output_list))
# display(p[varname])
p_val = @sprintf "%.3f" (p[varname])
p_val = @sprintf "%.4f" (p[varname])
# display(p_val)
plot!(plts[1], mean.(data.recorded_status_totals.S); ribbon = std.(data.recorded_status_totals.S),
label = "$varname = $(p_val)", line_z = i, color=:blues,
ylabel = "no. of people",
colorbar = false,
title = "Susceptible over time"
title = "Susceptible over time",
legend=:topright
)
plot!(plts[2], mean.(data.recorded_status_totals.I); ribbon = std.(data.recorded_status_totals.I),
label = "$varname = $(p_val)", line_z = i, color=:blues,
ylabel = "no. of people",
colorbar = false,
title = "Infected over time"
title = "Infected over time",
legend=false,
)
plot!(plts[3], mean.(data.recorded_status_totals.R); ribbon = std.(data.recorded_status_totals.R),
label = "$varname = $(p_val)", line_z = i, color=:blues,
ylabel = "no. of people",
colorbar = false,
title = "Recovered over time"
title = "Recovered over time",
legend=false
)
plot!(plts[4], mean.(data.recorded_status_totals.V); ribbon = std.(data.recorded_status_totals.V),
label = "$varname = $(p_val)", line_z = i, color=:blues,
ylabel = "no. of people",
colorbar = false,
title = "Vaccinated over time"
title = "Vaccinated over time",
legend=false
)
plot!(plts[5], mean.(data.total_vaccinators); ribbon = std.(data.total_vaccinators),
label = "$varname = $(p_val)", line_z = i, color=:blues,
ylabel = "no. of people",
colorbar = false,
title = "No. vaccinators over time"
title = "No. vaccinators over time",
legend=false
)
plot!(plts[6], mean.(data.mean_time_since_last_notification); ribbon = std.(data.mean_time_since_last_notification),
label = "$varname = $(p_val)", line_z = i, color=:blues,
ylabel = "Days",
colorbar = false,
title = "Mean time since last notification"
)
label = "$varname = $(p_val)", line_z = i, color=:blues,
ylabel = "Days",
colorbar = false,
title = "Mean time since last notification",
legend=false
)
end
return plot(plts...;layout = (6,1),size=(800,2500),leftmargin = 5Plots.mm)
return plot(plts...;layout = (6,1),size=(800,2500),leftmargin = 8Plots.mm)
end
CovidAlertVaccinationModel/src/ABM/solve.jl
View file @ 91c68770
#remove Base.@propagate_inbounds during tests if you get segfaults or change anything,
# as that might mean there is an indexing bug that is not being caught,
#and therefore accessing the wrong memory is causing background errors
function contact_weight(p, contact_time)
return 1 - (1-p)^contact_time
......@@ -10,35 +11,36 @@ function Φ(payoff,β)
return 1 / (exp(-1*β*payoff))
end
Base.@propagate_inbounds @views function update_alert_durations!(t,modelsol)
#remove Base.@propagate_inbounds if you get segfaults
@unpack notification_parameter = modelsol.params
@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 inf_network.graph_list[t]
Base.@propagate_inbounds @views function update_alert_durations!(t,modelsol) # Base.@propagate_inbounds
@unpack notification_parameter,notification_threshold = modelsol.params
@unpack time_of_last_alert, app_user_index,inf_network,covid_alert_notifications,app_user = modelsol
for (i,node) in enumerate(modelsol.app_user_index)
for j in 2:14
covid_alert_times[j-1,i] = covid_alert_times[j,i] #shift them all back
covid_alert_notifications[j-1,i] = covid_alert_notifications[j,i] #shift them all back
end
for j in neighbors(mixing_graph.g,node)
if app_user[j]
covid_alert_times[end,i] += get_weight(mixing_graph,GraphEdge(node,j)) #add the contact times for today to the back
total_weight_i = 0
for mixing_graph in inf_network.graph_list[t]
for j in neighbors(mixing_graph.g,node)
if app_user[j]
total_weight_i+= get_weight(mixing_graph,GraphEdge(node,j))
end
end
end
covid_alert_total_exposures = 1 - (1 - notification_parameter) ^ sum(covid_alert_times[:,i])
if rand(RNG) < covid_alert_total_exposures
coin_flip = 1 - (1 - notification_parameter)^total_weight_i
r = rand(RNG)
if r < coin_flip
covid_alert_notifications[end,i] = 1 #add the notifications for today
else
covid_alert_notifications[end,i] = 0
end
if sum(covid_alert_notifications[:,i])>=notification_threshold
time_of_last_alert[i] = t
end
end
# display(time_of_last_alert)
end
Base.@propagate_inbounds @views function update_infection_state!(t,modelsol)
#remove Base.@propagate_inbounds if you get segfaults
@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,status_totals = modelsol
......@@ -79,16 +81,12 @@ Base.@propagate_inbounds @views function update_infection_state!(t,modelsol)
end
end
end
# display(u_next_inf)
end
Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,modelsol,total_infections)
#remove Base.@propagate_inbounds if you get segfaults
@unpack π_base, η,γ, κ, ω, ρ, ω_en,ρ_en,γ,β = modelsol.params
@unpack demographics,time_of_last_alert, nodes, soc_network,u_vac,u_next_vac,app_user,app_user_list = modelsol
app_user_pointer = 0
for i in 1:nodes
vac_payoff = 0
soc_nbrs_vac = @MArray [0,0,0]
......
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