From: paradis <paradis@6e262413-ae40-0410-9e79-b911bd7a66b7>
Date: Fri, 15 Jul 2011 02:18:47 +0000 (+0000)
Subject: new compute.brtime() + minor change in dist.topo.R
X-Git-Url: https://git.donarmstrong.com/?p=ape.git;a=commitdiff_plain;h=1f8c4f7bf4f7100178bd11f274247ae9fac1b876

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 4ea8236..2089fab 100644
--- 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 9e66a45..cfd28ab 100644
--- 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
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 455107e..525178f 100644
--- 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
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}