diff --git a/R/computePopulationLogml.R b/R/computePopulationLogml.R index 32abb3f..9114a87 100644 --- a/R/computePopulationLogml.R +++ b/R/computePopulationLogml.R @@ -4,9 +4,13 @@ computePopulationLogml <- function(pops, adjprior, priorTerm) { # ======================================================== # # Limiting COUNTS size # # ======================================================== # - COUNTS <- COUNTS[ - seq_len(nrow(adjprior)), seq_len(ncol(adjprior)), pops, drop = FALSE - ] + if (!is.null(adjprior)) { + 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,25 +19,24 @@ computePopulationLogml <- function(pops, adjprior, priorTerm) { # ======================================================== # # Computation # # ======================================================== # - isarray <- length(dim(repmat(adjprior, c(1, 1, length(pops))))) > 2 - term1 <- squeeze( - sum( + term1 <- NULL + if (!is.null(adjprior)) { + isarray <- length(dim(repmat(adjprior, c(1, 1, length(pops))))) > 2 + term1 <- squeeze( sum( - reshape( - lgamma( - repmat(adjprior, c(1, 1, length(pops))) + - COUNTS[ - seq_len(nrow(adjprior)), seq_len(ncol(adjprior)), pops, - drop = !isarray - ] + sum( + reshape( + lgamma( + repmat(adjprior, c(1, 1, length(pops))) + COUNTS[nr, nc, pops, drop = !isarray] + ), + c(x, y, z) ), - c(x, y, z) + 1 ), - 1 - ), - 2 + 2 + ) ) - ) + } if (is.null(priorTerm)) priorTerm <- 0 popLogml <- term1 - sum(lgamma(1 + SUMCOUNTS[pops, ]), 2) - priorTerm return(popLogml) diff --git a/R/greedyMix.R b/R/greedyMix.R index 28c9815..cbc6912 100644 --- a/R/greedyMix.R +++ b/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, - adjprior, priorTerm, alleleCodesinp, popnames, fixedK = FALSE, - partition_compare = FALSE, verbose = TRUE + data, format, partitionCompare = NULL, ninds = NULL, rowsFromInd = NULL, + noalle = NULL, adjprior = NULL, npops = 1L, priorTerm = NULL, + 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 ========================================================== - return( - list( - "alleleCodes" = alleleCodes, "adjprior" = adjprior, "popnames" = popnames, - "rowsFromInd" = rowsFromInd, "data" = data, "npops" = npops, - "noalle" = noalle, "mixtureType" = "mix", "logml" = logml, - "changesInLogml" = changesInLogml - ) - ) - + # Updateing results ========================================================== + return(c(out, "changesInLogml" = changesInLogml)) } diff --git a/R/indMix.R b/R/indMix.R index 883d19a..3bfcd83 100644 --- a/R/indMix.R +++ b/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,25 +87,27 @@ indMix <- function(c, npops, dispText = TRUE) { ) } ninds <- size(rows, 1) - initialPartition <- admixture_initialization(initData, npops, Z) - sumcounts_counts_logml <- initialCounts( - initialPartition, data, npops, rows, noalle, adjprior - ) - sumcounts <- sumcounts_counts_logml$sumcounts - counts <- sumcounts_counts_logml$counts - logml <- sumcounts_counts_logml$logml - PARTITION <- zeros(ninds, 1) - for (i in 1:ninds) { - apu <- rows[i] - PARTITION[i] <- initialPartition[apu[1]] + if (!is.null(Z)) { + initialPartition <- admixture_initialization(initData, npops, Z) + sumcounts_counts_logml <- initialCounts( + initialPartition, data, npops, rows, noalle, adjprior + ) + sumcounts <- sumcounts_counts_logml$sumcounts + counts <- sumcounts_counts_logml$counts + logml <- sumcounts_counts_logml$logml + + PARTITION <- zeros(ninds, 1) + for (i in 1:ninds) { + 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 diff --git a/R/indMixWrapper.R b/R/indMixWrapper.R index 256b6c4..1153c2c 100644 --- a/R/indMixWrapper.R +++ b/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) } }