R/vulnerability.R

In brainGraph: Graph Theory Analysis of Brain MRI Data

#' Calculate graph vulnerability
#'
#' This function calculates the \emph{vulnerability} of the vertices of a graph.
#' Here, vulnerability is considered to be the proportional drop in global
#' efficiency when a given vertex is removed from the graph. The vulnerability
#' of the graph is considered the maximum across all vertices.
#'
#' @param g An \code{igraph} graph object
#' @param use.parallel Logical indicating whether or not to use \emph{foreach}
#'   (default: \code{TRUE})
#' @param weighted Logical indicating whether weighted efficiency should be
#'   calculated (default: \code{FALSE})
#' @export
#' @importFrom Matrix rowSums
#' @importFrom foreach getDoParRegistered
#' @importFrom doParallel registerDoParallel
#'
#' @return A numeric vector of length equal to the vertex count of \emph{g}
#'
#' @seealso \code{\link{efficiency}}
#' @author Christopher G. Watson, \email{cgwatson@@bu.edu}
#' @references Latora, V. and Marchiori, M. (2005) Variability and protection of
#'   infrastructure networks. \emph{Physical Review E}, \bold{71}, 015103.
#'   \doi{10.1103/physreve.71.015103}

vulnerability <- function(g, use.parallel=TRUE, weighted=FALSE) {
  i <- NULL
  stopifnot(is_igraph(g))
  if (isTRUE(weighted)) {
    weights <- NULL
    g.attr <- 'E.global.wt'
    e.attr <- 'weight'
  } else {
    weights <- NA
    g.attr <- 'E.global'
    e.attr <- NULL
    weighted <- NULL
  }
  A <- as_adj(g, sparse=FALSE, names=FALSE, attr=e.attr)
  calc <- g.attr %in% graph_attr_names(g)
  E.global <- if (calc) efficiency(g, 'global', weights) else graph_attr(g, g.attr)

  vuln <- rep.int(0, dim(A)[1L])
  verts <- which(rowSums((A > 0) + 0L) != 0L)
  `%d%` <- `%do%`
  if (isTRUE(use.parallel)) {
    if (!getDoParRegistered()) {
      cl <- makeCluster(getOption('bg.ncpus'))
      registerDoParallel(cl)
    }
    `%d%` <- `%dopar%`
  }
  vuln[verts] <- foreach(i=verts, .combine='c') %d% {
    g.tmp <- graph_from_adjacency_matrix(A[-i, -i, drop=FALSE], mode='undirected', weighted=weighted)
    E.global.tmp <- efficiency(g.tmp, 'global', weights)
    1 - (E.global.tmp / E.global)
  }
  return(vuln)
}