R/visualization.R

In ELAplus: Energy Landscape Analysis for Ecological Communities

#'PCplot:Scatter plot of community compositions colored by basin label
#'
#' @description Performs PCA on the community composition matrix (\code{ocmat}) and plots the first
#' two PCs. Points are colored by the stable state (basin) label estimated
#' from \link{Bi}. If \code{pruned = TRUE}, colors correspond to the
#' pruned/merged basin labels in \code{ssrep} (\code{stablestates2}), otherwise to the
#' original labels (\code{stablestates}). When \code{export = TRUE}, the plotting data
#' frame (stable-state labels + PCs + merged labels) is returned.
#'
#' @param ocmat Data frame (or matrix) of community compositions (rows = samples,
#'   columns = species/features).
#' @param sa Parameter set passed to \link{sa2params}.
#' @param env Optional environment object passed to \link{sa2params}.
#' @param ssrep Data frame-like object providing stable state label mapping with
#'   two columns: \code{stablestates} and \code{stablestates2}.
#' @param pruned Logical; if TRUE, use \code{stablestates2} for coloring.
#' @param export Logical; if TRUE, return the data used for plotting.
#' @importFrom dplyr select left_join bind_rows
#' @return NULL invisibly (or a data frame when \code{export = TRUE}).
#' @useDynLib ELAplus, .registration=TRUE
#' @export
PCplot <- function(ocmat, sa=sa, env=NULL, ssrep=NULL, pruned=TRUE, export=FALSE){
  ######################
  sa2p <- sa2params(sa, env)
  hestp <- sa2p[[1]]
  jestp <- sa2p[[2]]
  gestp <- sa2p[[3]]
  hgestp <- sa2p[[4]] 
  ####################
    if(length(ssrep)!=0){
    ssrep_ <- ssrep %>% as.data.frame() %>% 'colnames<-'(
    c('stablestates', 'stablestates2'))
    stablestates <- ssrep_$stablestates
    stablestates2 <- ssrep_$stablestates2
  PCA_ocmatrix <- as.data.frame(stats::prcomp(x=ocmat, scale=TRUE)[5]) %>%
    'colnames<-'(lapply(1: length(names(hgestp)), function(x){
      gsub(' ', '', paste('PC', as.character(x)), fixed=TRUE)}))
  states <- as.data.frame(t(as.data.frame(
    do.call(
      rbind,
      apply(ocmat, 1, function(x){
        t(Bi(x,
             hgestp, jestp))}))[, 1]))) %>% 'rownames<-'(
               rownames(ocmat)) %>% 'colnames<-'('stablestates')
  state_PC <- cbind(states, PCA_ocmatrix) %>% left_join(ssrep_, by='stablestates')
  if(pruned==FALSE){
    g <- ggplot(state_PC, aes(x=.data$PC1, y=.data$PC2, color=as.factor(stablestates)))}
  else{
    g <- ggplot(state_PC, aes(x=.data$PC1, y=.data$PC2, color=as.factor(stablestates2)))
  }
  g <- g + geom_point() + theme(aspect.ratio = 1)
  plot(g)
  if(export){return(state_PC)}else{return(NULL)}
    }else{
        cat("Please set ssrep")
    }}


#'temp.mat: supporting function for calc_ccc
#'
#' @description Builds an adjacency (0/1) matrix from an edge list \code{rips} with columns \code{V1}
#' and \code{V2}. The matrix size is determined by \code{max(rips)}, and entries
#' \code{matrix.p[V1, V2]} are set to 1.
#'
#' @param rips Edge list as a data frame with at least columns \code{V1} and \code{V2}.
#' @return A square adjacency matrix with dimnames \code{1:n}.
#' @useDynLib ELAplus, .registration=TRUE
#' @noRd
temp.mat <- function(rips){
  n <- max(rips)
  matrix.p <- matrix(0, nrow = n, ncol = n)
  for (i in 1:nrow(rips)) {
    matrix.p[rips$V1[i], rips$V2[i]] <- 1
  }
  dimnames(matrix.p) <- list(1:n, 1:n)
  return(matrix.p)}


#' calc_ccc: supporting function for DisconnectivityGraph
#'
#' @description Computes connected components (clusters) under an energy threshold \code{ww}.
#' Stable states with energy <= \code{ww} are included, and tipping connections with
#' energy <= \code{ww} are added as edges. The resulting graph is
#' decomposed into connected components, filtered to keep components of size >= 2
#' and (optionally) singleton components that correspond to a stable state.
#'
#' @param ww Numeric energy threshold.
#' @param stablestates List (or vector) of stable states.
#' @param ssenergies Energies for stable states.
#' @param tippingenergies List of energies in tipping points.
#' 
#' @importFrom dplyr group_by group_map filter
#' @importFrom igraph components graph_from_adjacency_matrix
#' @importFrom tidyr pivot_longer
#' @return A list of connected components (each component is a vector of node indices).
#' @useDynLib ELAplus, .registration=TRUE
#' @noRd
calc_ccc <- function(ww, stablestates, ssenergies, tippingenergies){
  sta <- data.frame(ssenergies) %>% subset(ssenergies <= ww) %>% rownames() %>%
    as.integer
  coordinates <- as.data.frame(gtools::permutations(n=length(
    stablestates), r=2, v=1:length(stablestates), repeats.allowed = TRUE))
  tip <- cbind(1: length(tippingenergies), as.data.frame(tippingenergies)) %>%
    #pivot_longer(cols=colnames(as.data.frame(tippingenergies))) %>% cbind(
    #  coordinates) %>% subset(value <= ww) %>% select(c('V1', 'V2'))
    pivot_longer(cols=colnames(as.data.frame(tippingenergies))) %>%
    dplyr::bind_cols(coordinates) %>%
    dplyr::filter(.data$value <= ww) %>%
    dplyr::select(c('V1', 'V2'))
  rips <- rbind(as.data.frame(cbind(sta, sta)) %>% 'colnames<-'(c('V1', 'V2')), tip)
  tempmat <- temp.mat(rips)
  ccc_pre1 <- as.vector(components(graph_from_adjacency_matrix(
    adjmatrix=as.matrix(tempmat) + t(tempmat), mode='undirected', diag=F))$membership)
  ccc_pre2 <- T(as.data.frame(rbind(ccc_pre1, 1: length(ccc_pre1)))) %>% 
    group_by(ccc_pre1) %>% group_map(as.list) %>% lapply(function(x){x$V2})
  ccc_pre2 <- ccc_pre2 %>% 'names<-'(1: length(ccc_pre2))
  ccc_pre2_length <- T(as.data.frame(lapply(
    ccc_pre2, function(x){length(x)}))) %>% 'rownames<-'(1: length(ccc_pre2))
  #ccc_final1 <- ccc_pre2[ccc_pre2_length %>% subset(V1 >= 2) %>% rownames]
  #ccc_final2_cond1 <- ccc_pre2_length %>% subset(V1 == 1) %>% rownames
  ccc_final1 <- ccc_pre2[rownames(ccc_pre2_length)[ccc_pre2_length[, "V1"] >= 2]]
  ccc_final2_cond1 <- rownames(ccc_pre2_length)[ccc_pre2_length[, "V1"] == 1]
  #print(ccc_final2_cond1)
  ccc_final2_cond2 <- cbind(ccc_pre2[ccc_final2_cond1], lapply(
    ccc_pre2[ccc_final2_cond1], function(x){length(intersect(x, sta))})) %>% 
    'colnames<-'(c('node', 'len')) %>% applymap(function(x){x[1]}) %>%
    as.data.frame
  ####################################################
  if(length(ccc_final2_cond2)>0){
      ccc_final2 <- ccc_pre2[rownames(ccc_final2_cond2['len'] == 1)]
      }else{ccc_final2 <- NULL}
  ####################################################
  ccc <- c(ccc_final1, ccc_final2)
  return(ccc)
}

#' Assort: supporting function for DisconnectivityGraph
#'
#' @description Sorts (assorts) a list of daughter sets according to their first matching
#' parent set, then by descending set size, then by the first element.
#'
#' @param parent List of parent sets.
#' @param daughter List of daughter sets.
#' @return A one-element list containing the reordered daughter sets.
#' @noRd
Assort <- function(parent, daughter){
  parent_ <- parent %>% unlist(recursive=FALSE)
  daughter_ <- daughter %>% unlist(recursive=FALSE)
  orderkey <- do.call(
    rbind, 
    daughter_  %>%
      lapply(function(x){list(
        match(TRUE, lapply(parent_, function(y){all(x %in% y)})), -length(x), x[1]
        )})) %>% as.data.frame %>% applymap(function(x){x[1]}) %>% as.data.frame %>%
    'colnames<-'(c('setloc', 'length', 'firstN')) %>% 
    cbind(index=1: length(daughter_))
  return(list(daughter_[order(orderkey$setloc, orderkey$length, orderkey$firstN)]))
}

#' calc_cee: supporting function for DisconnectivityGraph
#'
#' @description Iteratively constructs an ordered list of connected component sets
#' from hierarchical edge/component information, using \code{Assort} to keep
#' a consistent parent-daughter ordering.
#'
#' @param hee List of component sets (hierarchical sequence).
#' @return A list of ordered component sets.
#' @noRd
calc_cee <- function(hee){
  cee <- list()
  for(daughter in hee){
    daughter_ <- list(daughter)
    if(length(cee) == 0){
      cee <- append(cee, daughter_)
    }else{
      parent <- cee[length(cee)]
      cee <- append(cee, Assort(parent, daughter_))
    }
  }
  return(cee)
}

#' calc_xpos1: supporting function for DisconnectivityGraph
#'
#' @description Updates x-positions for components when moving from a previous level
#' to a refined partition. If the partition size is
#' unchanged, positions are kept; otherwise, new positions are computed as the
#' average of member positions.
#'
#' @param xpos_prev Previous x-position list.
#' @param rcee_ Current level component sets.
#' @importFrom tidyr all_of
#' @importFrom data.table transpose
#' @return A list containing component ids and their updated x-positions.
#' @noRd
calc_xpos1 <- function(xpos_prev, rcee_){
  rcee__ <- rcee_
  xpos_prev_ <- unlist(xpos_prev, recursive=FALSE)
  if(length(rcee__) == length(unlist(xpos_prev_[1], recursive=FALSE))){
    output <- xpos_prev
  }else{
    xx <- xpos_prev_
    yy <- rcee__
    # pp_pre <- lapply(yy, "length<-", max(lengths(yy))) %>% as.data.frame %>% T %>%
    #  'colnames<-'(1: max(unlist(lapply(yy, length)))) %>%
    #  cbind(index=1: length(yy)) %>%
    #  pivot_longer(cols=-index, names_to='columns', values_to='yy_') %>%
    #  stats::na.omit() %>% applymap(as.integer) %>% as.data.frame() %>%
    #  column_to_rownames(var='yy_')
    tmp <- lapply(yy, "length<-", max(lengths(yy))) %>% as.data.frame %>% T %>%
      'colnames<-'(1: max(unlist(lapply(yy, length)))) %>%
      cbind(index=1: length(yy))
    pp_pre <- pivot_longer(tmp,cols = setdiff(names(tmp),"index"),
                           names_to = "columns",values_to = "yy_") %>%
      stats::na.omit() %>% applymap(as.integer) %>% as.data.frame()
    
    rownames(pp_pre) <- pp_pre$yy_
    pp_pre$yy_ <- NULL
    
    pp <- pp_pre[as.character(unlist(xx[1])), 'index']
    pp2 <- cbind(pp, unlist(xx[2])) %>% as.data.frame %>% 'colnames<-'(c('i', 'j'))
    pp3 <- unique(pp) %>% map(function(x){pp2[pp2$i == x,]})
    output <- list(list(
      yy %>% map(function(x){x[1]}),
      pp3 %>% lapply(function(x){
        sum(unlist(do.call(cbind, x)[, 'j'])) / length(unlist(x$j))})))
  }
  return(output)
}

#' calc_xpos: supporting function for DisconnectivityGraph
#'
#'
#' @param rcee List of component partitions across energy levels.
#' @return A list of x-position specifications for each level.
#' @noRd
calc_xpos <- function(rcee){
  xpositions <- list()
  for(rcee_ in rcee){
    if(length(xpositions) == 0){
      xpositions <- append(
        xpositions, 
        list(list(rcee_ %>% lapply(function(x){x[1]}), 1: length(rcee_))))
    }else{
      xpos_prev <- xpositions[length(xpositions)]
      xpositions <- append(xpositions, calc_xpos1(xpos_prev, rcee_))
    }
  }
  return(xpositions)
}

#' GraphObj: Generating graph objects for DisconnectivityGraph
#'
#' @description Creates a data frame describing nodes/branches for a disconnectivity graph.
#' From ELA output (\code{ela} = stable states, their energies, tipping points, and
#' tipping energies), this function:
#'  - constructs connectivity under decreasing energy thresholds,
#'  - derives component hierarchy,
#'  - assigns x-positions to components for visualization,
#'  - maps energies to representative points,
#' and returns a tidy object used by \link{DisconnectivityGraph}.
#' If \code{bootstrap_ene} is provided, statistical significance level of 
#' energy gap (i.e., energy difference between stable state and tipping point) 
#' in each stable state is computed. 
#'
#' @param ela Tuple/list: (stablestates, ssenergies, tippingpoints, tippingenergies).
#' @param bootstrap_ene Optional matrix/data frame of bootstrap energies
#' @importFrom stats quantile
#' @return A data frame containing component sets, key, energy,
#'   x-position, and representative point,
#'   (plus \code{qt_min/qt_max} when \code{bootstrap_ene} is given).
#' @useDynLib ELAplus, .registration=TRUE
#' @noRd
GraphObj <- function(ela,bootstrap_ene=NULL){
  stablestates <- ela[[1]]
  ssenergies <- ela[[2]]
  tippingpoints <- ela[[3]]
  tippingenergies <- ela[[4]]
  energyall <- setdiff(union(ssenergies, unlist(tippingenergies)), Inf)
  rre_pre <- rbind(
    cbind(ssenergies, stablestates) %>% 'colnames<-'(c('energy', 'point')),
    cbind(
      tippingenergies %>% as.data.frame %>% 'rownames<-'(
        1: length(tippingenergies)) %>% 'colnames<-'(
        1: length(t(tippingenergies))) %>% pivot_longer(
          cols=1: length(t(tippingenergies))),
    tippingpoints %>% as.data.frame %>% 'rownames<-'(
      1: length(tippingpoints)) %>% 'colnames<-'(
        1: length(t(tippingpoints))) %>% pivot_longer(
          cols=1: length(t(tippingpoints)))) %>%
    'colnames<-'(c('V1.1', 'energy', 'V1.2', 'point')) %>%
    select(c('energy', 'point')))
  rre <- rre_pre[!duplicated(rre_pre),]
  eee <- cbind(energyall, lapply(energyall, function(ww){calc_ccc(
    ww, stablestates, ssenergies, tippingenergies)})) %>% 'colnames<-'(c('energy', 'ccc')) %>%
    as.data.frame
  eee$energy <- unlist(eee[,'energy'])
  eee <- eee[order(eee$energy, decreasing=TRUE),]
  ese <- eee[!duplicated(eee$ccc, fromLast=TRUE),]
  eps <- ese$ccc
  ori <- eps[1][[1]]$'1'
  hee <- apply(cbind(eps %>% map(function(x){t(t(x))}), eps %>% map(
    function(x){t(setdiff(ori, unlist(x)))})), 1, function(x){unlist(
      x, recursive=FALSE)})
  cee <- calc_cee(hee)
  rcee <- cee[length(cee): 1]
  xpositions <- do.call(rbind, calc_xpos(rcee))[, 2][length(rcee): 1]
  invalid_shape <- lapply(xpositions, function(x){dim(data.frame(x))[1]}) != 1
  xpositions[invalid_shape] <- lapply(
    xpositions[invalid_shape],
    function(x){lapply(unlist(x), function(y){y})})
  nodes2xposi_pre <- cbind(
    cee=unlist(cee, recursive=FALSE),
    xpositions=unlist(xpositions, recursive=FALSE)) %>% as.data.frame
  nodes2xposi_pre$xpositions <- nodes2xposi_pre$xpositions %>% as.double
  nodes2xposi <- nodes2xposi_pre[!duplicated(nodes2xposi_pre),]
  ff_pre <- do.call(
    rbind,
    ese %>% apply(1, function(x){cbind(ccc=x$ccc, energy=x$energy)})) %>%
    as.data.frame
  ff_pre2 <- cbind(
    ccc=ff_pre$ccc,
    ccc_str=ff_pre$ccc %>% as.character,
    energy=ff_pre$energy %>% as.double)
  ff <- ff_pre2 %>% as.data.frame %>% left_join(
    cbind(
      cee=nodes2xposi$cee,
      cee_str=nodes2xposi$cee %>% as.character,
      nodes2xposi=nodes2xposi$xpositions) %>% as.data.frame,
    by=c('ccc_str' = 'cee_str'))
  ff <- ff[, c('cee', 'ccc_str', 'energy', 'nodes2xposi')]
  ff$ccc_str <- ff$ccc_str %>% as.character
  ff[, c('energy', 'nodes2xposi')] <- ff[, c('energy', 'nodes2xposi')] %>%
    apply(1, as.double) %>% T
  graphinfo_pre <- ff[order(ff$ccc_str, ff$energy),]
  graphinfo_pre2 <- graphinfo_pre[!duplicated(graphinfo_pre$ccc_str),]
  graphinfo_pre2$len_cee <- lapply(
    graphinfo_pre2$cee, function(x){length(x)}) %>% as.integer
  graphinfo_pre2$sum_cee <- lapply(graphinfo_pre2$cee, sum) %>% as.integer
  graphinfo <- graphinfo_pre2[order(
    graphinfo_pre2$len_cee, graphinfo_pre2$sum_cee),]
  mm <- as.data.frame(as.character(graphinfo$energy))
  colnames(mm) <- "energy"
  mmm <- (mm  %>% left_join(as.data.frame(rre),by='energy'))
  point <- mmm$point
  mmm <- mmm[!duplicated(mmm),]
  rownames(mmm) <- mmm$point
  graphinfo <- cbind(graphinfo, point)
  if(!is.null(bootstrap_ene)){
    ssobj <- graphinfo[graphinfo$len_cee==1,]
    for(i in 1:dim(ssobj)[1]){
      ss_str <- ssobj[i,"ccc_str"]
      ssid <- ssobj[i,"point"]
      tipp_pre <-  graphinfo[sapply(graphinfo$ccc_str, function(x) {
        v <- eval(parse(text=x))
        ss_str %in% unlist(v)
      }), ]
      tipp_pre <- tipp_pre[tipp_pre$len_cee > 1,]
      tippid <- tipp_pre[order(tipp_pre$len_cee),][1,"point"]
      graphinfo[graphinfo$ccc_str==ss_str,"tippid"] <- tippid
      deltaE <- bootstrap_ene[tippid,] - bootstrap_ene[ssid,]
      graphinfo[graphinfo$ccc_str==ss_str,'pvalue'] <- length(deltaE[deltaE < 0])/length(deltaE)
    }
    graphinfo_final <- graphinfo[, c(
      'cee', 'ccc_str', 'energy', 'nodes2xposi', 'point', 'tippid', 'pvalue')]
  }else{
    graphinfo <- cbind(graphinfo, point)
    graphinfo_final <- graphinfo[, c(
      'cee', 'ccc_str', 'energy', 'nodes2xposi', 'point')]
  }
  return(graphinfo_final)}


#' DisconnectivityGraph: Generate a disconnectivity graph
#'
#' @description
#' Draws a disconnectivity graph from the data frame returned by \code{GraphObj}.
#' Nodes are plotted at their assigned x-positions and energy levels, and are
#' connected to their parent components by paths.
#'
#' If \code{show_pvalue = TRUE}, approximate p-values for the energy gaps of
#' stable states are displayed. The p-values are shown on a log10 scale using
#' the values specified in \code{scale_labels}.
#'
#' @param grobj A data frame returned by \code{GraphObj}.
#' @param s Integer; number of species/features (used for labeling by \link{id2bin}).
#' @param DG_sample Character; title of the plot.
#' @param annot_adj Numeric vector of length 2; horizontal and vertical adjustment for labels.
#' @param min_pval Numeric (must be > 0); minimum p-value shown in the visualization. This is
#' usually set to the minimum detection level in bootstrap sampling
#' (e.g., if the number of bootstrap iterations is 1000, the minimum level is
#' 1/1000 = 0.001).
#' @param scale_labels Numeric vector of values between 0 and 1; labels for the
#' p-value scale. These values are converted to log-scaled values for plotting.
#' @param show_pvalue Logical; if \code{TRUE}, plots approximate p-values based
#' on the energy gaps between each stable state and its corresponding tipping
#' point estimated by bootstrap resampling.
#' @param getGraphObj Logical; if \code{TRUE}, returns the data frame used for
#' plotting the disconnectivity graph, so that users can manually modify the
#' graphical settings.
#' @param plot Logical; if TRUE, show disconnectivity graph 
#' @return
#' If \code{getGraphObj = TRUE}, returns a data frame containing energy values,
#' x-positions, and p-values for the energy gaps of the identified stable states
#' and tipping points. Otherwise, returns \code{NULL}.
#'
#' @importFrom purrr map
#' @importFrom ggplot2 ggplot aes geom_rect geom_point geom_text geom_path xlab ggtitle theme coord_cartesian scale_size scale_color_viridis_c
#' @importFrom rlang .data
#' @useDynLib ELAplus, .registration=TRUE
#' @export
DisconnectivityGraph <- function(
    grobj, s, DG_sample='DG_sample', annot_adj=c(0.75, 2.00), min_pval = 5e-02, scale_labels = c(0.05,0.1,0.2,0.5),
    show_pvalue=FALSE, getGraphObj=FALSE, plot =TRUE){

  range_ <- c(min(grobj$energy), max(grobj$energy))
  grobj_pre <- grobj
  grobj_pre$len_cee <- lapply(
    grobj_pre$cee, function(x){length(x)}) %>% as.integer
  grobj_ <- cbind(
    grobj_pre,
    do.call(
      cbind,
      grobj_pre$point %>% map(function(x){id2bin(x, s)})) %>% t
      )
  jun_pre <- grobj_[grobj_$len_cee == 1,]
  jun <- jun_pre[order(jun_pre$nodes2xposi),]
  jen_pre <- grobj_[grobj_$len_cee > 1,]
  jen <- jen_pre[order(jen_pre$nodes2xposi),]

  # grobj to plot
  grobj_to_plot <- data.frame()
  if(nrow(grobj_) != 1){
    for(i in 1: (nrow(grobj_))){
      aa <- grobj_[i,]
      #aa$point_str <- paste0('C', aa$point)
      aa$point_str <- aa$point
      bb_pre <- grobj_[grobj_$cee %>% map(
        function(x){all(unlist(aa$cee) %in% unlist(x))}) %>% unlist(),]
      bb <- bb_pre[bb_pre$ccc_str != aa$ccc_str,][1,]
      bb$point_str <- rep("", length(bb$point))
      # bb$point_str <- paste0('C', bb$point)
      between_aa_bb <- aa
      between_aa_bb$energy <- bb$energy
      between_aa_bb$point_str <- ''
      line_break <- between_aa_bb
      line_break$nodes2xposi <- NA
      line_break$energy <- NA
      grobj_to_plot_ <- rbind(aa, between_aa_bb, bb, line_break)
      grobj_to_plot_$id <- i
      grobj_to_plot <- grobj_to_plot %>% rbind(grobj_to_plot_)
    }
  }

  # scatter with annot and line
  g <- ggplot(
    grobj_to_plot,
    aes(x=.data$nodes2xposi, y=.data$energy, label=.data$point_str))
  xrange <- range(grobj_to_plot$nodes2xposi, na.rm=TRUE)
  yrange <- range(grobj_to_plot$energy, na.rm=TRUE)
  xlim <- c(xrange[[1]]-0.05,xrange[[2]]+0.05)
  ylim <- c(yrange[[1]]-0.25,yrange[[2]]+0.25)
  
  g <- g + xlab("") + geom_text(
    hjust=annot_adj[1], vjust=annot_adj[2], aes(fontface=2)) + geom_path()
  
  grobj_to_plot <- grobj_to_plot[!is.na(grobj_to_plot$len_cee), ]
  grobj_to_plot[as.logical((grobj_to_plot$len_cee > 1)*(grobj_to_plot$point_str != "")),"type"] <- "tipp"
  grobj_to_plot[as.logical((grobj_to_plot$len_cee == 1)*(grobj_to_plot$point_str != "")),"type"] <- "ss"
  
  grobj_point <- grobj_to_plot[grobj_to_plot$point_str != "",]
  grobj_tipp <- grobj_point[grobj_point$len_cee > 1,]
  grobj_ss <- grobj_point[grobj_point$len_cee == 1,]
  
  if (show_pvalue){
    #log_scale_labels <- scale_labels
    log_scale_labels <- log10(scale_labels)
    nom_scale_labels <- (log_scale_labels - min(log_scale_labels))/(max(log_scale_labels) - min(log_scale_labels))*3 + 1
    
    min_pval <- max(min(scale_labels),min_pval)
    max_pval <- max(scale_labels)
    grobj_ss[grobj_ss$pvalue < min_pval,"pvalue"] <- min_pval
    grobj_ss[grobj_ss$pvalue > max_pval,"pvalue"] <- max_pval
    log_pvalue <- log10(grobj_ss$pvalue)
    #log_pvalue <- grobj_ss$pvalue
    nom_pvalue <- (log_pvalue - min(log_scale_labels))/(max(log_scale_labels) - min(log_scale_labels))*3 + 1
    print(nom_pvalue)
    
    log_scale_labels <- trunc(log_scale_labels * 100) / 100
    
    grobj_ss[,"pval_size"] <- nom_pvalue
    g <- g + 
      geom_point(data = grobj_tipp, color = "black",shape = 15,size=2)
    
    g <- g +
      geom_point(
        data = grobj_ss,
        mapping = aes(
          x = .data$nodes2xposi,
          y = .data$energy,
          color = .data$pval_size,
          size = .data$pval_size
        )
      ) +
      scale_size(range = c(3, 7),limits = c(1,4), name="",
                 breaks=nom_scale_labels,labels=as.character(log_scale_labels)) +
      scale_color_viridis_c(option = "viridis",alpha=0.8,
                            direction=-1,na.value = "grey",limits=c(1,4),
                            name="log10(Bootstrap-approx. p-val)",
                            breaks=nom_scale_labels,labels=as.character(log_scale_labels))
  
  }else{
    g <- g + 
      geom_point(data = grobj_tipp, color = "black",shape = 15, size=2.5)
    
    g <- g +
      geom_point(data = grobj_ss,color="black",shape = 19, size = 3)
  }
  # pie
  #g <- g + geom_scatterpie(
  #  aes(x=nodes2xposi, y=energy),
  #  data=grobj_, cols=d_list)

  # plot
  g <- g + ggtitle(DG_sample) + coord_cartesian(xlim=xlim, ylim=ylim)
  if (plot){
    plot(g)
  }
  if (getGraphObj){
    return(grobj_to_plot)}
  else{
    return(NULL)
  }
}

#' showDG: draw DisconnectivityGraph
#'
#' @description Computing graph objects from \code{ela} (and optionally
#' bootstrap energies) and draw the disconnectivity graph for a given community
#' composition matrix \code{oc}. More than two stable state are required.
#'
#' @param ela ELA output tuple/list used by \code{GraphObj}.
#' @param oc Community composition matrix/data frame (used to infer \code{s = ncol(oc)}).
#' @param label Character; plot title.
#' @param bootstrap_ene Optional bootstrap energy matrix/data frame for computing statistical significance levels of energy gaps.
#' @param min_pval Numeric; minimum p-value shown in the visualization. This is
#' usually set to the minimum detection level in bootstrap sampling
#' (e.g., if the number of bootstrap iterations is 1000, the minimum level is
#' 1/1000 = 0.001).
#' @param scale_labels Numeric vector of values between 0 and 1; labels for the
#' p-value scale. These values are converted to log-scaled values for plotting.
#' @param plot Logical; if TRUE, show disconnectivity graph 
#' @param getGraphObj Logical; if TRUE, returns data.frame object used for plotting DisconnectivityGraph
#' Users can manually change the graphical configurations.
#' @return A data frame containing energy values, xpositions, pvalues for energy gaps in identified stable states and tipping points.
#' @useDynLib ELAplus, .registration=TRUE
#' @export
showDG <- function(ela, oc, label="",bootstrap_ene=NULL,min_pval = 5e-02, plot=TRUE,
                   scale_labels = c(0.05,0.1,0.2,0.5), getGraphObj=FALSE){
  if(length(ela[[1]])>1){
  s <- ncol(oc)
  if(is.null(bootstrap_ene)){
    grobj <- GraphObj(ela)
    plotObj <- DisconnectivityGraph(grobj, s, label, annot_adj=c(0.75, 2.00), min_pval = min_pval,
                                    scale_labels = scale_labels, getGraphObj = TRUE, plot =plot)
  }else{
    grobj <- GraphObj(ela,bootstrap_ene=bootstrap_ene)
    plotObj <- DisconnectivityGraph(grobj, s, label, annot_adj=c(0.75, 2.00), min_pval = min_pval,
                                    scale_labels = scale_labels,show_pvalue = TRUE,
                                    getGraphObj = TRUE, plot =plot)
  }
  if (getGraphObj){
    return(plotObj)
  }else{
    return(NULL)
  }
  }else{
      return(cat("only one stable state found\n"))
  }}

#' Plot histograms of bootstrap deltaE values
#'
#' For each stable state, this function computes the difference between the energy values 
#' of the corresponding tipping point and the stable state across bootstrap samples,
#' then draws a histogram of the resulting deltaE values.
#'
#' @param bootstrap_ene A data frame containing energies in stable states and tipping points,
#' output of \link{bootstrap_ELA}.
#' @param grobj_to_plot A data frame containing stable state - tipping point relationship,
#' output of \link{showDG}
#' @param getData Logical. If `TRUE`, return the computed deltaE values as a list.
#'   If `FALSE`, return `NULL`.
#'
#' @return
#' If `getData = TRUE`, a list of numeric vectors, one per deltaE for each stable state.
#' If `getData = FALSE`, `NULL`.
#'
#' @export

deltaE_histogram <- function(bootstrap_ene, grobj_to_plot, getData = TRUE) {
  grobj_point <- grobj_to_plot[grobj_to_plot$point_str != "", ]
  grobj_tipp  <- grobj_point[grobj_point$len_cee > 1, ]
  ssobj       <- grobj_point[grobj_point$len_cee == 1, ]
  
  n <- nrow(ssobj)
  nr <- ceiling(sqrt(n))
  nc <- ceiling(n / nr)
  
  op <- graphics::par(mfrow = c(nr, nc), mar = c(4, 4, 2, 1))
  on.exit(graphics::par(op))
  
  deltaE_list <- vector("list", n)
  
  for (i in seq_len(n)) {
    ssid   <- ssobj[i, "point_str"]
    tippid <- ssobj[i, "tippid"]
    deltaE <- bootstrap_ene[tippid, ] - bootstrap_ene[ssid, ]
    
    graphics::hist(as.numeric(deltaE),
         main = paste0(ssid, "-", tippid),
         xlab = "deltaE")
    
    deltaE_list[[i]] <- deltaE
  }
  
  if (getData) {
    return(deltaE_list)
  } else {
    return(NULL)
  }
}


#' showIntrGraph: draw interaction matrix on a PCoA plot
#'
#' Computes a 2D PCoA embedding from the interaction matrix (\code{jest}) and overlays
#' interaction links above a threshold \code{th}. Positive and negative links are
#' drawn separately, and species names are annotated at their coordinates.
#'
#' @param ela \link{ELA} output tuple/list (used to obtain stable states etc.; currently
#'   mainly used for context).
#' @param sa Parameter set passed to \link{sa2params}.
#' @param th Numeric threshold for drawing links (absolute value).
#' @param annot_adj Numeric length-2; horizontal/vertical adjustment for labels.
#' @param env target environmental factor.
#' @return NULL.
#' @useDynLib ELAplus, .registration=TRUE
#' @importFrom ggplot2 scale_color_manual theme_classic
#' @export
showIntrGraph <- function(ela, sa, env=NULL, th=0., annot_adj=c(0.75, 2.00)){
  ######################
  sa2p <- sa2params(sa, env)
  hestp <- sa2p[[1]]
  jestp <- sa2p[[2]]
  gestp <- sa2p[[3]]
  hgestp <- sa2p[[4]] 
  ####################
  jest <- jestp
  s <- length(hestp)
  tss <- ela[[1]]
  tssen <- ela[[2]]
  tti <- ela[[3]]
  tien <- ela[[4]]
  #memberq <- do.call(cbind, tss %>% map(
  #  function(x){CIntegerDigits(ssid=x, n=s)})) %>%
  #  'colnames<-'(tss) %>%
  #  apply(2, function(x){names(x[x == 1]) %>% as.integer()}) %>%
  #  as.data.frame() %>% as.list() %>% 'names<-'(NULL)
  
  ############
  rn <- rownames(jest)
  colnames(jest) <- rn
  jj <- jest
  jest <- as.data.frame(jest)
  ############
  
  #fmq <- memberq %>% unlist() %>% unique() %>% sort()
  #tname <- colnames(jest)[fmq]
  tname <- colnames(jest)
  mat <- jest[tname, tname]
  mat$species1 <- rownames(mat)
  link_ <- mat %>% pivot_longer(cols=tname, names_to='species2') %>% as.data.frame()
  plink_pre <- link_[link_$value > th, c("species1", "species2"), drop = FALSE]
  qlink_pre <- link_[link_$value < -th, c("species1", "species2"), drop = FALSE]
  pcp <- as.data.frame(stats::cmdscale(jj, k=2)) %>%
    'colnames<-'(lapply(1:2, function(x){
      gsub(' ', '', paste('PCo', as.character(x)), fixed=TRUE)}))
  lmap <- cbind(
    link_[, c('species1', 'species2')],
    lmap=link_$value %>% map(function(x){min(max(-1, x), 1)}) %>% unlist())
  tmap <- cbind(
    link_[, c('species1', 'species2')],
    tmap=link_$value %>% map(function(x){min(max(-1, x), 1)}) %>% unlist())
  ltmap <- cbind(lmap, tmap=tmap$tmap)
  plink <- plink_pre %>% left_join(ltmap, by=c('species1', 'species2'))
  qlink <- qlink_pre %>% left_join(ltmap, by=c('species1', 'species2'))
  
  # grobj to plot
  grobj_to_plot <- data.frame()
  if(nrow(plink) != 1){
    for(i in 1: nrow(plink)){
      grobj_ <- pcp[unlist(plink[i, c('species1', 'species2')]),]
      grobj_$lmap <- plink[i, 'lmap']
      grobj_$tmap <- plink[i, 'tmap']
      grobj_$id <- 'positive'
      rownames(grobj_) <- NULL
      line_break <- grobj_[1,]
      line_break$PCo1 <- NA
      line_break$PCo2 <- NA
      line_break$tmap <- NA
      line_break$lmap <- NA
      grobj_to_plot <- rbind(grobj_to_plot, grobj_, line_break)
    }
  }
  if(nrow(qlink) != 1){
    for(i in 1: nrow(qlink)){
      grobj_ <- pcp[unlist(qlink[i, c('species1', 'species2')]),]
      grobj_$lmap <- qlink[i, 'lmap']
      grobj_$tmap <- qlink[i, 'tmap']
      grobj_$id <- 'negative'
      rownames(grobj_) <- NULL
      line_break <- grobj_[1,]
      line_break$PCo1 <- NA
      line_break$PCo2 <- NA
      line_break$tmap <- NA
      line_break$lmap <- NA
      grobj_to_plot <- rbind(grobj_to_plot, grobj_, line_break)
    }
  }
  grobj_to_plot$species <- NA
  for(i in 1: nrow(pcp)){
    pcp_ <- pcp[i,]
    grobj_ <- cbind(pcp_, lmap=Inf, tmap=1, id='species', species=rownames(pcp_))
    line_break <- as.list(grobj_ %>% sapply(function(x){NA}))
    line_break$id <- 'species'
    line_break$species <- rownames(pcp_)
    grobj_to_plot <- rbind(grobj_to_plot, grobj_, line_break)
  }
  grobj_to_plot$lmap2 <- grobj_to_plot$lmap %>% sapply(function(x){stats::pnorm(x)})
  grobj_to_plot$tmap2 <- grobj_to_plot$tmap %>% sapply(function(x){1 * abs(x)})
  
  g <- ggplot(
    grobj_to_plot,
    aes(x=.data$PCo1, y=.data$PCo2, label=.data$species, 
        color=.data$id, alpha=.data$lmap2))
  g <- g + geom_point() + geom_text(
    hjust=annot_adj[1], vjust=annot_adj[2], aes(fontface=2)) + geom_path() +
    scale_color_manual(values = c('red','blue',  'black'))
  g <- g + theme_classic()
  
  # plot
  plot(g)
  
  return(NULL)}


#' showSSD: draw stable state diagram
#'
#' @description Plots energy values of stable states across a factor (e.g., environmental or
#' control parameter) using the output from \code{gela}. Each unique stable state is
#' drawn as a separate line with points, and a legend is added.
#'
#' @param gela The output of \link{GradELA}. Object containing stable state lists, energies, the factor,
#'   and environment labels as used in the code.
#' @return NULL (produces a plot as side effect).
#' @importFrom RColorBrewer brewer.pal
#' @importFrom graphics lines
#' @useDynLib ELAplus, .registration=TRUE
#' @export
showSSD <- function(gela){
  mem <- foreach(i=gela[[1]]) %do% {
      i[[1]][[1]]
  }
  uqss <- unique(unlist(mem))

  memen <- foreach(i=gela[[1]]) %do% {
    i[[1]][[2]]
    }
  minen <- min(unlist(memen))
  maxen <- max(unlist(memen))

  cols <- rep(brewer.pal(8, "Set1"), times=100)

  yrange <- c(minen-abs(0.05*minen), maxen+abs(0.05*maxen))
  graphics::par(xpd=T)
  graphics::par(mar=c(6,6,3,12))

  cc <- 0
  factor <- gela[[2]]
  env <- gela[[3]]
  for(j in uqss){
    cc <- cc + 1
    
    po <- foreach(i = gela[[1]]) %do% {
      pp <- which(i[[1]][[1]] == j)
      if(identical(pp, integer(0))) NA else i[[1]][[2]][pp]
    }
    
    energy <- unlist(po)
    
    if(cc == 1){
      plot(factor, energy,
           type = "o",
           ylim = yrange,
           col = cols[cc],
           xlab = names(env[is.na(env)]),
           ylab = "Energy")
    } else {
      lines(factor, energy, type = "o", col = cols[cc])
    }
  }

  graphics::legend(graphics::par()$usr[2] + 0.2, graphics::par()$usr[4], legend = uqss, col = cols, pch=1, lty=1)}


#' showSSD_ggplot: draw stable state diagram using ggplot
#'
#' @description Plots energy values of stable states across a factor (e.g., environmental or
#' control parameter) using the output from \code{gela}. Each unique stable state is
#' drawn as a separate line with points, and a legend is added.
#'
#' @param gela The output of \link{GradELA}. Object containing stable state lists, energies, the factor,
#'   and environment labels as used in the code.
#' @param plot Logical; if TRUE, show SSD plot 
#' @param getSSDobj Logical; if TRUE, returns data.frame object used for plotting SSD diagram
#' Users can manually change the graphical configurations.
#' @return A data frame used for plotting SSD diagram.
#' @importFrom RColorBrewer brewer.pal
#' @useDynLib ELAplus, .registration=TRUE
#' @export
showSSD_ggplot <- function(gela, plot = TRUE, getSSDobj=TRUE) {
  egrad <- gela[[2]]
  ssd_df <- foreach(i = 1:length(egrad), .combine = "rbind") %do% {
    data.frame(
      env = egrad[[i]],
      SS_ID = gela[[1]][[i]][[1]][[1]],
      Energy = gela[[1]][[i]][[1]][[2]]
    )
  }
  
  if (plot) {
    g <- ggplot(ssd_df, aes(x = .data$env, y = .data$Energy, color = .data$SS_ID)) +
      geom_point(shape = 21) +
      geom_line(aes(group = .data$SS_ID))
    print(g)
  }
  if(getSSDobj){
    return(ssd_df)
  }else{return(NULL)}
}

#' GELA3D: draw 3D plot of energy landscape
#'
#' @description Draws a 3D surface of the energy landscape and overlays
#' line trajectories using \pkg{plot3D} function.
#'
#' @param gelsobj List with two elements: surface data frame and line data frame.
#' @param theta numeric value. plot3D option: Azimuth angle (in degrees) controlling rotation around the z-axis.
#' @param phi numeric value. plot3D option: Elevation angle (in degrees) controlling the vertical viewing position.
#' @importFrom stringr str_c
#' @return NULL (produces a 3D plot as side effect).
#' @useDynLib ELAplus, .registration=TRUE
#' @export
showGELA3D <- function(gelsobj,theta=25,phi=50) {
    # surface data
    surf3d <- gelsobj[[1]]
    # line data
    l3d <- gelsobj[[2]]
    GetMesh <- function (surf3d) {
        nx <- length(union(surf3d$x, c()))
        ny <- length(union(surf3d$y, c()))
        x <- matrix(0, nrow = nx, ncol = ny)
        y <- matrix(0, nrow = nx, ncol = ny)
        z <- matrix(0, nrow = nx, ncol = ny)
        for (i in 1:nx) {
            for (j in 1:ny) {
                x[i, j] <- surf3d[i + (j - 1) * nx, ]$x
                y[i, j] <- surf3d[i + (j - 1) * nx, ]$y
                z[i, j] <- surf3d[i + (j - 1) * nx, ]$z
            }
        }
        list(x, y, z)
    }
    x <- GetMesh(surf3d)[[1]]
    y <- GetMesh(surf3d)[[2]]
    z <- GetMesh(surf3d)[[3]]
    # surface plot
    plot3D::persp3D(
        x, y, z, 
        theta = theta, phi = phi, expand = 0.5, facets=NA, ticktype = "detailed",
        alpha = 0.1, 
        main = "GradELA with plot3D",
        xlab = "X", ylab = "Y", zlab = "Energy"
    )
    # lines plot
    for (a_l3d in split(l3d, l3d$Index)) {
        ifelse(length(a_l3d$x) == 1, type <- "p", type <- "l")
        x <- unlist(a_l3d$x)
        y <- unlist(a_l3d$y)
        z <- unlist(a_l3d$z)
        plot3D::scatter3D(
            x, y, z,
            type=type, add=TRUE, col='green', lwd=4
        )
        i <- a_l3d$Index[[1]]
        plot3D::text3D(
            mean(x), mean(y), mean(z)*0.95, 
            labels = str_c("C", i), add=TRUE, cex=1.2
        )
    }
}