.Rhistory

vecs <- list()
for (i in 1:n) {
vecs[[i]] <- tree.vec.match(trees[[i]],lambda,labelmatch[[i]],k,n)
centre <- centre + vecs[[i]]*likes[[i]]
}
centre <- centre/n
d <- list()
for (i in 1:n){
v <- vecs[[i]]-centre
d[[i]] <- sqrt(sum(v^2))
}
class(d) <- "numeric"
md <- min(d)
median <- which(d==md)
result <- list()
result$centre <- centre
result$median <- median
result$mindist <- md
return(result)
}
med.tree <- cmpfun(med.tree)
library(ape)
library(Rcpp)
library(inline)
computeKCMetricVectorCPP <- function(tree, lambda=0, return_lambda_function=F) {
if(lambda<0 || lambda>1) stop("Pick lambda in [0,1]")
num_leaves <- length(tree$tip.label)
num_edges <- nrow(tree$edge)
tip_order <- match(1:num_leaves, order(tree$tip.label))
edge_order <- order(tree$edge[,1], decreasing=T)
edges <- tree$edge[edge_order,]
edge_lengths <- tree$edge.length[edge_order]
annotated_nodes <- list()
for(i in 1:num_edges) {
parent <- edges[i,1]
child <- edges[i,2]
if(child <= num_leaves) {
child <- tip_order[child]
annotated_nodes[[child]] <- list(root_distance=NULL, edges_to_root=1, partitions=list(child))
}
aggregated_partitions <- annotated_nodes[[child]]$partitions[[1]]
if((child > num_leaves)) {
for(p in 2:length(annotated_nodes[[child]]$partitions))
aggregated_partitions <- c(aggregated_partitions, annotated_nodes[[child]]$partitions[[p]])
}
annotated_nodes[[child]]$root_distance <- edge_lengths[i]
if(parent > length(annotated_nodes) || is.null(annotated_nodes[[parent]])) {
annotated_nodes[[parent]] <- list(root_distance=NULL, edges_to_root=1, partitions=list(aggregated_partitions))
}
else {
annotated_nodes[[parent]]$partitions[[length(annotated_nodes[[parent]]$partitions)+1]] <- aggregated_partitions
}
}
annotated_nodes[[num_leaves+1]]$root_distance <- 0
annotated_nodes[[num_leaves+1]]$edges_to_root <- 0
for(i in num_edges:1) {
parent <- edges[i,1]
child <- edges[i,2]
if(child <= num_leaves)
child <- tip_order[child]
annotated_nodes[[child]]$root_distance <- annotated_nodes[[child]]$root_distance + annotated_nodes[[parent]]$root_distance
annotated_nodes[[child]]$edges_to_root <- annotated_nodes[[child]]$edges_to_root + annotated_nodes[[parent]]$edges_to_root
}
vector_length <- (num_leaves*(num_leaves-1)/2) + num_leaves
length_root_distances <- double(vector_length)
topological_root_distances <- integer(vector_length)
topological_root_distances[(vector_length-num_leaves+1):vector_length] <- 1
length_root_distances[(vector_length-num_leaves+1):vector_length] <- edge_lengths[match(1:num_leaves, edges[,2])][order(tree$tip.label)]
index_offsets <- c(0, cumsum((num_leaves-1):1))
sapply(annotated_nodes, function(node) {
if(length(node$partitions) > 1 && node$root_distance > 0) {
num_groups <- length(node$partitions)
for(group_a in 1:(num_groups-1)) {
for(group_b in (group_a+1):num_groups) {
CPP_update_combinations(length_root_distances, topological_root_distances, node$partitions[[group_a]],
node$partitions[[group_b]], index_offsets, node$root_distance, node$edges_to_root)
}
}
}
})
if(!return_lambda_function)
return(lambda * length_root_distances + (1-lambda) * topological_root_distances)
else {
return(function(l) {
if(l<0 || l>1) stop("Pick lambda in [0,1]")
return(l * length_root_distances + (1-l) * topological_root_distances) })
}
}
CPP_update_combinations <- cppFunction("void updateDistancesWithCombinations(NumericVector& length_root_distances,
NumericVector& topological_root_distances,
IntegerVector& left_partition,
IntegerVector& right_partition,
IntegerVector& index_offsets,
double distance_to_root,
int edges_to_root)
{
// Iterate through all combinations.
for(int i=0; i < left_partition.size(); ++i) {
for(int j=0; j < right_partition.size(); ++j) {
int first_leaf = left_partition[i];
int second_leaf = right_partition[j];
// Because of the symmetric distances.
if(left_partition[i] > right_partition[j]) {
first_leaf = right_partition[j];
second_leaf = left_partition[i];
}
// Roll the index (notice we take into account C++ indices here, starting at 0).
int combination_index = index_offsets[first_leaf-1] + (second_leaf - first_leaf) - 1;
// Update the vectors.
length_root_distances[combination_index] = distance_to_root;
topological_root_distances[combination_index] = edges_to_root;
}
}
}")
computeKCTreeDistance <- function(tree_a, tree_b, lambda=0, return_lambda_function=F) {
metric_a <- computeKCMetricVectorCPP(tree_a, lambda, return_lambda_function)
metric_b <- computeKCMetricVectorCPP(tree_b, lambda, return_lambda_function)
if(!return_lambda_function) {
return(sqrt(sum((metric_a - metric_b)^2)))
}
else {
return(function(l) {
return(sqrt(sum((metric_a(l) - metric_b(l))^2)))
})
}
}
computeKCTreeDistances <- function(trees, lambda=0, return_lambda_function=F, save_memory=F) {
num_trees <- length(trees)
if(!return_lambda_function) {
distances <- matrix(0.0, num_trees, num_trees)
if(!save_memory) {
tree_metrics <- t(sapply(trees, function(tree) {computeKCMetricVectorCPP(tree, lambda, F)}))
sapply(1:(num_trees-1), function(i) {
sapply((i+1):num_trees, function(j) {
distances[i,j] <<- distances[j,i] <<- sqrt(sum((tree_metrics[i,] - tree_metrics[j,])^2))
})
})
}
else {
sapply(1:(num_trees-1), function(i) {
sapply((i+1):num_trees, function(j) {
distances[i,j] <<- distances[j,i] <<- computeKCTreeDistance(trees[[i]], trees[[j]], lambda, F)
})
})
}
return(as.dist(distances))
}
else {
if(save_memory)
warning("save_memory=T is incompatible with return_lambda_function=T, setting save_memory=F")
tree_metric_functions <- sapply(trees, function(tree) {computeKCMetricVectorCPP(tree, lambda, T)})
compute_distance_matrix_function <- function(l) {
distances <- matrix(0.0, num_trees, num_trees)
sapply(1:(num_trees-1), function(i) {
sapply((i+1):num_trees, function(j) {
distances[i,j] <<- distances[j,i] <<- sqrt(sum((tree_metric_functions[[i]](l) - tree_metric_functions[[j]](l))^2))
})
})
return(as.dist(distances))
}
return(compute_distance_matrix_function)
}
}
testMetricVectorMethods <- function(num_trees=100, num_leaves=100, seed=NULL, lambda=NULL) {
if(!is.null(seed))
set.seed(seed)
else {
seed <- sample(1:10000,1)
set.seed(seed)
}
print(paste("Testing metric vector methods, seed:", seed), quote=F)
for(i in 1:num_trees) {
tree <- di2multi(rtree(num_leaves), runif(1, 0, 0.8))
if(is.null(lambda))
l <- runif(1)
else
l <- lambda
true_vector <- tree.vec(tree, l)
optimized_cpp_vector <- computeKCMetricVectorCPP(tree,l)
if(any(abs(true_vector - optimized_cpp_vector) > 10e-10))
stop("Incoherence in methods results")
}
print("Metric vector test passed with flying colors...", quote=F)
}
testDistanceMethods <- function(num_trees=100, num_leaves=100, seed=NULL, lambda=NULL) {
if(!is.null(seed))
set.seed(seed)
else {
seed <- sample(1:10000,1)
set.seed(seed)
}
print(paste("Testing distance methods, seed:", seed), quote=F)
for(i in 1:num_trees) {
tree_a <- di2multi(rtree(num_leaves), runif(1, 0, 0.8))
tree_b <- di2multi(rtree(num_leaves), runif(1, 0, 0.8))
if(is.null(lambda))
l <- runif(1)
else
l <- lambda
true_dist <- tree.dist(tree_a, tree_b, l)
optimized_cpp_dist <- computeKCTreeDistance(tree_a, tree_b, l)
if(abs(true_dist - optimized_cpp_dist) > 10e-10)
stop("Incoherence in methods results")
}
print("Distance test passed with flying colors...", quote=F)
}
testMultiDistanceMethods <- function(num_trees=100, num_leaves=100, seed=NULL, lambda=NULL) {
if(!is.null(seed))
set.seed(seed)
else {
seed <- sample(1:10000,1)
set.seed(seed)
}
print(paste("Testing multiple distance methods, seed:", seed), quote=F)
trees <- rmtree(num_trees, num_leaves)
if(is.null(lambda))
l <- runif(1)
else
l <- lambda
true_dist <- multi.dist(trees, l)
optimized_cpp_dist <- computeKCTreeDistances(trees, l, save_memory = F)
if(any(abs(true_dist - optimized_cpp_dist) > 10e-10))
stop("Incoherence in methods results")
print("Multipe Distance test passed with flying colors...", quote=F)
}
testMetricVectorMethods(num_trees=10, num_leaves=100)
testDistanceMethods(num_trees=100, num_leaves=100)
testMultiDistanceMethods(num_trees=100, num_leaves=50)
tree <- rtree(2500)
print("Time for original metric vector method: 2500 leaves, lambda 0.5", quote=F)
system.time(tree.vec(tree,0.5))
print("Time for C++ optimized metric vector method: 2500 leaves, lambda 0.5", quote=F)
system.time(computeKCMetricVectorCPP(tree,0.5))
tree_a <- di2multi(rtree(2500),runif(1, 0, 0.8))
tree_b <- di2multi(rtree(2500),runif(1, 0, 0.8))
print("Time for original distance method: 2500 leaves, lambda 0.5", quote=F)
system.time(tree.dist(tree_a, tree_b, 0.5))
print("Time for C++ optimized distance method: 2500 leaves, lambda 0.5", quote=F)
system.time(computeKCTreeDistance(tree_a, tree_b, 0.5))
trees <- rmtree(100, 500)
print("Time for original multi-distance method: 100 trees, 500 leaves, lambda 0.5", quote=F)
system.time(multi.dist(trees, lambda=0.5))
print("Time for the C++ optimized multi-distance method: 100 trees, 500 leaves, lambda 0.5", quote=F)
system.time(computeKCTreeDistances(trees, lambda = 0.5, return_lambda_function = F))
require(phangorn)
system.time(computeKCTreeDistance(tree_a, tree_b, 0.5))
system.time(tree.vec(tree,0.5))
tree.vec(tree,0.5)
require(ape)
tree.vec(tree,0.5)
tr1 <- rtree(6)
computeKCMetricVectorCPP(tr1)
tree.vec(tr1)
pen.edge.tree <- function(tree,k) {tree$edge[match(1:k, tree$edge[,2]),] }
pen.edge.tree <- cmpfun(pen.edge.tree)
pen.edge.treematch  <- function(tree,labelmatch) {tree$edge[match(labelmatch, tree$edge[,2]),] }
pen.edge.treematch <- cmpfun(pen.edge.treematch)
tree.vec <- function(tr1,lambda=0,type="number") { # allow output type to be number or function
if (type=="number"){
if (lambda<0) {stop("Pick lambda in [0,1]")}
if (lambda>1) {stop("Pick lambda in [0,1]")}
k <- length(tr1$tip.label)
if (lambda!=0) { # if lambda=0 then we don't need edge lengths to be defined, but if lambda!=0 then we do
if (is.null(tr1$edge.length)) {
stop("edge lengths not defined")
}
}
M1 <- linear.mrca(tr1,k); # kxk MRCA matrix for tr1
pairs <- combn2(1:k)
tiporder <- order(tr1$tip.label)
if (lambda!=1){ # make a copy with edge lengths = 1 because we need to know topological distances
TR1 <- tr1;       TR1$edge.length <- rep(1,length(tr1$edge.length))
D1 <- dist.nodes(TR1);
}
if (lambda!=0) { # if lambda!=0 we need to know branch length distances
d1 <- dist.nodes(tr1);
}
if (lambda==1) { vt <- rep(0,k*(k-1)/2)}
else {
vt <- apply(pairs, 1, function(x) D1[k+1,M1[[tiporder[[x[1]]],tiporder[[x[2]]]]]])
}
vt <- as.numeric(c(vt,rep(1,k)))
if (lambda==0) { vl <- rep(0,k*(k+1)/2) }
else {
vl <- apply(pairs, 1, function(x) d1[k+1,M1[[tiporder[[x[1]]],tiporder[[x[2]]]]]])
ep1 <- pen.edge.treematch(tr1,tiporder);
pen.length1 <- apply(ep1, 1, function(x) d1[x[1],x[2]])
vl <- as.numeric(c(vl,pen.length1))
}
v <- (1-lambda)*vt + lambda*vl
return(v)
}
if (type=="function") {
lambda <- integer()
k <- length(tr1$tip.label)
if (is.null(tr1$edge.length)) {
stop("edge lengths not defined")
}
M1 <- linear.mrca(tr1,k); # kxk MRCA matrix for tree 1
pairs <- combn2(1:k)
tiporder <- order(tr1$tip.label)
TR1 <- tr1
TR1$edge.length <- rep(1,length(tr1$edge.length));
D1 <- dist.nodes(TR1);
d1 <- dist.nodes(tr1);
vt <- apply(pairs, 1, function(x) D1[k+1,M1[[tiporder[[x[1]]],tiporder[[x[2]]]]]])
vl <- apply(pairs, 1, function(x) d1[k+1,M1[[tiporder[[x[1]]],tiporder[[x[2]]]]]])
ep1 <- pen.edge.treematch(tr1,tiporder);
pen.length1 <- apply(ep1, 1, function(x) d1[x[1],x[2]])
vlambda <- function(lambda) {
if (lambda<0) {stop("Pick lambda in [0,1]")}
if (lambda>1) {stop("Pick lambda in [0,1]")}
(c(((1-lambda)*vt + lambda*vl),(lambda*pen.length1))) }
return(vlambda)
}
}
tree.vec <- cmpfun(tree.vec)
tree.vec(tree,0.5)
tree
tree.vec(tr1)
tree <- rtree(2500)
tree.vec(tree,0.5)
tree.vec(tr1)
install_github("ThibautJombart/treescape",args='-l "C:/icnas3.cc.ic.ac.uk/mlkendal/R/win-library/3.1"')
library(devtools)
install_github("ThibautJombart/treescape",args='-l "C:/icnas3.cc.ic.ac.uk/mlkendal/R/win-library/3.1"')
library(treescape, lib.loc="C:/icnas3.cc.ic.ac.uk/mlkendal/R/win-library/3.1")
tree.vec(tr1)
library(shiny)
?shiny
install.packages(c("adegenet", "ape", "dendextend", "devtools", "e1071", "goftest", "httpuv", "httr", "igraph", "phangorn", "phylobase", "plyr", "polyclip", "pracma", "R6", "RCurl", "rgl", "scales", "shiny", "stringr", "XML"))
library(treescape, lib.loc="C:/icnas3.cc.ic.ac.uk/mlkendal/R/win-library/3.1")
treescapeServer()
suppressPackageStartupMessages(library(installr))
RStudio_CRAN_data_folder <- download_RStudio_CRAN_data(START = '2015-09-10',
END = '2015-10-19',
log_folder = "C:/Users/mlkendal/Dropbox/CRANlogs")
my_RStudio_CRAN_data <- read_RStudio_CRAN_data(RStudio_CRAN_data_folder)
my_RStudio_CRAN_data <- format_RStudio_CRAN_data(my_RStudio_CRAN_data)
barplot_package_users_per_day("treescape", my_RStudio_CRAN_data)
lineplot_package_downloads(pkg_names = c("treescape","phyloTop"),
dataset = my_RStudio_CRAN_data)
lineplot_package_downloads(pkg_names = c("treescape","kdetrees"),
dataset = my_RStudio_CRAN_data)
suppressPackageStartupMessages(library(installr))
RStudio_CRAN_data_folder <- download_RStudio_CRAN_data(START = '2015-09-10',
END = '2016-01-11',
log_folder = "C:/Users/mlkendal/Dropbox/CRANlogs")
my_RStudio_CRAN_data <- read_RStudio_CRAN_data(RStudio_CRAN_data_folder)
my_RStudio_CRAN_data <- format_RStudio_CRAN_data(my_RStudio_CRAN_data)
barplot_package_users_per_day("outbreaker", my_RStudio_CRAN_data)
barplot_package_users_per_day("treescape", my_RStudio_CRAN_data)
lineplot_package_downloads(pkg_names = c("treescape","kdetrees"),
dataset = my_RStudio_CRAN_data)
lineplot_package_downloads(pkg_names = c("treescape","rtwy"),
dataset = my_RStudio_CRAN_data)
lineplot_package_downloads(pkg_names = c("treescape","outbreaker"),
dataset = my_RStudio_CRAN_data)
setwd("C:/Users/mlkendal/Dropbox/treescape")
library(devtools)
install_github("juba/scatterD3",args='-l "C:/icnas3.cc.ic.ac.uk/mlkendal/R/win-library/3.2"')
library(scatterD3, lib.loc="C:/icnas3.cc.ic.ac.uk/mlkendal/R/win-library/3.2")
install.packages("digest")
install_github("juba/scatterD3",args='-l "C:/icnas3.cc.ic.ac.uk/mlkendal/R/win-library/3.2"')
library(devtools)
install_github("juba/scatterD3",args='-l "C:/icnas3.cc.ic.ac.uk/mlkendal/R/win-library/3.2"')
devtools::load_all(".")
devtools::load_all(".")
install.packages("roxygen2")
devtools::load_all(".")
install.packages("testthat")
devtools::load_all(".")
install.packages("Rcpp")
install.packages("ape")
check()
q()
n
check()
q()
n
knit("README.Rmd")
q()
n