#' @title Plotting clusters in a heatmap representation of the cell distances

#' @description  This functions plots a heatmap of the distance matrix grouped

#'   by clusters. Individual clusters are highlighted with rainbow colors along

#'   the x and y-axes.

#' @param object \code{DISCBIO} class object.

#' @param clustering_method either "k-means" or "model-based" ("k" and "mb" are

#' also accepted)

#' @param hmethod  Agglomeration method used for determining the cluster order

#'   from hierarchical clustering of the cluster medoids. This should be one of

#'   "ward.D", "ward.D2", "single", "complete", "average". Default is "single".

#' @param quiet if `TRUE`, intermediary output is suppressed

#' @param rseed Random integer to fix random results.

#' @param plot if `TRUE`, plots the heatmap; otherwise, just prints cclmo

#' @return Unless otherwise specified, a heatmap and a vector of the underlying

#'   cluster order.

#' @importFrom stats hclust as.dist cor

setGeneric(

  "clustheatmap",

  function(object,

           clustering_method = "k-means",

           hmethod = "single",

           rseed = NULL,

           quiet = FALSE,

           plot = TRUE) {

    standardGeneric("clustheatmap")

  }

)

#' @export

#' @rdname clustheatmap

setMethod(

  "clustheatmap",

  signature = "DISCBIO",

  definition = function(

    object, clustering_method, hmethod, rseed, quiet, plot

  ) {

    x <- object@fdata

    if (tolower(clustering_method) %in% c("k-means", "k")) {

      part <- object@kmeans$kpart

    } else if (tolower(clustering_method) %in% c("model-based", "mb")) {

      object@clusterpar$metric <- "pearson"

      y <- clustfun(

        object@fdata,

        clustnr = 20,

        bootnr = 50,

        metric = "pearson",

        do.gap = TRUE,

        SE.method = "Tibs2001SEmax",

        SE.factor = .25,

        B.gap = 50,

        cln = 0,

        rseed = rseed,

        quiet = quiet

      )

      object@distances <- as.matrix(y$di)

      part <- object@MBclusters$clusterid

    }

    na <- vector()

    j <- 0

    for (i in 1:max(part)) {

      if (sum(part == i) == 0) {

        next

      }

      j <- j + 1

      na <- append(na, i)

      d <- x[, part == i]

      if (sum(part == i) == 1) {

        cent <- d

      } else {

        cent <- apply(d, 1, mean)

      }

      if (j == 1) {

        tmp <- data.frame(cent)

      } else {

        tmp <- cbind(tmp, cent)

      }

    }

    names(tmp) <- paste("cl", na, sep = ".")

    if (max(part) > 1) {

      cclmo <-

        hclust(dist.gen(as.matrix(

          dist.gen(t(tmp), method = object@clusterpar$metric)

        )), method = hmethod)$order

    } else {

      cclmo <- 1

    }

    q <- part

    for (i in 1:max(part)) {

      q[part == na[cclmo[i]]] <- i

    }

    part <- q

    di <- as.data.frame(as.matrix(dist.gen(t(object@distances))))

    pto <- part[order(part, decreasing = FALSE)]

    ptn <- vector()

    for (i in 1:max(pto)) {

      pt <-

        names(pto)[pto == i]

      z <-

        if (length(pt) == 1) {

          pt

        } else {

          pt[hclust(as.dist(t(di[pt, pt])), method = hmethod)$order]

        }

      ptn <- append(ptn, z)

    }

    col <- c("black", "blue", "green", "red", "yellow", "gray")

    mi <- min(di, na.rm = TRUE)

    ma <- max(di, na.rm = TRUE)

    if (plot) {

      layout(

        matrix(

          data = c(1, 3, 2, 4),

          nrow = 2,

          ncol = 2

        ),

        widths = c(5, 1, 5, 1),

        heights = c(5, 1, 1, 1)

      )

      ColorRamp <- colorRampPalette(

        brewer.pal(n = 7, name = "RdYlBu")

      )(100)

      ColorLevels <- seq(mi, ma, length = length(ColorRamp))

      if (mi == ma) {

        ColorLevels <- seq(

          0.99 * mi, 1.01 * ma,

          length = length(ColorRamp)

        )

      }

      opar <- withr::local_par(mar = c(3, 5, 2.5, 2))

      on.exit(withr::local_par(opar))

      image(as.matrix(di[ptn, ptn]), col = ColorRamp, axes = FALSE)

      abline(0, 1)

      box()

      tmp <- vector()

      for (u in 1:max(part)) {

        ol <- (0:(length(part) - 1) /

          (length(part) - 1))[ptn %in% names(x)[part == u]]

        points(

          rep(0, length(ol)),

          ol,

          col = col[cclmo[u]],

          pch = 15,

          cex = .75

        )

        points(

          ol,

          rep(0, length(ol)),

          col = col[cclmo[u]],

          pch = 15,

          cex = .75

        )

        tmp <- append(tmp, mean(ol))

      }

      axis(1, at = tmp, labels = cclmo)

      axis(2, at = tmp, labels = cclmo)

      opar <- withr::local_par(mar = c(3, 2.5, 2.5, 2))

      on.exit(withr::local_par(opar))

      image(

        1,

        ColorLevels,

        matrix(

          data = ColorLevels,

          ncol = length(ColorLevels),

          nrow = 1

        ),

        col = ColorRamp,

        xlab = "",

        ylab = "",

        las = 2,

        xaxt = "n"

      )

      layout(1)

    }

    return(cclmo)

  }

)