Translated writeMixtureInfo()
This commit is contained in:
Waldir Leoncio
parent
d05fd94cfb
commit
e89c51de57
3 changed files with 396 additions and 0 deletions
|
|
@ -40,6 +40,7 @@ export(testaaPop)
|
|||
export(times)
|
||||
export(uigetfile)
|
||||
export(uiputfile)
|
||||
export(writeMixtureInfo)
|
||||
importFrom(methods,is)
|
||||
importFrom(stats,runif)
|
||||
importFrom(utils,read.delim)
|
||||
|
|
|
|||
337
R/writeMixtureInfo.R
Normal file
337
R/writeMixtureInfo.R
Normal file
|
|
@ -0,0 +1,337 @@
|
|||
#' @title Write Mixture Info
|
||||
#' @description Writes information about the mixture
|
||||
#' @param logml logml
|
||||
#' @param rowsFromInd rowsFromInd
|
||||
#' @param data data
|
||||
#' @param adjprior adjprior
|
||||
#' @param priorTerm priorTerm
|
||||
#' @param outPutFile outPutFile
|
||||
#' @param inputFile inputFile
|
||||
#' @param partitionSummary partitionSummary
|
||||
#' @param popnames popnames
|
||||
#' @param fixedK fixedK
|
||||
#' @param PARTITION PARTITION
|
||||
#' @param COUNTS COUNTS
|
||||
#' @param SUMCOUNTS SUMCOUNTS
|
||||
#' @param LOGDIFF LOGDIFF
|
||||
#' @return changesInLogml
|
||||
#' @export
|
||||
writeMixtureInfo <- function(
|
||||
logml, rowsFromInd, data, adjprior, priorTerm, outPutFile, inputFile, partitionSummary, popnames, fixedK, PARTITION, COUNTS, SUMCOUNTS,
|
||||
LOGDIFF
|
||||
) {
|
||||
|
||||
changesInLogml <- list()
|
||||
ninds <- size(data, 1) / rowsFromInd
|
||||
npops <- size(COUNTS, 3)
|
||||
# Check that the names refer to individuals
|
||||
names <- (size(popnames, 1) == ninds) #Tarkistetaan ett?nimet viittaavat yksil<69>ihin
|
||||
|
||||
if (length(outPutFile) > 0) {
|
||||
fid <- load(outPutFile)
|
||||
} else {
|
||||
fid <- -1
|
||||
message('Diverting output to baps4_output.baps')
|
||||
# TODO: replace sink with option that will record input and output
|
||||
sink('baps4_output.baps', split=TRUE) # save in text anyway.
|
||||
}
|
||||
|
||||
dispLine()
|
||||
cat('RESULTS OF INDIVIDUAL LEVEL MIXTURE ANALYSIS:')
|
||||
cat(c('Data file: ', inputFile))
|
||||
cat('Model: independent')
|
||||
cat(c('Number of clustered individuals: ', ownNum2Str(ninds)))
|
||||
cat(c('Number of groups in optimal partition: ', ownNum2Str(npops)))
|
||||
cat(c('Log(marginal likelihood) of optimal partition: ', ownNum2Str(logml)))
|
||||
cat(' ')
|
||||
if (fid != -1) {
|
||||
append(fid, 'RESULTS OF INDIVIDUAL LEVEL MIXTURE ANALYSIS:\n')
|
||||
append(fid, c('Data file: ', inputFile, '\n'))
|
||||
append(
|
||||
fid,
|
||||
c('Number of clustered individuals: ', ownNum2Str(ninds), '\n')
|
||||
)
|
||||
append(
|
||||
fid,
|
||||
c(
|
||||
'Number of groups in optimal partition: ',
|
||||
ownNum2Str(npops), '\n'
|
||||
)
|
||||
)
|
||||
append(
|
||||
fid,
|
||||
c(
|
||||
'Log(marginal likelihood) of optimal partition: ',
|
||||
ownNum2Str(logml),
|
||||
'\n'
|
||||
)
|
||||
)
|
||||
}
|
||||
|
||||
cluster_count <- length(unique(PARTITION))
|
||||
cat('Best Partition: ')
|
||||
if (fid != -1) {
|
||||
append(fid, c('Best Partition: ', '\n'))
|
||||
}
|
||||
for (m in 1:cluster_count) {
|
||||
indsInM <- find(PARTITION == m)
|
||||
length_of_beginning <- 11 + floor(log10(m))
|
||||
cluster_size <- length(indsInM)
|
||||
|
||||
if (names) {
|
||||
text <- c(
|
||||
'Cluster ',
|
||||
as.character(m),
|
||||
': {',
|
||||
as.character(popnames[[indsInM[1]]])
|
||||
)
|
||||
for (k in 2:cluster_size) {
|
||||
text <- c(text, ', ', as.character(popnames[[indsInM[k]]]))
|
||||
}
|
||||
} else {
|
||||
text <- c(
|
||||
'Cluster ', as.character(m), ': {', as.character(indsInM[1])
|
||||
)
|
||||
for (k in 2:cluster_size) {
|
||||
text <- c(text, ', ', as.character(indsInM[k]))
|
||||
}
|
||||
}
|
||||
text <- c(text, '}')
|
||||
while (length(text) > 58) {
|
||||
# Take one line and display it.
|
||||
new_line <- takeLine(text, 58)
|
||||
text <- (length(new_line) + 1):end
|
||||
cat(new_line)
|
||||
if (fid != -1) {
|
||||
append(fid, new_line)
|
||||
append(fid,'\n')
|
||||
}
|
||||
if (length(text) > 0) {
|
||||
text <- c(blanks(length_of_beginning), text)
|
||||
} else {
|
||||
text <- ""
|
||||
}
|
||||
}
|
||||
if (text != "") {
|
||||
cat(text)
|
||||
if (fid != -1) {
|
||||
append(fid, text)
|
||||
append(fid,'\n')
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (npops > 1) {
|
||||
cat(' ')
|
||||
cat(' ')
|
||||
cat(
|
||||
'Changes in log(marginal likelihood)',
|
||||
' if indvidual i is moved to group j:'
|
||||
)
|
||||
if (fid != -1) {
|
||||
append(fid, ' ')
|
||||
append(fid, '\n')
|
||||
append(fid, ' ')
|
||||
append(fid, '\n')
|
||||
append(
|
||||
fid,
|
||||
c(
|
||||
'Changes in log(marginal likelihood)',
|
||||
'if indvidual i is moved to group j:'
|
||||
)
|
||||
)
|
||||
append(fid, '\n')
|
||||
}
|
||||
|
||||
if (names) {
|
||||
nameSizes <- zeros(ninds, 1)
|
||||
for (i in 1:ninds) {
|
||||
nimi <- as.character(popnames[i])
|
||||
nameSizes[i] <- length(nimi)
|
||||
}
|
||||
maxSize <- max(nameSizes)
|
||||
maxSize <- max(maxSize, 5)
|
||||
erotus <- maxSize - 5
|
||||
alku <- blanks(erotus)
|
||||
ekarivi <- c(alku, ' ind', blanks(6 + erotus))
|
||||
} else {
|
||||
ekarivi <- ' ind '
|
||||
}
|
||||
for (i in 1:cluster_count) {
|
||||
ekarivi <- c(ekarivi, ownNum2Str(i), blanks(8 - floor(log10(i))))
|
||||
}
|
||||
cat(ekarivi)
|
||||
if (fid != -1) {
|
||||
append(fid, ekarivi)
|
||||
append(fid, '\n')
|
||||
}
|
||||
|
||||
# %ninds = size(data,1)/rowsFromInd;
|
||||
changesInLogml <- t(LOGDIFF)
|
||||
for (ind in 1:ninds) {
|
||||
muutokset <- changesInLogml[, ind]
|
||||
|
||||
if (names) {
|
||||
nimi <- as.character(popnames[ind])
|
||||
rivi <- c(blanks(maxSize - length(nimi)), nimi, ':')
|
||||
} else {
|
||||
rivi <- c(blanks(4 - floor(log10(ind))), ownNum2Str(ind), ':')
|
||||
}
|
||||
for (j in 1:npops) {
|
||||
rivi <- c(rivi, ' ', logml2String(omaRound(muutokset[j])))
|
||||
}
|
||||
cat(rivi)
|
||||
if (fid != -1) {
|
||||
append(fid, rivi)
|
||||
append(fid, '\n')
|
||||
}
|
||||
}
|
||||
|
||||
cat(' ')
|
||||
cat(' ')
|
||||
cat('KL-divergence matrix in PHYLIP format:')
|
||||
|
||||
dist_mat <- zeros(npops, npops)
|
||||
if (fid != -1) {
|
||||
append(fid, ' ')
|
||||
append(fid, ' ')
|
||||
append(fid, c('KL-divergence matrix in PHYLIP format:'))
|
||||
append(fid, '\n')
|
||||
}
|
||||
|
||||
maxnoalle <- size(COUNTS, 1)
|
||||
nloci <- size(COUNTS, 2)
|
||||
d <- zeros(maxnoalle, nloci, npops)
|
||||
prior <- adjprior
|
||||
prior[find(prior == 1)] <- 0
|
||||
nollia <- find(all(prior == 0)) # Loci in which only one allele was detected.
|
||||
prior[1, nollia] <- 1
|
||||
for (pop1 in 1:npops) {
|
||||
d[, , pop1] <- (squeeze(COUNTS[, , pop1]) + prior) /
|
||||
repmat(sum(squeeze(COUNTS[, , pop1]) + prior), c(maxnoalle, 1))
|
||||
}
|
||||
ekarivi <- as.character(npops)
|
||||
cat(ekarivi)
|
||||
if (fid != -1) {
|
||||
append(fid, ekarivi)
|
||||
append(fid, '\n')
|
||||
}
|
||||
|
||||
for (pop1 in 1:npops) {
|
||||
for (pop2 in 1:(pop1 - 1)) {
|
||||
dist1 <- d[, , pop1]
|
||||
dist2 <- d[, , pop2]
|
||||
div12 <- sum(
|
||||
sum(dist1 * log2((dist1 + 10 ^ -10) / (dist2 + 10 ^ -10)))
|
||||
) / nloci
|
||||
div21 <- sum(
|
||||
sum(dist2 * log2((dist2 + 10 ^ -10) / (dist1 + 10 ^ -10)))
|
||||
) / nloci
|
||||
div <- (div12 + div21) / 2
|
||||
dist_mat(pop1, pop2) <- div
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
dist_mat <- dist_mat + t(dist_mat) # make it symmetric
|
||||
for (pop1 in 1:npops) {
|
||||
rivi <- c('Cluster_', as.character(pop1), ' ')
|
||||
for (pop2 in 1:npops) {
|
||||
rivi <- c(rivi, kldiv2str(dist_mat[pop1, pop2]), ' ')
|
||||
}
|
||||
cat(rivi)
|
||||
if (fid != -1) {
|
||||
append(fid, rivi)
|
||||
append(fid, '\n')
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
cat(' ')
|
||||
cat(' ');
|
||||
cat(
|
||||
'List of sizes of 10 best visited partitions',
|
||||
'and corresponding log(ml) values'
|
||||
)
|
||||
|
||||
if (fid != -1) {
|
||||
append(fid, ' ')
|
||||
append(fid, '\n')
|
||||
append(fid, ' ')
|
||||
append(fid, '\n')
|
||||
append(
|
||||
fid,
|
||||
c(
|
||||
'List of sizes of 10 best visited partitions',
|
||||
'and corresponding log(ml) values'
|
||||
)
|
||||
)
|
||||
append(fid, '\n')
|
||||
}
|
||||
|
||||
partitionSummary <- sortrows(partitionSummary, 2)
|
||||
partitionSummary <- partitionSummary[size(partitionSummary, 1):1, ]
|
||||
partitionSummary <- partitionSummary[find(partitionSummary[, 2] > -1e49), ]
|
||||
if (size(partitionSummary, 1) > 10) {
|
||||
vikaPartitio <- 10
|
||||
} else {
|
||||
vikaPartitio <- size(partitionSummary, 1)
|
||||
}
|
||||
for (part in 1:vikaPartitio) {
|
||||
line <- c(
|
||||
as.character(partitionSummary[part, 1]),
|
||||
' ',
|
||||
as.character(partitionSummary(part, 2))
|
||||
)
|
||||
cat(line)
|
||||
if (fid != -1) {
|
||||
append(fid, line)
|
||||
append(fid, '\n')
|
||||
}
|
||||
}
|
||||
|
||||
if (!fixedK) {
|
||||
cat(' ')
|
||||
cat(' ')
|
||||
cat('Probabilities for number of clusters')
|
||||
|
||||
if (fid != -1) {
|
||||
append(fid, ' ')
|
||||
append(fid, '\n')
|
||||
append(fid, ' ')
|
||||
append(fid, '\n')
|
||||
append(fid, c('Probabilities for number of clusters'))
|
||||
append(fid, '\n')
|
||||
}
|
||||
|
||||
npopsTaulu <- unique(partitionSummary[, 1])
|
||||
len <- length(npopsTaulu)
|
||||
probs <- zeros(len, 1)
|
||||
partitionSummary[, 2] <- partitionSummary[, 2] -
|
||||
max(partitionSummary[, 2])
|
||||
sumtn <- sum(exp(partitionSummary[, 2]))
|
||||
for (i in 1:len) {
|
||||
npopstn <- sum(
|
||||
exp(
|
||||
partitionSummary[find(
|
||||
partitionSummary[, 1] == npopsTaulu[i]
|
||||
), 2]
|
||||
)
|
||||
)
|
||||
probs[i] <- npopstn / sumtn
|
||||
}
|
||||
for (i in 1:len) {
|
||||
if (probs[i] > 1e-5) {
|
||||
line <- c(
|
||||
as.character(npopsTaulu[i]), ' ', as.character(probs[i])
|
||||
)
|
||||
cat(line)
|
||||
if (fid != -1) {
|
||||
append(fid, line)
|
||||
append(fid, '\n')
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
return(changesInLogml)
|
||||
}
|
||||
58
man/writeMixtureInfo.Rd
Normal file
58
man/writeMixtureInfo.Rd
Normal file
|
|
@ -0,0 +1,58 @@
|
|||
% Generated by roxygen2: do not edit by hand
|
||||
% Please edit documentation in R/writeMixtureInfo.R
|
||||
\name{writeMixtureInfo}
|
||||
\alias{writeMixtureInfo}
|
||||
\title{Write Mixture Info}
|
||||
\usage{
|
||||
writeMixtureInfo(
|
||||
logml,
|
||||
rowsFromInd,
|
||||
data,
|
||||
adjprior,
|
||||
priorTerm,
|
||||
outPutFile,
|
||||
inputFile,
|
||||
partitionSummary,
|
||||
popnames,
|
||||
fixedK,
|
||||
PARTITION,
|
||||
COUNTS,
|
||||
SUMCOUNTS,
|
||||
LOGDIFF
|
||||
)
|
||||
}
|
||||
\arguments{
|
||||
\item{logml}{logml}
|
||||
|
||||
\item{rowsFromInd}{rowsFromInd}
|
||||
|
||||
\item{data}{data}
|
||||
|
||||
\item{adjprior}{adjprior}
|
||||
|
||||
\item{priorTerm}{priorTerm}
|
||||
|
||||
\item{outPutFile}{outPutFile}
|
||||
|
||||
\item{inputFile}{inputFile}
|
||||
|
||||
\item{partitionSummary}{partitionSummary}
|
||||
|
||||
\item{popnames}{popnames}
|
||||
|
||||
\item{fixedK}{fixedK}
|
||||
|
||||
\item{PARTITION}{PARTITION}
|
||||
|
||||
\item{COUNTS}{COUNTS}
|
||||
|
||||
\item{SUMCOUNTS}{SUMCOUNTS}
|
||||
|
||||
\item{LOGDIFF}{LOGDIFF}
|
||||
}
|
||||
\value{
|
||||
changesInLogml
|
||||
}
|
||||
\description{
|
||||
Writes information about the mixture
|
||||
}
|
||||
Loading…
Add table
Reference in a new issue