Nothing
# This file is part of the TimeMachine.
# Copyright (C) 2013 Gianluca Campanella <gianluca@campanella.org>
#
# The TimeMachine is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# The TimeMachine is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE.
# See the GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along with
# theTimeMachine. If not, see <http://www.gnu.org/licenses/>.
setGeneric("pophistory", def=function(x, probs=seq(0, 1, 0.25)) UseMethod("pophistory"))
setGeneric("ttc", def=function(x) UseMethod("ttc"))
setGeneric("ttm", def=function(x) UseMethod("ttm"))
full.transitions <- function(unitary.transitions, loci) {
unitary.transitions <- as.matrix(unitary.transitions)
loci <- as.integer(loci)
.Call("compute_full_transitions", unitary.transitions, loci)
}
stationary.dist <- function(transitions) {
transitions <- as.matrix(transitions)
.Call("compute_stationary_distribution", transitions)
}
sample.pop <- function(transitions, pi=NULL, pop.size, mu) {
transitions <- as.matrix(transitions)
if (!is.null(pi)) {
pi <- as.numeric(pi)
} else {
pi <- stationary.dist(transitions)
}
pop.size <- as.integer(pop.size)
mu <- as.numeric(mu)
types <- nrow(transitions)
population <- integer(types)
population[sample(types, 1, prob=pi)] <- 2
while (sum(population) < pop.size) {
# Sample ancestor type
ancestor.type <- sample(types, 1, prob=population)
# Sample type of event
p.mutation <- mu / (sum(population) - 1 + mu)
u <- runif(1)
if (u < p.mutation) {
# Mutation
offspring.type <- sample(types, 1, prob=transitions[,ancestor.type])
population[ancestor.type] <- population[ancestor.type] - 1
population[offspring.type] <- population[offspring.type] + 1
} else {
# Split
population[ancestor.type] <- population[ancestor.type] + 1
}
}
population
}
tm <- function(transitions, pi=NULL, population, n, mu, samples, threads=NULL) {
transitions <- as.matrix(transitions)
if (!is.null(pi)) {
pi <- as.numeric(pi)
} else {
pi <- stationary.dist(transitions)
}
population <- as.integer(population)
n <- as.integer(n)
mu <- as.numeric(mu)
samples <- as.integer(samples)
if (!is.null(threads)) {
threads <- as.integer(threads)
}
structure(.Call("estimate_loglik", transitions, pi, population, n, mu,
samples, threads), class="tm")
}
mu.mle <- function(transitions, pi=NULL, population, n, mu.int, samples,
threads=NULL, ...) {
transitions <- as.matrix(transitions)
if (!is.null(pi)) {
pi <- as.numeric(pi)
} else {
pi <- stationary.dist(transitions)
}
population <- as.integer(population)
n <- as.integer(n)
mu.int <- as.numeric(mu.int)
samples <- as.integer(samples)
threads <- as.integer(threads)
opt.res <- optimize(function(mu) {
mean(tm(transitions, pi, population, n, mu, samples, threads))
}, mu.int, maximum=TRUE, ...)
list(mu.hat=opt.res$maximum, loglik=opt.res$objective)
}
density.tm <- function(x, ...) {
density(x$logliks, ...)
}
mean.tm <- function(x, ...) {
mean(x$logliks)
}
pophistory.tm <- function(x, probs=seq(0, 1, 0.25)) {
h <- t(apply(x$coalescent.events, 1, quantile, probs=probs))
rownames(h) <- seq(sum(x$population) - 1, max(2, x$n))
h
}
print.tm <- function(x, ...) {
summary(x)
}
quantile.tm <- function(x, ...) {
quantile(x$logliks, ...)
}
summary.tm <- function(object, ..., digits=max(3, getOption("digits")-3)) {
n.percentage <- object$n / sum(object$population)# 100
mean.sim.time <- mean(object$simulation.times)
mean.loglik <- mean(object)
loglik.quantiles <- quantile(object, probs=c(0.05, 0.95))
cat("Parameters\n")
cat(sprintf(" Mutation rate: %.*f\n", digits, object$mu))
cat(sprintf(" Target population size: %d (%.*f%%)\n", object$n, digits,
n.percentage))
cat(sprintf(" Samples: %d\n", length(object$logliks)))
cat("Performance\n")
cat(sprintf(" Mean simulation time: %.*f\n", digits, mean.sim.time))
cat(sprintf(" Total computation time: %.*f\n", digits, object$total.time))
cat("Log-likelihood\n")
cat(sprintf(" Mean: %.*f\n", digits, mean(object)))
cat(sprintf(" 5%%-95%% quantiles: %.*f %.*f\n", digits, loglik.quantiles[1],
digits, loglik.quantiles[2]))
}
ttc.tm <- function(x) {
ttc.table <- tabulate(apply(x$coalescent.events, 2, diff))
structure(ttc.table / sum(ttc.table), names=1:length(ttc.table))
}
ttm.tm <- function(x) {
ttm.table <- tabulate(do.call(c, apply(x$coalescent.events, 2, function(y) {
diff(setdiff(0:y[length(y)], y))
})))
structure(ttm.table / sum(ttm.table), names=1:length(ttm.table))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.