Skip to content
GitLab
Commit 2ce4277d authored by Peter Jentsch's avatar Peter Jentsch
Browse files

basic tests

parent 688f2ea5
Expand all Hide whitespace changes
Inline Side-by-side
Manifest.toml
View file @ 2ce4277d
......@@ -15,10 +15,10 @@ version = "2.1.0"
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
[[ArnoldiMethod]]
deps = ["DelimitedFiles", "LinearAlgebra", "Random", "SparseArrays", "Test"]
git-tree-sha1 = "a5f9581eaf9371fe12b4a0b12c05690b45912e71"
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "f87e559f87a45bece9c9ed97458d3afe98b1ebb9"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.0.2"
version = "0.1.0"
[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
......@@ -58,9 +58,9 @@ version = "0.9.0"
[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "89a0b14325d0f02f9caed7c8ba91181a5d254874"
git-tree-sha1 = "53fed426c9af1eb68e63b3999e96454c2db79757"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "0.9.26"
version = "0.9.27"
[[ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random", "StaticArrays"]
......@@ -112,15 +112,15 @@ uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"
[[DataAPI]]
git-tree-sha1 = "ad84f52c0b8f05aa20839484dbaf01690b41ff84"
git-tree-sha1 = "25ccd31003243d2ce83e474cf11663dddf48035f"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.4.0"
version = "1.4.1"
[[DataFrames]]
deps = ["CategoricalArrays", "Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "b46e1deb4592a5df7416b10dfcd6b01fb194ab9a"
git-tree-sha1 = "c5f9f2385a52e35c75efd90911091f8a749177a5"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "0.22.2"
version = "0.22.4"
[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
......@@ -152,9 +152,9 @@ uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
[[Distributions]]
deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"]
git-tree-sha1 = "5c6a645ea65af6aac4f96b5dedaae660b7b56267"
git-tree-sha1 = "f0e06a5b5ccda38e2fb8f59d91316e657b67047d"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.24.10"
version = "0.24.12"
[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
......@@ -186,9 +186,9 @@ version = "4.3.1+4"
[[FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays"]
git-tree-sha1 = "ff537e5a3cba92fb48f30fec46723510450f2c0e"
git-tree-sha1 = "50eabdace27aa27b143f65b65e762bb0112a7708"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.10.2"
version = "0.11.1"
[[FixedPointNumbers]]
deps = ["Statistics"]
......@@ -244,9 +244,9 @@ version = "0.53.0+0"
[[GeometryBasics]]
deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
git-tree-sha1 = "876f77f0d3253e882ff588af1c95d0e4a86c9766"
git-tree-sha1 = "f574945bcabe9805b78292216279c1be910168bb"
uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
version = "0.3.5"
version = "0.3.8"
[[Gettext_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
......@@ -271,6 +271,12 @@ git-tree-sha1 = "c7ec02c4c6a039a98a15f955462cd7aea5df4508"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.8.19"
[[ImportAll]]
deps = ["Test"]
git-tree-sha1 = "5aa65f6204e1b68dfef024421e61fbcacecfbfa6"
uuid = "c65182e5-40f4-518f-8165-175b85689199"
version = "1.1.0"
[[Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
......@@ -422,9 +428,9 @@ version = "2.34.0+7"
[[LightGraphs]]
deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
git-tree-sha1 = "677464beb4a91772d151ee62ea5e181d4b7f50ec"
git-tree-sha1 = "432428df5f360964040ed60418dd5601ecd240b6"
uuid = "093fc24a-ae57-5d10-9952-331d41423f4d"
version = "1.3.4"
version = "1.3.5"
[[LinearAlgebra]]
deps = ["Libdl"]
......@@ -466,9 +472,9 @@ version = "0.1.0"
[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "ed61674a0864832495ffe0a7e889c0da76b0f4c8"
git-tree-sha1 = "f8c673ccc215eb50fcadb285f522420e29e69e1c"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "0.4.4"
version = "0.4.5"
[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"
......@@ -571,9 +577,9 @@ version = "1.0.10"
[[Plots]]
deps = ["Base64", "Contour", "Dates", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "JSON", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs"]
git-tree-sha1 = "4797acb266b8d9ff316f4581924e71c6709f152d"
git-tree-sha1 = "3acf7ee21b0c0ea99ef0815e7768c1c0fde82629"
uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
version = "1.10.1"
version = "1.10.2"
[[PooledArrays]]
deps = ["DataAPI"]
......@@ -583,9 +589,9 @@ version = "0.5.3"
[[PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "237170206bf38a66fee4d845f4ae57f63788eeb0"
git-tree-sha1 = "42126c4e2677cdc664baea004c98cc60a664fe40"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "0.10.1"
version = "0.11.0"
[[Printf]]
deps = ["Unicode"]
......@@ -629,10 +635,9 @@ uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c"
version = "0.2.1"
[[Reexport]]
deps = ["Pkg"]
git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0"
git-tree-sha1 = "57d8440b0c7d98fc4f889e478e80f268d534c9d5"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "0.2.0"
version = "1.0.0"
[[Referenceables]]
git-tree-sha1 = "4a32b1dd124a846580608eb347d4337f873c2499"
......@@ -712,9 +717,9 @@ version = "1.2.1"
[[SplittablesBase]]
deps = ["Setfield", "Test"]
git-tree-sha1 = "e6e5e266aa58c2e419d19e8eccbe6af55d5dc6e1"
git-tree-sha1 = "c2b4f3cd29827787a5f61afe12ea13222197ab98"
uuid = "171d559e-b47b-412a-8079-5efa626c420e"
version = "0.1.11"
version = "0.1.12"
[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
......@@ -746,9 +751,9 @@ version = "0.4.4"
[[StructTypes]]
deps = ["Dates", "UUIDs"]
git-tree-sha1 = "d94235fcdc4a09649f263365c5f7e4ed4ba6ed34"
git-tree-sha1 = "65a43f5218197bc7091b76bc273a5e323a1d7b0d"
uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
version = "1.2.1"
version = "1.2.3"
[[SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
......@@ -766,9 +771,9 @@ version = "1.0.0"
[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "240d19b8762006ff04b967bdd833269ad642d550"
git-tree-sha1 = "8dc2bb7d3548e315d890706547b24502ed79504f"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.2.2"
version = "1.3.1"
[[Tar]]
deps = ["ArgTools", "SHA"]
......@@ -786,9 +791,9 @@ version = "0.1.7"
[[Transducers]]
deps = ["ArgCheck", "BangBang", "CompositionsBase", "DefineSingletons", "Distributed", "InitialValues", "Logging", "Markdown", "MicroCollections", "Requires", "Setfield", "SplittablesBase", "Tables"]
git-tree-sha1 = "eb5d4ff48c25762319f1aae0e5b501394beaefd7"
git-tree-sha1 = "02413f5795ad272f9f9912e4e4c83d9b1572750c"
uuid = "28d57a85-8fef-5791-bfe6-a80928e7c999"
version = "0.4.56"
version = "0.4.58"
[[UUIDs]]
deps = ["Random", "SHA"]
......@@ -958,9 +963,9 @@ uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
[[Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "6f1abcb0c44f184690912aa4b0ba861dd64f11b9"
git-tree-sha1 = "2c1332c54931e83f8f94d310fa447fd743e8d600"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.4.5+2"
version = "1.4.8+0"
[[ZygoteRules]]
deps = ["MacroTools"]
......
Project.toml
View file @ 2ce4277d
......@@ -9,6 +9,7 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
ImportAll = "c65182e5-40f4-518f-8165-175b85689199"
LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
NamedTupleTools = "d9ec5142-1e00-5aa0-9d6a-321866360f50"
NetworkLayout = "46757867-2c16-5918-afeb-47bfcb05e46a"
......
covidabm.code-workspace 0 → 100644
View file @ 2ce4277d
{
"folders": [
{
"name": "CovidAlertVaccinationModel",
"path": "."
}
],
"settings": {}
}
\ No newline at end of file
data/csv/COVID_ontario_data.csv
View file @ 2ce4277d
This diff is collapsed.
data/csv/distancing_data.csv 0 → 100644
View file @ 2ce4277d
This diff is collapsed.
data/csv/home_compositions.csv 0 → 100644
View file @ 2ce4277d
This diff is collapsed.
fetch_data.sh
View file @ 2ce4277d
......@@ -2,3 +2,6 @@
echo "fetching full cases data..."
curl https://data.ontario.ca/dataset/f4112442-bdc8-45d2-be3c-12efae72fb27/resource/455fd63b-603d-4608-8216-7d8647f43350/download/conposcovidloc.csv > data/csv/COVID_ontario_data.csv
echo "done."
echo "fetching google distancing data..."
curl https://www.gstatic.com/covid19/mobility/Global_Mobility_Report.csv?cachebust=fbb340e43a0602e1 | { head -1; grep ",CA-ON,"; } > data/csv/distancing_data.csv
echo "done."
\ No newline at end of file
src/CovidAlertVaccinationModel.jl
View file @ 2ce4277d
......@@ -8,6 +8,7 @@ using Plots
using DataFrames
using Distributions
using StatsBase
import StatsBase:mean
using Dates
using ThreadsX
using DelimitedFiles
......@@ -28,12 +29,14 @@ include("graphs.jl")
const population = 14.57e6 #population of ontario
const color_palette = palette(:seaborn_bright) #color theme for the plots
const age_bins = [(0.0, 20.0),(20.0,60.0),(60.0,Inf)]
const age_bins = [(0.0, 25.0),(25.0,65.0),(65.0,Inf)]
const PACKAGE_FOLDER = dirname(dirname(pathof(CovidAlertVaccinationModel)))
const infection_data = parse_cases_data()
const demographic_distribution = get_canada_demographic_distribution()
#shift distribution,shift means, geometric changes to poisson if mean goes below 1
#ZW expectation is (1 - α)(expecation of base distribution)
#expectation/(1 - α)
# CUDA.allowscalar(false)
......@@ -42,26 +45,16 @@ default(framestyle = :box)
using BenchmarkTools
function main()
vaccination_rate_test(5000,0.1,vaccinate_uniformly!)
rng = Xoroshiro128Plus()
agent_model = AgentModel(rng,50,2)
u_0 = get_u_0(rng,length(agent_model.demographics))
steps = 10
sol1,graphs = solve!(rng,u_0,get_parameters(),steps,agent_model,vaccinate_uniformly!);
plot_model(agent_model.base_network,graphs,sol1)
# @btime solve!($u_0,$base_network,$params,$steps,$r,$pop_list,$index_vectors,vaccinate_uniformly!);
end
#returns the total infected given a certain num_households, vaccines per day, and strategy
function vaccination_rate_test(num_households,vac_per_day,vac_strategy; rng = Xoroshiro128Plus())
agent_model = AgentModel(rng,num_households,2)
agent_model = AgentModel(rng,500,2)
u_0 = get_u_0(rng,length(agent_model.demographics))
params = merge(get_parameters(),(vaccines_per_day = vac_per_day,))
steps = 300
sol1,graphs = solve!(rng,u_0,params,steps,agent_model,vaccinate_uniformly!);
total_infections = count(x->x == Recovered,sol1[end])
# sol1,graphs = solve!(rng,u_0,get_parameters(),steps,agent_model,vaccinate_uniformly!);
# plot_model(agent_model.base_network,graphs,sol1)
# @btime solve!($rng,$u_0,$get_parameters(),$steps,$agent_model,$vaccinate_uniformly!);
end
#returns the total infected given a certain num_households, vaccines per day, and strategy
end
\ No newline at end of file
src/agents.jl
View file @ 2ce4277d
......@@ -15,21 +15,21 @@ end
end
struct AgentModel{GraphType}
struct AgentModel{GraphType,DistMatrix1,DistMatrix2}
demographics::Vector{AgentDemographic}
demographic_index_vectors::Vector{Vector{Int}}
base_network::GraphType
workschool_contacts_mean_adjusted::DistMatrix1
rest_contacts_mean_adjusted::DistMatrix2
function AgentModel(rng,num_households,num_ltc)
pop_list,base_network,index_vectors = generate_population(rng,num_households,num_ltc)
return new{typeof(base_network)}(pop_list,index_vectors,base_network)
workschool_contacts_mean_adjusted = symmetrize_means(length.(index_vectors),initial_workschool_mixing_matrix)
rest_contacts_mean_adjusted = symmetrize_means(length.(index_vectors),initial_rest_mixing_matrix)
return new{typeof(base_network),typeof(workschool_contacts_mean_adjusted),typeof(rest_contacts_mean_adjusted)}(pop_list,index_vectors,base_network,workschool_contacts_mean_adjusted,rest_contacts_mean_adjusted)
end
end
# func(a) = println("agent")
function Base.show(io::IO, status::AgentStatus)
if status == Susceptible
print(io, "S")
......
src/data.jl
View file @ 2ce4277d
......@@ -8,7 +8,7 @@ function get_canada_demographic_distribution()::Vector{Float64}
df = DataFrame([f[1,i] => f[2:end, i] for i = 1:length(f[1,:])])
binned_data = df[:,:demographic_data]
source_bins = map(parse_string_as_float_pair, df[:,:demographic_data_bins])
return [sum(binned_data[1:4]),sum(binned_data[5:12]),sum(binned_data[13:end])]
return [sum(binned_data[1:5]),sum(binned_data[6:13]),sum(binned_data[14:end])]
end
......@@ -18,6 +18,31 @@ function get_canada_case_fatality()::Tuple{Vector{Tuple{Float64,Float64}},Vector
return map(parse_string_as_float_pair, df[:,:case_fatality_bins]), df[:,:case_fatality_data]
# https://www.publichealthontario.ca/-/media/documents/ncov/epi/covid-19-severe-outcomes-ontario-epi-summary.pdf?la=en
end
function find_household_composition(df_row)
age_resp_to_bin = Dict(
"M" => 2,
"O" => 3,
"Y" => 1
)
u25_bins = [:U15CHILD,:O15CHILD,:YSPOUSE]
m_bins = [:MPAR, :MCHILD,:MHHADULT]
o_bins = [:OPAR, :OSPOUSE,:OHHADULT]
age_distribution = [0,0,0]
age_distribution[1] += sum(df_row[field ] for field in u25_bins)
age_distribution[2] += sum(df_row[field] for field in m_bins)
age_distribution[3] += sum(df_row[field] for field in o_bins)
age_distribution[age_resp_to_bin[df_row[:AGERESP]]] += 1
return age_distribution
end
function sample_household_data(rng,n)
f = readdlm(joinpath(PACKAGE_FOLDER,"data/csv/home_compositions.csv"), ',')
df = DataFrame([f[1,i] => f[2:end, i] for i = 1:length(f[1,:])])
weight_vector::Vector{Float64} = df[!,:WGHT_PER]/sum(df[!,:WGHT_PER])
households = map(find_household_composition,eachrow(df))
return sample(rng,households, Weights(weight_vector),n)
# https://www.publichealthontario.ca/-/media/documents/ncov/epi/covid-19-severe-outcomes-ontario-epi-summary.pdf?la=en
end
function parse_string_as_float_pair(s)::Tuple{Float64,Float64}
parsed = strip(s, ['(',')']) |> s -> split(s, ",")
......
src/graphs.jl
View file @ 2ce4277d
#add a bipartite graph derived from mixing matrices onto g
#is there a better way to do this?
function generate_mixing_graph!(rng,g,index_vectors,mixing_matrix)
for i in 1:5, j in 1:i #diagonal
agent_contact_dist = mixing_matrix[i,j]
i_to_j_contacts = rand(rng,agent_contact_dist,length(index_vectors[i]))
agent_contact_dist_i = mixing_matrix[i,j]
agent_contact_dist_j = mixing_matrix[j,i]
i_to_j_contacts = rand(rng,agent_contact_dist_i,length(index_vectors[i]))
j_to_i_contacts = rand(rng,agent_contact_dist,length(index_vectors[j]))
equalize_degree_lists!(i_to_j_contacts,j_to_i_contacts)
j_to_i_contacts = rand(rng,agent_contact_dist_j,length(index_vectors[j]))
while sum(i_to_j_contacts) - sum(j_to_i_contacts) != 0
i_to_j_contacts[sample(rng,1:length(index_vectors[i]))] = rand(rng,agent_contact_dist_i)
j_to_i_contacts[sample(rng,1:length(index_vectors[j]))] = rand(rng,agent_contact_dist_j)
end
random_bipartite_graph_fast_CL!(rng,g,index_vectors[i],index_vectors[j],i_to_j_contacts,j_to_i_contacts)
end
return g
end
#Remove contacts from the degree list with the larger sum uniformly at random until they are equal
#needed to satisfy constraint that a bipartite graph exists with degree lists l1 and l2
function equalize_degree_lists!(l1,l2)
if sum(l1) > sum(l2)
map(i -> l1[i]-=1, sample(1:length(l1),sum(l1) - sum(l2)))
else
map(i -> l2[i]-=1, sample(1:length(l2),sum(l2) - sum(l1)))
end
end
#modify g so that nodes specified in anodes and bnodes are connected by a bipartite graph with expected degrees given by aseq and bseq
#implemented from Aksoy, S. G., Kolda, T. G., & Pinar, A. (2017). Measuring and modeling bipartite graphs with community structure
#simple algorithm, might produce parallel edges for small graphs
......@@ -40,62 +32,64 @@ function random_bipartite_graph_fast_CL!(rng,g,anodes,bnodes,aseq,bseq)
return g
end
function symmetrize_means(population_sizes,mixing_matrix::Matrix{<:ZWDist})
mixing_means = mean.(mixing_matrix)
symmetrized_mixing = similar(mixing_matrix)
for i in 1:length(population_sizes), j in 1:length(population_sizes)
symmetrized_mean = 0.5*(mixing_means[i,j] + population_sizes[j]/population_sizes[i] * mixing_means[j,i])
symmetrized_mixing[i,j] = from_mean(typeof(mixing_matrix[i,j]),mixing_matrix[i,j].α,symmetrized_mean)
end
return symmetrized_mixing
end
function symmetrize_means(population_sizes,mixing_matrix::Matrix{<:Distribution})
mixing_means = mean.(mixing_matrix)
symmetrized_mixing = similar(mixing_matrix)
for i in 1:length(population_sizes), j in 1:length(population_sizes)
symmetrized_mean = 0.5*(mixing_means[i,j] + population_sizes[j]/population_sizes[i] * mixing_means[j,i])
symmetrized_mixing[i,j] = from_mean(typeof(mixing_matrix[i,j]),symmetrized_mean)
end
return symmetrized_mixing
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(rng,num_households,num_LTC)
household_sizes = Int.(rand(rng,household_size_distribution,num_households))
total_household_population = sum(household_sizes)
# LTC_sizes = rand(LTC_distribution,num_LTC)
prob_adult_is_HCW = 0.01 #assume 1% of adults are HCW
prob_adult_is_HCW = 0.01#assume 1% of adults are HCW
prob_elderly_in_LTC = 0.05 #5% of elderly are in LTC
agent_distribution_in_households = Categorical([ #find a better way to do this maybe
demographic_distribution[1],
demographic_distribution[2] * (1 - prob_adult_is_HCW),
demographic_distribution[2] * prob_adult_is_HCW,
demographic_distribution[3] * (1 - prob_elderly_in_LTC),
demographic_distribution[3] * prob_elderly_in_LTC
])
households = Vector{Vector{AgentDemographic}}()
ltc_pop_list = Vector{AgentDemographic}()
households_composition = sample_household_data(rng,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)
for household_size in household_sizes
household = Vector{AgentDemographic}()
while length(household)<household_size
agent_sample = AgentDemographic(rand(rng,agent_distribution_in_households))
if agent_sample == ElderlyLTC
push!(ltc_pop_list,agent_sample)
else
push!(household,agent_sample)
end
end
push!(households,household)
end
ltc_pop = length(ltc_pop_list)
println("ltc: $ltc_pop")
ltc_home_size = Int(ceil(ltc_pop/num_LTC))
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, household_sizes)
# display(household_networks)
population_list = vcat(vcat(households...), vcat(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()))
# display(network)
index_vectors = [findall(x -> x == AgentDemographic(i), population_list) for i in 1:AgentDemographic.size+1]
return (;
population_list,
network,
......
src/mixing_distributions.jl
View file @ 2ce4277d
......@@ -10,6 +10,24 @@ function Base.rand(rng::AbstractRNG, s::ZWDist)
return Base.rand(rng,s.base_dist)
end
end
# function from_mean(::Type{Geometric{T}},μ) where T
# return Geometric(1/μ)
# end
function from_mean(::Type{Geometric{T}},μ) where T
if μ > 1.0
return Geometric(1/(μ+1))
else
return Poisson(μ)
end
end
function from_mean(::Type{Poisson{T}},μ) where T
return Poisson(μ)
end
function from_mean(::Type{ZWDist{DistType,T}},α,μ) where {DistType <: Distribution{Univariate,W} where W, T}
return ZWDist(α,from_mean(DistType,μ/(1-α)))
end
function ZeroPoisson(α,λ)
return ZWDist(α,Poisson(λ))
end
......@@ -17,37 +35,25 @@ function ZeroGeometric(α,p)
return ZWDist(α,Geometric(p))
end
StatsBase.mean(d::ZWDist{Dist,T}) where {Dist,T} = (1 - d.α)*StatsBase.mean(d.base_dist)
#constant for now anyway..
const home_mixing_matrix = [
Poisson(1.749) Poisson(2.629) Poisson(0.648) Poisson(2.629) Poisson(0.648)
Poisson(0.726) Poisson(1.909) Poisson(0.357) Poisson(1.909) Poisson(0.357)
Poisson(0.022) Poisson(0.120) Poisson(0.677) Poisson(0.120) Poisson(0.677)
Poisson(0.726) Poisson(1.909) Poisson(0.357) Poisson(1.909) Poisson(0.357)
Poisson(0.022) Poisson(0.120) Poisson(0.677) Poisson(0.120) Poisson(0.677)
]
const workschool_mixing_matrix = [
ZeroGeometric(0.439,0.126) ZeroGeometric(0.414,0.230) ZeroPoisson(0.907,0.243) ZeroGeometric(0.414,0.230) ZeroPoisson(0.907,0.243)
ZeroGeometric(0.767,0.457) ZeroGeometric(0.786,0.106) ZeroPoisson(0.842,0.122) ZeroGeometric(0.786,0.106) ZeroPoisson(0.842,0.122)
ZeroPoisson(0.898,0.015) ZeroPoisson(0.786,0.036) ZeroPoisson(0.941,0.373) ZeroPoisson(0.786,0.036) ZeroPoisson(0.941,0.373)
ZeroGeometric(0.767,0.457) ZeroGeometric(0.786,0.106) ZeroPoisson(0.842,0.122) ZeroGeometric(0.786,0.106) ZeroPoisson(0.842,0.122)
ZeroPoisson(0.898,0.015) ZeroPoisson(0.786,0.036) ZeroPoisson(0.941,0.373) ZeroPoisson(0.786,0.036) ZeroPoisson(0.941,0.373)
]
const initial_workschool_mixing_matrix = map(t->from_mean(t...),[
(ZWDist{Geometric{Float64},Float64}, 0.433835,4.104848) (ZWDist{Geometric{Float64},Float64},0.406326,2.568782) (ZWDist{Poisson{Float64},Float64},0.888015,0.017729) (ZWDist{Geometric{Float64},Float64},0.406326,2.568782) (ZWDist{Poisson{Float64},Float64},0.888015,0.017729)
(ZWDist{Geometric{Float64},Float64}, 0.600966,0.975688) (ZWDist{Geometric{Float64},Float64},0.4306,5.057572) (ZWDist{Poisson{Float64},Float64}, 0.84513,0.021307) (ZWDist{Geometric{Float64},Float64},0.4306,5.057572) (ZWDist{Poisson{Float64},Float64}, 0.84513,0.021307)
(ZWDist{Poisson{Float64},Float64},0.887392,0.001937) (ZWDist{Poisson{Float64},Float64},0.793004,0.00722) (ZWDist{Poisson{Float64},Float64},0.940473, 0.022134) (ZWDist{Poisson{Float64},Float64},0.793004,0.00722) (ZWDist{Poisson{Float64},Float64},0.940473, 0.022134)
(ZWDist{Geometric{Float64},Float64}, 0.600966,0.975688) (ZWDist{Geometric{Float64},Float64},0.4306,5.057572) (ZWDist{Poisson{Float64},Float64}, 0.84513,0.021307) (ZWDist{Geometric{Float64},Float64},0.4306,5.057572) (ZWDist{Poisson{Float64},Float64}, 0.84513,0.021307)
(ZWDist{Poisson{Float64},Float64},0.887392,0.001937) (ZWDist{Poisson{Float64},Float64},0.793004,0.00722) (ZWDist{Poisson{Float64},Float64},0.940473, 0.022134) (ZWDist{Poisson{Float64},Float64},0.793004,0.00722) (ZWDist{Poisson{Float64},Float64},0.940473, 0.022134)
])
const rest_mixing_matrix = [
Geometric(0.450) Geometric(0.784) Poisson(0.155) Geometric(0.784) Poisson(0.155)
Poisson(0.850) Geometric(0.269) Poisson(0.741) Geometric(0.269) Poisson(0.741)
Poisson(0.217) Poisson(0.783) Poisson(0.667) Poisson(0.783) Poisson(0.667)
Poisson(0.850) Geometric(0.269) Poisson(0.741) Geometric(0.269) Poisson(0.741)
Poisson(0.217) Poisson(0.783) Poisson(0.667) Poisson(0.783) Poisson(0.667)
]
const initial_rest_mixing_matrix = 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)
])
const contact_time_distribution_matrix = [Hypergeometric(5,2,4) for i in 1:AgentDemographic.size + 1, j in 1:AgentDemographic.size + 1]
const household_size_distribution = Truncated(Poisson(2.60),0,Inf) #https://www.researchgate.net/publication/226081704_Household_size_and_the_Poisson_distribution
src/model.jl
View file @ 2ce4277d
function get_u_0(rng,nodes)
u_0 = fill(Susceptible,nodes)
init_indices = rand(rng, 1:nodes, 20)
init_indices = rand(rng, 1:nodes, 5)
u_0[init_indices] .= Infected #start with two infected
return u_0