final commit for ape 3.0-8

[ape.git] / R / rtree.R

diff --git a/R/rtree.R b/R/rtree.R
index f331ca3101f0472575ef35b6e92f3abf84021764..b36f9d307a40c655cbaf85b2395fada04839d64b 100644 (file)
--- a/R/rtree.R
+++ b/R/rtree.R
@@ -1,8 +1,8 @@
-## rtree.R (2000-05-13)
+## rtree.R (2013-04-02)
 
-##   Generates Random Trees
+##   Generates Trees
 
-## Copyright 2004-2009 Emmanuel Paradis
+## Copyright 2004-2013 Emmanuel Paradis
 
 ## This file is part of the R-package `ape'.
 ## See the file ../COPYING for licensing issues.
@@ -10,7 +10,7 @@
 rtree <- function(n, rooted = TRUE, tip.label = NULL, br = runif, ...)
 {
     foo <- function(n, pos) {
-        n1 <- .Internal(sample(n - 1, 1, FALSE, NULL))
+        n1 <- sample.int(n - 1, 1, FALSE, NULL)
         n2 <- n - n1
         po2 <- pos + 2*n1 - 1
         edge[c(pos, po2), 1] <<- nod
@@ -54,39 +54,39 @@ rtree <- function(n, rooted = TRUE, tip.label = NULL, br = runif, ...)
             i <- which(is.na(edge[, 2]))
             edge[i, 2] <- 1:n
         } else { # n > 4
-            n1 <- .Internal(sample(n - 2, 1, FALSE, NULL))
-            if (n1 == n - 2) {
-                n2 <- n3 <- 1
+            n1 <- sample.int(n - 2L, 1L)
+            if (n1 == n - 2L) {
+                n2 <- n3 <- 1L
             } else {
-                n2 <- .Internal(sample(n - n1 - 1, 1, FALSE, NULL))
+                n2 <- sample.int(n - n1 - 1L, 1L)
                 n3 <- n - n1 - n2
             }
-            po2 <- 2*n1
-            po3 <- 2*(n1 + n2) - 1
-            edge[c(1, po2, po3), 1] <- nod
-            nod <- nod + 1
-            if (n1 > 2) {
-                edge[1, 2] <- nod
-                foo(n1, 2)
-            } else if (n1 == 2) {
-                edge[2:3, 1] <- edge[1, 2] <- nod
+            po2 <- 2L * n1
+            po3 <- 2L * (n1 + n2) - 1L
+            edge[c(1, po2, po3), 1L] <- nod
+            nod <- nod + 1L
+            if (n1 > 2L) {
+                edge[1L, 2L] <- nod
+                foo(n1, 2L)
+            } else if (n1 == 2L) {
+                edge[2:3, 1L] <- edge[1L, 2L] <- nod
                 nod <- nod + 1L
             }
-            if (n2 > 2) {
-                edge[po2, 2] <- nod
-                foo(n2, po2 + 1)
-            } else if (n2 == 2) {
-                edge[c(po2 + 1, po2 + 2), 1] <- edge[po2, 2] <- nod
+            if (n2 > 2L) {
+                edge[po2, 2L] <- nod
+                foo(n2, po2 + 1L)
+            } else if (n2 == 2L) {
+                edge[c(po2 + 1L, po2 + 2), 1L] <- edge[po2, 2L] <- nod
                 nod <- nod + 1L
             }
             if (n3 > 2) {
-                edge[po3, 2] <- nod
-                foo(n3, po3 + 1)
-            } else if (n3 == 2) {
-                edge[c(po3 + 1, po3 + 2), 1] <- edge[po3, 2] <- nod
+                edge[po3, 2L] <- nod
+                foo(n3, po3 + 1L)
+            } else if (n3 == 2L) {
+                edge[c(po3 + 1L, po3 + 2), 1L] <- edge[po3, 2L] <- nod
                 ## nod <- nod + 1L
             }
-            i <- which(is.na(edge[, 2]))
+            i <- which(is.na(edge[, 2L]))
             edge[i, 2] <- 1:n
         }
     }
@@ -94,9 +94,13 @@ rtree <- function(n, rooted = TRUE, tip.label = NULL, br = runif, ...)
     phy$tip.label <-
       if (is.null(tip.label)) paste("t", sample(n), sep = "")
       else sample(tip.label)
-    if (is.function(br)) phy$edge.length <- br(nbr, ...)
+    if (!is.null(br)) {
+        phy$edge.length <-
+            if (is.function(br)) br(nbr, ...) else rep(br, length.out = nbr)
+    }
     phy$Nnode <- n - 2L + rooted
     class(phy) <- "phylo"
+    attr(phy, "order") <- "cladewise"
     phy
 }
 
@@ -106,8 +110,8 @@ rcoal <- function(n, tip.label = NULL, br = "coalescent", ...)
     nbr <- 2*n - 2
     edge <- matrix(NA, nbr, 2)
     ## coalescence times by default:
-    x <- if (is.character(br)) 2*rexp(n - 1)/(n:2 * (n - 1):1)
-    else br(n - 1, ...)
+    x <- if (is.character(br)) 2*rexp(n - 1)/(as.double(n:2) * as.double((n - 1):1))
+    else if (is.numeric(br)) rep(br, length.out = n - 1) else br(n - 1, ...)
     if (n == 2) {
         edge[] <- c(3L, 3L, 1:2)
         edge.length <- rep(x, 2)
@@ -150,3 +154,60 @@ rmtree <- function(N, n, rooted = TRUE, tip.label = NULL, br = runif, ...)
     class(a) <- "multiPhylo"
     a
 }
+
+stree <- function(n, type = "star", tip.label = NULL)
+{
+    type <- match.arg(type, c("star", "balanced", "left", "right"))
+    n <- as.integer(n)
+    if (type == "star") {
+        N <- n
+        m <- 1L
+    } else {
+        m <- n - 1L
+        N <- n + m - 1L
+    }
+    edge <- matrix(0L, N, 2)
+
+    switch(type, "star" = {
+        edge[, 1] <- n + 1L
+        edge[, 2] <- 1:n
+    }, "balanced" = {
+        if (log2(n) %% 1)
+            stop("'n' is not a power of 2: cannot make a balanced tree")
+        foo <- function(node, size) {
+            if (size == 2) {
+                edge[c(i, i + 1L), 1L] <<- node
+                edge[c(i, i + 1L), 2L] <<- c(nexttip, nexttip + 1L)
+                nexttip <<- nexttip + 2L
+                i <<- i + 2L
+            } else {
+                for (k in 1:2) { # do the 2 subclades
+                    edge[i, ] <<- c(node, nextnode)
+                    nextnode <<- nextnode + 1L
+                    i <<- i + 1L
+                    foo(nextnode - 1L, size/2)
+                }
+            }
+        }
+        i <- 1L
+        nexttip <- 1L
+        nextnode <- n + 2L
+        foo(n + 1L, n)
+    }, "left" = {
+        edge[c(seq.int(from = 1, to = N - 1, by = 2), N), 2L] <- 1:n
+        nodes <- (n + 1L):(n + m)
+        edge[seq.int(from = 2, to = N - 1, by = 2), 2L] <- nodes[-1]
+        edge[, 1L] <- rep(nodes, each = 2)
+    }, "right" = {
+        nodes <- (n + 1L):(n + m)
+        edge[, 1L] <- c(nodes, rev(nodes))
+        edge[, 2L] <- c(nodes[-1], 1:n)
+    })
+
+    if (is.null(tip.label))
+        tip.label <- paste("t", 1:n, sep = "")
+    phy <- list(edge = edge, tip.label = tip.label, Nnode = m)
+    class(phy) <- "phylo"
+    attr(phy, "order") <- "cladewise"
+    phy
+}