updated phymltest() for PhyML 3.0 + changed rbind.DNAbin()

[ape.git] / R / rtree.R

diff --git a/R/rtree.R b/R/rtree.R
index 6d0172c1ef8b5029739dbd91eb7d027064f1e37f..9425dadf22f5df689990e1ffe6a3829e3ce99c6b 100644 (file)
--- a/R/rtree.R
+++ b/R/rtree.R
@@ -1,4 +1,4 @@
-## rtree.R (2008-01-08)
+## rtree.R (2008-05-07)
 
 ##   Generates Random Trees
 
@@ -14,39 +14,39 @@ rtree <- function(n, rooted = TRUE, tip.label = NULL, br = runif, ...)
         n2 <- n - n1
         po2 <- pos + 2*n1 - 1
         edge[c(pos, po2), 1] <<- nod
-        nod <<- nod + 1
+        nod <<- nod + 1L
         if (n1 > 2) {
             edge[pos, 2] <<- nod
             foo(n1, pos + 1)
         } else if (n1 == 2) {
             edge[c(pos + 1, pos + 2), 1] <<- edge[pos, 2] <<- nod
-            nod <<- nod + 1
+            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
-            nod <<- nod + 1
+            nod <<- nod + 1L
         }
     }
 
     if (n < 2) stop("a tree must have at least 2 tips.")
-    nbr <- 2 * n - 2
-    if (!rooted) nbr <- nbr - 1
+    nbr <- 2 * n - 3 + rooted
     edge <- matrix(NA, nbr, 2)
 
+    n <- as.integer(n)
     if (n == 2) {
-        if (rooted) edge[] <- c(3, 3, 1, 2)
+        if (rooted) edge[] <- c(3L, 3L, 1L, 2L)
         else stop("an unrooted tree must have at least 3 tips.")
     } else if (n == 3) {
         edge[] <-
-          if (rooted) c(4, 5, 5, 4, 5, 1:3)
-          else c(4, 4, 4, 1:3)
+          if (rooted) c(4L, 5L, 5L, 4L, 5L, 1:3)
+          else c(4L, 4L, 4L, 1:3)
     } else if (n == 4 && !rooted) {
-        edge[] <- c(5, 6, 6, 5, 5, 6, 1:4)
+        edge[] <- c(5L, 6L, 6L, 5L, 5L, 6L, 1:4)
     } else {
-        nod <- n + 1
+        nod <- n + 1L
         if (rooted) { # n > 3
             foo(n, 1)
             ## The following is slightly more efficient than affecting the
@@ -70,21 +70,21 @@ rtree <- function(n, rooted = TRUE, tip.label = NULL, br = runif, ...)
                 foo(n1, 2)
             } else if (n1 == 2) {
                 edge[2:3, 1] <- edge[1, 2] <- nod
-                nod <- nod + 1
+                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
-                nod <- nod + 1
+                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
-                ## nod <- nod + 1
+                ## nod <- nod + 1L
             }
             i <- which(is.na(edge[, 2]))
             edge[i, 2] <- 1:n
@@ -95,28 +95,31 @@ rtree <- function(n, rooted = TRUE, tip.label = NULL, br = runif, ...)
       if (is.null(tip.label)) paste("t", sample(n), sep = "")
       else sample(tip.label)
     if (is.function(br)) phy$edge.length <- br(nbr, ...)
-    phy$Nnode <- if (rooted) n - 1 else n - 2
+    phy$Nnode <- n - 2L + rooted
     class(phy) <- "phylo"
     phy
 }
 
-rcoal <- function(n, tip.label = NULL, br = rexp, ...)
+rcoal <- function(n, tip.label = NULL, br = "coalescent", ...)
 {
+    n <- as.integer(n)
     nbr <- 2*n - 2
     edge <- matrix(NA, nbr, 2)
-    x <- br(n - 1, ...) # coalescence times
+    ## coalescence times by default:
+    x <- if (is.character(br)) 2*rexp(n - 1)/(n:2 * (n - 1):1)
+    else br(n - 1, ...)
     if (n == 2) {
-        edge[] <- c(3, 3, 1:2)
+        edge[] <- c(3L, 3L, 1:2)
         edge.length <- rep(x, 2)
     } else if (n == 3) {
-        edge[] <- c(4, 5, 5, 4, 5, 1:3)
+        edge[] <- c(4L, 5L, 5L, 4L, 5L, 1:3)
         edge.length <- c(x[2], x[1], x[1], sum(x))
     } else {
         edge.length <- numeric(nbr)
         h <- numeric(2*n - 1) # initialized with 0's
         node.height <- cumsum(x)
         pool <- 1:n
-        nextnode <- 2*n - 1
+        nextnode <- 2L*n - 1L
         for (i in 1:(n - 1)) {
             y <- sample(pool, size = 2)
             ind <- (i - 1)*2 + 1:2
@@ -125,14 +128,14 @@ rcoal <- function(n, tip.label = NULL, br = rexp, ...)
             edge.length[ind] <- node.height[i] - h[y]
             h[nextnode] <- node.height[i]
             pool <- c(pool[! pool %in% y], nextnode)
-            nextnode <- nextnode - 1
+            nextnode <- nextnode - 1L
         }
     }
     phy <- list(edge = edge, edge.length = edge.length)
     phy$tip.label <-
       if (is.null(tip.label)) paste("t", 1:n, sep = "")
       else tip.label
-    phy$Nnode <- n - 1
+    phy$Nnode <- n - 1L
     class(phy) <- "phylo"
     ##reorder(phy)
     ## to avoid crossings when converting with as.hclust: