R/level_partitioning.R

In ISMtools: Interpretive Structural Modelling Analysis Tools

#' Level Partitioning for ISM
#'
#' Performs hierarchical level decomposition of system elements based on their
#' reachability and antecedent sets. This is a core step in Interpretive
#' Structural Modelling (ISM) analysis.
#'
#' @param reach_matrix A square reachability matrix (n x n) with 0/1 entries,
#'   typically computed using \code{\link{compute_reachability}}. The diagonal
#'   should be 1 (self-reachability).
#'
#' @return An object of class \code{ism_levels}, which is a list containing:
#'   \itemize{
#'     \item Each element is a vector of node indices belonging to that level
#'     \item Level 1 is the top level (outcomes/dependent variables)
#'     \item Higher numbered levels are lower in the hierarchy (drivers/independent variables)
#'     \item Attribute \code{labels}: node names if the input matrix has dimnames
#'   }
#'
#' @details
#' The algorithm implements the standard ISM level partitioning procedure:
#'
#' \enumerate{
#'   \item For each remaining element i, compute:
#'     \itemize{
#'       \item \strong{Reachability set R(i)}: elements that i can reach (within remaining set)
#'       \item \strong{Antecedent set A(i)}: elements that can reach i (within remaining set)
#'     }
#'   \item An element belongs to the current (top) level if:
#'     \eqn{R(i) \cap A(i) = R(i)}, i.e., the reachability set equals the intersection
#'   \item Remove top-level elements and repeat until all elements are assigned
#' }
#'
#' This implementation correctly operates on the \strong{remaining subset} at each
#' iteration, which is essential for correct level assignment.
#'
#' @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}
#'
#' Sage, A. P. (1977). \emph{Interpretive Structural Modeling: Methodology for
#' Large-scale Systems}. McGraw-Hill.
#'
#' @seealso
#' \code{\link{compute_reachability}} for computing reachability matrices,
#' \code{\link{plot_ism}} for visualization,
#' \code{\link{plot_interactive_ism}} for interactive visualization,
#' \code{\link{micmac_analysis}} for MICMAC analysis.
#'
#' @export
#' @examples
#' # Create 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)
#' rownames(adj_matrix) <- colnames(adj_matrix) <- c("A", "B", "C", "D")
#'
#' # Compute reachability matrix
#' reach_mat <- compute_reachability(adj_matrix)
#'
#' # Perform level partitioning
#' levels <- level_partitioning(reach_mat)
#' print(levels)
#'
#' # Access specific levels
#' levels[[1]]  # Top level elements (outcomes)
#' levels[[length(levels)]]  # Bottom level elements (root causes)
level_partitioning <- function(reach_matrix) {
  # Parameter validation
  if (!is.matrix(reach_matrix)) {
    stop("Input must be a matrix", call. = FALSE)
  }
  if (nrow(reach_matrix) != ncol(reach_matrix)) {
    stop("Matrix must be square", call. = FALSE)
  }
  if (!all(reach_matrix %in% c(0, 1))) {
    stop("Matrix must contain only 0s and 1s", call. = FALSE)
  }

  n <- nrow(reach_matrix)

  # Check if diagonal is all 1s (self-reachability)
  if (!all(diag(reach_matrix) == 1)) {
    warning("Diagonal elements should be 1 (self-reachability). ",
            "Consider using compute_reachability() with include_self = TRUE.",
            call. = FALSE)
  }

  # Get labels if available
  labels <- rownames(reach_matrix)
  if (is.null(labels)) {
    labels <- as.character(seq_len(n))
  }

  remaining <- seq_len(n)
  levels <- list()
  max_iterations <- n  # Safety limit

  while (length(remaining) > 0 && max_iterations > 0) {
    # Find top level elements in current remaining subset
    # An element i is at top level if R(i) ∩ A(i) = R(i)
    # i.e., its reachability set (within remaining) equals the intersection

    is_top_level <- vapply(remaining, function(i) {
      # Reachability set: elements in remaining that i can reach
      reach_set <- remaining[reach_matrix[i, remaining] == 1]

      # Antecedent set: elements in remaining that can reach i
      ante_set <- remaining[reach_matrix[remaining, i] == 1]

      # Check if R(i) ∩ A(i) = R(i)
      # Equivalent to: R(i) is a subset of A(i)
      intersection <- intersect(reach_set, ante_set)
      setequal(reach_set, intersection)
    }, logical(1))

    current_level <- remaining[is_top_level]

    if (length(current_level) == 0) {
      # This should not happen with a valid reachability matrix
      # If it does, there might be a cycle or invalid matrix
      warning("Could not identify top level elements. ",
              "Assigning remaining elements to final level. ",
              "Check if the reachability matrix is valid.",
              call. = FALSE)
      levels[[length(levels) + 1]] <- remaining
      break
    }

    levels[[length(levels) + 1]] <- current_level
    remaining <- setdiff(remaining, current_level)
    max_iterations <- max_iterations - 1
  }

  # Create ism_levels object with labels attribute
  structure(levels,
            class = "ism_levels",
            labels = labels)
}

#' Print ISM Levels
#'
#' Print method for objects of class \code{ism_levels} created by
#' \code{\link{level_partitioning}}.
#'
#' @param x An object of class \code{ism_levels}
#' @param ... Additional arguments (ignored)
#'
#' @return Invisibly returns the input object
#'
#' @export
#' @method print ism_levels
print.ism_levels <- function(x, ...) {
  labels <- attr(x, "labels")

  cat("ISM Hierarchical Levels\n")
  cat("=======================\n")
  cat("(Level 1 = Top/Outcomes, Higher levels = Drivers)\n\n")

  for (i in seq_along(x)) {
    indices <- x[[i]]
    if (!is.null(labels)) {
      node_names <- labels[indices]
      cat(sprintf("Level %d: %s\n", i, paste(node_names, collapse = ", ")))
    } else {
      cat(sprintf("Level %d: %s\n", i, paste(indices, collapse = ", ")))
    }
  }

  cat(sprintf("\nTotal: %d levels, %d elements\n",
              length(x), sum(lengths(x))))

  invisible(x)
}

#' Summary of ISM Levels
#'
#' @param object An object of class \code{ism_levels}
#' @param ... Additional arguments (ignored)
#'
#' @return A data frame with level information
#'
#' @export
#' @method summary ism_levels
summary.ism_levels <- function(object, ...) {
  labels <- attr(object, "labels")

  result <- data.frame(
    level = rep(seq_along(object), lengths(object)),
    index = unlist(object),
    stringsAsFactors = FALSE
  )

  if (!is.null(labels)) {
    result$label <- labels[result$index]
  }

  class(result) <- c("summary.ism_levels", "data.frame")
  result
}