1 ## rTrait.R (2010-01-11)
5 ## Copyright 2010 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
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)
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")
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)
26 if (length(rate) != k*(k - 1))
27 stop("`rate' must have k(k - 1) elements")
29 Q[col(Q) != row(Q)] <- rate
31 if (length(rate) != k*(k - 1)/2)
32 stop("`rate' must have k(k - 1)/2 elements")
34 sel <- col(Q) < row(Q)
40 if (is.matrix(model)) {
42 if (ncol(Q) != nrow(Q))
43 stop("the matrix given as `model' must be square")
46 phy <- reorder(phy, "pruningwise")
47 n <- length(phy$tip.label)
50 x <- integer(n + phy$Nnode)
51 x[ROOT] <- as.integer(root.value)
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])
61 diag(Q) <- -rowSums(Q)
63 p <- matexpo(Q * freq * el[i])[x[anc[i]], ]
64 x[des[i]] <- .Internal(sample(k, size = 1, FALSE, prob = p))
69 if (is.null(phy$node.label)) phy <- makeNodeLabel(phy)
70 names(x) <- c(phy$tip.label, phy$node.label)
73 names(x) <- phy$tip.label
81 function(phy, model = "BM", sigma = 0.1, alpha = 1,
82 theta = 0, ancestor = FALSE, root.value = 0)
84 if (is.null(phy$edge.length))
85 stop("tree has no branch length")
86 if (any(phy$edge.length < 0))
87 stop("at least one branch length negative")
89 phy <- reorder(phy, "pruningwise")
90 n <- length(phy$tip.label)
93 x <- numeric(n + phy$Nnode)
100 if (is.function(model)) {
101 environment(model) <- environment()
102 for (i in N:1) x[des[i]] <- model(x[anc[i]], el[i])
104 model <- pmatch(model, c("BM", "OU"))
105 if (length(sigma) == 1) sigma <- rep(sigma, N)
106 else if (length(sigma) != N)
107 stop("'sigma' must have one or Nedge(phy) elements")
108 if (model == 2) { # "OU"
109 if (length(alpha) == 1) alpha <- rep(alpha, N)
110 else if (length(alpha) != N)
111 stop("'alpha' must have one or Nedge(phy) elements")
112 if (length(theta) == 1) theta <- rep(theta, N)
113 else if (length(theta) != N)
114 stop("'theta' must have one or Nedge(phy) elements")
116 .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")
120 if (is.null(phy$node.label)) phy <- makeNodeLabel(phy)
121 names(x) <- c(phy$tip.label, phy$node.label)
124 names(x) <- phy$tip.label