R/rTrait.R

   1 ## rTrait.R (2010-02-03)
   2
   3 ##   Trait Evolution
   4
   5 ## Copyright 2010 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 rTraitDisc <-
  11     function(phy, model = "ER", k = if (is.matrix(model)) ncol(model) else 2,
  12              rate = 0.1, states = LETTERS[1:k], freq = rep(1/k, k),
  13              ancestor = FALSE, root.value = 1)
  14 {
  15     if (is.null(phy$edge.length))
  16         stop("tree has no branch length")
  17     if (any(phy$edge.length < 0))
  18         stop("at least one branch length negative")
  19
  20     if (is.character(model)) {
  21         switch(model, "ER" = {
  22                    if (length(rate) != 1)
  23                        stop("`rate' must have one element")
  24                    Q <- matrix(rate, k, k)
  25                }, "ARD" = {
  26                    if (length(rate) != k*(k - 1))
  27                        stop("`rate' must have k(k - 1) elements")
  28                    Q <- matrix(0, k, k)
  29                    Q[col(Q) != row(Q)] <- rate
  30                }, "SYM" = {
  31                    if (length(rate) != k*(k - 1)/2)
  32                        stop("`rate' must have k(k - 1)/2 elements")
  33                    Q <- matrix(0, k, k)
  34                    sel <- col(Q) < row(Q)
  35                    Q[sel] <- rate
  36                    Q <- t(Q)
  37                    Q[sel] <- rate
  38                })
  39     }
  40     if (is.matrix(model)) {
  41         Q <- model
  42         if (ncol(Q) != nrow(Q))
  43             stop("the matrix given as `model' must be square")
  44     }
  45
  46     phy <- reorder(phy, "pruningwise")
  47     n <- length(phy$tip.label)
  48     N <- dim(phy$edge)[1]
  49     ROOT <- n + 1L
  50     x <- integer(n + phy$Nnode)
  51     x[ROOT] <- as.integer(root.value)
  52
  53     anc <- phy$edge[, 1]
  54     des <- phy$edge[, 2]
  55     el <- phy$edge.length
  56
  57     if (is.function(model)) {
  58         environment(model) <- environment() # to find 'k'
  59         for (i in N:1) x[des[i]] <- model(x[anc[i]], el[i])
  60     } else {
  61         diag(Q) <- -rowSums(Q)
  62         freq <- rep(freq, each = k)
  63         for (i in N:1) {
  64             p <- matexpo(Q * freq * el[i])[x[anc[i]], ]
  65             x[des[i]] <- .Internal(sample(k, size = 1, FALSE, prob = p))
  66         }
  67     }
  68
  69     if (ancestor) {
  70         if (is.null(phy$node.label)) phy <- makeNodeLabel(phy)
  71         names(x) <- c(phy$tip.label, phy$node.label)
  72     } else {
  73         x <- x[1:n]
  74         names(x) <- phy$tip.label
  75     }
  76     class(x) <- "factor"
  77     levels(x) <- states
  78     x
  79 }
  80
  81 rTraitCont <-
  82     function(phy, model = "BM", sigma = 0.1, alpha = 1,
  83              theta = 0, ancestor = FALSE, root.value = 0)
  84 {
  85     if (is.null(phy$edge.length))
  86         stop("tree has no branch length")
  87     if (any(phy$edge.length < 0))
  88         stop("at least one branch length negative")
  89
  90     phy <- reorder(phy, "pruningwise")
  91     n <- length(phy$tip.label)
  92     N <- dim(phy$edge)[1]
  93     ROOT <- n + 1L
  94     x <- numeric(n + phy$Nnode)
  95     x[ROOT] <- root.value
  96
  97     anc <- phy$edge[, 1]
  98     des <- phy$edge[, 2]
  99     el <- phy$edge.length
 100
 101     if (is.function(model)) {
 102         environment(model) <- environment()
 103         for (i in N:1) x[des[i]] <- model(x[anc[i]], el[i])
 104     } else {
 105         model <- pmatch(model, c("BM", "OU"))
 106         if (length(sigma) == 1) sigma <- rep(sigma, N)
 107         else if (length(sigma) != N)
 108             stop("'sigma' must have one or Nedge(phy) elements")
 109         if (model == 2) { # "OU"
 110             if (length(alpha) == 1) alpha <- rep(alpha, N)
 111             else if (length(alpha) != N)
 112                 stop("'alpha' must have one or Nedge(phy) elements")
 113             if (length(theta) == 1) theta <- rep(theta, N)
 114             else if (length(theta) != N)
 115                 stop("'theta' must have one or Nedge(phy) elements")
 116         }
 117         .C("rTraitCont", as.integer(model), as.integer(N), as.integer(anc - 1L), as.integer(des - 1L), el, sigma, alpha, theta, x, DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")
 118     }
 119
 120     if (ancestor) {
 121         if (is.null(phy$node.label)) phy <- makeNodeLabel(phy)
 122         names(x) <- c(phy$tip.label, phy$node.label)
 123     } else {
 124         x <- x[1:n]
 125         names(x) <- phy$tip.label
 126     }
 127     x
 128 }