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

removed LTC/HCW stuff, contact generation needs fixing, use bytearrays instead...

removed LTC/HCW stuff, contact generation needs fixing, use bytearrays instead of bitarrays for weights
parent ba3e1ee0
No related branches found
No related tags found
No related merge requests found
......@@ -41,11 +41,15 @@ default(dpi = 300)
default(framestyle = :box)
using BenchmarkTools
function main()
agent_model = AgentModel(5000,5)
# display(size.(agent_model.demographic_index_vectors))
agent_model = AgentModel(5000)
u_0 = get_u_0(length(agent_model.demographics))
steps = 300
# Random.seed!(RNG,1)
Random.seed!(RNG,1)
@btime solve!($u_0,$get_parameters(),$steps,$agent_model,$vaccinate_uniformly!);
# agent_model = AgentModel(10_000)
@btime solve!($u_0,$get_parameters(),$steps,$agent_model,$vaccinate_uniformly!);
# sol,graphs = solve!(u_0,get_parameters(),steps,agent_model,vaccinate_uniformly!);
# return aggregate_timeseries(sol)
......
#Agent type definitions
#agent type definitions
@enum AgentDemographic begin
Young = 1
Adult
Elderly
AdultHCW
ElderlyLTC
end
@enum AgentStatus begin
Susceptible = 1
......@@ -22,8 +19,8 @@ struct AgentModel{GraphType,DistMatrixList1,DistMatrixList2,DistMatrixList3}
workschool_contacts_mean_adjusted_list::DistMatrixList1
rest_contacts_mean_adjusted_list::DistMatrixList2
home_contacts_duration_list::DistMatrixList3
function AgentModel(num_households,num_ltc)
pop_list,base_network,index_vectors = generate_population(num_households,num_ltc)
function AgentModel(num_households)
pop_list,base_network,index_vectors = generate_population(num_households)
pop_sizes = length.(index_vectors)
workschool_contacts_mean_adjusted_list = [symmetrize_means(pop_sizes,initial_workschool_mixing_matrix)]
......@@ -65,6 +62,8 @@ function Base.show(io::IO, status::AgentStatus)
print(io, "I")
elseif status == Recovered
print(io, "R")
elseif status == Immunized
print(io, "Vac")
end
end
......@@ -73,12 +72,7 @@ function Base.show(io::IO, status::AgentDemographic)
print(io, "<20")
elseif status == Adult
print(io, "20-60")
elseif status == AdultHCW
print(io, "20-60 HCW")
elseif status == Elderly
print(io, ">60")
elseif status == ElderlyLTC
print(io, ">60 LTC")
end
end
......@@ -86,7 +86,7 @@ function generate_contact_vectors!(ij_dist,ji_dist,i_to_j_contacts, j_to_i_conta
end
#add a bipartite graph derived from mixing matrices onto g
function generate_mixing_graph!(g,index_vectors,contacts)
for i in 1:5, j in 1:5
for i in 1:length(index_vectors), j in 1:length(index_vectors)
random_bipartite_graph_fast_CL!(g,index_vectors[i],index_vectors[j],contacts.contact_array[i,j],contacts.contact_array[j,i])
end
return g
......@@ -109,7 +109,7 @@ function random_bipartite_graph_fast_CL!(g::SimpleGraph,anodes,bnodes,aseq,bseq)
return g
end
#same algorithm but with a bitarray
#same algorithm but with a bitarray, this cannot produce parallel edges
function random_bipartite_graph_fast_CL!(g::T,anodes,bnodes,aseq,bseq) where T<:AbstractArray
lena = length(aseq)
lenb = length(bseq)
......@@ -124,41 +124,18 @@ function random_bipartite_graph_fast_CL!(g::T,anodes,bnodes,aseq,bseq) where T<:
end
#generate population with households distributed according to household_size_distribution
#Assume 5% of elderly are in LTC roughly, these are divide evenly among a specific number of LTC ltc_homes
#LTC homes and households are assumed to be complete static graphs.
function generate_population(num_households,num_LTC)
# LTC_sizes = rand(LTC_distribution,num_LTC)
prob_adult_is_HCW = 0.01#assume 1% of adults are HCW
prob_elderly_in_LTC = 0.05 #5% of elderly are in LTC
function generate_population(num_households)
households_composition = sample_household_data(num_households)
households = map( l -> [fill(AgentDemographic(i),n) for (i,n) in enumerate(l)],households_composition) |>
l -> reduce(vcat,l) |>
l -> reduce(vcat,l)
total_household_pop = sum(sum.(households_composition))
ltc_pop = Int(div(total_household_pop,(1 - prob_elderly_in_LTC))) - total_household_pop
hcw_pop = Int(ceil(prob_adult_is_HCW*total_household_pop))
ltc_pop_list = fill(ElderlyLTC,ltc_pop)
adults_indices = findall(x -> x == Adult, households)
households[sample(RNG,adults_indices,hcw_pop)] .= AdultHCW
ltc_home_size = Int(ceil(ltc_pop/num_LTC))
@assert ltc_home_size != 0 error("must have at least one LTC resident,make problem larger")
if ltc_home_size != 0
ltc_homes = [ltc_pop_list[i:min(i + ltc_home_size - 1,ltc_pop)] for i = 1:ltc_home_size:ltc_pop]
else
ltc_homes = Vector{Vector{AgentDemographic}}()
end
ltc_networks = map(LightGraphs.complete_graph,length.(ltc_homes))
household_networks = map(LightGraphs.complete_graph, sum.(households_composition))
population_list = vcat(reduce(vcat,households), reduce(vcat,ltc_homes))
network = blockdiag(reduce(blockdiag,household_networks; init = SimpleGraph()), reduce(blockdiag,ltc_networks; init = SimpleGraph()))
index_vectors = [findall(x -> x == AgentDemographic(i), population_list) for i in 1:AgentDemographic.size+1]
population_list = reduce(vcat,households)
network = reduce(blockdiag,household_networks; init = SimpleGraph())
index_vectors = [findall(x -> x == AgentDemographic(i), population_list) for i in 1:(AgentDemographic.size-1)]
return (;
population_list,
network,
......
......@@ -56,26 +56,24 @@ function from_mean(::Type{ZWDist{DistType,T}},α,μ) where {DistType <: Distribu
end
StatsBase.mean(d::ZWDist{Dist,T}) where {Dist,T} = (1 - d.α)*StatsBase.mean(d.base_dist)
const initial_workschool_type = Union{ZWDist{Geometric{Float64},Discrete},ZWDist{Poisson{Float64},Discrete}}
const initial_workschool_mixing_matrix = convert(Array{initial_workschool_type},map(t->from_mean(t...),[
(ZWDist{Geometric{Float64},Discrete}, 0.433835,4.104848) (ZWDist{Geometric{Float64},Discrete},0.406326,2.568782) (ZWDist{Poisson{Float64},Discrete},0.888015,0.017729) (ZWDist{Geometric{Float64},Discrete},0.406326,2.568782) (ZWDist{Poisson{Float64},Discrete},0.888015,0.017729)
(ZWDist{Geometric{Float64},Discrete}, 0.600966,0.975688) (ZWDist{Geometric{Float64},Discrete},0.4306,5.057572) (ZWDist{Poisson{Float64},Discrete}, 0.84513,0.021307) (ZWDist{Geometric{Float64},Discrete},0.4306,5.057572) (ZWDist{Poisson{Float64},Discrete}, 0.84513,0.021307)
(ZWDist{Poisson{Float64},Discrete},0.887392,0.001937) (ZWDist{Poisson{Float64},Discrete},0.793004,0.00722) (ZWDist{Poisson{Float64},Discrete},0.940473, 0.022134) (ZWDist{Poisson{Float64},Discrete},0.793004,0.00722) (ZWDist{Poisson{Float64},Discrete},0.940473, 0.022134)
(ZWDist{Geometric{Float64},Discrete}, 0.600966,0.975688) (ZWDist{Geometric{Float64},Discrete},0.4306,5.057572) (ZWDist{Poisson{Float64},Discrete}, 0.84513,0.021307) (ZWDist{Geometric{Float64},Discrete},0.4306,5.057572) (ZWDist{Poisson{Float64},Discrete}, 0.84513,0.021307)
(ZWDist{Poisson{Float64},Discrete},0.887392,0.001937) (ZWDist{Poisson{Float64},Discrete},0.793004,0.00722) (ZWDist{Poisson{Float64},Discrete},0.940473, 0.022134) (ZWDist{Poisson{Float64},Discrete},0.793004,0.00722) (ZWDist{Poisson{Float64},Discrete},0.940473, 0.022134)
(ZWDist{Geometric{Float64},Discrete}, 0.433835,4.104848) (ZWDist{Geometric{Float64},Discrete},0.406326,2.568782) (ZWDist{Poisson{Float64},Discrete},0.888015,0.017729)
(ZWDist{Geometric{Float64},Discrete}, 0.600966,0.975688) (ZWDist{Geometric{Float64},Discrete},0.4306,5.057572) (ZWDist{Poisson{Float64},Discrete}, 0.84513,0.021307)
(ZWDist{Poisson{Float64},Discrete},0.887392,0.001937) (ZWDist{Poisson{Float64},Discrete},0.793004,0.00722) (ZWDist{Poisson{Float64},Discrete},0.940473, 0.022134)
]))
const initial_rest_type = Union{Geometric{Float64},Poisson{Float64}}
const initial_rest_mixing_matrix = convert(Array{initial_rest_type},map(t->from_mean(t...),[
(Geometric{Float64},2.728177) (Geometric{Float64},1.382557) (Poisson{Float64},0.206362) (Geometric{Float64},1.382557) (Poisson{Float64},0.206362)
(Poisson{Float64},1.139072) (Geometric{Float64},3.245594) (Poisson{Float64},0.785297) (Geometric{Float64},3.245594) (Poisson{Float64},0.785297)
(Poisson{Float64},0.264822) (Poisson{Float64},0.734856) (Poisson{Float64},0.667099) (Poisson{Float64},0.734856) (Poisson{Float64},0.667099)
(Poisson{Float64},1.139072) (Geometric{Float64},3.245594) (Poisson{Float64},0.785297) (Geometric{Float64},3.245594) (Poisson{Float64},0.785297)
(Poisson{Float64},0.264822) (Poisson{Float64},0.734856) (Poisson{Float64},0.667099) (Poisson{Float64},0.734856) (Poisson{Float64},0.667099)
(Geometric{Float64},2.728177) (Geometric{Float64},1.382557) (Poisson{Float64},0.206362)
(Poisson{Float64},1.139072) (Geometric{Float64},3.245594) (Poisson{Float64},0.785297)
(Poisson{Float64},0.264822) (Poisson{Float64},0.734856) (Poisson{Float64},0.667099)
]))
const contact_time_distribution_matrix = convert(Array{initial_rest_type},[Geometric() for i in 1:AgentDemographic.size + 1, j in 1:AgentDemographic.size + 1])
const contact_time_distribution_matrix = convert(Array{initial_rest_type},[Geometric() for i in 1:(AgentDemographic.size-1), j in 1:(AgentDemographic.size-1)])
......@@ -3,12 +3,12 @@ using RandomNumbers.Xorshifts
using Test
using ThreadsX
import StatsBase.mean
model_sizes = [(100,1),(1000,1),(1000,3),(5000,3)]
model_sizes = [100,1000,1000,5000]
vaccination_strategies = [vaccinate_uniformly!]
vaccination_rates = [0.000,0.005,0.01,0.05]
infection_rates = [0.01,0.05,0.1]
agent_models = ThreadsX.map(model_size -> AgentModel(model_size...), model_sizes)
dem_cat = AgentDemographic.size + 1
dem_cat = AgentDemographic.size -1
......
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