some changes for ape 3.0-6

[ape.git] / R / compute.brtime.R

diff --git a/R/compute.brtime.R b/R/compute.brtime.R
index 727fcb6a37dc9fb261f8e4187d3f3d441a7293ed..24e9b83226cdd057dda1fa3113febedb09dce2c3 100644 (file)
--- a/R/compute.brtime.R
+++ b/R/compute.brtime.R
@@ -1,8 +1,8 @@
-## compute.brtime.R (2011-07-15)
+## compute.brtime.R (2012-03-02)
 
 ##   Compute and Set Branching Times
 
-## Copyright 2011 Emmanuel Paradis
+## Copyright 2011-2012 Emmanuel Paradis
 
 ## This file is part of the R-package `ape'.
 ## See the file ../COPYING for licensing issues.
@@ -16,12 +16,12 @@ compute.brtime <-
     m <- phy$Nnode
     N <- Nedge(phy)
 
-    ## x: branching times (aka, node ages or heights)
+    ## x: branching times (aka, node ages, depths, or heights)
 
     if (identical(method, "coalescent")) { # the default
-        x <- 2 * rexp(n - 1)/(as.double(n:2) * as.double((n - 1):1))
-        if (is.null(force.positive))
-            force.positive <- TRUE
+        x <- 2 * rexp(m)/(as.double((m + 1):2) * as.double(m:1))
+        ## x <- 2 * rexp(n - 1)/(as.double(n:2) * as.double((n - 1):1))
+        if (is.null(force.positive)) force.positive <- TRUE
     } else if (is.numeric(method)) {
         x <- as.vector(method)
         if (length(x) != m)
@@ -32,15 +32,24 @@ compute.brtime <-
 
     y <- c(rep(0, n), x) # for all nodes (terminal and internal)
 
-    phy <- reorder(phy, "pruningwise")
-    e1 <- phy$edge[, 1L] # local copies of the pointer
+    e1 <- phy$edge[, 1L] # local copies of the pointers
     e2 <- phy$edge[, 2L] #
-    el <- numeric(N)
-    if (force.positive) y[unique(e1)] <- sort(x)
 
-    for (i in 1:N) el[i] <- y[e1[i]] - y[e2[i]]
+    if (force.positive) {
+        o <- .Call("seq_root2tip", phy$edge, n, m, PACKAGE = "ape")
+        list.nodes <- list(n + 1L)
+        i <- 2L
+        repeat {
+            z <- sapply(o, "[", i)
+            z <- unique(z[!(z <= n | is.na(z))])
+            if (!length(z)) break
+            list.nodes[[i]] <- z
+            i <- i + 1L
+        }
+        nodes <- unlist(lapply(list.nodes, function(x) x[sample(length(x))]))
+        y[nodes] <- sort(x, decreasing = TRUE)
+    }
 
-    phy$edge.length <- el
-    reorder(phy)
+    phy$edge.length <- y[e1] - y[e2]
+    phy
 }
-