R/compute_reachability.R

In ISMtools: Interpretive Structural Modelling Analysis Tools

#' Compute Reachability Matrix
#'
#' Calculates the reachability matrix from an adjacency matrix using Warshall's
#' algorithm. The reachability matrix shows all direct and indirect relationships
#' between elements in a system.
#'
#' @param adj_matrix A square adjacency matrix (n x n) containing only 0s and 1s.
#'   A value of 1 at position (i,j) indicates that element i directly influences
#'   element j.
#' @param include_self Logical. If \code{TRUE} (default), the diagonal elements
#'   are set to 1, indicating that each element can reach itself (self-reachability).
#'   This is the standard ISM convention.
#'
#' @return A reachability matrix of the same dimension as the input.
#'   \itemize{
#'     \item A value of 1 at position (i,j) indicates that element i can reach
#'       element j either directly or through intermediate elements.
#'     \item When \code{include_self = TRUE}, diagonal elements are always 1.
#'   }
#'
#' @details
#' The function implements Warshall's algorithm for computing the transitive
#' closure of a directed graph. The time complexity is O(n^3) where n is the
#' number of elements.
#'
#' In standard ISM methodology, the reachability matrix is defined as:
#' \deqn{R = (A + I)^k}
#' where A is the adjacency matrix, I is the identity matrix, and k is the

#' smallest integer such that \eqn{(A + I)^k = (A + I)^{k+1}}.
#'
#' The diagonal elements being 1 (self-reachability) is essential for correct
#' level partitioning in ISM analysis.
#'
#' @references
#' Warfield, J. N. (1974). Developing interconnection matrices in structural
#' modeling. \emph{IEEE Transactions on Systems, Man, and Cybernetics},
#' SMC-4(1), 81-87. \doi{10.1109/TSMC.1974.5408524}
#'
#' @seealso
#' \code{\link{create_relation_matrix}} for creating adjacency matrices,
#' \code{\link{level_partitioning}} for hierarchical decomposition,
#' \code{\link{plot_ism}} for visualization.
#'
#' @export
#' @examples
#' # Create a 4x4 adjacency matrix
#' adj_matrix <- matrix(c(0, 1, 0, 0,
#'                        0, 0, 1, 1,
#'                        0, 0, 0, 0,
#'                        0, 0, 0, 0),
#'                      nrow = 4, byrow = TRUE)
#'
#' # Compute reachability matrix (with self-reachability)
#' reach_matrix <- compute_reachability(adj_matrix)
#' print(reach_matrix)
#'
#' # Note: diagonal elements are 1 (self-reachability)
#' diag(reach_matrix)
#'
#' # With named elements
#' rownames(adj_matrix) <- colnames(adj_matrix) <- c("A", "B", "C", "D")
#' reach_matrix <- compute_reachability(adj_matrix)
#' print(reach_matrix)
compute_reachability <- function(adj_matrix, include_self = TRUE) {
  # Parameter validation
  if (!is.matrix(adj_matrix)) {
    stop("Input must be a matrix", call. = FALSE)
  }
  if (nrow(adj_matrix) != ncol(adj_matrix)) {
    stop("Matrix must be square", call. = FALSE)
  }
  if (any(adj_matrix != 0 & adj_matrix != 1)) {
    stop("Matrix must contain only 0s and 1s", call. = FALSE)
  }

  n <- nrow(adj_matrix)

  # Preserve dimnames
  dn <- dimnames(adj_matrix)

  # Initialize reachability matrix
  # Standard ISM: R = (A + I)^k, so we start with A + I
  reach_matrix <- adj_matrix != 0
  if (include_self) {
    diag(reach_matrix) <- TRUE
  }

  # Warshall's algorithm for transitive closure
  for (k in seq_len(n)) {
    for (i in seq_len(n)) {
      if (reach_matrix[i, k]) {
        reach_matrix[i, ] <- reach_matrix[i, ] | reach_matrix[k, ]
      }
    }
  }

  # Convert back to numeric matrix
  storage.mode(reach_matrix) <- "integer"

  # Restore dimnames
  dimnames(reach_matrix) <- dn

  return(reach_matrix)
}