diff --git a/DESCRIPTION b/DESCRIPTION index 695be07..9d7ef04 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: rBAPS Title: Bayesian Analysis of Population Structure -Version: 0.0.0.9029 +Version: 0.0.0.9030 Date: 2020-11-09 Authors@R: c( diff --git a/NEWS.md b/NEWS.md index f09f421..421c424 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,7 @@ # rBAPS (development version) +* Fixed `greedyMix()` for BAPS and Genepop files. + # rBAPS 0.0.0.9021 * Added a `NEWS.md` file to track changes to the package. diff --git a/R/greedyMix.R b/R/greedyMix.R index cea602b..0ac6ff1 100644 --- a/R/greedyMix.R +++ b/R/greedyMix.R @@ -80,7 +80,7 @@ greedyMix <- function( # Generating partition summary =============================================== ekat <- seq(1L, ninds * c[["rowsFromInd"]], c[["rowsFromInd"]]) c[["rows"]] <- cbind(ekat, ekat + c[["rowsFromInd"]] - 1L) - logml_npops_partitionSummary <- indMixWrapper(c, npops, counts, sumcounts, max_iter, fixedK, verbose) # FIXME: not working for FASTA, GenePop + logml_npops_partitionSummary <- indMixWrapper(c, npops, counts, sumcounts, max_iter, fixedK, verbose) # FIXME: not working for FASTA logml <- logml_npops_partitionSummary[["logml"]] npops <- logml_npops_partitionSummary[["npops"]] partitionSummary <- logml_npops_partitionSummary[["partitionSummary"]] diff --git a/R/kldiv2str.R b/R/kldiv2str.R index a73ac78..4a21b4d 100644 --- a/R/kldiv2str.R +++ b/R/kldiv2str.R @@ -20,5 +20,5 @@ kldiv2str <- function(div, max_chars = 6L) { mjono[3] <- "." mjono[2] <- palautaYks(abs(div), suurinYks) } - return(mjono) + paste(mjono, collapse = "") } diff --git a/R/laskeLoggis.R b/R/laskeLoggis.R index bb6ced9..52eae31 100644 --- a/R/laskeLoggis.R +++ b/R/laskeLoggis.R @@ -2,8 +2,6 @@ laskeLoggis <- function(counts, sumcounts, adjprior) { npops <- size(counts, 3) replicated_adjprior <- array(adjprior, c(nrow(adjprior), ncol(adjprior), npops)) sum1 <- sum(sum(sum(lgamma(counts + replicated_adjprior)))) - sum3 <- sum(sum(lgamma(adjprior))) - sum(sum(lgamma(1 + sumcounts))) - logml2 <- sum1 - npops * sum3 - loggis <- logml2 - return(loggis) + sum2 <- npops * sum(sum(lgamma(adjprior))) + sum(sum(lgamma(1 + sumcounts))) + sum1 - sum2 } diff --git a/R/laskeMuutokset12345.R b/R/laskeMuutokset12345.R index 78c9f5c..e940324 100644 --- a/R/laskeMuutokset12345.R +++ b/R/laskeMuutokset12345.R @@ -423,7 +423,7 @@ greedyMix_muutokset <- R6Class( globals$SUMCOUNTS[i1, ] <- globals$SUMCOUNTS[i1, ] + diffInSumCounts if (i1 < npops) { - i2 <- c(1:(i1 - 1), (i1 + 1):npops) + i2 <- c(seq_len(i1 - 1), (i1 + 1):npops) } else { i2 <- 1:(i1 - 1) } @@ -458,12 +458,12 @@ greedyMix_muutokset <- R6Class( muutokset <- zeros(npops2, npops) i1_logml <- globals$POP_LOGML[i1] - for (pop2 in 1:npops2) { + for (pop2 in seq_len(npops2)) { inds <- inds2[matlab2r::find(T2 == pop2)] ninds <- length(inds) if (ninds > 0) { rows <- list() - for (i in 1:ninds) { + for (i in seq_len(ninds)) { ind <- inds[i] lisa <- globalRows[ind, 1]:globalRows[ind, 2] rows <- c(rows, t(lisa)) @@ -480,7 +480,7 @@ greedyMix_muutokset <- R6Class( globals$SUMCOUNTS[i1, ] <- globals$SUMCOUNTS[i1, ] + diffInSumCounts if (i1 < npops) { - i2 <- c(1:(i1 - 1), (i1 + 1):npops) + i2 <- c(seq_len(i1 - 1), (i1 + 1):npops) } else { i2 <- 1:(i1 - 1) } diff --git a/R/laskeOsaDist.R b/R/laskeOsaDist.R index b6c8570..e915276 100644 --- a/R/laskeOsaDist.R +++ b/R/laskeOsaDist.R @@ -11,7 +11,7 @@ laskeOsaDist <- function(inds2, dist, ninds) { ninds2 <- length(inds2) apu <- zeros(choose(ninds2, 2), 2) rivi <- 1 - for (i in 1:(ninds2 - 1)) { + for (i in seq_len(ninds2 - 1)) { for (j in (i + 1):ninds2) { apu[rivi, 1] <- inds2[i] apu[rivi, 2] <- inds2[j] @@ -20,6 +20,6 @@ laskeOsaDist <- function(inds2, dist, ninds) { } apu <- (apu[, 1] - 1) * ninds - apu[, 1] / 2 * (apu[, 1] - 1) + (apu[, 2] - apu[, 1]) - dist2 <- dist(apu) + dist2 <- dist[apu] return(dist2) } diff --git a/R/writeMixtureInfo.R b/R/writeMixtureInfo.R index 79294e5..6514454 100644 --- a/R/writeMixtureInfo.R +++ b/R/writeMixtureInfo.R @@ -65,9 +65,9 @@ writeMixtureInfo <- function( } cluster_count <- length(unique(globals$PARTITION)) - if (verbose) cat("Best Partition: ") + if (verbose) cat("\nBest Partition:\n") if (fid != -1) { - append(fid, c("Best Partition: ", "\n")) + append(fid, c("\nBest Partition:", "\n")) } for (m in 1:cluster_count) { indsInM <- matlab2r::find(globals$PARTITION == m) @@ -116,7 +116,7 @@ writeMixtureInfo <- function( cat("\n") cat( "Changes in log(marginal likelihood)", - " if indvidual i is moved to group j:\n" + "if indvidual i is moved to group j:\n" ) } if (fid != -1) { @@ -163,12 +163,12 @@ writeMixtureInfo <- function( muutokset <- changesInLogml[, ind] if (names) { nimi <- as.character(popnames[ind]) - rivi <- c(blanks(maxSize - length(nimi)), nimi, ":\n") + rivi <- c(blanks(maxSize - length(nimi)), nimi, ":") } else { - rivi <- c("\n", blanks(4 - floor(log10(ind))), ownNum2Str(ind), ":\n") + rivi <- c("\n", blanks(4 - floor(log10(ind))), ownNum2Str(ind), ":") } for (j in 1:npops) { - rivi <- c(rivi, logml2String(omaRound(muutokset[j]), "")) + rivi <- c(rivi, " ", logml2String(omaRound(muutokset[j]))) } if (verbose) cat(rivi, sep = "") if (fid != -1) { @@ -176,11 +176,11 @@ writeMixtureInfo <- function( append(fid, "\n") } } - if (verbose) cat("\n\nKL-divergence matrix in PHYLIP format: ") + if (verbose) cat("\n\nKL-divergence matrix in PHYLIP format:\n") dist_mat <- zeros(npops, npops) if (fid != -1) { - append(fid, " KL-divergence matrix in PHYLIP format: ") + append(fid, " KL-divergence matrix in PHYLIP format:\n") } globals$COUNTS <- globals$COUNTS[seq_len(nrow(adjprior)), seq_len(ncol(adjprior)), , drop = FALSE] @@ -196,7 +196,10 @@ writeMixtureInfo <- function( prior[1, nollia] <- 1 for (pop1 in 1:npops) { squeezed_COUNTS_prior <- squeeze(globals$COUNTS[, , pop1]) + prior - d[, , pop1] <- squeezed_COUNTS_prior / sum(squeezed_COUNTS_prior) + repmat_squeezed_COUNTS_prior <- matlab2r::repmat( + colSums(squeezed_COUNTS_prior), c(maxnoalle, 1L) + ) + d[, , pop1] <- squeezed_COUNTS_prior / repmat_squeezed_COUNTS_prior } ekarivi <- as.character(npops) if (verbose) cat(ekarivi) @@ -223,9 +226,9 @@ writeMixtureInfo <- function( dist_mat <- dist_mat + t(dist_mat) # make it symmetric for (pop1 in 1:npops) { - rivi <- c("\nCluster_", as.character(pop1), "\n") + rivi <- c("\nCluster_", as.character(pop1), " ") for (pop2 in 1:npops) { - rivi <- c(rivi, kldiv2str(dist_mat[pop1, pop2], 4L)) + rivi <- c(rivi, kldiv2str(dist_mat[pop1, pop2])) } if (verbose) cat(rivi, sep = "") if (fid != -1) { diff --git a/tests/testthat/test-greedyMix.R b/tests/testthat/test-greedyMix.R index 2c80ca3..402af13 100644 --- a/tests/testthat/test-greedyMix.R +++ b/tests/testthat/test-greedyMix.R @@ -65,13 +65,20 @@ test_that("greedyMix() fails successfully", { }) test_that("greedyMix() works when it should", { + tol <- 1e-4 baps_file <- file.path(path_inst, "BAPS_clustering_diploid.txt") - genepop_file <- file.path(path_inst, "GenePop.txt") - fasta_file <- file.path(path_inst, "FASTA_clustering_haploid.fasta") greedy_baps <- greedyMix(baps_file, "BAPS") - expect_error(greedy_genepop <- greedyMix(genepop_file, "GenePop")) # TEMP: fails expect_type(greedy_baps, "list") expect_length(greedy_baps, 10L) + expect_equal(greedy_baps[['logml']], -78.10151, tolerance = tol) + + genepop_file <- file.path(path_inst, "GenePop.txt") + greedy_genepop <- greedyMix(genepop_file, "GenePop") + expect_type(greedy_genepop, "list") + expect_length(greedy_genepop, 10L) + expect_equal(greedy_genepop[['logml']], -10.76224, tolerance = tol) + + fasta_file <- file.path(path_inst, "FASTA_clustering_haploid.fasta") }) context("Linkage")