# =============================================================================
# Appendix: How Fast Is the Network?
# Network Modeling for Epidemics (NME)
#
# The master variable of this module is the ratio of two clocks: how fast the
# network turns over, versus how fast the epidemic moves. This script develops
# three tools for reading that ratio:
#   (1) momentary vs cumulative degree, a diagnostic for where you sit on it;
#   (2) the epidemic threshold, the result that lives at the slow-network limit;
#   (3) temporal reachability, which shows the threshold dissolving at the fast limit.
# The movies write to standalone .html files in your working directory; open them
# in a browser. On a single-core machine, change every ncores = 5 to ncores = 1.
# =============================================================================

library("EpiModel")
library("ndtv")
library("sna")    # component.dist(): connected components of a static graph
library("tsna")   # tReach(): reachability by time-respecting paths

# ---- Animation helper (same four-step pipeline as the other tutorials) -------
# Writes a standalone HTML movie to `file` in the working directory.
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 fixed poker-chip SI epidemic used for the movies (certain transmission,
# 10 seeds, 25 steps), so every difference we see comes from the network.
movie.param   <- param.net(inf.prob = 1)
movie.init    <- init.net(i.num = 10)
movie.control <- control.net(type = "SI", nsteps = 25, nsims = 1)


# =============================================================================
# 1. Diagnostic: momentary vs cumulative degree
# =============================================================================
# network.collapse() flattens a dynamic network to a static one. At a single
# instant (at =) it gives the momentary network; over a window (onset, terminus)
# it gives the cumulative network, where an edge exists if the pair was EVER
# tied in that window. degree() then counts partners in each.

degree_profile <- function(duration, horizon = 50, seed = 1234) {
  set.seed(seed)
  nw  <- network_initialize(100)
  est <- netest(nw, ~edges, 50,
                dissolution_coefs(~offset(edges), duration = duration),
                verbose = FALSE)
  # We want the contact network's OWN dynamics with no epidemic on it, so we
  # simulate it with netdx (keep.tnetwork = TRUE keeps the full networkDynamic)
  # rather than netsim, then pull that dynamic network out to collapse below.
  dx  <- netdx(est, nsims = 1, nsteps = horizon, dynamic = TRUE,
               keep.tnetwork = TRUE, verbose = FALSE)
  nwd <- get_network(dx, sim = 1)
  data.frame(
    momentary  = degree(network.collapse(nwd, at = horizon), gmode = "graph"),
    cumulative = degree(network.collapse(nwd, onset = 1, terminus = horizon + 1),
                        gmode = "graph"))
}

# Same mean degree (1), same horizon; only the network clock differs.
# duration 200 = slow (partnerships outlast the horizon); duration 2 = fast
# (about 25 turnovers during it).
slow <- degree_profile(duration = 200)
fast <- degree_profile(duration = 2)

# Compare the two degree measures: are they the same number, and the same people?
# top15_overlap = how many of the 15 people you'd vaccinate agree between measures.
compare <- function(d) {
  top <- function(x) order(x, decreasing = TRUE)[1:15]
  data.frame(
    mean_momentary  = round(mean(d$momentary), 2),
    mean_cumulative = round(mean(d$cumulative), 2),
    correlation     = round(cor(d$momentary, d$cumulative), 2),
    top15_overlap   = length(intersect(top(d$momentary), top(d$cumulative))))
}
rbind(slow = compare(slow), fast = compare(fast))
# Slow: momentary ~ cumulative, correlation ~0.9, lists share ~13/15.
# Fast: momentary reads ~0.7 in this single snapshot (target 1) but cumulative
#   ~22, correlation ~0, lists share ~3/15.
#
# This is P = I x D applied to the network: momentary degree is the PREVALENCE
# of partnerships, duration is D, and their ratio is the INCIDENCE (turnover).
# Cumulative degree is what that incidence delivers over the horizon.


# =============================================================================
# 2. Slow limit: the epidemic threshold on a frozen network
# =============================================================================
# When the network clock is much slower than the epidemic, the network is frozen
# and reach is capped by the connected component a seed lands in. Sweep the mean
# degree and watch the giant component (and take-off) appear near mean degree 1.

threshold_row <- function(edges, n = 100, nsims = 100, nsteps = 80) {
  set.seed(1234)
  nw  <- network_initialize(n)
  est <- netest(nw, ~edges, edges,
                dissolution_coefs(~offset(edges), duration = 1e5), verbose = FALSE)
  # est$fit is the fitted cross-sectional ERGM stored inside the netest object;
  # simulate()ing from it draws static networks from the same ensemble the
  # epidemic runs on. Largest connected component (fraction), over 20 draws.
  giant <- mean(replicate(20,
    max(component.dist(simulate(est$fit, output = "network"),
                       connected = "weak")$csize) / n))
  # Epidemic final size from a single seed with certain transmission.
  sim <- netsim(est, param.net(inf.prob = 1), init.net(i.num = 1),
                control.net(type = "SI", nsteps = nsteps, nsims = nsims,
                            ncores = 5, verbose = FALSE))
  final <- as.data.frame(sim)$i.num[as.data.frame(sim)$time == nsteps]
  data.frame(mean_degree     = 2 * edges / n,
             giant_component = round(giant, 2),
             mean_final_size = round(mean(final), 1),
             p_takeoff       = round(mean(final > 10), 2))
}
sweep <- do.call(rbind, lapply(c(20, 40, 50, 60, 80, 100, 120), threshold_row))
sweep

# Both the giant component and the epidemic lift off around mean degree 1.
plot(sweep$mean_degree, sweep$giant_component, type = "b", pch = 19,
     col = "steelblue", ylim = c(0, 1),
     xlab = "mean degree", ylab = "fraction of the population")
lines(sweep$mean_degree, sweep$mean_final_size / 100, type = "b",
      pch = 17, col = "firebrick")
abline(v = 1, lty = 2, col = "grey50")
legend("topleft", bty = "n", pch = c(19, 17),
       col = c("steelblue", "firebrick"),
       legend = c("giant component", "mean final size / 100"))

# A frozen network is a movie where only the colors change, never the edges.
# At mean degree 0.8 (below the giant-component threshold) each seed saturates
# its small cluster and stops: the infection runs out of network.
set.seed(1234)
nw <- network_initialize(100)
est.frozen <- netest(nw, ~edges, 40,
                     dissolution_coefs(~offset(edges), duration = 1e5), verbose = FALSE)
sim.frozen <- netsim(est.frozen, movie.param, movie.init, movie.control)
animate(sim.frozen, "movie-4-frozen.html")


# =============================================================================
# 3. Fast limit: temporal reachability
# =============================================================================
# Same momentary mean degree (0.8), but duration 2 instead of frozen. The right
# question is now "whom can I reach by a time-respecting path" (contacts that
# move forward in time), which tReach() computes. Compare it to two static
# baselines: a single snapshot, and the cumulative graph (the reachability ceiling).
set.seed(1234)
nw <- network_initialize(100)
est.churn <- netest(nw, ~edges, 40,
                    dissolution_coefs(~offset(edges), duration = 2), verbose = FALSE)
dx  <- netdx(est.churn, nsims = 1, nsteps = 50, dynamic = TRUE,
             keep.tnetwork = TRUE, verbose = FALSE)
nwd <- get_network(dx, sim = 1)

snapshot_giant     <- max(component.dist(network.collapse(nwd, at = 25),
                                         connected = "weak")$csize)
temporal_reach     <- mean(tReach(nwd, direction = "fwd"))
cumulative_ceiling <- max(component.dist(network.collapse(nwd, onset = 1, terminus = 51),
                                         connected = "weak")$csize)
c(snapshot_giant     = snapshot_giant,       # ~8: what one static picture predicts
  temporal_reach     = round(temporal_reach),# ~100: what the churn delivers
  cumulative_ceiling = cumulative_ceiling)   # 100: the union-graph ceiling

# Same momentary degree as the frozen run, opposite reach: the epidemic threshold
# is a property of the slow-network limit, not a law of the network.
sim.churn <- netsim(est.churn, movie.param, movie.init, movie.control)
animate(sim.churn, "movie-4-churn.html")

# Takeaway: ask the clock question (and the density question) first. How fast
# does the network turn over relative to the epidemic, and how dense is it
# relative to the percolation threshold? That tells you whether you are in a
# percolation problem or a concurrency problem before you commit to either.