1 ## rtree.R (2009-11-03)
3 ## Generates Random Trees
5 ## Copyright 2004-2009 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)
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"