R/netclique.R

In EpiModel: Mathematical Modeling of Infectious Disease Dynamics

#' @title Clique Layer for Network Epidemic Models
#'
#' @description Builds a network layer in which every pair of nodes sharing a
#'              value of a grouping attribute is connected, so that each group
#'              (a household, classroom, hospital ward, or cabin) is a clique.
#'              [netsim()] accepts the layer anywhere in its layer list
#'              alongside the layers estimated with [netest()]. The edges of a
#'              clique layer never form or dissolve: there is no formation
#'              model, no dissolution model, and no resimulation. It is an
#'              addition to `netest` for structures that are groups by
#'              definition, not a substitute for an ERGM (see Details).
#'
#' @param nw An object of class `network` or `networkLite` holding the node
#'        set and any vertex attributes, as passed to [netest()] for the other
#'        layers of the same model. The network must be undirected and not bipartite. Any edges
#'        on `nw` are ignored.
#' @param group.attr Name of a vertex attribute on `nw` whose values partition
#'        the nodes into groups. Every pair of nodes sharing a non-missing
#'        value becomes an edge. Nodes with a missing (`NA`) value are in no
#'        group and are isolates on this layer.
#' @param arrivals Rule for placing nodes that arrive during a simulation with
#'        vital dynamics: `"isolate"` (the default; the new node has no edges
#'        on this layer and its group attribute is `NA`), `"new"` (each
#'        arrival starts a group of size one, with a fresh group id), or
#'        `"join"` (each arrival joins an existing group and is connected to
#'        all of its current members). See Details.
#' @param arrivals.FUN Optional function used with `arrivals = "join"` to
#'        choose the group each arrival joins. See Details.
#'
#' @details
#' A `netest` object carries a formation model, a dissolution model, and a
#' starting network; `netsim` resimulates its layer each time step. That is
#' the right representation whenever tie existence is something to model: a
#' network with a degree distribution, mixing by attribute, and clustering to
#' reproduce is an ERGM's job whether its ties turn over quickly, slowly, or
#' not at all, and a long partnership duration in [dissolution_coefs()] keeps
#' a `netest` layer close to fixed.
#'
#' Some contact structures are groups by construction instead: every pair of
#' co-residents is a household contact for the whole simulation, every pair
#' of pupils in a classroom is a classroom contact, every pair of cabin-mates
#' shares a cabin. There is no tie-formation process to estimate, and an ERGM
#' can reproduce such cliques only through a `nodematch` term on the group
#' attribute targeted at its maximum, the number of within-group pairs. A
#' target at the maximum of a statistic has no finite coefficient, so the fit
#' runs to its iteration limit, the coefficient it stops at is arbitrary, and
#' a few groups are typically left incomplete. Before this function, the
#' alternative was to carry the clique edgelist into `netsim` as a parameter
#' and walk it in a custom infection module. `netclique` is for these cases.
#' It builds a layer object that fills a slot in the layer list, is skipped
#' by the network resimulation and the edges correction, and is read by the
#' built-in infection modules through [discord_edgelist()] like any other
#' layer.
#'
#' A network whose edges are given rather than generated, such as a fully
#' observed contact network, is a different case again: it may carry its own
#' edge dynamics and it skips estimation for a different reason. That case is
#' handled by [netcensus()].
#'
#' The group attribute must be an integer, numeric, or character vector; a
#' factor is refused because arrivals under `"new"` need to create values the
#' factor does not have. It should also be set on the network passed to
#' [netest()] for the other layers, since `netsim` reads nodal attributes from
#' the first layer in its list; when it is not, `netsim` copies it from the
#' clique layer.
#'
#' @section Arrivals:
#' When a model has vital dynamics, each arriving node has to be placed on the
#' clique layer. `netsim` applies the rule chosen with `arrivals` inside
#' [arrive_nodes()], once per clique layer, after the arrivals module has
#' created the node and set its other attributes:
#'
#'  * `"isolate"`: the node gets no edges and its group attribute is `NA`. A
#'    custom module may place it later.
#'  * `"new"`: the node gets a fresh group id and no edges. Fresh ids
#'    continue from the largest id used so far in the simulation, so an id is
#'    never given to a second group, even after every member of the first has
#'    departed; for character ids they are new strings of the form
#'    `arrival_<i>`. Groups then grow only through further `"join"` arrivals,
#'    so this rule suits models in which arrivals are single-person
#'    households.
#'  * `"join"`: the node joins an existing group and is connected to every
#'    active member of that group, including other nodes joining the same
#'    group in the same time step. By default the group is drawn with
#'    probability proportional to its current size, which is the distribution
#'    of a randomly chosen existing member. `arrivals.FUN` replaces that draw:
#'    it is called as `arrivals.FUN(dat, at, new_ids, network)` and must return
#'    one group id per element of `new_ids` (an `NA` leaves that node as an
#'    isolate). Because the function receives the full `dat` object it can
#'    read any nodal attribute; the example below places each newborn in a
#'    household that already has a young child, and a function that reads the
#'    group attribute back for `new_ids` lets a custom arrivals module set the
#'    group itself. A function that draws from existing groups should leave
#'    `new_ids` and inactive nodes out of its pool: when an arrivals module
#'    gives new nodes a default value of the grouping attribute, such as `0`,
#'    a pool taken from all nodes contains that value, and every arrival that
#'    draws it joins the other arrivals in a spurious group.
#'
#' The layer keeps `arrivals.FUN` with its enclosing environment, and so does
#' every `netsim` object simulated from it. Define the function at the top
#' level of a script or in a package, not inside another function whose local
#' objects (fitted networks, population data) would then be saved with each
#' simulation.
#'
#' The rule decides the group attribute of every arriving node. A value given
#' by `attr.rules` in [control.net()] or by a custom arrivals module is
#' replaced, so that the layer always matches the groups, and `netsim` warns
#' when `attr.rules` names the group attribute. The exception is `"join"` with
#' an `arrivals.FUN` that reads that value back, as described above.
#'
#' @section Transmission over a clique layer:
#' The built-in infection modules treat every layer alike, so a clique layer
#' transmits with the model's `inf.prob` and `act.rate` unless those are given
#' per layer as [multilayer()] objects in [param.net()]:
#' `param.net(inf.prob = multilayer(0.45, 0.10), act.rate = multilayer(1, 2))`
#' with the household layer first and the community layer second. The
#' transmission matrix records the layer of each transmission in its
#' `network` column.
#'
#' Nodes departing the population are removed from a clique layer with their
#' edges, as on every other layer. Network statistics for the layer are
#' recorded through `nwstats.formula` in [control.net()], with the default
#' `"formation"` meaning `~edges` for a clique layer, and the cumulative
#' edgelist records its edges with a start time of 0. [netdx()] refuses a
#' `netclique` object, since there is no model to diagnose; the group size
#' distribution and mean degree are shown by `print`.
#'
#' @return
#' An object of class `netclique`, a list with elements:
#'
#'  * **newnetwork:** a `networkLite` object holding the node set, the vertex
#'    attributes of `nw`, and the clique edges. A `networkLite` stores the
#'    edges as an edgelist, so the layer stays small for large populations;
#'    [trim_netest()] has nothing further to remove from it.
#'  * **group.attr:** the grouping attribute name.
#'  * **arrivals**, **arrivals.FUN:** the arrival rule.
#'  * **target.stats**, **target.stats.names:** the edge count of the layer,
#'    named `"edges"`, so that the network statistics table printed by
#'    [print.netsim()] and plotted by [plot.netsim()] shows the layer's own
#'    edge count as the target.
#'  * **edapprox:** `FALSE`.
#'  * **summary:** a list with the node count, number of groups, edge count,
#'    mean degree, number of isolates, and the group size distribution.
#'
#' @seealso [sample_groups()] builds a population from a table of group
#'   types, and [assign_groups()] assigns group ids to an existing population
#'   with composition constraints. [netsim()] runs the simulation, and
#'   [multilayer()] sets per-layer parameters and controls.
#'
#' @export
#'
#' @examples
#' \donttest{
#' # Households as cliques under a community TERGM layer
#' hh <- sample_groups(500, c("adult" = 0.3, "adult adult" = 0.4,
#'                            "adult adult child" = 0.2,
#'                            "adult adult child child" = 0.1),
#'                     attr.name = "age")
#' nw <- network_initialize(n = 500)
#' nw <- set_vertex_attribute(nw, "age", hh$age)
#' nw <- set_vertex_attribute(nw, "hh_id", hh$group)
#'
#' est_hh <- netclique(nw, group.attr = "hh_id", arrivals = "join")
#' est_hh
#' print(est_hh, by = "age")
#'
#' est_com <- netest(nw, formation = ~edges + nodematch("age"),
#'                   target.stats = c(250, 175),
#'                   coef.diss = dissolution_coefs(~offset(edges), 20),
#'                   verbose = FALSE)
#'
#' param <- param.net(inf.prob = multilayer(0.3, 0.05), act.rate = 1)
#' init <- init.net(i.num = 10)
#' control <- control.net(type = "SI", nsteps = 50, nsims = 1,
#'                        tergmLite = TRUE, resimulate.network = TRUE,
#'                        verbose = FALSE)
#' sim <- netsim(list(est_hh, est_com), param, init, control)
#'
#' # Share of transmissions on the household layer
#' tm <- get_transmat(sim)
#' mean(tm$network == 1)
#'
#' # Newborns join a household that already has a child
#' est_hh2 <- netclique(nw, group.attr = "hh_id", arrivals = "join",
#'   arrivals.FUN = function(dat, at, new_ids, network) {
#'     active <- get_attr(dat, "active")
#'     hh <- get_attr(dat, "hh_id")
#'     age <- get_attr(dat, "age")
#'     existing <- active == 1
#'     existing[new_ids] <- FALSE
#'     pool <- hh[which(existing & age == "child" & !is.na(hh))]
#'     pool[sample.int(length(pool), length(new_ids), replace = TRUE)]
#'   })
#' }
#'
netclique <- function(nw, group.attr, arrivals = c("isolate", "new", "join"),
                      arrivals.FUN = NULL) {

  if (!inherits(nw, "network")) {
    stop("`nw` must be an object of class `network`.")
  }
  if (is.directed(nw)) {
    stop("`netclique` supports undirected networks only.")
  }
  if (is.bipartite(nw)) {
    stop("`netclique` does not support bipartite networks.")
  }
  arrivals <- match.arg(arrivals)
  if (!is.null(arrivals.FUN)) {
    if (!is.function(arrivals.FUN)) {
      stop("`arrivals.FUN` must be a function.")
    }
    if (arrivals != "join") {
      stop("`arrivals.FUN` applies only with `arrivals = \"join\"`.")
    }
  }
  if (missing(group.attr) || !is.character(group.attr) ||
        length(group.attr) != 1) {
    stop("`group.attr` must be a single attribute name.")
  }
  if (!group.attr %in% list.vertex.attributes(nw)) {
    stop("There is no vertex attribute called `", group.attr, "` on `nw`.")
  }
  group <- get_vertex_attribute(nw, group.attr)
  if (is.factor(group)) {
    stop("The `", group.attr, "` attribute is a factor. Convert it with ",
         "`as.character()` or `as.integer()` before calling `netclique`.")
  }
  if (!is.numeric(group) && !is.character(group)) {
    stop("The `", group.attr, "` attribute must be numeric or character.")
  }

  n <- network.size(nw)
  el <- group_edgelist(group)
  sizes <- tabulate(factor(group[!is.na(group)]))

  ## the layer network: node set and attributes of nw, clique edges only. It
  ## is a networkLite: a `network` object keeps a list for every edge, which
  ## for a household layer of a few hundred thousand edges is over a gigabyte
  attrs <- list.vertex.attributes(nw)
  attr_list <- lapply(attrs, function(a) {
    get.vertex.attribute(nw, a, null.na = FALSE, unlist = FALSE)
  })
  names(attr_list) <- attrs
  newnetwork <- networkLite(structure(el, n = n), attr_list, atomize = TRUE)

  deg <- tabulate(c(el[, 1], el[, 2]), nbins = n)

  out <- list()
  out$newnetwork <- newnetwork
  out$group.attr <- group.attr
  out$arrivals <- arrivals
  out$arrivals.FUN <- arrivals.FUN
  out$target.stats <- nrow(el)
  out$target.stats.names <- "edges"
  out$edapprox <- FALSE
  out$summary <- list(n = n,
                      groups = length(sizes),
                      edges = nrow(el),
                      mean.degree = mean(deg),
                      isolates = sum(deg == 0),
                      group.sizes = table(sizes, dnn = NULL))

  class(out) <- "netclique"
  return(out)
}

#' @export
print.netclique <- function(x, by = NULL, digits = 3, ...) {

  s <- x$summary

  cat("EpiModel Clique Network Layer")
  cat("\n=======================")
  cat("\nModel class:", class(x))
  cat("\nGrouping attribute:", x$group.attr)

  cat("\n\nLayer Summary")
  cat("\n-----------------------")
  cat("\nNodes:", s$n)
  cat("\nGroups:", s$groups)
  cat("\nEdges:", s$edges)
  cat("\nMean degree:", round(s$mean.degree, digits))
  cat("\nIsolates:", s$isolates)

  cat("\n\nGroup Size Distribution")
  cat("\n-----------------------\n")
  gs <- rbind(size = as.integer(names(s$group.sizes)),
              groups = as.integer(s$group.sizes))
  colnames(gs) <- rep("", ncol(gs))
  print(gs)

  if (!is.null(by)) {
    if (is.null(x$newnetwork)) {
      cat("\nMean degree by `", by, "` is not available: the layer ",
          "network is not stored on this object.\n", sep = "")
    } else {
      if (!by %in% list.vertex.attributes(x$newnetwork)) {
        stop("There is no vertex attribute called `", by, "` on the layer.")
      }
      el <- as.edgelist(x$newnetwork)
      deg <- tabulate(c(el[, 1], el[, 2]), nbins = network.size(x$newnetwork))
      by.val <- get_vertex_attribute(x$newnetwork, by)
      cat("\nMean Degree by `", by, "`", sep = "")
      cat("\n-----------------------\n")
      print(round(tapply(deg, by.val, mean), digits))
    }
  }

  cat("\nArrivals:", x$arrivals)
  if (x$arrivals == "join") {
    if (is.null(x$arrivals.FUN)) {
      cat(" (group drawn in proportion to its size)")
    } else {
      cat(" (group chosen by arrivals.FUN)")
    }
  }
  cat("\n")

  invisible()
}

# A layer with no network model behind it (a clique layer or an observed
# census), which the network simulation and the edges correction skip. The
# predicate is applied both to the layer object passed to netsim() and to its
# nwparam record on the dat object.
is_model_free_layer <- function(x) {
  inherits(x, c("netclique", "netcensus"))
}

# A clique layer specifically: the only model-free layer with a grouping
# attribute and arrival rules.
is_clique_layer <- function(x) {
  inherits(x, "netclique")
}
# All within-group pairs for a grouping vector, as a sorted two-column integer
# matrix with tail < head. Nodes with NA are in no group. Groups of the same
# size are processed together, so the cost is one combn() per distinct size.
group_edgelist <- function(group) {
  ids <- which(!is.na(group))
  empty <- matrix(integer(0), ncol = 2)
  if (length(ids) == 0) {
    return(empty)
  }
  g <- group[ids]
  o <- order(g, ids)
  ids <- ids[o]
  g <- g[o]
  sizes <- tabulate(factor(g, levels = unique(g)))
  starts <- cumsum(c(1L, sizes))[seq_along(sizes)]

  pieces <- list()
  for (k in sort(unique(sizes[sizes > 1]))) {
    gk <- which(sizes == k)
    members <- matrix(ids[rep(starts[gk], each = k) + rep(0:(k - 1), length(gk))],
                      ncol = k, byrow = TRUE)
    pairs <- combn(k, 2)
    pieces[[length(pieces) + 1]] <- cbind(
      as.vector(members[, pairs[1, ], drop = FALSE]),
      as.vector(members[, pairs[2, ], drop = FALSE])
    )
  }
  if (length(pieces) == 0) {
    return(empty)
  }
  el <- do.call(rbind, pieces)
  el <- cbind(pmin(el[, 1], el[, 2]), pmax(el[, 1], el[, 2]))
  el <- el[order(el[, 1], el[, 2]), , drop = FALSE]
  storage.mode(el) <- "integer"
  dimnames(el) <- NULL
  el
}

# Append rows to a tergmLite edgelist, keeping its attributes (n and the
# as.edgelist class and metadata) and its sorted order. A clique layer can
# hold hundreds of thousands of edges and gains a few each step, so the new
# rows are merged into place by position instead of sorting the whole
# edgelist again; the result is the same as rbind() followed by a stable sort.
add_edges_to_el <- function(el, new) {
  if (NROW(new) == 0) {
    return(el)
  }
  a <- attributes(el)
  old <- matrix(el, ncol = 2)
  new <- matrix(new, ncol = 2)
  new <- new[order(new[, 1], new[, 2]), , drop = FALSE]
  base <- as.numeric(max(old, new)) + 1
  key.old <- old[, 1] * base + old[, 2]
  if (is.unsorted(key.old)) {
    out <- rbind(old, new)
    out <- out[order(out[, 1], out[, 2]), , drop = FALSE]
  } else {
    pos <- findInterval(new[, 1] * base + new[, 2], key.old) +
      seq_len(nrow(new))
    out <- matrix(0L, nrow(old) + nrow(new), 2)
    out[pos, ] <- new
    out[-pos, ] <- old
  }
  a$dim <- dim(out)
  a$dimnames <- NULL
  attributes(out) <- a
  out
}

# Add edges to a clique layer in either storage mode.
add_clique_edges <- function(dat, network, el_new) {
  if (NROW(el_new) == 0) {
    return(dat)
  }
  at <- get_current_timestep(dat)
  el_new <- cbind(pmin(el_new[, 1], el_new[, 2]), pmax(el_new[, 1], el_new[, 2]))
  el_new <- unique(el_new)

  if (get_control(dat, "tergmLite") == TRUE) {
    dat$run$el[[network]] <- add_edges_to_el(dat$run$el[[network]], el_new)
    if (get_network_control(dat, network, "tergmLite.track.duration") == TRUE) {
      dat$run$net_attr[[network]][["lasttoggle"]] <- rbind(
        dat$run$net_attr[[network]][["lasttoggle"]],
        cbind(el_new, at)
      )
    }
  } else {
    nw <- get_network(dat, network = network)
    nw <- networkDynamic::add.edges.active(nw, tail = el_new[, 1],
                                           head = el_new[, 2],
                                           onset = at, terminus = Inf)
    dat <- set_network(dat, nw = nw, network = network)
  }
  return(dat)
}

# Fresh group ids for k arrivals, of the same type as the existing ids.
# `last` is the running maximum kept by record_group_ids(), so that an id is
# never given to a second group, even after every member of the first has
# departed: for numeric ids it is the largest id used so far, and for
# character ids the largest i of the "arrival_<i>" ids used so far.
new_group_ids <- function(group, k, last = NULL) {
  existing <- group[!is.na(group)]
  if (is.numeric(group) || length(existing) == 0) {
    return(max(c(0, last, existing)) + seq_len(k))
  }
  ids <- character(0)
  i <- NVL(last, 0)
  while (length(ids) < k) {
    i <- i + 1
    id <- paste0("arrival_", i)
    if (!id %in% existing) {
      ids <- c(ids, id)
    }
  }
  ids
}

# Raise the running maximum of the ids of a grouping attribute, kept on
# dat$run so that it carries over a restart, to cover the ids just assigned.
# Character ids count only in the "arrival_<i>" form that new_group_ids()
# creates.
record_group_ids <- function(dat, group.attr, ids) {
  ids <- ids[!is.na(ids)]
  if (is.character(ids)) {
    ids <- ids[grepl("^arrival_[0-9]+$", ids)]
    ids <- as.numeric(sub("^arrival_", "", ids))
  }
  if (length(ids) > 0) {
    dat$run$group_id_max[[group.attr]] <-
      max(c(dat$run$group_id_max[[group.attr]], ids))
  }
  return(dat)
}

# Apply a clique layer's arrival rule to the nodes new_ids, which arrive_nodes
# has already added to the layer as isolates. Sets the group attribute for
# the new nodes and, under "join", connects each of them to every active
# member of its group, including other nodes joining the same group in this
# time step.
clique_layer_arrivals <- function(dat, network, new_ids) {
  nwparam <- get_nwparam(dat, network = network)
  group.attr <- nwparam$group.attr
  rule <- nwparam$arrivals
  if (is.null(group.attr) || length(new_ids) == 0) {
    return(dat)
  }

  active <- get_attr(dat, "active")
  n <- length(active)
  group <- get_attr(dat, group.attr, override.null.error = TRUE)
  if (is.null(group)) {
    group <- rep(NA, n)
  } else if (length(group) < n) {
    group <- c(group, rep(NA, n - length(group)))
  }

  ## the rule sets the group of every arrival, replacing a value that
  ## attr.rules or a custom arrivals module gave it; an isolate has no group
  if (rule == "isolate") {
    group[new_ids] <- NA
    dat <- set_attr(dat, group.attr, group)
    return(dat)
  }

  last <- dat$run$group_id_max[[group.attr]]
  if (rule == "new") {
    group[new_ids] <- new_group_ids(group, length(new_ids), last)
    dat <- set_attr(dat, group.attr, group)
    dat <- record_group_ids(dat, group.attr, group[new_ids])
    return(dat)
  }

  ## the join rule: choose a group for each arrival, then connect it
  if (!is.null(nwparam$arrivals.FUN)) {
    at <- get_current_timestep(dat)
    new_group <- nwparam$arrivals.FUN(dat, at, new_ids, network)
    if (length(new_group) != length(new_ids)) {
      stop("`arrivals.FUN` for clique layer ", network, " returned ",
           length(new_group), " group ids for ", length(new_ids),
           " arriving nodes.")
    }
  } else {
    old_ids <- setdiff(which(active == 1), new_ids)
    pool <- group[old_ids]
    pool <- pool[!is.na(pool)]
    if (length(pool) == 0) {
      new_group <- new_group_ids(group, length(new_ids), last)
    } else {
      new_group <- pool[sample.int(length(pool), length(new_ids),
                                   replace = TRUE)]
    }
  }
  group[new_ids] <- new_group
  dat <- set_attr(dat, group.attr, group)
  dat <- record_group_ids(dat, group.attr, new_group)

  joining <- new_ids[!is.na(new_group)]
  if (length(joining) == 0) {
    return(dat)
  }
  ## split only the groups being joined, which keeps the cost from growing
  ## with the population
  in_group <- which(active == 1 & group %in% group[joining])
  members <- split(in_group, group[in_group])
  el_new <- lapply(joining, function(i) {
    m <- members[[as.character(group[i])]]
    m <- m[m != i]
    if (length(m) == 0) {
      return(NULL)
    }
    cbind(pmin(i, m), pmax(i, m))
  })
  el_new <- do.call(rbind, el_new)
  dat <- add_clique_edges(dat, network, el_new)

  return(dat)
}

# Remove every clique-layer edge that has an endpoint in ids, in either
# storage mode. Under tergmLite the edgelist keeps its attributes and sorted
# order, and the duration record, when tracked, drops the same dyads.
remove_clique_edges <- function(dat, network, ids) {
  if (length(ids) == 0) {
    return(dat)
  }
  at <- get_current_timestep(dat)

  if (get_control(dat, "tergmLite") == TRUE) {
    el <- dat$run$el[[network]]
    drop <- el[, 1] %in% ids | el[, 2] %in% ids
    if (any(drop)) {
      a <- attributes(el)
      out <- matrix(el, ncol = 2)[!drop, , drop = FALSE]
      a$dim <- dim(out)
      a$dimnames <- NULL
      attributes(out) <- a
      dat$run$el[[network]] <- out
    }
    if (get_network_control(dat, network, "tergmLite.track.duration") == TRUE) {
      lt <- dat$run$net_attr[[network]][["lasttoggle"]]
      keep <- !(lt[, 1] %in% ids | lt[, 2] %in% ids)
      dat$run$net_attr[[network]][["lasttoggle"]] <- lt[keep, , drop = FALSE]
    }
  } else {
    nw <- get_network(dat, network = network)
    eids <- unique(unlist(lapply(ids, function(i) {
      networkDynamic::get.edgeIDs.active(nw, v = i, at = at)
    })))
    if (length(eids) > 0) {
      nw <- networkDynamic::deactivate.edges(nw, onset = at, terminus = Inf,
                                             e = eids)
      dat <- set_network(dat, nw = nw, network = network)
    }
  }
  return(dat)
}

#' @title Move Nodes to Another Group of a Clique Layer
#'
#' @description Changes the group of existing nodes during a simulation and
#'              rewires a [netclique()] layer to match: each moved node loses
#'              its edges to the members of its old group and is connected to
#'              every active member of its new group. Custom modules use it to
#'              represent people changing groups, such as a young adult leaving
#'              home, a couple forming a household, or a child moving in with
#'              relatives.
#'
#' @param dat Main `netsim_dat` object passed through [netsim()] calls.
#' @param ids Positional ids of the nodes to move. Every node must be active.
#' @param group New group ids: one per element of `ids`, or a single value for
#'        all of them. An id in use joins that group; an id not in use starts
#'        a new group, which nodes given the same new id in one call form
#'        together; `NA` takes the node out of its group and leaves it as an
#'        isolate on the layer.
#' @param network Index of the clique layer in the list of layers passed to
#'        [netsim()]. It may be omitted when the model has exactly one clique
#'        layer. Every other clique layer built on the same grouping attribute
#'        is rewired as well.
#'
#' @details
#' The clique edges are built from the grouping attribute once, when the
#' simulation starts, and are afterwards changed only by arrivals, departures,
#' and this function. Setting the grouping attribute with [set_attr()] alone
#' therefore does not rewire the layer: the old edges remain and the layer no
#' longer matches the groups. `move_to_group()` sets the attribute and the
#' edges together, in either network storage mode (`tergmLite = TRUE` or
#' `FALSE`), so that the layer stays exactly the union of the cliques of the
#' current groups.
#'
#' Nodes that arrive in a time step join the layer when [nwupdate.net()]
#' runs, which places them by the layer's arrival rule (see [netclique()]), so
#' they cannot be moved before then. To place arriving nodes in chosen groups,
#' set the grouping attribute in the arrivals module and use
#' `arrivals = "join"` with an `arrivals.FUN` that returns that attribute for
#' the new nodes.
#'
#' A node moved to its own current group keeps its edges. When the moves leave
#' a group without members, the group ceases to exist, as when its last member
#' departs. For numeric group ids, one more than the largest id in use,
#' `max(get_attr(dat, group.attr), na.rm = TRUE) + 1`, is an unused id; ids of
#' groups whose members have all departed are also unused, so reusing one
#' starts a new group. The ids set here count toward the running maximum from
#' which the `"new"` arrival rule of [netclique()] draws, so a later arrival
#' is never given the id of a group formed by a move.
#'
#' Under `tergmLite = TRUE` with `tergmLite.track.duration = TRUE`, the moved
#' nodes' new edges are recorded as formed at the current time step. In the
#' cumulative edgelist (see [control.net()]), a move ends the node's old
#' edges and starts new ones.
#'
#' @return The updated `netsim_dat` object.
#'
#' @seealso [netclique()] for the clique layer and its rules for placing
#'   arriving nodes.
#'
#' @export
#'
#' @examples
#' \donttest{
#' # A module in which people aged 18 to 29 leave home with a probability of
#' # 0.05 per time step, each starting a household of their own
#' leave_home <- function(dat, at) {
#'   active <- get_attr(dat, "active")
#'   age <- get_attr(dat, "age")
#'   hh_id <- get_attr(dat, "hh_id")
#'   elig <- which(active == 1 & age >= 18 & age < 30)
#'   movers <- elig[runif(length(elig)) < 0.05]
#'   if (length(movers) > 0) {
#'     new_ids <- max(hh_id, na.rm = TRUE) + seq_along(movers)
#'     dat <- move_to_group(dat, ids = movers, group = new_ids)
#'   }
#'   return(dat)
#' }
#'
#' # 50 households of 4 under a community layer
#' set.seed(12)
#' nw <- network_initialize(n = 200)
#' nw <- set_vertex_attribute(nw, "age", sample(0:79, 200, replace = TRUE))
#' nw <- set_vertex_attribute(nw, "hh_id", rep(1:50, each = 4))
#' hh <- netclique(nw, group.attr = "hh_id")
#' est <- netest(nw, formation = ~edges, target.stats = 100,
#'               coef.diss = dissolution_coefs(~offset(edges), 20),
#'               verbose = FALSE)
#'
#' param <- param.net(inf.prob = multilayer(0.2, 0.05), act.rate = 1)
#' init <- init.net(i.num = 10)
#' control <- control.net(type = NULL, nsteps = 20, nsims = 1,
#'                        infection.FUN = infection.net,
#'                        leave_home.FUN = leave_home,
#'                        tergmLite = TRUE, resimulate.network = TRUE,
#'                        save.run = TRUE, verbose = FALSE)
#' sim <- netsim(list(hh, est), param, init, control)
#'
#' # More households at the end than at the start
#' length(unique(sim$run[[1]]$attr$hh_id))
#' }
#'
move_to_group <- function(dat, ids, group, network = NULL) {

  clique <- which(vapply(dat$nwparam, is_clique_layer, logical(1)))
  if (is.null(network)) {
    if (length(clique) != 1) {
      stop("The model has ", length(clique), " clique layers; give the ",
           "clique layer to change in `network`.")
    }
    network <- clique
  } else if (length(network) != 1 || !network %in% clique) {
    stop("Layer ", paste(network, collapse = ", "), " is not a clique layer.")
  }
  if (length(ids) == 0) {
    return(dat)
  }
  if (anyNA(ids) || anyDuplicated(ids) > 0) {
    stop("`ids` must not contain missing or duplicated values.")
  }
  if (length(group) == 1) {
    group <- rep(group, length(ids))
  }
  if (length(group) != length(ids)) {
    stop("`group` must have one value per element of `ids`, or a single ",
         "value.")
  }

  active <- get_attr(dat, "active")
  if (any(ids < 1 | ids > length(active)) || any(active[ids] != 1)) {
    stop("Only active nodes can be moved to another group.")
  }

  group.attr <- get_nwparam(dat, network = network)$group.attr
  grp <- get_attr(dat, group.attr)

  ## a node that an arrivals module has just created joins the layer when
  ## arrive_nodes() runs in nwupdate.net(), which places it by the layer's
  ## arrival rule; moving it before then would leave the layer out of step
  if (get_control(dat, "tergmLite") == TRUE) {
    n_layer <- dat$run$net_attr[[network]][["n"]]
  } else {
    n_layer <- network.size(get_network(dat, network = network))
  }
  if (any(ids > n_layer) || length(grp) < length(active)) {
    stop("Nodes that arrived in this time step cannot be moved before ",
         "nwupdate.net() has placed them on clique layer ", network, " by ",
         "its arrival rule. To place arriving nodes in chosen groups, set the `",
         group.attr, "` attribute in the arrivals module and use ",
         "arrivals = \"join\" with an arrivals.FUN that returns ",
         "get_attr(dat, \"", group.attr, "\")[new_ids]; or move them in a ",
         "later time step.", call. = FALSE)
  }

  ## every clique layer built on the same grouping attribute is rewired, so
  ## that they all stay the union of the cliques of the current groups
  layers <- clique[vapply(clique, function(k) {
    identical(get_nwparam(dat, network = k)$group.attr, group.attr)
  }, logical(1))]
  if (!all(is.na(group)) && is.numeric(grp) != is.numeric(group)) {
    stop("`group` must be of the same type as the `", group.attr,
         "` attribute (numeric or character).")
  }
  stay <- !is.na(group) & !is.na(grp[ids]) & grp[ids] == group
  ids <- ids[!stay]
  group <- group[!stay]
  if (length(ids) == 0) {
    return(dat)
  }
  grp[ids] <- group
  dat <- set_attr(dat, group.attr, grp)
  dat <- record_group_ids(dat, group.attr, group)

  ## edges from each moved node to every active member of its new group,
  ## including other nodes moved into the same group in this call
  joining <- ids[!is.na(group)]
  in_group <- which(active == 1 & grp %in% grp[joining])
  members <- split(in_group, grp[in_group])
  el_new <- do.call(rbind, lapply(joining, function(i) {
    m <- members[[as.character(grp[i])]]
    m <- m[m != i]
    if (length(m) == 0) {
      return(NULL)
    }
    cbind(pmin(i, m), pmax(i, m))
  }))

  if (!is.null(el_new)) {
    el_new <- unique(el_new)
  }
  for (layer in layers) {
    dat <- rewire_clique_layer(dat, layer, ids, el_new)
  }

  return(dat)
}

# Remove the clique-layer edges of the nodes ids and add the edges el_new, in
# either storage mode.
rewire_clique_layer <- function(dat, network, ids, el_new) {
  dat <- remove_clique_edges(dat, network, ids)
  if (is.null(el_new)) {
    return(dat)
  }
  if (get_control(dat, "tergmLite") == TRUE) {
    dat <- add_clique_edges(dat, network, el_new)
  } else {
    ## reactivate a dyad's earlier edge rather than add a second one
    nw <- get_network(dat, network = network)
    eid <- vapply(seq_len(nrow(el_new)), function(r) {
      e <- network::get.edgeIDs(nw, v = el_new[r, 1], alter = el_new[r, 2])
      if (length(e) == 0) NA_integer_ else as.integer(e[1])
    }, integer(1))
    if (any(!is.na(eid))) {
      at <- get_current_timestep(dat)
      nw <- networkDynamic::activate.edges(nw, onset = at, terminus = Inf,
                                           e = eid[!is.na(eid)])
      dat <- set_network(dat, nw = nw, network = network)
    }
    dat <- add_clique_edges(dat, network, el_new[is.na(eid), , drop = FALSE])
  }
  dat
}


#' @title Sample a Population of Groups from a Table of Group Types
#'
#' @description Builds a population one group at a time from a table of group
#'              types, such as household compositions, returning each node's
#'              group id and the attributes its group type implies. The result
#'              supplies both the grouping attribute for [netclique()] and the
#'              nodal attributes for the dynamic layers of the same model.
#'
#' @param n Number of nodes in the population.
#' @param types The group types. Either a character vector in which each
#'        element lists the members of one type separated by `sep`, such as
#'        `c("adult", "adult adult", "adult adult child")`; a named numeric
#'        vector whose names are those templates and whose values are the
#'        type weights, in which case `prob` is taken from the values; or a
#'        list of `data.frame`s, each with one row per member and any number of
#'        attribute columns, for types that carry more than one attribute per
#'        member.
#' @param prob Sampling weights for the types, one per type. They are
#'        normalized to sum to one. The default draws every type with equal
#'        probability.
#' @param attr.name Name of the attribute column in the result when `types` is
#'        a character vector.
#' @param sep Separator between member labels within a character template.
#'
#' @details
#' Types are drawn with replacement until the population reaches `n`. The
#' final group is drawn from the types whose size fits the remaining slots, so
#' that the population has exactly `n` nodes without cutting a group short.
#' When no type fits (no type has size one, and the remainder is smaller than
#' every type), the last group is truncated and a message says so.
#'
#' The realized distribution of group types, and so of group sizes and of the
#' attribute values, follows the `prob` weights up to sampling variation; the
#' expected mean group size is the weighted mean of the type sizes. The
#' expected share of nodes with a given attribute value is the weighted share
#' of that value among the members of all types. Tune the type table to hit
#' both.
#'
#' @return
#' A `data.frame` with one row per node and `n` rows, in group order:
#' `group` (integer group id, `1` to the number of groups), `type` (index of
#' the type the group was drawn from), and the attribute columns (one column
#' named `attr.name` for character templates, or the columns of the
#' `data.frame`s in `types`).
#'
#' @seealso [netclique()] to turn the `group` column into a clique
#'   layer, and [assign_groups()] for the reverse problem of assigning group
#'   ids to a population whose attributes already exist.
#'
#' @export
#'
#' @examples
#' # Household types by age group, weights chosen so that every child lives
#' # with at least one adult and about a quarter of older adults live alone
#' hh_types <- c("adult"                    = 0.15,
#'               "adult adult"              = 0.20,
#'               "elderly"                  = 0.12,
#'               "elderly elderly"          = 0.10,
#'               "adult elderly"            = 0.03,
#'               "adult child"              = 0.05,
#'               "adult adult child"        = 0.15,
#'               "adult adult child child"  = 0.15,
#'               "adult adult child elderly" = 0.05)
#' pop <- sample_groups(1000, hh_types, attr.name = "age")
#' head(pop, 10)
#' table(tabulate(pop$group))                 # household size distribution
#' prop.table(table(pop$age))                 # person-level age distribution
#'
#' # Every child shares a household with an adult
#' has_adult <- tapply(pop$age == "adult", pop$group, any)
#' all(has_adult[as.character(pop$group[pop$age == "child"])])
#'
#' # Types with more than one attribute per member
#' types <- list(
#'   data.frame(age = "adult", sex = "F"),
#'   data.frame(age = c("adult", "adult"), sex = c("F", "M")),
#'   data.frame(age = c("adult", "adult", "child"), sex = c("F", "M", "F"))
#' )
#' pop2 <- sample_groups(100, types, prob = c(0.3, 0.4, 0.3))
#' head(pop2)
#'
sample_groups <- function(n, types, prob = NULL, attr.name = "member",
                          sep = " ") {

  if (!is.numeric(n) || length(n) != 1 || n < 1 || n != round(n)) {
    stop("`n` must be a single positive integer.")
  }

  ## normalize types to a list of data.frames, one row per member
  if (is.character(types)) {
    if (!is.null(names(types))) {
      stop("`types` is a named character vector; give the weights in `prob` ",
           "or pass a named numeric vector.")
    }
    templates <- types
  } else if (is.numeric(types)) {
    if (is.null(names(types))) {
      stop("A numeric `types` must be named by the group-type templates.")
    }
    if (!is.null(prob)) {
      stop("Specify the type weights either as the values of `types` or in ",
           "`prob`, not both.")
    }
    prob <- unname(types)
    templates <- names(types)
  } else if (is.list(types)) {
    if (!all(vapply(types, is.data.frame, logical(1)))) {
      stop("A list `types` must contain one `data.frame` per group type.")
    }
    templates <- NULL
  } else {
    stop("`types` must be a character vector, a named numeric vector, or a ",
         "list of data.frames.")
  }
  if (!is.null(templates)) {
    members <- strsplit(trimws(templates), sep, fixed = TRUE)
    types <- lapply(members, function(m) {
      df <- data.frame(m, stringsAsFactors = FALSE)
      names(df) <- attr.name
      df
    })
  }
  K <- length(types)
  sizes <- vapply(types, nrow, integer(1))
  if (K == 0 || any(sizes == 0)) {
    stop("Every group type must have at least one member.")
  }
  if (any(c("group", "type") %in% names(types[[1]]))) {
    stop("The attribute columns of `types` cannot be named `group` or `type`.")
  }

  if (is.null(prob)) {
    prob <- rep(1, K)
  }
  if (length(prob) != K || any(is.na(prob)) || any(prob < 0) || sum(prob) == 0) {
    stop("`prob` must be one non-negative weight per group type.")
  }
  prob <- prob / sum(prob)

  ## draw whole groups until the population reaches n
  mean.size <- sum(prob * sizes)
  draw <- integer(0)
  total <- 0
  while (total < n) {
    n.draw <- ceiling(1.2 * (n - total) / mean.size) + 10
    batch <- sample.int(K, n.draw, replace = TRUE, prob = prob)
    cs <- total + cumsum(sizes[batch])
    keep <- which(cs <= n)
    draw <- c(draw, batch[keep])
    total <- if (length(keep) > 0) cs[max(keep)] else total
    if (length(keep) < n.draw) {
      break
    }
  }

  ## fill the remaining slots with types that fit; truncate as a last resort
  truncate <- 0L
  remaining <- n - total
  while (remaining > 0) {
    fit <- which(sizes <= remaining & prob > 0)
    if (length(fit) == 0) {
      message("No group type fits the ", remaining, " remaining slots; the ",
              "last group is truncated.")
      t <- sample.int(K, 1, prob = prob)
      truncate <- as.integer(sizes[t] - remaining)
      draw <- c(draw, t)
      remaining <- 0
    } else {
      t <- fit[sample.int(length(fit), 1, prob = prob[fit])]
      draw <- c(draw, t)
      remaining <- remaining - sizes[t]
    }
  }

  ## expand the drawn types into node rows, one block per type
  pieces <- lapply(seq_len(K), function(t) {
    gidx <- which(draw == t)
    if (length(gidx) == 0) {
      return(NULL)
    }
    k <- sizes[t]
    piece <- types[[t]][rep(seq_len(k), times = length(gidx)), , drop = FALSE]
    piece <- cbind(group = rep(gidx, each = k), type = t, piece)
    piece
  })
  out <- do.call(rbind, pieces)
  out <- out[order(out$group), , drop = FALSE]
  if (truncate > 0) {
    out <- out[seq_len(nrow(out) - truncate), , drop = FALSE]
  }
  rownames(out) <- NULL
  out
}


#' @title Assign Group Ids to an Existing Population
#'
#' @description Assigns nodes to groups of a target size distribution, subject
#'              to a composition rule that keeps dependent members (such as
#'              children) in groups that contain an anchor member (such as an
#'              adult). This is the reverse of [sample_groups()]: the nodal
#'              attributes already exist, from a census age distribution for
#'              example, and the group ids are built to match them.
#'
#' @param size.dist Distribution of group sizes: a numeric vector of weights
#'        named by group size, such as `c("1" = 0.28, "2" = 0.34, "3" = 0.15,
#'        "4" = 0.14, "5" = 0.09)`. An unnamed vector is taken over sizes
#'        `1, 2, ...`. Weights are normalized to sum to one.
#' @param role Optional vector of node roles, one per node, such as an age
#'        group. Its length sets the population size.
#' @param anchor Values of `role` that can anchor a group. Every group that
#'        contains a `dependent` member also contains an anchor, as long as
#'        there are enough anchors.
#' @param dependent Values of `role` that must share a group with an anchor.
#' @param n Population size when `role` is not given.
#'
#' @details
#' The algorithm has four steps:
#'
#'  1. Group sizes are drawn from `size.dist` until they sum to at least `n`;
#'     the last group is trimmed so that the sizes sum to exactly `n`.
#'  2. One anchor is placed in each group, in random order, until either the
#'     groups or the anchors run out. When there are fewer anchors than
#'     groups, the groups without an anchor cannot receive dependents.
#'  3. Dependents are placed in the open slots of anchored groups, each open
#'     slot being equally likely, so that a group receives dependents in
#'     proportion to its remaining size. When the dependents outnumber those
#'     slots, the surplus is added to random anchored groups, which then
#'     exceed their drawn size; this is the only case in which the realized
#'     size distribution departs from `size.dist`, and a message reports it.
#'  4. All remaining nodes (further anchors, and roles that are neither) fill
#'     the remaining open slots in the same way.
#'
#' Without `role`, or without `anchor`, nodes are assigned to the drawn group
#' sizes at random and no composition rule applies.
#'
#' The rule guarantees cross-generational contact for every dependent. It does
#' not otherwise control group composition: the number of dependents per
#' group, the pairing of anchors, and whether a third role (older adults, for
#' example) lives alone or with others all follow from the size distribution
#' and the population shares of the roles. When the composition itself is the
#' quantity to control, build the population from a table of group types with
#' [sample_groups()] instead.
#'
#' @return An integer vector of group ids, one per node, with ids running from
#'   `1` to the number of groups.
#'
#' @seealso [sample_groups()] and [netclique()].
#'
#' @export
#'
#' @examples
#' # Ages from a census-like distribution, then households around them
#' set.seed(1)
#' n <- 1000
#' age <- sample(c("child", "adult", "elderly"), n, replace = TRUE,
#'               prob = c(0.22, 0.60, 0.18))
#' size.dist <- c("1" = 0.28, "2" = 0.34, "3" = 0.15, "4" = 0.14, "5" = 0.09)
#' hh_id <- assign_groups(size.dist, role = age,
#'                        anchor = c("adult", "elderly"), dependent = "child")
#'
#' table(tabulate(hh_id))                          # realized size distribution
#' has_anchor <- tapply(age %in% c("adult", "elderly"), hh_id, any)
#' all(has_anchor[as.character(hh_id[age == "child"])])  # every child has one
#'
#' nw <- network_initialize(n)
#' nw <- set_vertex_attribute(nw, "age", age)
#' nw <- set_vertex_attribute(nw, "hh_id", hh_id)
#' est_hh <- netclique(nw, group.attr = "hh_id")
#' print(est_hh, by = "age")
#'
assign_groups <- function(size.dist, role = NULL, anchor = NULL,
                          dependent = NULL, n = NULL) {

  if (!is.numeric(size.dist) || length(size.dist) == 0 ||
        any(is.na(size.dist)) || any(size.dist < 0) || sum(size.dist) == 0) {
    stop("`size.dist` must be a vector of non-negative weights over group ",
         "sizes.")
  }
  sizes <- if (is.null(names(size.dist))) {
    seq_along(size.dist)
  } else {
    suppressWarnings(as.integer(names(size.dist)))
  }
  if (anyNA(sizes) || any(sizes < 1)) {
    stop("The names of `size.dist` must be positive integer group sizes.")
  }
  prob <- size.dist / sum(size.dist)

  if (!is.null(role)) {
    if (!is.null(n) && n != length(role)) {
      stop("`n` does not match the length of `role`.")
    }
    n <- length(role)
  }
  if (is.null(n) || !is.numeric(n) || length(n) != 1 || n < 1) {
    stop("Give the population size in `n` or through `role`.")
  }
  n <- as.integer(n)
  if (!is.null(dependent) && is.null(anchor)) {
    stop("`dependent` requires `anchor`.")
  }

  ## 1. group sizes summing to exactly n
  drawn <- integer(0)
  total <- 0L
  mean.size <- sum(prob * sizes)
  while (total < n) {
    k <- ceiling(1.2 * (n - total) / mean.size) + 10
    batch <- sizes[sample.int(length(sizes), k, replace = TRUE, prob = prob)]
    drawn <- c(drawn, batch)
    total <- total + sum(batch)
  }
  cs <- cumsum(drawn)
  last <- which(cs >= n)[1]
  drawn <- drawn[seq_len(last)]
  drawn[last] <- drawn[last] - (cs[last] - n)
  G <- length(drawn)

  ## no composition rule: random assignment to the drawn sizes
  if (is.null(role) || is.null(anchor)) {
    group <- integer(n)
    group[sample.int(n)] <- rep(seq_len(G), drawn)
    return(group)
  }

  is.anchor <- role %in% anchor
  is.dep <- if (is.null(dependent)) rep(FALSE, n) else role %in% dependent
  if (any(is.anchor & is.dep)) {
    stop("`anchor` and `dependent` must not share values.")
  }
  ## random order of each role; sample() on a single index k would permute
  ## 1:k instead
  shuffle <- function(x) x[sample.int(length(x))]
  anchors <- shuffle(which(is.anchor))
  deps <- shuffle(which(is.dep))
  others <- shuffle(which(!is.anchor & !is.dep))

  group <- rep(NA_integer_, n)
  open <- drawn

  ## 2. one anchor per group while both last
  n.head <- min(G, length(anchors))
  heads <- anchors[seq_len(n.head)]
  group[heads] <- seq_len(n.head)
  open[seq_len(n.head)] <- open[seq_len(n.head)] - 1L
  anchors <- anchors[-seq_len(n.head)]
  anchored <- seq_len(n.head)

  ## 3. dependents into the open slots of anchored groups
  if (length(deps) > 0) {
    if (length(anchored) == 0) {
      message("No anchors in the population; dependents are placed at random.")
      others <- c(others, deps)
    } else {
      slots <- rep(anchored, open[anchored])
      n.fit <- min(length(deps), length(slots))
      if (n.fit > 0) {
        pick <- slots[sample.int(length(slots), n.fit)]
        group[deps[seq_len(n.fit)]] <- pick
        used <- tabulate(pick, nbins = G)
        open <- open - used
      }
      if (length(deps) > n.fit) {
        surplus <- deps[seq.int(n.fit + 1L, length(deps))]
        message(length(surplus), " dependents exceed the open slots of ",
                "anchored groups and are added to random anchored groups.")
        group[surplus] <- anchored[sample.int(length(anchored),
                                              length(surplus), replace = TRUE)]
      }
    }
  }

  ## 4. everyone else into the remaining open slots
  rest <- c(anchors, others)
  if (length(rest) > 0) {
    slots <- rep(seq_len(G), pmax(open, 0L))
    n.fit <- min(length(rest), length(slots))
    if (n.fit > 0) {
      group[rest[seq_len(n.fit)]] <- slots[sample.int(length(slots), n.fit)]
    }
    if (length(rest) > n.fit) {
      surplus <- rest[seq.int(n.fit + 1L, length(rest))]
      group[surplus] <- sample.int(G, length(surplus), replace = TRUE)
    }
  }

  ## drop ids of groups that ended up empty (possible only after surplus)
  group <- as.integer(factor(group, levels = sort(unique(group))))
  group
}