R/compute.brtime.R

   1 ## compute.brtime.R (2012-03-02)
   2
   3 ##   Compute and Set Branching Times
   4
   5 ## Copyright 2011-2012 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 compute.brtime <-
  11     function(phy, method = "coalescent", force.positive = NULL)
  12 {
  13     if (!inherits(phy, "phylo"))
  14         stop('object "phy" is not of class "phylo"')
  15     n <- length(phy$tip.label)
  16     m <- phy$Nnode
  17     N <- Nedge(phy)
  18
  19     ## x: branching times (aka, node ages, depths, or heights)
  20
  21     if (identical(method, "coalescent")) { # the default
  22         x <- 2 * rexp(m)/(as.double((m + 1):2) * as.double(m:1))
  23         ## x <- 2 * rexp(n - 1)/(as.double(n:2) * as.double((n - 1):1))
  24         if (is.null(force.positive)) force.positive <- TRUE
  25     } else if (is.numeric(method)) {
  26         x <- as.vector(method)
  27         if (length(x) != m)
  28             stop("number of branching times given is not equal to the number of nodes")
  29         if (is.null(force.positive))
  30             force.positive <- FALSE
  31     }
  32
  33     y <- c(rep(0, n), x) # for all nodes (terminal and internal)
  34
  35     e1 <- phy$edge[, 1L] # local copies of the pointers
  36     e2 <- phy$edge[, 2L] #
  37
  38     if (force.positive) {
  39         o <- .Call("seq_root2tip", phy$edge, n, m, PACKAGE = "ape")
  40         list.nodes <- list(n + 1L)
  41         i <- 2L
  42         repeat {
  43             z <- sapply(o, "[", i)
  44             z <- unique(z[!(z <= n | is.na(z))])
  45             if (!length(z)) break
  46             list.nodes[[i]] <- z
  47             i <- i + 1L
  48         }
  49         nodes <- unlist(lapply(list.nodes, function(x) x[sample(length(x))]))
  50         y[nodes] <- sort(x, decreasing = TRUE)
  51     }
  52
  53     phy$edge.length <- y[e1] - y[e2]
  54     phy
  55 }