1 ## rtree.R (2008-01-13)
3 ## Generates Random Trees
5 ## Copyright 2004-2008 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 <- .Internal(sample(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.")
36 if (!rooted) nbr <- nbr - 1
37 edge <- matrix(NA, nbr, 2)
40 if (rooted) edge[] <- c(3, 3, 1, 2)
41 else stop("an unrooted tree must have at least 3 tips.")
44 if (rooted) c(4, 5, 5, 4, 5, 1:3)
46 } else if (n == 4 && !rooted) {
47 edge[] <- c(5, 6, 6, 5, 5, 6, 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 <- .Internal(sample(n - 2, 1, FALSE, NULL))
61 n2 <- .Internal(sample(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)
97 if (is.function(br)) phy$edge.length <- br(nbr, ...)
98 phy$Nnode <- if (rooted) n - 1 else n - 2
103 rcoal <- function(n, tip.label = NULL, br = "coalescent", ...)
106 edge <- matrix(NA, nbr, 2)
107 ## coalescence times by default:
108 x <- if (is.character(br)) 2*rexp(n - 1)/(n:2 * (n - 1):1)
111 edge[] <- c(3, 3, 1:2)
112 edge.length <- rep(x, 2)
114 edge[] <- c(4, 5, 5, 4, 5, 1:3)
115 edge.length <- c(x[2], x[1], x[1], sum(x))
117 edge.length <- numeric(nbr)
118 h <- numeric(2*n - 1) # initialized with 0's
119 node.height <- cumsum(x)
122 for (i in 1:(n - 1)) {
123 y <- sample(pool, size = 2)
124 ind <- (i - 1)*2 + 1:2
126 edge[ind, 1] <- nextnode
127 edge.length[ind] <- node.height[i] - h[y]
128 h[nextnode] <- node.height[i]
129 pool <- c(pool[! pool %in% y], nextnode)
130 nextnode <- nextnode - 1
133 phy <- list(edge = edge, edge.length = edge.length)
135 if (is.null(tip.label)) paste("t", 1:n, sep = "")
138 class(phy) <- "phylo"
140 ## to avoid crossings when converting with as.hclust:
141 read.tree(text = write.tree(phy))
144 rmtree <- function(N, n, rooted = TRUE, tip.label = NULL, br = runif, ...)
146 a <- replicate(N, rtree(n, rooted = rooted, tip.label = tip.label,
147 br = br, ...), simplify = FALSE)
148 class(a) <- "multiPhylo"