R/evolve.phylo.R

   1 ## evolve.tree.R (2005-12-04)
   2
   3 ##   Character Simulation under a Brownian Model
   4
   5 ## Copyright 2005 Julien Dutheil
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 evolve.phylo <- function(phy, value, var) {
  11   if (!("phylo" %in% class(phy)))
  12       stop("object \"phy\" is not of class \"phylo\"")
  13   if (is.null(phy$edge.length))
  14       stop("tree \" phy\" must have branch lengths.")
  15   nchar <- max(length(value), length(var))
  16   value <- rep(value, length=nchar)
  17   var   <- rep(var,   length=nchar)
  18   char.names <- names(value);
  19
  20   ## added by EP for the new coding of "phylo" (2006-10-04):
  21   phy <- new2old.phylo(phy)
  22   ## End
  23   edges <- phy$edge
  24   nodes <- unique(as.vector(edges))
  25   n <- length(nodes) # Number of nodes
  26   root <- match("-1", nodes)
  27   states<-list();
  28   states[["-1"]] <- value
  29   for(node in nodes[-root]) {
  30     edge.index <- match(node, edges[,2])
  31     edge.length <- phy$edge.length[edge.index]
  32     ancestor <- edges[edge.index, 1]
  33     ancestor.node.index <- match(ancestor, nodes)
  34     ancestor.states <- states[[ancestor.node.index]]
  35     index <- match(node, nodes)
  36     x <- numeric(nchar)
  37     for(i in 1:nchar) {
  38       x[i] <- rnorm(1, mean=ancestor.states[i], sd=sqrt(var[i]*edge.length))
  39     }
  40     states[[index]] <- x;
  41   }
  42   nodes.states <- as.data.frame(matrix(ncol=nchar, nrow=0))
  43   if(!is.null(char.names)) names(nodes.states) <- char.names
  44   count <- 1
  45   for(i in unique(edges[,1])) {
  46     nodes.states[i,] <- states[[match(i, nodes)]]
  47     count <- count + 1
  48   }
  49
  50   nl <- length(phy$tip.label) #Number of leaves
  51   leaves.states <- as.data.frame(matrix(ncol=nchar, nrow=0))
  52   if(!is.null(char.names)) names(leaves.states) <- char.names
  53   count <- 1
  54   for(i in 1:nl) {
  55     leaves.states[as.character(count),] <- states[[match(as.character(i), nodes)]]
  56     count <- count + 1
  57   }
  58
  59   phy[["node.character"]] <- nodes.states;
  60   phy[["tip.character"]]  <- leaves.states;
  61   if(! "ancestral" %in% class(phy)) class(phy) <- c("ancestral", class(phy));
  62   ## added by EP for the new coding of "phylo" (2006-10-04):
  63   phy <- old2new.phylo(phy)
  64   ## End
  65   return(phy)
  66 }