R/rtree.R

   1 ## rtree.R (2013-04-02)
   2
   3 ##   Generates Trees
   4
   5 ## Copyright 2004-2013 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 <- sample.int(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 <- sample.int(n - 2L, 1L)
  58             if (n1 == n - 2L) {
  59                 n2 <- n3 <- 1L
  60             } else {
  61                 n2 <- sample.int(n - n1 - 1L, 1L)
  62                 n3 <- n - n1 - n2
  63             }
  64             po2 <- 2L * n1
  65             po3 <- 2L * (n1 + n2) - 1L
  66             edge[c(1, po2, po3), 1L] <- nod
  67             nod <- nod + 1L
  68             if (n1 > 2L) {
  69                 edge[1L, 2L] <- nod
  70                 foo(n1, 2L)
  71             } else if (n1 == 2L) {
  72                 edge[2:3, 1L] <- edge[1L, 2L] <- nod
  73                 nod <- nod + 1L
  74             }
  75             if (n2 > 2L) {
  76                 edge[po2, 2L] <- nod
  77                 foo(n2, po2 + 1L)
  78             } else if (n2 == 2L) {
  79                 edge[c(po2 + 1L, po2 + 2), 1L] <- edge[po2, 2L] <- nod
  80                 nod <- nod + 1L
  81             }
  82             if (n3 > 2) {
  83                 edge[po3, 2L] <- nod
  84                 foo(n3, po3 + 1L)
  85             } else if (n3 == 2L) {
  86                 edge[c(po3 + 1L, po3 + 2), 1L] <- edge[po3, 2L] <- nod
  87                 ## nod <- nod + 1L
  88             }
  89             i <- which(is.na(edge[, 2L]))
  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     attr(phy, "order") <- "cladewise"
 104     phy
 105 }
 106
 107 rcoal <- function(n, tip.label = NULL, br = "coalescent", ...)
 108 {
 109     n <- as.integer(n)
 110     nbr <- 2*n - 2
 111     edge <- matrix(NA, nbr, 2)
 112     ## coalescence times by default:
 113     x <- if (is.character(br)) 2*rexp(n - 1)/(as.double(n:2) * as.double((n - 1):1))
 114     else if (is.numeric(br)) rep(br, length.out = n - 1) else br(n - 1, ...)
 115     if (n == 2) {
 116         edge[] <- c(3L, 3L, 1:2)
 117         edge.length <- rep(x, 2)
 118     } else if (n == 3) {
 119         edge[] <- c(4L, 5L, 5L, 4L, 5L, 1:3)
 120         edge.length <- c(x[c(2, 1, 1)], sum(x))
 121     } else {
 122         edge.length <- numeric(nbr)
 123         h <- numeric(2*n - 1) # initialized with 0's
 124         node.height <- cumsum(x)
 125         pool <- 1:n
 126         nextnode <- 2L*n - 1L
 127         for (i in 1:(n - 1)) {
 128             y <- sample(pool, size = 2)
 129             ind <- (i - 1)*2 + 1:2
 130             edge[ind, 2] <- y
 131             edge[ind, 1] <- nextnode
 132             edge.length[ind] <- node.height[i] - h[y]
 133             h[nextnode] <- node.height[i]
 134             pool <- c(pool[! pool %in% y], nextnode)
 135             nextnode <- nextnode - 1L
 136         }
 137     }
 138     phy <- list(edge = edge, edge.length = edge.length)
 139     if (is.null(tip.label))
 140         tip.label <- paste("t", 1:n, sep = "")
 141     phy$tip.label <- sample(tip.label)
 142     phy$Nnode <- n - 1L
 143     class(phy) <- "phylo"
 144     phy <- reorder(phy)
 145     ## to avoid crossings when converting with as.hclust:
 146     phy$edge[phy$edge[, 2] <= n, 2] <- 1:n
 147     phy
 148 }
 149
 150 rmtree <- function(N, n, rooted = TRUE, tip.label = NULL, br = runif, ...)
 151 {
 152     a <- replicate(N, rtree(n, rooted = rooted, tip.label =  tip.label,
 153                             br = br, ...), simplify = FALSE)
 154     class(a) <- "multiPhylo"
 155     a
 156 }
 157
 158 stree <- function(n, type = "star", tip.label = NULL)
 159 {
 160     type <- match.arg(type, c("star", "balanced", "left", "right"))
 161     n <- as.integer(n)
 162     if (type == "star") {
 163         N <- n
 164         m <- 1L
 165     } else {
 166         m <- n - 1L
 167         N <- n + m - 1L
 168     }
 169     edge <- matrix(0L, N, 2)
 170
 171     switch(type, "star" = {
 172         edge[, 1] <- n + 1L
 173         edge[, 2] <- 1:n
 174     }, "balanced" = {
 175         if (log2(n) %% 1)
 176             stop("'n' is not a power of 2: cannot make a balanced tree")
 177         foo <- function(node, size) {
 178             if (size == 2) {
 179                 edge[c(i, i + 1L), 1L] <<- node
 180                 edge[c(i, i + 1L), 2L] <<- c(nexttip, nexttip + 1L)
 181                 nexttip <<- nexttip + 2L
 182                 i <<- i + 2L
 183             } else {
 184                 for (k in 1:2) { # do the 2 subclades
 185                     edge[i, ] <<- c(node, nextnode)
 186                     nextnode <<- nextnode + 1L
 187                     i <<- i + 1L
 188                     foo(nextnode - 1L, size/2)
 189                 }
 190             }
 191         }
 192         i <- 1L
 193         nexttip <- 1L
 194         nextnode <- n + 2L
 195         foo(n + 1L, n)
 196     }, "left" = {
 197         edge[c(seq.int(from = 1, to = N - 1, by = 2), N), 2L] <- 1:n
 198         nodes <- (n + 1L):(n + m)
 199         edge[seq.int(from = 2, to = N - 1, by = 2), 2L] <- nodes[-1]
 200         edge[, 1L] <- rep(nodes, each = 2)
 201     }, "right" = {
 202         nodes <- (n + 1L):(n + m)
 203         edge[, 1L] <- c(nodes, rev(nodes))
 204         edge[, 2L] <- c(nodes[-1], 1:n)
 205     })
 206
 207     if (is.null(tip.label))
 208         tip.label <- paste("t", 1:n, sep = "")
 209     phy <- list(edge = edge, tip.label = tip.label, Nnode = m)
 210     class(phy) <- "phylo"
 211     attr(phy, "order") <- "cladewise"
 212     phy
 213 }