X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2Frtree.R;h=030b47e77885a70a4cc1b362d461878133ce8dc5;hb=b41db3e6053bfb0e4b9000d5de678d65a1af9670;hp=cba5d4e4a40c967dbbb443d30396ab7333913aef;hpb=d88302b4735b5b7c9132387090bb592d906fe1cb;p=ape.git

diff --git a/R/rtree.R b/R/rtree.R
index cba5d4e..030b47e 100644
--- a/R/rtree.R
+++ b/R/rtree.R
@@ -1,8 +1,8 @@
-## rtree.R (2008-01-13)
+## rtree.R (2009-11-03)
 
 ##   Generates Random Trees
 
-## Copyright 2004-2008 Emmanuel Paradis
+## Copyright 2004-2009 Emmanuel Paradis
 
 ## This file is part of the R-package `ape'.
 ## See the file ../COPYING for licensing issues.
@@ -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
@@ -94,31 +94,35 @@ 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, ...)
-    phy$Nnode <- if (rooted) n - 1 else n - 2
+    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"
     phy
 }
 
 rcoal <- function(n, tip.label = NULL, br = "coalescent", ...)
 {
+    n <- as.integer(n)
     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(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.length <- c(x[2], x[1], x[1], sum(x))
+        edge[] <- c(4L, 5L, 5L, 4L, 5L, 1:3)
+        edge.length <- c(x[c(2, 1, 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
@@ -127,18 +131,19 @@ rcoal <- function(n, tip.label = NULL, br = "coalescent", ...)
             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
+    if (is.null(tip.label))
+        tip.label <- paste("t", 1:n, sep = "")
+    phy$tip.label <- sample(tip.label)
+    phy$Nnode <- n - 1L
     class(phy) <- "phylo"
-    ##reorder(phy)
+    phy <- reorder(phy)
     ## to avoid crossings when converting with as.hclust:
-    read.tree(text = write.tree(phy))
+    phy$edge[phy$edge[, 2] <= n, 2] <- 1:n
+    phy
 }
 
 rmtree <- function(N, n, rooted = TRUE, tip.label = NULL, br = runif, ...)