#' @name getClustNum

#' @title Get estimation of optimal clustering number

#' @description This function provides two measurements (i.e., clustering prediction index [CPI] and Gap-statistics) and aims to search the optimal number for multi-omics integrative clustering. In short, the peaks reach by the red (CPI) and blue (Gap-statistics) lines should be referred to determine `N.clust`.

#' @param data List of matrices.

#' @param is.binary A logicial vector to indicate if the subdata is binary matrix of 0 and 1 such as mutation.

#' @param try.N.clust A integer vector to indicate possible choices of number of clusters.

#' @param center A logical value to indicate if the variables should be centered. TRUE by default.

#' @param scale A logical value to indicate if the variables should be scaled. FALSE by default.

#' @param fig.path A string value to indicate the output figure path.

#' @param fig.name A string value to indicate the name of the figure.

#' @export

#' @return A figure that helps to choose the optimal clustering number (argument of `N.clust`) for `get%algorithm_name%()` or `getMOIC()`, and a list contains the following components:

#'

#'         \code{CPI}   possible cluster number identified by clustering prediction index

#'

#'         \code{Gapk}  possible cluster number identified by Gap-statistics

#' @import IntNMF

#' @import mogsa

#' @import SNFtool

#' @importFrom ggplot2 alpha

#' @importFrom dplyr %>%

#' @importFrom grDevices dev.copy2pdf

#' @examples # There is no example and please refer to vignette.

#' @references Chalise P, Fridley BL (2017). Integrative clustering of multi-level omic data based on non-negative matrix factorization algorithm. PLoS One, 12(5):e0176278.

#'

#'Tibshirani, R., Walther, G., Hastie, T. (2001). Estimating the number of data clusters via the Gap statistic. J R Stat Soc Series B Stat Methodol, 63(2):411-423.

getClustNum <- function(data        = NULL,

                        is.binary   = rep(FALSE, length(data)),

                        try.N.clust = 2:8,

                        center      = TRUE,

                        scale       = TRUE,

                        fig.path    = getwd(),

                        fig.name    = "optimal_number_cluster") {

  # check data

  n_dat <- length(data)

  if(n_dat > 6){

    stop('current verision of MOVICS can support up to 6 datasets.')

  }

  if(n_dat < 2){

    stop('current verision of MOVICS needs at least 2 omics data.')

  }

  data.backup <- data # save a backup

  #--------------------------------------------#

  # Cluster Prediction Index (CPI) from IntNMF #

  # remove features that made of categories not equal to 2 otherwise Error in svd(X) : a dimension is zero

  if(!all(!is.binary)) {

    bindex <- which(is.binary == TRUE)

    for (i in bindex) {

      a <- which(rowSums(data[[i]]) == 0)

      b <- which(rowSums(data[[i]]) == ncol(data[[i]]))

      if(length(a) > 0) {

        data[[i]] <- data[[i]][which(rowSums(data[[i]]) != 0),] # remove all zero

      }

      if(length(b) > 0) {

        data[[i]] <- data[[i]][which(rowSums(data[[i]]) != ncol(data[[i]])),] # remove all one

      }

      if(length(a) + length(b) > 0) {

        message(paste0("--", names(data)[i],": a total of ",length(a) + length(b), " features were removed due to the categories were not equal to 2!"))

      }

    }

  }

  # In order to make the input data fit non-negativity constraint of intNMF,

  # the values of the data were shifted to positive direction by adding absolute value of the smallest negative number.

  # Further, each data was rescaled by dividing by maximum value of the data to make the magnitudes comparable (between 0 and 1) across the several datasets.

  dat <- lapply(data, function (dd){

    if (!all(dd >= 0)) dd <- pmax(dd + abs(min(dd)), 0) + .Machine$double.eps # .Machine$double.eps as The smallest positive floating-point number x

    dd <- dd/max(dd)

    return(dd %>% as.matrix)

  })

  #dat <- lapply(dat, t)

  dat <- lapply(dat, function(x) t(x) + .Machine$double.eps)

  message("calculating Cluster Prediction Index...")

  optk1 <- IntNMF::nmf.opt.k(dat      = dat,

                            n.runs    = 5,

                            n.fold    = 5,

                            k.range   = try.N.clust,

                            st.count  = 10,

                            maxiter   = 100,

                            make.plot = FALSE)

  optk1 <- as.data.frame(optk1)

  #-------------------------------#

  # Gap-statistics from MoCluster #

  message("calculating Gap-statistics...")

  moas <- data.backup %>% mogsa::mbpca(ncomp      = 2,

                                       k          = "all",

                                       method     = "globalScore",

                                       center     = center,

                                       scale      = scale,

                                       moa        = TRUE,

                                       svd.solver = "fast",

                                       maxiter    = 1000,

                                       verbose    = FALSE)

  gap <- mogsa::moGap(moas, K.max = max(try.N.clust), cluster = "hclust", plot = FALSE)

  optk2 <- as.data.frame(gap$Tab)[-1,] # remove k=1

  #---------------------#

  # Eigen-gaps from SNF #

  # message("Calculating Eigen-gaps...")

  # data <- lapply(data.backup, t)

  # dat <- lapply(data, function (dd){

  #   dd <- dd %>% as.matrix

  #   W <- dd %>% SNFtool::dist2(dd) %>% SNFtool::affinityMatrix(K = 30, sigma = 0.5)

  # })

  # W <-  SNFtool::SNF(Wall = dat,

  #                    K = 30,

  #                    t = 20)

  # optk3 <- SNFtool::estimateNumberOfClustersGivenGraph(W, NUMC = try.N.clust)

  # nemo.ag <- nemo.affinity.graph(data.backup, k = 20)

  # optk4 <- nemo.num.clusters(nemo.ag, NUMC = try.N.clust)

  #---------------------------#

  # calculate optimal N.clust #

  # N.clust <- as.numeric(names(which(table(c(as.numeric((which.max(apply(optk1, 1, mean)) + 1)),

  #                                           as.numeric((which.max(optk2$gap) + 1)),

  #                                           optk3$`Eigen-gap best`)) >= 2)))

  N.clust <- as.numeric(which.max(apply(optk1, 1, mean) + optk2$gap)) + 1

  if(length(N.clust) == 0) {

    message("--fail to define the optimal cluster number!")

    N.clust <- "null"

  }

  #---------------#

  # Visualization #

  outFig <- paste0(fig.name,".pdf")

  par(bty="o", mgp = c(1.9,.33,0), mar=c(3.1,3.1,2.1,3.1)+.1, las=1, tcl=-.25)

  plot(NULL, NULL,

       xlim = c(min(try.N.clust),max(try.N.clust)),

       #ylim = c(min(optk1), max(optk1)),

       ylim = c(0,1),

       xlab = "Number of Multi-Omics Clusters",ylab = "")

  rect(par("usr")[1],par("usr")[3],par("usr")[2],par("usr")[4],col = "#EAE9E9",border = FALSE)

  grid(col = "white", lty = 1, lwd = 1.5)

  # for (m in 1:n.runs) points(try.N.clust, optk1[, m], pch = 20, cex = 1.5, col = "#224A8D")

  points(try.N.clust, apply(optk1, 1, mean), pch = 19, col = ggplot2::alpha("#224A8D"), cex = 1.5)

  lines(try.N.clust, apply(optk1, 1, mean), col = "#224A8D", lwd = 2, lty = 4)

  mtext("Cluster Prediction Index", side = 2, line = 2, cex = 1.5, col = "#224A8D", las = 3)

  par(new = TRUE, xpd = FALSE)

  plot(NULL,NULL,

       xlim = c(min(try.N.clust),max(try.N.clust)),

       #ylim = c(min(optk2$gap), max(optk2$gap)),

       ylim = c(0,1),

       xlab = "",ylab = "",xaxt = "n",yaxt = "n")

  points(try.N.clust, optk2$gap, pch = 19, col = ggplot2::alpha("#E51718",0.8), cex = 1.5)

  lines(try.N.clust, optk2$gap, col = "#E51718", lwd = 2, lty = 4)

  #axis(side = 4, at = pretty(range(optk2$gap), n = 6))

  axis(side = 4, at = seq(0,1,0.2), labels = c("0.0","0.2","0.4","0.6","0.8","1.0"))

  mtext("Gap-statistics", side = 4, line = 2,las = 3, cex = 1.5, col = "#E51718")

  # abline(v = optk3$`Eigen-gap best`, col = "#008B8A", lwd = 2, lty = 4)

  # text(optk3$`Eigen-gap best`, par("usr")[3] + 0.1,"Eigen-gaps", cex = 1.5, col = "#008B8A", adj = -0.05)

  invisible(dev.copy2pdf(file = file.path(fig.path, outFig), width = 5.5, height = 5))

  message("visualization done...")

  if(N.clust > 1) {

    message(paste0("--the imputed optimal cluster number is ", N.clust, " arbitrarily, but it would be better referring to other priori knowledge."))

  }

  #return(list(N.clust = N.clust, CPI = optk1, Gapk = optk2, Eigen = optk3))

  return(list(N.clust = N.clust, CPI = optk1, Gapk = optk2))

}