1 ## rtree.R (2008-05-07)
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.")
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 <- .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 <- n - 2L + rooted
103 rcoal <- function(n, tip.label = NULL, br = "coalescent", ...)
107 edge <- matrix(NA, nbr, 2)
108 ## coalescence times by default:
109 x <- if (is.character(br)) 2*rexp(n - 1)/(n:2 * (n - 1):1)
112 edge[] <- c(3L, 3L, 1:2)
113 edge.length <- rep(x, 2)
115 edge[] <- c(4L, 5L, 5L, 4L, 5L, 1:3)
116 edge.length <- c(x[2], x[1], x[1], sum(x))
118 edge.length <- numeric(nbr)
119 h <- numeric(2*n - 1) # initialized with 0's
120 node.height <- cumsum(x)
122 nextnode <- 2L*n - 1L
123 for (i in 1:(n - 1)) {
124 y <- sample(pool, size = 2)
125 ind <- (i - 1)*2 + 1:2
127 edge[ind, 1] <- nextnode
128 edge.length[ind] <- node.height[i] - h[y]
129 h[nextnode] <- node.height[i]
130 pool <- c(pool[! pool %in% y], nextnode)
131 nextnode <- nextnode - 1L
134 phy <- list(edge = edge, edge.length = edge.length)
136 if (is.null(tip.label)) paste("t", 1:n, sep = "")
139 class(phy) <- "phylo"
141 ## to avoid crossings when converting with as.hclust:
142 read.tree(text = write.tree(phy))
145 rmtree <- function(N, n, rooted = TRUE, tip.label = NULL, br = runif, ...)
147 a <- replicate(N, rtree(n, rooted = rooted, tip.label = tip.label,
148 br = br, ...), simplify = FALSE)
149 class(a) <- "multiPhylo"