some basic stuff done, working on LTC home allocations

Project.toml

@@ -3,6 +3,19 @@ uuid = "9260c4ec-b5cf-4bc2-ad29-e1d23bf2bd6f"
 authors = ["pjentsch <pjentsch@uwaterloo.ca> and contributors"]
 version = "0.1.0"
 
+[deps]
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
+NetworkLayout = "46757867-2c16-5918-afeb-47bfcb05e46a"
+Pipe = "b98c9c47-44ae-5843-9183-064241ee97a0"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+RandomNumbers = "e6cf234a-135c-5ec9-84dd-332b85af5143"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
+
 [compat]
 julia = "1"
 

data/csv/case_fatality_data.csv

+case_fatality_bins,case_fatality_data
+"(0.0, 20.0)",0.0
+"(20.0, 29.0)",0.0011
+"(30.0, 39.0)",0.0027
+"(40.0, 49.0)",0.0067
+"(50.0, 59.0)",0.0203
+"(60.0, 69.0)",0.0652
+"(70.0, 79.0)",0.2089
+"(80.0, 125.0)",0.3700

data/csv/demographic_data.csv

+demographic_data_bins,demographic_data
+"(0, 5)",0.05169494947785886
+"(5, 10)",0.054253579120547776
+"(10, 15)",0.054051659753256126
+"(15, 20)",0.05625635853132738
+"(20, 25)",0.06588844441798299
+"(25, 30)",0.06984638325700568
+"(30, 35)",0.0692734536794045
+"(35, 40)",0.06863718154402712
+"(40, 45)",0.06440693089425378
+"(45, 50)",0.06375240886612778
+"(50, 55)",0.06657930661155305
+"(55, 60)",0.07314924139771618
+"(60, 65)",0.06682461602997154
+"(65, 70)",0.05577675347816086
+"(70, 75)",0.0454055203318437
+"(75, 80)",0.03097365944561508
+"(80, 85)",0.02092895572145045
+"(85, 90)",0.013589732088914116
+"(90, 95)",0.00645274706377582
+"(95, 100)",0.0019709352101672014
+"(100, 125)",0.00028718307903996624

fetch_data.sh

+#!/bin/bash
+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."

src/CovidAlertVaccinationModel.jl

 module CovidAlertVaccinationModel
+export get_u_0, solve, agents_step!,get_graph,plot_network,main,parse_cases_data,generate_population
+using LightGraphs
+using RandomNumbers.Xorshifts
+using UnPack
+using Plots
+using DataFrames
+using Distributions
+using StatsBase
+using Dates
+using DelimitedFiles
+# using CUDA
+using NetworkLayout:SFDP
+
+
+
+include("data.jl")
+include("model.jl")
+include("plotting.jl")
+include("types.jl")
+
+
+const population = 14.57e6
+
+const color_palette = palette(:seaborn_bright)
+const age_bins = [(0.0, 20.0),(20.0,60.0),(60.0,Inf)]
+const PACKAGE_FOLDER = dirname(dirname(pathof(CovidAlertVaccinationModel)))
+const infection_data = parse_cases_data()
+const demographic_distribution = get_canada_demographic_distribution()
+
+
+
+
+# CUDA.allowscalar(false)
+default(dpi = 300)
+default(framestyle = :box)
+
+using BenchmarkTools
+function main()
+    graph = get_graph()
+    u_0 = get_u_0(nv(graph))
+    params = get_parameters()
+    steps = 200
+    
+    r = Xoroshiro128Plus()
+    # @btime sol1 = solve!($u_0,$graph,$params,$steps,$r)
+    # display(sol1)
+    household_size_dist = Truncated(Poisson(2.60),0,Inf) #https://www.researchgate.net/publication/226081704_Household_size_and_the_Poisson_distribution
+    pop = generate_population(2,1,household_size_dist)
+    return pop
+end
+
+
+
 
-# Write your package code here.
 
 end

src/data.jl

+#Data parsing and import functions
+#These functions specify output types 
+#so the compiler is able to do type inference throughout the rest of the code
+#I don't think it is necessary in all cases, but it's good when importing data from files to specify because type instability makes things rly slow
+
+function get_canada_demographic_distribution()::Vector{Float64}
+    f = readdlm(joinpath(PACKAGE_FOLDER,"data/csv/demographic_data.csv"), ',')
+    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])]
+end
+
+
+function get_canada_case_fatality()::Tuple{Vector{Tuple{Float64,Float64}},Vector{Float64}}
+    f = readdlm(joinpath(PACKAGE_FOLDER,"data/csv/case_fatality_data.csv"), ',')
+    df = DataFrame([f[1,i] => f[2:end, i] for i = 1:length(f[1,:])])
+    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 parse_string_as_float_pair(s)::Tuple{Float64,Float64}
+    parsed = strip(s, ['(',')']) |> s -> split(s, ",")
+    if length(parsed) > 2
+        error("input tuple too long")
+    end
+    return (parse(Float64, parsed[begin]), parse(Float64, parsed[end]))
+end
+
+bin_str_to_int(bin_str) = occursin('<',bin_str) ? 1 : parse(Int, strip(bin_str, 's')) #this function parses the bin formatting in the input csv
+
+bin_to_index(bin) = findfirst(((lower,upper),) -> lower <= bin < upper, age_bins)
+
+
+
+function parse_cases_data()::Tuple{Vector{Date},Matrix{Float64}} #this type declaration probably not needed
+    f = readdlm(joinpath(PACKAGE_FOLDER,"data/csv/COVID_ontario_data.csv"), ',')
+    df = DataFrame([f[1,i] => f[2:end, i] for i = 1:length(f[1,:])])
+
+    dropmissing!(df,:Accurate_Episode_Date)
+    map!(x -> Date(x),df[!,:Accurate_Episode_Date],df[!,:Accurate_Episode_Date])
+    sort!(df,:Accurate_Episode_Date) 
+
+    #there are some cases before jan 1 2020 but they are neglible so we assume that cases start here
+    begin_date = Date(2020,1,1) 
+    cases_by_date = groupby(df,:Accurate_Episode_Date)
+
+    #Ontario health says that case data with 2 weeks of current date are subject to change
+
+    end_date = Date(maximum(df[!,:Accurate_Episode_Date])) - Day(14) 
+    dates = collect(begin_date:Day(1):end_date)
+    infection_data = zeros(Float64,length(dates), length(age_bins))
+
+
+    for (i,day) in enumerate(dates)
+        cases_on_date = get(cases_by_date, (Accurate_Episode_Date = day,),nothing)
+
+        if !isnothing(cases_on_date)
+            cases_by_age_group = groupby(cases_on_date,:Age_Group)
+
+            for age_group_key in keys(cases_by_age_group)
+                cases = cases_by_age_group[age_group_key]
+                bin_str = age_group_key[:Age_Group]
+                num_cases = nrow(cases)
+                if length(bin_str)>0 && !occursin("UNKNOWN",bin_str)
+                    bin = bin_str_to_int(bin_str)
+                    bin_index = bin_to_index(bin)
+                    infection_data[i,bin_index] = num_cases
+                else
+                    #account for missing age data?
+                end
+            end
+        end
+    end
+
+    return dates, infection_data
+end
\ No newline at end of file

src/model.jl

+
+function generate_population(num_households,num_LTC,household_size_distribution)
+    household_sizes = rand(household_size_distribution,num_households)
+    # LTC_sizes = rand(LTC_distribution,num_LTC)
+
+    prob_adult_is_HCW = 0.01
+    prob_elderly_in_LTC = 0.05
+    agent_distribution_in_households = [ #find a better way to do this
+        demographic_distribution[1],
+        demographic_distribution[2] * (1 - prob_adult_is_HCW),
+        demographic_distribution[2] * prob_adult_is_HCW,
+        demographic_distribution[3],
+        demographic_distribution[3] * 0.0
+    ]
+    display(agent_distribution_in_households)
+    agents_demographics = map(household_sizes) do household_size
+        return AgentDemographic.(rand(Categorical(agent_distribution_in_households),Int(household_size)))
+    end
+    display(agents_demographics)
+end
+
+function get_graph(population)
+    num_households = 100
+    household_sizes = rand(household_size_distribution,num_households)
+    households = map(nodes -> LightGraphs.complete_graph(nodes),household_sizes)
+    
+    num_households = 100
+    long_term_care_sizes = rand(household_size_distribution,num_households)
+    households = map(nodes -> LightGraphs.complete_graph(nodes),long_term_care_sizes)
+
+
+    return LightGraphs.grid([100,100]) 
+end
+function get_graph()
+    return LightGraphs.grid([100,100]) 
+end
+
+function get_u_0(nodes)
+    u_0 = fill(Susceptible,nodes)
+    u_0[1] = Infected
+
+    return u_0
+end
+
+
+function agents_step!(u_next,u,graph,params,t,rand_gen)
+    @unpack p, recovery_rate  = params  
+    u_next .= u
+    @inbounds for (v,agent) in enumerate(u)
+        if agent == Susceptible
+            for w in LightGraphs.neighbors(graph,v)
+                if u[w] == Infected && rand(rand_gen) < p
+                    u_next[v] = Infected
+                end
+            end
+        elseif agent == Infected
+            if rand(rand_gen) < recovery_rate
+                u_next[v] = Recovered
+            end
+        end
+    end
+end
+function solve!(u_0,graph,params,steps,rand_gen)
+    solution = vcat([u_0], [fill(Susceptible,length(u_0)) for i in 1:steps])
+    for t in 1:steps
+        agents_step!(solution[t+1],solution[t],graph,params,t,rand_gen)
+    end
+   
+    return solution
+end
+# function solve_2!(u_0,graph,params,steps,rand_gen)
+#     S_vec = findall(x -> x == Susceptible, u_0)
+#     solution = vcat((u_0,S_vec), [(fill(Susceptible,length(u_0)), deepcopy(S_vec)) for i in 1:steps])
+#     for t in 1:steps
+#         agents_step_2!(solution[t+1],solution[t],graph,params,t,rand_gen)
+#     end
+   
+#     return solution
+# end
+function get_parameters()
+    params = (
+        p = 0.9,
+        recovery_rate = 0.1,
+    )
+    return params
+end
+
+function fit_parameters()
+
+
+end
+
+

src/plotting.jl

+
+function plot_network(g,solution,weights) #it's shockingly hard to plot animated networks in julia so this function does that...
+    a = LightGraphs.adjacency_matrix(g) 
+    network_pts = SFDP.layout(a,SFDP.Point2f0;tol=0.1,C=1,K=1,iterations=10) # generate 2D layout
+    # display(network_pts)
+    xpts = [i[1] for i in network_pts]
+    ypts = [i[2] for i in network_pts]
+    plt = plot()
+    for edge in edges(g)
+        e = [src(edge), dst(edge)]
+        plot!(plt,xpts[e],ypts[e], color = :black)
+    end
+    anim = Animation()
+    for (t,u_t) in enumerate(solution)
+        frame_t = deepcopy(plt)
+        for (i,var) in enumerate(instances(AgentStatus))
+            indices_of_var = findall(x -> x == var, u_t)
+            scatter!(frame_t,xpts[indices_of_var],ypts[indices_of_var]; color = color_palette[i], legend = false,markersize = 5)
+        end
+        plot!(frame_t;title = "t = $t")
+        frame(anim,frame_t)
+    end
+        
+    gif(anim, "graphplot.gif", fps = 1.5)
+end
\ No newline at end of file

src/types.jl

+#Agent type definitions
+
+
+@enum AgentDemographic begin
+    Young = 1
+    Adult
+    AdultHCW
+    Elderly
+    ElderlyLTC
+end
+@enum AgentStatus begin
+    Susceptible
+    Infected
+    Recovered
+end
+
+struct AgentList
+    demographic::Vector{AgentDemographic}
+    status::Vector{AgentStatus}
+end
+
+function Base.show(io::IO, status::AgentStatus) 
+    if status == Susceptible
+        print(io, "S")
+    elseif status == Infected
+        print(io, "I")
+    elseif status == Recovered
+        print(io, "R")
+    end
+end
+
+function Base.show(io::IO, status::AgentDemographic) 
+    if status == Young
+        print(io, "<20")
+    elseif status == Adult
+        print(io, "20-60")
+    elseif status == Adult
+        print(io, "20-60,HCW")
+    elseif status == Elderly
+        print(io, ">60")
+    end
+end 
+