-## 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)
-}