R/vulnerability.R

In cograph: Analysis and Visualization of Complex Networks

# =============================================================================
# Node Vulnerability Analysis
# =============================================================================


#' Node Vulnerability
#'
#' Computes the vulnerability of each node, defined as the relative drop in
#' global efficiency when that node is removed from the network.
#'
#' \deqn{V(i) = \frac{E_{global} - E_{global \setminus i}}{E_{global}}}
#'
#' where \eqn{E_{global}} is the global efficiency of the full network and
#' \eqn{E_{global \setminus i}} is the global efficiency after removing node i
#' and all its edges.
#'
#' @param x Network input: matrix, igraph, network, cograph_network, or tna
#'   object.
#' @param directed Logical or NULL. If NULL (default), auto-detect from matrix
#'   symmetry. Set TRUE to force directed, FALSE to force undirected.
#' @param normalized Logical. If TRUE (default), return the proportional drop.
#'   If FALSE, return the raw efficiency difference.
#' @param weighted Logical. If TRUE, honor edge weights when computing
#'   shortest paths (Dijkstra); distance is \code{1/weight^alpha} per the
#'   usual qgraph/tna convention when \code{invert_weights = TRUE}. If FALSE
#'   (default, matches prior behavior), all edges are treated as unit length.
#' @param invert_weights Logical. If TRUE (default) and weights are present,
#'   invert weights to distances via \code{1/weight^alpha} so that higher
#'   weight = shorter path (matches \code{centrality()}'s default). Ignored
#'   when \code{weighted = FALSE}.
#' @param alpha Weight-to-distance exponent (default 1).
#' @param digits Integer or NULL. Round scores to this many decimal places.
#'   Default NULL (no rounding).
#' @param ... Additional arguments passed to \code{\link{to_igraph}}.
#'
#' @return A data frame of class \code{"cograph_vulnerability"} with columns:
#'   \describe{
#'     \item{node}{Node labels.}
#'     \item{vulnerability}{Vulnerability scores, sorted descending.}
#'   }
#'   The original input network and normalization mode are stored as
#'   attributes.
#'
#' @details
#' Global efficiency is defined as:
#'
#' \deqn{E_{global} = \frac{1}{n(n-1)} \sum_{i \neq j} \frac{1}{d(i,j)}}
#'
#' Nodes with high vulnerability are critical to the network's communication
#' efficiency. Removing them causes the greatest drop in global efficiency.
#'
#' \strong{Performance note:} This function computes all-pairs shortest paths
#' once for the full graph and once per node removal, giving O(n) calls to
#' the shortest-path algorithm. A warning is issued for networks with more
#' than 500 nodes.
#'
#' @references
#' Latora, V. & Marchiori, M. (2007). A measure of centrality based on
#' network efficiency. \emph{New Journal of Physics}, 9(6), 188.
#' \doi{10.1088/1367-2630/9/6/188}
#'
#' @seealso \code{\link{network_global_efficiency}}, \code{\link{robustness}},
#'   \code{\link{centrality}}
#'
#' @export
#' @examples
#' # Star network: hub is most vulnerable
#' star <- matrix(c(0,1,1,1, 1,0,0,0, 1,0,0,0, 1,0,0,0), 4, 4)
#' rownames(star) <- colnames(star) <- c("hub", "a", "b", "c")
#' cograph::vulnerability(star)
#'
#' # Complete graph: all nodes equally vulnerable
#' k4 <- matrix(1, 4, 4); diag(k4) <- 0
#' rownames(k4) <- colnames(k4) <- c("A", "B", "C", "D")
#' cograph::vulnerability(k4)
vulnerability <- function(x, directed = NULL, normalized = TRUE,
                          weighted = FALSE, invert_weights = TRUE,
                          alpha = 1, digits = NULL, ...) {

  stopifnot(is.logical(normalized), length(normalized) == 1L)
  stopifnot(is.logical(weighted), length(weighted) == 1L)

  g <- to_igraph(x, directed = directed, ...)
  is_dir <- igraph::is_directed(g)
  n <- igraph::vcount(g)

  # Resolve distance-weights once. weights = NA forces unweighted (matches
  # prior behavior); a numeric vector triggers Dijkstra. Inversion mirrors
  # centrality()'s path-based convention (higher weight = shorter path).
  raw_weights <- igraph::E(g)$weight
  dist_weights <- if (weighted && !is.null(raw_weights)) {
    if (invert_weights) {
      w <- 1 / (raw_weights ^ alpha)
      w[!is.finite(w)] <- .Machine$double.xmax
      w
    } else {
      raw_weights
    }
  } else {
    NA
  }

  labels <- if (!is.null(igraph::V(g)$name)) {
    igraph::V(g)$name
  } else {
    as.character(seq_len(n))
  }

  .make_result <- function(scores) {
    if (!is.null(digits)) scores <- round(scores, digits)
    df <- data.frame(
      node = names(scores),
      vulnerability = unname(scores),
      stringsAsFactors = FALSE
    )
    attr(df, "network") <- x
    attr(df, "normalized") <- normalized
    class(df) <- c("cograph_vulnerability", "data.frame")
    df
  }

  if (n <= 1) {
    return(.make_result(stats::setNames(rep(NA_real_, n), labels)))
  }

  if (n > 500) {
    warning("Network has ", n, " nodes. vulnerability() computes all-pairs ",
            "shortest paths n+1 times, which may be slow (O(n^3) total).",
            call. = FALSE)
  }

  e_global <- .compute_global_efficiency(g, is_dir, dist_weights)

  if (e_global == 0) {
    return(.make_result(stats::setNames(rep(0, n), labels)))
  }

  # Both full and reduced efficiency use the ORIGINAL n*(n-1) denominator
  # per Latora & Marchiori (2007): vulnerability must be non-negative
  orig_denom <- n * (n - 1)
  # Pre-compute the mapping from original-graph edge index to each reduced
  # subgraph's edge order. With per-vertex deletion, igraph reorders edges,
  # so we resolve weights per-iteration via E(g_red)$weight after deletion.
  vuln_scores <- vapply(seq_len(n), function(i) {
    g_red <- igraph::delete_vertices(g, i)
    red_weights <- if (is.numeric(dist_weights)) {
      # pull the aligned weights from the reduced graph
      w <- igraph::E(g_red)$weight
      if (is.null(w)) NA else {
        # The delete_vertices() path preserves $weight but the values are the
        # ORIGINAL edge weights (not our inverted vector). Re-invert here.
        if (weighted && !is.null(w) && invert_weights) {
          iw <- 1 / (w ^ alpha)
          iw[!is.finite(iw)] <- .Machine$double.xmax
          iw
        } else {
          w
        }
      }
    } else {
      NA
    }
    sp_red <- igraph::distances(g_red,
                                mode = if (is_dir) "out" else "all",
                                weights = red_weights)
    diag(sp_red) <- NA
    valid <- is.finite(sp_red) & sp_red > 0
    e_red <- sum(1 / sp_red[valid]) / orig_denom
    e_global - e_red
  }, numeric(1))

  if (normalized) vuln_scores <- vuln_scores / e_global

  scores <- stats::setNames(vuln_scores, labels)
  scores <- sort(scores, decreasing = TRUE)
  .make_result(scores)
}


#' Plot Node Vulnerability
#'
#' @param x A \code{cograph_vulnerability} object.
#' @param top Integer or NULL. Show only top N nodes. Default NULL (all).
#' @param col Bar color. Default \code{"steelblue"}.
#' @param ... Additional arguments passed to \code{\link[graphics]{barplot}}.
#'
#' @return Invisible \code{x}.
#' @method plot cograph_vulnerability
#' @export
#' @examples
#' star <- matrix(c(0,1,1,1, 1,0,0,0, 1,0,0,0, 1,0,0,0), 4, 4)
#' rownames(star) <- colnames(star) <- c("hub", "a", "b", "c")
#' v <- cograph::vulnerability(star)
#' plot(v)
plot.cograph_vulnerability <- function(x, top = NULL, col = "steelblue", ...) {
  scores <- stats::setNames(x$vulnerability, x$node)
  if (!is.null(top)) scores <- utils::head(scores, top)
  norm <- attr(x, "normalized") %||% TRUE
  ylab <- if (norm) "Efficiency Drop (proportion)" else "Efficiency Drop"
  graphics::barplot(scores, las = 2, col = col, border = "white",
                    main = "Node Vulnerability", ylab = ylab, ...)
  graphics::grid(nx = NA, ny = NULL,
                 col = grDevices::adjustcolor("gray50", 0.3), lty = 1)
  invisible(x)
}

#' @method print cograph_vulnerability
#' @export
#' @noRd
print.cograph_vulnerability <- function(x, ...) {
  norm <- attr(x, "normalized") %||% TRUE
  cat(sprintf("Node Vulnerability (%s)\n",
              if (norm) "normalized" else "raw"))
  print.data.frame(x, row.names = FALSE, ...)
  invisible(x)
}


#' Compute global efficiency of an igraph object
#'
#' Lightweight version for vulnerability computation. Uses same formula as
#' \code{network_global_efficiency()} but operates on igraph directly.
#' @param g An igraph object.
#' @param directed Logical.
#' @param weights Distance weights (numeric vector) or NA for unweighted.
#' @return Numeric scalar: global efficiency.
#' @keywords internal
#' @noRd
.compute_global_efficiency <- function(g, directed, weights = NA) {
  n <- igraph::vcount(g)
  if (n <= 1) return(0)
  sp <- igraph::distances(g, mode = if (directed) "out" else "all",
                          weights = weights)
  diag(sp) <- NA
  inv_sp <- 1 / sp
  inv_sp[is.infinite(sp)] <- 0
  sum(inv_sp, na.rm = TRUE) / (n * (n - 1))
}