diff --git a/CovidAlertVaccinationModel/Manifest.toml b/CovidAlertVaccinationModel/Manifest.toml
index c0aed9a0a9d3dfbbe9af4375b98dafa0ed3b1bc4..149934dba10ae26ebae41ff8d6d9db0d9ded0735 100644
--- a/CovidAlertVaccinationModel/Manifest.toml
+++ b/CovidAlertVaccinationModel/Manifest.toml
@@ -8,9 +8,9 @@ version = "1.0.1"
[[AbstractMCMC]]
deps = ["BangBang", "ConsoleProgressMonitor", "Distributed", "Logging", "LoggingExtras", "ProgressLogging", "Random", "StatsBase", "TerminalLoggers", "Transducers"]
-git-tree-sha1 = "29683bc1b52e1879ac0951253d8b0e2f60bf4cb4"
+git-tree-sha1 = "21279159f6be4b2fd00e1a4a1f736893100408fc"
uuid = "80f14c24-f653-4e6a-9b94-39d6b0f70001"
-version = "3.1.0"
+version = "3.2.0"
[[AbstractTrees]]
git-tree-sha1 = "03e0550477d86222521d254b741d470ba17ea0b5"
@@ -92,15 +92,15 @@ version = "1.16.0+6"
[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
-git-tree-sha1 = "e6b23566e025d3b0d9ccc397f5c7a134af552e27"
+git-tree-sha1 = "9b0375dc013ab0fc472b37cb8b18eed66b83f76b"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-version = "0.9.42"
+version = "0.9.43"
[[CheapThreads]]
deps = ["ArrayInterface", "IfElse", "Requires", "Static", "StrideArraysCore", "ThreadingUtilities", "VectorizationBase"]
-git-tree-sha1 = "6596dbdb8fadd45ef9dff0087ef3a6ec9c5473bc"
+git-tree-sha1 = "44899e10c2fd8a3c5b907c4b688c04f737d13959"
uuid = "b630d9fa-e28e-4980-896d-83ce5e2106b2"
-version = "0.2.3"
+version = "0.2.4"
[[ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random", "StaticArrays"]
@@ -128,9 +128,9 @@ version = "0.3.0"
[[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 = "0a817fbe51c976de090aa8c997b7b719b786118d"
+git-tree-sha1 = "0900bc19193b8e672d9cd477e6cd92d9e7c02f99"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "3.28.0"
+version = "3.29.0"
[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
@@ -322,9 +322,9 @@ version = "2.8.0"
[[FiniteDifferences]]
deps = ["ChainRulesCore", "LinearAlgebra", "Printf", "Random", "Richardson", "StaticArrays"]
-git-tree-sha1 = "80e1a7416cbf08fe80c8885e1834c45cfc399c61"
+git-tree-sha1 = "b97a0b1adc3f10482d704158073cf3ac77493e2c"
uuid = "26cc04aa-876d-5657-8c51-4c34ba976000"
-version = "0.12.5"
+version = "0.12.6"
[[FixedPointNumbers]]
deps = ["Statistics"]
@@ -678,9 +678,9 @@ version = "0.4.6"
[[LoopVectorization]]
deps = ["ArrayInterface", "CheapThreads", "DocStringExtensions", "IfElse", "LinearAlgebra", "OffsetArrays", "Requires", "SLEEFPirates", "Static", "StrideArraysCore", "ThreadingUtilities", "UnPack", "VectorizationBase"]
-git-tree-sha1 = "2ad016117de05750443ed219dada9df3f9b9fa8f"
+git-tree-sha1 = "8fb5abd2ac992013fa14b110ca1308c83f8a4098"
uuid = "bdcacae8-1622-11e9-2a5c-532679323890"
-version = "0.12.18"
+version = "0.12.19"
[[LsqFit]]
deps = ["Distributions", "ForwardDiff", "LinearAlgebra", "NLSolversBase", "OptimBase", "Random", "StatsBase"]
@@ -1025,9 +1025,9 @@ uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
[[SLEEFPirates]]
deps = ["IfElse", "Libdl", "VectorizationBase"]
-git-tree-sha1 = "91a650350dcf6e0fc1a014b59669e704d8f579ae"
+git-tree-sha1 = "e262a4909eda52b6d9f1d1f0eac6537d8267a4b5"
uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa"
-version = "0.6.17"
+version = "0.6.19"
[[Scratch]]
deps = ["Dates"]
@@ -1080,10 +1080,10 @@ deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
[[SpecialFunctions]]
-deps = ["ChainRulesCore", "OpenSpecFun_jll"]
-git-tree-sha1 = "5919936c0e92cff40e57d0ddf0ceb667d42e5902"
+deps = ["ChainRulesCore", "LogExpFunctions", "OpenSpecFun_jll"]
+git-tree-sha1 = "9146da51b38e9705b9f5ccfadc3ab10a482cae36"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "1.3.0"
+version = "1.4.0"
[[SplittablesBase]]
deps = ["Setfield", "Test"]
@@ -1099,9 +1099,9 @@ version = "0.2.4"
[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
-git-tree-sha1 = "fb46e45ef2cade8be20bb445b3ffeca3c6d6f7d3"
+git-tree-sha1 = "c635017268fd51ed944ec429bcc4ad010bcea900"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.1.3"
+version = "1.2.0"
[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
@@ -1126,9 +1126,9 @@ version = "0.9.8"
[[StrideArraysCore]]
deps = ["ArrayInterface", "Requires", "ThreadingUtilities", "VectorizationBase"]
-git-tree-sha1 = "f93118d367c8dec873c26a32ad2dea84989edd7d"
+git-tree-sha1 = "72429259dd8bbab3ee6ac78e3967e5409460cfc1"
uuid = "7792a7ef-975c-4747-a70f-980b88e8d1da"
-version = "0.1.7"
+version = "0.1.9"
[[StructArrays]]
deps = ["Adapt", "DataAPI", "Tables"]
@@ -1184,9 +1184,9 @@ uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
[[ThreadingUtilities]]
deps = ["VectorizationBase"]
-git-tree-sha1 = "063f52eee44ec303f1721cd59b4d7892cae9f1cc"
+git-tree-sha1 = "d9117912ec78dbba3294fcb2962b6826be6107a5"
uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5"
-version = "0.4.1"
+version = "0.4.2"
[[ThreadsX]]
deps = ["ArgCheck", "BangBang", "ConstructionBase", "InitialValues", "MicroCollections", "Referenceables", "Setfield", "SplittablesBase", "Transducers"]
@@ -1202,9 +1202,9 @@ version = "1.5.5"
[[Transducers]]
deps = ["Adapt", "ArgCheck", "BangBang", "Baselet", "CompositionsBase", "DefineSingletons", "Distributed", "InitialValues", "Logging", "Markdown", "MicroCollections", "Requires", "Setfield", "SplittablesBase", "Tables"]
-git-tree-sha1 = "c277f1190f76f108cfdb89b9d5da87d9602e5593"
+git-tree-sha1 = "34f27ac221cb53317ab6df196f9ed145077231ff"
uuid = "28d57a85-8fef-5791-bfe6-a80928e7c999"
-version = "0.4.64"
+version = "0.4.65"
[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
@@ -1225,15 +1225,17 @@ uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
[[VectorizationBase]]
deps = ["ArrayInterface", "Hwloc", "IfElse", "Libdl", "LinearAlgebra", "Static"]
-git-tree-sha1 = "c258467e1a3473e328c8a9109efcce845723593e"
+git-tree-sha1 = "24a5c54dab1907b0f2f086e5d15c65522c7e8ce1"
uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
-version = "0.19.37"
+version = "0.20.3"
[[VectorizedRNG]]
deps = ["Distributed", "Random", "UnPack", "VectorizationBase"]
-git-tree-sha1 = "f60638a5cdd839eba16fd9658909513a8291f862"
+git-tree-sha1 = "0a65169fbbf9c94d2136ba75946e73a88af3b491"
+repo-rev = "master"
+repo-url = "https://github.com/JuliaSIMD/VectorizedRNG.jl.git"
uuid = "33b4df10-0173-11e9-2a0c-851a7edac40e"
-version = "0.2.8"
+version = "0.2.11"
[[VersionParsing]]
git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f"
diff --git a/CovidAlertVaccinationModel/abm_parameter_fits/fit_all_parameters.dat b/CovidAlertVaccinationModel/abm_parameter_fits/fit_all_parameters.dat
new file mode 100644
index 0000000000000000000000000000000000000000..3e797eb7dd9642396b7242477e837c3fcb3b3dd5
Binary files /dev/null and b/CovidAlertVaccinationModel/abm_parameter_fits/fit_all_parameters.dat differ
diff --git a/CovidAlertVaccinationModel/src/ABM/contact_vectors.jl b/CovidAlertVaccinationModel/src/ABM/contact_vectors.jl
index 3937cae491e683eaf4ea6431777798b72b4e344c..412e1c1faa4725c76596e15924ccbcf0d5b619e3 100644
--- a/CovidAlertVaccinationModel/src/ABM/contact_vectors.jl
+++ b/CovidAlertVaccinationModel/src/ABM/contact_vectors.jl
@@ -27,17 +27,9 @@ Fills i_to_j_contacts and j_to_i_contacts with degrees sampled from ij_dist and
Given `μz_i = mean(ij_dist)` and `μ_j = mean(ji_dist)`, these must satisfy `μ_i* length(i_to_j_contacts) == μ_j* length(j_to_i_contacts)`
"""
function generate_contact_vectors!(ij_dist,ji_dist,i_to_j_contacts::Vector{T}, j_to_i_contacts::Vector{T}) where T
- # display(ij_dist)
-
- try
- rand!(Random.default_rng(Threads.threadid()),ij_dist,i_to_j_contacts)
- catch e
- display(ij_dist.p)
- throw(e)
-
- end
- rand!(Random.default_rng(Threads.threadid()),ji_dist,j_to_i_contacts)
+ rand!(local_rng(),ij_dist,i_to_j_contacts)
+ rand!(local_rng(),ji_dist,j_to_i_contacts)
l_i = length(i_to_j_contacts)
l_j = length(j_to_i_contacts)
@@ -50,10 +42,10 @@ function generate_contact_vectors!(ij_dist,ji_dist,i_to_j_contacts::Vector{T}, j
sample_list_j = similar(sample_list_i)
while csum != 0
- sample!(Random.default_rng(Threads.threadid()),1:l_i,index_list_i)
- sample!(Random.default_rng(Threads.threadid()),1:l_j,index_list_j)
- rand!(Random.default_rng(Threads.threadid()),ij_dist,sample_list_i)
- rand!(Random.default_rng(Threads.threadid()),ji_dist,sample_list_j)
+ sample!(local_rng(),1:l_i,index_list_i)
+ sample!(local_rng(),1:l_j,index_list_j)
+ rand!(local_rng(),ij_dist,sample_list_i)
+ rand!(local_rng(),ji_dist,sample_list_j)
@inbounds 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)
diff --git a/CovidAlertVaccinationModel/src/ABM/mixing_distributions.jl b/CovidAlertVaccinationModel/src/ABM/mixing_distributions.jl
index ee092e7d02b5802ba7d1aa260b2ba297c5e4e7a9..6a700dc82836d1646e5b1d5ae50de68df43f3527 100644
--- a/CovidAlertVaccinationModel/src/ABM/mixing_distributions.jl
+++ b/CovidAlertVaccinationModel/src/ABM/mixing_distributions.jl
@@ -57,5 +57,5 @@ const unemployment_matrix = alpha_matrix(
Sample initial_workschool_mixing_matrix, which is the workschool distributions symmetrized for the full Canadian population, rather than subsets (as used in the ABM). This is used in IntervalsModel.
"""
@views function ws_sample(age)
- return rand.(Random.default_rng(Threads.threadid()),initial_workschool_mixing_matrix[age,:]) * (rand(Random.default_rng(Threads.threadid())) < (5/7))
+ return rand.(local_rng(),initial_workschool_mixing_matrix[age,:]) * (rand(local_rng()) < (5/7))
end
diff --git a/CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl b/CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl
index 86c0676d0a86b0368af424bf89d126bf78b66872..3019c906fb6a69bd93873837c725d9d59b598275 100644
--- a/CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl
+++ b/CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl
@@ -36,7 +36,7 @@ function sample_mixing_graph!(mixing_graph)
# 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!(Random.default_rng(Threads.threadid()), mixing_edges.sampler_matrix[j,i],mixing_edges.sample_cache[j,i])
+ rand!(local_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)
@@ -97,8 +97,8 @@ function create_mixing_edges(demographic_index_vectors,mixing_matrix,weights_dis
stubs_i = Vector{Int}(undef,m)
stubs_j = Vector{Int}(undef,m)
if m>0
- sample!(Random.default_rng(Threads.threadid()),demographic_index_vectors[i],Weights(num_degrees_ij./m),stubs_i)
- sample!(Random.default_rng(Threads.threadid()),demographic_index_vectors[j],Weights(num_degrees_ji./m),stubs_j)
+ sample!(local_rng(),demographic_index_vectors[i],Weights(num_degrees_ij./m),stubs_i)
+ sample!(local_rng(),demographic_index_vectors[j],Weights(num_degrees_ji./m),stubs_j)
tot += m
end
@@ -175,7 +175,7 @@ function time_dep_mixing_graphs(len,base_network,demographics,index_vectors,ws_m
if !(day_of_week == 3 || day_of_week == 4) #simulation begins on thursday I guess
push!(l, ws_static_edges)
end
- if rand(Random.default_rng(Threads.threadid()))<5/7
+ if rand(local_rng())<5/7
push!(l, ws_weekly_edges)
push!(l, rest_weekly_edges)
end
diff --git a/CovidAlertVaccinationModel/src/ABM/model_setup.jl b/CovidAlertVaccinationModel/src/ABM/model_setup.jl
index 91ea11c955fd7a0d63b9c85d56013f87028b99d3..637e6e7442b3f67daa99a4ac85cfb08dd829c6bc 100644
--- a/CovidAlertVaccinationModel/src/ABM/model_setup.jl
+++ b/CovidAlertVaccinationModel/src/ABM/model_setup.jl
@@ -3,10 +3,13 @@ function get_parameters()
params = (
sim_length = 500,
num_households = 5000,
- I_0_fraction = 0.002,
- β_y = 0.0001,
- β_m = 0.0001,
- β_o = 1.0,
+ I_0_fraction = 0.005,
+ β_y = 0.1,
+ β_m = 0.1,
+ β_o = 0.1,
+ α_y = 0.0,
+ α_m = 0.0,
+ α_o = 0.0,
recovery_rate = 1/5,
Ï€_base_y = -5.0,
Ï€_base_m = -5.0,
@@ -30,7 +33,7 @@ function get_parameters()
end
function get_u_0(nodes,I_0_fraction,vaccinator_prob)
- is_vaccinator = rand(Random.default_rng(Threads.threadid()),nodes) .< vaccinator_prob
+ is_vaccinator = rand(local_rng(),nodes) .< vaccinator_prob
status = fill(Susceptible,nodes)
return status,is_vaccinator
end
@@ -44,7 +47,7 @@ function app_users(demographics,app_usage_prob)
is_app_user = Vector{Bool}(undef,length(demographics))
@inbounds for i in eachindex(demographics)
demo = demographics[i]
- is_app_user[i] = rand(Random.default_rng(Threads.threadid())) < app_usage_prob*ymo_usage[Int(demo)]
+ is_app_user[i] = rand(local_rng()) < app_usage_prob*ymo_usage[Int(demo)]
end
return is_app_user
end
diff --git a/CovidAlertVaccinationModel/src/ABM/parameter_optimization.jl b/CovidAlertVaccinationModel/src/ABM/parameter_optimization.jl
index d425167b29f06f23e12d3f47cf0580d1095a2777..dab46fda1299b247dcfa64a521bb4753f4650832 100644
--- a/CovidAlertVaccinationModel/src/ABM/parameter_optimization.jl
+++ b/CovidAlertVaccinationModel/src/ABM/parameter_optimization.jl
@@ -86,6 +86,50 @@ function fit_post_inf_behavioural_parameters(p_tuple)
return NamedTuple{p_names}(ntuple(i -> out.P[i].particles,length(p_names)))
end
+function fit_all_parameters(p_tuple)
+ p_names = (:ω,:β_y,:β_m,:β_o,:π_base_y,:π_base_m,:π_base_o)
+ priors = Factored(
+ Uniform(0.0,0.1),
+ Uniform(0.0,0.1),
+ Uniform(0.0,0.1),
+ Uniform(0.0,1.0),
+ Uniform(-5.0,5.0),
+ Uniform(-5.0,5.0),
+ Uniform(-5.0,5.0),
+ )
+ #simulation begins in july
+ #60 days for opinion dynamics to stabilize, then immunization begins in september,
+ #infection introduced at beginning of december
+ sim_length = 300
+ p_tuple_adjust = merge(p_tuple,
+ (
+ sim_length = sim_length,
+ I_0_fraction = 0.002,
+ immunization_begin_day =60,
+ infection_introduction_day = 180,
+ immunizing = true,
+ )
+ )
+ function cost(p)
+ output,model = solve_w_parameters(p_tuple_adjust, p_names,p)
+ target_ymo_vac = ymo_vac .* sum(vaccination_data[1:end]) .* length.(model.index_vectors)
+ ymo_vaccination_ts = output.daily_immunized_by_age
+ total_postinf_vaccination = sum.(eachrow(ymo_vaccination_ts[:,180:end]))
+
+ final_size = sum.(eachrow(output.daily_cases_by_age))
+
+ target_final_size = ymo_attack_rate .* length.(model.index_vectors)
+ target_ymo_vac = ymo_vac .* sum(vaccination_data[1:4]) .* length.(model.index_vectors)
+ ymo_vaccination_ts = output.daily_immunized_by_age
+ total_preinfection_vaccination = sum.(eachrow(ymo_vaccination_ts))
+ return 1e-1*sum((total_preinfection_vaccination .- target_ymo_vac).^2) + 1e-1*sum((total_postinf_vaccination .- target_ymo_vac).^2) + sum((final_size .- target_final_size).^2)
+ end
+
+ # display(cost([0.000,0.001,0.001,1.0]))
+ out =smc(priors,cost; verbose = true, nparticles = 800, parallel = true)# ABCDE(priors,cost,1e6; verbose=true, nparticles=300,generations=1000, parallel = true) #this one has better NaN handling
+ return NamedTuple{p_names}(ntuple(i -> out.P[i].particles,length(p_names)))
+end
+
function plot_behavioural_fit(particles,p_tuple)
p_names = (:Ï€_base_y,:Ï€_base_m,:Ï€_base_o)
sim_length = 210
@@ -122,25 +166,32 @@ end
function fit_parameters(default_parameters)
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")
+ fit_all_parameters_path = joinpath(PACKAGE_FOLDER,"abm_parameter_fits","fit_all_parameters.dat")
# pre_inf_behaviour_parameters = fit_pre_inf_behavioural_parameters(default_parameters)
# serialize(pre_inf_behaviour_parameters_path, pre_inf_behaviour_parameters)
- pre_inf_behaviour_parameters = deserialize(pre_inf_behaviour_parameters_path)
- display(map(mode,pre_inf_behaviour_parameters))
- 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)
+ # pre_inf_behaviour_parameters = deserialize(pre_inf_behaviour_parameters_path)
+ # display(map(mode,pre_inf_behaviour_parameters))
+ # 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(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)
+
+ # fitted_parameter_tuple = (;
+ # deserialize(pre_inf_behaviour_parameters_path)...,
+ # deserialize(post_inf_behaviour_parameters_path)...
+ # )
+ # display(map(mode,fitted_parameters_with_post_inf_behaviour))
+
+ output = fit_all_parameters(default_parameters)
+ serialize(fit_all_parameters_path,output)
+
+
+ fitted_parameter_tuple = deserialize(fit_all_parameters_path)
+ fitted_sol = plot_fitting_posteriors("post_inf_fitting",fitted_parameter_tuple,default_parameters)
return fitted_sol
end
@@ -190,20 +241,18 @@ function plot_fitting_posteriors(fname,particles_tuple,parameters)
return out
end
-#error function for inf
-#
-function visualize_Ï€_base(particles_tuple)
+# function visualize_Ï€_base(particles_tuple)
- param_keys = [:Ï€_base_y,:Ï€_base_m,:Ï€_base_o]
- # π_bases_array = map(f -> getproperty(particles_tuple,f),param_keys) |> l -> mapreduce(t-> [t...],hcat,zip(l...))
- # display(Ï€_bases_array)
+# param_keys = [:Ï€_base_y,:Ï€_base_m,:Ï€_base_o]
+# # π_bases_array = map(f -> getproperty(particles_tuple,f),param_keys) |> l -> mapreduce(t-> [t...],hcat,zip(l...))
+# # display(Ï€_bases_array)
- # p = cornerplot(Ï€_bases_array'; labels = string.(param_keys))
- params = NamedTuple{(param_keys...,)}(particles_tuple)
- # display(params)
- p = corner(params)
- display(p)
+# # p = cornerplot(Ï€_bases_array'; labels = string.(param_keys))
+# params = NamedTuple{(param_keys...,)}(particles_tuple)
+# # display(params)
+# p = corner(params)
+# display(p)
- # display(scatter(map(f -> getproperty(particles_tuple,f),param_keys)...))
-end
\ No newline at end of file
+# # display(scatter(map(f -> getproperty(particles_tuple,f),param_keys)...))
+# end
\ No newline at end of file
diff --git a/CovidAlertVaccinationModel/src/ABM/solve.jl b/CovidAlertVaccinationModel/src/ABM/solve.jl
index 2366d3eed085c0093d49bbb29d2ecb49107251b0..22a0e62f771f1686a2c2f04c68a453138e24958f 100644
--- a/CovidAlertVaccinationModel/src/ABM/solve.jl
+++ b/CovidAlertVaccinationModel/src/ABM/solve.jl
@@ -27,7 +27,7 @@ Base.@propagate_inbounds @views function update_alert_durations!(t,modelsol) # B
end
end
coin_flip = 1 - (1 - notification_parameter)^total_weight_i
- r = rand(Random.default_rng(Threads.threadid()))
+ r = rand(local_rng())
if r < coin_flip
covid_alert_notifications[end,i] = 1 #add the notifications for today
else
@@ -40,7 +40,7 @@ Base.@propagate_inbounds @views function update_alert_durations!(t,modelsol) # B
end
Base.@propagate_inbounds @views function update_infection_state!(t,modelsol)
- @unpack β_y,β_m,β_o,recovery_rate,immunizing,immunization_begin_day = modelsol.params
+ @unpack β_y,β_m,β_o,α_y,α_m,α_o,recovery_rate,immunizing,immunization_begin_day = modelsol.params
@unpack u_inf,u_vac,u_next_inf,u_next_vac,demographics,inf_network,status_totals, immunization_countdown = modelsol
modelsol.daily_cases_by_age .= 0
@@ -56,16 +56,20 @@ Base.@propagate_inbounds @views function update_infection_state!(t,modelsol)
agent_status = u_inf[i]
is_vaccinator = u_vac[i]
agent_demo = demographics[i]
- if agent_status == Susceptible
- if is_vaccinator && immunizing && immunization_countdown[i] == -1 && t> immunization_begin_day
- immunization_countdown[i] = 14
- else
- for mixing_graph in inf_network.graph_list[t]
- for j in neighbors(mixing_graph,i)
- if u_inf[j] == Infected && u_next_inf[i] != Infected
-
- β_vec = (β_y,β_m,β_o)
- if rand(Random.default_rng(Threads.threadid())) < contact_weight(β_vec[Int(agent_demo)],get_weight(mixing_graph,GraphEdge(i,j)))
+ if agent_status == Susceptible && is_vaccinator && immunizing && immunization_countdown[i] == -1 && t> immunization_begin_day
+ immunization_countdown[i] = 14
+ end
+ if agent_status == Susceptible || agent_status == Immunized
+ for mixing_graph in inf_network.graph_list[t]
+ for j in neighbors(mixing_graph,i)
+ if u_inf[j] == Infected && u_next_inf[i] != Infected
+ β_vec = (β_y,β_m,β_o)
+ α_vec = (α_y,α_m,α_o)
+ if rand(local_rng()) < contact_weight(β_vec[Int(agent_demo)],get_weight(mixing_graph,GraphEdge(i,j)))
+ if agent_status == Immunized && rand(local_rng()) < 1- α_vec[Int(agent_demo)]
+ modelsol.daily_cases_by_age[Int(agent_demo)]+=1
+ agent_transition!(i, Immunized,Infected)
+ elseif agent_status == Susceptible
modelsol.daily_cases_by_age[Int(agent_demo)]+=1
agent_transition!(i, Susceptible,Infected)
end
@@ -74,7 +78,7 @@ Base.@propagate_inbounds @views function update_infection_state!(t,modelsol)
end
end
elseif agent_status == Infected
- if rand(Random.default_rng(Threads.threadid())) < recovery_rate
+ if rand(local_rng()) < recovery_rate
agent_transition!(i, Infected,Recovered)
end
end
@@ -97,10 +101,10 @@ Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,mod
Ï€_base = @SVector [Ï€_base_y,Ï€_base_m,Ï€_base_o]
vac_payoff = 0
num_soc_nbrs = 0
- random_soc_network = sample(Random.default_rng(Threads.threadid()), soc_network.graph_list[t])
+ random_soc_network = sample(local_rng(), soc_network.graph_list[t])
if !isempty(neighbors(random_soc_network,i))
- random_neighbour = sample(Random.default_rng(Threads.threadid()), neighbors(random_soc_network.g,i))
+ random_neighbour = sample(local_rng(), neighbors(random_soc_network.g,i))
if u_vac[random_neighbour] == u_vac[i]
vac_payoff += π_base[Int(demographics[i])] + total_infections*ω
if app_user[i] && time_of_last_alert[app_user_list[i]]>=0
@@ -108,11 +112,11 @@ Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,mod
end
if u_vac[i]
- if rand(Random.default_rng(Threads.threadid())) < 1 - Φ(vac_payoff,ξ)
+ if rand(local_rng()) < 1 - Φ(vac_payoff,ξ)
u_next_vac[i] = false
end
else
- if rand(Random.default_rng(Threads.threadid())) < Φ(vac_payoff,ξ)
+ if rand(local_rng()) < Φ(vac_payoff,ξ)
u_next_vac[i] = true
end
end
@@ -144,7 +148,7 @@ function agents_step!(t,modelsol)
end
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))
+ init_indices = rand(local_rng(), 1:modelsol.nodes, round(Int,modelsol.nodes*modelsol.params.I_0_fraction))
modelsol.u_inf[init_indices] .= Infected
modelsol.status_totals[Int(Infected)] += length(init_indices)
end
diff --git a/CovidAlertVaccinationModel/src/IntervalsModel/hh_durations_model.jl b/CovidAlertVaccinationModel/src/IntervalsModel/hh_durations_model.jl
index 9704af083d6228315f4e7a221f0519cdf2c73d20..5e78911c341c91f80da9c328d1f1a65e9822d19a 100644
--- a/CovidAlertVaccinationModel/src/IntervalsModel/hh_durations_model.jl
+++ b/CovidAlertVaccinationModel/src/IntervalsModel/hh_durations_model.jl
@@ -42,7 +42,7 @@ function err_hh(p,dist)
@inbounds for age_j in YOUNG:OLD #for a given age_sample loop over possible contact ages
#get num_contacts from the HH data
- durs = trunc.(Int,rand(Random.default_rng(Threads.threadid()),age_dists[age_sample,age_j],num_contacts_subarray[i,age_j])) .% durmax
+ durs = trunc.(Int,rand(local_rng(),age_dists[age_sample,age_j],num_contacts_subarray[i,age_j])) .% durmax
#this returns the MEAN of the distribution of total durations, given we have durations given by durs, and we sample comparison_samples from that distribution.
expdur = tot_dur_sample(comparison_samples,durs)
diff --git a/CovidAlertVaccinationModel/src/IntervalsModel/interval_overlap_sampling.jl b/CovidAlertVaccinationModel/src/IntervalsModel/interval_overlap_sampling.jl
index c5e81477b3c05d233dad63f12712a901de64a848..d9a4db41fa515f292263efe4c0f6e8db7c921943 100644
--- a/CovidAlertVaccinationModel/src/IntervalsModel/interval_overlap_sampling.jl
+++ b/CovidAlertVaccinationModel/src/IntervalsModel/interval_overlap_sampling.jl
@@ -27,7 +27,7 @@ function tot_dur_sample(n,durlist)
int_list = Vector{Interval{Int,Closed,Closed}}()
sizehint!(int_list,numcontact*2)
- start_matrix = trunc.(Int,rand(Random.default_rng(Threads.threadid()),startdist,(numcontact,n)))
+ start_matrix = trunc.(Int,rand(local_rng(),startdist,(numcontact,n)))
@inbounds for i in 1:n
empty!(int_list)
@inbounds for j in 1:numcontact
@@ -59,7 +59,7 @@ function tot_dur_sample!(sample_list,durlist)
n = length(sample_list)
int_list = Vector{Interval{Int,Closed,Closed}}()
sizehint!(int_list,numcontact*2)
- start_matrix = trunc.(Int,rand(Random.default_rng(Threads.threadid()),startdist,(numcontact,n)))
+ start_matrix = trunc.(Int,rand(local_rng(),startdist,(numcontact,n)))
for i in 1:n
empty!(int_list)
for j in 1:numcontact
diff --git a/CovidAlertVaccinationModel/src/IntervalsModel/rest_durations_model.jl b/CovidAlertVaccinationModel/src/IntervalsModel/rest_durations_model.jl
index 61100a7027cc35936eb1eec127c8c611aac6b7a1..0d1c9cd4491ae0bf813c3dbfde2a9b8541d22d6b 100644
--- a/CovidAlertVaccinationModel/src/IntervalsModel/rest_durations_model.jl
+++ b/CovidAlertVaccinationModel/src/IntervalsModel/rest_durations_model.jl
@@ -35,11 +35,11 @@ function err_rest(p,dist)
@inbounds for _ in 1:sample_repeat
#sample a matrix of neighourhood sizes for distributions in rest_distributions
- neighourhoods = rand.(Random.default_rng(Threads.threadid()),rest_distributions)
+ neighourhoods = rand.(local_rng(),rest_distributions)
@inbounds for age_sample in YOUNG:OLD, age_j in YOUNG:OLD #for a given age_sample loop over possible contact ages
if neighourhoods[age_sample,age_j] > 0
#get durations from our candidate distribtions for each of the contacts in neighbourhood
- durs = trunc.(Int,rand(Random.default_rng(Threads.threadid()),age_dists[age_sample,age_j],neighourhoods[age_sample,age_j])) .% durmax
+ durs = trunc.(Int,rand(local_rng(),age_dists[age_sample,age_j],neighourhoods[age_sample,age_j])) .% durmax
#this MODIFIES sample_list to contain samples from the distribution of total_durations, given intervals of length dur.
tot_dur_sample!(sample_list,durs)
diff --git a/CovidAlertVaccinationModel/src/IntervalsModel/ws_durations_model.jl b/CovidAlertVaccinationModel/src/IntervalsModel/ws_durations_model.jl
index 794de2e59ed49aa90c617192a0a9446248b0d0db..0e30d32f3203e10b35881cdd5c58924c4e70c99d 100644
--- a/CovidAlertVaccinationModel/src/IntervalsModel/ws_durations_model.jl
+++ b/CovidAlertVaccinationModel/src/IntervalsModel/ws_durations_model.jl
@@ -39,7 +39,7 @@ function err_ws(p,dist)
for age_j in YOUNG:OLD #for a given age_sample loop over possible contact ages
if neighourhoods[age_j] > 0
#get durations from our candidate distribtions for each of the contacts in neighbourhood
- durs = trunc.(Int,rand(Random.default_rng(Threads.threadid()),age_dists[age_sample,age_j],neighourhoods[age_j])) .% durmax
+ durs = trunc.(Int,rand(local_rng(),age_dists[age_sample,age_j],neighourhoods[age_j])) .% durmax
#this MODIFIES sample_list to contain samples from the distribution of total_durations, given intervals of length dur.
tot_dur_sample!(sample_list,durs)
diff --git a/CovidAlertVaccinationModel/src/data.jl b/CovidAlertVaccinationModel/src/data.jl
index e238bf4ae1ff90a6da3bcc010b42d004ff9297ea..03fba52571cec40ad224eaab875d7a185babfc25 100644
--- a/CovidAlertVaccinationModel/src/data.jl
+++ b/CovidAlertVaccinationModel/src/data.jl
@@ -44,11 +44,11 @@ function read_household_data()
end
function sample_household_data(n)
- return sample(Random.default_rng(Threads.threadid()),household_data.households, Weights(household_data.weight_vector),n)
+ return sample(local_rng(),household_data.households, n)
end
function get_household_data_proportions()
- households_by_demographic_sum = sum.([map(((l,w),)-> l[i]*w,zip(household_data.households,household_data.weight_vector)) for i in 1:3])
+ households_by_demographic_sum = sum.([map(l -> l[i],household_data.households) for i in 1:3])
return households_by_demographic_sum./sum(households_by_demographic_sum)
# https://www.publichealthontario.ca/-/media/documents/ncov/epi/covid-19-severe-outcomes-ontario-epi-summary.pdf?la=en
end
diff --git a/post_inf_fitting.pdf b/post_inf_fitting.pdf
index 142fdb59fb513b713a13f1459af0a6bae237de1f..bb385cb9da521998480fb34f201f0d7384af44a5 100644
Binary files a/post_inf_fitting.pdf and b/post_inf_fitting.pdf differ
diff --git a/post_inf_fitting_posteriors.pdf b/post_inf_fitting_posteriors.pdf
index 29cf52fb7dea91c2cdd8f2253be437f00f799a86..f229d15e9fe61d08d9229d706f8ff9ac027dbe49 100644
Binary files a/post_inf_fitting_posteriors.pdf and b/post_inf_fitting_posteriors.pdf differ
diff --git a/timeseries.pdf b/timeseries.pdf
index 024049c519625a870b168c82773773c9a8b7be6e..a259a5eee94d30a346c39f883344a11108b61179 100644
Binary files a/timeseries.pdf and b/timeseries.pdf differ