1 ## rtree.R (2013-04-02)
5 ## Copyright 2004-2013 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 - 2L, 1L)
61 n2 <- sample.int(n - n1 - 1L, 1L)
65 po3 <- 2L * (n1 + n2) - 1L
66 edge[c(1, po2, po3), 1L] <- nod
71 } else if (n1 == 2L) {
72 edge[2:3, 1L] <- edge[1L, 2L] <- nod
78 } else if (n2 == 2L) {
79 edge[c(po2 + 1L, po2 + 2), 1L] <- edge[po2, 2L] <- nod
85 } else if (n3 == 2L) {
86 edge[c(po3 + 1L, po3 + 2), 1L] <- edge[po3, 2L] <- nod
89 i <- which(is.na(edge[, 2L]))
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"
103 attr(phy, "order") <- "cladewise"
107 rcoal <- function(n, tip.label = NULL, br = "coalescent", ...)
111 edge <- matrix(NA, nbr, 2)
112 ## coalescence times by default:
113 x <- if (is.character(br)) 2*rexp(n - 1)/(as.double(n:2) * as.double((n - 1):1))
114 else if (is.numeric(br)) rep(br, length.out = n - 1) else br(n - 1, ...)
116 edge[] <- c(3L, 3L, 1:2)
117 edge.length <- rep(x, 2)
119 edge[] <- c(4L, 5L, 5L, 4L, 5L, 1:3)
120 edge.length <- c(x[c(2, 1, 1)], sum(x))
122 edge.length <- numeric(nbr)
123 h <- numeric(2*n - 1) # initialized with 0's
124 node.height <- cumsum(x)
126 nextnode <- 2L*n - 1L
127 for (i in 1:(n - 1)) {
128 y <- sample(pool, size = 2)
129 ind <- (i - 1)*2 + 1:2
131 edge[ind, 1] <- nextnode
132 edge.length[ind] <- node.height[i] - h[y]
133 h[nextnode] <- node.height[i]
134 pool <- c(pool[! pool %in% y], nextnode)
135 nextnode <- nextnode - 1L
138 phy <- list(edge = edge, edge.length = edge.length)
139 if (is.null(tip.label))
140 tip.label <- paste("t", 1:n, sep = "")
141 phy$tip.label <- sample(tip.label)
143 class(phy) <- "phylo"
145 ## to avoid crossings when converting with as.hclust:
146 phy$edge[phy$edge[, 2] <= n, 2] <- 1:n
150 rmtree <- function(N, n, rooted = TRUE, tip.label = NULL, br = runif, ...)
152 a <- replicate(N, rtree(n, rooted = rooted, tip.label = tip.label,
153 br = br, ...), simplify = FALSE)
154 class(a) <- "multiPhylo"
158 stree <- function(n, type = "star", tip.label = NULL)
160 type <- match.arg(type, c("star", "balanced", "left", "right"))
162 if (type == "star") {
169 edge <- matrix(0L, N, 2)
171 switch(type, "star" = {
176 stop("'n' is not a power of 2: cannot make a balanced tree")
177 foo <- function(node, size) {
179 edge[c(i, i + 1L), 1L] <<- node
180 edge[c(i, i + 1L), 2L] <<- c(nexttip, nexttip + 1L)
181 nexttip <<- nexttip + 2L
184 for (k in 1:2) { # do the 2 subclades
185 edge[i, ] <<- c(node, nextnode)
186 nextnode <<- nextnode + 1L
188 foo(nextnode - 1L, size/2)
197 edge[c(seq.int(from = 1, to = N - 1, by = 2), N), 2L] <- 1:n
198 nodes <- (n + 1L):(n + m)
199 edge[seq.int(from = 2, to = N - 1, by = 2), 2L] <- nodes[-1]
200 edge[, 1L] <- rep(nodes, each = 2)
202 nodes <- (n + 1L):(n + m)
203 edge[, 1L] <- c(nodes, rev(nodes))
204 edge[, 2L] <- c(nodes[-1], 1:n)
207 if (is.null(tip.label))
208 tip.label <- paste("t", 1:n, sep = "")
209 phy <- list(edge = edge, tip.label = tip.label, Nnode = m)
210 class(phy) <- "phylo"
211 attr(phy, "order") <- "cladewise"