# Concurrency and Duration
#
# Accompanies the SISMID course "Network Modeling for Epidemics" (NME).
# Source page: mod5-Timing.qmd
#
# Two features of a network's TIMING shape an epidemic more than almost anything
# else: whether partnerships overlap (concurrency) and how long they last
# (duration). We change each in turn on the edges-only baseline from
# mod5-Movie, and close with a single 200-run comparison across all four
# specifications, because a movie shows a mechanism but cannot measure an effect.
#
# Each render writes an HTML movie to your working directory. Run top to bottom.


# ---- Setup ----

library("EpiModel")
library("ndtv")

# Animation settings, and a helper that wraps the four-step pipeline from
# mod5-Movie so each model is one call. (The tutorial page embeds movies inline;
# here we write each to its own HTML file so you can compare them.)
slice.par <- list(start = 1, end = 25, interval = 1, aggregate.dur = 1, rule = "any")
render.par <- list(tween.frames = 10, show.time = FALSE)
plot.par <- list(mar = c(0, 0, 0, 0))

animate <- function(sim, file) {
  nw <- get_network(sim)
  nw <- color_tea(nw, verbose = FALSE)
  compute.animation(nw, slice.par = slice.par, verbose = FALSE)
  render.d3movie(
    nw, render.par = render.par, plot.par = plot.par,
    vertex.cex = 0.9, vertex.col = "ndtvcol",
    edge.col = "darkgrey", vertex.border = "lightgrey",
    displaylabels = FALSE, filename = file.path(getwd(), file))
}

# The epidemic is held fixed throughout: the same poker-chip SI process as
# mod5-Movie, so every difference comes from the network.
param <- param.net(inf.prob = 1)
init <- init.net(i.num = 10)
control <- control.net(type = "SI", nsteps = 25, nsims = 1)


# ---- The baseline ----

# Re-fit the edges-only baseline; every change below is measured against it.
set.seed(1234)
nw <- network_initialize(n = 100)
coef.diss <- dissolution_coefs(dissolution = ~offset(edges), duration = 20)
est.base <- netest(nw, ~edges, 40, coef.diss)
sim.base <- netsim(est.base, param, init, control)
animate(sim.base, "movie-baseline.html")


# ---- Concurrency ----

# concurrent counts nodes with 2+ partners at once. It speeds transmission (an
# infection can cross an overlapping partnership without waiting). We drive it to
# 0: strict serial monogamy. First gauge the chance expectation: at mean degree
# 0.8, about 19 of 100 nodes would be concurrent at random.
ppois(1, lambda = 2 * 40 / 100, lower.tail = FALSE) * 100

set.seed(1234)
nw <- network_initialize(n = 100)
est.noconc <- netest(nw, ~edges + concurrent, c(40, 0), coef.diss)
sim.noconc <- netsim(est.noconc, param, init, control)

# With no concurrency the edges form a matching (disjoint pairs), so mean degree
# is capped at 1.0. At the target of 40 edges, ABOUT 80 of 100 nodes are paired
# and ABOUT 20 isolated at any moment: expected values, and the realized counts
# fluctuate (typically mid-30s to low-40s edges). Serial monogamy slows spread:
# an infection must wait for one tie to dissolve before travelling along the next.
animate(sim.noconc, "movie-noconcurrency.html")


# ---- Duration ----
#
# Duration only means something relative to two timescales:
#   - the TIME STEP: a duration of 1 is memoryless (edges redrawn each step);
#     longer durations carry memory (the same partners persist).
#   - the INFECTIOUS PERIOD: fast turnover relative to it approaches mass action;
#     durations that outlast it confine transmission to current partners.
# This SI model has no recovery, so the infectious period is the whole 25-step
# run, and we look at duration against that horizon.

# -- Very long duration: a static network --

# Mean duration 100,000 fixes every edge for the whole run. Little spreads even
# at 100% transmission: the infection is trapped in whatever component a seed
# lands in (a static-network-model assumption: a fixed graph, reachable set
# fixed at time zero, concurrency maximal by construction). Boundary: those are
# properties of THIS idealization, and the clean degree-distribution thresholds
# hold for specific ensembles, not universally.
set.seed(1234)
nw <- network_initialize(n = 100)
est.long <- netest(nw, ~edges, 40, dissolution_coefs(~offset(edges), duration = 1e5))
sim.long <- netsim(est.long, param, init, control)
animate(sim.long, "movie-longduration.html")

# -- One-step duration: rapid mixing --

# Duration 1 is the shortest possible: every edge lasts one step, so each step is
# an independent draw of momentary contacts. Persistence is 1 - 1/1 = 0, so the
# coefficient is logit(0) = -Inf (the correct answer, not an error).
dissolution_coefs(dissolution = ~offset(edges), duration = 1)

set.seed(1234)
nw <- network_initialize(n = 100)
est.short <- netest(nw, ~edges, 40, dissolution_coefs(~offset(edges), duration = 1))
sim.short <- netsim(est.short, param, init, control)

# The epidemic saturates: partner identity is memoryless, so an infected person
# keeps meeting fresh susceptibles. This is close to what a compartmental (ODE)
# model assumes (beta*S*I/N has no partnerships). But that headline is a stylistic
# generalization, not an identity: this is a stochastic, discrete-time network,
# an analogue of the rapid-mixing limit, and rapid mixing does not NECESSARILY
# overstate spread (it depends on the disease, degree distribution, and outcome).
animate(sim.short, "movie-shortduration.html")


# ---- Replicated evidence ----

# The movies above are single draws. Run each specification 200 times and report
# the mean number infected at the final step. These numbers are specific to THESE
# specifications; read them as evidence for the DIRECTIONS the movies suggested.
final_prev <- function(formation, target.stats, duration, nsims = 200) {
  set.seed(1234)
  nw <- network_initialize(n = 100)
  cd <- dissolution_coefs(~offset(edges), duration = duration)
  est <- netest(nw, formation, target.stats, cd, verbose = FALSE)
  sim <- netsim(est, param, init,
                control.net(type = "SI", nsteps = 25, nsims = nsims,
                            ncores = 5, verbose = FALSE))
  d <- as.data.frame(sim, out = "mean")
  round(d$i.num[d$time == 25], 1)
}

data.frame(
  specification = c("Baseline (edges only, duration 20)",
                    "No concurrency (duration 20)",
                    "Very long duration (100,000)",
                    "One-step duration (1)"),
  mean_infected_of_100 = c(
    final_prev(~edges, 40, 20),
    final_prev(~edges + concurrent, c(40, 0), 20),
    final_prev(~edges, 40, 1e5),
    final_prev(~edges, 40, 1)))