--- a
+++ b/R/7.stability.R
@@ -0,0 +1,545 @@
+# ========7.stability===========
+
+#' Evaluate the stability of a network
+#'
+#' @param go_ls an igraph object or igraph list.
+#' @param mode "robust_test", "vulnerability", "robustness"
+#' @param partial how much percent vertexes be removed in total (default: 0.5, only for robust_test)
+#' @param step how many nodes be removed each time? (default: 10, only for robust_test)
+#' @param reps simulation number (default: 9)
+#' @param threads threads
+#' @param verbose verbose
+#' @param keystone remove 70%% keystones instead of remove 50%% nodes (default: False, only for robustness)
+#'
+#' @return a data.frame
+#' @export
+#'
+#' @examples
+#' \donttest{
+#' data("c_net")
+#' if (requireNamespace("ggpmisc")) {
+#'   c_net_stability(co_net, mode = "robust_test", step = 20, reps = 9) -> robust_res
+#'   plot(robust_res, index = "Average_degree", mode = 2)
+#' }
+#'
+#' c_net_stability(co_net, mode = "vulnerability") -> vulnerability_res
+#' plot(vulnerability_res)
+#'
+#' robustness(co_net) -> robustness_res
+#' plot(robustness_res)
+#'
+#' module_detect(co_net) -> co_net_modu
+#' zp_analyse(co_net_modu, mode = 2) -> co_net_modu
+#'
+#' c_net_stability(co_net_modu, mode = "robustness", keystone = TRUE) -> robustness_res
+#' plot(robustness_res)
+#' }
+c_net_stability <- function(go_ls, mode = "robust_test", partial = 0.5, step = 10, reps = 9,
+                            threads = 1, verbose = TRUE, keystone = FALSE) {
+  mode <- match.arg(mode, c("robust_test", "vulnerability", "robustness"))
+  if (mode == "robust_test") {
+    res <- robust_test(go_ls, partial = partial, step = step, reps = reps, threads = threads, verbose = verbose)
+  } else if (mode == "vulnerability") {
+    res <- vulnerability(go_ls, threads = threads, verbose = verbose)
+  } else if (mode == "robustness") {
+    res <- robustness(go_ls, keystone = keystone, reps = reps, threads = threads, verbose = verbose)
+  }
+  return(res)
+}
+
+#' Robust_test for a network
+#' @rdname c_net_stability
+#'
+#' @return data.frame (robustness class)
+#' @export
+#'
+robust_test <- function(go_ls, partial = 0.5, step = 10, reps = 9, threads = 1, verbose = TRUE) {
+  if ("igraph" %in% class(go_ls)) {
+    robustness_res <- robust_test_in(go_ls, partial = partial, step = step, reps = reps, threads = threads, verbose = verbose)
+    robustness_res <- data.frame(robustness_res, group = "Network")
+  } else {
+    if (!is.igraph(go_ls[[1]])) stop("No igraph or igraph-list.")
+    robustness_res <- lapply(names(go_ls), \(i){
+      tmp <- robust_test_in(go_ls[[i]], partial = partial, step = step, reps = reps, threads = threads, verbose = verbose)
+      data.frame(tmp, group = i)
+    })
+    robustness_res <- do.call(rbind, robustness_res)
+    robustness_res$group <- factor(robustness_res$group, levels = names(go_ls))
+  }
+  class(robustness_res) <- c("robust", "data.frame")
+  return(robustness_res)
+}
+
+robust_test_in <- function(go, partial = 0.5, step = 10, reps = 9, threads = 1, verbose = TRUE) {
+  i <- NULL
+
+  cal_del <- \(go, partial, step, rep){
+    nodes <- length(igraph::V(go))
+    floor(nodes * partial) -> del_i
+    del_i_indexs <- data.frame()
+    sequ <- seq(0, del_i, step)
+    if (sequ[length(sequ)] < del_i) sequ <- c(sequ, del_i)
+    res <- lapply(sequ, \(i){
+      # remove i nodes in the network
+      remove_node <- sample(1:nodes, i)
+      dp <- igraph::delete.vertices(go, remove_node)
+      dp <- igraph::delete.vertices(dp, igraph::V(dp)[igraph::degree(dp) == 0])
+
+      # calculate network parameters
+      tmp_ind <- (MetaNet::net_par(dp, mode = "n")$n_index)
+      data.frame(tmp_ind, i = i)
+    })
+    del_i_indexs <- do.call(rbind, res)
+
+    # for (i in sequ) {
+    #   # remove i nodes in the network
+    #   remove_node <- sample(1:nodes, i)
+    #   dp <- igraph::delete.vertices(go, remove_node)
+    #   dp=igraph::delete.vertices(dp,igraph::V(dp)[igraph::degree(dp)==0])
+    #
+    #   # calculate network parameters
+    #   tmp_ind=(MetaNet::net_par(dp,mode = "n")$n_index)
+    #   del_i_indexs <- rbind(
+    #     del_i_indexs,
+    #     data.frame(tmp_ind, i = i)
+    #   )
+    # }
+
+    return(data.frame(del_i_indexs, "rep" = rep))
+  }
+
+  # parallel
+  # main function
+  loop <- function(i) {
+    cal_del(go, partial, step, i)
+  }
+  {
+    if (threads > 1) {
+      pcutils::lib_ps("foreach", "doSNOW", "snow", library = FALSE)
+      pb <- utils::txtProgressBar(max = reps, style = 3)
+      if (verbose) {
+        opts <- list(progress = function(n) utils::setTxtProgressBar(pb, n))
+        cl <- snow::makeCluster(threads)
+      } else {
+        opts <- NULL
+      }
+      doSNOW::registerDoSNOW(cl)
+      res <- foreach::`%dopar%`(
+        foreach::foreach(i = 1:reps, .options.snow = opts),
+        loop(i)
+      )
+      snow::stopCluster(cl)
+      gc()
+    } else {
+      res <- lapply(1:reps, loop)
+    }
+  }
+  # simplify method
+  del_i_indexs <- do.call(rbind, res)
+
+  del_i_indexs$i_ratio <- del_i_indexs$i / length(V(go)) %>% round(., 3)
+  class(del_i_indexs) <- c("robust", class(del_i_indexs))
+  return(del_i_indexs)
+}
+
+#' Plot robust
+#'
+#' @param x `robust_test()` result (robust object)
+#' @param mode plot mode, 1~3
+#' @param indexes indexes selected to show
+#' @param use_ratio use the delete nodes ratio rather than nodes number
+#' @param ... additional arguments for \code{\link[pcutils]{group_box}}
+#'
+#' @return a ggplot
+#' @method plot robust
+#' @exportS3Method
+#' @rdname robust_test
+plot.robust <- function(x, indexes = c("Natural_connectivity", "Average_degree"), use_ratio = FALSE, mode = 1, ...) {
+  i <- group <- variable <- value <- se <- eq.label <- adj.rr.label <- Natural_connectivity <- lm <- NULL
+  robust_res <- x
+  xlab <- "Removed_nodes"
+  if (use_ratio) {
+    robust_res$i <- robust_res$i_ratio
+    xlab <- "Removed_nodes_ratio"
+  }
+  robust_res %>%
+    dplyr::select(i, group, !!indexes) %>%
+    reshape2::melt(id.var = c("i", "group")) -> pdat
+
+  pdat %>%
+    dplyr::group_by(i, variable, group) %>%
+    dplyr::summarise(mean = mean(value), sd = sd(value), se = sd / sqrt(length(value))) -> sdd
+
+  if (mode == 1) {
+    p <- ggplot(sdd, aes(x = i, y = mean, col = group)) +
+      geom_line() +
+      geom_errorbar(data = sdd, aes(ymax = mean + se, ymin = mean - se)) +
+      # geom_smooth(se = FALSE,method = "loess",formula = 'y ~ x') +
+      facet_wrap(. ~ variable, scales = "free") +
+      labs(x = xlab, y = NULL) +
+      theme_bw()
+  }
+  if (mode == 2) {
+    lib_ps("ggpmisc", library = FALSE)
+    p <- ggplot(sdd, aes(x = i, y = mean, col = group)) +
+      geom_point(size = 0.2, alpha = 0.4) +
+      geom_smooth(se = FALSE, method = "loess", formula = "y ~ x") +
+      ggpmisc::stat_poly_eq(
+        aes(label = paste(after_stat(eq.label), after_stat(adj.rr.label), sep = "~~~~~")),
+        formula = y ~ x, parse = TRUE, label.x = "right",
+        size = 3,
+      ) +
+      facet_wrap(. ~ variable, scales = "free") +
+      labs(x = xlab, y = NULL) +
+      theme_bw()
+  }
+  if (mode == 3) {
+    robust_res %>% dplyr::select(i, rep, group, `Natural_connectivity`) -> pdat
+
+    pdat %>%
+      # filter(i<0.4)%>%
+      dplyr::group_by(rep, group) %>%
+      dplyr::summarise(slope = coefficients(lm(`Natural_connectivity` ~ i))[2]) -> slope
+
+    p <- pcutils::group_box(slope["slope"], group = "group", metadata = slope, alpha = TRUE, ...) + theme_bw()
+  }
+  return(p)
+}
+
+#' Vulnerability calculation
+#'
+#' @rdname c_net_stability
+#'
+#' @return a vector
+#' @export
+#' @description
+#' \deqn{Vi=\frac{E-Ei}{E}}
+#' E is the global efficiency and Ei is the global efficiency after the removal of the node i and its entire links.
+#'
+vulnerability <- function(go_ls, threads = 1, verbose = TRUE) {
+  if ("igraph" %in% class(go_ls)) {
+    vulnerability_res <- vul_max(go_ls, threads = threads, verbose = verbose)
+  } else {
+    if (!"igraph" %in% class(go_ls[[1]])) stop("No igraph or igraph-list.")
+    vulnerability_res <- lapply(names(go_ls), \(i){
+      vul_max(go_ls[[i]], threads = threads, verbose = verbose)
+    })
+    names(vulnerability_res) <- names(go_ls)
+  }
+  class(vulnerability_res) <- c("vulnerability", class(vulnerability_res))
+  return(vulnerability_res)
+}
+
+vul_max <- function(go, threads = 1, verbose = TRUE) {
+  i <- NULL
+  stopifnot(is.igraph(go))
+  if (is.null(V(go)$name)) V(go)$name <- V(go)
+
+  # parallel
+  # main function
+  loop <- function(i) igraph::global_efficiency(igraph::delete_vertices(go, i))
+  {
+    if (threads > 1) {
+      pcutils::lib_ps("foreach", "doSNOW", "snow", library = FALSE)
+      if (verbose) {
+        pb <- utils::txtProgressBar(max = length(V(go)$name), style = 3)
+        opts <- list(progress = function(n) utils::setTxtProgressBar(pb, n))
+      } else {
+        opts <- NULL
+      }
+      cl <- snow::makeCluster(threads)
+      doSNOW::registerDoSNOW(cl)
+      res <- foreach::`%dopar%`(
+        foreach::foreach(i = V(go)$name, .options.snow = opts),
+        loop(i)
+      )
+      snow::stopCluster(cl)
+      gc()
+    } else {
+      res <- lapply(V(go)$name, loop)
+    }
+  }
+  # simplify method
+  tmpv <- do.call(c, res)
+
+  (igraph::global_efficiency(go) - tmpv) / igraph::global_efficiency(go)
+}
+
+#' Plot vulnerability
+#'
+#' @param x `vulnerability()` result (vulnerability object)
+#' @param ... add
+#'
+#' @return a ggplot
+#' @method plot vulnerability
+#' @exportS3Method
+#' @rdname vulnerability
+plot.vulnerability <- function(x, ...) {
+  group <- NULL
+  vulnerability_res <- x
+  if (!is.list(vulnerability_res)) {
+    vulnerability_res <- list(Network = vulnerability_res)
+  }
+  pdat <- data.frame(
+    group = factor(names(vulnerability_res), levels = names(vulnerability_res)),
+    vulnerability = vapply(vulnerability_res, max, numeric(1))
+  )
+  ggplot(pdat, aes(group, vulnerability)) +
+    geom_col(aes(fill = group), ...) +
+    geom_text(aes(group, vulnerability * 1.05, label = round(vulnerability, 3))) +
+    theme_bw()
+}
+
+
+#' Robustness after remove 50%% nodes or some hubs, need the metanet contains "cor" attribute.
+#'
+#' @rdname c_net_stability
+#'
+#' @export
+robustness <- function(go_ls, keystone = FALSE, reps = 9, threads = 1, verbose = TRUE) {
+  if ("igraph" %in% class(go_ls)) {
+    robustness_res <- robustness_in(go_ls, keystone = keystone, reps = reps, threads = threads, verbose = verbose)
+  } else {
+    if (!"igraph" %in% class(go_ls[[1]])) stop("No igraph or igraph-list.")
+    robustness_res <- lapply(names(go_ls), \(i){
+      tmp <- robustness_in(go_ls[[i]], keystone = keystone, reps = reps, threads = threads, verbose = verbose)
+      data.frame(tmp, group = i)
+    })
+    robustness_res <- do.call(rbind, robustness_res)
+    robustness_res$group <- factor(robustness_res$group, levels = names(go_ls))
+  }
+  class(robustness_res) <- c("robustness", class(robustness_res))
+  return(robustness_res)
+}
+
+robustness_in <- function(go, keystone = FALSE, reps = 9, threads = 1, verbose = TRUE) {
+  roles <- name <- from <- to <- i <- NULL
+
+  nodes <- length(V(go))
+  floor(nodes * 0.5) -> del_i
+  if (keystone) {
+    get_v(go) -> tmp_v
+    if (!"module" %in% colnames(tmp_v)) stop("no modules, please `module_detect()` first")
+    if (!"roles" %in% colnames(tmp_v)) stop("no roles, please `zp_analyse()` first")
+    tmp_v %>%
+      dplyr::filter(roles == "Module hubs") %>%
+      dplyr::pull(name) -> hubs
+    floor(length(hubs) * 0.7) -> del_i
+  }
+
+  # parallel
+  # main function
+  loop <- \(i){
+    if (!keystone) remove_node <- sample(1:nodes, del_i)
+    if (keystone) remove_node <- sample(hubs, del_i)
+
+    dp <- igraph::delete.vertices(go, remove_node)
+    dead_s <- "init"
+    # calculated the abundance-weighted mean interaction strength (wMIS) of nodes,<=0 will dead
+    # repeat until all >0
+    while (length(dead_s) > 0) {
+      get_e(dp) -> edge_list
+      edge_list %>%
+        dplyr::select(from, cor) %>%
+        rbind(., dplyr::select(edge_list, to, cor) %>% dplyr::rename(from = to)) %>%
+        dplyr::group_by(from) %>%
+        dplyr::summarise(w_degree = sum(cor)) -> w_degree
+      w_degree %>%
+        dplyr::filter(w_degree <= 0) %>%
+        dplyr::pull(from) -> dead_s
+      dp <- igraph::delete.vertices(dp, dead_s)
+    }
+    tmp_ind <- (MetaNet::net_par(dp, mode = "n")$n_index)
+    data.frame(tmp_ind, reps = i)
+  }
+  {
+    if (threads > 1) {
+      pcutils::lib_ps("foreach", "doSNOW", "snow", library = FALSE)
+      if (verbose) {
+        pb <- utils::txtProgressBar(max = reps, style = 3)
+        opts <- list(progress = function(n) utils::setTxtProgressBar(pb, n))
+      } else {
+        opts <- NULL
+      }
+      cl <- snow::makeCluster(threads)
+      doSNOW::registerDoSNOW(cl)
+      res <- foreach::`%dopar%`(
+        foreach::foreach(i = 1:reps, .options.snow = opts),
+        loop(i)
+      )
+      snow::stopCluster(cl)
+      gc()
+    } else {
+      res <- lapply(1:reps, loop)
+    }
+  }
+  # simplify method
+  del_i_indexs <- do.call(rbind, res)
+  return(del_i_indexs)
+}
+
+#' Plot robustness
+#'
+#' @param x `robustness()` result (robustness object)
+#' @param indexes indexes selected to show
+#' @param ... additional arguments for \code{\link[pcutils]{group_box}}
+#'
+#' @return a ggplot
+#' @method plot robustness
+#' @exportS3Method
+#' @rdname robustness
+plot.robustness <- function(x, indexes = "Node_number", ...) {
+  robustness_res <- x
+  if (!"group" %in% colnames(robustness_res)) robustness_res$group <- "Network"
+  pcutils::group_box(robustness_res[, indexes, drop = FALSE],
+    group = "group",
+    metadata = robustness_res, ..., facet = FALSE
+  ) + theme_bw()
+}
+
+
+#' Cohesion calculation
+#'
+#' @param otutab otutab
+#' @param reps iteration time
+#' @param threads threads
+#' @param mycor a correlation matrix you want to use, skip the null model build when mycor is not NULL, default: NULL
+#' @param verbose verbose
+#'
+#' @return Cohesion object: a list with two dataframe
+#' @export
+#' @references
+#' Herren, C. M. & McMahon, K. (2017) Cohesion: a method for quantifying the connectivity of microbial communities. doi:10.1038/ismej.2017.91.
+#' @examples
+#' \donttest{
+#' data("otutab", package = "pcutils")
+#' # set reps at least 99 when you run.
+#' Cohesion(otutab[1:50, ], reps = 19) -> cohesion_res
+#' if (requireNamespace("ggpubr")) {
+#'   plot(cohesion_res, group = "Group", metadata = metadata, mode = 1)
+#'   plot(cohesion_res, group = "Group", metadata = metadata, mode = 2)
+#' }
+#' }
+Cohesion <- function(otutab, reps = 200, threads = 1, mycor = NULL, verbose = TRUE) {
+  i <- NULL
+  d <- t(otutab)
+  rel.d <- d / rowSums(d)
+
+  if (is.null(mycor)) {
+    # Create observed correlation matrix
+    cor.mat.true <- cor(rel.d)
+    nc <- ncol(rel.d)
+
+    loop <- \(which.taxon){
+      # create vector to hold correlations from every permutation for each single otu
+      ## perm.cor.vec.mat stands for permuted correlations vector matrix
+      perm.cor.vec.mat <- vector()
+
+      for (i in 1:reps) {
+        # Create empty matrix of same dimension as rel.d
+        perm.rel.d <- matrix(numeric(0), dim(rel.d)[1], dim(rel.d)[2])
+        rownames(perm.rel.d) <- rownames(rel.d)
+        colnames(perm.rel.d) <- colnames(rel.d)
+
+        # For each otu
+        for (j in seq_len(dim(rel.d)[2])) {
+          # Replace the original taxon vector with a permuted taxon vector
+          perm.rel.d[, j] <- sample(rel.d[, j])
+        }
+
+        # Do not randomize focal column
+        perm.rel.d[, which.taxon] <- rel.d[, which.taxon]
+
+        # Calculate correlation matrix of permuted matrix
+        cor.mat.null <- cor(perm.rel.d, perm.rel.d[, which.taxon])
+
+        # For each iteration, save the vector of null matrix correlations between focal taxon and other taxa
+        perm.cor.vec.mat <- cbind(perm.cor.vec.mat, cor.mat.null)
+      }
+
+      # Save the median correlations between the focal taxon and all other taxa
+      apply(perm.cor.vec.mat, 1, median)
+    }
+    {
+      if (threads > 1) {
+        pcutils::lib_ps("foreach", "doSNOW", "snow", library = FALSE)
+        pb <- utils::txtProgressBar(max = reps, style = 3)
+        if (verbose) {
+          opts <- list(progress = function(n) utils::setTxtProgressBar(pb, n))
+          cl <- snow::makeCluster(threads)
+        } else {
+          opts <- NULL
+        }
+        doSNOW::registerDoSNOW(cl)
+        res <- foreach::`%dopar%`(
+          foreach::foreach(i = 1:nc, .options.snow = opts),
+          loop(i)
+        )
+        snow::stopCluster(cl)
+        gc()
+      } else {
+        res <- lapply(1:nc, loop)
+      }
+    }
+
+    simplify2array(res) -> med.tax.cors
+
+    obs.exp.cors.mat <- cor.mat.true - med.tax.cors
+  } else {
+    obs.exp.cors.mat <- mycor
+  }
+
+  diag(obs.exp.cors.mat) <- 0
+
+  # Calculate connectedness by averaging positive and negative observed - expected correlations
+  pn.mean <- \(vector, mode = "p"){
+    if (mode == "p") {
+      pos.vals <- vector[which(vector > 0)]
+    } else if (mode == "n") pos.vals <- vector[which(vector < 0)]
+    p.mean <- mean(pos.vals)
+    if (length(pos.vals) == 0) p.mean <- 0
+    return(p.mean)
+  }
+  connectedness.pos <- apply(obs.exp.cors.mat, 2, pn.mean, mode = "p")
+  connectedness.neg <- apply(obs.exp.cors.mat, 2, pn.mean, mode = "n")
+
+  # Calculate cohesion by multiplying the relative abundance dataset by associated connectedness
+  cohesion.pos <- rel.d %*% connectedness.pos
+  cohesion.neg <- rel.d %*% connectedness.neg
+
+  #### Combine vectors into one list and print
+  output <- list(
+    data.frame(neg = connectedness.neg, pos = connectedness.pos),
+    data.frame(neg = cohesion.neg, pos = cohesion.pos)
+  )
+
+  names(output) <- c("Connectedness", "Cohesion")
+
+  class(output) <- c("cohesion", class(output))
+  return(output)
+}
+
+
+#' Plot cohesion
+#'
+#' @param x `Cohesion()` result (cohesion object)
+#' @param mode plot mode, 1~2
+#' @param group group name in colnames(metadata)
+#' @param metadata metadata
+#' @param ... additional arguments for \code{\link[pcutils]{group_box}} (mode=1) or \code{\link[pcutils]{group_box}} (mode=2)
+#'
+#' @return a ggplot
+#' @method plot cohesion
+#' @exportS3Method
+#' @rdname Cohesion
+plot.cohesion <- function(x, group, metadata, mode = 1, ...) {
+  neg <- pos <- NULL
+  cohesion_res <- x
+  if (mode == 1) p <- pcutils::stackplot(abs(t(cohesion_res$Cohesion)), group = group, metadata = metadata, ...)
+  if (mode == 2) {
+    co <- cohesion_res$Cohesion %>% dplyr::transmute(`neg:pos` = neg / pos)
+    p <- pcutils::group_box(co, group = group, metadata = metadata, p_value2 = TRUE, ...) +
+      ylab("neg:pos cohesion") + theme_bw()
+  }
+  p
+}