R/tergm.godfather.R

In statnet/tergm: Fit, Simulate and Diagnose Models for Network Evolution Based on Exponential-Family Random Graph Models

#  File R/tergm.godfather.R in package tergm, part of the
#  Statnet suite of packages for network analysis, https://statnet.org .
#
#  This software is distributed under the GPL-3 license.  It is free,
#  open source, and has the attribution requirements (GPL Section 7) at
#  https://statnet.org/attribution .
#
#  Copyright 2008-2023 Statnet Commons
################################################################################

#' A function to apply a given series of changes to a network.
#' 
#' Gives the network a series of timed proposals it can't refuse. Returns the
#' statistics of the network, and, optionally, the final network.
#' 
#' 
#' @param formula An \code{\link{summary.formula}}-style formula, with
#'   either a \code{\link{network}} or a \code{\link{networkDynamic}}
#'   as the LHS and statistics to be computed on the RHS. If LHS is a
#'   \code{\link{networkDynamic}}, it will be used to derive the
#'   changes to the network whose statistics are wanted. Otherwise,
#'   either \code{changes} or \code{toggles} must be specified, and
#'   the LHS \code{\link{network}} will be used as the starting
#'   network.
#' @param changes A matrix with four columns: time, tail, head, and
#'   new value, describing the changes to be made. Can only be used if
#'   LHS of \code{formula} is not a \code{\link{networkDynamic}}.
#' @param toggles A matrix with three columns: time, tail, and head,
#'   giving the dyads which had changed. Can only be used if LHS of
#'   \code{formula} is not a \code{\link{networkDynamic}}.
#' @param start Time from which to start applying changes.  Note that
#'   the first set of changes will take effect at \code{start + 1}.
#'   Defaults to the time point 1 before the earliest change passed.
#' @param end Time at which to finish applying changes.  Defaults to
#'   the last time point at which a change occurs.
#' @param end.network Whether to return the network that
#'   results. Defaults to \code{FALSE}.
#' @param stats.start Whether to return the network statistics at
#'   \code{start} (before any changes are applied) as the first row of
#'   the statistics matrix.  Defaults to \code{FALSE}, to produce
#'   output similar to that of \code{\link[=simulate.tergm]{simulate}} 
#'   for TERGMs when \code{output="stats"}, where initial network's 
#'   statistics are not returned.
#' @template verbose
#' @param control A control list generated by
#'   \code{\link{control.tergm.godfather}}.
#' @return If \code{end.network} is \code{FALSE} (the default), an
#'   \code{\link{mcmc}} object with the requested network statistics
#'   associated with the network series produced by applying the
#'   specified changes. Its \code{\link{mcmc}} attributes encode the
#'   timing information: so \code{\link{start}(out)} gives the time
#'   point associated with the first row returned, and
#'   \code{\link{end}(out)} out the last. The "thinning interval" is
#'   always 1.
#' 
#'   If \code{end.network} is \code{TRUE}, return a \code{\link{network}} object with
#'   \code{\link{lasttoggle}} "extension", representing the final network, with a
#'   matrix of statistics described in the previous paragraph attached to it as
#'   an \code{attr}-style attribute \code{"stats"}.
#'
#' @seealso [simulate.tergm()], [simulate_formula.network()], [simulate_formula.networkDynamic()]
#'
#' @export tergm.godfather
tergm.godfather <- function(formula, changes=NULL, toggles=changes[,-4,drop=FALSE],
                           start=NULL, end=NULL,
                           end.network=FALSE,
                           stats.start=FALSE,
                           verbose=FALSE,
                           control=control.tergm.godfather()){
  on.exit(ergm_Cstate_clear())

  check.control.class("tergm.godfather", "tergm.godfather")

  nw <- ergm.getnetwork(formula)
  
  formula <- nonsimp_update.formula(formula, nw~., from.new="nw")
  
  if(is.networkDynamic(nw)){
    if(!is.null(toggles)) stop("Network passed already contains change or toggle information.")

    toggles <- do.call(rbind, lapply(nw$mel, function(e) if(length(c(e$atl$active)[is.finite(c(e$atl$active))])) cbind(c(e$atl$active)[is.finite(c(e$atl$active))], e$outl,e$inl) else NULL))
    toggles[,1] <- ceiling(toggles[,1]) # Fractional times take effect at the end of the time step.
    
    net.obs.period<-nw%n%'net.obs.period'
    nwend <- if(!is.null(net.obs.period)) .get.last.obs.time(nw) else NULL
    nwstart<- if(!is.null(net.obs.period)) .get.first.obs.time(nw) else NULL
   
    start <- NVL(start,
                 nwstart,
                 suppressWarnings(min(toggles[,1]))-1
                 )
    if(start==Inf) stop("networkDynamic passed contains no change events or attributes. You must specify start explicitly.")

    end <- NVL(end,
               nwend,
               suppressWarnings(max(toggles[,1]))
               )
    if(end==-Inf) stop("networkDynamic passed contains no change events or attributes. You must specify end explicitly.")

    # The reason why it's > start is that the toggles that took effect
    # at start have already been applied to the network. Conversely,
    # it's <= end because we do "observe" the network at end, so we
    # need to apply the toggles that take effect then.
    toggles <- toggles[toggles[,1]>start & toggles[,1]<=end,,drop=FALSE]

    # This is important, since end is inclusive, but terminus is exclusive.
    if(!all(is.active(nw, onset=start, terminus=end+.Machine$double.eps*end*2, v=seq_len(network.size(nw)), rule="any")
            ==is.active(nw, onset=start, terminus=end+.Machine$double.eps*end*2, v=seq_len(network.size(nw)), rule="all")))
      stop("Network size and/or composition appears to change in the interval between start and end. This is not supported by ergm.godfather() at this time.")

    # Finally, we are ready to extract the network.
    nw <- network.extract.with.lasttoggle(nw, at=start)

  }else{
    if(is.null(toggles)) stop("Either pass a networkDynamic, or provide change or toggle information.")
      
    start <- NVL(start,
                 attr(toggles, "start"),
                 min(toggles[,1])-1
                 )
    end <- NVL(end,
               attr(toggles, "end"),
               max(toggles[,1])
               )

    # The reason why it's > start is that the toggles that took effect
    # at start have already been applied to the network. Conversely,
    # it's <= end because we do "observe" the network at end, so we
    # need to apply the toggles that take effect then.
    toggles <- toggles[toggles[,1]>start & toggles[,1]<=end,,drop=FALSE]

    nw %n% "time" <- start
  }

  if(!is.directed(nw)) {
    tails <- toggles[,2]
    heads <- toggles[,3]
    toggles[,2] <- pmin(tails, heads)
    toggles[,3] <- pmax(tails, heads)
  }
  
  toggles <- toggles[order(toggles[,1],toggles[,2],toggles[,3]),,drop=FALSE]

  formula <- nonsimp_update.formula(formula, nw~., from.new="nw")
  m <- ergm_model(formula, nw, dynamic=TRUE, term.options=control$term.options, extra.aux=list(system=~.lasttoggle))

  state <- ergm_state(nw, model=m)
  m$obs <- summary(state)

  if(verbose) message("Applying changes...")

  z <- .Call("godfather_wrapper",
             state,
             as.integer(nrow(toggles)),
             as.integer(toggles[,1]),
             as.integer(toggles[,2]),
             as.integer(toggles[,3]),
             as.integer(start),
             as.integer(end),
             as.integer(verbose),
             PACKAGE="tergm")

  if(z$status!=0) stop("tergm godfather errored with code ", z$status)

  z$state <- update(z$state)

  stats <- matrix(z$s + m$obs, ncol=nparam(state,canonical=TRUE), byrow=TRUE)
  colnames(stats) <- param_names(m, canonical = TRUE)
  if(!stats.start) stats <- stats[-1,,drop=FALSE]
  #' @importFrom coda mcmc
  stats <- mcmc(stats, start=if(stats.start) start else start+1)
  
  if(end.network){ 
    if(verbose) message("Creating new network...")
    newnetwork <- as.network(z$state)
    attr(newnetwork,"stats")<-stats
    newnetwork
  }else stats
}




####################################################################
# The <control.godfather> function allows for tuning of the
# <ergm.godfather> function
#
# --PARAMETERS--
#   maxedges          : the maximum number of edges to make space
#                       for for the new network; this is ignored
#                       if 5*Clist$nedges is greater; this is also
#                       ignored if 'return_new_network' is FALSE;
#                       default=100000
#
#
# --RETURNED--
#   a list of the above parameters
#
####################################################################
#' Control parameters for [tergm.godfather()].
#'
#' Returns a list of its arguments.
#'
#' @template term_options
#' 
#' @export control.tergm.godfather
control.tergm.godfather <- function(term.options = NULL) {
  control<-list()
  for(arg in names(formals(sys.function())))
    control[arg]<-list(get(arg))
  
  control <- set.control.class("control.tergm.godfather")
  control
}