# Building an Epidemic Movie
#
# Accompanies the SISMID course "Network Modeling for Epidemics" (NME).
# Source page: mod5-Movie.qmd
#
# The baseline and reference implementation for Module 5. We build one epidemic
# on one simple edges-only network, walk the whole pipeline from an empty
# network to an animated movie, and then look at the static transmission-tree
# views. Think of the 100 nodes as a small community (a residence hall, a
# workplace, a village) and the edges as the recurring close contacts among
# them: the specific pairs who spend enough time together for a close-contact
# infection to pass. Most people have only a handful of such contacts at a time.
#
# Each render below writes an HTML movie to your working directory (and, run
# interactively, opens it in your browser). Run top to bottom.


# ---- Setup ----

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

# A seed is set before each model in this module (not once at the top), so each
# model is reproducible on its own and inserting or re-running one cannot
# silently change the others. Skip the seed lines for different draws.


# ---- The network and epidemic model ----

# Edges-only: the only thing we specify is HOW MANY ties, 40 among 100 nodes,
# with a mean partnership duration of 20 steps. target.stats = 40 is an EXPECTED
# value: realized edge counts fluctuate around it, typically in a band a few
# edges wide. Mean degree is 2 * 40 / 100 = 0.8.
set.seed(1234)
nw <- network_initialize(n = 100)
formation <- ~edges
target.stats <- 40
coef.diss <- dissolution_coefs(dissolution = ~offset(edges), duration = 20)
est <- netest(nw, formation, target.stats, coef.diss)

# Check the network model hit its target before simulating. netdx takes its OWN
# nsims/nsteps, independent of the epidemic: the movie needs a small network
# over a short horizon, the diagnostic does not. ncores = 5 is safe anywhere;
# EpiModel caps it at your actual core count, so set ncores = 1 only if you want
# to force fully serial execution.
dx <- netdx(est, nsims = 5, nsteps = 500, ncores = 5)
dx

# The epidemic is the poker-chip SI process: a 100% per-act transmission
# probability, so every S-I contact transmits and the movie shows the NETWORK's
# limits on spread, not the disease's. One simulation, 25 steps, 10 seeds.
param <- param.net(inf.prob = 1)
init <- init.net(i.num = 10)
control <- control.net(type = "SI", nsteps = 25, nsims = 1)
sim <- netsim(est, param, init, control)


# ---- Building the movie: the four-step pipeline ----

# Step 1: extract the networkDynamic object, the network WITH timing attached
# (when each edge is active, how each node's status changes over the run).
nw <- get_network(sim)
nw

# Step 2: color nodes by disease status over time. color_tea() attaches a
# temporally extended attribute (TEA) named "ndtvcol", which is why every render
# passes vertex.col = "ndtvcol".
nw <- color_tea(nw, verbose = FALSE)

# Step 3: compute the layout. Only slice.par feeds compute.animation (which
# steps to show, how to bin dynamic edges into frames); render.par (playback)
# and plot.par (margins) are used later at render time. The default Kamada-Kawai
# layout has no meaningful coordinates (only relative position matters) and is
# chained frame to frame so nodes do not jump.
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))
compute.animation(nw, slice.par = slice.par, verbose = FALSE)

# Step 4: render. Here we WRITE an HTML file (the tutorial page embeds inline
# instead, with output.mode = "htmlWidget").
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(), "movie-1-baseline.html"))


# ---- Optional: transmission trees ----

# The static views of who infected whom. The transmission matrix records each
# infection event (who, by whom, when).
tm <- get_transmat(sim)
head(tm)

# Two whole-matrix plots. Left: a directed who-infected-whom network; with 10
# seeds it is NOT one component (each seed founds its own chain). Right: the
# transmission timeline. Its HORIZONTAL axis is model (calendar) time. Its
# VERTICAL axis, despite the "generation" label ndtv prints on it, is NOT a
# transmission generation: it is just a stacking counter that gives each
# infected node its own height (in ndtv's traversal order, 0 up to about the
# number infected) so points do not overlap. It carries no epidemiologic
# meaning: read only horizontal position (infection time) and the connecting
# lines (who infected whom).
par(mfrow = c(1, 2), mar = c(3, 3, 2, 1))
plot(tm, style = "network", displaylabels = TRUE)
plot(tm, style = "transmissionTimeline")

# Phylograms. This is the model's EXACT, KNOWN transmission tree (from
# get_transmat), not an estimated molecular phylogeny. Convert with
# as.phylo.transmat and index into the result. There are fewer trees than seeds
# (a seed that never transmits leaves no tree), and the trees vary enormously in
# size from identical starting conditions, purely by where each seed landed:
# the same reason a single run demonstrates but cannot measure.
tmPhylo <- as.phylo.transmat(tm)
sizes <- sapply(tmPhylo, function(x) length(x$tip.label))
sizes

# In each tree, only the TIPS are the distinct infected individuals (their count
# is the "N infected" in the panel title). The numbers at the internal branch
# points repeat the infector at each transmission event, so do not read the
# repeats as extra people.
biggest <- order(sizes, decreasing = TRUE)[1:min(4, length(sizes))]
par(mfrow = c(2, 2), mar = c(1, 1, 2, 1))
for (i in biggest) {
  plot(tmPhylo[[i]], show.node.label = TRUE, root.edge = TRUE, cex = 0.5,
       main = paste0(names(tmPhylo)[i], ": ", sizes[i], " infected"))
}