#' # Multi-Layer Networks
#' 
## ----setupHide, include = FALSE-----------------------------------------------
library(knitr)
knitr::opts_chunk$set(comment = NA)

#' 
#' Real contact processes rarely happen on a single network. People have main partners and casual partners; a pathogen might spread over sexual contacts and needle-sharing contacts; students share dorm rooms and classrooms. A **multi-layer network** captures this: one shared set of people (the nodes), with two or more separate edge sets (the layers), each with its own formation and dissolution model. An epidemic can spread over the edges in either layer.
#' 
#' We build a two-layer SI epidemic in two steps, easiest first:
#' 
#' - **Part 1: independent layers.** The two layers are fit and simulated separately and simply handed to `netsim()` together. If you have built a single-layer network model before, you already know almost everything here.
#' - **Part 2: dependent layers.** We add the one realistic complication that makes multi-layer models genuinely different: a person's activity in one layer constrains the other. This is where all of the new machinery lives.
#' 
#' ::: {.callout-note}
#' Download the R script to follow along [here](mod11-Tutorial.R). Part 2 sources a small helper file, [`mod11-helpers.R`](mod11-helpers.R), which holds the cross-layer plumbing.
#' :::
#' 
#' ::: {.callout-note .concept}
#' ## The one new idea: same people, two edge sets
#' 
#' A layer is just a network model of the kind you have already been fitting. A multi-layer network is nothing more than **two of them defined over the same nodes.** Everything you know about `netest`, `netdx`, and `netsim` carries over, layer by layer. The only genuinely new mechanic is that you hand `netsim()` a *list* of fitted models instead of one, and it runs an epidemic that can transmit on either layer's edges.
#' 
#' Part 1 does exactly that, and nothing more. Part 2 adds the single feature that couples the layers, and everything hard about multi-layer modeling lives in that one feature.
#' :::
#' 
#' The picture below is the whole idea. On the left, one set of people carries two separate edge sets, a sparse long-lasting main layer and a denser short-lived casual layer, with the dotted lines marking the same person in each. On the right is the one complication Part 2 adds: the layers can interact, because a person busy in one layer tends to be sparse in the other.
#' 
## ----layerDiagram, echo = FALSE, fig.width = 9, fig.height = 4.4, fig.cap = "A multi-layer network is one node set with two edge sets. Left: the same eight people carry a main layer and a casual layer. Right: a hub in the main layer is often an isolate in the casual layer, which EpiModel represents by letting a node's degree in one layer enter the other layer's formation model."----
op <- par(mfrow = c(1, 2), mar = c(0.4, 0.4, 2.2, 0.4))

## A plane is a light rectangle; nodes sit at the same (x, y) within it, so
## stacking two planes shows one node set carrying two separate edge sets.
draw_plane <- function(xs, ys, H, edges, col_edge, lwd_edge, fill, label, hi = NULL) {
  rect(-0.03, H, 1.03, H + 1, border = "grey75", col = fill)
  for (i in seq_len(nrow(edges))) {
    a <- edges[i, 1]; b <- edges[i, 2]
    segments(xs[a], ys[a] + H, xs[b], ys[b] + H, col = col_edge, lwd = lwd_edge)
  }
  bg <- rep("white", length(xs)); if (!is.null(hi)) bg[hi] <- "#e74c3c"
  points(xs, ys + H, pch = 21, bg = bg, col = "grey30", cex = 1.5)
  text(0, H + 0.93, label, adj = c(0, 1), font = 2, cex = 0.85, col = col_edge)
}
Hm <- 2.2; Hc <- 0

## Panel A: the same eight people carry two edge sets
ax <- c(0.15, 0.50, 0.85, 0.22, 0.58, 0.90, 0.35, 0.75)
ay <- c(0.20, 0.12, 0.28, 0.62, 0.70, 0.60, 0.92, 0.95)
main_e   <- rbind(c(1, 3), c(4, 6), c(7, 8), c(2, 5))
casual_e <- rbind(c(1, 2), c(2, 3), c(4, 5), c(5, 6), c(1, 4), c(3, 6), c(5, 7), c(6, 8))
plot(NA, xlim = c(-0.05, 1.05), ylim = c(-0.1, 3.35), axes = FALSE,
     xlab = "", ylab = "", main = "Same people, two edge sets")
segments(ax, ay + Hc, ax, ay + Hm, col = "grey85", lty = 3)
draw_plane(ax, ay, Hc, casual_e, "#c0392b", 1.3, "#fdecea", "casual layer")
draw_plane(ax, ay, Hm, main_e,   "#2c3e50", 2.6, "#eaf0f6", "main layer")

## Panel B: a hub in one layer tends to be an isolate in the other
bx <- c(0.50, 0.25, 0.50, 0.75)
by <- c(0.28, 0.80, 0.92, 0.80)
hub_e <- rbind(c(1, 2), c(1, 3), c(1, 4))
cas_e <- rbind(c(2, 3), c(3, 4))
plot(NA, xlim = c(-0.05, 1.05), ylim = c(-0.1, 3.35), axes = FALSE,
     xlab = "", ylab = "", main = "Why the layers interact")
segments(bx, by + Hc, bx, by + Hm, col = "grey85", lty = 3)
draw_plane(bx, by, Hc, cas_e, "#c0392b", 1.3, "#fdecea", "casual layer", hi = 1)
draw_plane(bx, by, Hm, hub_e, "#2c3e50", 2.6, "#eaf0f6", "main layer", hi = 1)
par(op)

#' 
#' ## Setup
#' 
## ----setup, message = FALSE, warning = FALSE----------------------------------
library(EpiModel)
set.seed(12345)

# Simulation sizes, reused throughout.
nsims  <- 5
ncores <- 5
nsteps <- 250

#' 
#' One shared node set carries both layers. We attach a binary `race` attribute for homophily in the main layer.
#' 
## ----nwInit-------------------------------------------------------------------
n <- 500
nw <- network_initialize(n)
nw <- set_vertex_attribute(nw, "race", rep(0:1, length.out = n))

#' 
#' ## Part 1: Independent Layers
#' 
#' Nothing here is new. A layer is just a network model; the only new idea is keeping two of them for the same people.
#' 
#' We fit two ordinary single-layer TERGMs:
#' 
#' - **Layer 1 (main):** steady, long-lasting ties (mean duration 200) with race homophily.
#' - **Layer 2 (casual):** short, high-turnover ties (mean duration 20), driven by a degree-1 target.
#' 
#' Layer 2's formation references only its own structure (`degree(1)`), not layer 1. That independence is exactly what lets us skip all of the cross-layer machinery until Part 2.
#' 
## ----p1est, results = "hide", message = FALSE, warning = FALSE----------------
est1 <- netest(nw, ~edges + nodematch("race"), target.stats = c(90, 60),
               coef.diss = dissolution_coefs(~offset(edges), duration = 200))

est2 <- netest(nw, ~edges + degree(1), target.stats = c(75, 120),
               coef.diss = dissolution_coefs(~offset(edges), duration = 20))

#' 
#' Each layer is diagnosed on its own, exactly as any single-layer model.
#' 
## ----p1dx, results = "hide", message = FALSE, warning = FALSE-----------------
dx1 <- netdx(est1, nsims = nsims, ncores = ncores, nsteps = nsteps, verbose = FALSE)
dx2 <- netdx(est2, nsims = nsims, ncores = ncores, nsteps = nsteps, verbose = FALSE)

#' 
## ----p1dxPlot-----------------------------------------------------------------
par(mfrow = c(1, 2), mar = c(3, 3, 2, 1), mgp = c(2, 1, 0))
plot(dx1, main = "Layer 1 (main)")
plot(dx2, main = "Layer 2 (casual)")

#' 
#' Now the epidemic. The **one** multi-layer step is `netsim(list(est1, est2), ...)`: passing a *list* of fitted models. The list order sets the layer numbering (`est1` is layer 1). `resimulate.network = TRUE` redraws both dynamic layers each step, and the optional `multilayer()` wrapper lets us monitor a separate set of network statistics for each layer.
#' 
## ----p1sim, results = "hide", message = FALSE, warning = FALSE----------------
param   <- param.net(inf.prob = 0.5, act.rate = 2)
init    <- init.net(i.num = 10)
control <- control.net(type = "SI", nsteps = nsteps, nsims = nsims, ncores = ncores,
                       tergmLite = TRUE, resimulate.network = TRUE,
                       nwstats.formula = multilayer(~edges + nodematch("race"),
                                                    ~edges + degree(0:3)),
                       verbose = FALSE)

sim.ind <- netsim(list(est1, est2), param, init, control)

#' 
#' And we plot the epidemic, which is the whole point, and which single-layer intuition already prepares you to read. We do not set `mean.col`: with more than one outcome, `plot.netsim` gives each its own coordinated mean-and-band color.
#' 
## ----p1plot-------------------------------------------------------------------
sim.ind <- mutate_epi(sim.ind, prev = i.num / num)
par(mfrow = c(1, 2), mar = c(3, 3, 2, 1), mgp = c(2, 1, 0))
plot(sim.ind, y = c("s.num", "i.num"), legend = TRUE, main = "State sizes (independent)")
plot(sim.ind, y = "prev", main = "Prevalence")

#' 
#' The short-duration casual layer drives early spread, while the long-duration main layer sustains transmission over time. Because either layer can transmit, the epidemic reflects both.
#' 
#' Because we passed a `multilayer()` `nwstats.formula`, `netsim` also monitored each layer's formation statistics over the run. Plotting them is the network-side payoff: it confirms that both dynamic layers held their target structure while the epidemic ran on top of them.
#' 
## ----p1netplot----------------------------------------------------------------
par(mfrow = c(1, 1), mar = c(3, 3, 2, 1), mgp = c(2, 1, 0))
plot(sim.ind, type = "formation")

#' 
#' ## Part 2: Dependent Layers
#' 
#' Part 1's layers ignore each other, which allows something unrealistic: a person can be a hub in *both* the main and the casual layer at once. In reality people have finite relational capacity, so being busy in one partnership type tends to reduce activity in the other. That produces a **negative cross-layer degree correlation**, and reproducing it is the single decision that separates a dependent multi-layer model from an independent one.
#' 
#' ::: {.callout-note .concept}
#' ## Coupling the layers creates two problems to solve
#' 
#' We couple the layers by adding, to each layer's formation model, a term that reads the *other* layer's degree: `nodefactor("deg1+.net2")` in layer 1 and `nodefactor("deg1+.net1")` in layer 2, each with a small target so that few nodes are active in both. That one-line addition creates two chicken-and-egg problems, and naming them separately is the key to not getting lost:
#' 
#' 1. **A cold start.** Layer 1's model wants each node's layer-2 degree, but layer 2 does not exist yet, and vice versa. We have to seed plausible starting degrees *before* fitting anything.
#' 2. **Staying current.** During the simulation each layer's degree changes every step, so the other layer's model must be shown the updated value each time, or the coupling silently freezes at its starting values.
#' 
#' We solve the first with a `san()` bootstrap, which is one-time setup and lives in [`mod11-helpers.R`](mod11-helpers.R), and the second with a short per-step callback that we write out below.
#' :::
#' 
## ----p2source-----------------------------------------------------------------
source("mod11-helpers.R")

#' 
#' ### The cold-start bootstrap
#' 
#' `init_cross_degree()` runs the `san()` (simulated annealing) seeding described above and returns the network with the two cross-layer attributes attached. Open the helper file to see exactly what it does; here we just call it.
#' 
## ----p2boot, results = "hide", message = FALSE, warning = FALSE---------------
nw <- init_cross_degree(nw)

#' 
#' ### Estimate the dependent layers
#' 
#' The formation models are Part 1's plus one cross-layer `nodefactor` term each, and the targets gain one small statistic (10 of the 90 or 75 edges) that forces the negative correlation. The durations are unchanged, so the *only* difference from Part 1 is the coupling.
#' 
## ----p2est, results = "hide", message = FALSE, warning = FALSE----------------
est.1 <- netest(nw, ~edges + nodematch("race") + nodefactor("deg1+.net2"),
                target.stats = c(90, 60, 10),
                coef.diss = dissolution_coefs(~offset(edges), duration = 200))

est.2 <- netest(nw, ~edges + degree(1) + nodefactor("deg1+.net1"),
                target.stats = c(75, 120, 10),
                coef.diss = dissolution_coefs(~offset(edges), duration = 20))

#' 
#' Read the fitted models with `summary()`, the same tool you have used on network models since earlier in the course. The cross-layer `nodefactor` row is the one that matters here:
#' 
## ----p2coef-------------------------------------------------------------------
summary(est.1)
summary(est.2)

#' 
#' The cross-layer coefficient (`nodefactor.deg1+.net2` in layer 1, `nodefactor.deg1+.net1` in layer 2) is negative in both: a partner in the other layer *suppresses* partnering here. That negative sign is the whole point of Part 2.
#' 
#' ### Keep the cross-layer degree current during the simulation
#' 
#' Each step, `netsim` resimulates the two layers in turn, and each layer's formation model reads the *other* layer's current degree. If that degree is not refreshed between the resimulations it freezes at its starting value and the coupling silently switches off. The fix is a short callback that `netsim` calls at set points each step (flagged by the `network` argument), pulling each layer's current edgelist with `get_edgelist()` and recomputing the "has a partner in the other layer" attribute:
#' 
## ----p2callback---------------------------------------------------------------
network_layer_updates <- function(dat, at, network) {
  if (network == 0) {          # before layer 1 is redrawn: refresh from layer 2's edges
    dat <- set_attr(dat, "deg1+.net2", pmin(get_degree(get_edgelist(dat, network = 2)), 1))
  } else if (network == 1) {   # before layer 2 is redrawn: refresh from layer 1's edges
    dat <- set_attr(dat, "deg1+.net1", pmin(get_degree(get_edgelist(dat, network = 1)), 1))
  } else if (network == 2) {   # after layer 2 is redrawn: refresh for the summary stats
    dat <- set_attr(dat, "deg1+.net2", pmin(get_degree(get_edgelist(dat, network = 2)), 1))
  }
  dat
}

#' 
#' ### Simulate with the cross-layer callback
#' 
#' The control settings are Part 1's plus `dat.updates = network_layer_updates`, which wires the callback into the run.
#' 
## ----p2sim, results = "hide", message = FALSE, warning = FALSE----------------
control <- control.net(type = "SI", nsteps = nsteps, nsims = nsims, ncores = ncores,
                       tergmLite = TRUE, resimulate.network = TRUE,
                       dat.updates = network_layer_updates,
                       nwstats.formula = multilayer("formation",
                         ~edges + nodefactor("deg1+.net1") + degree(0:4)),
                       verbose = FALSE)

sim.dep <- netsim(list(est.1, est.2), param, init, control)
sim.dep <- mutate_epi(sim.dep, prev = i.num / num)

#' 
#' ## The Payoff: Does the Coupling Matter?
#' 
#' We built the dependency to be more realistic, but does it change the epidemic? Overlay the two prevalence curves.
#' 
## ----overlay------------------------------------------------------------------
par(mfrow = c(1, 1), mar = c(3, 3, 2, 1), mgp = c(2, 1, 0))
plot(sim.ind, y = "prev", mean.col = "steelblue", mean.lwd = 3, qnts = FALSE,
     ylim = c(0, 1), main = "Prevalence: independent vs dependent layers",
     ylab = "Prevalence", xlab = "Time step")
plot(sim.dep, y = "prev", mean.col = "firebrick", mean.lwd = 3, qnts = FALSE,
     add = TRUE)
legend("bottomright", c("Independent layers", "Dependent layers"), lwd = 3,
       col = c("steelblue", "firebrick"), bty = "n")

#' 
#' Under dependence the high-activity nodes in the two layers no longer coincide, so there are fewer people who bridge both partnership types and carry the pathogen between them. The plumbing of Part 2 is therefore not decoration: it changes which nodes do the transmitting, and so it changes the epidemic. Tying the network mechanics back to an epidemiological difference like this is the reason to bother with the dependent model at all.
#' 
#' ## Extension: Which Layer Drives Transmission?
#' 
#' A natural research question for a multi-layer model is what *share* of infections happens over each layer: are most new cases coming from the long main partnerships or the churning casual ones? The built-in `SI` module does not record this, because it transmits over both layers at once and reports a single `si.flow`. To split the count we replace the transmission module with a small custom one that walks the layers separately and attributes each new infection to the layer whose edge reached it.
#' 
#' This is the extension API from Module 9 applied to the infection step: read the state, do the process (here, once per layer), write it back, and record two new summary statistics with `set_epi()`.
#' 
## ----blmod--------------------------------------------------------------------
infect_bylayer <- function(dat, at) {
  status    <- get_attr(dat, "status")
  infTime   <- get_attr(dat, "infTime")
  inf.prob  <- get_param(dat, "inf.prob")
  act.rate  <- get_param(dat, "act.rate")
  finalProb <- 1 - (1 - inf.prob)^act.rate

  already <- integer(0)                 # infected this step, credited to the first layer that reached them
  count   <- c(main = 0, casual = 0)

  for (k in 1:2) {
    del <- discord_edgelist(dat, at, network = k)      # susceptible-infected pairs on layer k
    if (is.null(del)) next
    del <- del[!(del$sus %in% already), , drop = FALSE]
    if (nrow(del) == 0) next
    newly <- setdiff(unique(del$sus[rbinom(nrow(del), 1, finalProb) == 1]), already)
    if (length(newly) > 0) {
      status[newly]  <- "i"
      infTime[newly] <- at
      already  <- c(already, newly)
      count[k] <- length(newly)
    }
  }

  dat <- set_attr(dat, "status", status)
  dat <- set_attr(dat, "infTime", infTime)
  dat <- set_epi(dat, "si.flow.main",   at, count[1])
  dat <- set_epi(dat, "si.flow.casual", at, count[2])
  dat
}

#' 
#' We run the independent model again, but pass `type = NULL` and our module in the `infection.FUN` slot so `netsim` uses ours instead of the built-in transmission:
#' 
## ----blsim, results = "hide", message = FALSE, warning = FALSE----------------
control.bl <- control.net(type = NULL, nsteps = nsteps, nsims = nsims, ncores = ncores,
                          tergmLite = TRUE, resimulate.network = TRUE,
                          infection.FUN = infect_bylayer, verbose = FALSE)
sim.bl <- netsim(list(est1, est2), param, init, control.bl)

#' 
#' Now the by-layer flows are in the output. Plot them and read off the overall share:
#' 
## ----blplot-------------------------------------------------------------------
par(mfrow = c(1, 1), mar = c(3, 3, 2, 1), mgp = c(2, 1, 0))
plot(sim.bl, y = c("si.flow.main", "si.flow.casual"), legend = TRUE,
     main = "New infections per step, by transmission layer")

d   <- as.data.frame(sim.bl)
tot <- c(main   = sum(d$si.flow.main,   na.rm = TRUE),
         casual = sum(d$si.flow.casual, na.rm = TRUE))
round(tot / sum(tot), 2)

#' 
#' The printed shares tell you where transmission concentrated. Which layer dominates depends on the interplay of mean degree and duration: a denser, faster-churning layer reaches more distinct susceptibles, while a long-lasting layer keeps re-exposing the same partners. The point is the recipe, not the specific split: once each infection is attributed to a layer, "what fraction of transmission happens where" is a one-line `set_epi()` away.
#' 
#' ## When Do I Need Which?
#' 
#' The whole module reduces to one decision and a short checklist.
#' 
#' | | Independent layers | Dependent layers |
#' |---|---|---|
#' | **Use when** | activity in one layer does not constrain the other | finite relational capacity couples the layers |
#' | **Formation** | each layer references only its own structure | add a `nodefactor("deg1+.netX")` cross-term per layer |
#' | **Before fitting** | nothing extra | `san()` bootstrap to seed the cross-layer degree (`init_cross_degree()`) |
#' | **During simulation** | nothing extra | per-step callback to refresh it (`dat.updates = network_layer_updates`) |
#' | **The multi-layer step** | `netsim(list(est1, est2), ...)` | the same: `netsim(list(est.1, est.2), ...)` |
#' 
#' One-line rule: **independent = a list of `netest` calls; dependent = add a cross-layer term, a `san()` seed, and a `dat.updates` callback.**
#' 
#' ## Key EpiModel Functions for Multi-Layer Models
#' 
#' | Function / argument | Purpose |
#' |---|---|
#' | `netsim(list(est1, est2), ...)` | a *list* of fitted models is what makes a model multi-layer; list order sets the layer number |
#' | `resimulate.network = TRUE` | redraw both dynamic layers each step |
#' | `tergmLite = TRUE` | lightweight edgelist representation; required for the current-edgelist access (`get_edgelist(dat, network = k)`) the callback uses |
#' | `nwstats.formula = multilayer(f1, f2)` | monitor a separate set of network statistics for each layer |
#' | `dat.updates = <fn>` | per-step callback (dependent layers only) that refreshes the cross-layer attributes |
#' 
#' ## Your Turn
#' 
#' Put these models to work in the [Multi-Layer Networks Lab](mod11-Lab.html). You will re-time a layer to show that duration, not just formation, drives the epidemic; predict and then check whether the cross-layer coupling changes peak prevalence; use the by-layer tracker above to ask *where* transmission happens under independence versus dependence; and add a third layer for the price of one more `netest`. A worked solution accompanies the lab.