Skip to content

Commit

Permalink
updated lectures
Browse files Browse the repository at this point in the history
  • Loading branch information
@annabellasd
annabellasd committed Sep 26, 2023
1 parent d5bdc9b commit 0dbf48d
Show file tree
Hide file tree
Showing 2 changed files with 103 additions and 99 deletions.
121 changes: 62 additions & 59 deletions lectures/continuous_time/covid_sde.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -230,40 +230,41 @@ First, construct our $F$ from {eq}`dfcvsde` and $G$ from {eq}`dG`

```{code-cell} julia
function F(x, p, t)
s, i, r, d, R₀, δ = x
(;γ, R̄₀, η, σ, ξ, θ, δ_bar) = p
return [-γ*R₀*s*i; # ds/dt
γ*R₀*s*i - γ*i; # di/dt
(1-δ)*γ*i; # dr/dt
δ*γ*i; # dd/dt
η*(R̄₀(t, p) - R₀);# dR₀/dt
θ*(δ_bar - δ); # dδ/dt
]
s, i, r, d, R_0, delta = x
(; gamma, R_bar_0, eta, sigma, xi, theta, delta_bar) = p
return [-gamma * R_0 * s * i; # ds/dt
gamma * R_0 * s * i - gamma * i; # di/dt
(1 - delta) * gamma * i; # dr/dt
delta * gamma * i; # dd/dt
eta * (R_bar_0(t, p) - R_0); # dR_0/dt
theta * (delta_bar - delta); # ddelta/dt
]
end
function G(x, p, t)
s, i, r, d, R₀, δ = x
(;γ, R̄₀, η, σ, ξ, θ, δ_bar) = p
s, i, r, d, R_0, delta = x
(; gamma, R_bar_0, eta, sigma, xi, theta, delta_bar) = p
return [0; 0; 0; 0; σ*sqrt(R₀); ξ*sqrt(δ * (1-δ))]
return [0; 0; 0; 0; sigma * sqrt(R_0); xi * sqrt(delta * (1 - delta))]
end
```

Next create a settings generator, and then define a [SDEProblem](https://docs.sciml.ai/stable/tutorials/sde_example/#Example-2:-Systems-of-SDEs-with-Diagonal-Noise-1) with Diagonal Noise.

```{code-cell} julia
p_gen = @with_kw (T = 550.0, γ = 1.0 / 18, η = 1.0 / 20,
R₀_n = 1.6, R̄₀ = (t, p) -> p.R₀_n, δ_bar = 0.01,
σ = 0.03, ξ = 0.004, θ = 0.2, N = 3.3E8)
p = p_gen() # use all defaults
p_gen(; T = 550.0, gamma = 1.0 / 18, eta = 1.0 / 20,
R_0_n = 1.6, R_bar_0 = (t, p) -> p.R_0_n, delta_bar = 0.01,
sigma = 0.03, xi = 0.004, theta = 0.2, N = 3.3E8) = (; T, gamma, eta, R_0_n, R_bar_0, delta_bar, sigma, xi, theta, N)
p = p_gen() # use all defaults
i_0 = 25000 / p.N
r_0 = 0.0
d_0 = 0.0
s_0 = 1.0 - i_0 - r_0 - d_0
R̄₀_0 = 0.5 # starting in lockdown
δ_0 = p.δ_bar
x_0 = [s_0, i_0, r_0, d_0, R̄₀_0, δ_0]
R_bar_0_0 = 0.5 # starting in lockdown
delta_0 = p.delta_bar
x_0 = [s_0, i_0, r_0, d_0, R_bar_0_0, delta_0]
prob = SDEProblem(F, G, x_0, (0, p.T), p)
```
Expand Down Expand Up @@ -370,29 +371,29 @@ We will shut down the shocks to the mortality rate (i.e. $\xi = 0$) to focus on
Consider $\eta = 1/50$ and $\eta = 1/20$, where we start at the same initial condition of $R_0(0) = 0.5$.

```{code-cell} julia
function generate_η_experiment(η; p_gen = p_gen, trajectories = 100,
saveat = 1.0, x_0 = x_0, T = 120.0)
p = p_gen(η = η, ξ = 0.0)
function generate_eta_experiment(eta; p_gen, trajectories = 100, saveat = 1.0, x_0, T = 120.0)
p = p_gen(eta = eta, xi = 0.0)
ensembleprob = EnsembleProblem(SDEProblem(F, G, x_0, (0, T), p))
sol = solve(ensembleprob, SOSRI(), EnsembleThreads(),
trajectories = trajectories, saveat = saveat)
sol = solve(ensembleprob, SOSRI(), EnsembleThreads(), trajectories = trajectories, saveat = saveat)
return EnsembleSummary(sol)
end
# Evaluate two different lockdown scenarios
η_1 = 1/50
η_2 = 1/20
summ_1 = generate_η_experiment(η_1)
summ_2 = generate_η_experiment(η_2)
plot(summ_1, idxs = (4,5),
eta_1 = 1 / 50
eta_2 = 1 / 20
summ_1 = generate_eta_experiment(eta_1)
summ_2 = generate_eta_experiment(eta_2)
plot(summ_1, idxs = (4, 5),
title = ["Proportion Dead" L"R_0"],
ylabel = [L"d(t)" L"R_0(t)"], xlabel = L"t",
legend = [:topleft :bottomright],
labels = L"Middle 95% Quantile, $\eta = %$η_1$",
labels = L"Middle 95% Quantile, $\eta = %$eta_1$",
layout = (2, 1), size = (900, 900), fillalpha = 0.5)
plot!(summ_2, idxs = (4,5),
plot!(summ_2, idxs = (4, 5),
legend = [:topleft :bottomright],
labels = L"Middle 95% Quantile, $\eta = %$η_2$", size = (900, 900), fillalpha = 0.5)
labels = L"Middle 95% Quantile, $\eta = %$eta_2$", size = (900, 900), fillalpha = 0.5)
```

While the the mean of the $d(t)$ increases, unsurprisingly, we see that the 95% quantile for later time periods is also much larger - even after the $R_0$ has converged.
Expand All @@ -411,24 +412,24 @@ Since empirical estimates of $R_0(t)$ discussed in {cite}`NBERw27128` and other
We start the model with 100,000 active infections.

```{code-cell} julia
R₀_L = 0.5 # lockdown
η_experiment = 1.0/10
σ_experiment = 0.04
R̄₀_lift_early(t, p) = t < 30.0 ? R₀_L : 2.0
R̄₀_lift_late(t, p) = t < 120.0 ? R₀_L : 2.0
R_0_L = 0.5 # lockdown
eta_experiment = 1.0 / 10
sigma_experiment = 0.04
p_early = p_gen(R̄₀ = R̄₀_lift_early, η = η_experiment, σ = σ_experiment)
p_late = p_gen(R̄₀ = R̄₀_lift_late, η = η_experiment, σ = σ_experiment)
R_bar_0_lift_early(t, p) = t < 30.0 ? R_0_L : 2.0
R_bar_0_lift_late(t, p) = t < 120.0 ? R_0_L : 2.0
p_early = p_gen(R_bar_0 = R_bar_0_lift_early, eta = eta_experiment, sigma = sigma_experiment)
p_late = p_gen(R_bar_0 = R_bar_0_lift_late, eta = eta_experiment, sigma = sigma_experiment)
# initial conditions
i_0 = 100000 / p_early.N
r_0 = 0.0
d_0 = 0.0
s_0 = 1.0 - i_0 - r_0 - d_0
δ_0 = p_early.δ_bar
delta_0 = p_early.delta_bar
x_0 = [s_0, i_0, r_0, d_0, R₀_L, δ_0] # start in lockdown
x_0 = [s_0, i_0, r_0, d_0, R_0_L, delta_0] # start in lockdown
prob_early = SDEProblem(F, G, x_0, (0, p_early.T), p_early)
prob_late = SDEProblem(F, G, x_0, (0, p_late.T), p_late)
```
Expand All @@ -438,13 +439,13 @@ Simulating for a single realization of the shocks, we see the results are qualit
```{code-cell} julia
sol_early = solve(prob_early, SOSRI())
sol_late = solve(prob_late, SOSRI())
plot(sol_early, vars = [5, 1,2,4],
plot(sol_early, vars = [5, 1, 2, 4],
title = [L"R_0" "Susceptible" "Infected" "Dead"],
layout = (2, 2), size = (900, 600),
ylabel = [L"R_0(t)" L"s(t)" L"i(t)" L"d(t)"], xlabel = L"t",
legend = [:bottomright :topright :topright :topleft],
label = ["Early" "Early" "Early" "Early"])
plot!(sol_late, vars = [5, 1,2,4],
plot!(sol_late, vars = [5, 1, 2, 4],
legend = [:bottomright :topright :topright :topleft],
label = ["Late" "Late" "Late" "Late"], xlabel = L"t")
```
Expand Down Expand Up @@ -535,24 +536,26 @@ We will redo the "Ending Lockdown" simulation from above, where the only differe

```{code-cell} julia
function F_reinfect(x, p, t)
s, i, r, d, R₀, δ = x
(;γ, R̄₀, η, σ, ξ, θ, δ_bar, ν) = p
return [-γ*R₀*s*i + ν*r; # ds/dt
γ*R₀*s*i - γ*i; # di/dt
(1-δ)*γ*i - ν*r # dr/dt
δ*γ*i; # dd/dt
η*(R̄₀(t, p) - R₀);# dR₀/dt
θ*(δ_bar - δ); # dδ/dt
]
s, i, r, d, R_0, delta = x
(;gamma, R_bar_0, eta, sigma, xi, theta, delta_bar, nu) = p
return [
-gamma * R_0 * s * i + nu * r; # ds/dt
gamma * R_0 * s * i - gamma * i; # di/dt
(1 - delta) * gamma * i - nu * r; # dr/dt
delta * gamma * i; # dd/dt
eta * (R_bar_0(t, p) - R_0); # dR_0/dt
theta * (delta_bar - delta); # ddelta/dt
]
end
p_re_gen = @with_kw ( T = 550.0, γ = 1.0 / 18, η = 1.0 / 20,
R₀_n = 1.6, R̄₀ = (t, p) -> p.R₀_n,
δ_bar = 0.01, σ = 0.03, ξ = 0.004, θ = 0.2, N = 3.3E8, ν = 1/360)
p_re_gen(; T = 550.0, gamma = 1.0 / 18, eta = 1.0 / 20,
R_0_n = 1.6, R_bar_0 = (t, p) -> p.R_0_n,
delta_bar = 0.01, sigma = 0.03, xi = 0.004, theta = 0.2,
N = 3.3E8, nu = 1/360) = (; T, gamma, eta, R_0_n, R_bar_0, delta_bar, sigma, xi, theta, N, nu)
p_re_early = p_re_gen(R̄₀ = R̄₀_lift_early, η = η_experiment, σ = σ_experiment)
p_re_late = p_re_gen(R̄₀ = R̄₀_lift_late, η = η_experiment, σ = σ_experiment)
p_re_early = p_re_gen(R_bar_0 = R_bar_0 _lift_early, eta = eta_experiment, sigma = sigma_experiment)
p_re_late = p_re_gen(R_bar_0 = R_bar_0 _lift_late, eta = eta_experiment, sigma = sigma_experiment)
trajectories = 400
saveat = 1.0
Expand Down
81 changes: 41 additions & 40 deletions lectures/continuous_time/seir_model.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -171,13 +171,13 @@ We begin by implementing a simple version of this model with a constant $R_0$ an
First, define the system of equations

```{code-cell} julia
function F_simple(x, p, t; γ = 1/18, R₀ = 3.0, σ = 1/5.2)
function F_simple(x, p, t; gamma = 1/18, R_0 = 3.0, sigma = 1/5.2)
s, e, i, r = x
return [-γ*R₀*s*i; # ds/dt = -γR₀si
γ*R₀*s*i - σ*e;# de/dt = γR₀si -σe
σ*e - γ*i; # di/dt = σe -γi
γ*i; # dr/dt = γi
return [-gamma * R_0 * s * i; # ds/dt
gamma * R_0 * s * i - sigma * e; # de/dt
sigma * e - gamma * i; # di/dt
gamma * i; # dr/dt
]
end
```
Expand All @@ -195,7 +195,7 @@ tspan = (0.0, 350.0) # ≈ 350 days
prob = ODEProblem(F_simple, x_0, tspan)
```

With this, choose an ODE algorithm and solve the initial value problem. A good default algorithm for non-stiff ODEs of this sort might be `Tsit5()`, which is the Tsitouras 5/4 Runge-Kutta method).
With this, choose an ODE algorithm and solve the initial value problem. A good default algorithm for non-stiff ODEs of this sort might be `Tsit5()`, which is the Tsitouras 5/4 Runge-Kutta method.

```{code-cell} julia
sol = solve(prob, Tsit5())
Expand Down Expand Up @@ -242,7 +242,7 @@ While we could integrate the deaths given the solution to the model ex-post, it

This is a common trick when solving systems of ODEs. While equivalent in principle to using the appropriate quadrature scheme, this becomes especially convenient when adaptive time-stepping algorithms are used to solve the ODEs (i.e. there is not a regular time grid). Note that when doing so, $d(0) = \int_0^0 \delta \gamma i(\tau) d \tau = 0$ is the initial condition.

The system {eq}`seir_system` and the supplemental equations can be written in vector form $x := [s, e, i, r, R₀, c, d]$ with parameter tuple $p := (\sigma, \gamma, \eta, \delta, \bar{R}_0(\cdot))$
The system {eq}`seir_system` and the supplemental equations can be written in vector form $x := [s, e, i, r, R_0, c, d]$ with parameter tuple $p := (\sigma, \gamma, \eta, \delta, \bar{R}_0(\cdot))$

Note that in those parameters, the targeted reproduction number, $\bar{R}_0(t)$, is an exogenous function.

Expand Down Expand Up @@ -289,35 +289,36 @@ First, construct our $F$ from {eq}`dfcv`

```{code-cell} julia
function F(x, p, t)
s, e, i, r, R₀, c, d = x
(;σ, γ, R̄₀, η, δ) = p
return [-γ*R₀*s*i; # ds/dt
γ*R₀*s*i - σ*e; # de/dt
σ*e - γ*i; # di/dt
γ*i; # dr/dt
η*(R̄₀(t, p) - R₀);# dR₀/dt
σ*e; # dc/dt
δ*γ*i; # dd/dt
s, e, i, r, R_0, c, d = x
(;sigma, gamma, R_0_bar, eta, delta) = p
return [-gamma * R_0 * s * i; # ds/dt
gamma * R_0 * s * i - sigma * e; # de/dt
sigma * e - γ*i; # di/dt
gamma * i; # dr/dt
eta * (R_bar_0(t, p) - R_0); # dR_0/dt
sigma * e; # dc/dt
delta * gamma * i; # dd/dt
]
end;
```

This function takes the vector `x` of states in the system and extracts the fixed parameters passed into the `p` object.

The only confusing part of the notation is the `R̄₀(t, p)` which evaluates the `p.R̄₀` at this time (and also allows it to depend on the `p` parameter).
The only confusing part of the notation is the `R_0(t, p)` which evaluates the `p.R_bar_0` at this time (and also allows it to depend on the `p` parameter).

### Parameters

#####HERE ENDED
The baseline parameters are put into a named tuple generator (see previous lectures using [Parameters.jl](https://github.com/mauro3/Parameters.jl)) with default values discussed above.

```{code-cell} julia
p_gen = @with_kw ( T = 550.0, γ = 1.0 / 18, σ = 1 / 5.2, η = 1.0 / 20,
R₀_n = 1.6, δ = 0.01, N = 3.3E8,
R̄₀ = (t, p) -> p.R₀_n);
p_gen(T = 550.0, gamma = 1.0 / 18, sigma = 1 / 5.2, eta = 1.0 / 20,
R_0_n = 1.6, delta = 0.01, N = 3.3E8,
R_0_bar = (t, p) -> p.R_0_n) = (T, gamma, sigma, eta, R_0_n, delta, N, R_0_bar)
```

Note that the default $\bar{R}_0(t)$ function always equals $R_{0n}$ -- a parameterizable natural level of $R_0$ used only by the `R̄₀` function
Note that the default $\bar{R}_0(t)$ function always equals $R_{0n}$ -- a parameterizable natural level of $R_0$ used only by the `R_hat_0` function

Setting initial conditions, we choose a fixed $s, i, e, r$, as well as $R_0(0) = R_{0n}$ and $m(0) = 0.01$

Expand All @@ -328,7 +329,7 @@ i_0 = 1E-7
e_0 = 4.0 * i_0
s_0 = 1.0 - i_0 - e_0
x_0 = [s_0, e_0, i_0, 0.0, p.R₀_n, 0.0, 0.0]
x_0 = [s_0, e_0, i_0, 0.0, p.R_0_n, 0.0, 0.0]
tspan = (0.0, p.T)
prob = ODEProblem(F, x_0, tspan, p)
```
Expand Down Expand Up @@ -372,9 +373,9 @@ Let's start with the case where $\bar{R}_0(t) = R_{0n}$ is constant.
We calculate the time path of infected people under different assumptions of $R_{0n}$:

```{code-cell} julia
R₀_n_vals = range(1.6, 3.0, length = 6)
sols = [solve(ODEProblem(F, x_0, tspan, p_gen(R₀_n = R₀_n)),
Tsit5(), saveat=0.5) for R₀_n in R₀_n_vals];
R_0_n_vals = range(1.6, 3.0, length = 6)
sols = [solve(ODEProblem(F, x_0, tspan, p_gen(R_0_n = R_0_n)),
Tsit5(), saveat=0.5) for R_0_n in R_0_n_vals];
```

Here we chose `saveat=0.5` to get solutions that were evenly spaced every `0.5`.
Expand All @@ -384,7 +385,7 @@ Changing the saved points is just a question of storage/interpolation, and does
Let's plot current cases as a fraction of the population.

```{code-cell} julia
labels = permutedims([L"R_0 = %$r" for r in R₀_n_vals])
labels = permutedims([L"R_0 = %$r" for r in R_0_n_vals])
infecteds = [sol[3,:] for sol in sols]
plot(infecteds, label=labels, legend=:topleft, lw = 2, xlabel = L"t",
ylabel = L"i(t)", title = "Current Cases")
Expand Down Expand Up @@ -413,21 +414,21 @@ In the simple case, where $\bar{R}_0(t) = R_{0n}$ is independent of the state, t

We will examine the case where $R_0(0) = 3$ and then it falls to $R_{0n} = 1.6$ due to the progressive adoption of stricter mitigation measures.

The parameter `η` controls the rate, or the speed at which restrictions are
The parameter `eta` controls the rate, or the speed at which restrictions are
imposed.

We consider several different rates:

```{code-cell} julia
η_vals = [1/5, 1/10, 1/20, 1/50, 1/100]
labels = permutedims([L"\eta = %$η" for η in η_vals]);
eta_vals = [1/5, 1/10, 1/20, 1/50, 1/100]
labels = permutedims([L"\eta = %$eta" for eta in eta_vals]);
```

Let's calculate the time path of infected people, current cases, and mortality

```{code-cell} julia
x_0 = [s_0, e_0, i_0, 0.0, 3.0, 0.0, 0.0]
sols = [solve(ODEProblem(F, x_0, tspan, p_gen(η=η)), Tsit5(), saveat=0.5) for η in η_vals];
sols = [solve(ODEProblem(F, x_0, tspan, p_gen(eta=eta)), Tsit5(), saveat=0.5) for eta in eta_vals];
```

Next, plot the $R_0$ over time:
Expand Down Expand Up @@ -468,21 +469,21 @@ The parameters considered here start the model with 25,000 active infections
and 75,000 agents already exposed to the virus and thus soon to be contagious.

```{code-cell} julia
R₀_L = 0.5 # lockdown
R̄₀_lift_early(t, p) = t < 30.0 ? R₀_L : 2.0
R̄₀_lift_late(t, p) = t < 120.0 ? R₀_L : 2.0
p_early = p_gen(R̄₀= R̄₀_lift_early, η = 10.0)
p_late = p_gen(R̄₀= R̄₀_lift_late, η = 10.0)
R_0_L = 0.5 # lockdown
R_bar_0_lift_early(t, p) = t < 30.0 ? R_0_L : 2.0
R_bar_0_lift_late(t, p) = t < 120.0 ? R_0_L : 2.0
p_early = p_gen(R_bar_0= R_bar_0_lift_early, eta = 10.0)
p_late = p_gen(R_bar_0= R_bar_0_lift_late, eta = 10.0)
# initial conditions
i_0 = 25000 / p_early.N
e_0 = 75000 / p_early.N
s_0 = 1.0 - i_0 - e_0
x_0 = [s_0, e_0, i_0, 0.0, R₀_L, 0.0, 0.0] # start in lockdown
x_0 = [s_0, e_0, i_0, 0.0, R_0_L, 0.0, 0.0] # start in lockdown
# create two problems, with rapid movement of R₀(t) towards R̄₀(t)
# create two problems, with rapid movement of R_0(t) towards R_hat_0(t)
prob_early = ODEProblem(F, x_0, tspan, p_early)
prob_late = ODEProblem(F, x_0, tspan, p_late)
```
Expand All @@ -501,7 +502,7 @@ plot!(sol_late, vars = [7], label = "Lift Late", xlabel = L"t")
Next we examine the daily deaths, $\frac{d D(t)}{dt} = N \delta \gamma i(t)$.

```{code-cell} julia
flow_deaths(sol, p) = p.N * p.δ * p.γ * sol[3,:]
flow_deaths(sol, p) = p.N * p.delta * p.gamma * sol[3,:]
plot(sol_early.t, flow_deaths(sol_early, p_early), title = "Flow Deaths", label = "Lift Early")
plot!(sol_late.t, flow_deaths(sol_late, p_late), label = "Lift Late", xlabel = L"t")
Expand Down

0 comments on commit 0dbf48d

Please sign in to comment.