R/rtree.R

   1 ## rtree.R (2000-05-13)
   2
   3 ##   Generates Random Trees
   4
   5 ## Copyright 2004-2009 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 rtree <- function(n, rooted = TRUE, tip.label = NULL, br = runif, ...)
  11 {
  12     foo <- function(n, pos) {
  13         n1 <- .Internal(sample(n - 1, 1, FALSE, NULL))
  14         n2 <- n - n1
  15         po2 <- pos + 2*n1 - 1
  16         edge[c(pos, po2), 1] <<- nod
  17         nod <<- nod + 1L
  18         if (n1 > 2) {
  19             edge[pos, 2] <<- nod
  20             foo(n1, pos + 1)
  21         } else if (n1 == 2) {
  22             edge[c(pos + 1, pos + 2), 1] <<- edge[pos, 2] <<- nod
  23             nod <<- nod + 1L
  24         }
  25         if (n2 > 2) {
  26             edge[po2, 2] <<- nod
  27             foo(n2, po2 + 1)
  28         } else if (n2 == 2) {
  29             edge[c(po2 + 1, po2 + 2), 1] <<- edge[po2, 2] <<- nod
  30             nod <<- nod + 1L
  31         }
  32     }
  33
  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)
  37
  38     n <- as.integer(n)
  39     if (n == 2) {
  40         if (rooted) edge[] <- c(3L, 3L, 1L, 2L)
  41         else stop("an unrooted tree must have at least 3 tips.")
  42     } else if (n == 3) {
  43         edge[] <-
  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)
  48     } else {
  49         nod <- n + 1L
  50         if (rooted) { # n > 3
  51             foo(n, 1)
  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]))
  55             edge[i, 2] <- 1:n
  56         } else { # n > 4
  57             n1 <- .Internal(sample(n - 2, 1, FALSE, NULL))
  58             if (n1 == n - 2) {
  59                 n2 <- n3 <- 1
  60             } else {
  61                 n2 <- .Internal(sample(n - n1 - 1, 1, FALSE, NULL))
  62                 n3 <- n - n1 - n2
  63             }
  64             po2 <- 2*n1
  65             po3 <- 2*(n1 + n2) - 1
  66             edge[c(1, po2, po3), 1] <- nod
  67             nod <- nod + 1
  68             if (n1 > 2) {
  69                 edge[1, 2] <- nod
  70                 foo(n1, 2)
  71             } else if (n1 == 2) {
  72                 edge[2:3, 1] <- edge[1, 2] <- nod
  73                 nod <- nod + 1L
  74             }
  75             if (n2 > 2) {
  76                 edge[po2, 2] <- nod
  77                 foo(n2, po2 + 1)
  78             } else if (n2 == 2) {
  79                 edge[c(po2 + 1, po2 + 2), 1] <- edge[po2, 2] <- nod
  80                 nod <- nod + 1L
  81             }
  82             if (n3 > 2) {
  83                 edge[po3, 2] <- nod
  84                 foo(n3, po3 + 1)
  85             } else if (n3 == 2) {
  86                 edge[c(po3 + 1, po3 + 2), 1] <- edge[po3, 2] <- nod
  87                 ## nod <- nod + 1L
  88             }
  89             i <- which(is.na(edge[, 2]))
  90             edge[i, 2] <- 1:n
  91         }
  92     }
  93     phy <- list(edge = edge)
  94     phy$tip.label <-
  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
  99     class(phy) <- "phylo"
 100     phy
 101 }
 102
 103 rcoal <- function(n, tip.label = NULL, br = "coalescent", ...)
 104 {
 105     n <- as.integer(n)
 106     nbr <- 2*n - 2
 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)
 110     else br(n - 1, ...)
 111     if (n == 2) {
 112         edge[] <- c(3L, 3L, 1:2)
 113         edge.length <- rep(x, 2)
 114     } else if (n == 3) {
 115         edge[] <- c(4L, 5L, 5L, 4L, 5L, 1:3)
 116         edge.length <- c(x[c(2, 1, 1)], sum(x))
 117     } else {
 118         edge.length <- numeric(nbr)
 119         h <- numeric(2*n - 1) # initialized with 0's
 120         node.height <- cumsum(x)
 121         pool <- 1:n
 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
 126             edge[ind, 2] <- y
 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
 132         }
 133     }
 134     phy <- list(edge = edge, edge.length = edge.length)
 135     if (is.null(tip.label))
 136         tip.label <- paste("t", 1:n, sep = "")
 137     phy$tip.label <- sample(tip.label)
 138     phy$Nnode <- n - 1L
 139     class(phy) <- "phylo"
 140     phy <- reorder(phy)
 141     ## to avoid crossings when converting with as.hclust:
 142     phy$edge[phy$edge[, 2] <= n, 2] <- 1:n
 143     phy
 144 }
 145
 146 rmtree <- function(N, n, rooted = TRUE, tip.label = NULL, br = runif, ...)
 147 {
 148     a <- replicate(N, rtree(n, rooted = rooted, tip.label =  tip.label,
 149                             br = br, ...), simplify = FALSE)
 150     class(a) <- "multiPhylo"
 151     a
 152 }