R/rtree.R

   1 ## rtree.R (2009-11-03)
   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.null(br)) {
  98         phy$edge.length <-
  99             if (is.function(br)) br(nbr, ...) else rep(br, length.out = nbr)
 100     }
 101     phy$Nnode <- n - 2L + rooted
 102     class(phy) <- "phylo"
 103     phy
 104 }
 105
 106 rcoal <- function(n, tip.label = NULL, br = "coalescent", ...)
 107 {
 108     n <- as.integer(n)
 109     nbr <- 2*n - 2
 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, ...)
 114     if (n == 2) {
 115         edge[] <- c(3L, 3L, 1:2)
 116         edge.length <- rep(x, 2)
 117     } else if (n == 3) {
 118         edge[] <- c(4L, 5L, 5L, 4L, 5L, 1:3)
 119         edge.length <- c(x[c(2, 1, 1)], sum(x))
 120     } else {
 121         edge.length <- numeric(nbr)
 122         h <- numeric(2*n - 1) # initialized with 0's
 123         node.height <- cumsum(x)
 124         pool <- 1:n
 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
 129             edge[ind, 2] <- y
 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
 135         }
 136     }
 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)
 141     phy$Nnode <- n - 1L
 142     class(phy) <- "phylo"
 143     phy <- reorder(phy)
 144     ## to avoid crossings when converting with as.hclust:
 145     phy$edge[phy$edge[, 2] <= n, 2] <- 1:n
 146     phy
 147 }
 148
 149 rmtree <- function(N, n, rooted = TRUE, tip.label = NULL, br = runif, ...)
 150 {
 151     a <- replicate(N, rtree(n, rooted = rooted, tip.label =  tip.label,
 152                             br = br, ...), simplify = FALSE)
 153     class(a) <- "multiPhylo"
 154     a
 155 }