1 ## rTrait.R (2010-02-03)
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)
62 freq <- rep(freq, each = k)
64 p <- matexpo(Q * freq * el[i])[x[anc[i]], ]
65 x[des[i]] <- .Internal(sample(k, size = 1, FALSE, prob = p))
70 if (is.null(phy$node.label)) phy <- makeNodeLabel(phy)
71 names(x) <- c(phy$tip.label, phy$node.label)
74 names(x) <- phy$tip.label
82 function(phy, model = "BM", sigma = 0.1, alpha = 1,
83 theta = 0, ancestor = FALSE, root.value = 0)
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")
90 phy <- reorder(phy, "pruningwise")
91 n <- length(phy$tip.label)
94 x <- numeric(n + phy$Nnode)
101 if (is.function(model)) {
102 environment(model) <- environment()
103 for (i in N:1) x[des[i]] <- model(x[anc[i]], el[i])
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")
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")
121 if (is.null(phy$node.label)) phy <- makeNodeLabel(phy)
122 names(x) <- c(phy$tip.label, phy$node.label)
125 names(x) <- phy$tip.label