diff --git a/CovidAlertVaccinationModel/Manifest.toml b/CovidAlertVaccinationModel/Manifest.toml
index eb8556d36e7ed83ee2d442ae96c37d6ab8db6a7f..e7787fd374febf3e457f49d2603b126aa0044c19 100644
--- a/CovidAlertVaccinationModel/Manifest.toml
+++ b/CovidAlertVaccinationModel/Manifest.toml
@@ -1229,6 +1229,12 @@ git-tree-sha1 = "293aa2c5cbf201e6b98810cb36d9eeafdafdafd1"
uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
version = "0.19.34"
+[[VectorizedRNG]]
+deps = ["Distributed", "Random", "UnPack", "VectorizationBase"]
+git-tree-sha1 = "f60638a5cdd839eba16fd9658909513a8291f862"
+uuid = "33b4df10-0173-11e9-2a0c-851a7edac40e"
+version = "0.2.8"
+
[[VersionParsing]]
git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f"
uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
diff --git a/CovidAlertVaccinationModel/Project.toml b/CovidAlertVaccinationModel/Project.toml
index 654a0383385af58661c7c754a06c26f15d89578a..0b775449a22244be70a6d01753f7f18fe29235f5 100644
--- a/CovidAlertVaccinationModel/Project.toml
+++ b/CovidAlertVaccinationModel/Project.toml
@@ -43,6 +43,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d"
UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
+VectorizedRNG = "33b4df10-0173-11e9-2a0c-851a7edac40e"
ZeroWeightedDistributions = "24733ad3-391a-4e41-8839-c7177de7dea4"
[compat]
diff --git a/CovidAlertVaccinationModel/abm_parameter_fits/outbreak_inf_parameters.dat b/CovidAlertVaccinationModel/abm_parameter_fits/outbreak_inf_parameters.dat
deleted file mode 100644
index 0a9c77acbdc26d0c06bb190a30a251756ea8e83b..0000000000000000000000000000000000000000
Binary files a/CovidAlertVaccinationModel/abm_parameter_fits/outbreak_inf_parameters.dat and /dev/null differ
diff --git a/CovidAlertVaccinationModel/abm_parameter_fits/seasonal_inf_parameters.dat b/CovidAlertVaccinationModel/abm_parameter_fits/seasonal_inf_parameters.dat
deleted file mode 100644
index 0d1edc9d05c43aba2ce2008755ad97be408348ce..0000000000000000000000000000000000000000
Binary files a/CovidAlertVaccinationModel/abm_parameter_fits/seasonal_inf_parameters.dat and /dev/null differ
diff --git a/CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl b/CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl
index dc9ebeb45643946d361f39ae4a1e837ed7e3742b..e6b4270b2e7632394495f00205da16137f7db32b 100644
--- a/CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl
+++ b/CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl
@@ -138,11 +138,11 @@ List of lists of graphs, one list for each day.
graph_list::Vector{Vector{G}}
"""
-struct TimeDepMixingGraph{N,G}
- resampled_graphs::NTuple{N,G}
+struct TimeDepMixingGraph{N,G,GNT}
+ resampled_graphs::NTuple{N,GNT}
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}(
+ function TimeDepMixingGraph(len,resampled_graphs::NTuple{N,GNT},base_graph_list::Vector{G}) where {GNT,G,N}
+ return new{N,G,GNT}(
resampled_graphs,
[copy(base_graph_list) for i in 1:len]
)
@@ -156,17 +156,20 @@ Assumes the simulation begins on Thursday arbitrarily.
"""
function time_dep_mixing_graphs(len,base_network,demographics,index_vectors,ws_matrix_tuple,rest_matrix_tuple)
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)
+
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)
+ infected_mixing_graph = TimeDepMixingGraph(len,((ws_daily_edges,ws_matrix_tuple.daily),(rest_daily_edges,rest_matrix_tuple.daily)),inf_network_list)
+ soc_mixing_graph = TimeDepMixingGraph(len,((ws_daily_edges,ws_matrix_tuple.daily),(rest_daily_edges,rest_matrix_tuple.daily)),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)
@@ -187,8 +190,8 @@ end
"""
Completely remake all the graphs in `time_dep_mixing_graph.resampled_graphs`.
"""
-function remake!(time_dep_mixing_graph,demographic_index_vectors,mixing_matrix)
- for weighted_graph in time_dep_mixing_graph.resampled_graphs
+function remake!(time_dep_mixing_graph,demographic_index_vectors)
+ for (weighted_graph,mixing_matrix) in time_dep_mixing_graph.resampled_graphs
empty!.(weighted_graph.g.fadjlist) #empty all the vector edgelists
weighted_graph.g.ne = 0
weighted_graph.mixing_edges = create_mixing_edges(demographic_index_vectors,mixing_matrix,weighted_graph.mixing_edges.sampler_matrix)
diff --git a/CovidAlertVaccinationModel/src/ABM/model_setup.jl b/CovidAlertVaccinationModel/src/ABM/model_setup.jl
index 901bd33dfc260e82b1bf810e5f29bd1acbeb7cce..2806d3c5dc7f29b10644e30a7dfe6098c5601dc4 100644
--- a/CovidAlertVaccinationModel/src/ABM/model_setup.jl
+++ b/CovidAlertVaccinationModel/src/ABM/model_setup.jl
@@ -1,7 +1,7 @@
function get_parameters()
params = (
- sim_length = 120,
+ sim_length = 500,
num_households = 5000,
I_0_fraction = 0.002,
base_transmission_probability = 0.001,
@@ -21,8 +21,8 @@ function get_parameters()
notification_threshold = 2,
immunizing = true,
immunization_delay = 14,
- immunization_begin_day = 60,
- infection_introduction_day = 120
+ immunization_begin_day =60,
+ infection_introduction_day = 180,
)
return params
end
diff --git a/CovidAlertVaccinationModel/src/ABM/output.jl b/CovidAlertVaccinationModel/src/ABM/output.jl
index 20e30853b781cdd577ec26599d4d4362c7747807..52d1b7d4a12d41cebdb5dad72fa9747943a45146 100644
--- a/CovidAlertVaccinationModel/src/ABM/output.jl
+++ b/CovidAlertVaccinationModel/src/ABM/output.jl
@@ -121,12 +121,12 @@ function plot_model(varname,univariate_series, output_list::Vector{T},infection_
sim_length = length(output_list[1].recorded_status_totals.S)
ts_list(data) = [
(data.recorded_status_totals.S, "Susceptible over time"),
- (data.recorded_status_totals.I, "Infected over time"),
(data.recorded_status_totals.R, "Recovered over time"),
- (data.recorded_status_totals.V, "Vaccinated over time"),
(data.total_vaccinators, "No. vaccinators over time"),
(data.mean_time_since_last_notification, "Mean time since last notification"),
- (sum.(eachcol(mean.(data.daily_cases_by_age))),"Daily (incident) cases"),
+ (data.daily_cases_by_age.Y,"Daily (incident) Y cases"),
+ (data.daily_cases_by_age.M,"Daily (incident) M cases"),
+ (data.daily_cases_by_age.O,"Daily (incident) O cases"),
(data.daily_immunized_by_age.Y, "new Y vaccinations each day"),
(data.daily_immunized_by_age.M, "new M vaccinations each day"),
(data.daily_immunized_by_age.O, "new O vaccinations each day"),
diff --git a/CovidAlertVaccinationModel/src/ABM/parameter_optimization.jl b/CovidAlertVaccinationModel/src/ABM/parameter_optimization.jl
index e866872c8f16c9c508fb88d8a12ab29bb6db7c97..55d0ce63ffd35a9e22d4de24c5c073ae2893f66f 100644
--- a/CovidAlertVaccinationModel/src/ABM/parameter_optimization.jl
+++ b/CovidAlertVaccinationModel/src/ABM/parameter_optimization.jl
@@ -3,10 +3,11 @@ using KissABC
const target_cumulative_vac_proportion = 0.33
const vaccination_data = [0.0,0.043,0.385,0.424,0.115,0.03,0.005] #by month starting in august
const ymo_vac = [0.255,0.278,0.602]
-const ymo_attack_rate = []
+const ymo_attack_rate = [10.376,5.636,7.2]
function solve_w_parameters(default_p, p_names, new_p_list)
+
new_params = merge(default_p, NamedTuple{p_names}(ntuple(i -> new_p_list[i],length(p_names))))
out = DebugRecorder(0,default_p.sim_length)
model = abm(new_params,out)
@@ -36,17 +37,17 @@ function fit_pre_inf_behavioural_parameters(p_tuple)
immunizing = true,
)
)
-
function cost(p)
output,model = solve_w_parameters(p_tuple_adjust, p_names,p)
target_cumulative_vaccinations = target_cumulative_vac_proportion*model.nodes
target_ymo_vac = ymo_vac .* sum(vaccination_data[1:4]) .* target_cumulative_vaccinations
- ymo_vaccination_ts = output.new_ymo_immunization
+ ymo_vaccination_ts = output.daily_immunized_by_age
total_preinfection_vaccination = sum.(eachrow(ymo_vaccination_ts))
# display(target_ymo_vac)
# display(total_preinfection_vaccination)
return sum((total_preinfection_vaccination .- target_ymo_vac).^2)
end
+
out = smc(priors,cost; verbose = true, nparticles = 800, parallel = true)#
return NamedTuple{p_names}(ntuple(i -> out.P[i].particles,length(p_names)))
end
@@ -68,22 +69,25 @@ function fit_post_inf_behavioural_parameters(p_tuple)
immunizing = true,
)
)
- target_cumulative_vac_proportion = 0.33
- vaccination_data = @SVector [0.0,0.043,0.385,0.424,0.115,0.03,0.005] #by month starting in august
- ymo_vac = @SVector [0.255,0.278,0.602]
function cost(p)
output,model = solve_w_parameters(p_tuple_adjust, p_names,p)
target_cumulative_vaccinations = target_cumulative_vac_proportion*model.nodes
target_ymo_vac = ymo_vac .* sum(vaccination_data[5:end]) .* target_cumulative_vaccinations
- ymo_vaccination_ts = output.new_ymo_immunization
+ ymo_vaccination_ts = output.daily_immunized_by_age
total_postinf_vaccination = sum.(eachrow(ymo_vaccination_ts[:,180:end]))
- final_size = output.recorded_status_totals.R[end]
- return sum((total_postinf_vaccination .- target_ymo_vac).^2)
+ final_size = sum.(eachrow(output.daily_cases_by_age))
+ target_final_size = ymo_attack_rate .* model.nodes
+
+ display(final_size)
+ display((sum((total_postinf_vaccination .- target_ymo_vac).^2) , 1e-3*sum((final_size .- target_final_size).^2)))
+ return sum((total_postinf_vaccination .- target_ymo_vac).^2) + sum((final_size .- target_final_size).^2)
end
- out =smc(priors,cost; verbose = true, nparticles = 200, parallel = true)# ABCDE(priors,cost,1e6; verbose=true, nparticles=300,generations=1000, parallel = true) #this one has better NaN handling
+
+ display(cost([0.0001,0.0005]))
+ # out =smc(priors,cost; verbose = true, nparticles = 200, parallel = true)# ABCDE(priors,cost,1e6; verbose=true, nparticles=300,generations=1000, parallel = true) #this one has better NaN handling
return NamedTuple{p_names}((out.P.particles,))
end
@@ -96,6 +100,7 @@ function plot_behavioural_fit(particles,p_tuple)
sim_length = sim_length,
I_0_fraction = 0.000,
immunization_begin_day =60,
+ infection_introduction_day = 180,
immunizing = true,
)
)
@@ -106,7 +111,7 @@ function plot_behavioural_fit(particles,p_tuple)
target_cumulative_vac_proportion = 0.33
vaccination_data = @SVector [0.0,0.043,0.385,0.424,0.115,0.03,0.005] #by month starting in august
ymo_vac = @SVector [0.255,0.278,0.602]
- ymo_vaccination_ts = mean.(out.new_ymo_immunization)
+ ymo_vaccination_ts = mean.(out.daily_immunized_by_age)
total_preinfection_vaccination = sum.(eachrow(ymo_vaccination_ts))
display(total_preinfection_vaccination)
p = [plot(),plot(),plot()]
@@ -120,34 +125,25 @@ end
# plot_max_posterior("outbreak", outbreak_transmission_dist,default_parameters)
function fit_parameters(default_parameters)
- seasonal_parameters_path = joinpath(PACKAGE_FOLDER,"abm_parameter_fits","outbreak_inf_parameters.dat")
pre_inf_behaviour_parameters_path =joinpath(PACKAGE_FOLDER,"abm_parameter_fits","pre_inf_behaviour_parameters.dat")
post_inf_behaviour_parameters_path = joinpath(PACKAGE_FOLDER,"abm_parameter_fits","post_inf_behaviour_parameters.dat")
- # seasonal_transmission_dist = CovidAlertVaccinationModel.fit_epi_parameters(default_parameters,0.073) ##seasonal
- # plot_max_posterior("seasonal", seasonal_transmission_dist,default_parameters)
- # pre_inf_behaviour_parameters = CovidAlertVaccinationModel.fit_pre_inf_behavioural_parameters(default_parameters)
- # fitted_parameters = map(mode,((;seasonal_transmission_dist...,pre_inf_behaviour_parameters...)))
- # fitted_parameter_tuple = merge(default_parameters,fitted_parameters)
- # serialize(seasonal_parameters_path,seasonal_transmission_dist)
+ # pre_inf_behaviour_parameters = fit_pre_inf_behavioural_parameters(default_parameters)
# serialize(pre_inf_behaviour_parameters_path, pre_inf_behaviour_parameters)
- # post_inf_behaviour_parameters = fit_post_inf_behavioural_parameters(fitted_parameter_tuple)
- # serialize(post_inf_behaviour_parameters_path, post_inf_behaviour_parameters)
-
+
+ pre_inf_behaviour_parameters = deserialize(pre_inf_behaviour_parameters_path)
+ post_inf_behaviour_parameters = fit_post_inf_behavioural_parameters(merge(default_parameters,map(mode,pre_inf_behaviour_parameters)))
+ serialize(post_inf_behaviour_parameters_path, post_inf_behaviour_parameters)
# visualize_Ï€_base(deserialize(pre_inf_behaviour_parameters_path))
-
-
fitted_parameters_with_post_inf_behaviour = (;
- deserialize(seasonal_parameters_path)...,
deserialize(pre_inf_behaviour_parameters_path)...,
deserialize(post_inf_behaviour_parameters_path)...
)
display(map(mode,fitted_parameters_with_post_inf_behaviour))
-
fitted_sol = plot_fitting_posteriors("post_inf_fitting",fitted_parameters_with_post_inf_behaviour,default_parameters)
return fitted_sol
end
diff --git a/CovidAlertVaccinationModel/src/ABM/solve.jl b/CovidAlertVaccinationModel/src/ABM/solve.jl
index b49c8083eb113c771d5461e4dedebf8e964895f0..aef858346ad639aa85cb67d181fadb18143bc7fa 100644
--- a/CovidAlertVaccinationModel/src/ABM/solve.jl
+++ b/CovidAlertVaccinationModel/src/ABM/solve.jl
@@ -122,14 +122,33 @@ Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,mod
end
+function weighted_degree(node,network::TimeDepMixingGraph)
+ weighted_degree = 0
+ for g_list in network.graph_list
+ for g in g_list
+ for j in neighbors(g.g,node)
+ weighted_degree += get_weight(g,GraphEdge(node,j))
+ end
+ end
+ end
+ return weighted_degree
+end
function agents_step!(t,modelsol)
- remake!(modelsol.inf_network,modelsol.index_vectors,modelsol.ws_matrix_tuple.daily)
- remake!(modelsol.soc_network,modelsol.index_vectors,modelsol.rest_matrix_tuple.daily)
- for network in modelsol.inf_network.graph_list[t]
+ remake!(modelsol.inf_network,modelsol.index_vectors)
+ for network in modelsol.inf_network.graph_list[t] #this also resamples the soc network weights since they point to the same objects, but those are never used
sample_mixing_graph!(network) #get new contact weights
end
+ # @show modelsol.inf_network.graph_list[t]
+ # @show weighted_degree(1,modelsol.inf_network)
+ # @show weighted_degree(1,modelsol.soc_network)
+
+
+
+
+
+
if t == modelsol.params.infection_introduction_day
init_indices = rand(Random.default_rng(Threads.threadid()), 1:modelsol.nodes, round(Int,modelsol.nodes*modelsol.params.I_0_fraction))
modelsol.u_inf[init_indices] .= Infected
diff --git a/timeseries.pdf b/timeseries.pdf
index 8d4307ce48f21b64e04bbda1880b3ea4bdb684f3..0c08e5343c5f98a5d66a2802b97baff66bb2eac9 100644
Binary files a/timeseries.pdf and b/timeseries.pdf differ