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

multivariate sims worked

parent 0254b735
......@@ -6,7 +6,7 @@ using CovidAlertVaccinationModel:vaccination_data,ymo_vac,ymo_attack_rate
const samples = 25
function solve_and_plot_parameters()
p = CovidAlertVaccinationModel.get_parameters()
p = CovidAlertVaccinationModel.get_app_parameters()
display(p)
out,avg_populations = mean_solve(samples, p,DebugRecorder)
p = plot_model(nothing,[nothing],[out],p.infection_introduction_day,p.immunization_begin_day)
......
......@@ -13,10 +13,10 @@ function get_parameters()#(0.0000,0.00048,0.0005,0.16,-1.30,-1.24,-0.8,0.35,0.35
π_base_y = -1.37,
π_base_m = -1.46,
π_base_o = -0.95,
η = 0.0,
η = 0.5,
κ = 0.0,
ω = 0.0055,
ω_en = 0.0,
ω_en = 0.01,
Γ = 1/7,
ξ = 5.0,
notification_parameter = 0.0005,
......@@ -78,6 +78,8 @@ mutable struct ModelSolution{T,InfNet,SocNet,WSMixingDist,RestMixingDist}
daily_cases_by_age::Vector{Int}
daily_unvac_cases_by_age::Vector{Int}
daily_immunized_by_age::Vector{Int}
daily_total_notifications::Int
daily_total_notified_agents::Int
avg_weighted_degree_of_vaccinators::Variance{Float64,Float64,EqualWeight}
avg_weighted_degree::Variance{Float64,Float64,EqualWeight}
ws_matrix_tuple::WSMixingDist
......@@ -133,6 +135,8 @@ mutable struct ModelSolution{T,InfNet,SocNet,WSMixingDist,RestMixingDist}
[0,0,0],
[0,0,0],
[0,0,0],
0,
0,
Variance(),
Variance(),
ws_matrix_tuple,
......
......@@ -19,6 +19,8 @@ struct DebugRecorder{ElType,ArrT1<:AbstractArray{ElType},ArrT2<:AbstractVector{E
recorded_status_totals::ArrT1
daily_cases_by_age::ArrT3
total_vaccinators::ArrT2
daily_total_notifications::ArrT2
daily_total_notified_agents::ArrT2
unvac_final_size_by_age::ArrT2
total_postinf_vaccination::ArrT2
total_preinf_vaccination::ArrT2
......@@ -46,6 +48,10 @@ function DebugRecorder(val::T,sim_length) where T
totals = [copy(val) for i in 1:4, j in 1:sim_length]
state_totals = @LArray totals (S = (1,:),I = (2,:),R = (3,:), V = (4,:))
total_vaccinators = [copy(val) for j in 1:sim_length]
daily_total_notifications = [copy(val) for j in 1:sim_length]
daily_total_notified_agents = [copy(val) for j in 1:sim_length]
mean_time_since_last_notification = [copy(val) for j in 1:sim_length]
totals_ymo = [copy(val) for i in 1:3, j in 1:sim_length]
daily_cases_by_age = @LArray deepcopy(totals_ymo) (Y = (1,:),M = (2,:),O = (3,:))
......@@ -59,6 +65,8 @@ function DebugRecorder(val::T,sim_length) where T
state_totals,
daily_cases_by_age,
total_vaccinators,
daily_total_notifications,
daily_total_notified_agents,
unvac_final_size_by_age,
total_postinf_vaccination,
total_preinf_vaccination,
......@@ -88,6 +96,9 @@ function record!(t,modelsol, recorder::DebugRecorder)
recorder.total_vaccinators[t] = modelsol.daily_vaccinators
recorder.daily_cases_by_age[:,t] .= modelsol.daily_cases_by_age
recorder.recorded_status_totals[:,t] .= modelsol.status_totals
recorder.daily_total_notifications[t] = modelsol.daily_total_notifications
recorder.daily_total_notified_agents[t] = modelsol.daily_total_notified_agents
alerts = filter(>(0),modelsol.time_of_last_alert)
if !isempty(alerts)
mean_alert_time = mean(t .- alerts)
......@@ -126,7 +137,7 @@ end
function mean_solve(samples,parameter_tuple,recorder;progmeter = nothing)
stat_recorder = recorder(Variance(), parameter_tuple.sim_length)
avg_populations = [0.0,0.0,0.0]
sol_list = ThreadsX.map(1:samples) do _
sol_list = map(1:samples) do _
output_recorder = recorder(0.0,parameter_tuple.sim_length)
sol = abm(parameter_tuple,output_recorder)
isnothing(progmeter) || next!(progmeter)
......@@ -143,7 +154,7 @@ end
function mean_solve(samples,parameter_tuple,recorder::Type{HeatmapRecorder};progmeter = nothing)
stat_recorder = recorder(parameter_tuple.sim_length)
avg_populations = [0.0,0.0,0.0]
sol_list = ThreadsX.map(1:samples) do _
sol_list = map(1:samples) do _
output_recorder = recorder(parameter_tuple.sim_length)
sol = abm(parameter_tuple,output_recorder)
isnothing(progmeter) || next!(progmeter)
......
const samples = 50
const univariate_path = "CovidAlertVaccinationModel/plots/univariate/"
const bivariate_path = "CovidAlertVaccinationModel/plots/univariate/"
function univarate_test(variable, variable_range,samples; progmeter = nothing)
function univarate_test(variable, variable_range; progmeter = nothing)
default_parameters = get_app_parameters()
parameter_range_list = [merge(default_parameters,NamedTuple{(variable,)}((value,))) for value in variable_range]
solve_fn(p) = mean_solve(samples, p,DebugRecorder;progmeter)[1]
......@@ -17,23 +17,22 @@ if !ispath(univariate_path)
mkdir(univariate_path)
end
function univariate_simulations()
len = 5
samples = 10
len = 8
univarate_test_list = (
# (:I_0_fraction, range(0.0, 0.05; length = len)),
# (:base_transmission_probability, range(0.0002, 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(-4.5, -3.5; length = len)),
(:η, range(0.0, 0.005; length = len)),
(:η, range(0.0,2.0; length = len)),
# (:κ, range(0.5, 1.5; length = len)),
# (:ω, range(0.0, 0.01; length = len)),
(:ω_en, range(0.0, 0.001; length = len)),
(:ω_en, range(0.0, 1e-1; length = len)),
# (:γ, range(0.0, 0.5; length = len)),
# (:ξ, range(1, 15; length = len)),
(:notification_parameter, range(0.000, 0.002; length = len)),
(:app_user_fraction, range(0.05, 0.8; length = len)),
(:notification_threshold, (1:len)),
(:notification_parameter, range(0.000, 0.001; length = len)),
# (:app_user_fraction, range(0.05, 0.8; length = len)),
# (:notification_threshold, (1:len)),
# (:immunization_delay, [7,10,14,20]),
)
......@@ -41,8 +40,9 @@ function univariate_simulations()
numsim = sum(map(t -> length(t[2]), univarate_test_list))
display(numsim)
progmeter = Progress(numsim*samples)
display(get_app_parameters())
plt_list = ThreadsX.map(univarate_test_list) do ur
out = univarate_test(ur...,samples;progmeter)
out = univarate_test(ur...;progmeter)
return out
end
for ((varname,_),p) in zip(univarate_test_list,plt_list)
......@@ -52,14 +52,13 @@ end
using AxisKeys
function multivariate_simulations()
len = 10
samples = 25
len = 20
app_simulations = (
(:η, range(0.0, 0.005; length = len)),
(:ω_en, range(0.0, 0.05; length = len)),
(:η, range(0.0, 2.0; length = len)),
(:ω_en, range(0.0, 1e-1; length = len)),
# (:notification_threshold, (1:len)),
)
run_multivariate_sims(app_simulations,samples)
run_multivariate_sims(app_simulations)
# for ((varname,_),p) in zip(univarate_test_list,plt_list)
......@@ -67,26 +66,27 @@ function multivariate_simulations()
# end
end
using ProgressMeter
function run_multivariate_sims(sims,samples)
function run_multivariate_sims(sims)
varnames, sim_ranges = zip(sims...)
without_app_future = @spawn mean_solve(samples,get_parameters(),HeatmapRecorder)
simvars = Iterators.product(sim_ranges...)
progmeter = Progress((length(simvars)+1)*samples)
without_app, _ = mean_solve(samples,get_parameters(),HeatmapRecorder; progmeter)
next!(progmeter)
default_parameters = get_app_parameters()
simvars = Iterators.product(sim_ranges...)
progmeter = Progress(length(simvars)*samples)
app_output = ThreadsX.map(simvars) do vars
vars_with_names = NamedTuple{varnames}(vars)
parameters = merge(default_parameters,vars_with_names)
out,_ = mean_solve(samples, parameters,HeatmapRecorder;progmeter)
return out
end
display(length(simvars))
fname = join(string.(varnames),"_")
keyed_output = KeyedArray(app_output;NamedTuple{varnames}(sim_ranges)...)
path = joinpath(PACKAGE_FOLDER,"abm_output","$fname.dat")
serialize(path,(fetch(without_app_future)[1],keyed_output))
serialize(path,(without_app,keyed_output))
return fname
end
using ColorSchemes
......@@ -114,7 +114,7 @@ function plot_parameter_plane(input_fname)
"Change in average final size of infection outbreak",
]
for (fname,title,datamap) in zip(fnames,titles,datamaps)
p = heatmap(var_ranges[1],var_ranges[2],datamap; title, xlabel = vars[1], ylabel = vars[2], seriescolor=cs, size = (800,600))
p = heatmap(var_ranges[1],var_ranges[2],datamap; title, xlabel = vars[1], ylabel = vars[2], seriescolor=cs, size = (600,400))
savefig(p,joinpath(PACKAGE_FOLDER,"plots","app_heatmaps","$fname"))
end
end
\ No newline at end of file
......@@ -19,6 +19,8 @@ function plot_model(varname,univariate_series, output_list::Vector{T},infection_
(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"),
(data.daily_total_notifications, "total notifications each day"),
(data.daily_total_notified_agents, "total notified agents each day"),
]
l = length(ts_list(output_list[1]))
......
......@@ -14,6 +14,9 @@ end
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
modelsol.daily_total_notifications = 0
for (i,node) in enumerate(app_user_index)
for j in 2:14
covid_alert_notifications[j-1,i] = covid_alert_notifications[j,i] #shift them all back
......@@ -34,6 +37,7 @@ Base.@propagate_inbounds @views function update_alert_durations!(t,modelsol) # B
covid_alert_notifications[end,i] = 0
end
if sum(covid_alert_notifications[:,i])>=notification_threshold
modelsol.daily_total_notifications += 1
time_of_last_alert[i] = t
end
end
......@@ -105,6 +109,8 @@ end
Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,modelsol,total_infections)
@unpack infection_introduction_day, π_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
modelsol.daily_total_notified_agents = 0
for i in 1:nodes
π_base = t<infection_introduction_day ?
(π_base_y,π_base_m,π_base_o) :
......@@ -112,14 +118,14 @@ Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,mod
random_soc_network = sample(Random.default_rng(Threads.threadid()), 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))
app_vac_payoff = 0.0
if app_user[i] && time_of_last_alert[app_user_list[i]]>=0
modelsol.daily_total_notified_agents += 1
app_vac_payoff = Γ^((t - time_of_last_alert[app_user_list[i]])) * (η + total_infections*ω_en)
# display(t - time_of_last_alert[app_user_list[i]])
end
if u_vac[random_neighbour] == u_vac[i]
app_vac_payoff = 0.0
if app_user[i] && time_of_last_alert[app_user_list[i]]>=0
app_vac_payoff = Γ^(-1*(t - time_of_last_alert[app_user_list[i]])) * (η + total_infections*ω_en)
end
vac_payoff = π_base[Int(demographics[i])] + total_infections*ω + app_vac_payoff
if u_vac[i]
# display(1 - Φ(vac_payoff,ξ))
if rand(Random.default_rng(Threads.threadid())) < 1 - Φ(vac_payoff,ξ)
......@@ -134,7 +140,6 @@ Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,mod
end
end
end
modelsol.daily_vaccinators = count(==(true),u_vac)
end
......@@ -188,10 +193,11 @@ function solve!(modelsol,recording::T) where T
update_vaccination_opinion_state!(t,modelsol,modelsol.status_totals[Int(Infected)])
#advance agent states based on the new network
modelsol.u_vac .= modelsol.u_next_vac
modelsol.u_inf .= modelsol.u_next_inf
record!(t,modelsol,recording)
modelsol.u_vac .= modelsol.u_next_vac
modelsol.u_inf .= modelsol.u_next_inf
end
end
......
No preview for this file type
Markdown is supported
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