evolve.phylo.R
## evolve.tree.R (2005-12-04)
## Character Simulation under a Brownian Model
## Copyright 2005 Julien Dutheil
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
evolve.phylo <- function(phy, value, var) {
if (!("phylo" %in% class(phy)))
stop("object \"phy\" is not of class \"phylo\"")
if (is.null(phy$edge.length))
stop("tree \" phy\" must have branch lengths.")
nchar <- max(length(value), length(var))
value <- rep(value, length=nchar)
var <- rep(var, length=nchar)
char.names <- names(value);
## added by EP for the new coding of "phylo" (2006-10-04):
phy <- new2old.phylo(phy)
## End
edges <- phy$edge
nodes <- unique(as.vector(edges))
n <- length(nodes) # Number of nodes
root <- match("-1", nodes)
states<-list();
states[["-1"]] <- value
for(node in nodes[-root]) {
edge.index <- match(node, edges[,2])
edge.length <- phy$edge.length[edge.index]
ancestor <- edges[edge.index, 1]
ancestor.node.index <- match(ancestor, nodes)
ancestor.states <- states[[ancestor.node.index]]
index <- match(node, nodes)
x <- numeric(nchar)
for(i in 1:nchar) {
x[i] <- rnorm(1, mean=ancestor.states[i], sd=sqrt(var[i]*edge.length))
}
states[[index]] <- x;
}
nodes.states <- as.data.frame(matrix(ncol=nchar, nrow=0))
if(!is.null(char.names)) names(nodes.states) <- char.names
count <- 1
for(i in unique(edges[,1])) {
nodes.states[i,] <- states[[match(i, nodes)]]
count <- count + 1
}
nl <- length(phy$tip.label) #Number of leaves
leaves.states <- as.data.frame(matrix(ncol=nchar, nrow=0))
if(!is.null(char.names)) names(leaves.states) <- char.names
count <- 1
for(i in 1:nl) {
leaves.states[as.character(count),] <- states[[match(as.character(i), nodes)]]
count <- count + 1
}
phy[["node.character"]] <- nodes.states;
phy[["tip.character"]] <- leaves.states;
if(! "ancestral" %in% class(phy)) class(phy) <- c("ancestral", class(phy));
## added by EP for the new coding of "phylo" (2006-10-04):
phy <- old2new.phylo(phy)
## End
return(phy)
}