Improved handling of NULL objects (#25)
This commit is contained in:
Waldir Leoncio
parent
843a929938
commit
b034158e2b
4 changed files with 69 additions and 50 deletions
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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))
|
||||
}
|
||||
|
|
|
|||
33
R/indMix.R
33
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
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
}
|
||||
}
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue