From 92a2cad8712eedefa03fd21e567273eca55e3f72 Mon Sep 17 00:00:00 2001
From: pjentsch <pjentsch@uwaterloo.ca>
Date: Fri, 16 Apr 2021 00:25:28 -0400
Subject: [PATCH] app dynamics and opinion dynamics done, I think? people seem
to rapidly switch opinions which is probably wrong
---
CovidAlertVaccinationModel/src/ABM/abm.jl | 15 +-
CovidAlertVaccinationModel/src/ABM/agents.jl | 34 ----
.../src/ABM/contact_vectors.jl | 37 ----
.../src/ABM/{graphs.jl => mixing_graphs.jl} | 59 ++++--
CovidAlertVaccinationModel/src/ABM/model.jl | 170 ------------------
.../src/ABM/model_setup.jl | 117 ++++++++++++
CovidAlertVaccinationModel/src/ABM/output.jl | 5 +
CovidAlertVaccinationModel/src/ABM/solve.jl | 130 ++++++++++++++
.../src/CovidAlertVaccinationModel.jl | 7 +-
9 files changed, 311 insertions(+), 263 deletions(-)
delete mode 100644 CovidAlertVaccinationModel/src/ABM/contact_vectors.jl
rename CovidAlertVaccinationModel/src/ABM/{graphs.jl => mixing_graphs.jl} (54%)
delete mode 100644 CovidAlertVaccinationModel/src/ABM/model.jl
create mode 100644 CovidAlertVaccinationModel/src/ABM/model_setup.jl
create mode 100644 CovidAlertVaccinationModel/src/ABM/output.jl
create mode 100644 CovidAlertVaccinationModel/src/ABM/solve.jl
diff --git a/CovidAlertVaccinationModel/src/ABM/abm.jl b/CovidAlertVaccinationModel/src/ABM/abm.jl
index 9174035..250c1d0 100644
--- a/CovidAlertVaccinationModel/src/ABM/abm.jl
+++ b/CovidAlertVaccinationModel/src/ABM/abm.jl
@@ -7,14 +7,15 @@ default(dpi = 300)
default(framestyle = :box)
-function abm()
- agent_model = AgentModel(5000)
+function bench()
+
steps = 100
- # return agent_model
- sol,graphs = solve!(get_parameters(),steps,agent_model);
- b1 = @benchmark solve!($get_parameters(),$steps,$agent_model) seconds = 20
+ model_sol = ModelSolution(steps,get_parameters(),5000);
+ solve!(model_sol,identity)
+end
+
+function abm()
+ b1 = @benchmark bench() seconds = 20
display(b1)
println("done")
- # sol,graphs = solve_record!(u_0,get_parameters(),steps,agent_model);
- # return aggregate_timeseries(sol)
end
diff --git a/CovidAlertVaccinationModel/src/ABM/agents.jl b/CovidAlertVaccinationModel/src/ABM/agents.jl
index 28774ea..555add2 100644
--- a/CovidAlertVaccinationModel/src/ABM/agents.jl
+++ b/CovidAlertVaccinationModel/src/ABM/agents.jl
@@ -11,40 +11,6 @@ end
Immunized
end
-struct AgentModel{GraphType,L1,L2}
- demographics::Vector{AgentDemographic}
- demographic_index_vectors::Vector{Vector{Int}}
- app_user::Vector{Bool}
- app_user_index::Vector{Int}
- base_network::GraphType
- ws_matrix_list::L1
- rest_matrix_list::L2
- function AgentModel(num_households)
- pop_list,base_network,index_vectors = generate_population(num_households)
- pop_sizes = length.(index_vectors)
- map_symmetrize(m_tuple) = map(md -> symmetrize_means(pop_sizes,md), m_tuple)
-
- workschool_contacts_mean_adjusted = map_symmetrize(workschool_mixing)
- rest_contacts_mean_adjusted = map_symmetrize(rest_mixing)
-
- app_user = rand(RNG,length(pop_list)) .< 0.1 #replace with parameter
- app_user_index = findall(==(true),app_user)
-
- return new{
- typeof(base_network),
- typeof(workschool_contacts_mean_adjusted),
- typeof(rest_contacts_mean_adjusted)
- }(
- pop_list,
- index_vectors,
- app_user,
- app_user_index,
- base_network,
- workschool_contacts_mean_adjusted,
- rest_contacts_mean_adjusted
- )
- end
-end
#generate population with households distributed according to household_size_distribution
function generate_population(num_households)
diff --git a/CovidAlertVaccinationModel/src/ABM/contact_vectors.jl b/CovidAlertVaccinationModel/src/ABM/contact_vectors.jl
deleted file mode 100644
index d23ff4d..0000000
--- a/CovidAlertVaccinationModel/src/ABM/contact_vectors.jl
+++ /dev/null
@@ -1,37 +0,0 @@
-@inline function reindex!(k,csum,index_list_i,index_list_j,j_to_i_contacts,i_to_j_contacts,sample_list_i,sample_list_j)
- i_index = index_list_i[k]
- j_index = index_list_j[k]
- csum += j_to_i_contacts[j_index] - i_to_j_contacts[i_index] + sample_list_i[k] - sample_list_j[k]
- i_to_j_contacts[i_index] = sample_list_i[k]
- j_to_i_contacts[j_index] = sample_list_j[k]
- return csum
-end
-using StaticArrays
-function generate_contact_vectors!(ij_dist,ji_dist,i_to_j_contacts::Vector{T}, j_to_i_contacts::Vector{T}) where T
- rand!(RNG,ij_dist,i_to_j_contacts)
- rand!(RNG,ji_dist,j_to_i_contacts)
- l_i = length(i_to_j_contacts)
- l_j = length(j_to_i_contacts)
-
- csum = sum(i_to_j_contacts) - sum(j_to_i_contacts)
-
- inner_iter = 10
- index_list_i = @MVector zeros(Int,inner_iter)
- index_list_j = similar(index_list_i)
- sample_list_i = @MVector zeros(T,inner_iter)
- sample_list_j = similar(sample_list_i)
-
- while csum != 0
- sample!(RNG,1:l_i,index_list_i)
- sample!(RNG,1:l_j,index_list_j)
- rand!(RNG,ij_dist,sample_list_i)
- rand!(RNG,ji_dist,sample_list_j)
- @inbounds @simd 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)
- end
- end
- end
-
- return nothing
-end
\ No newline at end of file
diff --git a/CovidAlertVaccinationModel/src/ABM/graphs.jl b/CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl
similarity index 54%
rename from CovidAlertVaccinationModel/src/ABM/graphs.jl
rename to CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl
index 3eff7ef..1c26f96 100644
--- a/CovidAlertVaccinationModel/src/ABM/graphs.jl
+++ b/CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl
@@ -1,12 +1,10 @@
-using LinearAlgebra
-using LightGraphs
#A type that defines a matrix of vectors, such that the sum of the ijth vector is equal to the sum of the jith vector
struct MixingContacts{V}
total_edges::Int
contact_array::V
- function MixingContacts(index_vectors,mixing_matrix)
+ function MixingContacts(demographic_index_vectors,mixing_matrix)
contacts = map(CartesianIndices(mixing_matrix)) do ind
- zeros(Int,length(index_vectors[ind[1]]))
+ zeros(Int,length(demographic_index_vectors[ind[1]]))
end
tot = 0
for i in 1:length(mixing_matrix[:,1]), j in 1:i #diagonal
@@ -38,15 +36,15 @@ struct MixingGraph{G, V, M}
g::G
weights::V
covid_alert_time::M
- function MixingGraph(population_list,index_vectors,mixing_matrix)
+ function MixingGraph(demographics,demographic_index_vectors,mixing_matrix)
- contacts = MixingContacts(index_vectors,mixing_matrix)
+ contacts = MixingContacts(demographic_index_vectors,mixing_matrix)
- g = Graph(length(population_list))
+ g = Graph(length(demographics))
weights = RobinDict{Tuple{Int,Int},UInt8}()
- for i in 1:length(index_vectors), j in 1:i #diagonal
- random_bipartite_graph_fast_CL!(g,index_vectors[i],index_vectors[j],contacts.contact_array[i,j],contacts.contact_array[j,i],weights)
+ for i in 1:length(demographic_index_vectors), j in 1:i #diagonal
+ random_bipartite_graph_fast_CL!(g,demographic_index_vectors[i],demographic_index_vectors[j],contacts.contact_array[i,j],contacts.contact_array[j,i],weights)
end
covid_alert_time = zeros(nv(g))
@@ -67,12 +65,12 @@ struct MixingGraph{G, V, M}
end
end
-function MixingGraph(g::SimpleGraph,population_list,contact_time_distribution_matrix)
+function MixingGraph(g::SimpleGraph,demographics,contact_time_distribution_matrix)
mg = MixingGraph(g)
return mg
end
-function MixingGraph(population_list,index_vectors,mixing_matrix,contact_time_distribution_matrix)
- mg = MixingGraph(population_list,index_vectors,mixing_matrix)
+function MixingGraph(demographics,demographic_index_vectors,mixing_matrix,contact_time_distribution_matrix)
+ mg = MixingGraph(demographics,demographic_index_vectors,mixing_matrix)
return mg
end
@@ -89,4 +87,41 @@ function sample_mixing_graph!(mixing_graph,population_demographics, contact_time
mixing_graph.weights[(i,j)] = contact_time
mixing_graph.weights[(j,i)] = contact_time
end
+end
+@inline function reindex!(k,csum,index_list_i,index_list_j,j_to_i_contacts,i_to_j_contacts,sample_list_i,sample_list_j)
+ i_index = index_list_i[k]
+ j_index = index_list_j[k]
+ csum += j_to_i_contacts[j_index] - i_to_j_contacts[i_index] + sample_list_i[k] - sample_list_j[k]
+ i_to_j_contacts[i_index] = sample_list_i[k]
+ j_to_i_contacts[j_index] = sample_list_j[k]
+ return csum
+end
+using StaticArrays
+function generate_contact_vectors!(ij_dist,ji_dist,i_to_j_contacts::Vector{T}, j_to_i_contacts::Vector{T}) where T
+ rand!(RNG,ij_dist,i_to_j_contacts)
+ rand!(RNG,ji_dist,j_to_i_contacts)
+ l_i = length(i_to_j_contacts)
+ l_j = length(j_to_i_contacts)
+
+ csum = sum(i_to_j_contacts) - sum(j_to_i_contacts)
+
+ inner_iter = 10
+ index_list_i = @MVector zeros(Int,inner_iter)
+ index_list_j = similar(index_list_i)
+ sample_list_i = @MVector zeros(T,inner_iter)
+ sample_list_j = similar(sample_list_i)
+
+ while csum != 0
+ sample!(RNG,1:l_i,index_list_i)
+ sample!(RNG,1:l_j,index_list_j)
+ rand!(RNG,ij_dist,sample_list_i)
+ rand!(RNG,ji_dist,sample_list_j)
+ @inbounds @simd 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)
+ end
+ end
+ end
+
+ return nothing
end
\ No newline at end of file
diff --git a/CovidAlertVaccinationModel/src/ABM/model.jl b/CovidAlertVaccinationModel/src/ABM/model.jl
deleted file mode 100644
index 3445351..0000000
--- a/CovidAlertVaccinationModel/src/ABM/model.jl
+++ /dev/null
@@ -1,170 +0,0 @@
-function get_u_0(nodes)
- #vaccinator with 50% probability
- is_vaccinator = rand(nodes) .> 0.5
-
- u_0 = (
- status = fill(Susceptible,nodes),
- is_vaccinator = is_vaccinator
- )
-
- init_indices = rand(RNG, 1:nodes, 5)
- u_0.status[init_indices] .= Infected #start with five infected
- return u_0
-end
-
-function contact_weight(p, contact_time)
- return 1 - (1-p)^contact_time
-end
-
-function EN_payoff()
- return 0.0
-end
-
-function vac_switch_probability(u, population_list,vac, params,total_infected,soc_nbrs)
- @unpack Ï€_base, η, κ, ω, Ï, ω_en,Ï_en,γ,β = params
- if vac
- soc_nbrs_vac = [0,0,0]
- for nbr in soc_nbrs
- if u.is_vaccinator[nbr]
- soc_nbrs_vac[Int(population_list[nbr])] += 1
- end
- end
- return Φ(Ï€_base + κ * (sum(soc_nbrs_vac)/length(soc_nbrs)) + dot(Ï,soc_nbrs_vac) + 0,β)
- else
- soc_nbrs_nonvac = 0
- for nbr in soc_nbrs
- if u.is_vaccinator[nbr]
- soc_nbrs_nonvac += 1
- end
- end
- return Φ(κ * soc_nbrs_nonvac/length(soc_nbrs),β)
- end
-end
-
-function agents_step!(t,u_next,u,population_list,soc_network,network_list,params,covid_alert_times,app_user_index)
- @unpack p, recovery_rate, immunization_loss_prob = params
-
- u_next.status .= u.status
- u_next.is_vaccinator .= u.is_vaccinator
-
- for (i,node) in enumerate(app_user_index), mixing_graph in network_list
- for j in 2:14
- covid_alert_times[i,j-1] = covid_alert_times[i,j] #shift them all back
- end
- for j in neighbors(mixing_graph.g,node)
- covid_alert_times[i,end] += mixing_graph.weights[(node,j)] #add the contact times for today to the back
- end
- end
-
- for i in 1:length(population_list)
- if rand(RNG)<vac_switch_probability(u, population_list,u.is_vaccinator[i],params,count(==(Infected),u),neighbors(soc_network.g,i))
- u_next.is_vaccinator[i] = !u.is_vaccinator[i]
- end
- end
- for i in 1:length(population_list)
- agent_status = u.status[i]
- is_vaccinator = u.is_vaccinator[i]
- agent_demo = population_list[i]
-
- if agent_status == Susceptible
- if is_vaccinator
- u_next.status[i] = Immunized
- else
- for mixing_graph in network_list
- for j in neighbors(mixing_graph.g,i)
- if u.status[j] == Infected && rand(RNG) < contact_weight(p,mixing_graph.weights[(i,j)])
- u_next.status[i] = Infected
- end
- end
- end
- end
- elseif agent_status == Infected
- if rand(RNG) < recovery_rate
- u_next.status[i] = Recovered
- end
- elseif agent_status == Immunized
- if rand(RNG) < immunization_loss_prob
- u_next.status[i] = Susceptible
- end
- end
- end
-end
-
-
-
-function solve!(params,steps,agent_model; record = false)
- population_list = agent_model.demographics #list of demographic status for each agent
- popsize = length(population_list)
- index_vectors = agent_model.demographic_index_vectors #lists of indices of each agent with a given demographic
- ws_contacts = agent_model.ws_matrix_list
- rest_contacts = agent_model.rest_matrix_list
- app_user_index = agent_model.app_user_index
-
- p = params.p
- u_0 = get_u_0(length(agent_model.demographics))
-
- home_static_edges = MixingGraph(agent_model.base_network,population_list,contact_time_distribution_matrix) #network with households and LTC homes
- ws_static_edges = MixingGraph(population_list,index_vectors,ws_contacts.daily,contact_time_distribution_matrix)
- ws_weekly_edges = MixingGraph(population_list,index_vectors,ws_contacts.twice_a_week,contact_time_distribution_matrix)
- ws_daily_edges_vector = [MixingGraph(population_list,index_vectors,ws_contacts.otherwise,contact_time_distribution_matrix) for i in 1:steps]
-
- rest_static_edges = MixingGraph(population_list,index_vectors,rest_contacts.daily,contact_time_distribution_matrix)
- rest_weekly_edges = MixingGraph(population_list,index_vectors,rest_contacts.twice_a_week,contact_time_distribution_matrix)
- rest_daily_edges_vector = [MixingGraph(population_list,index_vectors,rest_contacts.otherwise,contact_time_distribution_matrix) for i in 1:steps]
-
-
-
- solution = vcat([u_0], [deepcopy(u_0) for i in 1:steps])
- # @show solution
- network_lists = [
- [home_static_edges,rest_static_edges] for i in 1:steps
- ]
-
- covid_alert_times = zeros(length(app_user_index),14) #two weeks worth of values
-
- for (t,l) in enumerate(network_lists)
- day_of_week = mod(t,7)
- if !(day_of_week == 3 || day_of_week == 4) #simulation begins on thursday I guess
- push!(l, ws_static_edges)
- end
- if rand(RNG)<5/7
- push!(l, ws_weekly_edges)
- push!(l, rest_weekly_edges)
- end
- push!(l,ws_daily_edges_vector[t])
- push!(l,rest_daily_edges_vector[t])
- end
-
- for t in 1:steps
- network_list = network_lists[t]
- for network in network_list
- sample_mixing_graph!(network,population_list, contact_time_distribution_matrix) #get new contact weights
- end
- #advance agent states based on the new network
- agents_step!(t,solution[t+1],solution[t],population_list,home_static_edges,network_list,params,covid_alert_times,app_user_index)
- end
- return solution, network_lists
-end
-function get_parameters()
- params = (
- p = 0.05,
- recovery_rate = 0.1,
- immunization_loss_prob = 0.01,
- π_base = 0.1,
- η = 0.0,
- κ = 0.0,
- ω = 0.0,
- Ï = [0.0,0.0,0.0],
- ω_en = 0.0,
- Ï_en = 0.0,
- γ = 0.0,
- β = 0.1
- )
- return params
-end
-
-function Φ(payoff,β)
- return 1 / (exp(-1*β*payoff))
-end
-
-
diff --git a/CovidAlertVaccinationModel/src/ABM/model_setup.jl b/CovidAlertVaccinationModel/src/ABM/model_setup.jl
new file mode 100644
index 0000000..6bde0ba
--- /dev/null
+++ b/CovidAlertVaccinationModel/src/ABM/model_setup.jl
@@ -0,0 +1,117 @@
+
+function get_parameters()
+ params = (
+ I_0_fraction = 0.01,
+ base_transmission_probability = 0.5,
+ recovery_rate = 0.1,
+ immunization_loss_prob = 0.01,
+ π_base = 1.0,
+ η = 0.0,
+ κ = 0.0,
+ ω = 0.0,
+ Ï = [0.0,0.0,0.0],
+ ω_en = 0.0,
+ Ï_en = [0.0,0.0,0.0],
+ γ = 0.0,
+ β = 10.0,
+ notification_parameter = 0.0,
+ vaccinator_prob = 0.5,
+ app_user_fraction = 0.5,
+ )
+ return params
+end
+
+
+function get_u_0(nodes,I_0_fraction,vaccinator_prob)
+ #vaccinator with 50% probability
+ is_vaccinator = rand(nodes) .> vaccinator_prob
+ status = fill(Susceptible,nodes)
+ init_indices = rand(RNG, 1:nodes, round(Int,nodes*I_0_fraction))
+ status[init_indices] .= Infected #start with five infected
+ return status,is_vaccinator
+end
+
+struct ModelSolution{T,G}
+ sim_length::Int
+ nodes::Int
+ params::T
+ u_next_inf::Vector{AgentStatus}
+ u_next_vac::Vector{Bool}
+ u_inf::Vector{AgentStatus}
+ u_vac::Vector{Bool}
+ covid_alert_times::Array{Int,2}
+ time_of_last_alert::Vector{Int}
+ inf_network_lists::Vector{Vector{G}}
+ soc_networks::Vector{G}
+ index_vectors::Vector{Vector{Int}}
+ demographics::Vector{AgentDemographic}
+ app_user::Vector{Bool}
+ app_user_list::Vector{Int}
+ app_user_index::Vector{Int}
+ function ModelSolution(sim_length,params::T,num_households) where T
+ demographics,base_network,index_vectors = generate_population(num_households)
+ pop_sizes = length.(index_vectors)
+ map_symmetrize(m_tuple) = map(md -> symmetrize_means(pop_sizes,md), m_tuple)
+
+ ws_matrix_list = map_symmetrize(workschool_mixing)
+ rest_matrix_list = map_symmetrize(rest_mixing)
+
+ app_user_list = zeros(length(demographics))
+ is_app_user = rand(RNG,length(demographics)) .< params.app_user_fraction
+ app_user_index = findall(==(true),is_app_user)
+ app_user_list[is_app_user] .= collect(1:length(app_user_index))
+ # display(app)
+ nodes = length(demographics)
+
+ u_0_inf,u_0_vac = get_u_0(nodes,params.I_0_fraction,params.vaccinator_prob)
+
+ home_static_edges = MixingGraph(base_network,demographics,contact_time_distribution_matrix) #network with households and LTC homes
+ ws_static_edges = MixingGraph(demographics,index_vectors,ws_matrix_list.daily,contact_time_distribution_matrix)
+ ws_weekly_edges = MixingGraph(demographics,index_vectors,ws_matrix_list.twice_a_week,contact_time_distribution_matrix)
+ ws_daily_edges_vector = [MixingGraph(demographics,index_vectors,ws_matrix_list.otherwise,contact_time_distribution_matrix) for i in 1:sim_length]
+
+ rest_static_edges = MixingGraph(demographics,index_vectors,rest_matrix_list.daily,contact_time_distribution_matrix)
+ rest_weekly_edges = MixingGraph(demographics,index_vectors,rest_matrix_list.twice_a_week,contact_time_distribution_matrix)
+ rest_daily_edges_vector = [MixingGraph(demographics,index_vectors,rest_matrix_list.otherwise,contact_time_distribution_matrix) for i in 1:sim_length]
+
+ inf_network_lists = [
+ [home_static_edges,rest_static_edges] for i in 1:sim_length
+ ]
+ soc_network_list = [home_static_edges,rest_static_edges,ws_static_edges]
+
+ for (t,l) in enumerate(inf_network_lists)
+ day_of_week = mod(t,7)
+ if !(day_of_week == 3 || day_of_week == 4) #simulation begins on thursday I guess
+ push!(l, ws_static_edges)
+ end
+ if rand(RNG)<5/7
+ push!(l, ws_weekly_edges)
+ push!(l, rest_weekly_edges)
+ end
+ push!(l,ws_daily_edges_vector[t])
+ push!(l,rest_daily_edges_vector[t])
+ end
+
+ covid_alert_times = zeros(Int,length(app_user_index),14) #two weeks worth of values
+ time_of_last_alert = fill(-1,length(app_user_index)) #two weeks worth of values
+
+ return new{T,typeof(home_static_edges)}(
+ sim_length,
+ nodes,
+ params,
+ u_0_inf,
+ u_0_vac,
+ copy(u_0_inf),
+ copy(u_0_vac),
+ covid_alert_times,
+ time_of_last_alert,
+ inf_network_lists,
+ soc_network_list,
+ index_vectors,
+ demographics,
+ is_app_user,
+ app_user_list,
+ app_user_index
+ )
+ end
+end
diff --git a/CovidAlertVaccinationModel/src/ABM/output.jl b/CovidAlertVaccinationModel/src/ABM/output.jl
new file mode 100644
index 0000000..dbd0bbf
--- /dev/null
+++ b/CovidAlertVaccinationModel/src/ABM/output.jl
@@ -0,0 +1,5 @@
+abstract type AbstractOutputData end
+
+struct DebugOutputData <: AbstractOutputData
+ TotalInfected::
+end
\ No newline at end of file
diff --git a/CovidAlertVaccinationModel/src/ABM/solve.jl b/CovidAlertVaccinationModel/src/ABM/solve.jl
new file mode 100644
index 0000000..46c8764
--- /dev/null
+++ b/CovidAlertVaccinationModel/src/ABM/solve.jl
@@ -0,0 +1,130 @@
+
+function contact_weight(p, contact_time)
+ return 1 - (1-p)^contact_time
+end
+
+
+function EN_payoff()
+ return 0.0
+end
+
+function update_alert_durations!(t,modelsol)
+
+ @unpack notification_parameter = modelsol.params
+ @unpack time_of_last_alert, app_user_index,inf_network_lists,covid_alert_times,app_user = modelsol
+
+
+ for (i,node) in enumerate(modelsol.app_user_index), mixing_graph in modelsol.inf_network_lists[t]
+ for j in 2:14
+ covid_alert_times[i,j-1] = covid_alert_times[i,j] #shift them all back
+ end
+ for j in neighbors(mixing_graph.g,node)
+ if app_user[j]
+ covid_alert_times[i,end] += mixing_graph.weights[(node,j)] #add the contact times for today to the back
+ end
+ end
+ if rand(RNG) < 1 - (1- notification_parameter)^sum(covid_alert_times[i,:])
+ time_of_last_alert[i] = t
+ end
+ end
+end
+function update_infection_state!(t,modelsol)
+ @unpack base_transmission_probability,immunization_loss_prob,recovery_rate = modelsol.params
+ @unpack u_inf,u_vac,u_next_inf,u_next_vac,demographics,inf_network_lists = modelsol
+
+ for i in 1:modelsol.nodes
+ agent_status = u_inf[i]
+ is_vaccinator = u_vac[i]
+ agent_demo = demographics[i]
+ if agent_status == Susceptible
+ if is_vaccinator
+ u_next_inf[i] = Immunized
+ else
+ for mixing_graph in inf_network_lists[t]
+ for j in neighbors(mixing_graph.g,i)
+ if u_inf[j] == Infected && rand(RNG) < contact_weight(base_transmission_probability,mixing_graph.weights[(i,j)])
+ u_next_inf[i] = Infected
+ end
+ end
+ end
+ end
+ elseif agent_status == Infected
+ if rand(RNG) < recovery_rate
+ u_next_inf[i] = Recovered
+ end
+ elseif agent_status == Immunized
+ if rand(RNG) < immunization_loss_prob
+ u_next_inf[i] = Susceptible
+ end
+ end
+ end
+end
+function update_vaccination_opinion_state!(t,modelsol,total_infections)
+
+ @unpack Ï€_base, η,γ, κ, ω, Ï, ω_en,Ï_en,γ,β = modelsol.params
+ @unpack demographics,time_of_last_alert, nodes, soc_networks,u_vac,u_next_vac,app_user,app_user_list = modelsol
+ app_user_pointer = 0
+ for i in 1:nodes
+ vac_payoff = 0
+ soc_nbrs_vac = [0,0,0]
+ soc_nbrs_nonvac = 0
+ num_soc_nbrs = 0
+ for soc_network in soc_networks
+ soc_nbrs = neighbors(soc_network.g,i)
+ num_soc_nbrs += length(soc_nbrs)
+ for nbr in soc_nbrs
+ if u_vac[nbr]
+ soc_nbrs_vac[Int(demographics[nbr])] += 1
+ else
+ soc_nbrs_nonvac += 1
+ end
+ end
+ end
+
+ vac_payoff += Ï€_base + κ * (sum(soc_nbrs_vac)/num_soc_nbrs) + dot(Ï,soc_nbrs_vac) + total_infections*ω- κ * soc_nbrs_nonvac/num_soc_nbrs
+ if app_user[i] && time_of_last_alert[app_user_list[i]]>=0
+ vac_payoff += γ^(-1*(t - time_of_last_alert[app_user_list[i]]))* (η + dot(Ï_en,soc_nbrs_vac) + total_infections*ω_en)
+ end
+ if rand(RNG)<Φ(vac_payoff,β)
+ u_next_vac[i] = !u_vac[i]
+ end
+ if u_vac[i]
+ if rand(RNG)< 1- Φ(vac_payoff,β)
+ u_next_vac[i] = !u_vac[i]
+ end
+ else
+ if rand(RNG)< Φ(vac_payoff,β)
+ u_next_vac[i] = !u_vac[i]
+ end
+ end
+ end
+end
+
+
+function agents_step!(t,modelsol)
+ for network in modelsol.inf_network_lists[t]
+ sample_mixing_graph!(network,modelsol.demographics, contact_time_distribution_matrix) #get new contact weights
+ end
+ update_alert_durations!(t,modelsol)
+ update_vaccination_opinion_state!(t,modelsol,count(==(Infected),modelsol.u_inf))
+ update_infection_state!(t,modelsol)
+ modelsol.u_vac .= modelsol.u_next_vac
+ modelsol.u_inf .= modelsol.u_next_inf
+end
+
+
+
+function solve!(modelsol,recording_function)
+ for t in 1:modelsol.sim_length
+ #advance agent states based on the new network
+ agents_step!(t,modelsol)
+ display((count(==(true),modelsol.u_vac),count(==(Immunized),modelsol.u_inf)))
+ end
+ # return solution, network_lists
+end
+
+
+function Φ(payoff,β)
+ return 1 / (exp(-1*β*payoff))
+end
+
diff --git a/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl b/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl
index de68843..20f401c 100644
--- a/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl
+++ b/CovidAlertVaccinationModel/src/CovidAlertVaccinationModel.jl
@@ -7,6 +7,7 @@ using Plots
using Distributions
using StatsBase
using Dates
+using LinearAlgebra
using ThreadsX
using DelimitedFiles
using KernelDensity
@@ -33,10 +34,10 @@ include("data.jl")
include("ABM/abm.jl")
include("ABM/agents.jl")
include("ABM/mixing_distributions.jl")
-include("ABM/contact_vectors.jl")
-include("ABM/model.jl")
+include("ABM/mixing_graphs.jl")
include("ABM/plotting.jl")
-include("ABM/graphs.jl")
+include("ABM/model_setup.jl")
+include("ABM/solve.jl")
include("IntervalsModel/intervals_model.jl")
include("IntervalsModel/interval_overlap_sampling.jl")
include("IntervalsModel/hh_durations_model.jl")
--
GitLab