Skip to content

Instantly share code, notes, and snippets.

@affans

affans/covid19 - basic dynamical system

Last active March 9, 2021 17:00
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save affans/f529f8cd727989143c979afd9174b984 to your computer and use it in GitHub Desktop.
Save affans/f529f8cd727989143c979afd9174b984 to your computer and use it in GitHub Desktop.
Download ZIP
COVID-19 Model Implemention
Raw
covid19 - basic dynamical system
# Affan Shoukat
# Center for Infectious Disease Modelling and Analysis
# Yale University, Dec 2020
# Repository: https://gist.github.com/affans/f529f8cd727989143c979afd9174b984
module covid19ode
using OrdinaryDiffEq
using Parameters
using LinearAlgebra
using ForwardDiff
using ProgressMeter
using DelimitedFiles
using Distributions
using Bootstrap
using JLD2
using FileIO
using Distributed
using Gnuplot
#using DataFrames
#using Interpolations
#using RCall
#using ComponentArrays
# julia> using Distributed; addprocs(4, exeflags="--project=.")
# julia> @everywhere include("model.jl")
# julia> using Revise; Revise.track("model.jl")
# julia> @everywhere using .covid19ode
# par init
# using ClusterManagers, Distributed;
# addprocs(SlurmManager(288), N=9, topology=:master_worker, exeflags="--project=.")
# include("model.jl");
# @everywhere include("model.jl")
# define helper functions
∑(x) = round(sum(x); digits=0)
H(x) = x <= 0 ? 0 : 1
pinc(final, orig) = (final - orig) / abs(orig)
pdec(final, orig) = (orig - final) / abs(orig)
sevenaverage(X) = movingaverage(X, 7)
function movingaverage(X,numofele::Int)
BackDelta = div(numofele,2)
ForwardDelta = isodd(numofele) ? div(numofele,2) : div(numofele,2) - 1
len = length(X)
Y = zeros(Float64, len)
for n = 1:len
lo = max(1,n - BackDelta)
hi = min(len,n + ForwardDelta)
Y[n] = mean(X[lo:hi])
end
return Int.(round.(Y))
end
export ∑, H, pinc, pdec, sevenaverage, movingaverage
@with_kw mutable struct ModelParameters
β::Float64 = 0.0267 # transmission parameter, calibrated
β₊::Float64 = 0.0 # step size to increase beta by, set in model setup
α::Float64 = 0.75 # rel asymptomatic transmission, Sayampanathan Lancet, see new paper: doi:10.1001/jamanetworkopen.2020.35057
ϵ₁::NTuple{5, Float64} = 0.46 .* (0.0, 0.0, 1, 1, 1) # vaccine efficacy, pfizer and moderna papers
ϵ₂::NTuple{5, Float64} = 0.60 .* (0.0, 0.0, 1, 1, 1) # vaccine efficacy, pfizer and moderna papers
λ::Array{Int64,1} = [0, 0, 0, 0, 0] # see vaccine function
s_vac::Float64 = 0 # flags for susceptible -> V1 or V1 -> V2 with rate lambda, also controls vaccine on and off
v_vac::Float64 = 0
covreached::Array{Bool,1} = [0, 0, 1, 1, 1] # works "opposite". 1 means coverage NOT reached, probably better to rename
p::NTuple{5, Float64} = (0.3, 0.377, 0.328, 0.328, 0.188) # proportion of non-vaccinated asymptomatic, Buitrago-Garcia PLOS Medicine
ρ₁::NTuple{5, Float64} = 0.66 .* (1.0, 1.0, 1.0, 1.0, 1.0) # see vaccine function
ρ₂::NTuple{5, Float64} = 0.94 .* (1.0, 1.0, 1.0, 1.0, 1.0)
σ::Float64 = 1/3.35 # average latent period incubation period - the presymptomatic period
η::Float64 = 1/5 # average asymptomatic period # Li, Moghadas
θ::Float64 = 1/1.85 # average presymptomatic period
γ::Float64 = 1/3.4 # average symptomatic period
δ::Float64 = 1/1.325 # days before identification (in silent infectious stage) (0.8 to 2.8 days), calculated using min of asymp and pre (then Uniform from 0.8 to this minimum, see sample function)
τ::Float64 = 1/1.0 # days before identification (in symptomatic period) CDC https://www.cdc.gov/coronavirus/2019-ncov/hcp/planning-scenarios.html
κ::NTuple{5, Float64} = (1/6.0, 1/6.0, 1/6.0, 1/6.0, 1/4.0) # days before hospitalization for symptomatic patients CDC https://www.cdc.gov/coronavirus/2019-ncov/hcp/planning-scenarios.html
νᵣ::NTuple{5, Float64} = (1/3, 1/3, 1/3, 1/4, 1/6) # days spent in hospitalization before death
νₓ::NTuple{5, Float64} = (1/11, 1/11, 1/11, 1/14, 1/12) # days spent in hospitalization before death (assuming ICU) https://www.cdc.gov/coronavirus/2019-ncov/hcp/planning-scenarios.html
h::NTuple{5, Float64} = (0.006, 0.006, 0.217, 0.257, 0.512) # proportion of symptomatic individuals that will be hospitalized, https://gis.cdc.gov/grasp/COVIDNet/COVID19_5.html, data downloaded december 3
d::NTuple{5, Float64} = (0.009, 0.007, 0.031, 0.105, 0.253) # proportion of deaths out of the hospitalized individuals, see downloaded screenshot and folder
q::NTuple{5, Float64} = (0.0, 0.0, 0.0, 0.0, 0.0) # isolation during latent period (source, see function)
q₊::Float64 = 0.0 # the increase to q tuple at specified time
g::NTuple{5, Float64} = (0.0, 0.0, 0.0, 0.0, 0.0) # isolation during silent infectious stage (source, see function)
g₊::Float64 = 0.0 # the increase to g tuple at a specified time
f::NTuple{5, Float64} = (0.0, 0.0, 0.0, 0.0, 0.0) # isolation during symptomatic period
# state specific parameters
pop::NTuple{5, Float64} = (0, 0, 0, 0, 0) # age groups: 0-5, 6-19, 20-49, 50-64, 65+
ic_inf::NTuple{5, Float64} = (0, 0, 0, 0, 0) # this fits to texas data
preexist_imm::Float64 = 0.143 # default value see JAMA Swerdlow paper
flucov::NTuple{5, Float64} = (0, 0, 0, 0, 0)
vacslopeval::Float64 = 0.0
vacagedist::NTuple{5, Float64} = (0, 0, 0, 0, 0)
# interal system params
kffcov::NTuple{5, Float64} = (0.755, 0.602, 0.72, 0.74, 0.85)
npi_start::Int16 = -1
npi_end::Int16 = -1
scr_start::Int16 = -1
end
export ModelParameters
struct PlotData_t5
# the tuples are for avg, lo, hi
t::Vector{Float64} # time
W::NTuple{3, Vector{Float64}} # cumulative incidence
inc::NTuple{3, Vector{Float64}}
Z::NTuple{3, Vector{Float64}} # asymptomatic
Y::NTuple{3, Vector{Float64}} # symptomatic
U::NTuple{3, Vector{Float64}} # hospitalized
D::NTuple{3, Vector{Float64}} # deaths
end
PlotData_t5(n) = PlotData_t5(zeros(Float64, n), [Tuple([zeros(Float64, n) for _= 1:3]) for _=1:(fieldcount(PlotData_t5) - 1)]...)
export PlotData_t5
function contact_matrix()
## contact matrix for general population and in household.
## from Mossong
# 0 - 4, 5 - 19, 20 - 49, 50 - 65, 65 +
Moss = ones(Float64, 5, 5)
Moss[1, :] = [0.2287 0.1839 0.4219 0.1116 0.0539]
Moss[2, :] = [0.0276 0.5964 0.2878 0.0591 0.0291]
Moss[3, :] = [0.0376 0.1454 0.6253 0.1423 0.0494]
Moss[4, :] = [0.0242 0.1094 0.4867 0.2723 0.1074]
Moss[5, :] = [0.0207 0.1083 0.4071 0.2193 0.2446]
r = [10.21, 16.793, 13.795, 11.2669, 8.0027];
r̃ = [2.86, 4.70, 3.86, 3.15, 2.24];
M = Moss .* r
M̃ = Moss .* r̃
return M, M̃
end
function aModel!(du, u, params, t)
# internal compartments: W
S₁, S₂, S₃, S₄, S₅,
V₁₁, V₁₂, V₁₃, V₁₄, V₁₅,
V₂₁, V₂₂, V₂₃, V₂₄, V₂₅,
E₁, E₂, E₃, E₄, E₅,
F₁₁, F₁₂, F₁₃, F₁₄, F₁₅,
F₂₁, F₂₂, F₂₃, F₂₄, F₂₅,
Ẽ₁, Ẽ₂, Ẽ₃, Ẽ₄, Ẽ₅,
F̃₁₁, F̃₁₂, F̃₁₃, F̃₁₄, F̃₁₅,
F̃₂₁, F̃₂₂, F̃₂₃, F̃₂₄, F̃₂₅,
A₁, A₂, A₃, A₄, A₅,
H̃₁, H̃₂, H̃₃, H̃₄, H̃₅,
Ã₁, Ã₂, Ã₃, Ã₄, Ã₅,
P₁, P₂, P₃, P₄, P₅,
B̃₁, B̃₂, B̃₃, B̃₄, B̃₅,
P̃₁, P̃₂, P̃₃, P̃₄, P̃₅,
I₁, I₂, I₃, I₄, I₅,
Q̃₁, Q̃₂, Q̃₃, Q̃₄, Q̃₅,
Ĩ₁, Ĩ₂, Ĩ₃, Ĩ₄, Ĩ₅,
X₁, X₂, X₃, X₄, X₅,
D₁, D₂, D₃, D₄, D₅,
R₁, R₂, R₃, R₄, R₅,
W₁, W₂, W₃, W₄, W₅,
Z₁, Z₂, Z₃, Z₄, Z₅,
Y₁, Y₂, Y₃, Y₄, Y₅,
V₁, V₂, V₃, V₄, V₅,
U₁, U₂, U₃, U₄, U₅ = u
S = [S₁, S₂, S₃, S₄, S₅]
V₁ = [V₁₁, V₁₂, V₁₃, V₁₄, V₁₅]
V₂ = [V₂₁, V₂₂, V₂₃, V₂₄, V₂₅]
E = [E₁, E₂, E₃, E₄, E₅]
F₁ = [F₁₁, F₁₂, F₁₃, F₁₄, F₁₅]
F₂ = [F₂₁, F₂₂, F₂₃, F₂₄, F₂₅]
Ẽ = [Ẽ₁, Ẽ₂, Ẽ₃, Ẽ₄, Ẽ₅]
F̃₁ = [F̃₁₁, F̃₁₂, F̃₁₃, F̃₁₄, F̃₁₅]
F̃₂ = [F̃₂₁, F̃₂₂, F̃₂₃, F̃₂₄, F̃₂₅]
A = [A₁, A₂, A₃, A₄, A₅]
H̃ = [H̃₁, H̃₂, H̃₃, H̃₄, H̃₅]
à = [Ã₁, Ã₂, Ã₃, Ã₄, Ã₅]
P = [P₁, P₂, P₃, P₄, P₅]
B̃ = [B̃₁, B̃₂, B̃₃, B̃₄, B̃₅]
P̃ = [P̃₁, P̃₂, P̃₃, P̃₄, P̃₅]
I = [I₁, I₂, I₃, I₄, I₅]
Q̃ = [Q̃₁, Q̃₂, Q̃₃, Q̃₄, Q̃₅]
Ĩ = [Ĩ₁, Ĩ₂, Ĩ₃, Ĩ₄, Ĩ₅]
X = [X₁, X₂, X₃, X₄, X₅]
D = [D₁, D₂, D₃, D₄, D₅]
R = [R₁, R₂, R₃, R₄, R₅]
W = [W₁, W₂, W₃, W₄, W₅]
Z = [Z₁, Z₂, Z₃, Z₄, Z₅]
Y = [Y₁, Y₂, Y₃, Y₄, Y₅]
V = [V₁, V₂, V₃, V₄, V₅]
U = [U₁, U₂, U₃, U₄, U₅]
S′ = view(du, 1:5)
V₁′ = view(du, 6:10)
V₂′ = view(du, 11:15)
E′ = view(du, 16:20)
F₁′ = view(du, 21:25)
F₂′ = view(du, 26:30)
Ẽ′ = view(du, 31:35)
F̃₁′ = view(du, 36:40)
F̃₂′ = view(du, 41:45)
A′ = view(du, 46:50)
H̃′ = view(du, 51:55)
Ã′ = view(du, 56:60)
P′ = view(du, 61:65)
B̃′ = view(du, 66:70)
P̃′ = view(du, 71:75)
I′ = view(du, 76:80)
Q̃′ = view(du, 81:85)
Ĩ′ = view(du, 86:90)
X′ = view(du, 91:95)
D′ = view(du, 96:100)
R′ = view(du, 101:105)
W′ = view(du, 106:110)
Z′ = view(du, 111:115)
Y′ = view(du, 116:120)
V′ = view(du, 121:125)
U′ = view(du, 126:130)
M, M̃ = contact_matrix()
@unpack β, α, ϵ₁, ϵ₂, λ, s_vac, v_vac, covreached, p, ρ₁, ρ₂, σ, η, θ, γ, δ, τ, κ, νᵣ, νₓ, q, g, f, h, d, pop, flucov, npi_start = params
for a = 1:5
# define force of infection
ℱ = β * (α*dot(M[a, :], A./pop) + dot(M[a, :], P./pop) + dot(M[a, :], I./pop) + α*dot(M̃[a, :], Ã./pop) + α*dot(M̃[a, :], H̃./pop) +
dot(M̃[a, :], P̃./pop) + dot(M̃[a, :], B̃./pop) + dot(M̃[a, :], Ĩ./pop) + dot(M̃[a, :], Q̃./pop))
S′[a] = -S[a]*ℱ - H(S[a] - λ[a])*λ[a]*s_vac*covreached[a]
V₁′[a] = H(S[a] - λ[a])*λ[a]*s_vac*covreached[a] - H(V₁[a] - λ[a])*λ[a]*v_vac - (1 - ϵ₁[a])*V₁[a]*ℱ
V₂′[a] = H(V₁[a] - λ[a])*λ[a]*v_vac - (1 - ϵ₂[a])*V₂[a]*ℱ
E′[a] = (1 - q[a])*S[a]*ℱ - σ*E[a]
F₁′[a] = (1 - q[a])*(1 - ϵ₁[a])*V₁[a]*ℱ - σ*F₁[a]
F₂′[a] = (1 - q[a])*(1 - ϵ₂[a])*V₂[a]*ℱ - σ*F₂[a]
#F′[a] = (1 - q[a])*((1 - ϵ₁[a])*V₁[a]*ℱ + (1 - ϵ₂[a])*V₂[a]*ℱ) - σ*F[a]
Ẽ′[a] = q[a]*S[a]*ℱ - σ*Ẽ[a]
F̃₁′[a] = q[a]*(1 - ϵ₁[a])*V₁[a]*ℱ - σ*F̃₁[a]
F̃₂′[a] = q[a]*(1 - ϵ₂[a])*V₂[a]*ℱ - σ*F̃₂[a]
#F̃′[a] = q[a]*((1 - ϵ₁[a])*V₁[a]*ℱ + (1 - ϵ₂[a])*V₂[a]*ℱ) - σ*F̃[a]
A′[a] = p[a]*σ*E[a] + ρ₁[a]*σ*F₁[a] + ρ₂[a]*σ*F₂[a] - (1 - g[a])*η*A[a] - g[a]*δ*A[a]
H̃′[a] = p[a]*σ*Ẽ[a] + ρ₁[a]*σ*F̃₁[a] + ρ₂[a]*σ*F̃₂[a] - η*H̃[a]
Ã′[a] = g[a]*δ*A[a] - (δ*η / (δ - η))*Ã[a]
P′[a] = (1 - p[a])*σ*E[a] + (1 - ρ₁[a])*σ*F₁[a] + (1 - ρ₂[a])*σ*F₂[a] - (1 - g[a])*θ*P[a] - g[a]*δ*P[a]
B̃′[a] = (1 - p[a])*σ*Ẽ[a] + (1 - ρ₁[a])*σ*F̃₁[a] + (1 - ρ₂[a])*σ*F̃₂[a] - θ*B̃[a]
P̃′[a] = g[a]*δ*P[a] - (δ*θ / (δ - θ))*P̃[a]
I′[a] = (1 - g[a])*θ*P[a] - h[a]*κ[a]*I[a] - (1 - h[a])*(1 - f[a])*γ*I[a] - (1 - h[a])*f[a]*τ*I[a]
Q̃′[a] = θ*B̃[a] + (δ*θ / (δ - θ))*P̃[a] - h[a]*κ[a]*Q̃[a] - (1 - h[a])*γ*Q̃[a]
Ĩ′[a] = (1 - h[a])*f[a]*τ*I[a] - (τ*γ / (τ - γ))*Ĩ[a]
X′[a] = h[a]*κ[a]*I[a] + h[a]*κ[a]*Q̃[a] - (1 - d[a])*νᵣ[a]*X[a] - d[a]*νₓ[a]*X[a]
D′[a] = d[a]*νₓ[a]*X[a]
R′[a] = (1 - g[a])*η*A[a] + η*H̃[a] + (δ*η / (δ - η))*Ã[a] + (1 - h[a])*(1 - f[a])*γ*I[a] +
(1 - h[a])*γ*Q̃[a] + (τ*γ / (τ - γ))*Ĩ[a] + (1 - d[a])*νᵣ[a]*X[a]
# internal compartments
W′[a] = S[a]*ℱ + (1 - ϵ₁[a])*V₁[a]*ℱ + (1 - ϵ₂[a])*V₂[a]*ℱ # incidence
Z′[a] = p[a]*σ*E[a] + ρ₁[a]*σ*F₁[a] + ρ₂[a]*σ*F₂[a] + p[a]*σ*Ẽ[a] + ρ₁[a]*σ*F̃₁[a] + ρ₂[a]*σ*F̃₂[a] # asymptomatic
Y′[a] = (1 - p[a])*σ*E[a] + (1 - ρ₁[a])*σ*F₁[a] + (1 - ρ₂[a])*σ*F₂[a] + (1 - p[a])*σ*Ẽ[a] + (1 - ρ₁[a])*σ*F̃₁[a] + (1 - ρ₂[a])*σ*F̃₂[a] # symptomatic
V′[a] = H(S[a] - λ[a])*λ[a]*s_vac*covreached[a] # vaccinated
U′[a] = h[a]*κ[a]*I[a] # hospitalizations
end
end
# doses per day
# we have daily number of doses given up till february 11 (or 55 data points)
# so while this condition is true, change the number of doses per day
dpd_cond(u, t, integrator) = (t <= 55) ? true : false
function dpd_change!(integrator)
prms = integrator.p
slopevalue = prms.vacslopeval
doses_per_day = 9.9 + slopevalue*integrator.t
pop = prms.pop
totaladults = sum(pop[3:5])
props = round.(doses_per_day .* prms.vacagedist; digits = 0)
#integrator.p.λ = [0, 0, props...]
integrator.p.λ = [props...]
end
# flu vax coverage
# if an age group reaches flu vax conditions, the available doses are shifted to the next age group
function fluvax_condition(u, t, integrator)
# check if any of age groups recieved their maximum coverage
# if so, change the lambda allocation. Note that when t < 55,
# lambda is overwritten by the dpd_change! condition.
# but from my empirical tests, given the maximum doses per day, we don't hit maximum numbers till after day 55.
fvp = integrator.p.flucov
for i = 3:5 ## at most everything allocates to age group 3
# u[i] is the i'th compartment of the u solution vector, so u[121] is the internal V1 class for the first age group.
hval = H(fvp[i] - u[i+120]) == 0
covval = integrator.p.covreached[i] == 1
retval = hval && covval
if retval
#println("condition hit --t: $t, ag: $i, max:$(fvp[i]), current: $(u[i+120])")
return true
end
#retval && return true
end
return false
# can also use `any`: return any(calc_true_false_value() for i in 3:5)
# returns as soon as there is a true value
end
function fluvax_change!(integrator)
# lets figure out what age group has recieved their maximum vaccine coverage, and change the allocation
fvp = integrator.p.flucov
for i = 3:5 ## at most everything allocates to age group 3
#println("ag: $i $(integrator.p.λ), time: $(integrator.t), u: $(integrator.u[i+120]), flucov: $(fvp[i])")
hval = H(fvp[i] - integrator.u[i+120]) == 0
if hval
integrator.p.covreached[i] = 0
#integrator.p.λ[i - 1] += integrator.p.λ[i]
end
end
end
# switching between first and second doses
function dose_switch_condition(u,t,integrator)
t == 21.0 # at time 21, start the continous dosage
end
function dose_switch!(integrator)
integrator.p.s_vac = 0.70
integrator.p.v_vac = 0.30
end
# when to start lifting behavioural NPI
function npilift_condition(u, t, integrator)
integrator.p.npi_start <= t <= (integrator.p.npi_start + integrator.p.npi_end)
end
function npilift_start!(integrator)
# this condition corresponds to behavioural NPI - increase beta by x% over npi_end days
# calculation of beta+ happens at the modelparams setup (`run_model`)
integrator.p.β += integrator.p.β₊
end
function screen_condition(u, t, integrator)
integrator.p.scr_start == t
end
function screen_start!(integrator)
integrator.p.g = integrator.p.g .+ integrator.p.g₊
end
function _simulate(params)
#set up initial conditions and run the model
tspan = (0.0, 365.0)
total_pop = params.pop
preexist_imm = params.preexist_imm
u0 = zeros(Float64, 130)
u0[1:5] .= (1 - preexist_imm) .* (total_pop) ## susceptible
u0[16:20] = [params.ic_inf[1], params.ic_inf[2], params.ic_inf[3], params.ic_inf[4], params.ic_inf[5]] # one exposed individual in each of the E[a] class.
u0[101:105] .= preexist_imm .* (total_pop)
cb_dpd = DiscreteCallback(dpd_cond, dpd_change!, save_positions=(false,false))
cb_fluvaxcov = DiscreteCallback(fluvax_condition, fluvax_change!, save_positions=(false,false))
cb_doseswitch = DiscreteCallback(dose_switch_condition, dose_switch!, save_positions=(false,false))
cb_npi = DiscreteCallback(npilift_condition, npilift_start!, save_positions=(false,false))
cb_screen = DiscreteCallback(screen_condition, screen_start!, save_positions=(false,false))
# add model callbacks
cbset = CallbackSet()
params.scr_start > 0 && (cbset = CallbackSet(cbset, cb_screen))
if params.npi_start > 0 && params.npi_end > 0
cbset = CallbackSet(cbset, cb_npi)
end
params.s_vac > 0 && (cbset = CallbackSet(cbset, cb_dpd, cb_fluvaxcov, cb_doseswitch))
prob = ODEProblem(aModel!, u0, tspan, params)
sol = solve(prob, Rodas4(autodiff=false), dt=1, adaptive=false, callback=cbset)
return sol
end
export _simulate
@inline function sample_dis_params(params)
# sample relevant disease specific parameters
# to get the latent period: first sample the incubation period, then sample the presymptomatic period
# the difference between them is the latent period.
# other sample points: asymptomatic and symptomatic periods
# potential other sample points, f, g, q.
# also other potential sample points are days until moved out of compartment
# note: tau has to be smaller than gamma (which is already true by truncation)
# note: delta has to be smaller than theta and eta (but also in the range 0.8 to 2.8)
inc_per_dist = truncated(LogNormal(1.434, 0.661), 4, 7) # gives mean 5.2, as reference
pre_per_dist = truncated(Gamma(1.058, 2.174), 1, 3) # gives mean 1.85 (but the actual distribution gives mean 2.3)
asy_per_dist = truncated(Gamma(5, 1), 1, Inf) # mean 5 days as reference
sym_per_dist = truncated(Gamma(2.768, 1.1563), 1.1, Inf) # mean 3.4 days (original distribution 3.2)
# plot the histogram for inc_period
# incs = [rand(truncated(LogNormal(1.434, 0.661), 4, 7)) for _=1:10000]
# println(mean(incs))
# h = hist(incs, bs=0.5)
# @gp h.bins h.counts "w impulses t 'Data' lc rgb 'red'"
# # plot the histogram for presympt period
# pres = [rand() for i = 1:10000]
# println(mean(pres))
# h = hist(pres, bs=0.05)
# @gp h.bins h.counts "w impulses t 'Data' lc rgb 'red'"
inc = rand(inc_per_dist)
pre = rand(pre_per_dist)
lat = inc - pre
asy = rand(asy_per_dist)
sym = rand(sym_per_dist)
min_pre_asy = min(2.8, min(pre, asy) - 0.1) # the -0.1 is for stability, with 2.8 days max
delta = rand(Uniform(0.8, min_pre_asy))
params.δ = 1/delta
params.σ = 1/lat # average latent period
params.η = 1/asy # average asymptomatic period # Li, Moghadas
params.θ = 1/pre # average presymptomatic period
params.γ = 1/sym # average symptomatic period
#println("inc:$inc, \n lat: $lat_sigma ($difflat), \n asymp: $asy_eta ($diffeta),\n pre: $pre_theta ($difftheta) (delta: $delta, diff: $diffdelta) \n symp: $sym_gamma ($diffgamma)")
#println("rates: lat: $(1/lat_sigma), asymp: $(1/asy_eta), pre: $(1/pre_theta), delta: $(1/delta), symp: $(1/sym_gamma)")
end
export sample_dis_params
function default_qg_params(params)
# according to https://www.nature.com/articles/d41586-020-03518-4
# In New Jersey, just 49% of cases between July and November were contacted; only 31% of those provided any contact details.
# so 0.49 contact * 0.31 actual gave information * 0.80 of individuals adhered
# split this between those quarantine right away and those captured at silent infection stage.
# note: same calculation is done in ngm2
qg_total = 0.49 * 0.31 * 0.80
params.q = Tuple((qg_total / 2 for _ = 1:5))
params.f = (0.80, 0.80, 0.80, 0.80, 0.80)
params.g = (0, 0, 0, 0, 0)
end
export default_qg_params
function setup_state_params(params, st)
# to do: setup default parameters
# tx slopevalue = 9.9, props [0.3138, 0.2408, 0.4478] 55 stopping time
# slopevalue ct = 9.9, props [0.1789, 0.2025, 0.6185] 55 stopping time
if st == "TX"
params.pop = 100000 .* (0.068661166, 0.214502674, 0.41644308, 0.171608271, 0.128784809) #TX pop
params.ic_inf = (0, 80, 130, 100, 40) .* 1.25 #80/120/120/80
params.preexist_imm = 0.163 # see the news article in notes.md
params.flucov = params.kffcov .* params.pop
params.vacslopeval = 9.9
params.vacagedist = (0, 0, 0.3138, 0.2408, 0.4478)
elseif st == "CT"
params.pop = 100000 .* (0.050966444, 0.181716647, 0.376013207, 0.214531397, 0.176772305)
params.flucov = params.kffcov .* params.pop
params.ic_inf = (0, 5, 10, 7, 4) .* 20 #80/120/120/80
params.preexist_imm = 287396/3565000
params.vacslopeval = 9.9
params.vacagedist = (0, 0, 0.1789, 0.2025, 0.6185)
else
error("Model calibration parameters not set")
end
end
export setup_state_params
function calculate_beta_ss(params, perinc, over_time)
# over_time is basically p.npi_end but maybe this function runs before npi_end is set
# calculates the step size that is added to beta
ss = 0
if over_time > 0 # can also check npi_start
ss = ((params.β*perinc) - params.β)/over_time #since half time step
end
params.β₊ = ss
end
export calculate_beta_ss
function run_model(sim_id, s_vac, rvalue, rincrease, npi_start, npi_dur, scr_start, gplus, sample_p)
# run a single model realization with vaccine, R, npi, and sampling
params = ModelParameters()
setup_state_params(params, "TX")
reff = Inf
# upper limit reff for TX: 1.85, CT: 2.2
while sample_p && reff > 1.85
sample_dis_params(params)
reff = ngm2(rvalue, params)
end
# if sample_p
# while true
# sample_dis_params(params)
# reff = ngm2(rvalue, params)
# if reff <= 1.85 #1.90 worked
# break
# end
# end
# end
params.s_vac = s_vac
params.β = rvalue
calculate_beta_ss(params, rincrease, npi_dur) # function has to run after npi_start, end set
default_qg_params(params)
params.g₊ = gplus
params.scr_start = scr_start
params.npi_start = npi_start
params.npi_end = npi_dur
sol = _simulate(params)
return sol
end
export run_model
function run_par_model(beta)
println("updated with 2.4 effective")
#baseline status quo with vaccine on
p = Progress(500, barglyphs=BarGlyphs("[=> ]"))
bl = progress_map(1:500, progress = p, mapfun = pmap) do x
run_model(x, 1, beta, 0.0, -1, -1, -1, 0.0, true)
end
bpd = create_pd(bl)
end
export run_par_model
function create_pd(sols)
# creates the plot data object from the stochastic sols
# create the PD otherwise the data file is many many gigabytes
nsims = length(sols)
tlens = [length(sols[i].t) for i = 1:nsims]
!any(diff(tlens) .== 0) && error("time vectors do not match in multiple simulations")
xvals = sols[1].t
pd = PlotData_t5(length(xvals))
pd.t .= xvals
for comp in (:W, :inc, :Z, :Y, :U, :D)
sv = zeros(Float64, length(xvals), nsims)
for i = 1:nsims
gc = get_compartments(sols[i])
cdata = getfield(gc, comp)
sv[:, i] .= cdata
end
bmean = zeros(Float64, length(xvals))
qlo = zeros(Float64, length(xvals))
qhi = zeros(Float64, length(xvals))
for i = 1:length(xvals)
## basic bootstrap
bs1 = bootstrap(mean, sv[i, :], BasicSampling(1000))
bci1 = confint(bs1, NormalConfInt(0.95));
bmean[i] = bci1[1][1]
qlo[i] = bci1[1][2]
qhi[i] = bci1[1][3]
end
getfield(pd, comp)[1] .= bmean
getfield(pd, comp)[2] .= qlo
getfield(pd, comp)[3] .= qhi
end
println("sum data incidence all: $(sum(pd.inc[1]))")
println("sum data incidence 1:80: $(sum(pd.inc[1][1:80]))")
return pd
end
export create_pd
function bin_pd(pd::PlotData_t5, bins)
# multiplies the plotdata by bin count to get true population count
for comp in fieldnames(PlotData_t5)
comp == :t && continue
getfield(pd, comp)[1] .*= bins
getfield(pd, comp)[2] .*= bins
getfield(pd, comp)[3] .*= bins
end
end
export bin_pd
function plot_pd(pd::PlotData_t5)
# internal plotting for the create_pd stochastic simulations object
println("plotting inc")
# @gp :inc "set grid" "set key right"
# @gp :- :inc "set xrange [1:365]"
# @gp :- :inc "set title 'Cumulative Incidence'"
# @gp :- :inc pd.t pd.W[2] pd.W[3] "with filledcurves notitle fc 'grey' fs solid 0.5 border lc 'blue'"
# @gp :- :inc pd.t pd.W[1] "with lines t 'cum. incidence' lw 3 lc rgb 'grey'"
# @gp :- :inc pd.t pd.Z[2] pd.Z[3] "with filledcurves notitle fc '#8dd3c7' fs solid 0.5 border lc '#8dd3c7'"
# @gp :- :inc pd.t pd.Z[1] "with lines t 'asymptomatic' lw 3 lc rgb '#8dd3c7'"
# @gp :- :inc pd.t pd.Y[2] pd.Y[3] "with filledcurves notitle fc '#bebada' fs solid 0.5 border lc '#8dd3c7'"
# @gp :- :inc pd.t pd.Y[1] "with lines t 'symptomatic' lw 3 lc rgb '#bebada'"
# @gp :- :inc "set label 'overall AR: $((pd.W[1])[end])' at graph 0.1, 0.9"
# @gp :agp "set grid" "set key right"
# @gp :- :agp "set title 'Hospitalization and Death'"
# @gp :- :agp "set xrange [1:365]"
# @gp :- :agp pd.t pd.U[2] pd.U[3] "with filledcurves notitle fc 'grey' fs solid 0.5 border lc 'blue'"
# @gp :- :agp pd.t pd.U[1] "with lines t 'hospitalized' lw 3 lc rgb 'grey'"
# @gp :- :agp pd.t pd.D[2] pd.D[3] "with filledcurves notitle fc '#8dd3c7' fs solid 0.5 border lc '#8dd3c7'"
# @gp :- :agp pd.t pd.D[1] "with lines t 'deaths' lw 3 lc rgb '#8dd3c7'"
#sn = Symbol(rand(String, 3)) # gnuplot session name
@gp :sn "set grid" :-
@gp :- :sn "set xrange [1:365]" :-
@gp :- :sn pd.t pd.inc[1] "with lines notitle lw 3 lc rgb '#66c2a5'"
Gnuplot.save(:mdp, term="pdfcairo enhanced font 'Arial,14' size 8.5,4.5", output="vaccine_fitting.pdf")
#Gnuplot.save(:mdp, term="pngcairo enhanced font 'Arial,10' size 700,400", output="vaccine_fitting.png")
println("plotting finished")
#@gp :- sn "set xtics ('Dec 14' 1, 'Jan' 31, 'Feb' 61, 'Mar' 91, 'Apr' 121, 'May' 151, 'Jun' 181, 'Jul' 211, 'Aug' 241, 'Sep' 271, 'Oct' 301, 'Nov' 301, 'Dec' 331)" :-
#@gp :- sn bl.t[1:60] bl.inc[2] bl.inc[3] "with filledcurves notitle fc '#66c2a5' fs solid 0.5 border"
end
export plot_pd
function plot_txdata(sessname)
ct_data = reverse([494, 502, 2680, 0, 0, 787, 975, 1493, 1357, 2233, 0, 0, 1198, 547, 534, 580, 2905, 0, 0, 838, 1003, 888, 869, 4367, 0, 0, 1431, 937, 482, 2568, 3931, 0, 0, 1258, 1426, 2440, 1267, 5817, 0, 0, 2019, 1662, 1915, 2094, 6703, 0, 0, 1878, 968, 3529, 3689, 7364, 0, 0, 3236, 3304, 2486, 2332, 4516, 0, 4412, 0, 2045, 1696, 767, 8457, 0, 0, 0, 2038, 1745, 1583, 4595, 0, 0, 2680, 2321, 2319, 1470, 7231])
tbindata = (ct_data ./ 35.65)
tx_data = reverse([7822, 7747, 3821, 3815, 11073, 7955, 7389, 7517, 11809, 6365, 4484, 6486, 2937, 3131, 3766, 3348, 3889, 6933, 11282, 12502, 11890, 12897, 13329, 7485, 6959, 13897, 14495, 15281, 17620, 23047, 31811, 11370, 19234, 19076, 18220, 19613, 26274, 6319, 11565, 17672, 22646, 22360, 31255, 9476, 11590, 16402, 24657, 27204, 23064, 27343, 26052, 14834, 15855, 23290, 23520, 24578, 24010, 31630, 18182, 16095, 4763, 16311, 18725, 21469, 32552, 14829, 8046, 2694, 4335, 17064, 23363, 21147, 10280, 7780, 17907, 16792, 19849, 18802, 18926, 8901])
tbindata = (tx_data ./ 289.96)
@gp :- sessname 1:80 tbindata "with line title 'Texas Data (Dec 14 - Mar 3)' lw 1 dashtype 2 lc rgb 'black'"
println("sum data incidence: $(sum(tbindata))")
end
export plot_txdata
function test_vaccine_allocation()
# tests the vaccine allocation and generates the plots of fitting and coverage
params = ModelParameters()
setup_state_params(params, "CT")
params.s_vac = 1 # turn on vaccine.
params.β = 0
sol = _simulate(params)
totalpop = sum((1 - params.preexist_imm) .* (params.pop))
# plot varaibles (only look at t = 0 to 150, so 300 time points)
cols = ["black", "black", "red", "green", "blue"]
v1ag = [get_ode_class(i, sol) for i = 6:10]
v2ag = [get_ode_class(i, sol) for i = 11:15]
# v1_int is atleast vaccinated (since they dont move out of this internal class)
v1_int = [get_ode_class(i, sol) for i = 121:125]
println("by end of 365 days: total vaccinated: $(sum(v1_int)[end])")
println("by first 80 days: total vaccinated: $(sum(v1_int)[80])")
s = sum([get_ode_class(i, sol) for i = 1:5])
v1 = sum(v1ag) # total of all age groups
v2 = sum(v2ag)
v1int = sum(v1_int)
totalvac = v1 + v2 # should equal v1int
xvals = sol.t
# plot one, plots the vaccine doses given to each age group
println("setting qt")
gpexec("set term x11 size 89,89 font 'Arial,10'")
@gp :agp "set grid" "set key right"
@gp :- :agp "set title 'Vaccination by Age Groups'"
@gp :- :agp "set xlabel 'Time (in days)'" "set ylabel 'Vaccinated Individuals'"
for i = 1:5
max_vax = params.flucov[i]
@gp :- :agp "set arrow from 0,$max_vax to 731,$max_vax nohead lw 1.4 dt 3 lc rgb '$(cols[i])'"
@gp :- :agp xvals v1_int[i] "with lines title 'ag: $i' lt 1 lw 2 dt 1 lc rgb '$(cols[i])' "
end
Gnuplot.save(:agp, term="pngcairo enhanced font 'Arial,10' size 700,400", output="vaccine_agegroups.png")
Gnuplot.save(:agp, term="pdfcairo enhanced font 'Arial,14' size 8.5,4.5", output="vaccine_agegroups.pdf")
# plot two, overall vaccine coverage
@gp :fsp "set grid" "set key right"
#@gp :- "set yrange [0:10000]" "set ytics 0,1000,10000 nomirror"
@gp :- :fsp "set title 'Vaccine Distribution (first, second dose)'"
@gp :- :fsp "set xlabel 'Time (in days)'" "set ylabel 'Vaccinated (proportion)"
@gp :- :fsp xvals s./totalpop "with lines title 'Susceptible' lw 3 dt 3 lc rgb '#666666' "
@gp :- :fsp xvals v1./totalpop "with lines title 'Vaccinated (first dose)' lw 3 dt 1 lc rgb 'black' "
@gp :- :fsp xvals v2./totalpop "with lines title 'Vaccinated (second dose)' lw 3 dt 1 lc rgb 'blue' "
#@gp :- :fsp "set obj rect from 0, 0 to 84, 10000 fc lt 2 fillstyle pattern 5 noborder"
#@gp :- :fsp "set arrow from 84,0 to 84,10000 nohead lw 1 lc rgb 'black' lt 3"
Gnuplot.save(:fsp, term="pngcairo enhanced font 'Arial,10' size 700,400", output="vaccine_coverage.png")
Gnuplot.save(:fsp, term="pdfcairo enhanced font 'Arial,14' size 8.5,4.5", output="vaccine_coverage.pdf")
# length of texas data is 80 days, december 14 + 80 days = march 3 (including dec 14)
tx_data = [1.520899434, 21.15809077, 85.87736239, 169.1681611, 299.5102773, 362.160298, 396.2787971, 492.8334943, 585.9428887, 721.3374259, 780.9801352, 789.6296041, 831.6250517, 862.881087, 1028.924679, 1256.097393, 1546.909919, 1715.553869, 1734.05987, 1821.054628, 1873.182508, 2094.395779, 2336.405021, 2635.949786, 2980.538695, 3342.023038, 3509.911712, 3575.19313, 3850.365568, 4162.080977, 4478.214236, 4786.763692, 5069.206097, 5261.011864, 5305.811146, 5439.512347, 5689.812388, 6033.997793, 6363.363912, 6708.252862, 6916.802318, 7000.279349, 7331.335357, 7765.69182, 8262.501724, 8743.547386, 9197.275486, 9424.720651, 9509.339219, 9882.804525, 10348.35495, 10895.26486, 11488.91571, 12121.77542, 12515.36074, 12617.85419, 12949.75169, 13420.04414, 13994.27162, 14462.58794, 15012.99145, 15319.22334, 15400.91047, 15440.0538, 15507.63209, 15569.63374, 15687.95006, 15894.42337, 16155.3904, 16334.91171, 16688.41909, 17080.86633, 17634.522, 18349.13091, 19225.34143, 19772.30997, 20003.84191, 20394.71651, 20683.99434, 20893.93365]
# ct data starts at january 12, which is 29 days after december 14.
ct_data = [1580, 1650, 4390, 4680, 4680, 4680, 4680, 5660, 5980, 6480, 6980, 7410, 7720, 7930, 8170, 8490, 8840, 9230, 9560, 9810, 10030, 10140, 10260, 10540, 10870, 11100, 11280, 11410, 11570, 11870, 12260, 12630, 12980, 13270, 13270, 13870, 14280, 14730, 15670, 16050, 16250, 16700, 16900, 17340, 17660, 18310, 18910, 19450, 20070, 20270, 20760]
# plot the first 80 days
@gp :mdp "set grid" "set key left"
@gp :- :mdp "set title 'Fitting to vaccination data'"
#@gp :- "set yrange [0:10000]" "set ytics 0,1000,10000 nomirror"
@gp :- :mdp "set xlabel 'Time (in days)'" "set ylabel 'Vaccinated Individuals"
@gp :- :mdp 1:80 v1int[1:80] "with lines title 'Model Output' lw 3 dt 1 lc rgb 'black' "
# for texas
#@gp :- :mdp xvals[1:length(tx_data)] tx_data.*289.96 "with boxes title 'MD data' fill pattern 1 lc rgb '#666666'"
# for CT
@gp :- :mdp 30:80 ct_data "with boxes title 'vaccine data' fill pattern 1 lc rgb '#666666'"
Gnuplot.save(:mdp, term="pdfcairo enhanced font 'Arial,14' size 8.5,4.5", output="vaccine_fitting.pdf")
Gnuplot.save(:mdp, term="pngcairo enhanced font 'Arial,10' size 700,400", output="vaccine_fitting.png")
return sol
end
export test_vaccine_allocation
function run_ct()
nsims = 500
betaval = 0.0210
#jldfile = jldopen("tx_sims.jld", "a+")
end
function run_tx()
nsims = 500
#CT: 0.0210, TX: 0.0243
betaval = 0.0210
jldfile = jldopen("ct_sims.jld", "a+")
println("First Day of Vaccination in CT: December 14")
lift_days = (90, 180)
trans_inc = (1.30, 1.60, 1.90)
#baseline status quo with vaccine on
p = Progress(nsims, barglyphs=BarGlyphs("[=> ]"))
bl = progress_map(1:nsims, progress = p, mapfun = pmap) do x
run_model(x, 1, betaval, 0.0, -1, -1, -1, 0.0, true)
end
bpd = create_pd(bl)
# scenarios
spd = Array{PlotData_t5, 2}(undef, length(lift_days), length(trans_inc))
for (i, ld) in enumerate(lift_days)
for (j, tr) in enumerate(trans_inc)
println("($i, $j) $ld $tr")
p = Progress(nsims, barglyphs=BarGlyphs("[=> ]"))
s = progress_map(1:nsims, progress = p, mapfun = pmap) do x
run_model(x, 1, betaval, tr, ld, 1, -1, 0.0, true)
end
spd[i, j] = create_pd(s)
end
end
println("saving to jld")
jldfile["bl"] = bpd
jldfile["spd"] = spd
close(jldfile)
return nothing
end
export run_tx
function plot_tx()
jldfile = jldopen("tx_sims.jld", "r")
bl = jldfile["bl"]
sp = jldfile["spd"]
close(jldfile)
bin_pd(bl, 28996000 / 100000)
for _sp in sp
bin_pd(_sp, 28996000 / 100000)
end
# size(sp) = 2, 3
# texas data from csv file
tx_data = reverse([7822, 7747, 3821, 3815, 11073, 7955, 7389, 7517, 11809, 6365, 4484, 6486, 2937, 3131, 3766, 3348, 3889, 6933, 11282, 12502, 11890, 12897, 13329, 7485, 6959, 13897, 14495, 15281, 17620, 23047, 31811, 11370, 19234, 19076, 18220, 19613, 26274, 6319, 11565, 17672, 22646, 22360, 31255, 9476, 11590, 16402, 24657, 27204, 23064, 27343, 26052, 14834, 15855, 23290, 23520, 24578, 24010, 31630, 18182, 16095, 4763, 16311, 18725, 21469, 32552, 14829, 8046, 2694, 4335, 17064, 23363, 21147, 10280, 7780, 17907, 16792, 19849, 18802, 18926, 8901])
tbindata = (tx_data ./ 289.96)
println("""fitting proces
sum total incidence tx: $(sum(tx_data))
sum total incidence model (time 1:80): $(sum(bl.inc[1][1:80]))
sum total incidence tx binned: $(sum(tbindata)))
""")
# define the date range of the model
drange = Date(2020, 12, 14):Day(1):Date(2021, 12, 14)
(length(drange) != length(bl.W[1])) && error("plotting date range not matching simulation time") # date range must match model output.
# _xtics = map(d -> Dates.format(d, "dd"), drange)
_myxtics_idx = findall(x -> Dates.day(x) == 1, drange)
_myxtics_months = Dates.monthabbr.(drange[_myxtics_idx])
myxtics = ""
for i = 1:length(_myxtics_idx)
myxtics = string(myxtics, "'$(_myxtics_months[i])' $(_myxtics_idx[i]), ")
end
println(myxtics) # copy paste into the gnuplot with minor string formatting (brackets, comma, etc)
_insetxtics = filter(x -> x <= Date(2021, 03, 03), drange)
insetxtics = ""
for i = 1:5:length(_insetxtics)
insetxtics = string(insetxtics, "'$(monthabbr(_insetxtics[i])) $(day(_insetxtics[i]))' $i, ")
end
insetxtics = string(insetxtics, "'Mar 3' 80") # add march 3d
println(insetxtics) # copy paste into the gnuplot with minor string formatting (brackets, comma, etc)
# figure for plotting the baseline and fitting process
@gp "reset session"
gpexec("set term qt size 650,450 font 'Arial,10'")
@gp "set grid" "set key right"
@gp :- "set multiplot" # rowsfirst margins 0.10,0.95,0.05,0.9 spacing 0.05,0.0" #l r b t :-
@gp :- "set xtics noenhanced nomirror" # no escaping for strings :-
@gp :- "set ytics nomirror" :-
@gp :- 1 bl.t bl.inc[2] bl.inc[3] "with filledcurves notitle fc '#4daf4a' fs solid 0.4 border" :-
@gp :- 1 bl.t bl.inc[1] "with lines notitle lw 3 lc rgb '#4daf4a'" :-
#@gp :- 1 1:80 tx_data "with line notitle lw 1 dashtype 3 lc rgb 'black'" :-
@gp :- "set xtics noenhanced nomirror" # no escaping for strings :-
@gp :- "set ytics nomirror" :-
@gp :- "set xlabel 'Time'" "set ylabel 'Incidence'" :-
@gp :- "set xrange [1:365]" :-
@gp :- "set xtics ('Jan' 19, 'Feb' 50, 'Mar' 78, 'Apr' 109, 'May' 139, 'Jun' 170, 'Jul' 200, 'Aug' 231, 'Sep' 262, 'Oct' 292, 'Nov' 323, 'Dec' 353)" :-
@gp :- "set xtics out" :-
@gp :- 1 "set obj 1 rect from 0, graph 0 to 80, graph 1 fillcolor lt 3 fillstyle pattern 5 noborder" :-
#@gp :- "set obj rect from 0, 0 to 80, 10000 fc lt 2 fillstyle pattern 5 noborder" :-
# add inset plot
@gp :- 2 1:80 tx_data "with line notitle lw 1 dashtype 1 lc rgb 'black'" :-
@gp :- 2 1:80 bl.inc[1][1:80] "with line notitle lw 1 dashtype 1 lc rgb 'green'" :-
#@gp :- "unset object 1" :-
@gp :- 2 "unset grid" "unset key" "unset xlabel" :-
@gp :- 2 "set xrange[1:80]" :-
@gp :- 2 "set size .50,.45" :-
@gp :- 2 "set xtics ('Dec 14' 1, 'Dec 19' 6, 'Dec 24' 11, 'Dec 29' 16, 'Jan 3' 21, 'Jan 8' 26, 'Jan 13' 31, 'Jan 18' 36, 'Jan 23' 41, 'Jan 28' 46, 'Feb 2' 51, 'Feb 7' 56, 'Feb 12' 61, 'Feb 17' 66, 'Feb 22' 71, 'Feb 27' 76, 'Mar 3' 80) rotate by 45 right in font ', 9'" :-
@gp :- 2 "set ytics font ', 9'" :-
@gp :- 2 "unset xlabel" "unset ylabel" :-
@gp :- 2 "set origin 0.47, 0.50"
Gnuplot.save(term="pdfcairo enhanced font 'Arial,14' size 6.5,3.5", output="sn.pdf")
Gnuplot.save(term="pngcairo enhanced font 'Arial,10' size 650,350", output="sn.png")
Gnuplot.save(term="svg enhanced font 'Arial,10' size 650,350", output="sn.svg")
# The svg terminal doesn’t seem to accept the ‘in’ argument to the size option, but in SVG-land 100=1in.
# main figures
# first have to filter the date from march 3 onwards, and get their index (this matches the sol.t index)
fidx = findall(x -> x > Date(2021, 03, 03), drange)
fidx_month = findall(x -> x >= Date(2021, 03, 03) && day(x) == 1, drange)
# the fidx_month just gives you the time index of the first of the months
# println(fidx_month) # to generate the xtics
# Incidence Curves
@gp "reset session";
@gp "set grid"
@gp :- "set multiplot layout 2,2" # rowsfirst margins 0.10,0.95,0.05,0.9 spacing 0.05,0.0" #l r b t
@gp :- "set format y '%.0f'"
@gp :- "set xrange [80:365]"
@gp :- "set yrange [0:*]"
@gp :- "set format y '%.f'"
@gp :- "set xtics ('Mar' 80, 'Apr' 109, 'May' 139, 'Jun' 170, 'Jul' 200, 'Aug' 231, 'Sep' 262, 'Oct' 292, 'Nov' 323, 'Dec' 353)"
@gp :- 1 bl.t[fidx] bl.inc[2][fidx] bl.inc[3][fidx] "with filledcurves notitle fc '#66c2a5' fs solid 0.5 border"
@gp :- 1 bl.t[fidx] bl.inc[1][fidx] "with lines notitle lw 3 lc rgb '#66c2a5'"
@gp :- "set title 'Base Line'"
@gp :- 2 sp[1, 1].t[fidx] sp[1, 1].inc[2][fidx] sp[1, 1].inc[3][fidx] "with filledcurves notitle fc '#e41a1c' fs solid 0.5 border"
@gp :- 2 sp[1, 1].t[fidx] sp[1, 1].inc[1][fidx] "with lines title 'Mar 10' lw 3 lc rgb '#e41a1c'"
@gp :- 2 sp[2, 1].t[fidx] sp[2, 1].inc[2][fidx] sp[2, 1].inc[3][fidx] "with filledcurves notitle fc '#8da0cb' fs solid 0.5 border"
@gp :- 2 sp[2, 1].t[fidx] sp[2, 1].inc[1][fidx] "with lines title 'Jun 12' lw 3 lc rgb '#8da0cb'"
@gp :- "set title '30% Increase'"
@gp :- 3 sp[1, 2].t[fidx] sp[1, 2].inc[2][fidx] sp[1, 2].inc[3][fidx] "with filledcurves notitle fc '#e41a1c' fs solid 0.5 border"
@gp :- 3 sp[1, 2].t[fidx] sp[1, 2].inc[1][fidx] "with lines title 'Mar 10' lw 3 lc rgb '#e41a1c'"
@gp :- 3 sp[2, 2].t[fidx] sp[2, 2].inc[2][fidx] sp[2, 2].inc[3][fidx] "with filledcurves notitle fc '#8da0cb' fs solid 0.5 border"
@gp :- 3 sp[2, 2].t[fidx] sp[2, 2].inc[1][fidx] "with lines title 'Jun 12' lw 3 lc rgb '#8da0cb'"
@gp :- "set title '60% Increase'"
@gp :- 4 sp[1, 3].t[fidx] sp[1, 3].inc[2][fidx] sp[1, 3].inc[3][fidx] "with filledcurves notitle fc '#e41a1c' fs solid 0.5 border"
@gp :- 4 sp[1, 3].t[fidx] sp[1, 3].inc[1][fidx] "with lines title 'Mar 10' lw 3 lc rgb '#e41a1c'"
@gp :- 4 sp[2, 3].t[fidx] sp[2, 3].inc[2][fidx] sp[2, 3].inc[3][fidx] "with filledcurves notitle fc '#8da0cb' fs solid 0.5 border"
@gp :- 4 sp[2, 3].t[fidx] sp[2, 3].inc[1][fidx] "with lines title 'Jun 12' lw 3 lc rgb '#8da0cb'"
@gp :- "set title '90% Increase'"
Gnuplot.save(term="pdfcairo enhanced font 'Arial,14' size 8.5,8.5", output="tx_fig1_incidence.pdf")
Gnuplot.save(term="pngcairo enhanced font 'Arial,10' size 850,850", output="tx_fig1_incidence.png")
Gnuplot.save(term="svg enhanced font 'Arial,10' size 850,850", output="tx_fig1_incidence.svg")
#cumulative with bars
@gp "reset" :-
@gp :- "set multiplot layout 1,3" # rowsfirst margins 0.10,0.95,0.05,0.9 spacing 0.05,0.0" #l r b t
@gp :- "set style data histograms" :-
@gp :- "set key left"
@gp :- "set style histogram cluster errorbars gap 1" :-
# color them
#@gp :- "set style histogram errorbars linewidth 3" :-
@gp :- "set errorbars linecolor black linewidth 2" :-
@gp :- "set bars front" :-
@gp :- "set boxwidth 0.9 abs" :-
@gp :- "set xrange [-0.5:2.5]" :-
@gp :- "set xtics nomirror"
@gp :- "set style fill solid border -1" :-
@gp :- "set format y '%.0f'" :-
@gp :- "set xtics ('30%%' 0, '60%%' 1, '90%%' 2)"
@gp :- "set xtics format '' nomirror"
@gp :- "set format x '%g%%'"
#@gp :- "set yrange[2000000:20000000]"
#@gp :- "set ytics ('2,000,000' 2000000, '4,000,000' 4000000, '6,000,000' 6000000, '8,000,000' 8000000, '10,000,000' 10000000, '12,000,000' 12000000, '14,000,000' 14000000, '16,000,000' 16000000, '18,000,000' 18000000,)"
@gp :- "set xtics out "
@gp :- "unset xlabel"
@gp :- "set ylabel 'Total Infections' offset 0,0,0"
#@gp :- "unset key"
# @gp :- "set key at screen 0.50,screen 0.96 #for example"
# @gp :- "set key box lt -1 lw 1"
# #@gp :- "set key spacing 1 font 'Helvetica, 14'"
sub_fidx(arr, idx) = arr[end] - arr[idx[1] - 1]
# to get filtered cumulative data, can substract W[end] - W[80] or sum up icndeice
mean_incs = [sub_fidx(sp[1, i].W[1], fidx) for i = 1:3]
lo_incs = [sub_fidx(sp[1, i].W[2], fidx) for i = 1:3]
hi_incs = [sub_fidx(sp[1, i].W[3], fidx) for i = 1:3]
@gp :- 1 mean_incs lo_incs hi_incs "notitle lc rgb '#e41a1c'"
mean_incs = [sub_fidx(sp[2, i].W[1], fidx) for i = 1:3]
lo_incs = [sub_fidx(sp[2, i].W[2], fidx) for i = 1:3]
hi_incs = [sub_fidx(sp[2, i].W[3], fidx) for i = 1:3]
@gp :- 1 mean_incs lo_incs hi_incs "notitle lc rgb '#8da0cb'"
@gp :- "unset xlabel"
mean_incs = [sub_fidx(sp[1, i].U[1], fidx) for i = 1:3]
lo_incs = [sub_fidx(sp[1, i].U[2], fidx) for i = 1:3]
hi_incs = [sub_fidx(sp[1, i].U[3], fidx) for i = 1:3]
@gp :- 2 mean_incs lo_incs hi_incs "notitle lc rgb '#e41a1c'"
mean_incs = [sub_fidx(sp[2, i].U[1], fidx) for i = 1:3]
lo_incs = [sub_fidx(sp[2, i].U[2], fidx) for i = 1:3]
hi_incs = [sub_fidx(sp[2, i].U[3], fidx) for i = 1:3]
@gp :- 2 mean_incs lo_incs hi_incs "notitle lc rgb '#8da0cb'"
@gp :- "set ylabel 'Hospitalizations'"
@gp :- "set xlabel 'Transmission Increase'"
mean_incs = [sub_fidx(sp[1, i].D[1], fidx) for i = 1:3]
lo_incs = [sub_fidx(sp[1, i].D[2], fidx) for i = 1:3]
hi_incs = [sub_fidx(sp[1, i].D[3], fidx) for i = 1:3]
@gp :- 3 mean_incs lo_incs hi_incs "title 'Mar 3' lc rgb '#e41a1c'"
mean_incs = [sub_fidx(sp[2, i].D[1], fidx) for i = 1:3]
lo_incs = [sub_fidx(sp[2, i].D[2], fidx) for i = 1:3]
hi_incs = [sub_fidx(sp[2, i].D[3], fidx) for i = 1:3]
@gp :- 3 mean_incs lo_incs hi_incs "title 'Jun 12' lc rgb '#8da0cb'"
@gp :- "set ylabel 'Deaths'"
@gp :- "unset xlabel"
Gnuplot.save(term="pdfcairo enhanced font 'Arial,14' size 8.5,2.5", output="tx_fig2_bars.pdf")
Gnuplot.save(term="pngcairo enhanced font 'Arial,10' size 850,250", output="tx_fig2_bars.png")
Gnuplot.save(term="svg enhanced font 'Arial,10' size 850,250", output="tx_fig2_bars.svg")
# lets print all the information/data for all scenarios
println("""
Baseline
Inf: $(sub_fidx(bl.W[1], fidx)) CrI: $(sub_fidx(bl.W[2], fidx)) - $(sub_fidx(bl.W[3], fidx))
Hos: $(sub_fidx(bl.U[1], fidx)) CrI: $(sub_fidx(bl.U[2], fidx)) - $(sub_fidx(bl.U[3], fidx))
Ded: $(sub_fidx(bl.D[1], fidx)) CrI: $(sub_fidx(bl.D[2], fidx)) - $(sub_fidx(bl.D[3], fidx))
""")
lift_days = (90, 180)
trans_inc = (1.30, 1.60, 1.90)
for (j, ti) in enumerate(trans_inc)
for (i, ld) in enumerate(lift_days)
_pdojb = sp[i, j]
println("""
Index: $i ($ld), $j ($ti)
Inf: $(sub_fidx(_pdojb.W[1], fidx)) CrI: $(sub_fidx(_pdojb.W[2], fidx)) - $(sub_fidx(_pdojb.W[3], fidx))
Hos: $(sub_fidx(_pdojb.U[1], fidx)) CrI: $(sub_fidx(_pdojb.U[2], fidx)) - $(sub_fidx(_pdojb.U[3], fidx))
Ded: $(sub_fidx(_pdojb.D[1], fidx)) CrI: $(sub_fidx(_pdojb.D[2], fidx)) - $(sub_fidx(_pdojb.D[3], fidx))
""")
end
end
# Cumulative, Hospitalizations, Deaths
@gp "reset session";
@gp "set grid"
@gp :- "set multiplot layout 3,2 margins 0.1,0.99,0.02,0.90 spacing 0.1,0.03" # rowsfirst margins 0.10,0.95,0.05,0.9 spacing 0.05,0.0" #l r b t
@gp :- "set format y '%.0f'"
#@gp :- "unset key"
@gp :- "set xrange [80:365]"
@gp :- 1 "set key left"
@gp :- 1 sp[1, 1].t[fidx] sp[1, 1].W[2][fidx] sp[1, 1].W[3][fidx] "with filledcurves notitle fc '#fc8d62' fs solid 0.5 border"
@gp :- 1 sp[1, 1].t[fidx] sp[1, 1].W[1][fidx] "with lines title '30% Increase' lw 3 lc rgb '#fc8d62'"
@gp :- 1 sp[1, 2].t[fidx] sp[1, 2].W[2][fidx] sp[1, 2].W[3][fidx] "with filledcurves notitle fc '#8da0cb' fs solid 0.5 border"
@gp :- 1 sp[1, 2].t[fidx] sp[1, 2].W[1][fidx] "with lines title '60% Increase' lw 3 lc rgb '#8da0cb'"
@gp :- 1 sp[1, 3].t[fidx] sp[1, 3].W[2][fidx] sp[1, 3].W[3][fidx] "with filledcurves notitle fc '#e78ac3' fs solid 0.5 border"
@gp :- 1 sp[1, 3].t[fidx] sp[1, 3].W[1][fidx] "with lines title '90% Increase' lw 3 lc rgb '#e78ac3'"
@gp :- "set ytics 5000000"
@gp :- "set ylabel 'Cumulative Incidence'"
@gp :- 2 "set key left"
@gp :- 2 bl.t bl.U[2] bl.U[3] "with filledcurves notitle fc '#66c2a5' fs solid 0.5 border"
@gp :- 2 bl.t bl.U[1] "with lines title 'Base Line' lw 3 lc rgb '#66c2a5'"
@gp :- 2 s1.t s1.U[2] s1.U[3] "with filledcurves notitle fc '#fc8d62' fs solid 0.5 border"
@gp :- 2 s1.t s1.U[1] "with lines title '30% Increase' lw 3 lc rgb '#fc8d62'"
@gp :- 2 s2.t s2.U[2] s2.U[3] "with filledcurves notitle fc '#8da0cb' fs solid 0.5 border"
@gp :- 2 s2.t s2.U[1] "with lines title '60% Increase' lw 3 lc rgb '#8da0cb'"
@gp :- 2 s3.t s3.U[2] s3.U[3] "with filledcurves notitle fc '#e78ac3' fs solid 0.5 border"
@gp :- 2 s3.t s3.U[1] "with lines title '90% Increase' lw 3 lc rgb '#e78ac3'"
@gp :- "set ylabel 'Hospitalizations'"
@gp :- 3 "set key left"
@gp :- 3 bl.t bl.D[2] bl.D[3] "with filledcurves notitle fc '#66c2a5' fs solid 0.5 border"
@gp :- 3 bl.t bl.D[1] "with lines title 'Base Line' lw 3 lc rgb '#66c2a5'"
@gp :- 3 s1.t s1.D[2] s1.D[3] "with filledcurves notitle fc '#fc8d62' fs solid 0.5 border"
@gp :- 3 s1.t s1.D[1] "with lines title '30% Increase' lw 3 lc rgb '#fc8d62'"
@gp :- 3 s2.t s2.D[2] s2.D[3] "with filledcurves notitle fc '#8da0cb' fs solid 0.5 border"
@gp :- 3 s2.t s2.D[1] "with lines title '60% Increase' lw 3 lc rgb '#8da0cb'"
@gp :- 3 s3.t s3.D[2] s3.D[3] "with filledcurves notitle fc '#e78ac3' fs solid 0.5 border"
@gp :- 3 s3.t s3.D[1] "with lines title '90% Increase' lw 3 lc rgb '#e78ac3'"
@gp :- "set ylabel 'Deaths'"
Gnuplot.save(term="pdfcairo enhanced font 'Arial,14' size 8.5,11", output="tx_fig1.pdf")
Gnuplot.save(term="pngcairo enhanced font 'Arial,10' size 900,1200", output="tx_fig1.png")
end
export plot_tx
function get_ode_class(x, sol)
p = [sol.u[i][x] for i = 1:length(sol.u)]
return p
end
function get_compartments(sol, tidx=0)
# the tidx can subset the array
# note the step size so time = 80 correponds to index of 160 due to half step size
last_idx = tidx == 0 ? length(sol.t) : tidx
S = round.(sum([get_ode_class(i, sol) for i = 1:5]); digits = 3)[1:last_idx]
V₁ = round.(sum([get_ode_class(i, sol) for i = 6:10]); digits = 3)[1:last_idx]
V₂ = round.(sum([get_ode_class(i, sol) for i = 11:15]); digits = 3)[1:last_idx]
E = round.(sum([get_ode_class(i, sol) for i = 16:20]); digits = 3)[1:last_idx]
F₁ = round.(sum([get_ode_class(i, sol) for i = 21:25]); digits = 3)[1:last_idx]
F₂ = round.(sum([get_ode_class(i, sol) for i = 26:30]); digits = 3)[1:last_idx]
Ẽ = round.(sum([get_ode_class(i, sol) for i = 31:35]); digits = 3)[1:last_idx]
F̃₁ = round.(sum([get_ode_class(i, sol) for i = 36:40]); digits = 3)[1:last_idx]
F̃₂ = round.(sum([get_ode_class(i, sol) for i = 41:45]); digits = 3)[1:last_idx]
A = round.(sum([get_ode_class(i, sol) for i = 46:50]); digits = 3)[1:last_idx]
H̃ = round.(sum([get_ode_class(i, sol) for i = 51:55]); digits = 3)[1:last_idx]
à = round.(sum([get_ode_class(i, sol) for i = 56:60]); digits = 3)[1:last_idx]
P = round.(sum([get_ode_class(i, sol) for i = 61:65]); digits = 3)[1:last_idx]
B̃ = round.(sum([get_ode_class(i, sol) for i = 66:70]); digits = 3)[1:last_idx]
P̃ = round.(sum([get_ode_class(i, sol) for i = 71:75]); digits = 3)[1:last_idx]
I = round.(sum([get_ode_class(i, sol) for i = 76:80]); digits = 3)[1:last_idx]
Q̃ = round.(sum([get_ode_class(i, sol) for i = 81:85]); digits = 3)[1:last_idx]
Ĩ = round.(sum([get_ode_class(i, sol) for i = 86:90]); digits = 3)[1:last_idx]
X = round.(sum([get_ode_class(i, sol) for i = 91:95]); digits = 3)[1:last_idx]
D = round.(sum([get_ode_class(i, sol) for i = 96:100]); digits = 3)[1:last_idx]
R = round.(sum([get_ode_class(i, sol) for i = 101:105]); digits = 3)[1:last_idx]
W = round.(sum([get_ode_class(i, sol) for i = 106:110]); digits = 3)[1:last_idx]
Z = round.(sum([get_ode_class(i, sol) for i = 111:115]); digits = 3)[1:last_idx]
Y = round.(sum([get_ode_class(i, sol) for i = 116:120]); digits = 3)[1:last_idx]
V = round.(sum([get_ode_class(i, sol) for i = 121:125]); digits = 3)[1:last_idx]
U = round.(sum([get_ode_class(i, sol) for i = 126:130]); digits = 3)[1:last_idx]
inc_w = W - circshift(W, 1)
#inc_w[1] = 0 # no incidence on day 1.
inc_w[1] = inc_w[2] # setting it to zero causes crazy jump
inc_a = Z - circshift(Z, 1)
#inc_a[1] = 0
inc_a[1] = inc_a[2]
return (S = S, V₁ = V₁, V₂ = V₂, E = E, F₁ = F₁, F₂ = F₂, Ẽ = Ẽ, F̃₁ = F̃₁, F̃₂ = F̃₂, A = A, H̃ = H̃, Ã = Ã,
P = P, B̃ = B̃, P̃ = P̃, I = I, Q̃ = Q̃, Ĩ = Ĩ, X = X, D = D, R = R, W = W, Z = Z, Y = Y, inc = inc_w, inc_a = inc_a, V=V, U=U)
end
export get_compartments
# sol.t[end] = 365 (out of 731 total points due to half step)
# but we take care of the half step by multiplying by 2
# so it's more easier to interpret plot_model(sol, 50) ie 50 days of data
plot_model(sol) = plot_model(sol, length(sol.t))
function plot_model(sol, tidx)
xvals = sol.t[1:tidx]
yvals = get_compartments(sol, tidx)
tx_data = reverse([7822, 7747, 3821, 3815, 11073, 7955, 7389, 7517, 11809, 6365, 4484, 6486, 2937, 3131, 3766, 3348, 3889, 6933, 11282, 12502, 11890, 12897, 13329, 7485, 6959, 13897, 14495, 15281, 17620, 23047, 31811, 11370, 19234, 19076, 18220, 19613, 26274, 6319, 11565, 17672, 22646, 22360, 31255, 9476, 11590, 16402, 24657, 27204, 23064, 27343, 26052, 14834, 15855, 23290, 23520, 24578, 24010, 31630, 18182, 16095, 4763, 16311, 18725, 21469, 32552, 14829, 8046, 2694, 4335, 17064, 23363, 21147, 10280, 7780, 17907, 16792, 19849, 18802, 18926, 8901])
tbindata = (tx_data ./ 289.96)
ct_data = reverse([494, 502, 2680, 0, 0, 787, 975, 1493, 1357, 2233, 0, 0, 1198, 547, 534, 580, 2905, 0, 0, 838, 1003, 888, 869, 4367, 0, 0, 1431, 937, 482, 2568, 3931, 0, 0, 1258, 1426, 2440, 1267, 5817, 0, 0, 2019, 1662, 1915, 2094, 6703, 0, 0, 1878, 968, 3529, 3689, 7364, 0, 0, 3236, 3304, 2486, 2332, 4516, 0, 4412, 0, 2045, 1696, 767, 8457, 0, 0, 0, 2038, 1745, 1583, 4595, 0, 0, 2680, 2321, 2319, 1470, 7231])
tbindata = sevenaverage(ct_data ./ 35.65)
# plot non_infected classes (S, V, R, D)
@gp "reset session"
@gp "set grid" "set key right"
@gp :- "unset colorbox"
@gp :- "set multiplot layout 3,1" # rowsfirst margins 0.10,0.95,0.05,0.9 spacing 0.05,0.0" #l r b t
@gp :- 1 xvals yvals.S "with lines title 'susceptible' lw 2 dt 1 lc rgb 'black' " :-
@gp :- 1 xvals yvals.V₁ "with lines title 'vaccinated 1' lw 2 dt 7 lc rgb 'black'" :-
@gp :- 1 xvals yvals.V₂ "with lines title 'vaccinated 2' lw 2 dt 7 lc rgb 'green'" :-
@gp :- "set label 'overall AR: $((yvals.W)[end])' at graph 0.1, 0.5" :-
@gp :- "set label 'sum incidence: $(sum(yvals.inc[1:80]))' at graph 0.1, 0.6"
@gp :- "set label 'sum tx: $(sum(tbindata))' at graph 0.1, 0.8"
@gp :- 2 xvals yvals.inc "with lines title 'incidence' lw 1 dt 1 lc rgb 'black'" :-
@gp :- 2 tbindata "with lines title 'tx d14 - mar3"
@gp :- "unset label"
#yrangemax = maximum(yvals.inc) + 1
#@gp :- "set yrange [0:$yrangemax]"
@gp :- "set ytics auto" #tx incidence
# silent infections
@gp :- 3 xvals (yvals.W) "with lines title 'cumulative incidence' lw 2 dt 1 lc rgb 'black' " :-
@gp :- 3 xvals (yvals.Z) "with lines title 'asymp' lw 2 dt 1 lc rgb 'green'" :-
@gp :- 3 xvals (yvals.Y) "with lines title 'symp' lw 2 dt 1 lc rgb 'red'" :-
@gp :- "set yrange [*:*]"
@gp :- "set ytics auto"
# @gp "set grid" "set key right"
# @gp :- "unset colorbox"
# @gp :- "set yrange [0:10000]"
# @gp :- "set y2range [0:1600]"
# @gp :- "set y2tics 0, 100"
# @gp :- "set xtics 0,2,150"
# @gp :- "set ytics nomirror"
# @gp :- "set xlabel 'Time'" "set ylabel 'Prevalence'"
@gp
#save(term="pdfcairo enhanced size 8.5in,8.5in fontscale 0.8", output="basic_results.pdf")
end
export plot_model
ngm2(beta) = ngm2(beta, false)
function ngm2(beta, st::String, sample_p::Bool)
# model parameters
params = ModelParameters()
setup_state_params(params, st)
sample_p && sample_dis_params(params)
params.β = beta
params.s_vac = 0
params.v_vac = 0
default_qg_params(params)
ngm2(beta, params)
end
function ngm2(beta, params::ModelParameters)
# ngm https://mtbi.asu.edu/sites/ult/files/brn_mtbi_2016.pdf (see downloaded pdf)
## ix: DFE for the infect subsystem (passes to subfunction _f)
@unpack β, α, ϵ₁, ϵ₂, λ, s_vac, v_vac, p, ρ₁, ρ₂, σ, η, θ, γ, δ, τ, κ, νᵣ, νₓ, q, g, f, h, d, pop, npi_start = params
# calculate the number of susceptibles
total_pop = params.pop
preexist_imm = params.preexist_imm
susp = (1 - preexist_imm) .* (total_pop)
S₁ = susp[1]; S₂ = susp[2]; S₃ = susp[3]; S₄ = susp[4]; S₅ = susp[5];
R₁ = 0; R₂ = 0; R₃ = 0; R₄ = 0; R₅ = 0;
M, M̃ = contact_matrix()
# the disease free equilibrium
dfe_x = [0 for _ = 1:55]
function _f(x)
# x = DFE for each compartment
E = [x[1], x[2], x[3], x[4], x[5]]
Ẽ = [x[6], x[7], x[8], x[9], x[10]]
A = [x[11], x[12], x[13], x[14], x[15]]
H̃ = [x[16], x[17], x[18], x[19], x[20]]
à = [x[21], x[22], x[23], x[24], x[25]]
P = [x[26], x[27], x[28], x[29], x[30]]
B̃ = [x[31], x[32], x[33], x[34], x[35]]
P̃ = [x[36], x[37], x[38], x[39], x[40]]
I = [x[41], x[42], x[43], x[44], x[45]]
Q̃ = [x[46], x[47], x[48], x[49], x[50]]
Ĩ = [x[51], x[52], x[53], x[54], x[55]]
foi = [β * (α*dot(M[a, :], A./pop) + dot(M[a, :], P./pop) + dot(M[a, :], I./pop) + α*dot(M̃[a, :], Ã./pop) + α*dot(M̃[a, :], H̃./pop) +
dot(M̃[a, :], P̃./pop) + dot(M̃[a, :], B̃./pop) + dot(M̃[a, :], Ĩ./pop) + dot(M̃[a, :], Q̃./pop)) for a = 1:5]
[(1 - q[1])*S₁*foi[1], (1 - q[2])*S₂*foi[2], (1 - q[3])*S₃*foi[3], (1 - q[4])*S₄*foi[4], (1 - q[5])*S₅*foi[5],
q[1]*S₁*foi[1], q[2]*S₂*foi[2], q[3]*S₃*foi[3], q[4]*S₄*foi[4], q[5]*S₅*foi[5],
[0 for _=1:45]...]
end
function _y(x)
E = [x[1], x[2], x[3], x[4], x[5]]
Ẽ = [x[6], x[7], x[8], x[9], x[10]]
A = [x[11], x[12], x[13], x[14], x[15]]
H̃ = [x[16], x[17], x[18], x[19], x[20]]
à = [x[21], x[22], x[23], x[24], x[25]]
P = [x[26], x[27], x[28], x[29], x[30]]
B̃ = [x[31], x[32], x[33], x[34], x[35]]
P̃ = [x[36], x[37], x[38], x[39], x[40]]
I = [x[41], x[42], x[43], x[44], x[45]]
Q̃ = [x[46], x[47], x[48], x[49], x[50]]
Ĩ = [x[51], x[52], x[53], x[54], x[55]]
# strs = ["E", "Et", "A", "H", "At", "P", "B", "Pt", "I", "Q", "It"]
# retvals = map(strs) do x
# return ["$x$(i)" for i = 1:5]
# end
# println(vcat(retvals...))
Eₓ = σ .* E
Ẽₓ = σ .* Ẽ
Aₓ = -(p .* σ .* E) + ((1 .- g) .* η .* A) + (g .* δ .* A)
H̃ₓ = -(p .* σ .* Ẽ) + (η .* H̃)
Ãₓ = -(g .* δ .* A) + ((δ*η / (δ - η)) .* Ã)
Pₓ = -((1 .- p) .* σ .* E) + ((1 .- g) .* θ .* P) + (g .* δ .* P)
B̃ₓ = -((1 .- p) .* σ .* Ẽ) + (θ .* B̃)
P̃ₓ = -(g .* δ .* P) + ((δ*θ / (δ - θ)) .* P̃)
Iₓ = -((1 .- g) .* θ .* P) + ((1 .- f) .* γ .* I) + (f .* τ .* I)
Q̃ₓ = -(θ .* B̃) - ((δ*θ / (δ - θ)) .* P̃) + (γ .* Q̃)
Ĩₓ = -(f .* τ .* I) + ((γ*τ / (τ - γ)) .* Ĩ)
[Eₓ..., Ẽₓ..., Aₓ..., H̃ₓ..., Ãₓ..., Pₓ..., B̃ₓ..., P̃ₓ..., Iₓ..., Q̃ₓ..., Ĩₓ...]
end
gf = x -> ForwardDiff.jacobian(_f, x)
gy = x -> ForwardDiff.jacobian(_y, x)
fv = gf(dfe_x) * inv(gy(dfe_x))
#gv = gy(dfe_x)
#spec_radius = eigen(fv * gv)
#return maximum(spec_radius.values)
#sfv * gv
#gy(dfe_x)
#inv(gy(dfe_x))
maximum(abs.(eigen(fv).values))
end
export ngm2
function plot_ngm2(beta)
## samples ngm2 for sampled parameters.
## we want to adjust beta to get to the R we want
rvals = [ngm2(0.0240, true) for _ = 1:1000]
println("mean: $(mean(rvals))")
@gp "reset" :-
@gp :- "set style fill solid 0.5 border -1 " :-
@gp :- "set style data boxplot" :-
@gp :- "set boxwidth 0.5" :-
@gp :- "set pointsize 0.5" :-
#@gp :- "set xtics ('Using Next Generation Matrix (1000 replicates)' 1.00000, 'B' 2.00000)"
#@gp :- "set xtics nomirror"
@gp :- "unset xtics" :-
@gp :- "set ylabel 'Effective Reproduction Number'" :-
@gp :- "unset key" :-
@gp :- "set ytics nomirror" :-
@gp :- [1 for _i=1:length(rvals)] rvals "lc rgb '#fc8d62'"
Gnuplot.save(term="pdfcairo enhanced font 'Arial,14' size 4.5,4.5", output="r_boxplot.pdf")
Gnuplot.save(term="pngcairo enhanced font 'Arial,10' size 400,400", output="r_boxplot.png")
end
## comment when uploading to server
# using RCall
# @rlibrary EpiEstim
# #@rlibrary R0
# function r_estim(inc)
# # maybe @rlibrary in global scope.
# # serial interval from https://bmcinfectdis.biomedcentral.com/articles/10.1186/s12879-020-05577-4
# re_zero = filter(x -> x != 0 && x > 0, inc)
# estR0 = R"""
# res = EpiEstim::estimate_R($inc, method = "parametric_si", config = EpiEstim::make_config(list(mean_si = 5.30, std_si = 0.26)))
# """
# rval_jl = rcopy(estR0)
# return rval_jl
# end
end
# ####################################################################
# # DEPRECATED FUNCTSION.. includes good commands for heatmap plotting
# # deprecated function for 3d plotting only
# function plot_time_to_coverage()
# dda = readdlm("vaccinecoverages.dat")
# dda = convert(Matrix, dda)
# # lets plot the data as well
# @gp "reset session"; gpexec("set term qt size 700,400 font 'Arial,10'"); @gp "set grid"
# @gp :- "set size ratio 0.5"
# @gp :- "set lmargin at screen 0.0"
# @gp :- "set rmargin at screen 0.9"
# #@gp :- "set bmargin at screen 0.1"
# #@gp :- "set tmargin at screen 0.95"
# # @gp :- "set view map" # bird eye view of 3d plot or #set pm3d map
# # if i don't put "with pm3d" here, I have to have "set pm3d map interpolate" ..
# # don't know why the 'set view map' dosn't work
# @gp :- "set style increment user"
# #@gp :- "set dgrid3d 100,100,100"
# @gp :- "set grid x y lc 'black' front"
# @gp :- "unset border"
# # Note that grid lines only appear at axis tic locations.
# # But black tics would hide the white grid lines, so scale them to 0 length
# @gp :- "set tics scale 0"
# @gp :- "set pm3d map interpolate 1,1" # where "x" and "y" are the number of additional points in the x and y axes. So if "x" is set to 4 then Gnuplot will interpolate with 4 times as many points in the x-axis.
# #@gp :- "set format cb '%g%%'"
# @gp :- "set palette rgbformulae 33,13,10"
# @gp :- "set xtics 20,2,40"
# @gp :- "set xtics nomirror"
# @gp :- "set xlabel 'Dosage'"
# @gp :- "set ytics offset 1.5,0,0"
# #@gp :- "set label '' at screen 0.05,0.5 center front rotate"
# @gp :- "set ylabel 'Coverages' offset 1,0,0"
# @gp :- "set cblabel 'Time to coverages' offset 1,0,0"
# @gp :- "set style textbox opaque noborder"
# @gsp :- 20:40 5:35 dda "notitle ls 4" # if i don't put "with pm3d" here, I have to have "set pm3d map interpolate" .. don't know why the 'set view map' dosn't work
# # don't really need the ls 4, can turn ls 0 since we arn't drawing the surface lines since pm3d is set to map mode.
# #@gsp :- "set grid ytics"
# save(term="pdfcairo enhanced font 'Arial,14' size 8.5,4.5", output="time_to_coverage.pdf")
# save(term="pngcairo enhanced font 'Arial,10' size 700,400", output="time_to_coverage.png")
# @gp
# end
# # 3d print commands
# function plot_figure3()
# arfile_fast = readdlm("vaccinecoverages_attackrates_60.dat") .* 100
# arfile_med = readdlm("vaccinecoverages_attackrates_90.dat") .* 100
# arfile_slow = readdlm("vaccinecoverages_attackrates_120.dat") .* 100
# filetouse = arfile_fast
# #rows = coverage 5:60 (56 elements), cols = dosage level 20:80 (61 elements)
# #baseline = 0.28
# #aadat = reduce(hcat, aa)'
# #aadat = pinc.(aadat, baseline)
# aa = findfirst(x -> x < 38.59, filetouse)
# # add +4 to the row to get the coverage value
# # rows are vaccine coverages. columns are dosage values.
# # so row 22 col 2 is the first value when attack rate is less than 4.9
# # this corresponds to vaccination coverage of 22%
# @gp "reset session"
# gpexec("set term qt size 700,400 font 'Arial,10'")
# @gp "set grid"
# #@gp :- "set multiplot layout 3,1 margins 0.05,0.98,0.05,0.98 spacing 0,0.030" # rowsfirst margins 0.10,0.95,0.05,0.9 spacing 0.05,0.0" #l r b t
# @gp :- "set size ratio 0.5"
# @gp :- "set lmargin at screen 0.0"
# @gp :- "set rmargin at screen 0.9"
# #@gp :- "set bmargin at screen 0.1"
# #@gp :- "set tmargin at screen 0.95"
# @gp :- "set style increment user"
# @gp :- "set style line 1 lw 3" # the width of the contour line uses line style 1 by ult
# @gp :- "set style line 2 lc rgb 'black' " # contour lines uses line style 2 by ult.
# @gp :- "set contour"
# @gp :- "set cntrparam levels discrete 38.59" # Plot the selected contours
# @gp :- "set view map" # bird eye view of 3d plot
# @gp :- "set dgrid3d 10,10,10"
# @gp :- "set grid x y lc 'black' front"
# @gp :- "unset border"
# # Note that grid lines only appear at axis tic locations.
# # But black tics would hide the white grid lines, so scale them to 0 length
# @gp :- "set tics scale 0"
# @gp :- "set pm3d interpolate 0,0" # where "x" and "y" are the number of additional points in the x and y axes. So if "x" is set to 4 then Gnuplot will interpolate with 4 times as many points in the x-axis.
# @gp :- "set format cb '%g%%'"
# @gp :- "set palette rgbformulae 33,13,10"
# @gp :- "set xtics 5,5,60"
# @gp :- "set xtics nomirror"
# @gp :- "set xlabel 'Population Coverage when NPI is lifted'"
# @gp :- "set ytics offset 1.5,0,0"
# #@gp :- "set label '' at screen 0.05,0.5 center front rotate"
# @gp :- "set ylabel 'Vaccine Dose Rate per 10,000 individuals' offset 1,0,0"
# @gp :- "set cblabel 'Attack Rate' offset 1,0,0"
# @gp :- "set style textbox opaque noborder"
# @gsp :- 5:35 20:40 filetouse "w pm3d notitle ls 1" :-
# #@gsp :- 2 5:35 20:40 arfile_med "w pm3d notitle ls 1" :-
# #@gsp :- 3 5:35 20:40 arfile_fast "w pm3d notitle ls 1" :-
# #@gsp :- "set grid ytics"
# save(term="pdfcairo enhanced font 'Arial,14' size 8.5,4.5", output="Figure3_high.pdf")
# save(term="pngcairo enhanced font 'Arial,10' size 700,400", output="Figure3_high.png")
# @gp
# end
# # 3d print commands
# function plot_npi_panel()
# comb = readdlm("mfdata/attackrates_24.dat")
# comb = comb .* 10000
# comb = comb ./ 9000
# qvals = collect(1:size(comb)[1]) .+ 29
# gvals = collect(1:size(comb)[2]) .+ 29
# # reduce qvals
# qvals = qvals[qvals .<= 150]
# comb = comb[1:length(qvals), :]
# # calculate percentage increase from the baseline (i.e. vaccine on and no reduction in NPI)
# baseline_ar = get_compartments(results_one(9999, -1, 20, 0.0215)).W[end] / 9000
# println("baseline: $baseline_ar")
# comb = round.((comb .- baseline_ar) ./ baseline_ar * 100, digits=1)
# @gp "reset"
# @gp :- "set size ratio 0.5"
# @gp :- "set lmargin at screen 0.0"
# @gp :- "set rmargin at screen 0.9"
# #@gp :- "set bmargin at screen 0.1"
# #@gp :- "set tmargin at screen 0.95"
# @gp :- "set style increment user"
# @gp :- "set style line 1 lw 3" # the width of the contour line uses line style 1 by ult
# @gp :- "set style line 2 lc rgb 'white' " # contour lines uses line style 2 by ult.
# #@gp :- "set contour"
# #@gp :- "set cntrparam cubicspline" # Smooth out the lines"
# #@gp :- "set cntrparam levels incremental 0,100,2000"
# #@gp :- "set cntrlabel onecolor"
# #@gp :- "set cntrparam levels discrete 5" # Plot the selected contours
# #@gp :- "set cntrlabel start 5 interval 20"
# #set cntrlabel format '%5.3g' font ',7'
# #@gp :- "set xrange [1:5]"
# #@gp :- "set yrange [1:5]"
# #set urange [ 5 : 35 ] noreverse nowriteback
# #set vrange [ 5 : 35 ] noreverse nowriteback
# @gp :- "set view map" # bird eye view of 3d plot
# @gp :- "set xtics 30,10,150"
# @gp :- "set dgrid3d 10,10,10"
# @gp :- "set grid x y lc 'black' front"
# @gp :- "unset border"
# # Note that grid lines only appear at axis tic locations.
# # But black tics would hide the white grid lines, so scale them to 0 length
# @gp :- "set tics scale 0"
# @gp :- "set pm3d interpolate 0,0" # where "x" and "y" are the number of additional points in the x and y axes. So if "x" is set to 4 then Gnuplot will interpolate with 4 times as many points in the x-axis.
# ## 0,0 is automatic interpolation
# #@gp :- "set palette defined (0 0 0 0.5, 1 0 0 1, 2 0 0.5 1, 3 0 1 1, 4 0.5 1 0.5, 5 1 1 0, 6 1 0.5 0, 7 1 0 0, 8 0.5 0 0)"
# #@gp :- xlab = "contact tracing" ylab = "y label" palette(:thermal)
# #@gp :- "set format y '%.f%%'"
# #@gp :- "set format x '%.f%%'"
# @gp :- "set format cb '%g%%'"
# @gp :- "set palette rgbformulae 33,13,10"
# @gp :- "set ylabel 'Rate of Reduction (in days)' offset 1,0,0"
# @gp :- "set xlabel 'Starting Day of NPI Reduction'"
# @gp :- "set ytics offset 1.5,0,0"
# @gp :- "set cblabel 'Percent Increase in Attack Rate' offset 1,0,0"
# @gp :- "set style textbox opaque noborder"
# #@gp :- "set arrow from 80,30 to 80,90 nohead lw 2 linetype 2 dashtype 2 lc rgb 'black' front";
# @gsp :- qvals gvals comb "w pm3d notitle ls 1"
# #@gsp :- qvals gvals comb "w labels notitle"
# #@gsp :- qvals gvals comb " with lines lw 4 nosurface"
# # dump( Gnuplot.options.term)
# #save(term="pdfcairo enhanced font 'Arial,12' size 8.5,4.5", output="plot_heatmap.pdf")
# save(term="pngcairo enhanced font 'Arial,10' size 700,400", output="plot_heatmap_deaths.png")
# end
# # 3d print commands
# function plot_comb()
# comb = readdlm("comb.dat") ./ 100
# qvals = collect(2:50)
# gvals = collect(1:100)
# maxqval = qvals[qvals .<= 50]
# maxgval = gvals[gvals .<= 100]
# #filter the comb
# comb = comb[1:length(maxqval), 1:length(maxgval)]
# longcomb = zeros(Float64, length(vec(comb)), 3)
# longcomb[:, 1] .= repeat(maxqval, Int(length(vec(comb))/length(maxqval)))
# longcomb[:, 2] .= sort(repeat(maxgval, Int(length(vec(comb))/length(maxgval))))
# longcomb[:, 3] .= vec(comb)
# # we don't really use longcomb anymore.
# @gp "reset"
# #@gp :- "set term qt font 'Arial,12'"
# #@gp :- "set terminal pdfcairo font 'Helvetica,12' size 1cm,1cm fontscale 0.75"
# @gp :- "set size ratio 0.5"
# # gnuplot and pm3d don't play well with margins
# @gp :- "set lmargin at screen 0.0"
# #@gp :- "set rmargin at screen 0.9"
# #@gp :- "set bmargin at screen 0.1"
# #@gp :- "set tmargin at screen 0.95"
# # see https://stackoverflow.com/questions/18878163/gnuplot-contour-line-color-set-style-line-and-set-linetype-not-working
# # simple heatmap example gnuplot http://www.phyast.pitt.edu/~zov1/gnuplot/html/contour.html
# # more splot example http://www.gnuplotting.org/tag/splot/
# # plot line over heatmap https://stackoverflow.com/questions/43917389/how-to-plot-new-line-given-by-math-formula-over-pm3d-map
# @gp :- "set style increment user"
# @gp :- "set style line 1 lw 3" # the width of the contour line uses line style 1 by ult
# @gp :- "set style line 2 lc rgb 'white' " # contour lines uses line style 2 by ult.
# @gp :- "set contour"
# #@gp :- "set cntrparam cubicspline" # Smooth out the lines"
# #@gp :- "set cntrparam levels incremental 0,100,2000"
# #@gp :- "set cntrlabel onecolor"
# @gp :- "set cntrparam levels discrete 5" # Plot the selected contours
# #@gp :- "set cntrlabel start 5 interval 20"
# #set cntrlabel format '%5.3g' font ',7'
# #@gp :- "set xrange [1:5]"
# #@gp :- "set yrange [1:5]"
# #set urange [ 5 : 35 ] noreverse nowriteback
# #set vrange [ 5 : 35 ] noreverse nowriteback
# @gp :- "set view map" # bird eye view of 3d plot
# @gp :- "set xtics 0,5,50"
# @gp :- "set dgrid3d 10,10,10"
# @gp :- "set grid x y lc 'black' front"
# @gp :- "unset border"
# # Note that grid lines only appear at axis tic locations.
# # But black tics would hide the white grid lines, so scale them to 0 length
# @gp :- "set tics scale 0"
# @gp :- "set pm3d interpolate 0,0" # where "x" and "y" are the number of additional points in the x and y axes. So if "x" is set to 4 then Gnuplot will interpolate with 4 times as many points in the x-axis.
# ## 0,0 is automatic interpolation
# #@gp :- "set palette defined (0 0 0 0.5, 1 0 0 1, 2 0 0.5 1, 3 0 1 1, 4 0.5 1 0.5, 5 1 1 0, 6 1 0.5 0, 7 1 0 0, 8 0.5 0 0)"
# #@gp :- xlab = "contact tracing" ylab = "y label" palette(:thermal)
# #@gp :- "set format y '%.f%%'"
# #@gp :- "set format x '%.f%%'"
# @gp :- "set palette rgbformulae 33,13,10"
# @gp :- "set ylabel '% of silent infection identified' offset 2,0,0"
# @gp :- "set xlabel '% of exposed contacts traced'"
# @gp :- "set ytics offset 1.5,0,0"
# @gp :- "set cblabel 'Attack Rate'"
# @gp :- "set style textbox opaque noborder"
# @gsp :- qvals gvals comb "w pm3d notitle ls 1"
# #@gsp :- qvals gvals comb "w labels notitle"
# #@gsp :- qvals gvals comb " with lines lw 4 nosurface"
# # dump( Gnuplot.options.term)
# #save(term="pdfcairo enhanced font 'Arial,12' size 8.5,4.5", output="plot_heatmap.pdf")
# save(term="pngcairo enhanced font 'Arial,12' size 700,400", output="plot_heatmap.png")
# end
# # cubic spline
# function plot_comb_gp()
# # same function as plot_comb but uses @gp (problem with @gp is lack of contour lines)
# # but this function makes use of Interpolation package.
# comb = readdlm("comb.dat")
# qvals = collect(2:50)
# gvals = collect(1:100)
# #longcomb = zeros(Float64, length(vec(comb)), 3)
# #longcomb[:, 1] .= repeat(qvals, length(gvals))
# #longcomb[:, 2] .= repeat(gvals, length(qvals))
# #longcomb[:, 3] .= vec(comb)
# #@gp "reset session"
# #@gp qvals gvals comb "w image notit" "set view map" "set auto fix" #"set size ratio -1"
# #@gp :- xlab = "x" ylab = "y" "set cblabel 'attack rate'" palette(:thermal)
# cubint = CubicSplineInterpolation((1:49,1:100),comb)
# xi, yi = 1.0:0.2:49, 1.0:0.2:100
# fxyi = [cubint(u,v) for u = xi, v = yi]
# @gp xi yi fxyi "w image notit" "set view map" "set auto fix"
# @gp :- xlab = "contact tracing" ylab = "silent infection" palette(:thermal)
# end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment