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

new simulations with all parameters at once, reinstane VectorizedRNG, remove...

new simulations with all parameters at once, reinstane VectorizedRNG, remove weights for household sampling, need to redo intervalsmodel again
parent c241302e
......@@ -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"
......
......@@ -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)
......
......@@ -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
......@@ -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
......
......@@ -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
......
......@@ -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
......@@ -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
......
......@@ -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)
......
......@@ -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
......
......@@ -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)
......
......@@ -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)
......
......@@ -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
......
No preview for this file type
No preview for this file type
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