tests/testthat/test-newmodules.R

In EpiModel: Mathematical Modeling of Infectious Disease Dynamics

context("New Network Models")

test_that("New network models vignette example", {
  skip_on_cran()

  ## New Aging Module
  aging <- function(dat, at) {

    age <- get_attr(dat, "age", override.null.error = TRUE)
    if (is.null(age)) {
      active <- get_attr(dat, "active")
      n <- sum(active == 1)
      age <- sample(18:49, n, replace = TRUE)
    } else {
      age <- get_attr(dat, "age") + 1 / 12
    }
    dat <- set_attr(dat, "age", age)

    return(dat)
  }


  ## Replacement Departure Module
  dfunc <- function(dat, at) {
    active <- get_attr(dat, "active")
    exitTime <- get_attr(dat, "exitTime")
    idsElig <- which(active == 1)
    nElig <- length(idsElig)

    nDepartures <- 0

    if (nElig > 0) {
      ages <- get_attr(dat, "age")[idsElig]
      life.expt <- get_param(dat, "life.expt")
      departure.rates <- pmin(1, 1 / (life.expt * 12 - ages * 12))
      vecDepartures <- which(rbinom(nElig, 1, departure.rates) == 1)
      idsDepartures <- idsElig[vecDepartures]
      nDepartures <- length(idsDepartures)
      if (nDepartures > 0) {
        active[idsDepartures] <- 0
        exitTime[idsDepartures] <- at
        dat <- set_attr(dat, "active", active)
        dat <- set_attr(dat, "exitTime", exitTime)
      }
    }

    # Output
    dat <- set_epi(dat, "d.flow", nDepartures)
    return(dat)
  }


  ## Replacement Arrival Module
  afunc <- function(dat, at) {

    # Variables
    growth.rate <- get_param(dat, "growth.rate")
    exptPopSize <- get_epi(dat, "num", 1) * (1 + growth.rate * at)
    active <- get_attr(dat, "active")
    numNeeded <- exptPopSize - sum(active == 1)

    if (numNeeded > 0) {
      nArrivals <- rpois(1, numNeeded)
    } else {
      nArrivals <- 0
    }

    dat <- append_core_attr(dat, at, nArrivals)
    dat <- append_attr(dat, "status", "s", nArrivals)
    dat <- append_attr(dat, "infTime", NA, nArrivals)
    dat <- append_attr(dat, "age", 0, nArrivals)

    # Output
    dat <- set_epi(dat, "a.flow", nArrivals)

    return(dat)
  }


  ## Network Model
  nw <- network.initialize(50, directed = FALSE)
  est <- netest(nw, formation = ~edges, target.stats = 15,
                coef.diss = dissolution_coefs(~offset(edges), 60, 0.000274),
                verbose = FALSE)


  ## EpiModel
  param <- param.net(inf.prob = 0.35, growth.rate = 0.00083, life.expt = 70)
  init <- init.net(i.num = 10)
  control <- control.net(type = NULL, nsims = 1, nsteps = 5,
                         departures.FUN = dfunc,
                         arrivals.FUN = afunc, aging.FUN = aging,
                         infection.FUN = infection.net,
                         tergmLite = FALSE, resimulate.network = TRUE, verbose = FALSE)
  mod1 <- netsim(est, param, init, control)
  capture_output(
    mod1
  )

  expect_is(mod1, "netsim")
  expect_output(print(mod1), "resim_nets.FUN")
  expect_output(print(mod1), "infection.FUN")
  expect_output(print(mod1), "departures.FUN")
  expect_output(print(mod1), "arrivals.FUN")
  expect_output(print(mod1), "aging.FUN")

  ## Test module reordering
  control <- control.net(type = NULL, nsims = 1, nsteps = 10,
                         departures.FUN = dfunc,
                         arrivals.FUN = afunc, aging.FUN = aging,
                         infection.FUN = infection.net,
                         module.order = c("resim_nets.FUN", "summary_nets.FUN",
                                          "infection.FUN", "aging.FUN",
                                          "arrivals.FUN", "departures.FUN",
                                          "nwupdate.FUN", "prevalence.FUN"),
                         tergmLite = FALSE, resimulate.network = TRUE, verbose = FALSE)
  mod2 <- netsim(est, param, init, control)
  expect_is(mod2, "netsim")

  ### tergmLite replication
  param <- param.net(inf.prob = 0.35, growth.rate = 0.00083, life.expt = 70)
  init <- init.net(i.num = 10)
  control <- control.net(type = NULL, nsims = 1, nsteps = 10,
                         infection.FUN = infection.net,
                         departures.FUN = dfunc,
                         arrivals.FUN = afunc, aging.FUN = aging,
                         tergmLite = TRUE, verbose = FALSE,
                         resimulate.network = TRUE)
  mod3 <- netsim(est, param, init, control)
  expect_is(mod3, "netsim")

  ## Test module reordering
  control <- control.net(type = NULL, nsims = 1, nsteps = 10,
                         departures.FUN = dfunc,
                         arrivals.FUN = afunc, aging.FUN = aging,
                         infection.FUN = infection.net,
                         module.order = c("resim_nets.FUN", "summary_nets.FUN",
                                          "infection.FUN", "aging.FUN",
                                          "arrivals.FUN", "departures.FUN",
                                          "nwupdate.FUN", "prevalence.FUN"),
                         tergmLite = TRUE, resimulate.network = TRUE, verbose = FALSE)
  mod4 <- netsim(est, param, init, control)
  expect_is(mod4, "netsim")

  ## "updated" infection module
  infect <- infection.net
  control <- control.net(type = NULL, nsims = 1, nsteps = 10,
                         departures.FUN = dfunc,
                         arrivals.FUN = afunc, aging.FUN = aging,
                         infection.FUN = infect,
                         module.order = c("resim_nets.FUN", "summary_nets.FUN",
                                          "infection.FUN", "aging.FUN",
                                          "arrivals.FUN", "departures.FUN",
                                          "nwupdate.FUN", "prevalence.FUN"),
                         tergmLite = TRUE, resimulate.network = TRUE, verbose = FALSE)
  mod5 <- netsim(est, param, init, control)
  expect_is(mod5, "netsim")

  expect_output(print(mod5), "resim_nets.FUN")
  expect_output(print(mod5), "infection.FUN")
  expect_output(print(mod5), "departures.FUN")
  expect_output(print(mod5), "arrivals.FUN")
  expect_output(print(mod5), "aging.FUN")

})

test_that("module.order is independent of resimulate.network", {
  skip_on_cran()

  nw <- network_initialize(n = 50)
  est <- netest(nw, formation = ~edges, target.stats = 15,
                coef.diss = dissolution_coefs(~offset(edges), 60),
                verbose = FALSE)
  param <- param.net(inf.prob = 0.3)
  init <- init.net(i.num = 5)

  custom_order <- c("resim_nets.FUN", "summary_nets.FUN", "infection.FUN",
                     "nwupdate.FUN", "prevalence.FUN")

  # module.order preserved with resimulate.network = TRUE
  ctrl1 <- control.net(type = "SI", nsims = 1, nsteps = 5,
                        module.order = custom_order,
                        resimulate.network = TRUE, verbose = FALSE)
  expect_equal(ctrl1$module.order, custom_order)
  mod1 <- netsim(est, param, init, ctrl1)
  expect_s3_class(mod1, "netsim")

  # module.order preserved with resimulate.network = FALSE
  ctrl2 <- control.net(type = "SI", nsims = 1, nsteps = 5,
                        module.order = custom_order,
                        resimulate.network = FALSE, verbose = FALSE)
  expect_equal(ctrl2$module.order, custom_order)
  mod2 <- netsim(est, param, init, ctrl2)
  expect_s3_class(mod2, "netsim")
})

test_that("tergmLite = TRUE warns when overriding resimulate.network", {
  expect_warning(
    control.net(type = "SI", nsteps = 10,
                tergmLite = TRUE, resimulate.network = FALSE),
    "resimulate.network"
  )

  # confirm the override took effect
  ctrl <- suppressWarnings(
    control.net(type = "SI", nsteps = 10,
                tergmLite = TRUE, resimulate.network = FALSE)
  )
  expect_true(ctrl$resimulate.network)

  # no warning when resimulate.network = TRUE (no conflict)
  expect_no_warning(
    control.net(type = "SI", nsteps = 10,
                tergmLite = TRUE, resimulate.network = TRUE)
  )
})

context("Network Model with Param Updater")

test_that("netsim with param updater", {
  skip_on_cran()
  # Create the list.param.updaters
  list.param.updaters <- list(
    # this is one updater
    list(
      at = 10,
      verbose = TRUE,
      param = list(
        inf.prob = 0.3,
        act.rate = 0.3
      )
    ),
    # this is another updater
    list(
      at = 20,
      verbose = TRUE,
      param = list(
        # inf.prob = function(x) plogis(qlogis(x) - log(10)),
        # act.rate = function(x) plogis(qlogis(x) - log(10))
        inf.prob = 0.01
      )
    )
  )

  # Create the list.control.updaters
  list.control.updaters <- list(
    # this is one updater
    list(
      at = 30,
      verbose = TRUE,
      control = list(
        resimulate.network = FALSE
      )
    )
  )

  # Do not forget to add it to `param`
  param <- param.net(
    inf.prob = 0.1,
    act.rate = 0.1,
    .param.updater.list = list.param.updaters
  )

  # Enable the module in `control`
  control <- control.net(
    type = NULL, # must be NULL as we use a custom module
    nsims = 1,
    nsteps = 50,
    verbose = FALSE,
    infection.FUN = infection.net,
    .control.updater.list = list.control.updaters,
    resimulate.network = TRUE
  )

  nw <- network_initialize(n = 50)
  nw <- set_vertex_attribute(nw, "race", rbinom(50, 1, 0.5))
  est <- netest(
    nw,
    formation = ~edges,
    target.stats = 25,
    coef.diss = dissolution_coefs(~offset(edges), 10, 0),
    verbose = FALSE
  )

  init <- init.net(i.num = 10)

  expect_message(mod <- netsim(est, param, init, control))

  # `resimulate.network` is turned of at step 30. We check that the number of
  # observations in the "networkDynamic" object is < than 31 and not 50 (the
  # number of timestep in the simulation)
  n_obs <- length(
    get.network.attribute(mod$network[[1]][[1]], 'net.obs.period')$observations
  )
  expect_lt(n_obs, 31)
})

test_that("update_list merges named lists and replaces lists without names", {
  old <- list(inf.prob = multilayer(0.3, 0.1), act.rate = 1,
              nested = list(a = 1, b = 2), unnamed = list(1, 2), x = 2)
  new <- list(inf.prob = multilayer(0.05, 0.1), act.rate = multilayer(1, 2),
              nested = list(b = 3), unnamed = list(5), x = function(v) v * 3)
  out <- update_list(old, new)
  expect_identical(out$inf.prob, multilayer(0.05, 0.1))
  expect_identical(out$act.rate, multilayer(1, 2))
  expect_identical(out$nested, list(a = 1, b = 3))
  expect_identical(out$unnamed, list(5))
  expect_equal(out$x, 6)
})

context("Network Model with Scenarios")

test_that("SIS with scenarios", {
  skip_on_cran()
  set.seed(10)

  nw <- network_initialize(n = 200)
  est <- netest(nw,
    formation = ~edges, target.stats = 60,
    coef.diss = dissolution_coefs(~offset(edges), 10, 0),
    verbose = FALSE
  )

  param <- param.net(inf.prob = 0.9, rec.rate = 0.01, act.rate = 2)
  control <- control.net(type = "SIS", nsims = 1, nsteps = 50, verbose = FALSE)
  init <- init.net(i.num = 10)

  scenarios.df <- dplyr::tribble(
    ~.scenario.id, ~.at, ~inf.prob, ~rec.rate,
    "base", 0, 0.9, 0.01,
    "multiple_changes", 0, 0.1, 0.04,
    "multiple_changes", 20, 0.9, 0.01,
    "multiple_changes", 40, 0.1, 0.1
  )

  scenarios.list <- create_scenario_list(scenarios.df)
  expect_length(scenarios.list, 2)

  sc.param <- use_scenario(param, scenarios.list[[1]])
  expect_silent(netsim(est, sc.param, init, control))
  sc.param <- use_scenario(param, scenarios.list[[2]])
  expect_message(netsim(est, sc.param, init, control))

  # .at not a integer
  scenarios.df <- dplyr::tribble(
    ~.scenario.id, ~.at, ~inf.prob, ~rec.rate,
    "multiple_changes", "text", 0.1, 0.1
  )
  expect_error(create_scenario_list(scenarios.df))

  # inf_prob with an underscore
  scenarios.df <- dplyr::tribble(
    ~.scenario.id, ~.at, ~inf_prob, ~rec.rate,
    "multiple_changes", 0, 0.1, 0.1
  )
  expect_error(scenarios.list <- create_scenario_list(scenarios.df))

  # rec.rate2 not in param
  scenarios.df <- dplyr::tribble(
    ~.scenario.id, ~.at, ~inf.prob, ~rec.rate2,
    "multiple_changes", 0, 0.1, 0.1
  )
  scenarios.list <- create_scenario_list(scenarios.df)
  expect_error(sc.param <- use_scenario(param, scenarios.list[[1]]))
})

test_that("create_scenario_list handles a single parameter column (#1045)", {
  # Selecting a row of a one-parameter table dropped it to an unnamed scalar,
  # so `unflatten_params` was handed NULL names and `strsplit` failed with
  # "non-character argument". No netsim here, so this runs on CRAN too.
  sc.df <- data.frame(
    .scenario.id = c("a", "b"),
    .at          = 0,
    inf.prob     = c(0.2, 0.3),
    stringsAsFactors = FALSE
  )

  sc <- create_scenario_list(sc.df)
  expect_length(sc, 2)
  expect_equal(names(sc), c("a", "b"))
  expect_equal(sc[["a"]][[".param.updater.list"]][[1]]$param,
               list(inf.prob = 0.2))
  expect_equal(sc[["b"]][[".param.updater.list"]][[1]]$param,
               list(inf.prob = 0.3))

  # The value reaches `param` through `use_scenario` as well.
  param <- param.net(inf.prob = 0.9, act.rate = 2)
  expect_equal(use_scenario(param, sc[["b"]])$inf.prob, 0.3)

  # Several `.at` rows for one scenario, still one parameter column.
  sc.df2 <- data.frame(
    .scenario.id = "ramp",
    .at          = c(0, 10, 20),
    inf.prob     = c(0.1, 0.2, 0.3),
    stringsAsFactors = FALSE
  )
  updaters <- create_scenario_list(sc.df2)[["ramp"]][[".param.updater.list"]]
  expect_length(updaters, 3)
  expect_equal(vapply(updaters, function(u) u$param$inf.prob, numeric(1)),
               c(0.1, 0.2, 0.3))

  # A lone column carrying a position suffix rebuilds a length-1 vector.
  sc.df3 <- data.frame(.scenario.id = "s", .at = 0, d.rate_1 = 0.5,
                       check.names = FALSE, stringsAsFactors = FALSE)
  expect_equal(
    create_scenario_list(sc.df3)[["s"]][[".param.updater.list"]][[1]]$param,
    list(d.rate = 0.5)
  )

  # Two parameter columns were never affected; kept as the contrast case.
  sc.df4 <- data.frame(.scenario.id = "a", .at = 0, inf.prob = 0.2,
                       act.rate = 1, stringsAsFactors = FALSE)
  expect_equal(create_scenario_list(sc.df4)[["a"]][[".param.updater.list"]][[1]]$param,
               list(inf.prob = 0.2, act.rate = 1))
})

test_that("a scenario sets single positions of a vector parameter", {
  # Setting `x_1` alone replaced `x` with a vector of length 1, and `x_2`
  # alone with `c(NA, value)`. No netsim here, so this runs on CRAN too.
  param <- param.net(inf.prob = 0.3, x = c(0.1, 0.2, 0.3))
  sc_param <- function(..., p = param) {
    sc.df <- data.frame(.scenario.id = "s", .at = 0, ..., check.names = FALSE)
    use_scenario(p, create_scenario_list(sc.df)[["s"]])
  }

  # The positions set change and the others keep their current value.
  expect_equal(sc_param(x_1 = 0.5)$x, c(0.5, 0.2, 0.3))
  expect_equal(sc_param(x_2 = 0.5)$x, c(0.1, 0.5, 0.3))
  expect_equal(sc_param(x_3 = 0.7, x_1 = 0.5)$x, c(0.5, 0.2, 0.7))
  expect_equal(sc_param(x_1 = 1, x_2 = 2, x_3 = 3)$x, c(1, 2, 3))
  expect_equal(sc_param(x_2 = 0.5)$inf.prob, 0.3)

  # A column without a position suffix replaces the whole parameter.
  expect_equal(sc_param(x = 0.5)$x, 0.5)

  # Positions past the end extend the vector when they leave no gap.
  expect_equal(sc_param(x_4 = 0.4)$x, c(0.1, 0.2, 0.3, 0.4))
  expect_equal(sc_param(x_5 = 0.5, x_4 = 0.4)$x, c(0.1, 0.2, 0.3, 0.4, 0.5))
  expect_error(sc_param(x_5 = 0.5), "no value for position 4")
  expect_error(sc_param(x_6 = 0.6, x_4 = 0.4), "no value for position 5")

  # Names are kept and other atomic types are merged the same way.
  p.chr <- param.net(grp = c(a = "low", b = "mid", c = "high"))
  expect_equal(sc_param(grp_2 = "none", p = p.chr)$grp,
               c(a = "low", b = "none", c = "high"))

  # Rows with .at < 2 are applied in turn, each to the result of the last.
  sc.df <- data.frame(.scenario.id = "s", .at = c(0, 1), x_1 = c(0.5, 0.6),
                      x_3 = c(0.7, 0.7), check.names = FALSE)
  expect_equal(use_scenario(param, create_scenario_list(sc.df)[["s"]])$x,
               c(0.6, 0.2, 0.7))

  # The scenario records the positions next to the values.
  sc <- create_scenario_list(
    data.frame(.scenario.id = "s", .at = 5, x_2 = 0.5, inf.prob = 0.2,
               check.names = FALSE)
  )[["s"]]
  updater <- sc[[".param.updater.list"]][[1]]
  expect_equal(updater$param, list(x = c(NA, 0.5), inf.prob = 0.2))
  expect_equal(updater$positions, list(x = 2))
  # Scalar columns alone add no positions.
  sc.scalar <- create_scenario_list(
    data.frame(.scenario.id = "s", .at = 5, inf.prob = 0.2)
  )[["s"]]
  expect_null(sc.scalar[[".param.updater.list"]][[1]]$positions)

  # Parameters that are not plain vectors cannot be set by position.
  p.ml <- param.net(inf.prob = 0.3, act.rate = multilayer(1, 2))
  expect_error(sc_param(act.rate_2 = 3, p = p.ml),
               "`act.rate` is a multilayer object rather than a plain vector")
  p.bid <- param.net(inf.prob = by_infection_duration(c(0.5, 0.1)))
  expect_error(sc_param(inf.prob_2 = 0.2, p = p.bid),
               "`inf.prob` is a by_infection_duration object")
})

test_that("a scenario updater sets single positions of a vector parameter", {
  skip_on_cran()
  set.seed(10)

  nw <- network_initialize(n = 50)
  est <- netest(nw,
                formation = ~edges,
                target.stats = 25,
                coef.diss = dissolution_coefs(~offset(edges), 10, 0),
                verbose = FALSE)

  # Records the value of the vector parameter `x` at each step
  record_x <- function(dat, at) {
    x <- get_param(dat, "x")
    dat <- set_epi(dat, "x.len", length(x))
    for (i in 1:3) dat <- set_epi(dat, paste0("x", i), x[i])
    dat
  }

  param <- param.net(inf.prob = 0.3, x = c(0.1, 0.2, 0.3))
  init <- init.net(i.num = 10)
  control <- control.net(type = NULL, nsims = 1, nsteps = 6,
                         infection.FUN = infection.net,
                         record.FUN = record_x, verbose = FALSE)

  # The updater sets the second position at step 4
  sc.df <- data.frame(.scenario.id = "s", .at = 4, x_2 = 0.7,
                      check.names = FALSE)
  sc.param <- use_scenario(param, create_scenario_list(sc.df)[["s"]])
  expect_equal(sc.param$x, c(0.1, 0.2, 0.3))

  sim <- suppressMessages(netsim(est, sc.param, init, control))
  d <- as.data.frame(sim)[2:6, ]
  expect_equal(d$x.len, rep(3, 5))
  expect_equal(d$x1, rep(0.1, 5))
  expect_equal(d$x2, c(0.2, 0.2, 0.7, 0.7, 0.7))
  expect_equal(d$x3, rep(0.3, 5))
  expect_equal(sim$param$x, c(0.1, 0.7, 0.3))

  # An updater position that would leave a gap stops the run
  sc.gap <- data.frame(.scenario.id = "s", .at = 3, x_5 = 0.5,
                       check.names = FALSE)
  gap.param <- use_scenario(param, create_scenario_list(sc.gap)[["s"]])
  expect_error(suppressMessages(netsim(est, gap.param, init, control)),
               "no value for position 4")
})

context("Records: attr_history and Raw Objects")

test_that("Time varying elements", {
  skip_on_cran()
  test_logger <- function(dat, at) {
    nodes <- get_posit_ids(dat)

    some_nodes <- sample(nodes, 5)
    dat <- record_attr_history(
      dat, "attr_norm",
      rnorm(length(some_nodes)),
      posit_ids = some_nodes
    )

    some_nodes <- sample(nodes, 5)
    dat <- record_attr_history(
      dat, "attr_unif",
      runif(length(some_nodes)),
      posit_ids = some_nodes
    )

    some_nodes <- sample(nodes, 5)
    dat <- record_attr_history(
      dat, "attr_fix",
      at,
      posit_ids = some_nodes
    )

    # test when 0 nodes selected
    some_nodes <- integer(0)
    dat <- record_attr_history(
      dat, "attr_none",
      at,
      posit_ids = some_nodes
    )

    return(dat)
  }

  param <- param.net(
    inf.prob = 0.1,
    act.rate = 0.1
  )

  # Enable the module in `control`
  control <- control.net(
    type = NULL, # must be NULL as we use a custom module
    nsims = 1,
    nsteps = 20,
    verbose = FALSE,
    infection.FUN = infection.net,
    logger.FUN = test_logger
  )

  nw <- network_initialize(n = 50)
  nw <- set_vertex_attribute(nw, "race", rbinom(50, 1, 0.5))
  est <- netest(
    nw,
    formation = ~edges,
    target.stats = 25,
    coef.diss = dissolution_coefs(~ offset(edges), 10, 0),
    verbose = FALSE
  )

  init <- init.net(i.num = 10)

  mod <- netsim(est, param, init, control)
  attr_history <- get_attr_history(mod)
  expect_is(attr_history, "list")
  expect_is(attr_history[[1]], "data.frame")
  expect_setequal(
    names(attr_history),
    c("attr_norm", "attr_unif", "attr_fix", "attr_none"))
  # `attr_fix` records the time step as its value, at the default `at`
  expect_equal(attr_history$attr_fix$values, attr_history$attr_fix$time)
  expect_setequal(attr_history$attr_fix$time, 2:20)
})

test_that("the module `at` matches the current step used by set_epi", {
  skip_on_cran()
  step_logger <- function(dat, at) {
    dat <- set_epi(dat, "at_arg", at)
    dat <- set_epi(dat, "flow_like", at * 10)
    return(dat)
  }
  nw <- network_initialize(n = 50)
  est <- netest(nw, formation = ~edges, target.stats = 20,
                coef.diss = dissolution_coefs(~offset(edges), 10, 0),
                verbose = FALSE)
  param <- param.net(inf.prob = 0.3, act.rate = 1)
  init <- init.net(i.num = 5)
  control <- control.net(type = NULL, nsteps = 5, nsims = 1, verbose = FALSE,
                         tergmLite = TRUE, resimulate.network = TRUE,
                         save.run = TRUE, infection.FUN = infection.net,
                         logger.FUN = step_logger,
                         .tracker.list = list(cur = get_current_timestep))
  sim <- netsim(est, param, init, control)
  expect_equal(sim$epi$at_arg[, 1], c(NA, 2:5))
  expect_equal(sim$epi$flow_like[, 1], c(NA, 2:5) * 10)
  expect_equal(sim$epi$cur[, 1], c(NA, 2:5))

  control$start <- 6
  control$nsteps <- 8
  rs <- netsim(sim, param, init, control)
  expect_equal(rs$epi$at_arg[, 1], c(NA, 2:8))
  expect_equal(rs$epi$cur[, 1], c(NA, 2:8))
})

context("Custom Trackers")

test_that("netsim, SI, custom trackers", {
  skip_on_cran()
  nw <- network_initialize(n = 50)
  nw <- set_vertex_attribute(nw, "race", rbinom(50, 1, 0.5))
  est <- netest(
    nw,
    formation = ~edges,
    target.stats = 25,
    coef.diss = dissolution_coefs(~ offset(edges), 10, 0),
    verbose = FALSE
  )

  init <- init.net(i.num = 10)

  epi_s_num <- function(dat, at) {
    needed_attributes <- c("status")
    output <- with(get_attr_list(dat, needed_attributes), {
      sum(status == "s", na.rm = TRUE)
    })
    return(output)
  }

  epi_prop_infected <- function(dat, at) {
    needed_attributes <- c("status", "active")
    output <- with(get_attr_list(dat, needed_attributes), {
      pop <- active == 1
      cond <- status == "i"
      sum(cond & pop, na.rm = TRUE) / sum(pop, na.rm = TRUE)
    })
    return(output)
  }

  some.trackers <- list(
    prop_infected = epi_prop_infected,
    s_num = epi_s_num
  )

  control <- control.net(
    type = "SI",
    nsims = 1,
    nsteps = 50,
    verbose = FALSE,
    infection.FUN = infection.net,
    .tracker.list = some.trackers
  )

  param <- param.net(
    inf.prob = 0.3,
    act.rate = 0.1
  )

  mod <- netsim(est, param, init, control)

  d <- as.data.frame(mod)

  expect_true(all(c("prop_infected", "s_num") %in% names(d)))
  expect_is(d[["prop_infected"]], "numeric")
  # the custom epi trackers are not run during intialization so the first value
  # is always NA
  expect_true(all(d$s_num[2:50] == d$s.num[2:50]))
})

context("Load Parameters from data.frame")

test_that("Load parameters from data.frame", {
  skip_on_cran()
  params.df <- dplyr::tribble(
    ~param, ~value, ~type, ~detail,
    "p1", "10", "numeric", "foo",
    "p2", "TRUE", "logical", "bar",
    "p3_1", "1", "numeric", "baz",
    "p3_2", "3", "numeric", "foobar",
    "p4", "tsa", "character", "foobaz"
  )

  expect_silent(param <- param.net_from_table(params.df))
  expect_s3_class(param, "param.net")
  expect_type(param, "list")

  expect_silent(param <- param.net(data.frame.params = params.df))
  expect_s3_class(param, "param.net")
  expect_type(param, "list")
  # Verify the table was actually unpacked (not just stored as a dot arg)
  expect_equal(param$p1, 10)
  expect_true(param$p2)
  expect_equal(param$p3, c(1, 3))
  expect_equal(param$p4, "tsa")

  # convert back to a `long.param.df`
  param.df_back <- param.net_from_table(params.df) |> param.net_to_table()
  expect_true(all(params.df[c("param", "value", "type")] == param.df_back))

  # wrong column name
  params.df <- dplyr::tribble(~name, ~value, ~type, "p1", "10", "numeric")
  expect_error(param <- param.net_from_table(params.df))
  params.df <- dplyr::tribble(~param, ~val, ~type, "p1", "10", "numeric")
  expect_error(param <- param.net_from_table(params.df))
  params.df <- dplyr::tribble(~param, ~value, ~class, "p1", "10", "numeric")
  expect_error(param <- param.net_from_table(params.df))

  # wrong "type" value
  params.df <- dplyr::tribble(
    ~param, ~value, ~type, ~detail,
    "p1", "10", "numeric", "foo",
    "p2", "TRUE", "logical", "bar",
    "p3_1", "1", "factor", "baz",
    "p3_2", "3", "numeric", "foobar",
    "p4", "tsa", "character", "foobaz"
  )
  expect_error(param <- param.net_from_table(params.df))

  # wrong "param" format
  params.df <- dplyr::tribble(
    ~param, ~value, ~type, ~detail,
    ".p1", "10", "numeric", "foo",
    "p2", "TRUE", "logical", "bar",
    "p_3_1", "1", "numeric", "baz",
    "p3_2", "3", "numeric", "foobar",
    "p4", "tsa", "character", "foobaz"
  )
  expect_error(param <- param.net_from_table(params.df))
})

test_that("parameter vectors with ten or more elements round trip", {
  skip_on_cran()

  # The position suffix pattern was `(_[1-9]+)?`, which rejected any position
  # containing a zero. `x_10` failed while `x_11` passed.
  p <- param.net(inf.prob = 0.3, x = 1:10)
  tbl <- param.net_to_table(p)
  expect_true(all(paste0("x_", 1:10) %in% tbl$param))

  back <- param.net_from_table(tbl)
  expect_equal(back$x, as.numeric(1:10))
  expect_equal(back$inf.prob, 0.3)

  # Three-digit positions, and positions that are multiples of ten.
  p100 <- param.net(y = seq_len(100) / 10)
  expect_equal(param.net_from_table(param.net_to_table(p100))$y,
               seq_len(100) / 10)

  # Non-numeric vectors longer than nine elements.
  p.chr <- param.net(z = letters[1:12])
  expect_equal(param.net_from_table(param.net_to_table(p.chr))$z, letters[1:12])

  # Zero and leading-zero positions are still malformed.
  for (nm in c("x_0", "x_01", "x_007")) {
    bad <- data.frame(param = nm, value = "1", type = "numeric",
                      stringsAsFactors = FALSE)
    expect_error(param.net_from_table(bad), "malformed")
  }

  # Names that merely end in digits are unaffected.
  ok <- data.frame(param = c("beta2", "paramSet10"), value = c("1", "2"),
                   type = rep("numeric", 2), stringsAsFactors = FALSE)
  expect_silent(ok.param <- param.net_from_table(ok))
  expect_equal(ok.param$beta2, 1)
  expect_equal(ok.param$paramSet10, 2)

  # The scenario path shares `unflatten_params`.
  sc.df <- data.frame(.scenario.id = "s1", .at = 1, inf.prob = 0.3,
                      stringsAsFactors = FALSE)
  for (i in 1:12) sc.df[[paste0("v_", i)]] <- i
  sc <- create_scenario_list(sc.df)
  expect_equal(sc[[1]][[".param.updater.list"]][[1]]$param$v, as.numeric(1:12))
})

test_that("param.net_to_table refuses multilayer and duration parameters", {
  p <- param.net(inf.prob = multilayer(0.3, 0.05), act.rate = 1)
  expect_error(param.net_to_table(p), "parameters: `inf.prob`")
  p2 <- param.net(inf.prob = multilayer(0.3, c(0.1, 0.05)),
                  act.rate = multilayer(1, 2))
  expect_error(param.net_to_table(p2), "`inf.prob`, `act.rate`")
  p3 <- param.net(inf.prob = by_infection_duration(c(0.5, 0.1)), act.rate = 1)
  expect_error(param.net_to_table(p3), "by_infection_duration")

  # without the multilayer parameters the rest converts as before
  p2[c("inf.prob", "act.rate")] <- NULL
  expect_silent(tbl <- param.net_to_table(p2))
  expect_false(any(grepl("^(inf.prob|act.rate)", tbl$param)))
})

test_that("dot-prefixed parameter names are rejected as malformed", {
  skip_on_cran()

  # Dot-prefixed internals cannot come in through a parameter table, since a
  # flat parameter name must begin with a letter.
  for (nm in c(".param.updater.list", ".param.updater.list_1")) {
    df <- data.frame(param = c("inf.prob", nm), value = c("0.3", "1"),
                     type = rep("numeric", 2), stringsAsFactors = FALSE)
    expect_error(param.net_from_table(df), "malformed")
  }
})

test_that("param.net preserves data.frame.params values against constructor defaults", {
  skip_on_cran()

  # Regression for #1029: vector `act.rate` from a long.param.df must not be
  # silently overwritten by the scalar default of 1.
  params.df <- data.frame(
    param = c("act.rate_1", "act.rate_2", "act.rate_3", "act.rate_4",
              "inf.prob"),
    value = c("5", "1", "1", "1", "0.3"),
    type  = rep("numeric", 5),
    stringsAsFactors = FALSE
  )
  p <- param.net(data.frame.params = params.df)
  expect_equal(p$act.rate, c(5, 1, 1, 1))
  expect_equal(p$inf.prob, 0.3)

  # Scalar default still applies when neither formal arg nor table provides it.
  p_default <- param.net(inf.prob = 0.3)
  expect_equal(p_default$act.rate, 1)

  # Regression for #1031: `vital` must reflect table-supplied vital parameters
  # after a round trip through `param.net_to_table` -> `param.net`.
  p1 <- param.net(inf.prob = 0.3, a.rate = 0.1, ds.rate = 0.1, di.rate = 0.1)
  expect_true(p1$vital)
  p2 <- param.net(data.frame.params = param.net_to_table(p1))
  expect_true(p2$vital)

  # `vital` is FALSE when no vital params are supplied via either path.
  p_no_vital <- param.net(
    data.frame.params = data.frame(
      param = "inf.prob", value = "0.3", type = "numeric",
      stringsAsFactors = FALSE
    )
  )
  expect_false(p_no_vital$vital)
})

test_that("param.net rejects the misdocumented data.frame.parameters name", {
  skip_on_cran()

  # Regression for #1031: releases through v2.6.1 documented the table argument
  # as `data.frame.parameters`, which the constructor never accepted. It used to
  # fall through to `...` and be stored as a parameter of that name, leaving the
  # table unpacked and the model silently misparameterized.
  params.df <- data.frame(
    param = c("inf.prob", "act.rate"),
    value = c("0.3", "2"),
    type  = rep("numeric", 2),
    stringsAsFactors = FALSE
  )
  expect_error(
    param.net(data.frame.parameters = params.df),
    "Use data.frame.params instead"
  )

  # The correct name still works.
  p <- param.net(data.frame.params = params.df)
  expect_equal(p$inf.prob, 0.3)
  expect_equal(p$act.rate, 2)
})

# `module.order` validation (run at control.net() construction time) -----------

test_that("module.order with an entry that has no matching .FUN errors", {
  expect_error(
    control.net(type = "SI", nsteps = 10,
                module.order = c("resim_nets.FUN", "summary_nets.FUN",
                                 "infction.FUN", "nwupdate.FUN")),
    "no matching `\\.FUN` argument"
  )
})

test_that("module.order containing initialize.FUN or verbose.FUN errors", {
  expect_error(
    control.net(type = "SI", nsteps = 10,
                module.order = c("initialize.FUN", "resim_nets.FUN",
                                 "summary_nets.FUN", "infection.FUN",
                                 "nwupdate.FUN", "prevalence.FUN")),
    "`initialize\\.FUN` runs once at simulation start"
  )

  expect_error(
    control.net(type = "SI", nsteps = 10,
                module.order = c("resim_nets.FUN", "summary_nets.FUN",
                                 "infection.FUN", "nwupdate.FUN",
                                 "prevalence.FUN", "verbose.FUN")),
    "outside the per-step module loop"
  )
})

test_that("module.order omitting a critical built-in warns", {
  expect_warning(
    control.net(type = "SI", nsteps = 10,
                module.order = c("resim_nets.FUN", "summary_nets.FUN",
                                 "infection.FUN", "prevalence.FUN")),
    "nwupdate\\.FUN"
  )

  expect_warning(
    control.net(type = "SI", nsteps = 10,
                module.order = c("summary_nets.FUN", "infection.FUN",
                                 "nwupdate.FUN", "prevalence.FUN")),
    "resim_nets\\.FUN"
  )

  expect_warning(
    control.net(type = "SI", nsteps = 10,
                module.order = c("resim_nets.FUN", "infection.FUN",
                                 "nwupdate.FUN", "prevalence.FUN")),
    "summary_nets\\.FUN"
  )
})

test_that("module.order with the recommended full ordering emits no warning", {
  expect_no_warning(
    control.net(type = "SI", nsteps = 10,
                module.order = c("resim_nets.FUN", "summary_nets.FUN",
                                 "infection.FUN", "nwupdate.FUN",
                                 "prevalence.FUN"))
  )
})

test_that("module.order is not required (NULL skips validation)", {
  expect_no_warning(
    control.net(type = "SI", nsteps = 10)
  )
})

test_that("module.order does not warn for built-ins explicitly set to NULL", {
  # When a user disables a critical built-in by passing NULL, that module
  # leaves `bi.mods` and so should not trigger the missing-critical warning
  # if it is also absent from `module.order`.
  expect_no_warning(
    control.net(type = "SI", nsteps = 10,
                nwupdate.FUN = NULL,
                module.order = c("resim_nets.FUN", "summary_nets.FUN",
                                 "infection.FUN", "prevalence.FUN"))
  )
})

context("Custom Module Return Value Validation")

test_that("netsim errors clearly when a custom module omits return(dat)", {
  skip_on_cran()

  nw <- network_initialize(n = 50)
  est <- netest(nw,
                formation = ~edges,
                target.stats = 25,
                coef.diss = dissolution_coefs(~offset(edges), 10, 0),
                verbose = FALSE)

  # Custom module that forgets to return dat (returns NULL implicitly)
  bad_module <- function(dat, at) {
    status <- get_attr(dat, "status")
    invisible(NULL)
  }

  param <- param.net(inf.prob = 0.3)
  init <- init.net(i.num = 10)
  control <- control.net(type = NULL,
                         nsims = 1,
                         nsteps = 5,
                         infection.FUN = infection.net,
                         my_bad.FUN = bad_module,
                         verbose = FALSE)

  expect_error(
    suppressMessages(netsim(est, param, init, control)),
    "Module 'my_bad.FUN' did not return"
  )
})

test_that("netsim errors when a custom module returns a plain list", {
  skip_on_cran()

  nw <- network_initialize(n = 50)
  est <- netest(nw,
                formation = ~edges,
                target.stats = 25,
                coef.diss = dissolution_coefs(~offset(edges), 10, 0),
                verbose = FALSE)

  # Returns dat with its netsim_dat class stripped, mimicking a user who
  # rebuilds the object via list() rather than returning the input.
  unclass_module <- function(dat, at) {
    unclass(dat)
  }

  param <- param.net(inf.prob = 0.3)
  init <- init.net(i.num = 10)
  control <- control.net(type = NULL,
                         nsims = 1,
                         nsteps = 5,
                         infection.FUN = infection.net,
                         my_unclass.FUN = unclass_module,
                         verbose = FALSE)

  expect_error(
    suppressMessages(netsim(est, param, init, control)),
    "Module 'my_unclass.FUN' did not return"
  )
})

test_that("netsim_cond_msg names the condition, module, and step", {
  expect_equal(netsim_cond_msg("error", "my.FUN", 3),
               "\n\tAn error occurred in module 'my.FUN' at step 3")
  expect_equal(netsim_cond_msg("warning", "my.FUN", 3),
               "\n\tA warning occurred in module 'my.FUN' at step 3")
  expect_equal(netsim_cond_msg("message", "my.FUN", 3),
               "\n\tA message occurred in module 'my.FUN' at step 3")
})

test_that("netsim reports the module and step where an error occurred", {
  skip_on_cran()

  nw <- network_initialize(n = 50)
  est <- netest(nw,
                formation = ~edges,
                target.stats = 25,
                coef.diss = dissolution_coefs(~offset(edges), 10, 0),
                verbose = FALSE)

  failing_module <- function(dat, at) {
    if (at == 3) stop("boom")
    dat
  }

  param <- param.net(inf.prob = 0.3)
  init <- init.net(i.num = 10)
  control <- control.net(type = NULL,
                         nsims = 1,
                         nsteps = 5,
                         infection.FUN = infection.net,
                         my_failing.FUN = failing_module,
                         verbose = FALSE)

  msgs <- character(0)
  expect_error(
    withCallingHandlers(
      netsim(est, param, init, control),
      message = function(m) {
        msgs <<- c(msgs, conditionMessage(m))
        invokeRestart("muffleMessage")
      }
    ),
    "boom"
  )
  expect_true(any(grepl("An error occurred in module 'my_failing.FUN' at step 3",
                        msgs, fixed = TRUE)))
})

test_that("an attribute first appended mid-simulation covers every node", {
  skip_on_cran()

  ## Arrivals module creating a custom attribute that does not exist at
  ## initialization, so the first arrivals are the first nodes to hold it
  afunc <- function(dat, at) {
    n.new <- 2
    dat <- append_core_attr(dat, at, n.new)
    dat <- append_attr(dat, "status", "s", n.new)
    dat <- append_attr(dat, "infTime", NA, n.new)
    dat <- append_attr(dat, "arrival.step", at, n.new)
    dat <- set_epi(dat, "a.flow", n.new)
    return(dat)
  }

  nw <- network_initialize(n = 20)
  est <- netest(nw, formation = ~edges, target.stats = 5,
                coef.diss = dissolution_coefs(~offset(edges), 10),
                verbose = FALSE)
  param <- param.net(inf.prob = 0.3)
  init <- init.net(i.num = 2)
  control <- control.net(type = NULL, nsims = 1, nsteps = 3,
                         arrivals.FUN = afunc, infection.FUN = infection.net,
                         resimulate.network = TRUE, tergmLite = FALSE,
                         save.run = TRUE, verbose = FALSE)
  mod <- netsim(est, param, init, control)

  ## the modules run from the second time step, so 2 arrival steps of 2 nodes
  attrs <- mod$run[[1]]$attr
  n.nodes <- length(attrs$active)
  expect_equal(n.nodes, 20 + 2 * 2)
  expect_true(all(vapply(attrs, length, integer(1)) == n.nodes))
  ## the 20 initial nodes have no arrival step, the 4 arrivals do
  expect_equal(attrs$arrival.step, c(rep(NA, 20), rep(2:3, each = 2)))
})