Commit 800d7af5 authored by Peter Jentsch's avatar Peter Jentsch
Browse files

found bug in VectorizedRNG, back to default. Also reran intervals model, and...

found bug in VectorizedRNG, back to default. Also reran intervals model, and now should rerun model fitting.
parent f1adfdc0
......@@ -39,9 +39,9 @@ version = "0.1.0"
[[ArrayInterface]]
deps = ["IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "b08be763d0b8ddee6b162016dad746a69980616d"
git-tree-sha1 = "b09fe16aa9dc587cccce838e6cb6d6e1f4831d7f"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "3.1.11"
version = "3.1.12"
[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
......
......@@ -11,4 +11,8 @@ function solve_and_plot_parameters()
end
out = solve_and_plot_parameters()
println(sum.(eachrow(mean.(out.daily_immunized_by_age))))
\ No newline at end of file
println(sum.(eachrow(mean.(out.daily_immunized_by_age))))
final_size = sum.(eachrow(mean.(out.daily_cases_by_age)))
display(final_size)
\ No newline at end of file
......@@ -43,9 +43,10 @@ function abm(parameters, recorder)
model_sol = ModelSolution(parameters.sim_length,parameters,5000)
output = solve!(model_sol,recorder )
total_weighted_degree = map(modelsol.index_vectors) do age_group_indices
return mean(map(i -> weighted_degree(i,modelsol.inf_network),age_group_indices))
end
@show total_weighted_degree
# total_weighted_degree = map(modelsol.index_vectors) do age_group_indices
# return mean(map(i -> weighted_degree(i,modelsol.inf_network),age_group_indices))
# end
# @show total_weighted_degree
return model_sol
end
......@@ -27,8 +27,17 @@ 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
rand!(local_rng(),ij_dist,i_to_j_contacts)
rand!(local_rng(),ji_dist,j_to_i_contacts)
# 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)
l_i = length(i_to_j_contacts)
l_j = length(j_to_i_contacts)
......@@ -41,10 +50,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!(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)
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)
@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.(initial_workschool_mixing_matrix[age,:]) * (rand(local_rng()) < (5/7))
return rand.(Random.default_rng(Threads.threadid()),initial_workschool_mixing_matrix[age,:]) * (rand(Random.default_rng(Threads.threadid())) < (5/7))
end
......@@ -25,7 +25,6 @@ function Base.isequal(e1::GraphEdge,e2::GraphEdge)
return isequal(minmax(e1.a,e1.b),minmax(e2.a,e2.b))
end
"""
sample_mixing_graph!(mixing_graph::Graph)
......@@ -37,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!(local_rng(), mixing_edges.sampler_matrix[j,i],mixing_edges.sample_cache[j,i])
rand!(Random.default_rng(Threads.threadid()), 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)
......@@ -78,7 +77,7 @@ struct MixingEdges{M}
sampler_matrix::M
function MixingEdges(total_edges,contact_array,sampler_matrix::Matrix{M}) where M<:Sampleable{Univariate, Discrete}
sample_cache = map(v-> Vector{Int}(undef,length(v)),contact_array)
weights_dict = Dictionary{GraphEdge,UInt8}()
weights_dict = Dictionary{GraphEdge,UInt8}(;sizehint = total_edges)
new{typeof(sampler_matrix)}(total_edges,contact_array,sample_cache,weights_dict,sampler_matrix)
end
end
......@@ -98,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!(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)
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)
tot += m
end
......@@ -176,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(local_rng())<5/7
if rand(Random.default_rng(Threads.threadid()))<5/7
push!(l, ws_weekly_edges)
push!(l, rest_weekly_edges)
end
......
......@@ -4,17 +4,19 @@ function get_parameters()
sim_length = 500,
num_households = 5000,
I_0_fraction = 0.002,
base_transmission_probability = 0.001,
β_y = 0.0001,
β_m = 0.0001,
β_o = 1.0,
recovery_rate = 1/5,
π_base_y = -5.0,
π_base_m = -5.0,
π_base_o = 1.0,
π_base_o = 0.0,
η = 0.0,
κ = 0.0,
ω = 0.00005,
ω_en = 0.00,
γ = 0.0,
β = 5.0,
ξ = 5.0,
notification_parameter = 0.001,
vaccinator_prob = 0.2,
app_user_fraction = 0.0,
......@@ -28,7 +30,7 @@ function get_parameters()
end
function get_u_0(nodes,I_0_fraction,vaccinator_prob)
is_vaccinator = rand(local_rng(),nodes) .< vaccinator_prob
is_vaccinator = rand(Random.default_rng(Threads.threadid()),nodes) .< vaccinator_prob
status = fill(Susceptible,nodes)
return status,is_vaccinator
end
......@@ -42,7 +44,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(local_rng()) < app_usage_prob*ymo_usage[Int(demo)]
is_app_user[i] = rand(Random.default_rng(Threads.threadid())) < app_usage_prob*ymo_usage[Int(demo)]
end
return is_app_user
end
......
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 = [10.376,5.636,7.2]
const ymo_attack_rate = [10.376,5.636,7.2]./100
function solve_w_parameters(default_p, p_names, new_p_list)
......@@ -39,8 +38,7 @@ function fit_pre_inf_behavioural_parameters(p_tuple)
)
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
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))
# display(target_ymo_vac)
......@@ -54,8 +52,8 @@ end
function fit_post_inf_behavioural_parameters(p_tuple)
samples = 1
p_names = (:ω,:base_transmission_probability)
priors = Factored(Uniform(0.0,0.1),Uniform(0.0,0.1))
p_names = (:ω,:β_y,:β_m,:β_o)
priors = Factored(Uniform(0.0,0.1),Uniform(0.0,0.1),Uniform(0.0,0.1),Uniform(0.0,1.0))
#simulation begins in july
#60 days for opinion dynamics to stabilize, then immunization begins in september,
#infection introduced at beginning of december
......@@ -74,20 +72,22 @@ function fit_post_inf_behavioural_parameters(p_tuple)
target_cumulative_vaccinations = target_cumulative_vac_proportion*model.nodes
target_ymo_vac = ymo_vac .* sum(vaccination_data[5:end]) .* target_cumulative_vaccinations
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 .* 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)
target_final_size = ymo_attack_rate .* length.(model.index_vectors)
# display( length.(model.index_vectors))
# display((final_size,target_final_size))
# display((total_postinf_vaccination,target_ymo_vac))
# display((1e-1*sum((total_postinf_vaccination .- target_ymo_vac).^2) , sum((final_size .- target_final_size).^2)))
return 1e-1*sum((total_postinf_vaccination .- target_ymo_vac).^2) + sum((final_size .- target_final_size).^2)
end
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
# 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}((out.P.particles,))
end
......@@ -129,8 +129,8 @@ function fit_parameters(default_parameters)
post_inf_behaviour_parameters_path = joinpath(PACKAGE_FOLDER,"abm_parameter_fits","post_inf_behaviour_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 = 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)
......
......@@ -3,12 +3,12 @@
# 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
function contact_weight(β, contact_time)
return 1 - (1-β)^contact_time
end
function Φ(payoff,β)
return 1 / (exp(-1*β*payoff))
function Φ(payoff,ξ)
return 1 / (exp(-1*ξ*payoff))
end
Base.@propagate_inbounds @views function update_alert_durations!(t,modelsol) # Base.@propagate_inbounds
......@@ -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(local_rng())
r = rand(Random.default_rng(Threads.threadid()))
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 base_transmission_probability,recovery_rate,immunizing,immunization_begin_day = modelsol.params
@unpack β_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
......@@ -62,8 +62,10 @@ Base.@propagate_inbounds @views function update_infection_state!(t,modelsol)
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
if rand(local_rng()) < contact_weight(base_transmission_probability,get_weight(mixing_graph,GraphEdge(i,j)))
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)))
modelsol.daily_cases_by_age[Int(agent_demo)]+=1
agent_transition!(i, Susceptible,Infected)
end
......@@ -72,7 +74,7 @@ Base.@propagate_inbounds @views function update_infection_state!(t,modelsol)
end
end
elseif agent_status == Infected
if rand(local_rng()) < recovery_rate
if rand(Random.default_rng(Threads.threadid())) < recovery_rate
agent_transition!(i, Infected,Recovered)
end
end
......@@ -88,17 +90,17 @@ end
Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,modelsol,total_infections)
@unpack π_base_y,π_base_m,π_base_o, η,γ, κ, ω, ω_en,γ,β = modelsol.params
@unpack π_base_y,π_base_m,π_base_o, η,γ, κ, ω, ω_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
π_base = @SVector [π_base_y,π_base_m,π_base_o]
vac_payoff = 0
num_soc_nbrs = 0
random_soc_network = sample(local_rng(), soc_network.graph_list[t])
random_soc_network = sample(Random.default_rng(Threads.threadid()), soc_network.graph_list[t])
if !isempty(neighbors(random_soc_network,i))
random_neighbour = sample(local_rng(), neighbors(random_soc_network.g,i))
random_neighbour = sample(Random.default_rng(Threads.threadid()), 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
......@@ -106,11 +108,11 @@ Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,mod
end
if u_vac[i]
if rand(local_rng()) < 1 - Φ(vac_payoff,β)
if rand(Random.default_rng(Threads.threadid())) < 1 - Φ(vac_payoff,ξ)
u_next_vac[i] = false
end
else
if rand(local_rng()) < Φ(vac_payoff,β)
if rand(Random.default_rng(Threads.threadid())) < Φ(vac_payoff,ξ)
u_next_vac[i] = true
end
end
......@@ -142,7 +144,7 @@ function agents_step!(t,modelsol)
end
if t == modelsol.params.infection_introduction_day
init_indices = rand(local_rng(), 1:modelsol.nodes, round(Int,modelsol.nodes*modelsol.params.I_0_fraction))
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
modelsol.status_totals[Int(Infected)] += length(init_indices)
end
......
......@@ -62,7 +62,7 @@ function bayesian_estimation(fname, err_func, priors_list, dists, particles; alp
data_pairs = map(zip(dists,priors_list)) do (dist,priors)
init = rand(priors)
# @btime $err_func($init,$dist)
out = smc(priors,p -> err_func(p, dist), verbose=true, nparticles=particles, alpha=alpha, parallel = true)
out = smc(priors,p -> err_func(p, dist), verbose=true, nparticles=particles, parallel = true)
return string(dist) => out
end |> OrderedDict
......
#testing some GPU implementations of the ABM model
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
end
function agents_step_cuda!(t,u_next,u,population_list,graph,params,index_vectors,vaccination_algorithm!)
@unpack p, recovery_rate, vaccines_per_day,vaccination_start_day = params
u_next .= u#keep state the same if nothing else happens
if t >= vaccination_start_day
vaccination_algorithm!(Int(floor(vaccines_per_day*length(u))),u_next,u,graph,population_list,index_vectors)
end
@inbounds for i in 1:length(u)
agent_status = u[i]
agent_demo = population_list[i]
if agent_status == Susceptible
@inbounds for j in 1:length(u)
if graph[i,j] != 0 && u[j] == Infected
if rand(RNG) < contact_weight(p,graph[i,j])
(u_next[i] = Infected)
end
end
end
elseif agent_status == Infected
if rand(RNG) < recovery_rate
u_next[i] = Recovered
end
end
end
end
function agents_step_gpu!(t,u_next,u,graph,params)
u_next .= u
function kernel(state, _, (t,u_next,u,g,params,randstate))
@unpack p, recovery_rate = params
i = linear_index(state)
if i <= length(u)
agent_status = u[i]
if agent_status == Susceptible
for j in 1:length(u)
if g[i,j] == 1 && u[j] == Infected && GPUArrays.gpu_rand(Float64, state, randstate) < p
u_next[i] = Infected
end
end
elseif agent_status == Infected
if GPUArrays.gpu_rand(Float64, state, randstate) < recovery_rate
u_next[i] = Recovered
end
end
end
return nothing
end
gpu_call(kernel, u, (t,u_next,u,graph,params, GPUArrays.default_rng(typeof(u)).state))
end
function solve_cuda!(u_0,params,steps,agent_model,vaccination_algorithm!; record = false)
solution = vcat([u_0], [fill(Susceptible,length(u_0)) for i in 1:steps])
population_list = cu(agent_model.demographics) #list of demographic status for each agent
index_vectors = cu.(agent_model.demographic_index_vectors) #lists of indices of each agent with a given demographic
base_network = cu(agent_model.base_network) #network with households and LTC homes
ws_contacts = cu(agent_model.ws_matrix_list)
rest_contacts = cu(agent_model.rest_matrix_list)
ws_daily_graph = MixingGraph(index_vectors,ws_contacts.daily,contact_time_distribution_matrix)
rest_daily_graph = MixingGraph(index_vectors,rest_contacts.daily,contact_time_distribution_matrix)
ws_weekly_graph = MixingGraph(index_vectors,ws_contacts.twice_a_week,contact_time_distribution_matrix)
rest_weekly_graph = MixingGraph(index_vectors,rest_contacts.twice_a_week,contact_time_distribution_matrix)
graph = similar(base_network)
graphlist = typeof(graph)[]
for t in 1:steps
day_of_week = mod(t,7)
graph .= 0
graph .= base_network
if !(day_of_week == 3 || day_of_week == 4)
apply_mixing_graph(graph, ws_daily_graph)
end
apply_mixing_graph(graph, rest_daily_graph)
if rand(RNG)<5/7
apply_mixing_graph(graph, ws_weekly_graph)
apply_mixing_graph(graph, rest_weekly_graph)
end
ws_once_graph = MixingGraph(index_vectors,ws_contacts.otherwise,contact_time_distribution_matrix)
rest_once_graph = MixingGraph(index_vectors,rest_contacts.otherwise,contact_time_distribution_matrix)
apply_mixing_graph(graph, ws_once_graph)
apply_mixing_graph(graph, rest_once_graph)
if record
push!(graphlist, deepcopy(graph))
end
agents_step!(t,solution[t+1],solution[t],population_list,graph,params,index_vectors,vaccination_algorithm!)
# advance agent states based on the new network, vaccination process given by vaccination_algorithm!, which is just a function defined as above
end
return solution, graphlist
end
\ No newline at end of file
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