.initializeCluster <- function(N,

     len,
     z = NULL,
     initial = NULL,

     fixed = NULL) {

     # If initial values are given, then they will not be randomly initialized
     if (!is.null(initial)) {

         initValues <- sort(unique(initial))

         if (length(unique(initial)) != N || length(initial) != len ||

             !all(initValues %in% seq(N))) {

                 stop("'initial' needs to be a vector of length 'len'",
                     " containing N unique values.")

         }

         z <- as.integer(as.factor(initial))

     } else {
         z <- rep(NA, len)
     }

     # Set any values that need to be fixed during sampling
     if (!is.null(fixed)) {

         fixedValues <- sort(unique(fixed))
         if (length(fixed) != len || !all(fixedValues %in% seq(N))) {

             stop("'fixed' to be a vector of length 'len' where each entry is",
                 " one of N unique values or NA.")

         }

         fixedIx <- !is.na(fixed)
         z[fixedIx] <- fixed[fixedIx]
         zNotUsed <- setdiff(seq(N), unique(fixed[fixedIx]))

     } else {

         zNotUsed <- seq(N)
         fixedIx <- rep(FALSE, len)

     }

 
     # Randomly sample remaining values

     zNa <- which(is.na(z))
     if (length(zNa) > 0) {
         z[zNa] <- sample(zNotUsed, length(zNa), replace = TRUE)

     }

 
     # Check to ensure each value is in the vector at least once

     missing <- setdiff(seq(N), z)

     for (i in missing) {

         ta <- sort(table(z[!fixedIx]), decreasing = TRUE)

         if (ta[1] == 1) {
             stop("'len' is not long enough to accomodate 'N' unique values")
         }

         ix <- which(z == as.integer(names(ta))[1] & !fixedIx)

         z[sample(ix, 1)] <- i

     }

 
     return(z)

 }
 
 

 .initializeSplitZ <- function(counts,

     K,

     KSubcluster = NULL,

     alpha = 1,
     beta = 1,

     minCell = 3) {

 
     s <- rep(1, ncol(counts))

     if (is.null(KSubcluster))
         KSubcluster <- ceiling(sqrt(K))

     # Initialize the model with KSubcluster clusters

     res <- .celda_C(
         counts,

         K = KSubcluster,
         maxIter = 20,
         zInitialize = "random",

         alpha = alpha,
         beta = beta,

         splitOnIter = -1,
         splitOnLast = FALSE,

         verbose = FALSE,
         reorder = FALSE
     )

     overallZ <- as.integer(as.factor(clusters(res)$z))

     currentK <- max(overallZ)

     while (currentK < K) {

         # Determine which clusters are split-able

         # KRemaining <- K - currentK

         KPerCluster <- min(ceiling(K / currentK), KSubcluster)
         KToUse <- ifelse(KPerCluster < 2, 2, KPerCluster)

         zTa <- tabulate(overallZ, max(overallZ))

         zToSplit <- which(zTa > minCell & zTa > KToUse)
         if (length(zToSplit) > 1) {
             zToSplit <- sample(zToSplit)
         } else if (length(zToSplit) == 0) {

             break

         }
 

         # Cycle through each splitable cluster and split it up into
         # K.sublcusters
         for (i in zToSplit) {
             clustLabel <- .celda_C(counts[, overallZ == i, drop = FALSE],
                 K = KToUse,
                 zInitialize = "random",

                 alpha = alpha,
                 beta = beta,

                 maxIter = 20,
                 splitOnIter = -1,
                 splitOnLast = FALSE,

                 verbose = FALSE)

             tempZ <- as.integer(as.factor(clusters(clustLabel)$z))

             # Reassign clusters with label > 1

             splitIx <- tempZ > 1
             ix <- overallZ == i
             newZ <- overallZ[ix]
             newZ[splitIx] <- currentK + tempZ[splitIx] - 1

             overallZ[ix] <- newZ
             currentK <- max(overallZ)

 
             # Ensure that the maximum number of clusters does not get too large'

             if (currentK > K + 10) {

                 break

             }
         }

     }
 

     # Decompose counts for likelihood calculation

     p <- .cCDecomposeCounts(counts, s, overallZ, currentK)

     nS <- p$nS
     nG <- p$nG
     nM <- p$nM

     mCPByS <- p$mCPByS
     nGByCP <- p$nGByCP
     nCP <- p$nCP
     nByC <- p$nByC

 
     # Remove clusters 1-by-1 until K is reached

     while (currentK > K) {

         # Find second best assignment give current assignments for each cell

         probs <- .cCCalcEMProbZ(counts,

                 s = s,

                 z = overallZ,
                 K = currentK,
                 mCPByS = mCPByS,
                 nGByCP = nGByCP,
                 nByC = nByC,
                 nCP = nCP,

                 nG = nG,
                 nM = nM,
                 alpha = alpha,
                 beta = beta,

                 doSample = FALSE)
         zProb <- t(probs$probs)
         zProb[cbind(seq(nrow(zProb)), overallZ)] <- NA
         zSecond <- apply(zProb, 1, which.max)

         zTa <- tabulate(overallZ, currentK)
         zNonEmpty <- which(zTa > 0)

 
         # Find worst cluster by logLik to remove

         previousZ <- overallZ
         llShuffle <- rep(NA, currentK)
         for (i in zNonEmpty) {
             ix <- overallZ == i
             newZ <- overallZ
             newZ[ix] <- zSecond[ix]
 
             p <- .cCReDecomposeCounts(counts,

                 s,

                 newZ,
                 previousZ,
                 nGByCP,
                 currentK)
             nGByCP <- p$nGByCP
             mCPByS <- p$mCPByS
             llShuffle[i] <- .cCCalcLL(mCPByS,
                     nGByCP,

                     s,

                     newZ,
                     currentK,

                     nS,
                     nG,
                     alpha,
                     beta)

             previousZ <- newZ

         }

         # Remove the cluster which had the the largest likelihood after removal

         zToRemove <- which.max(llShuffle)

         ix <- overallZ == zToRemove
         overallZ[ix] <- zSecond[ix]

         p <- .cCReDecomposeCounts(counts,

             s,

             overallZ,
             previousZ,
             nGByCP,
             currentK)
         nGByCP <- p$nGByCP[, -zToRemove, drop = FALSE]
         mCPByS <- p$mCPByS[-zToRemove, , drop = FALSE]
         overallZ <- as.integer(as.factor(overallZ))
         currentK <- currentK - 1

     }

     return(overallZ)

 }
 

 .initializeSplitY <- function(counts,

     L,

     LSubcluster = NULL,
     tempK = 100,

     beta = 1,
     delta = 1,
     gamma = 1,

     minFeature = 3) {

 
     if (is.null(LSubcluster))
         LSubcluster <- ceiling(sqrt(L))

 
     # Collapse cells to managable number of clusters

     if (!is.null(tempK) && ncol(counts) > tempK) {

         z <- .initializeSplitZ(counts, K = tempK)

         counts <- .colSumByGroup(counts, z, length(unique(z)))

     }
 

     # Initialize the model with KSubcluster clusters

     res <- .celda_G(counts,

         L = LSubcluster,
         maxIter = 10,
         yInitialize = "random",

         beta = beta,
         delta = delta,
         gamma = gamma,

         splitOnIter = -1,
         splitOnLast = FALSE,

         verbose = FALSE,
         reorder = FALSE
     )

     overallY <- as.integer(as.factor(clusters(res)$y))

     currentL <- max(overallY)

     while (currentL < L) {

         # Determine which clusters are split-able

         yTa <- tabulate(overallY, max(overallY))
         yToSplit <- sample(which(yTa > minFeature & yTa > LSubcluster))

         if (length(yToSplit) == 0) {

             break

         }

         # Cycle through each splitable cluster and split it up into
         # LSublcusters

         for (i in yToSplit) {

             # make sure the colSums of subset counts is not 0
             countsY <- counts[overallY == i, , drop = FALSE]
             countsY <- countsY[, !(colSums(countsY) == 0)]
 
             if (ncol(countsY) == 0) {
                 next
             }
 

             clustLabel <- .celda_G(

                 countsY,

                 L = LSubcluster,
                 yInitialize = "random",

                 beta = beta,
                 delta = delta,
                 gamma = gamma,

                 maxIter = 20,
                 splitOnIter = -1,
                 splitOnLast = FALSE,

                 verbose = FALSE)

             tempY <- as.integer(as.factor(clusters(clustLabel)$y))

 
             # Reassign clusters with label > 1

             splitIx <- tempY > 1
             ix <- overallY == i
             newY <- overallY[ix]
             newY[splitIx] <- currentL + tempY[splitIx] - 1

             overallY[ix] <- newY
             currentL <- max(overallY)

 
             # Ensure that the maximum number of clusters does not get too large

             if (currentL > L + 10) {

                 break

             }
         }
     }
 
     ## Decompose counts for likelihood calculation

     p <- .cGDecomposeCounts(counts = counts, y = overallY, L = currentL)
     nTSByC <- p$nTSByC
     nByG <- p$nByG
     nByTS <- p$nByTS
     nGByTS <- p$nGByTS

     nM <- p$nM
     nG <- p$nG
     rm(p)
 
     # Pre-compute lgamma values

     lgbeta <- lgamma((seq(0, max(.colSums(counts, nrow(counts),
         ncol(counts))))) + beta)
     lggamma <- lgamma(seq(0, nrow(counts) + L) + gamma)
     lgdelta <- c(NA, lgamma(seq(nrow(counts) + L) * delta))

 
     # Remove clusters 1-by-1 until L is reached

     while (currentL > L) {

         # Find second best assignment give current assignments for each cell

         probs <- .cGCalcGibbsProbY(

             counts = counts,

             y = overallY,
             L = currentL,
             nTSByC = nTSByC,
             nByTS = nByTS,
             nGByTS = nGByTS,
             nByG = nByG,

             nG = nG,
             beta = beta,
             delta = delta,
             gamma = gamma,
             lgbeta = lgbeta,
             lggamma = lggamma,
             lgdelta = lgdelta,

             doSample = FALSE)
         yProb <- t(probs$probs)
         yProb[cbind(seq(nrow(yProb)), overallY)] <- NA
         ySecond <- apply(yProb, 1, which.max)

         yTa <- tabulate(overallY, currentL)
         yNonEmpty <- which(yTa > 0)

 
         # Find worst cluster by logLik to remove

         previousY <- overallY
         llShuffle <- rep(NA, currentL)
         for (i in yNonEmpty) {
             ix <- overallY == i
             newY <- overallY
             newY[ix] <- ySecond[ix]

 
             # Move arounds counts for likelihood calculation

             p <- .cGReDecomposeCounts(counts,
                     newY,
                     previousY,
                     nTSByC,
                     nByG,
                     currentL)
             nTSByC <- p$nTSByC
             nGByTS <- p$nGByTS
             nByTS <- p$nByTS
             llShuffle[i] <- .cGCalcLL(nTSByC,
                     nByTS,
                     nByG,
                     nGByTS,

                     nM,
                     nG,

                     currentL,

                     beta,
                     delta,
                     gamma)

             previousY <- newY

         }
 
         # Remove the cluster which had the the largest likelihood after removal

         yToRemove <- which.max(llShuffle)

         ix <- overallY == yToRemove
         overallY[ix] <- ySecond[ix]

 
         # Move around counts and remove module

         p <- .cGReDecomposeCounts(counts,
             overallY,
             previousY,
             nTSByC,
             nByG,
             currentL)
         nTSByC <- p$nTSByC[-yToRemove, , drop = FALSE]
         nGByTS <- p$nGByTS[-yToRemove]
         nByTS <- p$nByTS[-yToRemove]
         overallY <- as.integer(as.factor(overallY))
         currentL <- currentL - 1

     }

     return(overallY)

 }