From 203b7a3d5bc891cdd934cfcc5c266558d4588686 Mon Sep 17 00:00:00 2001
From: pjentsch <pjentsch@uwaterloo.ca>
Date: Thu, 25 Mar 2021 21:07:18 -0400
Subject: [PATCH] fixed a breakage moving mixing distributions to
IntervalsModel from CovidAlterVaccinationModel
---
CovidAlertVaccinationModel/src/data.jl | 65 -------------------
.../src/mixing_distributions.jl | 43 ------------
IntervalsModel/src/IntervalsModel.jl | 5 +-
IntervalsModel/src/rest_durations_model.jl | 10 ++-
IntervalsModel/src/ws_durations_model.jl | 29 ++++++++-
5 files changed, 38 insertions(+), 114 deletions(-)
diff --git a/CovidAlertVaccinationModel/src/data.jl b/CovidAlertVaccinationModel/src/data.jl
index f472779..abebc66 100644
--- a/CovidAlertVaccinationModel/src/data.jl
+++ b/CovidAlertVaccinationModel/src/data.jl
@@ -10,15 +10,6 @@ function get_canada_demographic_distribution()::Vector{Float64}
source_bins = map(parse_string_as_float_pair, df[:,:demographic_data_bins])
return [sum(binned_data[1:5]),sum(binned_data[6:13]),sum(binned_data[14:end])]
end
-function get_distancing_data_sets()::NTuple{3,Vector{Float64}}
- f = readdlm(joinpath(PACKAGE_FOLDER,"data/csv/distancing_data.csv"), ',')
- df = filter(row -> row[:country_region_code] == "CA" && row[:sub_region_1] == "Ontario", DataFrame([f[1,i] => f[2:end, i] for i = 1:length(f[1,:])]))
- distancing_data = df[:,:retail_and_recreation_percent_change_from_baseline] ./ 100
- workplace_data = df[:,:workplaces_percent_change_from_baseline] ./ 100
- residental_data = df[:,:residential_percent_change_from_baseline] ./ 100
-
- return (workplace_data, distancing_data,residental_data)
-end
function get_canada_case_fatality()::Tuple{Vector{Tuple{Float64,Float64}},Vector{Float64}}
f = readdlm(joinpath(PACKAGE_FOLDER,"data/csv/case_fatality_data.csv"), ',')
@@ -62,62 +53,6 @@ function get_household_data_proportions()
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
using CSV
function load_mixing_matrices()
diff --git a/CovidAlertVaccinationModel/src/mixing_distributions.jl b/CovidAlertVaccinationModel/src/mixing_distributions.jl
index 3fc3aba..60ad4a0 100644
--- a/CovidAlertVaccinationModel/src/mixing_distributions.jl
+++ b/CovidAlertVaccinationModel/src/mixing_distributions.jl
@@ -13,24 +13,6 @@ function adjust_distributions_mean(distribution_matrix::AbstractArray{T},mean_sh
end
end
-
-function ZeroPoisson(α,λ)
- return ZWDist(α,Poisson(λ))
-end
-function ZeroGeometric(α,p)
- return ZWDist(α,Geometric(p))
-end
-
-function symmetrize_means(population_sizes,mixing_matrix::Matrix{<:ZWDist})
- mixing_means = mean.(mixing_matrix)
- symmetrized_means = zeros((length(population_sizes),length(population_sizes)))
- for i in 1:length(population_sizes), j in 1:length(population_sizes)
- symmetrized_means[i,j] = 0.5*(mixing_means[i,j] + population_sizes[j]/population_sizes[i] * mixing_means[j,i])
- end
- return map(((dst,sym_mean),)->from_mean(typeof(dst),dst.α,sym_mean), zip(mixing_matrix,symmetrized_means))
-end
-
-
function symmetrize_means(population_sizes,mixing_matrix::Matrix{<:Distribution})
mixing_means = mean.(mixing_matrix)
symmetrized_means = zeros((length(population_sizes),length(population_sizes)))
@@ -72,31 +54,6 @@ const unemployment_matrix = alpha_matrix(
)
)
-function make_workschool_mixing_matrix()
-
- #all geometric with means given
- ws_mixing = map(t->from_mean(t...),[
- (Geometric{Float64}, 4.104848) (Geometric{Float64},2.568782) (Geometric{Float64},0.017729)
- (Geometric{Float64}, 0.975688) (Geometric{Float64},5.057572) (Geometric{Float64},0.021307)
- (Geometric{Float64},0.001937) (Geometric{Float64},0.00722) (Geometric{Float64}, 0.022134)
- ])
-
- #symmetrize WS mixing with respect to the population proportions in the data
- ws_mixing_w_unemployment_symmetrized = symmetrize_means(get_household_data_proportions(),ws_mixing)
-
- #define a function that adjusts the means of W according to the unemployment_matrix
- ws_adjust_mean(W) = (5/7) .* (1 .- unemployment_matrix) ./ ( mean.(W) + (5/7) .* (1 .- unemployment_matrix))
-
- #create a zero weighted distribution where the zero-weight is given by the unemployment_matrix, and the non-zero weight is given by ws_mixing, symmetrized, with means adjust upwards
- ws_mixing_w_unemployment_symmetrized_weekday_adjusted = ZWDist.(unemployment_matrix,Geometric.(ws_adjust_mean(ws_mixing_w_unemployment_symmetrized)))
-
- return ws_mixing_w_unemployment_symmetrized_weekday_adjusted
-end
-
-@views function ws_sample(age)
- return rand.(initial_workschool_mixing_matrix[age,:]) * (rand(RNG) < (5/7))
-end
-
const workschool_mixing, rest_mixing = load_mixing_matrices()
diff --git a/IntervalsModel/src/IntervalsModel.jl b/IntervalsModel/src/IntervalsModel.jl
index d0811f8..936a9d8 100644
--- a/IntervalsModel/src/IntervalsModel.jl
+++ b/IntervalsModel/src/IntervalsModel.jl
@@ -10,7 +10,6 @@ using Random
using KernelDensity
using StatsBase
using Distributions
-using CovidAlertVaccinationModel
using ZeroWeightedDistributions
const PACKAGE_FOLDER = dirname(dirname(pathof(IntervalsModel)))
using DataStructures:OrderedDict
@@ -58,7 +57,7 @@ end
priors_list::Vector{Distribution},
dists::Vector{Distribution},
particles::Integer;
- alpha = 0.95
+ alpha = 0.98
)
# Arguments
@@ -85,7 +84,7 @@ Number of particles to obtain representing posterior distribution of the distrib
# Optional Arguments
-- alpha = 0.95
+- alpha = 0.98
Threshold adaptive rate, see `KissABC.smc` for more details.
"""
diff --git a/IntervalsModel/src/rest_durations_model.jl b/IntervalsModel/src/rest_durations_model.jl
index c0f3203..8b3c7c2 100644
--- a/IntervalsModel/src/rest_durations_model.jl
+++ b/IntervalsModel/src/rest_durations_model.jl
@@ -1,5 +1,3 @@
-
-const rest_distributions = CovidAlertVaccinationModel.initial_rest_mixing_matrix
const rest_data = get_rest_data()
const kerneldensity_data_rest = (
Y = InterpKDE(kde(rest_data.Y.durs,weights = rest_data.Y.weights)),
@@ -11,6 +9,14 @@ const data_rest = hcat(
map(i->pdf(kerneldensity_data_rest.M,i), 0:143),
map(i->pdf(kerneldensity_data_rest.O,i), 0:143),
)
+using CovidAlertVaccinationModel: symmetrize_means, get_household_data_proportions, from_mean
+const rest_distributions = symmetrize_means(get_household_data_proportions(),map(t->from_mean(t...),[
+ (Geometric{Float64},2.728177) (Geometric{Float64},1.382557) (Geometric{Float64},0.206362)
+ (Geometric{Float64},1.139072) (Geometric{Float64},3.245594) (Geometric{Float64},0.785297)
+ (Geometric{Float64},0.264822) (Geometric{Float64},0.734856) (Geometric{Float64},0.667099)
+]))
+
+
"""
err_rest(p::Vector{Number},dist::Distribution)
diff --git a/IntervalsModel/src/ws_durations_model.jl b/IntervalsModel/src/ws_durations_model.jl
index c9d8de2..404e127 100644
--- a/IntervalsModel/src/ws_durations_model.jl
+++ b/IntervalsModel/src/ws_durations_model.jl
@@ -33,7 +33,7 @@ function err_ws(p,dist)
for _ in 1:sample_repeat
for age_sample in YOUNG:OLD
#Sample neighbourhood sizes. Note that this includes the weekend/weekday coinflip, implemented in CovidAlertVaccinationModel.
- neighourhoods = CovidAlertVaccinationModel.ws_sample(age_sample)
+ neighourhoods = ws_sample(age_sample)
for age_j in YOUNG:OLD #for a given age_sample loop over possible contact ages
if neighourhoods[age_j] > 0
@@ -57,3 +57,30 @@ function err_ws(p,dist)
end
return errsum./sample_repeat
end
+
+
+
+function make_workschool_mixing_matrix()
+
+ #all geometric with means given
+ ws_mixing = map(t->from_mean(t...),[
+ (Geometric{Float64}, 4.104848) (Geometric{Float64},2.568782) (Geometric{Float64},0.017729)
+ (Geometric{Float64}, 0.975688) (Geometric{Float64},5.057572) (Geometric{Float64},0.021307)
+ (Geometric{Float64},0.001937) (Geometric{Float64},0.00722) (Geometric{Float64}, 0.022134)
+ ])
+
+ #symmetrize WS mixing with respect to the population proportions in the data
+ ws_mixing_w_unemployment_symmetrized = symmetrize_means(get_household_data_proportions(),ws_mixing)
+
+ #define a function that adjusts the means of W according to the unemployment_matrix
+ ws_adjust_mean(W) = (5/7) .* (1 .- unemployment_matrix) ./ ( mean.(W) + (5/7) .* (1 .- unemployment_matrix))
+
+ #create a zero weighted distribution where the zero-weight is given by the unemployment_matrix, and the non-zero weight is given by ws_mixing, symmetrized, with means adjust upwards
+ ws_mixing_w_unemployment_symmetrized_weekday_adjusted = ZWDist.(unemployment_matrix,Geometric.(ws_adjust_mean(ws_mixing_w_unemployment_symmetrized)))
+
+ return ws_mixing_w_unemployment_symmetrized_weekday_adjusted
+end
+
+@views function ws_sample(age)
+ return rand.(initial_workschool_mixing_matrix[age,:]) * (rand(RNG) < (5/7))
+end
--
GitLab