Skip to content
Snippets Groups Projects
Commit 7316cefb authored by Peter Jentsch's avatar Peter Jentsch
Browse files

I think it works

parent 144ab904
No related branches found
No related tags found
No related merge requests found
......@@ -29,7 +29,7 @@ function solve_and_plot_parameters()
display(sum.(eachrow(ymo_vaccination_ts)) ./avg_populations)
l = length(p)
# savefig(plot(p...; layout = (l,1), size=(800,400*l),leftmargin = 12Plots.mm),"timeseries.pdf")
savefig(plot(p...; layout = (l,1), size=(800,400*l),leftmargin = 12Plots.mm),"timeseries.pdf")
return out
end
......
......@@ -19,7 +19,7 @@
π_base_o = -0.95,
η = 0.5,
κ = 0.0,
ω = 0.0053,
ω = 0.053,
ω_en = 0.01,
Γ = 1/7,
ξ = 5.0,
......@@ -74,7 +74,6 @@ mutable struct ModelSolution{T,InfNet,SocNet,WSMixingDist,RestMixingDist,Recorde
index_vectors::Vector{Vector{Int}}
demographics::Vector{AgentDemographic}
is_app_user::Vector{Bool}
status_totals::Vector{Int}
ws_matrix_tuple::WSMixingDist
rest_matrix_tuple::RestMixingDist
immunization_countdown::Vector{Int}
......@@ -101,9 +100,9 @@ mutable struct ModelSolution{T,InfNet,SocNet,WSMixingDist,RestMixingDist,Recorde
covid_alert_notifications = zeros(Bool,14,nodes) #two weeks worth of values
time_of_last_alert = fill(-1,nodes) #time of last alert is negative if no alert has been recieved
status_totals = [count(==(AgentStatus(i)), u_0_inf) for i in 1:AgentStatus.size]
immunization_countdown = fill(-1, nodes) #immunization countdown is negative if not counting down
output_data = Recorder(0,sim_length; record_degrees)
output_data.recorded_status_totals[:,1] = [count(==(AgentStatus(i)), u_0_inf) for i in 1:AgentStatus.size]
return new{T,typeof(infected_mixing_graph),
typeof(soc_mixing_graph),typeof(ws_matrix_tuple),
......@@ -122,7 +121,6 @@ mutable struct ModelSolution{T,InfNet,SocNet,WSMixingDist,RestMixingDist,Recorde
index_vectors,
demographics,
is_app_user,
status_totals,
ws_matrix_tuple,
rest_matrix_tuple,
immunization_countdown,
......
......@@ -73,7 +73,6 @@ Initialize a Recorder filled with (copies) of val.
function record!(t,modelsol)
recorder = modelsol.output_data
recorder.recorded_status_totals[:,t] .= modelsol.status_totals
alerts = filter(>(0),modelsol.time_of_last_alert)
if !isempty(alerts)
mean_alert_time = mean(t .- alerts)
......
......@@ -33,14 +33,16 @@ end
time_of_last_alert[agent] = t
end
end
function agent_transition!(modelsol,agent, from::AgentStatus,to::AgentStatus)
@unpack u_next_inf,status_totals, immunization_countdown = modelsol
function agent_transition!(t,modelsol,agent, from::AgentStatus,to::AgentStatus)
@unpack u_next_inf,output_data, immunization_countdown = modelsol
immunization_countdown[agent] = -1
status_totals[Int(from)] -= 1
status_totals[Int(to)] += 1
if t<size(output_data.recorded_status_totals,2)-1
output_data.recorded_status_totals[Int(from), t+1] -= 1
output_data.recorded_status_totals[Int(to), t+1] += 1
end
u_next_inf[agent] = to
end
function infect_susceptible_agent!(t, modelsol, agent,agent_inf_status,agent_demo,α_vec,β_vec)
function is_susceptible_infected!(t, modelsol, agent,agent_inf_status,agent_demo,α_vec,β_vec)
@unpack u_inf,inf_network, output_data = modelsol
for mixing_graph in inf_network.graph_list[t], neighbor in neighbors(mixing_graph,agent)
......@@ -50,14 +52,14 @@ function infect_susceptible_agent!(t, modelsol, agent,agent_inf_status,agent_dem
if rand(Random.default_rng(Threads.threadid())) < infection_threshold && agent_inf_status == Susceptible
output_data.daily_cases_by_age[Int(agent_demo),t]+=1
output_data.daily_unvac_cases_by_age[Int(agent_demo),t]+=1
agent_transition!(modelsol, agent, Susceptible,Infected)
agent_transition!(t,modelsol, agent, Susceptible,Infected)
return true
end
end
end
return false
end
function infect_immunized_agent!(t, modelsol, agent,agent_inf_status,agent_demo,α_vec,β_vec)
function is_immunized_infected!(t, modelsol, agent,agent_inf_status,agent_demo,α_vec,β_vec)
@unpack u_inf,inf_network, output_data = modelsol
for mixing_graph in inf_network.graph_list[t], neighbor in neighbors(mixing_graph,agent)
if u_inf[neighbor] == Infected
......@@ -65,7 +67,7 @@ function infect_immunized_agent!(t, modelsol, agent,agent_inf_status,agent_demo,
if rand(Random.default_rng(Threads.threadid())) < infection_threshold && agent_inf_status == Immunized
immunization_ineffective = rand(Random.default_rng(Threads.threadid())) < 1- α_vec[Int(agent_demo)]
if immunization_ineffective
agent_transition!(modelsol,agent, Immunized,Infected)
agent_transition!(t,modelsol,agent, Immunized,Infected)
output_data.daily_cases_by_age[Int(agent_demo),t]+=1
return true
end
......@@ -80,10 +82,9 @@ function vaccinate_agent!(t,agent,immunization_countdown)
end
end
function update_vaccination_opinion_state!(t,agent,modelsol,total_infections,π_base_vec)
@unpack immunization_intervals, η,Γ,ζ, ω, ω_en,ξ = modelsol.params
@unpack demographics,time_of_last_alert, soc_network,u_vac,u_next_vac,is_app_user,output_data = modelsol
function update_vaccination_opinion_state!(t,agent,modelsol,vaccinate_today,total_infections,π_base_vec)
@unpack infection_introduction_day,immunization_intervals, η,Γ,ζ, ω, ω_en,ξ = modelsol.params
@unpack immunization_countdown,demographics,time_of_last_alert, soc_network,u_vac,u_next_vac,is_app_user,output_data = modelsol
random_soc_network = sample(Random.default_rng(Threads.threadid()), soc_network.graph_list[t])
if !isempty(neighbors(random_soc_network,agent))
random_neighbour = sample(Random.default_rng(Threads.threadid()), neighbors(random_soc_network.g,agent))
......@@ -101,6 +102,9 @@ function update_vaccination_opinion_state!(t,agent,modelsol,total_infections,π_
else
if rand(Random.default_rng(Threads.threadid())) < Φ(vac_payoff,ξ)
u_next_vac[agent] = true
if vaccinate_today && t>infection_introduction_day
vaccinate_agent!(t, agent, immunization_countdown)
end
end
end
return vac_payoff
......@@ -136,7 +140,7 @@ function sample_initial_nodes(nodes,graphs,I_0_fraction)
end
function agents_step!(t,modelsol,vaccinate_today)
@unpack recovery_rate,β_y,β_m,β_o,α_y,α_m,α_o,π_base_y,π_base_m,π_base_o,ζ,immunization_intervals = modelsol.params
@unpack infection_introduction_day,recovery_rate,β_y,β_m,β_o,α_y,α_m,α_o,π_base_y,π_base_m,π_base_o,ζ,immunization_intervals = modelsol.params
@unpack u_inf,u_vac,u_next_vac,u_next_inf,demographics,inf_network, is_app_user, immunization_countdown, output_data = modelsol
u_next_inf .= u_inf
......@@ -153,16 +157,16 @@ function agents_step!(t,modelsol,vaccinate_today)
for (agent,(agent_inf_status,agent_is_vaccinator,agent_demo)) in enumerate(zip(u_inf,u_vac,demographics))
if agent_inf_status == Susceptible
infect_susceptible_agent!(t, modelsol, agent,agent_inf_status,agent_demo,α_vec,β_vec)
if vaccinate_today && agent_is_vaccinator
is_susceptible_infected!(t, modelsol, agent,agent_inf_status,agent_demo,α_vec,β_vec)
if vaccinate_today && agent_is_vaccinator && t<(infection_introduction_day-1)
vaccinate_agent!(t, agent, immunization_countdown)
end
elseif agent_inf_status == Immunized
infect_immunized_agent!(t, modelsol, agent,agent_inf_status,agent_demo,α_vec,β_vec)
is_immunized_infected!(t, modelsol, agent,agent_inf_status,agent_demo,α_vec,β_vec)
elseif agent_inf_status == Infected
if rand(Random.default_rng(Threads.threadid())) < recovery_rate
agent_transition!(modelsol, agent, Infected,Recovered)
agent_transition!(t,modelsol, agent, Infected,Recovered)
end
end
if t>first(minimum(modelsol.params.immunization_intervals))
......@@ -174,16 +178,15 @@ function agents_step!(t,modelsol,vaccinate_today)
if immunization_countdown[agent] == 0
output_data.daily_immunized_by_age[Int(agent_demo),t] += 1
fit!(output_data.avg_weighted_degree_of_vaccinators[Int(agent_demo)],weighted_degree_of_i)
agent_transition!(modelsol, agent, Susceptible,Immunized)
agent_transition!(t,modelsol, agent, Susceptible,Immunized)
elseif immunization_countdown[agent]>0
fit!(output_data.avg_weighted_degree[Int(agent_demo)],weighted_degree_of_i)
immunization_countdown[agent] -= 1
else
fit!(output_data.avg_weighted_degree[Int(agent_demo)],weighted_degree_of_i)
end
update_vaccination_opinion_state!(t, agent, modelsol, modelsol.status_totals[Int(Infected)],π_base_vec) #very bad pls change
total_infected = modelsol.output_data.recorded_status_totals[Int(Infected), t]
update_vaccination_opinion_state!(t, agent, modelsol,vaccinate_today, total_infected,π_base_vec)
end
output_data.total_vaccinators[t] = count(==(true),u_vac)
......@@ -209,7 +212,7 @@ function solve!(modelsol)
if !isempty(init_indices)
inf_index = pop!(init_indices)
modelsol.u_inf[inf_index] = Infected
modelsol.status_totals[Int(Infected)] += 1
modelsol.output_data.recorded_status_totals[1,Int(Infected)] += 1
end
end
......@@ -224,7 +227,7 @@ function solve!(modelsol)
vaccinate_today = false
end
total_weight += agents_step!(t,modelsol,vaccinate_today)
agents_step!(t,modelsol,vaccinate_today)
#advance agent states based on the new network
record!(t,modelsol)
......
No preview for this file type
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment