R/runELA.R

In ELAplus: Energy Landscape Analysis for Ecological Communities

#' Energy landscape analysis in parallel mode
#' @description Identifying stable states and tipping points.
#'
#' This function estimates stable states by repeated steepest-descent searches
#' \code{SSestimate} and then checks tipping points between all pairs of stable
#' states using \code{FindingTippingpoint_cpp}. 
#' The result is summarized and converted by \code{elaconvert}
#' into a list of stable state IDs, their energies, tipping point IDs, and
#' tipping energies.
#'
#' @param sa output of runSA function.
#' @param env vector of environmental factors for specifying a single environmental condition (optional).
#' @param SS.itr number of iterations of the process to identify the stable states.
#' @param FindingTip.itr number of iterations of the process to identify the tipping points.
#' @param threads number of cores (threads) used for calculation.
#' @param reporting enable/disable the reporting function.
#' @param fork logical; if TRUE and \code{threads > 1}, a cluster is created internally for parallel execution.
#'
#' @return A list returned by \code{elaconvert}: (1) stable state IDs, (2) stable state energies,
#'   (3) tipping point ID matrix/list, and (4) tipping energies matrix/list.
#'
#' @useDynLib ELAplus, .registration=TRUE
#' @importFrom Rcpp sourceCpp
#' @importFrom foreach foreach
#' @export
ELA <- function(sa=sa, env=NULL, SS.itr=20000, FindingTip.itr=10000,
                 threads=1, reporting=TRUE, fork=TRUE){
  ######################
  sa2p <- sa2params(sa, env)
  hestp <- sa2p[[1]]
  jestp <- sa2p[[2]]
  gestp <- sa2p[[3]]
  hgestp <- sa2p[[4]]
  if(threads>1&&fork==TRUE){cluster <- parallel::makeCluster(threads, type = if (.Platform$OS.type == "windows") "PSOCK" else "FORK")
      parallel::clusterCall(cluster, function(x) .libPaths(x), .libPaths())
      parallel::clusterSetRNGStream(cluster)
      doParallel::registerDoParallel(cluster)
      on.exit(parallel::stopCluster(cluster), add = TRUE)}
	    
  ####################
  ## ||||||||||||||||||||||||||||||||||||| ##
  if(reporting){cat('Start ELA:\n')}
  start <- proc.time()[3]
	speciesName <- rownames(jestp)

	## ||||||||||||||||||||||||||||||||||||| ##
	## -- Stable state estimatin
	minsets <- as.data.frame(SSestimate(hgestp, jestp, itr = SS.itr))
    uqm <- unique(minsets[,-length(minsets[1,])])
    energy <- apply(uqm,1,function(x){Energy(x,hgestp,jestp)})
    minsets <- cbind(uqm, energy)
    minsets <- as.matrix(minsets[order(minsets[,ncol(minsets)]),])
    colnames(minsets) <- c(speciesName, "energy")
    if(reporting){cat(nrow(minsets),"stable states were found.\n")}
	sampleSS <- data.frame(minsets)
	ssid <- sprintf("SS_%s", formatC(1:nrow(minsets), width = nchar(ncol(minsets)),
	                                 flag = "0"))
	rownames(minsets) <- ssid

	## ||||||||||||||||||||||||||||||||||||| ##
	## -- Tipping point estimation

	if(length(ssid)>1){
	comb <- expand.grid(1:nrow(minsets), 1:nrow(minsets))
	comb <- comb[comb[, 1] <comb[, 2], 2:1]

	tpnodeID <- sprintf("TPnode_%s", formatC(1:nrow(comb), width = nchar(nrow(comb)),
	                             flag = "0"))
  if(reporting){cat('Checking', nrow(comb), 'tipping points.\n')}
    if(threads>1){
	FindTip <- foreach (idx = 1:nrow(comb), .combine=rbind, .packages="ELAplus") %dopar% {

	    minsetsub <- minsets[as.integer(comb[idx, ]), ]
	    ss1 = minsetsub[1, -ncol(minsets)]
	    ss2 = minsetsub[2, -ncol(minsets)]

	    tippoint.tmp = FindingTippingpoint_cpp(s1 = ss1,  s2 = ss2,
	                                           alpha = hgestp, jj = jestp,
	                                           tmax = FindingTip.itr)

	    colnames(tippoint.tmp) <- c(speciesName, "energy")

	    ss <- data.frame(SS1=rownames(minsetsub)[1], SS1.energy=minsetsub[1,ncol(minsetsub)],
	               SS2=rownames(minsetsub)[2], SS2.energy=minsetsub[2,ncol(minsetsub)])
	    tp <- data.frame(TP=tpnodeID[idx], tippoint.tmp)
	    return( cbind(ss, tp))
	}}else{
		FindTip <- foreach (idx = 1:nrow(comb), .combine=rbind) %do% {

	    minsetsub <- minsets[as.integer(comb[idx, ]), ]
	    ss1 = minsetsub[1, -ncol(minsets)]
	    ss2 = minsetsub[2, -ncol(minsets)]

	    tippoint.tmp = FindingTippingpoint_cpp(s1 = ss1,  s2 = ss2,
	                                           alpha = hgestp, jj = jestp,
	                                           tmax = FindingTip.itr)

	    colnames(tippoint.tmp) <- c(speciesName, "energy")

	    ss <- data.frame(SS1=rownames(minsetsub)[1], SS1.energy=minsetsub[1,ncol(minsetsub)],
	               SS2=rownames(minsetsub)[2], SS2.energy=minsetsub[2,ncol(minsetsub)])
	    tp <- data.frame(TP=tpnodeID[idx], tippoint.tmp)
	    return( cbind(ss, tp))
	}}
	}

	## ||||||||||||||||||||||||||||||||||||| ##
	## -- Summarize
	if(length(ssid)>1){
	tpState <- (FindTip[,-c(1:5)]);
    uqt <- unique(tpState[,-length(tpState[1,])])
    energy <- apply(uqt,1,function(x){Energy(x,hgestp,jestp)})
    tpStateUni <- cbind(uqt, energy)
	tpid <- sprintf("TP_%s", formatC(1:nrow(tpStateUni), width = nchar(nrow(tpStateUni)),
	                                 flag = "0"))
	rownames(tpStateUni) <- tpid

	stateInfo <- rbind(data.frame(stateID=rownames(minsets), state=rep('Stable state', nrow(minsets)), minsets),
	                   data.frame(stateID=rownames(tpStateUni), state=rep('Tipping point', nrow(tpStateUni)), tpStateUni))
	                   }
	                   else{
	                   stateInfo <- data.frame(stateID=rownames(minsets), state=rep('Stable state', nrow(minsets)), minsets)
	                   }

	# ||||||||||||||||||||||||||||||||||||| ##
	if(length(ssid)>1){
	Stable <- FindTip[,c(1:4)]

	tmp <- rownames(tpStateUni); names(tmp) <-apply(tpStateUni[,-ncol(tpStateUni)], 1, paste, collapse='')
	ordered <- apply(tpState[,-ncol(tpState)], 1, paste, collapse='')
	tpidVec <- tmp[ordered]

	network <- data.frame(Stable, TP=tpidVec, Tp.energy=FindTip[,ncol(FindTip)])
	elasummary <- list(stateInfo = stateInfo, network =network)
	}else{
	elasummary <- list(stateInfo = stateInfo, network =NULL)
	}

	# ||||||||||||||||||||||||||||||||||||| ##
	if(reporting==TRUE){cat("converting...\n")}
  elaconv <- elaconvert(elasummary)
  end <- proc.time()[3]
	if(reporting==TRUE){cat(sprintf("ELA: elapsed time %.2f sec\n", end - start))}
	return(elaconv)}


#' ELA convert
#' @description Converting output of the ELA function.
#'
#' Converts the raw ELA summary into the standard
#' \code{ela} tuple used in downstream visualization/analysis functions:
#' stable state IDs (integer IDs converted from binary states), their energies,
#' tipping point IDs (with "Inf" for missing pairs), and tipping energies.
#'
#' @param ela output of the ELA function.
#' @return A list: (stablestates, ssenergy, tippingpoints, tpenergy).
#' @importFrom magrittr %>%
#' @noRd
elaconvert <- function(ela){
  stablestates <- apply(ela$stateInfo[ela$stateInfo$state=='Stable state', 3:(length(colnames(ela$stateInfo))-1)],
                        1,
                        bin2id) |> as.vector()

  ssenergy <- ela$stateInfo[ela$stateInfo$state=='Stable state','energy']

  # number of stable state
  num_ss <- nrow(ela$stateInfo[grepl('SS', ela$stateInfo$stateID),])
  num_d <- nchar(ncol(ela$stateInfo)-2)
  get_tp <- function(i, j) {
      ss1 <- sprintf("SS_%s", formatC(i, width = num_d,flag="0"))
      ss2 <- sprintf("SS_%s", formatC(j, width = num_d,flag="0"))
      return(ela$network[ela$network$SS1==ss1 & ela$network$SS2==ss2,
                        "TP"])}

  # tippingpoints
  if(num_ss>1){
  tippingpoints <- lapply(matrix(1:num_ss), function(i)
  apply(matrix(1:num_ss), 1, function(j) {
      tp <- get_tp(i, j)
      if(length(tp)==0){return("Inf")}
      bin <- ela$stateInfo[ela$stateInfo$stateID==tp, 3:(length(colnames(ela$stateInfo))-1)]
                                          return(apply(bin,1, bin2id))
                                          })) %>% 'names<-'(sapply(1:num_ss, \(i) paste(i, sep="")))
  }else{
      tippingpoints <- lapply(matrix(1:1),function(i) apply(matrix(1:1),1, function(j) {
      return("Inf")
      })) %>% 'names<-'(sapply(1:1, \(i) paste(i, sep="")))
  }

  ## tpenergy
  if(num_ss>1){
  tpenergy <- lapply(matrix(1:num_ss),function(i) apply(matrix(1:num_ss),1, function(j) {
      tp <- get_tp(i, j)
      if(length(tp)==0){return(Inf)}
      return(ela$stateInfo[ela$stateInfo$stateID==tp, 'energy'])
      })) %>% 'names<-'(sapply(1:num_ss, \(i) paste(i, sep="")))
  }else{
      tpenergy <- lapply(matrix(1:1),function(i) apply(matrix(1:1),1, function(j) {
      return(Inf)
      })) %>% 'names<-'(sapply(1:1, \(i) paste(i, sep="")))
  }

  return(list(stablestates, ssenergy, tippingpoints, tpenergy))}


#' ELPruning
#' @description Pruning shallow stable states of energy landscape.
#'
#' Iteratively prunes stable states (merges basins). 
#' If \code{type="relative"},  stable state whose energy depth is smaller than \code{th} times the deepest
#' basin depth are pruned and merged to the deeper basins. 
#' If \code{type="bootstrap"}, stable state whose bootstrap confidence interval (default: 95% quantile) of energy gap between tipping point
#' cross zero is pruned. 
#' The function updates stable states and corresponding tipping
#' structures, and returns both the pruned ELA object and a mapping table
#' from stable states before pruning to those after pruning.
#'
#' @param ela output of ELA function.
#' @param th ratio of the depth of the basin to be pruned with respect to the deepest basin.
#' @param type the type of pruning \code{"relative"}(use with \code{th}) or \code{"bootstrap"} (use with \code{lower_qr} and \code{bootstrap_ene}) 
#' @param lower_qr Numeric value. Quantile used for constructing confidence intervals from bootstrap samples to check 
#' whether the estimated energy gaps include zero. Default is \code{lower_qr = 0.05}.
#' @param bootstrap_ene a data frame of bootstrap energies generated by \link{bootstrap_ELA}. Necessary when \code{type="quantile"}.
#' @param threads number of cores (threads) used for calculation.
#' @param reporting enable/disable the reporting function.
#' @param fork logical; if TRUE and \code{threads > 1}, a cluster is created internally for parallel execution.
#'
#' @return A list of two elements: (1) pruned ELA tuple (same format as \link{ELA} output),
#'   and (2) a data frame with columns \code{ss.before.pruning} and \code{ss.after.pruning}.
#'
#' @export
ELPruning <- function(ela, type="relative",th=0.05, lower_qr=0.05, bootstrap_ene = NULL,
                      threads=1, reporting=TRUE, fork=TRUE){
  start <- proc.time()[3]
  if(reporting){cat('Start pruning:\n')}
  
  if(threads>1&&fork==TRUE){cluster <- parallel::makeCluster(threads, type = if (.Platform$OS.type == "windows") "PSOCK" else "FORK")
  parallel::clusterCall(cluster, function(x) .libPaths(x), .libPaths())
  parallel::clusterSetRNGStream(cluster)
  doParallel::registerDoParallel(cluster)
  on.exit(parallel::stopCluster(cluster), add = TRUE)}
  
  tss <- ela[[1]]
  tssen <- ela[[2]]
  tti <- ela[[3]]
  ttien <- ela[[4]]
  if(length(tss)>1){
    rat <- as.data.frame(gtools::permutations(
      n=length(tss), r=2, v=1:length(tss), repeats.allowed = TRUE))
    rat <- rat[rat$V1 < rat$V2, , drop = FALSE]
    tss0 <- tss
    ssrep <- tss
    if(type=="relative"){
      #%>% subset(V1 < V2)
      if(reporting==TRUE){cat("*")}
      cmpax <- calc_mist_pax(rat, ttien, tssen, threads)
      if(reporting==TRUE){cat(".")}
      mist <- as.numeric(unlist(cmpax[,(dim(cmpax)[2]-1):(dim(cmpax)[2])]))
      pax <- cmpax[,-dim(cmpax)[2]:-(dim(cmpax)[2]-1)]
      paxmax <- max(mist)
      tt <- 0
      while(min(pax['val1']) < th * paxmax){
          tt <- tt + 1
          rat <- as.data.frame(gtools::permutations(
              n=length(tss), r=2, v=1:length(tss), repeats.allowed = TRUE))
          #%>% subset(V1 < V2)
          rat <- rat[rat$V1 < rat$V2, , drop = FALSE]
          cmpax <- calc_mist_pax(rat, ttien, tssen, threads)
          if(reporting==TRUE){if(tt%%10 == 0){cat("*")}else{cat(".")}}
          mist <- as.numeric(unlist(cmpax[,(dim(cmpax)[2]-1):(dim(cmpax)[2])]))
          pax <- cmpax[,-dim(cmpax)[2]:-(dim(cmpax)[2]-1)]
          paxmax <- max(mist)
        if(min(pax['val1']) >= th * paxmax){break}else{
            pp <- pax[order(pax$val1),][1, c('val2_1', 'val2_2')] %>% as.list() %>% unlist()
            ssrep <- replace_value(tss[pp[2]], ssrep, tss[pp[1]])
            tss <- tss[tss != tss[pp[2]]]
            tssen <- tssen[tssen != tssen[pp[2]]]
            tti_pre <- tti %>% as.data.frame()
            tti <- tti_pre[rownames(tti_pre) != rownames(tti_pre)[pp[2]],
                           colnames(tti_pre) != colnames(tti_pre)[pp[2]]] %>% as.list()
            ttien_pre <- ttien %>% as.data.frame()
            ttien <- ttien_pre[rownames(ttien_pre) != rownames(ttien_pre)[pp[2]],
                               colnames(ttien_pre) != colnames(ttien_pre)[pp[2]]] %>% as.list()
          }
        if(length(unique(ssrep)) == 1){break}
      }
    }else if (type=="bootstrap"){
      if(is.null(bootstrap_ene)){cat("bootstrap results of ela is missing. Please specify by 'bootstrap_ene='")}
      # Using quantile range in bootstrap energy calculation for pruning
      done <- FALSE
      while(done != TRUE){
        rat <- as.data.frame(gtools::permutations(
          n=length(tss), r=2, v=1:length(tss), repeats.allowed = TRUE))
        rat <- rat[rat$V1 < rat$V2, , drop = FALSE]
        #qr_overlap <- calc_overlap_qr(tss, rat, tssen, tti, ttien, bootstrap_ene, lower_qr, threads)
        qr_overlap <- calc_deltaE_range(tss, rat, tssen, tti, ttien, bootstrap_ene, lower_qr, threads)
        qr_overlap <- qr_overlap[order(qr_overlap$val1),]
        #print(qr_overlap)
        for (i in 1:dim(qr_overlap)[1]) {
          #if (qr_overlap[i,"overlap1"] || qr_overlap[i,"overlap2"]){
          if (qr_overlap[i,"pvalue1"] > lower_qr || qr_overlap[i,"pvalue2"] > lower_qr){
            pp <- qr_overlap[i, c('val2_1', 'val2_2')]%>% as.list() %>% unlist()
            ssrep <- replace_value(tss[pp[2]], ssrep, tss[pp[1]])
            tss <- tss[tss != tss[pp[2]]]
            tssen <- tssen[tssen != tssen[pp[2]]]
            tti_pre <- tti %>% as.data.frame()
            tti <- tti_pre[rownames(tti_pre) != rownames(tti_pre)[pp[2]],
                           colnames(tti_pre) != colnames(tti_pre)[pp[2]]] %>% as.list()
            ttien_pre <- ttien %>% as.data.frame()
            ttien <- ttien_pre[rownames(ttien_pre) != rownames(ttien_pre)[pp[2]],
                               colnames(ttien_pre) != colnames(ttien_pre)[pp[2]]] %>% as.list()
            done <- FALSE
            break
          }
          done <- TRUE
        }
        if(length(unique(ssrep)) == 1){
          done <- TRUE
          break
        }
      }
    }
  tti <- tti %>% as.data.frame() %>% 'colnames<-'(1: length(tti)) %>% as.list()
  ttien <- ttien %>% as.data.frame() %>% 'colnames<-'(1: length(ttien)) %>% as.list()
  end <- proc.time()[3]
  if(reporting==TRUE){cat(sprintf("\nELPruning: elapsed time %.2f sec\n", end - start))}
    states <- as.data.frame(cbind(tss0, ssrep))
    colnames(states) <- c("ss.before.pruning", "ss.after.pruning")
    return(list(list(tss, tssen, tti, ttien), states))
  }else{
    end <- proc.time()[3]
    if(reporting==TRUE){cat(sprintf("\nELPruning: elapsed time %.2f sec\n", end - start))}
    states <- as.data.frame(cbind(tss, tss))
    colnames(states) <- c("ss.before.pruning", "ss.after.pruning")
    return(list(list(tss, tssen, tti, ttien), states))
  }}


#'Calc_overlap_qr
#'@description supporting function for ELpruning with absolute value
#' @param tss A string vector (stable states) 
#' @param rat A data.frame
#' @param tti A matix (tipping points).
#' @param bootstrap_ene A data.frame bootstrap results of estimated energies
#' @param qrange A vector (length: 2) quantile range 
#' @param tssen A numeric vector (energy of stable states).
#' @param ttien A matrix-like object convertible to a numeric matrix.
#' @param threads Integer. If greater than 1, parallel execution will be run.
#' 
#' @return A data frame
#' @importFrom stats quantile
#' @noRd
calc_overlap_qr <- function(tss, rat, tssen, tti, ttien, bootstrap_ene, qrange, threads){
  
  overlap_interval <- function(a_min, a_max, b_min, b_max) {
    v <- c(a_min, a_max, b_min, b_max)
    if (anyNA(v) || any(!is.finite(v))) return(FALSE)
    (a_min <= b_max) && (b_min <= a_max)
  }
  qr_ela <- as.matrix(t(apply(bootstrap_ene, 1, function(x) quantile(x, probs = qrange,na.rm=TRUE))))
  qr_col <- sapply(qrange,function(x){paste0(as.character(x*100),"%")})
  pax <- rat
  npax <- max(pax['V2'])
  tti.m <- as.matrix(tti)
  ttien.m <- as.matrix(ttien)
  ### requires starting clusters out side the function ###
  if(threads>1){a <- foreach(i = pax[,'V1'], j = pax[,'V2'], .combine="rbind", .packages=c("tidyverse","ELAplus")) %dopar% {
    row_ <- as.character(npax*(i-1) + j)
    tipp_id <- tti.m[[j]][i]
    te <- ttien.m[[j]][i]
    tipp_qr <- c(qr_ela[tipp_id,qr_col[1]],qr_ela[tipp_id,qr_col[2]])
    se <- c(tssen[i], tssen[j])
    id1 <- tss[[i]]
    id2 <- tss[[j]]
    # check overlap in the quantile range
    pax[row_, 'val1'] <- min(te - se[1], te - se[2])
    val2_pre <- pax[row_, c('V1', 'V2')] %>% T %>% cbind(se)
    val2_pre <- val2_pre[order(val2_pre$se)][,row_]
    pax[row_, 'val2_1'] <- val2_pre[1]
    pax[row_, 'val2_2'] <- val2_pre[2]
    pax[row_, 'overlap1'] <- overlap_interval(tipp_qr[1],tipp_qr[2],qr_ela[id1,qr_col[1]], qr_ela[id1,qr_col[2]])
    pax[row_, 'overlap2'] <- overlap_interval(tipp_qr[1],tipp_qr[2],qr_ela[id2,qr_col[1]], qr_ela[id2,qr_col[2]])
    pax[row_,]
  }}else{a <- foreach(i = pax[,'V1'], j = pax[,'V2'], .combine="rbind") %do% {
    row_ <- as.character(npax*(i-1) + j)
    tipp_id <- tti.m[[j]][i]
    te <- ttien.m[[j]][i]
    tipp_qr <- c(qr_ela[tipp_id,qr_col[1]],qr_ela[tipp_id,qr_col[2]])
    se <- c(tssen[i], tssen[j])
    id1 <- tss[[i]]
    id2 <- tss[[j]]
    # check overlap in the quantile range
    pax[row_, 'val1'] <- min(te - se[1], te - se[2])
    val2_pre <- pax[row_, c('V1', 'V2')] %>% T %>% cbind(se)
    val2_pre <- val2_pre[order(val2_pre$se)][,row_]
    pax[row_, 'val2_1'] <- val2_pre[1]
    pax[row_, 'val2_2'] <- val2_pre[2]
    pax[row_, 'overlap1'] <- overlap_interval(tipp_qr[1],tipp_qr[2],qr_ela[id1,qr_col[1]], qr_ela[id1,qr_col[2]])
    pax[row_, 'overlap2'] <- overlap_interval(tipp_qr[1],tipp_qr[2],qr_ela[id2,qr_col[1]], qr_ela[id2,qr_col[2]])
    pax[row_,]
  }}
  return(a)}



#'Calc mist pax
#'@description supporting function for ELpruning
#' @param rat A data.frame 
#' @param ttien A matrix-like object convertible to a numeric matrix.
#' @param tssen A numeric vector.
#' @param threads Integer. If greater than 1, parallel execution will be run.
#'
#' @return A data frame
#' @noRd
calc_mist_pax <- function(rat, ttien, tssen, threads){
  pax <- rat
  npax <- max(pax['V2'])
  ttien.m <- as.matrix(ttien)
  ### requires starting clusters out side the function ###
  if(threads>1){a <- foreach(i = pax[,'V1'], j = pax[,'V2'], .combine="rbind", .packages=c("tidyverse","ELAplus")) %dopar% {
    row_ <- as.character(npax*(i-1) + j)
    te <- ttien.m[[j]][i]
    se <- c(tssen[i], tssen[j])
    pax[row_, 'val1'] <- min(te - se[1], te - se[2])
    val2_pre <- pax[row_, c('V1', 'V2')] %>% T %>% cbind(se)
    val2_pre <- val2_pre[order(val2_pre$se)][,row_]
    pax[row_, 'val2_1'] <- val2_pre[1]
    pax[row_, 'val2_2'] <- val2_pre[2]
    pax[row_, 'mist1'] <- te - se[1]
    pax[row_, 'mist2'] <- te - se[2]
    pax[row_,]
    }}else{a <- foreach(i = pax[,'V1'], j = pax[,'V2'], .combine="rbind") %do% {
    row_ <- as.character(npax*(i-1) + j)
    te <- ttien.m[[j]][i]
    se <- c(tssen[i], tssen[j])
    pax[row_, 'val1'] <- min(te - se[1], te - se[2])
    val2_pre <- pax[row_, c('V1', 'V2')] %>% T %>% cbind(se)
    val2_pre <- val2_pre[order(val2_pre$se)][,row_]
    pax[row_, 'val2_1'] <- val2_pre[1]
    pax[row_, 'val2_2'] <- val2_pre[2]
    pax[row_, 'mist1'] <- te - se[1]
    pax[row_, 'mist2'] <- te - se[2]
    pax[row_,]
    }}
  return(a)}


#'Calc_deltaE_range
#'@description supporting function for ELpruning with absolute value
#' @param tss A string vector (stable states) 
#' @param rat A data.frame
#' @param tti A matix (tipping points).
#' @param bootstrap_ene A data.frame bootstrap results of estimated energies
#' @param qrange A numeric value for quantile range (confidence interval).
#' @param tssen A numeric vector (energy of stable states).
#' @param ttien A matrix-like object convertible to a numeric matrix.
#' @param threads Integer. If greater than 1, parallel execution will be run.
#' 
#' @return A data frame
#' @importFrom stats quantile
#' @noRd
calc_deltaE_range <- function(tss, rat, tssen, tti, ttien, bootstrap_ene, qrange, threads){
  
  pax <- rat
  npax <- max(pax['V2'])
  tti.m <- as.matrix(tti)
  ttien.m <- as.matrix(ttien)
  ### requires starting clusters out side the function ###
  if(threads>1){a <- foreach(i = pax[,'V1'], j = pax[,'V2'], .combine="rbind", .packages=c("tidyverse","ELAplus")) %dopar% {
    row_ <- as.character(npax*(i-1) + j)
    tipp_id <- tti.m[[j]][i]
    te <- ttien.m[[j]][i]
    se <- c(tssen[i], tssen[j])
    id1 <- tss[[i]]
    id2 <- tss[[j]]
    # check overlap in the quantile range
    pax[row_, 'val1'] <- min(te - se[1], te - se[2])
    val2_pre <- pax[row_, c('V1', 'V2')] %>% T %>% cbind(se)
    val2_pre <- val2_pre[order(val2_pre$se)][,row_]
    pax[row_, 'val2_1'] <- val2_pre[1]
    pax[row_, 'val2_2'] <- val2_pre[2]
    deltaE1 <- bootstrap_ene[tipp_id,] - bootstrap_ene[id1,]
    deltaE2 <- bootstrap_ene[tipp_id,] - bootstrap_ene[id2,]
    #pax[row_, 'overlap1'] <- as.numeric(quantile(deltaE1,probs=qrange,na.rm=TRUE)) < 0
    pax[row_, 'pvalue1'] <- length(deltaE1[deltaE1 < 0])/length(deltaE1)
    pax[row_, 'pvalue2'] <- length(deltaE1[deltaE2 < 0])/length(deltaE2)
    pax[row_,]
  }}else{a <- foreach(i = pax[,'V1'], j = pax[,'V2'], .combine="rbind") %do% {
    row_ <- as.character(npax*(i-1) + j)
    tipp_id <- tti.m[[j]][i]
    te <- ttien.m[[j]][i]
    se <- c(tssen[i], tssen[j])
    id1 <- tss[[i]]
    id2 <- tss[[j]]
    # check overlap in the quantile range
    pax[row_, 'val1'] <- min(te - se[1], te - se[2])
    val2_pre <- pax[row_, c('V1', 'V2')] %>% T %>% cbind(se)
    val2_pre <- val2_pre[order(val2_pre$se)][,row_]
    pax[row_, 'val2_1'] <- val2_pre[1]
    pax[row_, 'val2_2'] <- val2_pre[2]
    deltaE1 <- bootstrap_ene[tipp_id,] - bootstrap_ene[id1,]
    deltaE2 <- bootstrap_ene[tipp_id,] - bootstrap_ene[id2,]
    pax[row_, 'pvalue1'] <- length(deltaE1[deltaE1 < 0])/length(deltaE1)
    pax[row_, 'pvalue2'] <- length(deltaE1[deltaE2 < 0])/length(deltaE2)
    pax[row_,]
  }}
  return(a)}


#'Replace value
#'@description supporting function for ELpruning
#' @param value_before A value to be replaced.
#' @param list_before A vector (or list) to be copied and modified.
#' @param value_after A value to replace with.
#' @return A modified copy of input with replacements applied.
#' @noRd
replace_value <- function(value_before, list_before, value_after){
  list_after <- list_before
  for(i in 1: length(list_after)){
    if(list_after[i] == value_before){list_after[i] <- value_after}
  }
  return(list_after)
}

#' SteepestDescent
#' @description Perform steepest descent on the energy landscape.
#'
#' Starting from an initial binary state, the steepest descent algorithm 
#' will find a local minimum of the energy landscape. 
#' The computation is delegated to the C++ implementation
#' \code{SteepestDescent_cpp}.
#'
#' @param state Binary vector representing the initial state.
#' @param alpha vector of implicit/explicit environmental effects.
#' @param beta Matrix of species interactions.
#'
#' @return A numeric vector of the reached binary stable state and its energy.
#'
#' @useDynLib ELAplus, .registration=TRUE
#' @export
SteepestDescent <- function(state, alpha, beta){
	res <- SteepestDescent_cpp(state, alpha, beta)
	return(res)
}

#' Bi
#' @description Identify basin (stable state) and its energy for a given state.
#'
#' Applies \link{SteepestDescent} to the input state and returns the basin (stable state) reached,
#' together with its energy. The stable state is encoded as an integer ID using \link{bin2id}.
#'
#' @param xi Binary vector representing the initial state.
#' @param h vector of implicit/explicit environmental effects on species.
#' @param j Matrix of species interactions.
#'
#' @return A list with two elements:
#'   (1) integer stable state ID (basin assignment),
#'   (2) energy of the corresponding stable state.
#' @useDynLib ELAplus, .registration=TRUE
#' @export
Bi <- function(xi, h, j){
  sd <- SteepestDescent(xi,h,j)
  stateyen <- list(as.vector(sd[-length(sd)]), sd[length(sd)])
  state <- stateyen[[1]]
  energy <- stateyen[[2]]
  return(
    list(bin2id(state),
      energy))}


#' Energy
#' @description Calculate energy of a given state.
#'
#' Computing energy in \code{state} using given parameters,\code{alpha} and \code{beta}
#' based on the (extended) maximum pairwise entropy model.
#'  
#' @param state Binary vector representing a state (community composition).
#' @param alpha vector of implicit/explicit environmental effects on species.
#' @param beta Matrix of species interactions.
#'
#' @return A numeric scalar giving the energy of the state.
#'
#' @useDynLib ELAplus, .registration=TRUE
#' @export
Energy <- function(state, alpha, beta){
	res <- cEnergy(state, alpha, beta)
	return(res)
}


#' GradELA
#' @description Estimate energy landscape across environmental gradients.
#'
#' Runs \link{ELA} and \link{ELPruning} comprehensively across specified environmental parameters.
#' One environmental parameter can be changed by specifying \code{eid} over \code{steps} within \code{range}. 
#' If \code{env} is not given, the mean vector of \code{enmat} is used as the reference environmental condition.
#'
#' @param sa output of runSA function.
#' @param eid target environmental factor: the specified environmental factor to be changed.
#' @param enmat array/matrix of environmental factors for each sample (normalization required).
#' @param pruning_type the type of pruning \code{"relative"}(use with \code{th}) or \code{"bootstrap"} (use with \code{lower_qr},\code{ocmat} and \code{bs_params}) 
#' @param env array of environmental factor values. 
#' For the environmental parameters that are not specified for \code{eid}, these values are used as references.
#' If not specified, then the mean of values in \code{enmat} are used as references.
#' @param range range over which the environmental factor specified by eid is changed.
#' @param steps number of steps in the range.
#' @param th ratio of the depth of the basin to be pruned with respect to the deepest basin.
#' @param threads number of cores (threads) used for calculation.
#' @param SS.itr number of iterations to identify stable states.
#' @param FindingTip.itr number of iterations to identify tipping points.
#' @param bs_params output of \link{bootstrap_SA}. list of estimated parameters using bootstrap resampled data set.
#' necessary when \code{type="bootstrap"}.
#' @param ocmat community composition matrix.
#' @param lower_qr Numeric value. Quantile used for constructing confidence intervals from bootstrap samples to check 
#' whether the estimated energy gaps include zero. Default is \code{lower_qr = 0.05},necessary when \code{type="bootstrap"}.
#' @param threads number of cores (threads) used for calculation.
#'
#' @return A list of three elements: (1) ELA results for each step (pruned),
#'   (2) the sequence of environmental values, and (3) the reference env vector.
#'
#' @useDynLib ELAplus, .registration=TRUE
#' @export
GradELA <- function(sa=sa, pruning_type="relative",eid=NULL, enmat=enmat, env=NULL, 
                    range=NULL, steps=16, 
                    th=0.05, threads=1, SS.itr=20000, FindingTip.itr=10000,
                    lower_qr=0.05,bs_params=NULL,ocmat=NULL){
  start <- proc.time()[3]
  
  if(threads>1){cluster <- parallel::makeCluster(threads, type = if (.Platform$OS.type == "windows") "PSOCK" else "FORK")
  parallel::clusterCall(cluster, function(x) .libPaths(x), .libPaths())
  parallel::clusterSetRNGStream(cluster)
  doParallel::registerDoParallel(cluster)
  on.exit(parallel::stopCluster(cluster), add = TRUE)}
  
  if(is.null(eid)){return(cat("eid not specified prease set eid=position/colname\n"))}
  if(is.character(eid)){ei <- which(colnames(enmat)==eid)}else{ei<-eid}
  if(is.null(env)){refe <- apply(enmat, 2, mean)
  cat("env not specified, the mean of enmat is used\n")
  }else{refe <- env}
  if(is.null(range)){mi <- min(enmat[,ei])}else{mi <- range[1]}
  if(is.null(range)){ma <- max(enmat[,ei])}else{ma <- range[2]}
  dm <- (ma - mi)/(steps - 1)
  de <- seq(mi, ma, dm)
  cat("processing(")
  cat(length(de))
  cat(") |")
  els <- foreach(i=de) %dopar% {
      cat("=")
      ee <- replace(refe,ei,i)
      elanp <- ELA(sa, ee, SS.itr=SS.itr, FindingTip.itr=FindingTip.itr, threads=1, reporting=FALSE, fork=FALSE)
      if(pruning_type=="relative"){
        elap <- ELPruning(elanp, type = "relative", th=th, threads=1, reporting=FALSE, fork=FALSE)
      }
      else if(pruning_type=="bootstrap"){
        bootstrap <- bootstrap_ELA(bs_params,ocmat,ela=elanp,enmat=enmat,env=ee,threads=1,reporting=FALSE)
        bootstrap_ene <- bootstrap[[2]]
        elap <- ELPruning(elanp, type = "bootstrap",bootstrap_ene=bootstrap_ene,threads=1,
                          lower_qr=lower_qr,reporting = FALSE,fork=FALSE)
      }
      return(elap)
  }
  cat("|\n")
  end <- proc.time()[3]
  refe[eid] <- NA
  cat(sprintf("Elapsed time %.2f sec\n", end - start))
  return(list(els, de, refe))
}


#' SSentropy
#' @description Estimate stable state entropy for a given state.
#'
#' Evaluating the entropy of transitions among stable states
#' starting from a given community state. 
#'
#' The returned value includes the estimated entropy, the convergence time,
#' and the number of successful convergences within the given iteration limits.
#'
#' @param state Binary vector (or matrix with one row) representing the current state.
#' @param ss Matrix of stable states.
#' @param alpha Vector of implicit or explicit environmental effects.
#' @param beta Matrix of species interactions.
#' @param seitr Number of stochastic search iterations.
#' @param convTime Maximum number of iterations for each search.
#'
#' @return A numeric matrix with three columns:
#'   (1) stable state entropy,
#'   (2) convergence time,
#'   (3) number of converged searches.
#'
#' @useDynLib ELAplus, .registration=TRUE
#' @export
SSentropy <- function(state, ss, alpha, beta, seitr=1000, convTime=10000){ 
  res <- SSentropy_cpp(uoc= state, ss= ss, alpha= alpha, beta=beta, 
                       seitr=seitr, convTime=convTime) 
  return(res) 
  }


#' stability
#' @description Calculate energy gap and stable state entropy as stability indices and returns a dataframe.
#'
#' For each sample state, this function identifies the basin (stable state) by \link{Bi},
#' then computes (i) the energy gap and (ii) stable state entropy
#' estimated by \link{SSentropy}. 
#' If \code{enmat} is provided and \code{sa} was estimated with environmental
#' factors, the basin assignment is computed for each sample-specific \code{env}. 
#'
#' Note: the returned stable state id corresponds to the basin (stable state) the sample is assigned to.
#'
#' @param sa : output of runSA function.
#' @param ocmat : array of the vectors of presence/absence status for each sample. 
#' @param enmat : array of the vectors of environmental factors for each sample. Normalization required.
#' @param threads : number of cores (threads) used for calculation.
#' @param seitr : Number of stochastic search iterations.
#' @param seconvTime : Maximum number of iterations for each search.
#' @param reporting : enable/disable the reporting function. 
#' @return A data frame with stability indices. Columns include at least
#'   \code{energy.gap}, \code{ss.entropy}, \code{e.realize}, \code{e.stable}, \code{state.id}, and \code{stable.state.id}
#'   (and may include \code{conv_time} and \code{convergence} depending on the branch).
#'
#' @useDynLib ELAplus, .registration=TRUE
#' @export

Stability <- function(sa, ocmat, enmat=NULL, seitr=1000, seconvTime=10000, threads=1, reporting=TRUE){
  if(reporting){cat('Start Stability:\n')}
  start <- proc.time()[3]
  sadim <- dim(sa[[1]])
  if(reporting && sadim[2] > sadim[1]+1 && is.null(enmat)){cat("sa obtained with enmat, but enmat is not given.\nassigned a zero vector as the environmental condition for all data points.\n")}
  if(reporting && sadim[2] == sadim[1]+1 && !is.null(enmat)){cat("sa obtained without enmat, but enmat is given. enmat is omitted.\n")}
  
  no <- ncol(ocmat)
  if(threads>1){cluster <- parallel::makeCluster(threads, type = if (.Platform$OS.type == "windows") "PSOCK" else "FORK")
  parallel::clusterCall(cluster, function(x) .libPaths(x), .libPaths())
  parallel::clusterSetRNGStream(cluster)
  doParallel::registerDoParallel(cluster)
  on.exit(parallel::stopCluster(cluster), add = TRUE)}
  
  if(!is.null(enmat) && sadim[2] > sadim[1]+1){
      ne <- ncol(enmat)
      
      if(threads>1){stab <- foreach(x=t(cbind(as.data.frame(ocmat), enmat)), .combine="rbind", .packages=c("ELAplus", "tidyverse", "foreach")) %dopar%{
        env <- x[(no+1):(no+ne)]
        sa2p <- sa2params(sa, env)
        je <- sa2p[[2]]
        hge <- sa2p[[4]]
        state <- bin2id(x[1:no]) 
        # basin #
        bi <- Bi(x[1:no], hge, je)
        ssid <- bi[[1]]
        e.stable <- bi[[2]] 
        # energy gap #
        e.realize <- Energy(x[1:no], hge, je)
        energy.gap <- as.numeric(e.realize - e.stable)
        # ss entropy #
        minsets <- SSestimate(hge, je, itr = 10000)
        minsets <- unique(minsets[,-ncol(minsets)])
        sse <- SSentropy(
          as.matrix(t(as.data.frame(x[1:no]))),
          ss=minsets, alpha= hge, beta=je, seitr=seitr, convTime=seconvTime)
        ss.entropy <- sse[,1]
        conv_time <- sse[,2]
        convergence <- sse[,3]
        cbind(as.data.frame(t(as.matrix(c(energy.gap, ss.entropy, conv_time, convergence, e.realize, e.stable)))),
              as.data.frame(t(as.matrix(c(state, ssid))))) %>% "colnames<-" (c("energy.gap", "ss.entropy","conv_time","convergence","e.realize","e.stable","state.id","stable.state.id"))
      }
      
      }else{
        stab <- foreach(x=t(cbind(as.data.frame(ocmat), enmat)), .combine="rbind") %do% {
        env <- x[(no+1):(no+ne)]
        sa2p <- sa2params(sa, env)
        je <- sa2p[[2]]
        hge <- sa2p[[4]]
        state <- bin2id(x[1:no]) 
        # basin #
        bi <- Bi(x[1:no], hge, je)
        ssid <- bi[[1]]
        e.stable <- bi[[2]] 
        # energy gap #
        e.realize <- Energy(x[1:no], hge, je)
        energy.gap <- as.numeric(e.realize - e.stable)
        # ss entropy #
        minsets <- SSestimate(hge, je, itr = 10000)
        minsets <- unique(minsets[,-ncol(minsets)])
        sse <- SSentropy(
            as.matrix(t(as.data.frame(x[1:no]))),
                         ss=minsets, alpha= hge, beta=je, seitr=seitr, convTime=seconvTime)
        ss.entropy <- sse[,1]
        conv_time <- sse[,2]
        convergence <- sse[,3]
        cbind(as.data.frame(t(as.matrix(c(energy.gap, ss.entropy, conv_time, convergence, e.realize, e.stable)))),
              as.data.frame(t(as.matrix(c(state, ssid))))) %>% "colnames<-" (c("energy.gap", "ss.entropy","conv_time","convergence","e.realize","e.stable","state.id","stable.state.id"))
      }
      }
     
      # ss entropy check #
      if(reporting){cat(sprintf("Convergence rate: %.2f \n",mean(stab$convergence)/seitr))}
      if(mean(stab$convergence)/seitr < 0.9){cat("more than 10% searches were not converged. convTime may be too small, try larger seconvTime.\n")}
      
      rownames(stab) <- rownames(ocmat)
     
     }else{
      sa2p <- sa2params(sa)
      je <- sa2p[[2]]
      hge <- sa2p[[4]]
      sampleSS <- t(apply(ocmat, 1, SteepestDescent, alpha=hge, beta=je))
      minsets <- unique(sampleSS[,-ncol(sampleSS)])
      if(threads>1){
        sse <- foreach(i = 1:dim(ocmat)[1], .combine = 'rbind', .packages = c("ELAplus", "tidyverse", "foreach")) %dopar% {
          SSentropy(t(as.matrix(ocmat[i,])), ss=minsets, alpha= hge, beta=je, seitr=seitr, convTime=seconvTime)
        }
      }
      else{
        sse <- SSentropy(ocmat, ss=minsets, alpha= hge, beta=je, seitr=seitr, convTime=seconvTime)
      }
      if(reporting){cat(sprintf("Convergence rate: %.2f \n",mean(sse[,3])/seitr))}
      if(mean(sse[,3])/seitr < 0.9){cat("more than 10% searches were not converged. convTime may be too small, try larger convTime.\n")}
      ss.entropy <- sse[,1]
      names(ss.entropy) <- rownames(ocmat)
      e.stable <- sampleSS[, ncol(sampleSS)]
      e.realize <- apply(ocmat, 1, Energy, hge, je)
      energy.gap <- e.realize - e.stable
      state.id <- apply(ocmat, 1, bin2id)
      stable.state.id <- apply(sampleSS[, -ncol(sampleSS)], 1, bin2id)
      stab <- cbind(as.data.frame(cbind(energy.gap, ss.entropy, e.realize, e.stable)),
      as.data.frame(cbind(state.id, stable.state.id)))}
  end <- proc.time()[3]
  if(reporting==TRUE){cat(sprintf("Stability: elapsed time %.2f sec\n\n", end - start))}
  return(stab)}

#' gstability
#' @description Calculate energy gap, energy barrier and stable state entropy.
#'
#' This function evaluates stability indices for both (i) the pruned and (ii) the non-pruned
#' energy landscape. For each sample, it assigns the basin by \link{Bi} and computes:
#'  - energy gap (sample energy minus basin energy),
#'  - energy barrier (minimum tipping energy connected to the basin minus basin energy),
#'  - stable state entropy (\code{SSentropy}).
#'
#' For basin connectivity: when only one stable state exists, tipping energy is treated as \code{Inf}
#' and therefore energy barrier is returned as \code{Inf}.
#'
#' @param sa : output of runSA function.
#' @param ocmat : matrix of presence/absence status for each sample. 
#' @param enmat : matrix of environmental factors for each sample. Normalization required.
#' @param th : Numerical (0~1). Ratio of the depth of the basin to be pruned with respect to the deepest basin in each pair.
#' @param pruning_type the type of pruning \code{"relative"}(use with \code{th}) or \code{"bootstrap"} (use with \code{lower_qr} and \code{bootstrap_ene}) 
#' @param lower_qr Numeric value. Quantile used for constructing confidence intervals from bootstrap samples to check 
#' whether the estimated energy gaps include zero. Default is \code{lower_qr = 0.05},necessary when \code{type="bootstrap"}.
#' @param bs_params output of \link{bootstrap_SA}. list of estimated parameters using bootstrap resampled data set,
#' necessary when \code{type="bootstrap"}.
#' @param threads number of cores (threads) used for calculation.#' @param threads : number of cores (threads) used for calculation.
#' @param seitr : Number of stochastic search iterations.
#' @param seconvTime : Maximum number of iterations for each search.
#' @param reporting : enable/disable the reporting function. 
#' 
#' @return A list of four elements: (1) pruned stability data frame, (2) non-pruned stability data frame,
#'   (3) auxiliary object for pruned landscape (parameters + stable states table), and
#'   (4) auxiliary object for sample-wise evaluation (inputs and ELA objects).
#'
#' @useDynLib ELAplus, .registration=TRUE
#' @export
gStability <- function(sa, ocmat, enmat=NULL, seitr=1000, seconvTime=10000, pruning_type="relative",
                       th=0.1,lower_qr=0.05,threads=1, reporting=TRUE,bs_params=NULL){
  if(reporting){cat('Start gStability:\n')}
  start <- proc.time()[3]
  sadim <- dim(sa[[1]])
  if(reporting && sadim[2] > sadim[1]+1 && is.null(enmat)){cat("sa obtained with enmat, but enmat is not given.\nassigned a zero vector as the environmental condition for all data points.\n")}
  if(reporting && sadim[2] == sadim[1]+1 && !is.null(enmat)){cat("sa obtained without enmat, but enmat is given. enmat is omitted.\n")}
  cluster <- parallel::makeCluster(threads, type = if (.Platform$OS.type == "windows") "PSOCK" else "FORK")
  parallel::clusterCall(cluster, function(x) .libPaths(x), .libPaths())
  parallel::clusterSetRNGStream(cluster)
  doParallel::registerDoParallel(cluster)
  on.exit(parallel::stopCluster(cluster), add = TRUE)
  if(!is.null(enmat) && sadim[2] > sadim[1]+1){
    no <- ncol(ocmat)
    mo <- nrow(ocmat)
    ne <- ncol(enmat)
    ### ela
    elall <- foreach(x=t(cbind(as.data.frame(ocmat), enmat)), .combine="rbind", .packages=c("ELAplus", "tidyverse", "foreach", "gtools")) %dopar%{
        env <- x[(no+1):(no+ne)]
        elanp <- ELA(sa, env=env, threads=1, reporting=FALSE, fork=FALSE)
        if(pruning_type=="relative"){
          elap <- ELPruning(elanp, type = "relative", th=th, threads=1, reporting=FALSE, fork=FALSE)
        }
        else if(pruning_type=="bootstrap"){
          if(is.null(bs_params)){stop("'bs_params' is missing. Please get bootstrap SA results by running 'bootstrap_SA()'.")}
          bootstrap <- bootstrap_ELA(bs_params,ocmat,ela=elanp,enmat=enmat,env=env,threads=1,reporting=FALSE)
          bootstrap_ene <- bootstrap[[2]]
          elap <- ELPruning(elanp, type = "bootstrap",bootstrap_ene=bootstrap_ene,threads=1,
                            lower_qr=lower_qr,reporting = FALSE,fork=FALSE)
        }
        #ela <- ELPruning(elanp, th=th, threads=1, reporting=FALSE, fork=FALSE)
        list(elanp, elap)}
    #print("Finish pruning.")
    ### stability evaluation (1)
    stab <- foreach(x=t(cbind(as.data.frame(ocmat), enmat)), ela=elall[,2], elanp=elall[,1], .combine="rbind") %do% {
        env <- x[(no+1):(no+ne)]
        sa2p <- sa2params(sa, env)
        je <- sa2p[[2]]
        hge <- sa2p[[4]]
        state <- bin2id(x[1:no])
        #####
        stablestates <- ela[[1]][[1]]
        stablen <- ela[[1]][[2]]
        tippingpoints <- ela[[1]][[3]]
        tippingen <- ela[[1]][[4]]
        stablestates.np <- elanp[[1]]
        stablen.np <- elanp[[2]]
        tippingpoints.np <- elanp[[3]]
        tippingen.np <- elanp[[4]]
        # basin #
        bi <- Bi(x[1:no], hge, je)
        ssid.np <- bi[[1]]
        ssid <- ela[[2]][,2][ela[[2]][,1]==bi[[1]]]
        e.stable <- stablen[stablestates==ssid]
        e.stable.np <- bi[[2]]
        # energy gap #
        e.realize <- Energy(x[1:no], hge, je)
        energy.gap <- as.numeric(e.realize - e.stable)
        energy.gap.np <- as.numeric(e.realize - e.stable.np)
        # energy barrier #
        if(length(stablestates) > 1){
            e.tipping <- min(unlist(c(as.data.frame(tippingen)[stablestates==ssid,],as.data.frame(tippingen)[,stablestates==ssid])))
            energy.barrier <- e.tipping - e.stable}else{
            e.tipping <- Inf
            energy.barrier <- Inf}
        if(length(stablestates.np) > 1){
            e.tipping.np <- min(unlist(c(as.data.frame(tippingen.np)[stablestates.np==ssid.np,],as.data.frame(tippingen.np)[,stablestates.np==ssid.np])))
            energy.barrier.np <- e.tipping.np - e.stable.np}else{
            e.tipping.np <- Inf
            energy.barrier.np <- Inf}
        # ss entropy #
        minsets <- as.matrix(t(sapply(stablestates, function(x){id2bin(x, no)})))
        sse <- SSentropy(as.matrix(t(as.data.frame(x[1:no]))), ss=minsets, alpha= hge, beta=je, seitr=seitr, convTime=seconvTime)
        ss.entropy <- sse[,1]
        conv_time <- sse[,2]
        convergence <- sse[,3]
        minsets.np <- as.matrix(t(sapply(stablestates.np, function(x){ELAplus::id2bin(x, no)})))
        sse.np <- SSentropy(as.matrix(t(as.data.frame(x[1:no]))), ss=minsets.np, alpha= hge, beta=je, seitr=seitr, convTime=seconvTime)
        ss.entropy.np <- sse.np[,1]
        conv_time.np <- sse.np[,2]
        convergence.np <- sse.np[,3]
        cbind(as.data.frame(t(as.matrix(c(energy.gap, ss.entropy, energy.barrier, e.realize, e.stable, e.tipping, conv_time, convergence)))),
              as.data.frame(t(as.matrix(c(state, ssid)))),
              as.data.frame(t(as.matrix(c(energy.gap.np, ss.entropy.np, energy.barrier.np, e.realize, e.stable.np, e.tipping.np, conv_time.np, convergence.np)))),
              as.data.frame(t(as.matrix(c(state, ssid.np))))) %>% "colnames<-" (
                  c("energy.gap", "ss.entropy", "energy.barrier", "e.realize", "e.stable", "e.tipping",
                    "conv_time","convergence","state.id", "stable.state.id",
                    "energy.gap.np", "ss.entropy.np", "energy.barrier.np", "e.realize", "e.stable.np", "e.tipping.np", 
                    "conv_time.np","convergence.np","state.id.np", "stable.state.id.np"))}
     rownames(stab) <- rownames(ocmat)
     if(reporting){cat(sprintf("Convergence rate: %.2f \n",mean(stab$convergence)/seitr))}
     if(reporting){cat(sprintf("Convergence rate (no pruning): %.2f \n",mean(stab$convergence.np)/seitr))}
     if(mean(stab$convergence)/seitr < 0.9){cat("more than 10% searches were not converged. convTime may be too small, try larger seconvTime.\n")}
     if(mean(stab$convergence.np)/seitr < 0.9){cat("more than 10% searches were not converged(no pruning). convTime may be too small, try larger seconvTime.\n")}

     ### Stability evaluation (2)
     stab2 <- foreach(x=t(cbind(as.data.frame(ocmat), enmat)),  ela=elall[,2], elanp=elall[,1]) %do% {
        env <- x[(no+1):(no+ne)]
        sa2p <- sa2params(sa, env)
        stablestates <- ela[[1]][[1]]
        stablen <- ela[[1]][[2]]
        bin <- as.list(lapply(stablestates, function(x){id2bin(x, ncol(ocmat))}))
        sstable <- as.data.frame(cbind(stablestates,
                                       stablen, t(as.data.frame(bin)))) %>% 'colnames<-'(c('ID', 'Energy', colnames(ocmat))) %>% 'rownames<-'(1: length(stablestates))
        list(list(sa2p, sstable), list(ocmat, env, sa, ela, elanp))}
     stab2 <- list(lapply(stab2, function(x) x[[1]]), lapply(stab2, function(x) x[[2]]))
     }else{
      no <- ncol(ocmat)
      sa2p <- sa2params(sa)
      je <- sa2p[[2]]
      hge <- sa2p[[4]]
      ### ela
      elanp <- ELA(sa, env=NULL, threads=threads, reporting=FALSE, fork=FALSE)
      if(pruning_type=="relative"){
        ela <- ELPruning(elanp, type = "relative", th=th, threads=threads, reporting=FALSE, fork=FALSE)
      }
      else if(pruning_type=="bootstrap"){
        bootstrap <- bootstrap_ELA(bs_params,ocmat,ela=elanp,enmat=enmat,env=NULL,threads=1,reporting=FALSE)
        bootstrap_ene <- bootstrap[[2]]
        ela <- ELPruning(elanp, type = "bootstrap",bootstrap_ene=bootstrap_ene,threads=1,
                          lower_qr=lower_qr,reporting = FALSE,fork=FALSE)
      }
      
      ### Stability evaluation (1)
      stablestates <- ela[[1]][[1]]
      stablen <- ela[[1]][[2]]
      tippingpoints <- ela[[1]][[3]]
      tippingen <- ela[[1]][[4]]
      stablestates.np <- elanp[[1]]
      stablen.np <- elanp[[2]]
      tippingpoints.np <- elanp[[3]]
      tippingen.np <- elanp[[4]]
      stab <- foreach(x=t(as.data.frame(ocmat)), .combine="rbind", .packages=c("ELAplus", "tidyverse", "foreach")) %dopar% {
        state <- bin2id(x[1:no])
        # basin #
        bi <- Bi(x[1:no], hge, je)
        ssid.np <- bi[[1]]
        ssid <- ela[[2]][,2][ela[[2]][,1]==bi[[1]]]
        e.stable <- stablen[stablestates==ssid]
        e.stable.np <- bi[[2]]
        # energy gap #
        e.realize <- Energy(x[1:no], hge, je)
        energy.gap <- as.numeric(e.realize - e.stable)
        energy.gap.np <- as.numeric(e.realize - e.stable.np)
        # energy barrier #
        if(length(stablestates) > 1){
            e.tipping <- min(unlist(c(as.data.frame(tippingen)[stablestates==ssid,],as.data.frame(tippingen)[,stablestates==ssid])))
            energy.barrier <- e.tipping - e.stable}else{
            e.tipping <- Inf
            energy.barrier <- Inf}
        if(length(stablestates.np) > 1){
            e.tipping.np <- min(unlist(c(as.data.frame(tippingen.np)[stablestates.np==ssid.np,],as.data.frame(tippingen.np)[,stablestates.np==ssid.np])))
            energy.barrier.np <- e.tipping.np - e.stable.np}else{
            e.tipping.np <- Inf
            energy.barrier.np <- Inf}
        # ss entropy #
        minsets <- as.matrix(t(sapply(stablestates, function(x){id2bin(x, no)})))
        sse <- SSentropy(as.matrix(t(as.data.frame(x[1:no]))), ss=minsets, alpha= hge, beta=je, seitr=seitr, convTime=seconvTime)
        ss.entropy <- sse[,1]
        conv_time <- sse[,2]
        convergence <- sse[,3]
        minsets.np <- as.matrix(t(sapply(stablestates.np, function(x){ELAplus::id2bin(x, no)})))
        sse.np <- SSentropy(as.matrix(t(as.data.frame(x[1:no]))), ss=minsets.np, alpha= hge, beta=je, seitr=seitr, convTime=seconvTime)
        ss.entropy.np <- sse.np[,1]
        conv_time.np <- sse.np[,2]
        convergence.np <- sse.np[,3]
        data <- cbind(as.data.frame(t(as.matrix(c(energy.gap, ss.entropy, energy.barrier, e.realize, e.stable, e.tipping, conv_time, convergence)))),
              as.data.frame(t(as.matrix(c(state, ssid)))),
              as.data.frame(t(as.matrix(c(energy.gap.np, ss.entropy.np, energy.barrier.np, e.realize, e.stable.np, e.tipping.np, conv_time.np, convergence.np)))),
              as.data.frame(t(as.matrix(c(state, ssid.np))))) %>% "colnames<-" (
                c("energy.gap", "ss.entropy", "energy.barrier", "e.realize", "e.stable", "e.tipping",
                  "conv_time","convergence","state.id", "stable.state.id",
                  "energy.gap.np", "ss.entropy.np", "energy.barrier.np", "e.realize", "e.stable.np", "e.tipping.np", 
                  "conv_time.np","convergence.np","state.id.np", "stable.state.id.np"))}
      rownames(stab) <- rownames(ocmat)
      if(reporting){cat(sprintf("Convergence rate: %.2f \n",mean(stab$convergence)/seitr))}
      if(reporting){cat(sprintf("Convergence rate (no pruning): %.2f \n",mean(stab$convergence.np)/seitr))}
      if(mean(stab$convergence)/seitr < 0.9){cat("more than 10% searches were not converged. convTime may be too small, try larger seconvTime.\n")}
      if(mean(stab$convergence.np)/seitr < 0.9){cat("more than 10% searches were not converged(no pruning). convTime may be too small, try larger seconvTime.\n")}

      ### Stability evaluation (2)
      bin <- as.list(lapply(stablestates, function(x){id2bin(x, ncol(ocmat))}))
      sstable <- as.data.frame(cbind(stablestates,
                                      stablen, t(as.data.frame(bin)))) %>% 'colnames<-'(c('ID', 'Energy', colnames(ocmat))) %>% 'rownames<-'(1: length(stablestates))
      stab2 <- list(list(sa2p, sstable), list(ocmat, NULL, sa, ela, elanp))
     }
  end <- proc.time()[3]
  if(reporting==TRUE){cat(sprintf("gStability: elapsed time %.2f sec\n\n", end - start))}
  return(list(as.data.frame(stab)[1:8], as.data.frame(stab)[9:16], stab2[[1]], stab2[[2]]))
  }


#' sstable
#' @description Summary of stable states.
#'
#' Builds a table of stable states and their energies from an ELA output.
#' The returned data frame contains stable state IDs and the corresponding binary
#' presence/absence vector (per species). This function is typically used for
#' inspecting stable states after pruning (if \code{ela} is pruned) or before pruning.
#'
#' @param ela output of ELA or the first return value of ELPruning.
#' @param ocmat community composition matrix.
#'
#' @return A data frame with columns \code{ID}, \code{Energy}.
#'
#' @useDynLib ELAplus, .registration=TRUE
#' @export
sstable <- function(ela,ocmat){
  stablestates  <- ela[[1]]
  stablen <- ela[[2]]
  tippingpoints <- ela[[3]]
  tippingen <- ela[[4]]
  bin <- as.list(lapply(stablestates, function(x){id2bin(x, ncol(ocmat))}))
  names(bin) <- stablestates
  sstable <- as.data.frame(cbind(stablestates, stablen, t(as.data.frame(bin)))) %>%
    'colnames<-'(c('ID', 'Energy', colnames(ocmat))) %>%
    'rownames<-'(1: length(stablestates))
  return(sstable)}

#' tptable
#' @description Summary of tipping points.
#'
#' Builds a table of tipping points (transition states) and their energies from an ELA output.
#' This function lists only finite tipping points
#' (i.e., entries not equal to \code{Inf}), and returns the corresponding binary state vectors.
#'
#' Two stable states are required at least.
#'
#' @param ela  output of \link{ELA} or the first return value of \link{ELPruning}.
#' @param ocmat community composition matrix.
#'
#' @return A data frame with columns \code{TP}, \code{SS1}, \code{SS2}, \code{Energy}, 
#' and one column per species (binary state vector).
#' If no tipping points exist, returns NULL (invisibly) after printing a message.
#'
#' @useDynLib ELAplus, .registration=TRUE
#' @export
tptable <- function(ela,ocmat){
  stablestates  <- ela[[1]]
  stablen <- ela[[2]]
  tippingpoints <- ela[[3]]
  tippingen <- ela[[4]]
  tpid.df <- as.data.frame(tippingpoints)
  tpid <- unlist(tpid.df[upper.tri(tpid.df, diag = FALSE)], use.names = FALSE)
  tpen.df <- as.data.frame(tippingen)
  tpen <- unlist(tpen.df[upper.tri(tpen.df, diag = FALSE)], use.names = FALSE)
  tpbin <- as.list(lapply(tpid, function(x){id2bin(x, ncol(ocmat))}))
  if(length(tpbin) == 0){
    cat(sprintf('only one stable state is found.\nNo tipping points exist.\n'))
    return()
  }
  tp.ss1 <- stablestates[unlist(lapply(seq_along(tippingpoints), function(i) {
      non_inf_count <- sum(tippingpoints[[i]] != "Inf")
      return(rep(i, non_inf_count))}))]
  tp.ss2 <-stablestates[unlist(lapply(tippingpoints, function(x) which(x != "Inf")))]

  tptable <- as.data.frame(cbind(tpid, tp.ss1, tp.ss2, tpen, t(as.data.frame(tpbin)))) %>%
    'colnames<-'(c('TP','SS1','SS2', 'Energy', colnames(ocmat))) %>%
    'rownames<-'(1: length(tpid))
  return(tptable)
}

####################################
####################################

#'Calculation stability indices.
#'@description Calculating stable state energy, difference between sample's energy and stable state energy and stable state entropy.
#'
#'@param data data is a matrix.
#'@param alpha alpha is a vector which explicit/implicit preference. Use output of runSA function.
#'@param J J is a matrix which shows species preference. Use output of runSA function.
#'@param prun Numeric 0 ~ 1. ratio of the depth of the basin to be pruned with respect to the deepest basin.
#'. If 0, then pruning will be skipped (default = 0).
#'@param seitr serials for searching stable states
#'@param pruning_type the type of pruning \code{"relative"}(use with \code{th}) or \code{"bootstrap"} (use with \code{lower_qr} and \code{bootstrap_ene}) 
#' @param lower_qr Numeric value. Quantile used for constructing confidence intervals from bootstrap samples to check 
#' whether the estimated energy gaps include zero. Default is \code{lower_qr = 0.05},necessary when \code{type="bootstrap"}.
#'@param bootstrap_ene a data frame of bootstrap energies generated by \link{bootstrap_ELA}. Necessary when \code{type="bootstrap"}.
#'@param threads the number of cores (threads) for parallel computing
#'@param convTime total iterations for searching stable points
#'@param threads the number of cores used for parallel computing.
#'@useDynLib ELAplus, .registration=TRUE
#'@noRd
calcStability <- function(data=NULL, alpha, J,
						  prun=0,threads=1,pruning_type="relative",
						  lower_qr=0.05,bs_params=NULL,seitr=1000, convTime=10000){
						  	
    start <- proc.time()[3]
    sampleSS <- t(apply(data, 1, SteepestDescent_cpp, alpha=alpha, beta=J))
   	
   	uniqueSS <- data.frame(unique(sampleSS[,-ncol(sampleSS)]))
	uniqueSS$energy <- NA
	
	for(i in 1:nrow(uniqueSS)){
	    
	    state <- uniqueSS[i,-ncol(uniqueSS)]
	    detect <- apply(sampleSS[,-ncol(sampleSS)], 1, function(x){ all(x==state) })
	    
	    uniqueSS[i, ncol(uniqueSS)] <- mean(sampleSS[detect, ncol(sampleSS)])
	}
	sampleSS <- uniqueSS[order(uniqueSS$energy),]
	
    ## ================================ ##
	if(prun>0 || !is.null(bs_params)){
		## |||||||||||||||||||||||||||| ##
		## pruned landscape
	    sa_dummy = list(cbind(alpha, J),NULL)
	    elasummary <- ELA(sa = sa_dummy, threads = threads,FindingTip.itr = 10000,
	                      SS.itr = 10000)
		#redela <- ELPruning(elasummary, th=prun)
		if(pruning_type=="relative"){
		  redela <- ELPruning(elasummary, type="relative",th=prun, reporting=FALSE, fork=FALSE)
		}
	  else if(pruning_type=="bootstrap"){
	    print("Pruning starts")
	    bootstrap <- bootstrap_ELA(bs_params,data,ela=elasummary,
	                               enmat=NULL,env=NULL,threads=1,reporting=FALSE)
	    print("bootstrap finished")
	    bootstrap_ene <- bootstrap[[2]]
	    redela <- ELPruning(elasummary, type = "bootstrap",bootstrap_ene=bootstrap_ene,threads=1,
	                     lower_qr=lower_qr,reporting = FALSE,fork=FALSE)
	  }
	 print("Pruning Finished")
		stateInfo <- elasummary[[1]]
		ssInfo <- stateInfo[stateInfo$state=='Stable state',]
		
		## |||||||||||||||||||||||||||| ##
		## -- Convert to deeper landscape
		
		sampleSSpruned <- as.matrix(sampleSS)
		for(i in 1:nrow(ssInfo)){
		    
		    state <- ssInfo[i,-c(1, 2, ncol(ssInfo))]
		    detect <- apply(sampleSS[,-ncol(sampleSS)], 1, function(x){ all(x==state) })
		    
		    if(sum(detect)>1){
		        id <- ssInfo[i,1]
		        
		        if( !any(redela$log[,2]==id)){
		            log <- redela$log
		            change <- log[id,2]
		            
		            sampleSSpruned[detect,] <- matrix( unlist(ssInfo[change, -c(1,2)]), 
		            								   ncol=ncol(sampleSSpruned), nrow=sum(detect), 
		            								   byrow=TRUE)
		        }
		        
		    }
		
		}
		sampleSS <- sampleSSpruned
        ## |||||||||||||||||||||||||||| ##
	}
    
    stablestate=t(apply(sampleSS, 1, 
                    function(x){ 
                        sep <- round(length(x)/20)
                        line <- c()
                        for(i in 1:sep){
                            y <- x[-length(x)]
                            bi <- stats::na.omit(y[((1+20*(i-1)):(20*i))])
                            line <- c(line, str2i( paste0(bi, collapse='')))
                            
                        }
                        
                        ssid=paste(line, collapse='-')
                        c(ssid, x[length(x)]) }))
                  
    ssenergy <- sampleSS[, ncol(sampleSS)]
    energy <- apply(data, 1, Energy, alpha= alpha, beta=J)
    
    energy.gap <- energy-ssenergy
    ssent <- SSentropy_cpp(uoc=data, ss=unique(sampleSS[,-ncol(sampleSS)]),
          				alpha= alpha, beta=J, 
          				seitr=seitr, convTime=convTime)
          	  
    return( data.frame(StableState.id=stablestate[,1], SSenergy=ssenergy,
                       energy.gap=energy.gap, SSentropy= ssent[,1]))
                       
    ## ================================ ##
    end <- proc.time()[3]
    cat(sprintf('Stability calculation done.\nElapsed time %.2f sec\n\n', end-start))                   
}


#' Energy landscape analysis with bootstrap resampled SA results
#' @description Running energy landscape analysis with bootstrap sampling of estimated saresults
#'
#' Performs \link{ELA} and \link{ELPruning} using SA fitting results based on 
#' the bootstrap resampled occurrence and environmental data using \link{bootstrap_SA}.
#' 
#' Basically, representative ELA results with original occurrence and environmental data should be given by \code{ela=ela},
#' alternatively the function can select a representative parameter where commonly appearing stable states are most frequent.
#'
#' The function computes bootstrap energy distributions for all stable states and tipping points
#' appearing in the specified pruned landscape. Energies that are not present in a given trial
#' are stored as \code{NaN} in the summary table.
#'
#' @param bs_params output of \link{bootstrap_SA}. list of estimated parameters using bootstrap resampled data set.
#' @param ocmat community composition matrix.
#' @param env  target environmental factor.
#' @param ela representative ELA result (output of \link{ELA}). If \code{NULL}, the function returns the most representative one.
#' @param enmat array of the vectors of environmental factors for each sample. Normalization required.
#' @param threads number of cores (threads) used for calculation.
#' @param reporting enable/disable the reporting function.
#'
#' @return A list with two elements:
#'   (1) representative pruned ELA result (stable states, energies, tipping points, tipping energies),
#'   (2) a data frame of bootstrap energies (rows = state IDs, columns = trials) with \code{NaN} for missing states.
#'
#' @useDynLib ELAplus, .registration=TRUE
#' @export
bootstrap_ELA <- function(bs_params,ocmat,ela=NULL,env=NULL,enmat=NULL,threads=1,reporting=FALSE){
  iterations <- length(bs_params)
  nspecies <- dim(ocmat)[[2]]
  if(reporting==TRUE){cat("Run bootstrap sampling of sa params.\n")}
  if(threads>1){cluster <- parallel::makeCluster(threads, type = if (.Platform$OS.type == "windows") "PSOCK" else "FORK")
  parallel::clusterCall(cluster, function(x) .libPaths(x), .libPaths())
  parallel::clusterSetRNGStream(cluster)
  doParallel::registerDoParallel(cluster)
  on.exit(parallel::stopCluster(cluster), add = TRUE)}
  
  # summarize stable states after pruning
  if(is.null(ela)){
    if(threads>1){
    elaresults <- foreach(i = 1:iterations, 
                          .packages = c("ELAplus","doParallel","magrittr","gtools")) %dopar% {
                            rand_saparams <- list(bs_params[[i]],enmat)
                            # run ELA
                            elanp <- ELA(rand_saparams, env, threads=1, 
                                         FindingTip.itr = 5000, reporting=FALSE, fork=FALSE)
                            return(list(elanp,rand_saparams))
                          }
    }else{
      elaresults <- foreach(i = 1:iterations, 
                            .packages = c("ELAplus","doParallel","magrittr","gtools")) %do% {
                              rand_saparams <- list(bs_params[[i]],enmat)
                              # run ELA
                              elanp <- ELA(rand_saparams, env, threads=1, 
                                           FindingTip.itr = 5000, reporting=FALSE, fork=FALSE)
                              return(list(elanp,rand_saparams))
                            }
    }
    rand_ela_summary <- data.frame(matrix(NaN, nrow = 1, ncol = iterations))
    colnames(rand_ela_summary) <- paste0("trial", 1:iterations)
    rownames(rand_ela_summary) <- "0"
    default_row <- as.list(rep(NaN, ncol(rand_ela_summary)))
    for(ss in 1:length(elaresults)){
      ela_i <- elaresults[[ss]]
      sstates <- ela_i[[1]][[1]]
      ssenergy <- ela_i[[1]][[2]]
      
      for (s in 1:length(sstates)) {
        if (!(sstates[s] %in% rownames(rand_ela_summary))) {
          rand_ela_summary[sstates[s],] <- default_row
          rand_ela_summary[sstates[s],paste0("trial", ss)] <- ssenergy[s]}
        else{
          rand_ela_summary[sstates[s],paste0("trial", ss)] <- ssenergy[s]
        }
      }
    }
    rand_ela_summary <- rand_ela_summary[-1,]
    
    # select the most representative sample
    # (i.e., where the commonly appaeared stable states appeared most frequently)
    count_non_nan <- apply(rand_ela_summary,1,function(row) sum(!is.nan(row)))
    total_ss_score <- sapply(colnames(rand_ela_summary), function(x) 
    {si <- rand_ela_summary[,x] 
    si <- rownames(rand_ela_summary)[!is.nan(si)]
    return(sum(as.numeric(count_non_nan[si])))})
    rep_trial <- which.max(total_ss_score)
    rep_ela <- elaresults[[rep_trial]][[1]] # representative ELA results
    
  }else{
    rep_ela <- ela
  }
  # listing all tipping points and stable states in the representative resampling
  tipp_all <- unlist(rep_ela[[3]], use.names = FALSE)
  tipp_all <- unique(tipp_all[!(tipp_all=="Inf")])
  allpossible <- c(rep_ela[[1]],tipp_all)
  
  # Calculate qunatile distribution of energy
  if(reporting==TRUE){cat("calculate Energy.\n")}
  if(threads>1){
    energy_bootstrap <- foreach(i = 1:iterations,.packages = c("ELAplus","doParallel","magrittr","gtools")) %dopar% {
    rand_saparams <- sa2params(list(bs_params[[i]],enmat),env = env)
    hg_rand <- rand_saparams[[4]]
    J_rand <- rand_saparams[[2]]
    rand_cenergy <- sapply(allpossible, function(x){Energy(id2bin(x,nspecies),hg_rand,J_rand)})
    return(rand_cenergy)
    }
    }
  else{
      energy_bootstrap <- foreach(i = 1:iterations,.packages = c("ELAplus","doParallel","magrittr","gtools")) %do% {
        rand_saparams <- sa2params(list(bs_params[[i]],enmat),env = env)
        hg_rand <- rand_saparams[[4]]
        J_rand <- rand_saparams[[2]]
        rand_cenergy <- sapply(allpossible, function(x){Energy(id2bin(x,nspecies),hg_rand,J_rand)})
        return(rand_cenergy)}
  }
  rand_ene_summary <- data.frame(matrix(NaN, nrow = length(allpossible), ncol = iterations))
  colnames(rand_ene_summary) <- paste0("trial", 1:iterations)
  rownames(rand_ene_summary) <- allpossible
  for(r_ene in 1:length(energy_bootstrap)){
    energy_i <- energy_bootstrap[[r_ene]][allpossible]
    rand_ene_summary[,r_ene] <- energy_i
  }
  return(list(rep_ela,rand_ene_summary))
}