new compute.brtime() + minor change in dist.topo.R

git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@165 6e262413-ae40-0410-9e79-b911bd7a66b7

diff --git a/DESCRIPTION b/DESCRIPTION
index 4ea823664ab744ebfb3f630ec4acaa8f2f20915b..2089fab9c3acda403717513b610090c647b25430 100644 (file)
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: ape
 Version: 2.7-3
-Date: 2011-07-04
+Date: 2011-07-15
 Title: Analyses of Phylogenetics and Evolution
 Author: Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong, Richard Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel, Christoph Heibl, Daniel Lawson, Vincent Lefort, Pierre Legendre, Jim Lemon, Yvonnick Noel, Johan Nylander, Rainer Opgen-Rhein, Klaus Schliep, Korbinian Strimmer, Damien de Vienne
 Maintainer: Emmanuel Paradis <Emmanuel.Paradis@ird.fr>

diff --git a/NEWS b/NEWS
index 9e66a45982269e849ce15198fc94c2813a598e50..cfd28ab9e05dea767f55e37a6d18113b4b6a4c43 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -3,6 +3,8 @@
 
 NEW FEATURES
 
+    o The new function compute.brtime computes and sets branching times.
+
     o mantel.test() has a new argument 'alternative' which is
       "two-sided" by default. Previously, this test was one-tailed
       with no possibility to change.
@@ -21,6 +23,10 @@ BUG FIXES
     o Cross-validation in chronopl() did not work when 'age.max' was
       used.
 
+    o consensus(, p = 0.5) could return an incorrect tree if some
+      incompatible splits occur in 50% of the trees (especially with
+      small number of trees).
+
 
 
                CHANGES IN APE VERSION 2.7-2

diff --git a/R/compute.brtime.R b/R/compute.brtime.R
new file mode 100644 (file)
index 0000000..727fcb6
--- /dev/null
+++ b/R/compute.brtime.R
@@ -0,0 +1,46 @@
+## compute.brtime.R (2011-07-15)
+
+##   Compute and Set Branching Times
+
+## Copyright 2011 Emmanuel Paradis
+
+## This file is part of the R-package `ape'.
+## See the file ../COPYING for licensing issues.
+
+compute.brtime <-
+    function(phy, method = "coalescent", force.positive = NULL)
+{
+    if (!inherits(phy, "phylo"))
+        stop('object "phy" is not of class "phylo"')
+    n <- length(phy$tip.label)
+    m <- phy$Nnode
+    N <- Nedge(phy)
+
+    ## x: branching times (aka, node ages 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
+    } else if (is.numeric(method)) {
+        x <- as.vector(method)
+        if (length(x) != m)
+            stop("number of branching times given is not equal to the number of nodes")
+        if (is.null(force.positive))
+            force.positive <- FALSE
+    }
+
+    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
+    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]]
+
+    phy$edge.length <- el
+    reorder(phy)
+}
+

diff --git a/R/dist.topo.R b/R/dist.topo.R
index 455107ec924c43eab1527bbc52503b3fb3e72643..525178f2528ed1ba66321e178d8c79d2cf0c91ed 100644 (file)
--- a/R/dist.topo.R
+++ b/R/dist.topo.R
@@ -1,4 +1,4 @@
-## dist.topo.R (2011-06-14)
+## dist.topo.R (2011-07-13)
 
 ##      Topological Distances, Tree Bipartitions,
 ##   Consensus Trees, and Bootstrapping Phylogenies
@@ -268,6 +268,7 @@ consensus <- function(..., p = 1, check.labels = TRUE)
     ## Get all observed partitions and their frequencies:
     pp <- prop.part(obj, check.labels = FALSE)
     ## Drop the partitions whose frequency is less than 'p':
+    if (p == 0.5) p <- 0.5000001 # avoid incompatible splits
     pp <- pp[attr(pp, "number") >= p * ntree]
     ## Get the order of the remaining partitions by decreasing size:
     ind <- sort(unlist(lapply(pp, length)), decreasing = TRUE,

diff --git a/man/compute.brtime.Rd b/man/compute.brtime.Rd
new file mode 100644 (file)
index 0000000..8d0c8c4
--- /dev/null
+++ b/man/compute.brtime.Rd
@@ -0,0 +1,47 @@
+\name{compute.brtime}
+\alias{compute.brtime}
+\title{Compute and Set Branching Times}
+\description{
+  This function computes the branch lengths of a tree giving its
+  branching times (aka node ages or heights).
+}
+\usage{
+compute.brtime(phy, method = "coalescent", force.positive = NULL)
+}
+\arguments{
+  \item{phy}{an object of class \code{"phylo"}.}
+  \item{method}{either \code{"coalescent"} (the default), or a numeric
+    vector giving the branching times.}
+  \item{force.positive}{a logical value (see details).}
+}
+\details{
+  By default, a set of random branching times is generated from a simple
+  coalescent, and the option \code{force.positive} is set to \code{TRUE}
+  so that no branch length is negative.
+
+  If a numeric vector is passed to \code{method}, it is taken as the
+  branching times of the nodes with respect to their numbers (i.e., the
+  first element of \code{method} is the branching time of the node
+  numbered \eqn{n + 1} [= the root], the second element of the node
+  numbered \eqn{n + 2}, and so on), so \code{force.positive} is set to
+  \code{FALSE}. This may result in negative branch lengths. To avoid
+  this, one should use \code{force.positive = TRUE} in which case the
+  branching times are eventually reordered.
+}
+\value{
+  An object of class \code{"phylo"} with branch lengths and ultrametric.
+}
+\author{Emmanuel Paradis}
+\seealso{
+  \code{\link{compute.brlen}}, \code{\link{branching.times}}
+}
+\examples{
+tr <- rtree(10)
+layout(matrix(1:4, 2))
+plot(compute.brtime(tr)); axisPhylo()
+plot(compute.brtime(tr, force.positive = FALSE)); axisPhylo()
+plot(compute.brtime(tr, 1:9)); axisPhylo() # a bit nonsense
+plot(compute.brtime(tr, 1:9, TRUE)); axisPhylo()
+layout(1)
+}
+\keyword{manip}