changed rtree and rcoal

[ape.git] / R / rtree.R

diff --git a/R/rtree.R b/R/rtree.R
index 9425dadf22f5df689990e1ffe6a3829e3ce99c6b..030b47e77885a70a4cc1b362d461878133ce8dc5 100644 (file)
--- a/R/rtree.R
+++ b/R/rtree.R
@@ -1,8 +1,8 @@
-## rtree.R (2008-05-07)
+## 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.
@@ -94,7 +94,10 @@ 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"
     phy
@@ -106,14 +109,14 @@ 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)
     } else if (n == 3) {
         edge[] <- c(4L, 5L, 5L, 4L, 5L, 1:3)
-        edge.length <- c(x[2], x[1], x[1], sum(x))
+        edge.length <- c(x[c(2, 1, 1)], sum(x))
     } else {
         edge.length <- numeric(nbr)
         h <- numeric(2*n - 1) # initialized with 0's
@@ -132,14 +135,15 @@ rcoal <- function(n, tip.label = NULL, br = "coalescent", ...)
         }
     }
     phy <- list(edge = edge, edge.length = edge.length)
-    phy$tip.label <-
-      if (is.null(tip.label)) paste("t", 1:n, sep = "")
-      else tip.label
+    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, ...)