R/rtree.R

   1 ## rtree.R (2007-12-22)
   2
   3 ##   Generates Random Trees
   4
   5 ## Copyright 2004-2007 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 + 1
  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 + 1
  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 + 1
  31         }
  32     }
  33
  34     if (n < 2) stop("a tree must have at least 2 tips.")
  35     nbr <- 2 * n - 2
  36     if (!rooted) nbr <- nbr - 1
  37     edge <- matrix(NA, nbr, 2)
  38
  39     if (n == 2) {
  40         if (rooted) edge[] <- c(3, 3, 1, 2)
  41         else stop("an unrooted tree must have at least 3 tips.")
  42     } else if (n == 3) {
  43         edge[] <-
  44           if (rooted) c(4, 5, 5, 4, 5, 1:3)
  45           else c(4, 4, 4, 1:3)
  46     } else if (n == 4 && !rooted) {
  47         edge[] <- c(5, 6, 6, 5, 5, 6, 1:4)
  48     } else {
  49         nod <- n + 1
  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 + 1
  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 + 1
  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 + 1
  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 <- if (rooted) n - 1 else n - 2
  99     class(phy) <- "phylo"
 100     phy
 101 }
 102
 103 rcoal <- function(n, tip.label = NULL, br = rexp, ...)
 104 {
 105     nbr <- 2*n - 2
 106     edge <- matrix(NA, nbr, 2)
 107     x <- br(n - 1, ...) # coalescence times
 108     if (n == 2) {
 109         edge[] <- c(3, 3, 1:2)
 110         edge.length <- rep(x, 2)
 111     } else if (n == 3) {
 112         edge[] <- c(4, 5, 5, 4, 5, 1:3)
 113         edge.length <- c(x[2], x[1], x[1], sum(x))
 114     } else {
 115         edge.length <- numeric(nbr)
 116         h <- numeric(2*n - 1) # initialized with 0's
 117         node.height <- cumsum(x)
 118         pool <- 1:n
 119         nextnode <- 2*n - 1
 120         for (i in 1:(n - 1)) {
 121             y <- sample(pool, size = 2)
 122             ind <- (i - 1)*2 + 1:2
 123             edge[ind, 2] <- y
 124             edge[ind, 1] <- nextnode
 125             edge.length[ind] <- node.height[i] - h[y]
 126             h[nextnode] <- node.height[i]
 127             pool <- c(pool[! pool %in% y], nextnode)
 128             nextnode <- nextnode - 1
 129         }
 130     }
 131     phy <- list(edge = edge, edge.length = edge.length)
 132     phy$tip.label <-
 133       if (is.null(tip.label)) paste("t", 1:n, sep = "")
 134       else tip.label
 135     phy$Nnode <- n - 1
 136     class(phy) <- "phylo"
 137     ##reorder(phy)
 138     ## to avoid crossings when converting with as.hclust:
 139     read.tree(text = write.tree(phy))
 140 }