#' # Epidemics over Observed Networks
#' 
#' Everything in the course so far has followed one workflow: observe a network egocentrically (a sample of nodes and their partnerships), fit a temporal ERGM with `netest`, and simulate an epidemic from that fitted model. The fit is what lets a small sample stand in for a whole population, and what lets the network keep evolving during the run.
#' 
#' Sometimes you do not need the fit. If you have a **dynamic network census**, every node and every partnership directly observed over a series of time steps, you already have the object `netsim` would otherwise have to generate. This chapter shows the alternative workflow: place the observed `networkDynamic` object straight into the simulation and run the epidemic over it, with no ERGM estimation at all.
#' 
#' ::: {.callout-note .concept}
#' ## The one idea: skip the fit when the network is fully observed
#' 
#' `netest` exists to turn partial, egocentric data into a generative model you can simulate from. A full dynamic census is not partial, so there is nothing to estimate. You hand `netsim` the observed network, tell it not to regenerate the network each step (`resimulate.network = FALSE`), and supply two small custom modules so EpiModel reads from the observed object instead of from a model fit. Everything else, the transmission process, the trackers, the plots, is the EpiModel you already know.
#' :::
#' 
#' This is the one workflow in the course that bypasses the ERGM pipeline entirely. It is the right tool when the network really is observed in full over time (some contact-tracing, sensor, or animal-tracking datasets are), and the wrong tool when you have a sample and want to generalize from it, which is the far more common case and the reason the rest of the course teaches estimation.
#' 
#' ::: {.callout-note}
#' Download the R script to follow along [here](mod11-ObservedNets.R).
#' :::
#' 
#' ## Setup
#' 
#' We use a simulated dynamic network from the `networkDynamicData` package and treat it as if it were an observed census. The network is 1000 nodes with edge spells (onset and terminus times for each partnership) recorded over roughly 100 discrete time steps.
#' 
## ----setup, message = FALSE, warning = FALSE----------------------------------
library(EpiModel)
library(networkDynamicData)
set.seed(12345)

nsims  <- 5
ncores <- 5
nsteps <- 100

#' 
#' ## Two Custom Modules
#' 
#' Because there is no model fit, EpiModel's built-in initialization and infection modules, which both expect a `netest` object, do not apply. We replace them with two functions that work from the observed network directly. Both follow the read-process-write module contract from Module 9.
#' 
#' ### Initialization
#' 
#' `init_obsnw` replaces the built-in `initialize.net`. It builds the master data object, drops the observed `networkDynamic` straight into it, seeds the initial infections, and (optionally) stores disease status on the network as a temporally extended attribute so we can color network plots by status later.
#' 
## ----modInit------------------------------------------------------------------
init_obsnw <- function(x, param, init, control, s) {

  dat <- create_dat_object(param, init, control)

  dat$num.nw <- 1L
  dat$run$nw[[1]] <- x                              # the observed network goes in as-is
  dat <- set_param(dat, "groups", 1)

  i.num <- get_init(dat, "i.num")
  n <- network.size(dat$run$nw[[1]])
  dat <- append_core_attr(dat, 1, n)

  status <- rep("s", n)
  status[sample(1:n, i.num)] <- "i"                 # seed i.num random infections
  dat <- set_attr(dat, "status", status)

  infTime <- rep(NA, n)
  infTime[which(status == "i")] <- 1
  dat <- set_attr(dat, "infTime", infTime)

  track.nw.attr <- get_param(dat, "track.nw.attr", override.null.error = TRUE)
  if (!is.null(track.nw.attr) && track.nw.attr) {
    dat$run$nw[[1]] <- networkDynamic::activate.vertex.attribute(
      dat$run$nw[[1]], prefix = "testatus",
      value = get_attr(dat, "status"), onset = 1, terminus = Inf
    )
  }

  dat <- prevalence.net(dat, 1)
  return(dat)
}

#' 
#' ### Infection
#' 
#' `infect_obsnw` handles transmission. It reads the discordant (susceptible-infectious) partnerships active at the current step off the observed network with `discord_edgelist()`, exactly as a normal infection module would, and draws transmission. It supports two modes: a constant per-act probability (`inf.prob`), or a two-stage, duration-dependent probability (`inf.prob.stage1` / `inf.prob.stage2`), selected automatically by which parameters you supplied.
#' 
## ----modInfect----------------------------------------------------------------
infect_obsnw <- function(dat, at) {

  active  <- get_attr(dat, "active")
  status  <- get_attr(dat, "status")
  infTime <- get_attr(dat, "infTime")

  act.rate <- get_param(dat, "act.rate")
  inf.prob.stage1 <- get_param(dat, "inf.prob.stage1", override.null.error = TRUE)

  time_varying <- !is.null(inf.prob.stage1)
  if (time_varying) {
    inf.prob.stage2 <- get_param(dat, "inf.prob.stage2")
    dur.stage1 <- get_param(dat, "dur.stage1")
  } else {
    inf.prob <- get_param(dat, "inf.prob")
  }

  idsInf <- which(active == 1 & status == "i")
  nActive <- sum(active == 1)
  nElig <- length(idsInf)

  totInf <- 0

  if (nElig > 0 && nElig < nActive) {

    del <- discord_edgelist(dat, at)                # SI pairs active on the observed network

    if (!is.null(del)) {

      if (time_varying) {
        infDur.del <- at - infTime[del$inf]
        del$transProb <- ifelse(infDur.del <= dur.stage1,
                                inf.prob.stage1, inf.prob.stage2)
      } else {
        del$transProb <- inf.prob
      }

      del$actRate <- act.rate
      del$finalProb <- 1 - (1 - del$transProb) ^ del$actRate
      transmit <- rbinom(nrow(del), 1, del$finalProb)
      del <- del[which(transmit == 1), ]

      idsNewInf <- unique(del$sus)
      totInf <- length(idsNewInf)

      if (totInf > 0) {
        status[idsNewInf] <- "i"
        infTime[idsNewInf] <- at
        dat <- set_attr(dat, "status", status)
        dat <- set_attr(dat, "infTime", infTime)

        track.nw.attr <- get_param(dat, "track.nw.attr", override.null.error = TRUE)
        if (!is.null(track.nw.attr) && track.nw.attr) {
          dat$run$nw[[1]] <- networkDynamic::activate.vertex.attribute(
            dat$run$nw[[1]], prefix = "testatus",
            value = "i", onset = at, terminus = Inf, v = idsNewInf
          )
        }
      }
    }
  }

  dat <- set_epi(dat, "si.flow", at, totInf)
  return(dat)
}

#' 
#' ## Loading the Observed Network
#' 
#' The `concurrencyComparisonNets` dataset provides an observed dynamic network as its `base` object. We remove its pre-recorded disease status attribute, since we simulate our own epidemic on top of the observed contact structure.
#' 
## ----loadNw-------------------------------------------------------------------
data(concurrencyComparisonNets)
nw <- base
nw <- network::delete.vertex.attribute(nw, "status.active")
print(nw)

#' 
#' ## Example 1: A Basic SI Epidemic
#' 
#' The first run is a plain SI model: a constant 50% per-act transmission probability and one act per partnership per step.
#' 
## ----ex1param-----------------------------------------------------------------
param <- param.net(inf.prob = 0.5, act.rate = 1)
init  <- init.net(i.num = 10)

#' 
#' The control settings are where the observed-network workflow shows itself. We pass our two custom modules, keep the built-in `prevalence.net` recorder, and set `resimulate.network = FALSE` so the observed edges are used as-is. There is no ERGM fit to compute network statistics from, so we turn off `save.nwstats`. Critically, `nsteps` matches the number of observed time steps in the network.
#' 
## ----ex1control---------------------------------------------------------------
control <- control.net(
  type = NULL,
  nsteps = nsteps,
  nsims = nsims,
  ncores = ncores,
  initialize.FUN = init_obsnw,
  infection.FUN = infect_obsnw,
  prevalence.FUN = prevalence.net,
  resimulate.network = FALSE,
  save.nwstats = FALSE,
  verbose = FALSE
)

#' 
#' Note that the first argument to `netsim` is the observed network `nw`, not a fitted `netest` object. That substitution is the whole point.
#' 
## ----ex1sim, results = "hide", message = FALSE, warning = FALSE---------------
sim <- netsim(nw, param, init, control)

#' 
## ----ex1print-----------------------------------------------------------------
sim

#' 
#' ## A Caveat: Match `nsteps` to the Observation Window
#' 
#' The observed network has edge activity recorded over a finite window of about 100 steps. By a `networkDynamic` convention, edges active at the last observed time stay active indefinitely, so nothing stops you from simulating past the window, but the results stop meaning anything: the edge set is frozen, so it no longer reflects real turnover.
#' 
#' You can see the frozen edges directly. The dyads active at step 100 and step 200 are identical, because nothing dissolves after the observation ends:
#' 
## ----caveatDyads--------------------------------------------------------------
head(get.dyads.active(nw, at = 1),   5)
head(get.dyads.active(nw, at = 100), 5)
head(get.dyads.active(nw, at = 200), 5)

#' 
#' Running 200 steps on the roughly 100-step network shows the symptom: prevalence plateaus and incidence falls to zero once the frozen structure takes over.
#' 
## ----caveatSim, results = "hide", message = FALSE, warning = FALSE------------
control_caveat <- control.net(
  type = NULL, nsteps = 200, nsims = nsims, ncores = ncores,
  initialize.FUN = init_obsnw, infection.FUN = infect_obsnw,
  prevalence.FUN = prevalence.net,
  resimulate.network = FALSE, save.nwstats = FALSE, verbose = FALSE
)
sim_caveat <- netsim(nw, param, init, control_caveat)

#' 
## ----caveatPlot---------------------------------------------------------------
par(mfrow = c(1, 2), mar = c(3, 3, 2, 1), mgp = c(2, 1, 0))
plot(sim_caveat, y = "i.num", main = "Prevalence past the window",
     ylab = "Infected count", xlab = "Time step", legend = FALSE)
plot(sim_caveat, y = "si.flow", main = "Incidence past the window",
     ylab = "New infections", xlab = "Time step", legend = FALSE)

#' 
#' The rule is simple: **match `nsteps` to the number of observed time steps.**
#' 
#' ## Example 2: Time-Varying Transmission
#' 
#' The second run uses the duration-dependent mode of the infection module. Transmission is low in the primary stage (5% per act for the first 5 steps after infection) and higher afterward (15% per act), a common shape for infections with an early low-shedding period. Setting `track.nw.attr = TRUE` stores disease status on the network so we can visualize it.
#' 
## ----ex2param-----------------------------------------------------------------
param <- param.net(
  inf.prob.stage1 = 0.05,
  inf.prob.stage2 = 0.15,
  dur.stage1 = 5,
  act.rate = 1,
  track.nw.attr = TRUE
)
init <- init.net(i.num = 10)

control <- control.net(
  type = NULL, nsteps = nsteps, nsims = nsims, ncores = ncores,
  initialize.FUN = init_obsnw, infection.FUN = infect_obsnw,
  prevalence.FUN = prevalence.net,
  resimulate.network = FALSE, save.nwstats = FALSE, verbose = FALSE
)

#' 
## ----ex2sim, results = "hide", message = FALSE, warning = FALSE---------------
sim2 <- netsim(nw, param, init, control)

#' 
#' Because we tracked status on the network, we can plot the observed network colored by disease status at an early and a late step.
#' 
## ----ex2net-------------------------------------------------------------------
par(mfrow = c(1, 2), mar = c(1, 1, 1, 1))
plot(sim2, type = "network", col.status = TRUE, at = 2,      sims = 1)
plot(sim2, type = "network", col.status = TRUE, at = nsteps, sims = 1)

#' 
#' The same stored attribute can be read back at the individual level with the `networkDynamic` API, which is how you would extract who was infected when for a downstream analysis:
#' 
## ----ex2attr------------------------------------------------------------------
nwd <- get_network(sim2, 1)
head(get.vertex.attribute.active(nwd, "testatus", at = nsteps), 10)

#' 
#' ## Comparing the Two Runs
#' 
#' The constant, high transmission of Example 1 saturates quickly; the lower, staged transmission of Example 2 climbs more slowly.
#' 
## ----cmpPlot------------------------------------------------------------------
sim  <- mutate_epi(sim,  prev = i.num / num)
sim2 <- mutate_epi(sim2, prev = i.num / num)
par(mfrow = c(1, 2), mar = c(3, 3, 2, 1), mgp = c(2, 1, 0))
plot(sim,  y = "i.num", main = "Ex 1: constant transmission",
     ylab = "Infected count", xlab = "Time step", legend = FALSE)
plot(sim2, y = "i.num", main = "Ex 2: time-varying transmission",
     ylab = "Infected count", xlab = "Time step", legend = FALSE)

#' 
#' ## When to Use This
#' 
#' | | Fit and simulate (the rest of the course) | Observed network (this chapter) |
#' |---|---|---|
#' | **You have** | an egocentric sample of the network | a full dynamic census of nodes and edges over time |
#' | **Network step** | `netest` fits a TERGM you simulate from | the observed `networkDynamic` object is used directly |
#' | **`resimulate.network`** | `TRUE` (the model regenerates the network) | `FALSE` (the observed edges are fixed) |
#' | **Custom modules** | usually not needed for the network itself | a custom `initialize.FUN` and `infection.FUN` |
#' | **Run length** | any `nsteps` the model supports | `nsteps` capped at the observation window |
#' 
#' The estimation workflow is the general-purpose one, because full network censuses are rare and samples are the norm. Reach for the observed-network approach only when you genuinely have the whole network over time and want to run transmission on exactly that structure without a modeling layer in between.
#' 
#' ## Your Turn
#' 
#' - Convert the model to SIS or SIR by adding a recovery module (the Module 9 pattern), and see how the fixed observation window interacts with recovery.
#' - Swap in a different observed dynamic network from `networkDynamicData` and confirm you reset `nsteps` to its own observation length.
#' - Bridge the two workflows: take a *static* observed network, assume an edge dissolution rate, and fit a TERGM with `netest`, the approach from the Module 6 network-census lab, then contrast it with using a dynamic census directly here.