Improved handling of NULL objects (#25)

R/computePopulationLogml.R

@@ -4,9 +4,13 @@ computePopulationLogml <- function(pops, adjprior, priorTerm) {
   # ======================================================== #
   # Limiting COUNTS size                                     #
   # ======================================================== #
-  COUNTS <- COUNTS[
+  if (!is.null(adjprior)) {
-    seq_len(nrow(adjprior)), seq_len(ncol(adjprior)), pops, drop = FALSE
+    nr <- seq_len(nrow(adjprior))
-  ]
+    nc <- seq_len(ncol(adjprior))
+    COUNTS <- COUNTS[nr, nc, pops, drop = FALSE]
+  } else {
+    COUNTS <- NA
+  }
 
   x <- size(COUNTS, 1)
   y <- size(COUNTS, 2)
@@ -15,17 +19,15 @@ computePopulationLogml <- function(pops, adjprior, priorTerm) {
   # ======================================================== #
   # Computation                                              #
   # ======================================================== #
+  term1 <- NULL
+  if (!is.null(adjprior)) {
     isarray <- length(dim(repmat(adjprior, c(1, 1, length(pops))))) > 2
     term1 <- squeeze(
       sum(
         sum(
           reshape(
             lgamma(
-            repmat(adjprior, c(1, 1, length(pops))) +
+              repmat(adjprior, c(1, 1, length(pops))) + COUNTS[nr, nc, pops, drop = !isarray]
-              COUNTS[
-                seq_len(nrow(adjprior)), seq_len(ncol(adjprior)), pops,
-                drop = !isarray
-              ]
             ),
             c(x, y, z)
           ),
@@ -34,6 +36,7 @@ computePopulationLogml <- function(pops, adjprior, priorTerm) {
         2
       )
     )
+  }
   if (is.null(priorTerm)) priorTerm <- 0
   popLogml <- term1 - sum(lgamma(1 + SUMCOUNTS[pops, ]), 2) - priorTerm
   return(popLogml)

R/greedyMix.R

@@ -13,14 +13,24 @@
 #' data <- system.file("extdata", "FASTA_clustering_haploid.fasta", package = "rBAPS")
 #' greedyMix(data, "fasta")
 greedyMix <- function(
-  data, format, partitionCompare, ninds, rowsFromInd, noalle,
+  data, format, partitionCompare = NULL, ninds = NULL, rowsFromInd = NULL,
-  adjprior, priorTerm, alleleCodesinp, popnames, fixedK = FALSE,
+  noalle = NULL, adjprior = NULL, npops = 1L, priorTerm = NULL,
-  partition_compare = FALSE, verbose = TRUE
+  alleleCodes = NULL, inp = NULL, popnames = NULL, fixedK = FALSE, verbose = TRUE
 ) {
   # Importing and handling data ================================================
   data <- importFile(data, format, verbose)
+  c <- list(
+    # TODO: get elements from handleData()?
+    noalle = noalle,
+    rows = NA,
+    data = data,
+    adjprior = adjprior,
+    priorTerm = priorTerm,
+    rowsFromInd = rowsFromInd
+  )
 
-  if (partition_compare) {
+  # Comparing partitions =======================================================
+  if (!is.null(partitionCompare)) {
     logmls <- comparePartitions(
       data, nrows(data), partitionCompare[["partitions"]], ninds, rowsFromInd,
       noalle, adjprior
@@ -29,11 +39,20 @@ greedyMix <- function(
 
 
   # Generating partition summary ===============================================
-  logml_npops_partitionSummary <- indMixWrapper(data, npops, fixedK);
+  logml_npops_partitionSummary <- indMixWrapper(c, npops, fixedK, verbose);
   logml <- logml_npops_partitionSummary[["logml"]]
   npops <- logml_npops_partitionSummary[["npops"]]
   partitionSummary <- logml_npops_partitionSummary[["partitionSummary"]]
-  stopifnot(logml != 1)
+
+  # Generating output object ===================================================
+  out <- list(
+      "alleleCodes" = alleleCodes, "adjprior" = adjprior, "popnames" = popnames,
+      "rowsFromInd" = rowsFromInd, "data" = data, "npops" = npops,
+      "noalle" = noalle, "mixtureType" = "mix", "logml" = logml
+    )
+  if (logml == 1) {
+    return(out)
+  }
 
   # Writing mixture info =======================================================
   changesInLogml <- writeMixtureInfo(
@@ -41,14 +60,6 @@ greedyMix <- function(
     popnames, fixedK
   )
 
-  # Returning results ==========================================================
+  # Updateing results ==========================================================
-  return(
+  return(c(out, "changesInLogml" = changesInLogml))
-    list(
-      "alleleCodes" = alleleCodes, "adjprior" = adjprior, "popnames" = popnames,
-      "rowsFromInd" = rowsFromInd, "data" = data, "npops" = npops,
-      "noalle" = noalle, "mixtureType" = "mix", "logml" = logml,
-      "changesInLogml" = changesInLogml
-    )
-  )
-
 }

R/indMix.R

@@ -1,4 +1,4 @@
-indMix <- function(c, npops, dispText = TRUE) {
+indMix <- function(c, npops, counts = NULL, sumcounts = NULL, max_iter = 100L, dispText = FALSE) {
   # Greedy search algorithm with unknown number of classes for regular
   # clustering.
   # Input npops is not used if called by greedyMix or greedyPopMix.
@@ -17,8 +17,11 @@ indMix <- function(c, npops, dispText = TRUE) {
   if (isfield(c, "dist")) {
     dist <- c$dist
     Z <- c$Z
+  } else {
+    Z <- NULL
   }
 
+
   rm(c)
   nargin <- length(as.list(match.call())) - 1
   if (nargin < 2) {
@@ -84,6 +87,8 @@ indMix <- function(c, npops, dispText = TRUE) {
       )
     }
     ninds <- size(rows, 1)
+
+    if (!is.null(Z)) {
       initialPartition <- admixture_initialization(initData, npops, Z)
       sumcounts_counts_logml <- initialCounts(
         initialPartition, data, npops, rows, noalle, adjprior
@@ -97,12 +102,12 @@ indMix <- function(c, npops, dispText = TRUE) {
         apu <- rows[i]
         PARTITION[i] <- initialPartition[apu[1]]
       }
+    }
 
     COUNTS <- counts
     SUMCOUNTS <- sumcounts
-    POP_LOGML <- computePopulationLogml(1:npops, adjprior, priorTerm)
+    POP_LOGML <- computePopulationLogml(seq_len(npops), adjprior, priorTerm)
     LOGDIFF <- repmat(-Inf, c(ninds, npops))
-    rm(initialPartition, counts, sumcounts)
 
     # PARHAAN MIXTURE-PARTITION ETSIMINEN
     nRoundTypes <- 7

R/indMixWrapper.R

@@ -1,7 +1,7 @@
-indMixWrapper <- function(c, npops, fixedK = FALSE) {
+indMixWrapper <- function(c, npops, fixedK = FALSE, verbose = FALSE) {
   if (fixedK) {
     stop("indMix_fixK() not yet implemented.") # TODO: translate indMix_fixK.m
   } else {
-    indMix(c, npops, TRUE)
+    indMix(c, npops, verbose)
   }
 }