@@ -28,3 +28,4 @@ celda.Rproj

 .DS_Store

 src/*.o

 src/*.dll

+src/*.so

@@ -21,7 +21,7 @@ Imports:

     gtable,

     grDevices,

     graphics,

-    Rmpfr,

+    matrixStats,

     gplots,

     doParallel,

     cluster,

@@ -79,4 +79,3 @@ import(grid)

 import(gtable)

 import(scales)

 useDynLib(celda)

-importFrom(Rcpp, sourceCpp)

@@ -91,21 +91,21 @@ celda_C = function(counts, sample.label=NULL, K, alpha=1, beta=1,

   num.iter.without.improvement = 0L

   do.cell.split = TRUE

   while(iter <= max.iter & num.iter.without.improvement <= stop.iter) {

-    message(paste("Iteration: ", iter))

 #   next.z = cC.calcGibbsProbZ(counts=counts, m.CP.by.S=m.CP.by.S, n.G.by.CP=n.G.by.CP, n.by.C=n.by.C, n.CP=n.CP, z=z, s=s, K=K, nG=nG, nM=nM, alpha=alpha, beta=beta)

 #	next.z = cC.calcEMProbZ(counts=counts, m.CP.by.S=m.CP.by.S, n.G.by.CP=n.G.by.CP, n.by.C=n.by.C, n.CP=n.CP, z=z, s=s, K=K, nG=nG, nM=nM, alpha=alpha, beta=beta)

     next.z = do.call(algorithm.fun, list(counts=counts, m.CP.by.S=m.CP.by.S, n.G.by.CP=n.G.by.CP, n.by.C=n.by.C, n.CP=n.CP, z=z, s=s, K=K, nG=nG, nM=nM, alpha=alpha, beta=beta))

+

     m.CP.by.S = next.z$m.CP.by.S

     n.G.by.CP = next.z$n.G.by.CP

     n.CP = next.z$n.CP

     z = next.z$z

-    

+

     ## Perform split on i-th iteration of no improvement in log likelihood

     if(K > 2 & (((iter == max.iter | num.iter.without.improvement == stop.iter) & isTRUE(split.on.last)) | (split.on.iter > 0 & iter %% split.on.iter == 0 & isTRUE(do.cell.split)))) {

       logMessages(date(), " ... Determining if any cell clusters should be split.", logfile=logfile, append=TRUE, sep="")

-      res = split.each.z(counts=counts, z=z, K=K, z.prob=t(next.z$probs), alpha=alpha, beta=beta, s=s, LLFunction="calculateLoglikFromVariables.celda_C")

+	  res = cC.splitZ(counts, m.CP.by.S, n.G.by.CP, s, z, K, nS, nG, alpha, beta, z.prob=t(next.z$probs), max.clusters.to.try=10, min.cell=3)

       logMessages(res$message, logfile=logfile, append=TRUE)

 	  # Reset convergence counter if a split occured

@@ -118,9 +118,12 @@ celda_C = function(counts, sample.label=NULL, K, alpha=1, beta=1,

       ## Re-calculate variables

       z = res$z

-      m.CP.by.S = matrix(as.integer(table(factor(z, levels=1:K), s)), ncol=nS)

-      n.G.by.CP = colSumByGroup(counts, group=z, K=K)

-      n.CP = as.integer(colSums(n.G.by.CP))

+      #m.CP.by.S = matrix(as.integer(table(factor(z, levels=1:K), s)), ncol=nS)

+      #n.G.by.CP = colSumByGroup(counts, group=z, K=K)

+      #n.CP = as.integer(colSums(n.G.by.CP))

+      m.CP.by.S = res$m.CP.by.S

+      n.G.by.CP = res$n.G.by.CP

+      n.CP = res$n.CP

     }

@@ -225,14 +228,16 @@ cC.calcEMProbZ = function(counts, m.CP.by.S, n.G.by.CP, n.by.C, n.CP, z, s, K, n

   ## Maximization to find best label for each cell

   probs = eigenMatMultInt(phi, counts) + theta[, s]  

   #probs = (t(phi) %*% counts) + theta[, s]  

+  

+  z.previous = z

   z = apply(probs, 2, which.max)

   ## Recalculate counts based on new label

-  p = cC.decomposeCounts(counts, s, z, K)

+  #p = cC.decomposeCounts(counts, s, z, K)

+  p = cC.reDecomposeCounts(counts, s, z, z.previous, n.G.by.CP, K)

   m.CP.by.S = p$m.CP.by.S

   n.G.by.CP = p$n.G.by.CP

   n.CP = p$n.CP

-  n.by.C = p$n.by.C

   return(list(m.CP.by.S=m.CP.by.S, n.G.by.CP=n.G.by.CP, n.CP=n.CP, z=z, probs=probs))

 }

@@ -400,10 +405,21 @@ cC.decomposeCounts = function(counts, s, z, K) {

   n.G.by.CP = colSumByGroup(counts, group=z, K=K)

   n.CP = as.integer(colSums(n.G.by.CP))

   n.by.C = as.integer(colSums(counts))

-  

+

   return(list(m.CP.by.S=m.CP.by.S, n.G.by.CP=n.G.by.CP, n.CP=n.CP, n.by.C=n.by.C, nS=nS, nG=nG, nM=nM))

 }

+cC.reDecomposeCounts = function(counts, s, z, previous.z, n.G.by.CP, K) {

+

+  ## Recalculate counts based on new label

+  n.G.by.CP = colSumByGroupChange(counts, n.G.by.CP, z, previous.z, K)

+  nS = length(unique(s))

+  m.CP.by.S = matrix(as.integer(table(factor(z, levels=1:K), s)), ncol=nS)

+  n.CP = as.integer(colSums(n.G.by.CP))

+

+  return(list(m.CP.by.S=m.CP.by.S, n.G.by.CP=n.G.by.CP, n.CP=n.CP))  

+}

+

 #' Calculates the conditional probability of each cell belong to each cluster given all other cluster assignments

 #'

@@ -439,22 +455,13 @@ clusterProbability.celda_C = function(celda.mod, counts, log=FALSE, ...) {

 #' @export

 calculatePerplexity.celda_C = function(counts, celda.mod, precision=128) {

-  # TODO Can try to turn into a single giant matrix multiplication by duplicating

-  #     phi / theta / sl

-  # TODO Cast to sparse matrices?

   factorized = factorizeMatrix(counts = counts, celda.mod = celda.mod, "posterior")

   theta = log(factorized$posterior$sample.states)

   phi = log(factorized$posterior$gene.states)

   sl = celda.mod$sample.label

   inner.log.prob = (t(phi) %*% counts) + theta[, sl]  

-  inner.log.prob = Rmpfr::mpfr(inner.log.prob, precision)

-  inner.log.prob.exp = exp(inner.log.prob)

-  

-  log.px = 0

-  for(i in 1:ncol(inner.log.prob.exp)) {

-    log.px = log.px + Rmpfr::asNumeric(log(sum(inner.log.prob.exp[, i])))

-  }

+  log.px = sum(apply(inner.log.prob, 2, matrixStats::logSumExp))

   perplexity = exp(-(log.px/sum(counts)))

   return(perplexity)

@@ -84,18 +84,19 @@ celda_CG = function(counts, sample.label=NULL, K, L,

   p = cCG.decomposeCounts(counts, s, z, y, K, L)

   m.CP.by.S = p$m.CP.by.S

   n.TS.by.C = p$n.TS.by.C

-  n.CP.by.TS = p$n.CP.by.TS

+  n.TS.by.CP = p$n.TS.by.CP

   n.CP = p$n.CP

   n.by.G = p$n.by.G

   n.by.C = p$n.by.C

   n.by.TS = p$n.by.TS

   nG.by.TS = p$nG.by.TS

+  n.G.by.CP = p$n.G.by.CP

   nM = p$nM

   nG = p$nG

   nS = p$nS

   rm(p)

-  ll = cCG.calcLL(K=K, L=L, m.CP.by.S=m.CP.by.S, n.CP.by.TS=n.CP.by.TS, n.by.G=n.by.G, n.by.TS=n.by.TS, nG.by.TS=nG.by.TS, nS=nS, nG=nG, alpha=alpha, beta=beta, delta=delta, gamma=gamma)

+  ll = cCG.calcLL(K=K, L=L, m.CP.by.S=m.CP.by.S, n.TS.by.CP=n.TS.by.CP, n.by.G=n.by.G, n.by.TS=n.by.TS, nG.by.TS=nG.by.TS, nS=nS, nG=nG, alpha=alpha, beta=beta, delta=delta, gamma=gamma)

   set.seed(seed)

   logMessages(date(), "... Starting celda_CG to cluster cells and genes.", logfile=logfile, append=FALSE)

@@ -106,29 +107,27 @@ celda_CG = function(counts, sample.label=NULL, K, L,

   do.gene.split = TRUE  

   while(iter <= max.iter & num.iter.without.improvement <= stop.iter) {

-    ## Gibbs sampling for each cell

-    n.TS.by.C = rowSumByGroup(counts, group=y, L=L)

-    n.TS.by.CP = t(n.CP.by.TS)

-	#next.z = cC.calcGibbsProbZ(counts=n.TS.by.C, m.CP.by.S=m.CP.by.S, n.G.by.CP=n.TS.by.CP, n.CP=n.CP, n.by.C=n.by.C, z=z, s=s, K=K, nG=L, nM=nM, alpha=alpha, beta=beta)

+    ## Gibbs or EM sampling for each cell

 	next.z = do.call(algorithm.fun, list(counts=n.TS.by.C, m.CP.by.S=m.CP.by.S, n.G.by.CP=n.TS.by.CP, n.CP=n.CP, n.by.C=n.by.C, z=z, s=s, K=K, nG=L, nM=nM, alpha=alpha, beta=beta))

     m.CP.by.S = next.z$m.CP.by.S

     n.TS.by.CP = next.z$n.G.by.CP

     n.CP = next.z$n.CP

+    n.G.by.CP = colSumByGroupChange(counts, n.G.by.CP, next.z$z, z, K)

     z = next.z$z

-

+        

     ## Gibbs sampling for each gene

-    n.CP.by.G = t(colSumByGroup(counts, group=z, K=K))

-    n.CP.by.TS = t(n.TS.by.CP) 

- 	next.y = cG.calcGibbsProbY(counts.t=n.CP.by.G, n.C.by.TS=n.CP.by.TS, n.by.TS=n.by.TS, nG.by.TS=nG.by.TS, n.by.G=n.by.G, y=y, L=L, nG=nG, beta=beta, delta=delta, gamma=gamma)

-	n.CP.by.TS = next.y$n.C.by.TS

+ 	next.y = cG.calcGibbsProbY(counts=n.G.by.CP, n.TS.by.C=n.TS.by.CP, n.by.TS=n.by.TS, nG.by.TS=nG.by.TS, n.by.G=n.by.G, y=y, L=L, nG=nG, beta=beta, delta=delta, gamma=gamma)

+	n.TS.by.CP = next.y$n.TS.by.C

 	nG.by.TS = next.y$nG.by.TS

 	n.by.TS = next.y$n.by.TS

+	n.TS.by.C = rowSumByGroupChange(counts, n.TS.by.C, next.y$y, y, L)

 	y = next.y$y

+        

     ## Perform split on i-th iteration defined by split.on.iter

 	if(K > 2 & (((iter == max.iter | num.iter.without.improvement == stop.iter) & isTRUE(split.on.last)) | (split.on.iter > 0 & iter %% split.on.iter == 0 & isTRUE(do.cell.split)))) {

 	  logMessages(date(), " ... Determining if any cell clusters should be split.", logfile=logfile, append=TRUE, sep="")

-	  res = split.each.z(counts=counts, z=z, y=y, z.prob=t(next.z$probs), K=K, L=L, alpha=alpha, delta=delta, beta=beta, gamma=gamma, s=s, LLFunction="calculateLoglikFromVariables.celda_CG")

+	  res = cCG.splitZ(counts, m.CP.by.S, n.TS.by.C, n.TS.by.CP, n.by.G, n.by.TS, nG.by.TS, s, z, K, L, nS, nG, alpha, beta, delta, gamma, z.prob=t(next.z$probs), max.clusters.to.try=10, min.cell=3)

 	  logMessages(res$message, logfile=logfile, append=TRUE)

 	  # Reset convergence counter if a split occured

@@ -141,15 +140,14 @@ celda_CG = function(counts, sample.label=NULL, K, L,

 	  ## Re-calculate variables

 	  z = res$z      

-	  m.CP.by.S = matrix(as.integer(table(factor(z, levels=1:K), s)), ncol=nS)

-	  n.TS.by.C = rowSumByGroup(counts, group=y, L=L)

-	  n.CP.by.TS = t(colSumByGroup(n.TS.by.C, group=z, K=K))

-	  n.CP = as.integer(rowSums(n.CP.by.TS))

-	  n.CP.by.G = t(colSumByGroup(counts, group=z, K=K))

+	  m.CP.by.S = res$m.CP.by.S

+	  n.TS.by.CP = res$n.TS.by.CP

+	  n.CP = res$n.CP

+	  n.G.by.CP = colSumByGroup(counts, group=z, K=K)

 	}  

 	if(L > 2 & (((iter == max.iter | num.iter.without.improvement == stop.iter) & isTRUE(split.on.last)) | (split.on.iter > 0 & iter %% split.on.iter == 0 & isTRUE(do.gene.split)))) {

 	  logMessages(date(), " ... Determining if any gene clusters should be split.", logfile=logfile, append=TRUE, sep="")

-	  res = split.each.y(counts=counts, z=z, y=y, y.prob=t(next.y$probs), K=K, L=L, alpha=alpha, beta=beta, delta=delta, gamma=gamma, s=s, LLFunction="calculateLoglikFromVariables.celda_CG")

+	  res = cCG.splitY(counts, y, m.CP.by.S, n.G.by.CP, n.TS.by.C, n.TS.by.CP, n.by.G, n.by.TS, nG.by.TS, n.CP, s, z, K, L, nS, nG, alpha, beta, delta, gamma, y.prob=t(next.y$probs), max.clusters.to.try=10, min.cell=3)

 	  logMessages(res$message, logfile=logfile, append=TRUE)

 	  # Reset convergence counter if a split occured	    

@@ -162,15 +160,14 @@ celda_CG = function(counts, sample.label=NULL, K, L,

 	  ## Re-calculate variables

 	  y = res$y        

-	  n.TS.by.C = rowSumByGroup(counts, group=y, L=L)

-	  n.CP.by.TS = t(colSumByGroup(n.TS.by.C, group=z, K=K))

-	  n.CP = as.integer(rowSums(n.CP.by.TS))

-	  n.by.TS = as.integer(rowSumByGroup(matrix(n.by.G,ncol=1), group=y, L=L))

-	  nG.by.TS = as.integer(table(factor(y, levels=1:L)))

+	  n.TS.by.CP = res$n.TS.by.CP

+	  n.by.TS = res$n.by.TS

+	  nG.by.TS = res$nG.by.TS

+	  n.TS.by.C = rowSumByGroup(counts, group=y, L=L)	  

 	}      

     ## Calculate complete likelihood

-    temp.ll = cCG.calcLL(K=K, L=L, m.CP.by.S=m.CP.by.S, n.CP.by.TS=n.CP.by.TS, n.by.G=n.by.G, n.by.TS=n.by.TS, nG.by.TS=nG.by.TS, nS=nS, nG=nG, alpha=alpha, beta=beta, delta=delta, gamma=gamma)

+    temp.ll = cCG.calcLL(K=K, L=L, m.CP.by.S=m.CP.by.S, n.TS.by.CP=n.TS.by.CP, n.by.G=n.by.G, n.by.TS=n.by.TS, nG.by.TS=nG.by.TS, nS=nS, nG=nG, alpha=alpha, beta=beta, delta=delta, gamma=gamma)

     if((all(temp.ll > ll)) | iter == 1) {

       z.best = z

       y.best = y

@@ -223,7 +220,7 @@ celda_CG = function(counts, sample.label=NULL, K, L,

 #' @param ... Other arguments

 #' @export

 simulateCells.celda_CG = function(model, S=10, C.Range=c(50,100), N.Range=c(500,5000), 

-                                  G=1000, K=3, L=10, alpha=1, beta=1, gamma=1, 

+                                  G=1000, K=3, L=10, alpha=1, beta=1, gamma=5, 

                                   delta=1, seed=12345, ...) {

   set.seed(seed)

@@ -327,7 +324,7 @@ factorizeMatrix.celda_CG = function(counts, celda.mod, type=c("counts", "proport

   nM = p$nM

   m.CP.by.S = p$m.CP.by.S

   n.TS.by.C = p$n.TS.by.C

-  n.CP.by.TS = p$n.CP.by.TS

+  n.TS.by.CP = p$n.TS.by.CP

   n.by.G = p$n.by.G

   n.by.TS = p$n.by.TS

   nG.by.TS = p$nG.by.TS

@@ -345,8 +342,8 @@ factorizeMatrix.celda_CG = function(counts, celda.mod, type=c("counts", "proport

   rownames(n.G.by.TS) = celda.mod$names$row

   rownames(m.CP.by.S) = K.names

   colnames(m.CP.by.S) = celda.mod$names$sample

-  colnames(n.CP.by.TS) = L.names

-  rownames(n.CP.by.TS) = K.names

+  colnames(n.TS.by.CP) = K.names

+  rownames(n.TS.by.CP) = L.names

   counts.list = c()

   prop.list = c()

@@ -355,7 +352,7 @@ factorizeMatrix.celda_CG = function(counts, celda.mod, type=c("counts", "proport

   if(any("counts" %in% type)) {

     counts.list = list(sample.states = m.CP.by.S,

-            				   population.states = t(n.CP.by.TS), 

+            				   population.states = n.TS.by.CP, 

             				   cell.states = n.TS.by.C,

             				   gene.states = n.G.by.TS,

             				   gene.distribution = nG.by.TS)

@@ -365,8 +362,8 @@ factorizeMatrix.celda_CG = function(counts, celda.mod, type=c("counts", "proport

     ## Need to avoid normalizing cell/gene states with zero cells/genes

     unique.z = sort(unique(z))

-    temp.n.CP.by.TS = t(n.CP.by.TS)

-    temp.n.CP.by.TS[,unique.z] = normalizeCounts(temp.n.CP.by.TS[,unique.z], scale.factor=1)

+    temp.n.TS.by.CP = n.TS.by.CP

+    temp.n.TS.by.CP[,unique.z] = normalizeCounts(temp.n.TS.by.CP[,unique.z], scale.factor=1)

     unique.y = sort(unique(y))

     temp.n.G.by.TS = n.G.by.TS

@@ -374,7 +371,7 @@ factorizeMatrix.celda_CG = function(counts, celda.mod, type=c("counts", "proport

     temp.nG.by.TS = nG.by.TS/sum(nG.by.TS)

     prop.list = list(sample.states =  normalizeCounts(m.CP.by.S, scale.factor=1),

-    				   population.states = temp.n.CP.by.TS, 

+    				   population.states = temp.n.TS.by.CP, 

     				   cell.states = normalizeCounts(n.TS.by.C, scale.factor=1),

     				   gene.states = temp.n.G.by.TS, 

     				   gene.distribution = temp.nG.by.TS)

@@ -388,7 +385,7 @@ factorizeMatrix.celda_CG = function(counts, celda.mod, type=c("counts", "proport

 	temp.nG.by.TS = (nG.by.TS + gamma)/sum(nG.by.TS + gamma)

     post.list = list(sample.states = normalizeCounts(m.CP.by.S + alpha, scale.factor=1),

-          				   population.states = normalizeCounts(t(n.CP.by.TS + beta), scale.factor=1), 

+          				   population.states = normalizeCounts(n.TS.by.CP + beta, scale.factor=1), 

           				   gene.states = gs,

           				   gene.distribution = temp.nG.by.TS)

     res = c(res, posterior = list(post.list))						    

@@ -398,7 +395,7 @@ factorizeMatrix.celda_CG = function(counts, celda.mod, type=c("counts", "proport

 }

 # Calculate the loglikelihood for the celda_CG model

-cCG.calcLL = function(K, L, m.CP.by.S, n.CP.by.TS, n.by.G, n.by.TS, nG.by.TS, nS, nG, alpha, beta, delta, gamma) {

+cCG.calcLL = function(K, L, m.CP.by.S, n.TS.by.CP, n.by.G, n.by.TS, nG.by.TS, nS, nG, alpha, beta, delta, gamma) {

   nG.by.TS[nG.by.TS == 0] = 1

   nG = sum(nG.by.TS)

@@ -413,9 +410,9 @@ cCG.calcLL = function(K, L, m.CP.by.S, n.CP.by.TS, n.by.G, n.by.TS, nG.by.TS, nS

   ## Calculate for "Phi" component

   a = K * lgamma(L * beta)

-  b = sum(lgamma(n.CP.by.TS + beta))

+  b = sum(lgamma(n.TS.by.CP + beta))

   c = -K * L * lgamma(beta)

-  d = -sum(lgamma(rowSums(n.CP.by.TS + beta)))

+  d = -sum(lgamma(colSums(n.TS.by.CP + beta)))

   phi.ll = a + b + c + d

@@ -457,7 +454,7 @@ cCG.calcLL = function(K, L, m.CP.by.S, n.CP.by.TS, n.by.G, n.by.TS, nG.by.TS, nS

 calculateLoglikFromVariables.celda_CG = function(counts, sample.label, z, y, K, L, alpha, beta, delta, gamma) {  

   s = processSampleLabels(sample.label, ncol(counts))

   p = cCG.decomposeCounts(counts, s, z, y, K, L)

-  final = cCG.calcLL(K=K, L=L, m.CP.by.S=p$m.CP.by.S, n.CP.by.TS=p$n.CP.by.TS, n.by.G=p$n.by.G, n.by.TS=p$n.by.TS, nG.by.TS=p$nG.by.TS, nS=p$nS, nG=p$nG, alpha=alpha, beta=beta, delta=delta, gamma=gamma)

+  final = cCG.calcLL(K=K, L=L, m.CP.by.S=p$m.CP.by.S, n.TS.by.CP=p$n.TS.by.CP, n.by.G=p$n.by.G, n.by.TS=p$n.by.TS, nG.by.TS=p$nG.by.TS, nS=p$nS, nG=p$nG, alpha=alpha, beta=beta, delta=delta, gamma=gamma)

   return(final)

 }

@@ -472,18 +469,18 @@ calculateLoglikFromVariables.celda_CG = function(counts, sample.label, z, y, K,

 cCG.decomposeCounts = function(counts, s, z, y, K, L) {

   nS = length(unique(s))

   m.CP.by.S = matrix(as.integer(table(factor(z, levels=1:K), s)), ncol=nS)

-    n.TS.by.C = rowSumByGroup(counts, group=y, L=L)

-  n.CP.by.TS = t(colSumByGroup(n.TS.by.C, group=z, K=K))  # TODO Refactor celda_CG

-  n.CP = as.integer(rowSums(n.CP.by.TS))

+  n.TS.by.C = rowSumByGroup(counts, group=y, L=L)

+  n.TS.by.CP = colSumByGroup(n.TS.by.C, group=z, K=K)

+  n.CP = as.integer(colSums(n.TS.by.CP))

   n.by.G = as.integer(rowSums(counts))

   n.by.C = as.integer(colSums(counts))

   n.by.TS = as.integer(rowSumByGroup(matrix(n.by.G,ncol=1), group=y, L=L))

-  nG.by.TS = as.integer(table(factor(y, 1:L)))

-  n.CP.by.G = colSumByGroup(counts, group=z, K=K)

+  nG.by.TS = tabulate(y, L)

+  n.G.by.CP = colSumByGroup(counts, group=z, K=K)

   nG = nrow(counts)

   nM = ncol(counts)

-  return(list(m.CP.by.S=m.CP.by.S, n.TS.by.C=n.TS.by.C, n.CP.by.TS=n.CP.by.TS, n.CP=n.CP, n.by.G=n.by.G, n.by.C=n.by.C, n.by.TS=n.by.TS, nG.by.TS=nG.by.TS, n.CP.by.G=n.CP.by.G, nM=nM, nG=nG, nS=nS))

+  return(list(m.CP.by.S=m.CP.by.S, n.TS.by.C=n.TS.by.C, n.TS.by.CP=n.TS.by.CP, n.CP=n.CP, n.by.G=n.by.G, n.by.C=n.by.C, n.by.TS=n.by.TS, nG.by.TS=nG.by.TS, n.G.by.CP=n.G.by.CP, nM=nM, nG=nG, nS=nS))

 }  

@@ -511,11 +508,11 @@ clusterProbability.celda_CG = function(celda.mod, counts, log=FALSE, ...) {

   p = cCG.decomposeCounts(counts, s, z, y, K, L)

   ## Gibbs sampling for each cell

-  next.z = cC.calcGibbsProbZ(counts=p$n.TS.by.C, m.CP.by.S=p$m.CP.by.S, n.G.by.CP=t(p$n.CP.by.TS), n.CP=p$n.CP, n.by.C=p$n.by.C, z=z, s=s, K=K, nG=L, nM=p$nM, alpha=alpha, beta=beta, do.sample=FALSE)

+  next.z = cC.calcGibbsProbZ(counts=p$n.TS.by.C, m.CP.by.S=p$m.CP.by.S, n.G.by.CP=p$n.TS.by.CP, n.CP=p$n.CP, n.by.C=p$n.by.C, z=z, s=s, K=K, nG=L, nM=p$nM, alpha=alpha, beta=beta, do.sample=FALSE)

   z.prob = t(next.z$probs)

   ## Gibbs sampling for each gene

-  next.y = cG.calcGibbsProbY(counts.t=p$n.CP.by.G, n.C.by.TS=p$n.CP.by.TS, n.by.TS=p$n.by.TS, nG.by.TS=p$nG.by.TS, n.by.G=p$n.by.G, y=y, L=L, nG=nG, beta=beta, delta=delta, gamma=gamma, do.sample=FALSE)

+  next.y = cG.calcGibbsProbY(counts=p$n.CP.by.G, n.TS.by.C=p$n.TS.by.CP, n.by.TS=p$n.by.TS, nG.by.TS=p$nG.by.TS, n.by.G=p$n.by.G, y=y, L=L, nG=nG, beta=beta, delta=delta, gamma=gamma, do.sample=FALSE)

   y.prob = t(next.y$probs)

   if(!isTRUE(log)) {

@@ -528,7 +525,7 @@ clusterProbability.celda_CG = function(celda.mod, counts, log=FALSE, ...) {

 #' @export

-calculatePerplexity.celda_CG = function(counts, celda.mod, precision=128) {

+calculatePerplexity.celda_CG = function(counts, celda.mod) {

   factorized = factorizeMatrix(counts = counts, celda.mod = celda.mod, "posterior")

   theta = log(factorized$posterior$sample.states)

@@ -541,14 +538,8 @@ calculatePerplexity.celda_CG = function(counts, celda.mod, precision=128) {

   eta.prob = log(eta) * nG.by.TS

   gene.by.pop.prob = log(psi %*% phi)

   inner.log.prob = (t(gene.by.pop.prob) %*% counts) + theta[, sl]  

-  inner.log.prob = Rmpfr::mpfr(inner.log.prob, precision)

-  inner.log.prob.exp = exp(inner.log.prob)

-  

-  log.px = sum(eta.prob)

-  for(i in 1:ncol(inner.log.prob.exp)) {

-    log.px = log.px + Rmpfr::asNumeric(log(sum(inner.log.prob.exp[, i])))

-  }

+  log.px = sum(apply(inner.log.prob, 2, matrixStats::logSumExp)) + sum(eta.prob)

   perplexity = exp(-(log.px/sum(counts)))

   return(perplexity)

 }

@@ -58,7 +58,6 @@ celda_G = function(counts, L, beta=1, delta=1, gamma=1,

   ## Error checking and variable processing

   counts = processCounts(counts)  

-  counts.t = t(counts)

   if(is.null(count.checksum)) {

     count.checksum = digest::digest(counts, algo="md5")

@@ -70,7 +69,7 @@ celda_G = function(counts, L, beta=1, delta=1, gamma=1,

   ## Calculate counts one time up front

   p = cG.decomposeCounts(counts=counts, y=y, L=L)

-  n.C.by.TS = p$n.C.by.TS

+  n.TS.by.C = p$n.TS.by.C

   n.by.G = p$n.by.G

   n.by.TS = p$n.by.TS

   nG.by.TS = p$nG.by.TS

@@ -82,15 +81,15 @@ celda_G = function(counts, L, beta=1, delta=1, gamma=1,

   logMessages(date(), "... Starting celda_G to cluster genes.", logfile=logfile, append=FALSE)

   ## Calculate initial log likelihood

-  ll <- cG.calcLL(n.C.by.TS=n.C.by.TS, n.by.TS=n.by.TS, n.by.G=n.by.G, nG.by.TS=nG.by.TS, nM=nM, nG=nG, L=L, beta=beta, delta=delta, gamma=gamma)

+  ll <- cG.calcLL(n.TS.by.C=n.TS.by.C, n.by.TS=n.by.TS, n.by.G=n.by.G, nG.by.TS=nG.by.TS, nM=nM, nG=nG, L=L, beta=beta, delta=delta, gamma=gamma)

   iter <- 1L

   num.iter.without.improvement = 0L

   do.gene.split = TRUE

   while(iter <= max.iter & num.iter.without.improvement <= stop.iter) {

-	next.y = cG.calcGibbsProbY(counts.t=counts.t, n.C.by.TS=n.C.by.TS, n.by.TS=n.by.TS, nG.by.TS=nG.by.TS, n.by.G=n.by.G, y=y, nG=nG, L=L, beta=beta, delta=delta, gamma=gamma)

-	n.C.by.TS = next.y$n.C.by.TS

+	next.y = cG.calcGibbsProbY(counts=counts, n.TS.by.C=n.TS.by.C, n.by.TS=n.by.TS, nG.by.TS=nG.by.TS, n.by.G=n.by.G, y=y, nG=nG, L=L, beta=beta, delta=delta, gamma=gamma)

+	n.TS.by.C = next.y$n.TS.by.C

 	nG.by.TS = next.y$nG.by.TS

 	n.by.TS = next.y$n.by.TS

 	y = next.y$y

@@ -98,7 +97,7 @@ celda_G = function(counts, L, beta=1, delta=1, gamma=1,

     ## Perform split on i-th iteration of no improvement in log likelihood

     if(L > 2 & (((iter == max.iter | num.iter.without.improvement == stop.iter) & isTRUE(split.on.last)) | (split.on.iter > 0 & iter %% split.on.iter == 0 & isTRUE(do.gene.split)))) {

       logMessages(date(), " ... Determining if any gene clusters should be split.", logfile=logfile, append=TRUE, sep="")

-      res = split.each.y(counts=counts, y=y, L=L, y.prob=t(next.y$probs), beta=beta, delta=delta, gamma=gamma, LLFunction="calculateLoglikFromVariables.celda_G")

+      res = cG.splitY(counts, y, n.TS.by.C, n.by.TS, n.by.G, nG.by.TS, nM, nG, L, beta, delta, gamma, y.prob=t(next.y$probs), min=3, max.clusters.to.try=10)

       logMessages(res$message, logfile=logfile, append=TRUE)

 	  # Reset convergence counter if a split occured	    

@@ -111,13 +110,13 @@ celda_G = function(counts, L, beta=1, delta=1, gamma=1,

       ## Re-calculate variables

       y = res$y

-      n.C.by.TS = t(rowSumByGroup(counts, group=y, L=L))

-      n.by.TS = as.integer(rowSumByGroup(matrix(n.by.G,ncol=1), group=y, L=L))

-      nG.by.TS = as.integer(table(factor(y, 1:L)))

+      n.TS.by.C = res$n.TS.by.C

+      n.by.TS = res$n.by.TS

+      nG.by.TS = res$nG.by.TS

     }

     ## Calculate complete likelihood

-    temp.ll <- cG.calcLL(n.C.by.TS=n.C.by.TS, n.by.TS=n.by.TS, n.by.G=n.by.G, nG.by.TS=nG.by.TS, nM=nM, nG=nG, L=L, beta=beta, delta=delta, gamma=gamma)

+    temp.ll <- cG.calcLL(n.TS.by.C=n.TS.by.C, n.by.TS=n.by.TS, n.by.G=n.by.G, nG.by.TS=nG.by.TS, nM=nM, nG=nG, L=L, beta=beta, delta=delta, gamma=gamma)

     if((all(temp.ll > ll)) | iter == 1) {

       y.best = y

       ll.best = temp.ll

@@ -158,13 +157,13 @@ celda_G = function(counts, L, beta=1, delta=1, gamma=1,

 # @param delta The Dirichlet distribution parameter for Eta; adds a gene pseudocount to the numbers of genes each state.

 # @param beta Vector of non-zero concentration parameters for cluster <-> gene assignment Dirichlet distribution

 # @keywords log likelihood

-cG.calcGibbsProbY = function(counts.t, n.C.by.TS, n.by.TS, nG.by.TS, n.by.G, y, L, nG, beta, delta, gamma, do.sample=TRUE) {

+cG.calcGibbsProbY = function(counts, n.TS.by.C, n.by.TS, nG.by.TS, n.by.G, y, L, nG, beta, delta, gamma, do.sample=TRUE) {

   ## Set variables up front outside of loop

   probs = matrix(NA, ncol=nG, nrow=L)

   temp.nG.by.TS = nG.by.TS 

   temp.n.by.TS = n.by.TS 

-#  temp.n.C.by.TS = n.C.by.TS

+#  temp.n.TS.by.C = n.TS.by.C

   ix <- sample(1:nG)

   for(i in ix) {

@@ -173,25 +172,25 @@ cG.calcGibbsProbY = function(counts.t, n.C.by.TS, n.by.TS, nG.by.TS, n.by.G, y,

 	nG.by.TS[y[i]] = nG.by.TS[y[i]] - 1L

 	n.by.TS[y[i]] = n.by.TS[y[i]] - n.by.G[i]

-	n.C.by.TS_1 = n.C.by.TS

-	n.C.by.TS_1[,y[i]] = n.C.by.TS_1[,y[i]] - counts.t[,i]

-    n.C.by.TS_1 = .colSums(lgamma(n.C.by.TS_1 + beta), nrow(n.C.by.TS), ncol(n.C.by.TS))

+	n.TS.by.C_1 = n.TS.by.C

+	n.TS.by.C_1[y[i],] = n.TS.by.C_1[y[i],] - counts[i,]

+    n.TS.by.C_1 = .rowSums(lgamma(n.TS.by.C_1 + beta), nrow(n.TS.by.C), ncol(n.TS.by.C))

-	n.C.by.TS_2 = n.C.by.TS

+	n.TS.by.C_2 = n.TS.by.C

 	other.ix = (1:L)[-y[i]]

-	n.C.by.TS_2[,other.ix] = n.C.by.TS_2[,other.ix] + counts.t[,i]

-    n.C.by.TS_2 = .colSums(lgamma(n.C.by.TS_2 + beta), nrow(n.C.by.TS), ncol(n.C.by.TS))

-    	

+	n.TS.by.C_2[other.ix,] = n.TS.by.C_2[other.ix,] + counts[rep(i, length(other.ix)),]

+    n.TS.by.C_2 = .rowSums(lgamma(n.TS.by.C_2 + beta), nrow(n.TS.by.C), ncol(n.TS.by.C))

+    

 	## Calculate probabilities for each state

 	for(j in 1:L) {

 	  temp.nG.by.TS = nG.by.TS 

 	  temp.n.by.TS = n.by.TS 

-#	  temp.n.C.by.TS = n.C.by.TS

+#	  temp.n.TS.by.C = n.TS.by.C

 	  temp.nG.by.TS[j] = temp.nG.by.TS[j] + 1L      

 	  temp.n.by.TS[j] = temp.n.by.TS[j] + n.by.G[i]

-#	  temp.n.C.by.TS[,j] = temp.n.C.by.TS[,j] + counts.t[,i]

+#	  temp.n.TS.by.C[j,] = temp.n.TS.by.C[j,] + counts[i,]

       pseudo.nG.by.TS = temp.nG.by.TS

@@ -201,8 +200,8 @@ cG.calcGibbsProbY = function(counts.t, n.C.by.TS, n.by.TS, nG.by.TS, n.by.G, y,

 	  other.ix = (1:L)[-j]

       probs[j,i] <- 	sum(lgamma(pseudo.nG.by.TS + gamma)) -

 						sum(lgamma(sum(pseudo.nG.by.TS + gamma))) +

-						sum(n.C.by.TS_1[other.ix]) + n.C.by.TS_2[j] +

-#						sum(lgamma(temp.n.C.by.TS + beta)) +

+						sum(n.TS.by.C_1[other.ix]) + n.TS.by.C_2[j] +

+#						sum(lgamma(temp.n.TS.by.C + beta)) +

 						sum(lgamma(pseudo.nG.by.TS * delta)) -

 						(pseudo.nG * lgamma(delta)) -

 						sum(lgamma(temp.n.by.TS + (pseudo.nG.by.TS * delta)))

@@ -214,16 +213,16 @@ cG.calcGibbsProbY = function(counts.t, n.C.by.TS, n.by.TS, nG.by.TS, n.by.G, y,

 	if(isTRUE(do.sample)) y[i] = sample.ll(probs[,i])

 	if(prev.y != y[i]) {

-	  n.C.by.TS[,prev.y] = n.C.by.TS[,prev.y] - counts.t[,i]

-	  n.C.by.TS[,y[i]] = n.C.by.TS[,y[i]] + counts.t[,i]

+	  n.TS.by.C[prev.y,] = n.TS.by.C[prev.y,] - counts[i,]

+	  n.TS.by.C[y[i],] = n.TS.by.C[y[i],] + counts[i,]

 	}

 	nG.by.TS[y[i]] = nG.by.TS[y[i]] + 1L

 	n.by.TS[y[i]] = n.by.TS[y[i]] + n.by.G[i]

-	#n.C.by.TS[,y[i]] = n.C.by.TS[,y[i]] + counts.t[,i]

+	#n.TS.by.C[y[i],] = n.TS.by.C[y[i],] + counts[i,]

   }

-  return(list(n.C.by.TS=n.C.by.TS, nG.by.TS=nG.by.TS, n.by.TS=n.by.TS, y=y, probs=probs))

+  return(list(n.TS.by.C=n.TS.by.C, nG.by.TS=nG.by.TS, n.by.TS=n.by.TS, y=y, probs=probs))

 }

@@ -244,7 +243,7 @@ cG.calcGibbsProbY = function(counts.t, n.C.by.TS, n.by.TS, nG.by.TS, n.by.G, y,

 #' @param ... Other arguments

 #' @export

 simulateCells.celda_G = function(model, C=100, N.Range=c(500,5000),  G=1000, 

-                                 L=5, beta=1, gamma=1, delta=1, seed=12345, ...) {

+                                 L=5, beta=1, gamma=5, delta=1, seed=12345, ...) {

   set.seed(seed)

   eta = rdirichlet(1, rep(gamma, L))

@@ -314,7 +313,7 @@ factorizeMatrix.celda_G = function(counts, celda.mod, type=c("counts", "proporti

   gamma = celda.mod$gamma

   p = cG.decomposeCounts(counts=counts, y=y, L=L)

-  n.TS.by.C = t(p$n.C.by.TS)

+  n.TS.by.C = p$n.TS.by.C

   n.by.G = p$n.by.G

   n.by.TS = p$n.by.TS

   nG.by.TS = p$nG.by.TS

@@ -370,16 +369,16 @@ factorizeMatrix.celda_G = function(counts, celda.mod, type=c("counts", "proporti

 # Calculate log-likelihood of celda_CG model

-cG.calcLL = function(n.C.by.TS, n.by.TS, n.by.G, nG.by.TS, nM, nG, L, beta, delta, gamma) {

+cG.calcLL = function(n.TS.by.C, n.by.TS, n.by.G, nG.by.TS, nM, nG, L, beta, delta, gamma) {

   nG.by.TS[nG.by.TS == 0] = 1

   nG = sum(nG.by.TS)

   ## Calculate for "Phi" component

   a <- nM * lgamma(L * beta)

-  b <- sum(lgamma(n.C.by.TS + beta))

+  b <- sum(lgamma(n.TS.by.C + beta))

   c <- -nM * L * lgamma(beta)

-  d <- -sum(lgamma(rowSums(n.C.by.TS + beta)))

+  d <- -sum(lgamma(colSums(n.TS.by.C + beta)))

   phi.ll <- a + b + c + d

@@ -422,7 +421,7 @@ cG.calcLL = function(n.C.by.TS, n.by.TS, n.by.G, nG.by.TS, nM, nG, L, beta, delt

 calculateLoglikFromVariables.celda_G = function(counts, y, L, beta, delta, gamma) {

   p = cG.decomposeCounts(counts=counts, y=y, L=L)

-  final <- cG.calcLL(n.C.by.TS=p$n.C.by.TS, n.by.TS=p$n.by.TS, n.by.G=p$n.by.G, nG.by.TS=p$nG.by.TS, nM=p$nM, nG=p$nG, L=L, beta=beta, delta=delta, gamma=gamma)

+  final <- cG.calcLL(n.TS.by.C=p$n.TS.by.C, n.by.TS=p$n.by.TS, n.by.G=p$n.by.G, nG.by.TS=p$nG.by.TS, nM=p$nM, nG=p$nG, L=L, beta=beta, delta=delta, gamma=gamma)

   return(final)

 }

@@ -434,14 +433,24 @@ calculateLoglikFromVariables.celda_G = function(counts, y, L, beta, delta, gamma

 #' @param L The number of clusters being considered

 cG.decomposeCounts = function(counts, y, L) {

-  n.C.by.TS = t(rowSumByGroup(counts, group=y, L=L))

-  n.by.G = as.integer(rowSums(counts))

+  n.TS.by.C = rowSumByGroup(counts, group=y, L=L)

+  n.by.G = as.integer(.rowSums(counts, nrow(counts), ncol(counts)))

   n.by.TS = as.integer(rowSumByGroup(matrix(n.by.G,ncol=1), group=y, L=L))

-  nG.by.TS = as.integer(table(factor(y, 1:L)))

+  nG.by.TS = tabulate(y, L)

   nM = ncol(counts)

   nG = nrow(counts)

-  return(list(n.C.by.TS=n.C.by.TS, n.by.G=n.by.G, n.by.TS=n.by.TS, nG.by.TS=nG.by.TS, nM=nM, nG=nG))

+  return(list(n.TS.by.C=n.TS.by.C, n.by.G=n.by.G, n.by.TS=n.by.TS, nG.by.TS=nG.by.TS, nM=nM, nG=nG))

+}

+

+

+cG.reDecomposeCounts = function(counts, y, previous.y, n.TS.by.C, n.by.G, L) {

+  ## Recalculate counts based on new label

+  n.TS.by.C = rowSumByGroupChange(counts, n.TS.by.C, y, previous.y, L)

+  n.by.TS = as.integer(rowSumByGroup(matrix(n.by.G,ncol=1), group=y, L=L))

+  nG.by.TS = tabulate(y, L)

+

+  return(list(n.TS.by.C=n.TS.by.C, n.by.TS=n.by.TS, nG.by.TS=nG.by.TS))  

 }

@@ -463,7 +472,7 @@ clusterProbability.celda_G = function(celda.mod, counts, log=FALSE, ...) {

   ## Calculate counts one time up front

   p = cG.decomposeCounts(counts=counts, y=y, L=L)

-  next.y = cG.calcGibbsProbY(counts.t=t(counts), n.C.by.TS=p$n.C.by.TS, n.by.TS=p$n.by.TS, nG.by.TS=p$nG.by.TS, n.by.G=p$n.by.G, y=y, nG=p$nG, L=L, beta=beta, delta=delta, gamma=gamma, do.sample=FALSE)  

+  next.y = cG.calcGibbsProbY(counts=counts, n.TS.by.C=p$n.TS.by.C, n.by.TS=p$n.by.TS, nG.by.TS=p$nG.by.TS, n.by.G=p$n.by.G, y=y, nG=p$nG, L=L, beta=beta, delta=delta, gamma=gamma, do.sample=FALSE)  

   y.prob = t(next.y$probs)

   if(!isTRUE(log)) {

@@ -475,7 +484,7 @@ clusterProbability.celda_G = function(celda.mod, counts, log=FALSE, ...) {

 #' @export

-calculatePerplexity.celda_G = function(counts, celda.mod, precision=128) {

+calculatePerplexity.celda_G = function(counts, celda.mod) {

   factorized = factorizeMatrix(counts = counts, celda.mod = celda.mod, "posterior")

   phi <- factorized$posterior$gene.states

@@ -252,8 +252,8 @@ logMessages = function(..., sep = " ", logfile = NULL, append = FALSE) {

 #' @export

 distinct_colors = function(n,

 						   hues = c("red", "cyan", "orange", "blue", "yellow", "purple", "green", "magenta"),

-                           saturation.range = c(0.4, 1),

-                           value.range = c(0.4, 1)) {

+                           saturation.range = c(0.7, 1),

+                           value.range = c(0.7, 1)) {

   if(!(all(hues %in% grDevices::colors()))) {

     stop("Only color names listed in the 'color' function can be used in 'hues'")

@@ -11,10 +11,9 @@

 #' 

 #' @param counts The count matrix modeled in the celdaRun parameter

 #' @param celda.mod A single celda run (usually from the _res.list_ property of a celda_list)

-#' @param precision The amount of bits of precision to pass to Rmpfr

-#' @return The perplexity for the provided chain as an mpfr number

+#' @return The perplexity for the provided data and model

 #' @export

-calculatePerplexity = function(counts, celda.mod, precision=128) {

+calculatePerplexity = function(counts, celda.mod) {

   UseMethod("calculatePerplexity", celda.mod)

 }

@@ -115,7 +115,7 @@ plotDrCluster <- function(dim1, dim2, cluster, size = 1, xlab = "Dimension_1", y

 #' @param max.iter Numeric vector; determines iterations for tsne. Default 1000.

 #' @param distance Character vector; determines which distance metric to use for tsne. Options: cosine, hellinger, spearman.

 #' @export

-celdaTsne = function(counts, celda.mod, states=NULL, perplexity=20, max.iter=2500, distance="hellinger") {

+celdaTsne = function(counts, celda.mod, states=NULL, perplexity=20, max.iter=2500, distance="hellinger", seed=12345) {

   if (!isTRUE(class(celda.mod) %in% c("celda_CG","celda_C","celda_G"))) {

     stop("celda.mod argument is not of class celda_C, celda_G or celda_CG")

   } 

@@ -148,6 +148,7 @@ celdaTsne = function(counts, celda.mod, states=NULL, perplexity=20, max.iter=250

   }

   do.pca = class(celda.mod) == "celda_C"

+  set.seed(seed)

   res = Rtsne::Rtsne(d, pca=do.pca, max_iter=max.iter, perplexity = perplexity, 

                      check_duplicates = FALSE, is_distance = TRUE)$Y

   colnames(res) = c("tsne_1", "tsne_2")

@@ -1,50 +1,146 @@

-split.each.z = function(counts, z, K, LLFunction, z.prob, min.cell=3, max.clusters.to.try = 10, ...) { 

+# cC.calcLL = function(m.CP.by.S, n.G.by.CP, s, z, K, nS, nG, alpha, beta) 

+cC.splitZ = function(counts, m.CP.by.S, n.G.by.CP, s, z, K, nS, nG, alpha, beta, z.prob, max.clusters.to.try=10, min.cell=3) {

-  dot.args = list(...)

-  

   ## Identify clusters to split

-  z.ta = table(factor(z, levels=1:K))

+  z.ta = tabulate(z, K)

   z.to.split = which(z.ta >= min.cell)

   z.non.empty = which(z.ta > 0)

+  if(length(z.to.split) == 0) {

+    m = paste0(date(), " ... Cluster sizes too small. No additional splitting was performed.") 

+    return(list(z=z, m.CP.by.S, n.G.by.CP, n.CP=n.CP, message=m))  

+  }

+  

   ## Loop through each split-able Z and perform split

-  clust.split = list()

-  for(i in 1:K) { 

-    if(i %in% z.to.split) {

-      clustLabel = suppressMessages(celda_C(counts[,z == i], K=2, max.iter=5, split.on.iter=-1, split.on.last=FALSE))

-      clust.split = c(clust.split, list(clustLabel$z))

-    } else {

-      clust.split = c(clust.split, list(NA))

-    }  

+  clust.split = vector("list", K)

+  for(i in z.to.split) { 

+    clustLabel = suppressMessages(celda_C(counts[,z == i], K=2, max.iter=5, split.on.iter=-1, split.on.last=FALSE))

+    clust.split[[i]] = clustLabel$z

   }

   ## Find second best assignment give current assignments for each cell

   z.prob[cbind(1:nrow(z.prob), z)] = NA

   z.second = apply(z.prob, 1, which.max)

+  ## Set up initial variables

+  z.split = matrix(NA, nrow=length(z), ncol=length(z.to.split) * max.clusters.to.try)

+  z.split.ll = rep(NA, ncol=length(z.to.split) * max.clusters.to.try)  

+  z.split.ll[1] = cC.calcLL(m.CP.by.S, n.G.by.CP, s, z, K, nS, nG, alpha, beta) 

+  z.split[,1] = z

+

+  ## Select worst clusters to test for reshuffling  

+  previous.z = z

   ll.shuffle = rep(NA, K)

   for(i in z.non.empty) {

     ix = z == i

     new.z = z

     new.z[ix] = z.second[ix]

-    params = c(list(counts=counts, z=new.z, K=K), dot.args)

-    ll.shuffle[i] = do.call(LLFunction, params)

-  }

+    

+    p = cC.reDecomposeCounts(counts, s, new.z, previous.z, n.G.by.CP, K)

+    n.G.by.CP = p$n.G.by.CP

+    m.CP.by.S = p$m.CP.by.S

+    ll.shuffle[i] = cC.calcLL(m.CP.by.S, n.G.by.CP, s, z, K, nS, nG, alpha, beta) 

+    previous.z = new.z

+  } 

   z.to.shuffle = head(order(ll.shuffle, decreasing = TRUE, na.last=NA), n = max.clusters.to.try)

-  if(length(z.to.split) == 0 | length(z.to.shuffle) == 0) {

+  

+  pairs = c(NA, NA)

+  split.ix = 2

+  for(i in z.to.shuffle) {

+  

+    other.clusters = setdiff(z.to.split, i)

+   

+	for(j in other.clusters) {

+	  new.z = z

+	  

+      ## Assign cluster i to the next most similar cluster (excluding cluster j) 

+      ## as defined above by the correlation      

+      ix.to.move = z == i

+      new.z[ix.to.move] = z.second[ix.to.move]

+            

+      ## Split cluster j according to the clustering defined above

+      ix.to.split = z == j

+      new.z[ix.to.split] = ifelse(clust.split[[j]] == 1, j, i)

+

+      p = cC.reDecomposeCounts(counts, s, new.z, previous.z, n.G.by.CP, K)

+      n.G.by.CP = p$n.G.by.CP

+      m.CP.by.S = p$m.CP.by.S

+      

+	  ## Calculate likelihood of split

+	  z.split.ll[split.ix] = cC.calcLL(m.CP.by.S, n.G.by.CP, s, z, K, nS, nG, alpha, beta) 

+	  z.split[,split.ix] = new.z

+	  split.ix = split.ix + 1L

+      previous.z = new.z

+      

+	  pairs = rbind(pairs, c(i, j))

+	}  

+  }

+

+  select = which.max(z.split.ll) 

+

+  if(select == 1) {

+    m = paste0(date(), " ... No additional splitting was performed.") 

+  } else {

+    m = paste0(date(), " ... Cluster ", pairs[select,1], " was reassigned and cluster ", pairs[select,2], " was split in two.")

+  } 

+

+  p = cC.reDecomposeCounts(counts, s, z.split[,select], previous.z, n.G.by.CP, K)

+  return(list(z=z.split[,select], m.CP.by.S=p$m.CP.by.S, n.G.by.CP=p$n.G.by.CP, n.CP=p$n.CP, message=m))

+}

+

+

+# cCG.calcLL = function(K, L, m.CP.by.S, n.TS.by.CP, n.by.G, n.by.TS, nG.by.TS, nS, nG, alpha, beta, delta, gamma) 

+cCG.splitZ = function(counts, m.CP.by.S, n.TS.by.C, n.TS.by.CP, n.by.G, n.by.TS, nG.by.TS, s, z, K, L, nS, nG, alpha, beta, delta, gamma, z.prob, max.clusters.to.try=10, min.cell=3) {

+

+  ## Identify clusters to split

+  z.ta = tabulate(z, K)

+  z.to.split = which(z.ta >= min.cell)

+  z.non.empty = which(z.ta > 0)

+  

+  if(length(z.to.split) == 0) {

     m = paste0(date(), " ... Cluster sizes too small. No additional splitting was performed.") 

-    return(list(z=z,message=m))

+    return(list(z=z, m.CP.by.S=m.CP.by.S, n.TS.by.CP=n.TS.by.CP, n.CP=n.CP, message=m))  

   }

-  ## Set up variables for holding results with current iteration of z

-  params = c(list(counts=counts, z=z, K=K), dot.args)

-  z.split.ll = do.call(LLFunction, params)

-  z.split = z

+  ## Loop through each split-able Z and perform split

+  clust.split = vector("list", K)

+  for(i in z.to.split) { 

+    clustLabel = suppressMessages(celda_C(counts[,z == i], K=2, max.iter=5, split.on.iter=-1, split.on.last=FALSE))

+    clust.split[[i]] = clustLabel$z

+  }

+  ## Find second best assignment give current assignments for each cell

+  z.prob[cbind(1:nrow(z.prob), z)] = NA

+  z.second = apply(z.prob, 1, which.max)

+

+  ## Set up initial variables

+  z.split = matrix(NA, nrow=length(z), ncol=length(z.to.split) * max.clusters.to.try)

+  z.split.ll = rep(NA, ncol=length(z.to.split) * max.clusters.to.try)  

+  z.split.ll[1] = cCG.calcLL(K, L, m.CP.by.S, n.TS.by.CP, n.by.G, n.by.TS, nG.by.TS, nS, nG, alpha, beta, delta, gamma) 

+  z.split[,1] = z

+

+  ## Select worst clusters to test for reshuffling  

+  previous.z = z

+  ll.shuffle = rep(NA, K)

+  for(i in z.non.empty) {

+    ix = z == i

+    new.z = z

+    new.z[ix] = z.second[ix]

+    

+    p = cC.reDecomposeCounts(n.TS.by.C, s, new.z, previous.z, n.TS.by.CP, K)

+    n.TS.by.CP = p$n.G.by.CP

+    m.CP.by.S = p$m.CP.by.S

+    ll.shuffle[i] = cCG.calcLL(K, L, m.CP.by.S, n.TS.by.CP, n.by.G, n.by.TS, nG.by.TS, nS, nG, alpha, beta, delta, gamma) 

+    previous.z = new.z

+  } 

+  z.to.shuffle = head(order(ll.shuffle, decreasing = TRUE, na.last=NA), n = max.clusters.to.try)

+

+  

   pairs = c(NA, NA)

+  split.ix = 2

   for(i in z.to.shuffle) {

-    

+  

     other.clusters = setdiff(z.to.split, i)

 	for(j in other.clusters) {

@@ -59,18 +155,116 @@ split.each.z = function(counts, z, K, LLFunction, z.prob, min.cell=3, max.cluste

       ix.to.split = z == j

       new.z[ix.to.split] = ifelse(clust.split[[j]] == 1, j, i)

+      p = cC.reDecomposeCounts(n.TS.by.C, s, new.z, previous.z, n.TS.by.CP, K)

+      n.TS.by.CP = p$n.G.by.CP

+      m.CP.by.S = p$m.CP.by.S

+      

 	  ## Calculate likelihood of split

-	  params = c(list(counts=counts, z=new.z, K=K), dot.args)

-	  new.ll = do.call(LLFunction, params)

+	  z.split.ll[split.ix] = cCG.calcLL(K, L, m.CP.by.S, n.TS.by.CP, n.by.G, n.by.TS, nG.by.TS, nS, nG, alpha, beta, delta, gamma) 

+	  z.split[,split.ix] = new.z

+	  split.ix = split.ix + 1L

+      previous.z = new.z

+      

+	  pairs = rbind(pairs, c(i, j))

+	}  

+  }

+

+  select = which.max(z.split.ll) 

-	  z.split.ll = c(z.split.ll, new.ll)

-	  z.split = cbind(z.split, new.z)

+  if(select == 1) {

+    m = paste0(date(), " ... No additional splitting was performed.") 

+  } else {

+    m = paste0(date(), " ... Cluster ", pairs[select,1], " was reassigned and cluster ", pairs[select,2], " was split in two.")

+  } 

+

+  p = cC.reDecomposeCounts(n.TS.by.C, s, z.split[,select], previous.z, n.TS.by.CP, K)

+  return(list(z=z.split[,select], m.CP.by.S=p$m.CP.by.S, n.TS.by.CP=p$n.G.by.CP, n.CP=p$n.CP, message=m))

+}

+

+

+

+cCG.splitY = function(counts, y, m.CP.by.S, n.G.by.CP, n.TS.by.C, n.TS.by.CP, n.by.G, n.by.TS, nG.by.TS, n.CP, s, z, K, L, nS, nG, alpha, beta, delta, gamma, y.prob, max.clusters.to.try=10, min.cell=3) {

+

+  ## Identify clusters to split

+  y.ta = tabulate(y, L)

+  y.to.split = which(y.ta >= min.cell)