#' @title Inference of outlier cells

#' @description This functions performs the outlier identification for k-means

#' and model-based clustering

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

#' @param outminc minimal transcript count of a gene in a clusters to be tested

#'   for being an outlier gene. Default is 5.

#' @param outlg Minimum number of outlier genes required for being an outlier

#'   cell. Default is 2.

#' @param probthr outlier probability threshold for a minimum of \code{outlg}

#'   genes to be an outlier cell. This probability is computed from a negative

#'   binomial background model of expression in a cluster. Default is 0.001.

#' @param thr probability values for which the number of outliers is computed in

#'   order to plot the dependence of the number of outliers on the probability

#'   threshold. Default is 2**-(1:40).set

#' @param outdistquant Real number between zero and one. Outlier cells are

#'   merged to outlier clusters if their distance smaller than the

#'   outdistquant-quantile of the distance distribution of  pairs of cells in

#'   the orginal clusters after outlier removal. Default is 0.75.

#' @param K Number of clusters to be used.

#' @param plot if `TRUE`, produces a plot of -log10prob per K

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

#' @importFrom stats coef pnbinom

#' @return A named vector of the genes containing outlying cells and the number

#'   of cells on each.

#' @examples

#' sc <- DISCBIO(valuesG1msTest)

#' sc <- Clustexp(sc, cln = 2) # K-means clustering

#' FindOutliers(sc, K = 2)

#'

setGeneric(

  "FindOutliers",

  function(object, K, outminc = 5, outlg = 2, probthr = 1e-3, thr = 2**-(1:40),

           outdistquant = .75, plot = TRUE, quiet = FALSE) {

    standardGeneric("FindOutliers")

  }

)

#' @rdname FindOutliers

#' @export

setMethod(

  "FindOutliers",

  signature = "DISCBIO",

  definition = function(

    object, K, outminc, outlg, probthr, thr, outdistquant, plot, quiet

  ) {

    # ======================================================================

    # Validating

    # ======================================================================

    ran_k <- length(object@kmeans) > 0

    ran_m <- length(object@MBclusters) > 0

    if (ran_k) {

      clusters <- object@kmeans$kpart

    } else if (ran_m) {

      object <- Clustexp(

        object,

        clustnr = 20,

        bootnr = 50,

        metric = "pearson",

        do.gap = TRUE,

        SE.method = "Tibs2001SEmax",

        SE.factor = .25,

        B.gap = 50,

        cln = K,

        rseed = 17000,

        quiet = quiet

      )

      clusters <- object@MBclusters$clusterid

    } else {

      stop("run Clustexp before FindOutliers")

    }

    if (!is.numeric(outminc)) {

      stop("outminc has to be a non-negative integer")

    } else if (round(outminc) != outminc | outminc < 0) {

      stop("outminc has to be a non-negative integer")

    }

    if (!is.numeric(outlg)) {

      stop("outlg has to be a non-negative integer")

    } else if (round(outlg) != outlg | outlg < 0) {

      stop("outlg has to be a non-negative integer")

    }

    if (!is.numeric(probthr)) {

      stop("probthr has to be a number between 0 and 1")

    } else if (probthr < 0 | probthr > 1) {

      stop("probthr has to be a number between 0 and 1")

    }

    if (!is.numeric(thr)) {

      stop("thr hast to be a vector of numbers between 0 and 1")

    } else if (min(thr) < 0 | max(thr) > 1) {

      stop("thr hast to be a vector of numbers between 0 and 1")

    }

    if (!is.numeric(outdistquant)) {

      stop("outdistquant has to be a number between 0 and 1")

    } else if (outdistquant < 0 | outdistquant > 1) {

      stop("outdistquant has to be a number between 0 and 1")

    }

    object@outlierpar <- list(

      outminc = outminc,

      outlg = outlg,

      probthr = probthr,

      thr = thr,

      outdistquant = outdistquant

    )

    ### calibrate background model

    EXP <- object@expdata + 0.1

    m <- log2(apply(EXP, 1, mean))

    v <- log2(apply(EXP, 1, var))

    f <- m > -Inf & v > -Inf

    m <- m[f]

    v <- v[f]

    mm <- -8

    repeat {

      fit <- lm(v ~ m + I(m^2))

      if (coef(fit)[3] >= 0 | mm >= 3) {

        break

      }

      mm <- mm + .5

      f <- m > mm

      m <- m[f]

      v <- v[f]

    }

    object@background <- list()

    object@background$vfit <- fit

    object@background$lvar <-

      function(x, object) {

        2**(

          coef(object@background$vfit)[1] +

            log2(x) * coef(object@background$vfit)[2] +

            coef(object@background$vfit)[3] * log2(x)**2

        )

      }

    object@background$lsize <-

      function(x, object) {

        x**2 / (max(x + 1e-6, object@background$lvar(x, object)) - x)

      }

    ### identify outliers

    out <- vector()

    stest <- rep(0, length(thr))

    cprobs <- vector()

    for (n in 1:max(clusters)) {

      if (sum(clusters == n) == 1) {

        cprobs <-

          append(cprobs, .5)

        names(cprobs)[length(cprobs)] <-

          names(clusters)[clusters == n]

        next

      }

      x <- object@fdata[, clusters == n]

      x <- x[apply(x, 1, max) > outminc, ]

      z <-

        t(apply(x, 1, function(x) {

          apply(cbind(

            pnbinom(

              round(x, 0),

              mu = mean(x),

              size = object@background$lsize(mean(x), object)

            ),

            1 - pnbinom(

              round(x, 0),

              mu = mean(x),

              size = object@background$lsize(mean(x), object)

            )

          ), 1, min)

        }))

      cp <-

        apply(z, 2, function(x) {

          y <-

            p.adjust(x, method = "BH")

          y <- y[order(y, decreasing = FALSE)]

          return(y[outlg])

        })

      f <- cp < probthr

      cprobs <- append(cprobs, cp)

      if (sum(f) > 0) {

        out <- append(out, names(x)[f])

      }

      for (j in seq_len(length(thr))) {

        stest[j] <- stest[j] + sum(cp < thr[j])

      }

    }

    object@out <-

      list(

        out = out,

        stest = stest,

        thr = thr,

        cprobs = cprobs

      )

    ### cluster outliers

    clp2p.cl <- vector()

    cols <- names(object@fdata)

    di <- as.data.frame(object@distances)

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

      tcol <- cols[clusters == i]

      if (sum(!(tcol %in% out)) > 1) {

        clp2p.cl <- append(

          clp2p.cl,

          as.vector(

            t(di[tcol[!(tcol %in% out)], tcol[!(tcol %in% out)]])

          )

        )

      }

    }

    clp2p.cl <- clp2p.cl[clp2p.cl > 0]

    cpart <- clusters

    cadd <- list()

    if (length(out) > 0) {

      if (length(out) == 1) {

        cadd <- list(out)

      } else {

        n <- out

        m <- as.data.frame(di[out, out])

        for (i in seq_len(length(out))) {

          if (length(n) > 1) {

            o <-

              order(

                apply(

                  cbind(m, seq_len(dim(m)[1])),

                  1,

                  function(x) {

                    min(x[1:(length(x) - 1)][-x[length(x)]])

                  }

                ),

                decreasing = FALSE

              )

            m <- m[o, o]

            n <- n[o]

            f <- m[, 1] < quantile(clp2p.cl, outdistquant) |

              m[, 1] == min(clp2p.cl)

            ind <- 1

            if (sum(f) > 1) {

              for (j in 2:sum(f)) {

                comp1 <- m[f, f][j, c(ind, j)]

              }

            }

            comp2 <- quantile(clp2p.cl, outdistquant)

            if (apply(comp1 > comp2, 1, sum) == 0) {

              ind <- append(ind, j)

            }

            cadd[[i]] <- n[f][ind]

            g <- !n %in% n[f][ind]

            n <- n[g]

            m <- m[g, g]

            if (sum(g) == 0) {

              break

            }

          } else if (length(n) == 1) {

            cadd[[i]] <- n

            break

          }

        }

      }

      for (i in seq_len(length(cadd))) {

        cpart[cols %in% cadd[[i]]] <- max(cpart) + 1

      }

    }

    ### determine final clusters

    object@cpart <- cpart

    object@fcol <- rainbow(max(cpart))

    p <- clusters[order(clusters, decreasing = FALSE)]

    x <- object@out$cprobs[names(p)]

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

    if (plot) {

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

        y <- -log10(x + 2.2e-16)

        y[p != i] <- 0

        if (i == 1) {

          b <-

            barplot(

              y,

              ylim = c(0, max(-log10(

                x + 2.2e-16

              )) * 2.1),

              col = fcol[i],

              border = fcol[i],

              names.arg = FALSE,

              ylab = "-log10prob"

            )

        } else {

          barplot(

            y,

            add = TRUE,

            col = fcol[i],

            border = fcol[i],

            names.arg = FALSE,

            axes = FALSE

          )

        }

      }

      abline(

        -log10(object@outlierpar$probthr),

        0,

        col = "black",

        lty = 2

      )

      d <- b[2, 1] - b[1, 1]

      y <- 0

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

        y <- append(y, b[sum(p <= i), 1] + d / 2)

      }

      axis(1, at = (y[1:(length(y) - 1)] + y[-1]) / 2, labels = 1:max(p))

      box()

    }

    if (!quiet) {

      message(

        "The following cells are considered outliers: ",

        which(object@cpart > K),

        "\n"

      )

      print(which(object@cpart > K))

    }

    LL <- which(object@cpart > K)

    return(LL)

  }

)