1 ## rtree.R (2012-02-09)
5 ## Copyright 2004-2012 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 rtree <- function(n, rooted = TRUE, tip.label = NULL, br = runif, ...)
12 foo <- function(n, pos) {
13 n1 <- sample.int(n - 1, 1, FALSE, NULL)
16 edge[c(pos, po2), 1] <<- nod
22 edge[c(pos + 1, pos + 2), 1] <<- edge[pos, 2] <<- nod
29 edge[c(po2 + 1, po2 + 2), 1] <<- edge[po2, 2] <<- nod
34 if (n < 2) stop("a tree must have at least 2 tips.")
35 nbr <- 2 * n - 3 + rooted
36 edge <- matrix(NA, nbr, 2)
40 if (rooted) edge[] <- c(3L, 3L, 1L, 2L)
41 else stop("an unrooted tree must have at least 3 tips.")
44 if (rooted) c(4L, 5L, 5L, 4L, 5L, 1:3)
45 else c(4L, 4L, 4L, 1:3)
46 } else if (n == 4 && !rooted) {
47 edge[] <- c(5L, 6L, 6L, 5L, 5L, 6L, 1:4)
52 ## The following is slightly more efficient than affecting the
53 ## tip numbers in foo(): the gain is 0.006 s for n = 1000.
54 i <- which(is.na(edge[, 2]))
57 n1 <- sample.int(n - 2, 1, FALSE, NULL)
61 n2 <- sample.int(n - n1 - 1, 1, FALSE, NULL)
65 po3 <- 2*(n1 + n2) - 1
66 edge[c(1, po2, po3), 1] <- nod
72 edge[2:3, 1] <- edge[1, 2] <- nod
79 edge[c(po2 + 1, po2 + 2), 1] <- edge[po2, 2] <- nod
86 edge[c(po3 + 1, po3 + 2), 1] <- edge[po3, 2] <- nod
89 i <- which(is.na(edge[, 2]))
93 phy <- list(edge = edge)
95 if (is.null(tip.label)) paste("t", sample(n), sep = "")
96 else sample(tip.label)
99 if (is.function(br)) br(nbr, ...) else rep(br, length.out = nbr)
101 phy$Nnode <- n - 2L + rooted
102 class(phy) <- "phylo"
106 rcoal <- function(n, tip.label = NULL, br = "coalescent", ...)
110 edge <- matrix(NA, nbr, 2)
111 ## coalescence times by default:
112 x <- if (is.character(br)) 2*rexp(n - 1)/(as.double(n:2) * as.double((n - 1):1))
113 else if (is.numeric(br)) rep(br, length.out = n - 1) else br(n - 1, ...)
115 edge[] <- c(3L, 3L, 1:2)
116 edge.length <- rep(x, 2)
118 edge[] <- c(4L, 5L, 5L, 4L, 5L, 1:3)
119 edge.length <- c(x[c(2, 1, 1)], sum(x))
121 edge.length <- numeric(nbr)
122 h <- numeric(2*n - 1) # initialized with 0's
123 node.height <- cumsum(x)
125 nextnode <- 2L*n - 1L
126 for (i in 1:(n - 1)) {
127 y <- sample(pool, size = 2)
128 ind <- (i - 1)*2 + 1:2
130 edge[ind, 1] <- nextnode
131 edge.length[ind] <- node.height[i] - h[y]
132 h[nextnode] <- node.height[i]
133 pool <- c(pool[! pool %in% y], nextnode)
134 nextnode <- nextnode - 1L
137 phy <- list(edge = edge, edge.length = edge.length)
138 if (is.null(tip.label))
139 tip.label <- paste("t", 1:n, sep = "")
140 phy$tip.label <- sample(tip.label)
142 class(phy) <- "phylo"
144 ## to avoid crossings when converting with as.hclust:
145 phy$edge[phy$edge[, 2] <= n, 2] <- 1:n
149 rmtree <- function(N, n, rooted = TRUE, tip.label = NULL, br = runif, ...)
151 a <- replicate(N, rtree(n, rooted = rooted, tip.label = tip.label,
152 br = br, ...), simplify = FALSE)
153 class(a) <- "multiPhylo"
157 stree <- function(n, type = "star", tip.label = NULL)
159 type <- match.arg(type, c("star", "balanced", "left", "right"))
161 if (type == "star") {
168 edge <- matrix(0L, N, 2)
170 switch(type, "star" = {
175 stop("'n' is not a power of 2: cannot make a balanced tree")
176 foo <- function(node, size) {
178 edge[c(i, i + 1L), 1L] <<- node
179 edge[c(i, i + 1L), 2L] <<- c(nexttip, nexttip + 1L)
180 nexttip <<- nexttip + 2L
183 for (k in 1:2) { # do the 2 subclades
184 edge[i, ] <<- c(node, nextnode)
185 nextnode <<- nextnode + 1L
187 foo(nextnode - 1L, size/2)
196 edge[c(seq.int(from = 1, to = N - 1, by = 2), N), 2L] <- 1:n
197 nodes <- (n + 1L):(n + m)
198 edge[seq.int(from = 2, to = N - 1, by = 2), 2L] <- nodes[-1]
199 edge[, 1L] <- rep(nodes, each = 2)
201 nodes <- (n + 1L):(n + m)
202 edge[, 1L] <- c(nodes, rev(nodes))
203 edge[, 2L] <- c(nodes[-1], 1:n)
206 if (is.null(tip.label))
207 tip.label <- paste("t", 1:n, sep = "")
208 phy <- list(edge = edge, tip.label = tip.label, Nnode = m)
209 class(phy) <- "phylo"
210 attr(phy, "order") <- "cladewise"