change nlm to nlminb in ace() + new makeNodeLabel + fixed drop.tip

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

diff --git a/ChangeLog b/ChangeLog
index b6bffe8687b5bf9a7dc248c22ae7d5b7ed99725f..d659113a7cd989de0357b1303000d0a47674b66b 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -10,13 +10,22 @@ NEW FEATURES
     o The new function yule.time fits a user-defined time-dependent
       Yule model by maximum likelihood.
 
     o The new function yule.time fits a user-defined time-dependent
       Yule model by maximum likelihood.
 

+    o The new function makeNodeLabel creates and/or modifies node
+      labels in a flexible way.
+

 
 BUG FIXES
 
     o drop.tip() shuffled tip labels in some cases.
 
 
 BUG FIXES
 
     o drop.tip() shuffled tip labels in some cases.
 

+    o drop.tip() did not handle node.label correctly.
+

     o is.ultrametric() now checks the ordering of the edge matrix.
 
     o is.ultrametric() now checks the ordering of the edge matrix.
 

+    o ace() sometimes returned negative values of likelihoods of
+      ancestral states (thanks to Dan Rabosky for solving this long
+      lasting bug).
+

 
 OTHER CHANGES
 
 
 OTHER CHANGES
 

diff --git a/DESCRIPTION b/DESCRIPTION
index 1cc4f4a5db33950648020404c37519a60709e038..51c4cfc7521bf5690602170f13938196171210d9 100644 (file)
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: ape
 Version: 2.3
 Package: ape
 Version: 2.3

-Date: 2009-03-10
+Date: 2009-03-22

 Title: Analyses of Phylogenetics and Evolution
 Author: Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong,
   Richard Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel,
 Title: Analyses of Phylogenetics and Evolution
 Author: Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong,
   Richard Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel,

diff --git a/R/ace.R b/R/ace.R
index c38faeb6bdb2ef2b8b489139a841b3b5a440adaf..0439600832a0f1018af4b75974ea11703f2a6e02 100644 (file)
--- a/R/ace.R
+++ b/R/ace.R
@@ -1,4 +1,4 @@
-## ace.R (2009-01-19)
+## ace.R (2009-03-22)

 
 ##     Ancestral Character Estimation
 
 
 ##     Ancestral Character Estimation
 


@@ -170,13 +170,18 @@ as the number of categories in `x'")
             if (output.liks) return(liks[-(1:nb.tip), ])
             - 2 * log(sum(liks[nb.tip + 1, ]))
         }
             if (output.liks) return(liks[-(1:nb.tip), ])
             - 2 * log(sum(liks[nb.tip + 1, ]))
         }

-        out <- nlm(function(p) dev(p), p = rep(ip, length.out = np),
-                   hessian = TRUE)
-        obj$loglik <- -out$minimum / 2
-        obj$rates <- out$estimate
-        if (any(out$gradient == 0))
+        out <- nlminb(rep(ip, length.out = np), function(p) dev(p),
+                      lower = rep(0, np), upper = rep(Inf, np))
+        obj$loglik <- -out$objective/2
+        obj$rates <- out$par
+        oldwarn <- options("warn")
+        options(warn = -1)
+        h <- nlm(function(p) dev(p), p = obj$rates, iterlim = 1,
+                 stepmax = 0, hessian = TRUE)$hessian
+        options(oldwarn)
+        if (any(h == 0))

           warning("The likelihood gradient seems flat in at least one dimension (gradient null):\ncannot compute the standard-errors of the transition rates.\n")
           warning("The likelihood gradient seems flat in at least one dimension (gradient null):\ncannot compute the standard-errors of the transition rates.\n")

-        else obj$se <- sqrt(diag(solve(out$hessian)))
+        else obj$se <- sqrt(diag(solve(h)))

         obj$index.matrix <- index.matrix
         if (CI) {
             lik.anc <- dev(obj$rates, TRUE)
         obj$index.matrix <- index.matrix
         if (CI) {
             lik.anc <- dev(obj$rates, TRUE)

diff --git a/R/drop.tip.R b/R/drop.tip.R
index 69c957629e6533598488d0f7b3c4ede773faede7..c0bda8bc15042186f101ac7f548e7fc1fb20dd34 100644 (file)
--- a/R/drop.tip.R
+++ b/R/drop.tip.R
@@ -1,4 +1,4 @@
-## drop.tip.R (2009-03-04)
+## drop.tip.R (2009-03-22)

 
 ##   Remove Tips in a Phylogenetic Tree
 
 
 ##   Remove Tips in a Phylogenetic Tree
 


@@ -76,6 +76,7 @@ drop.tip <-
 {
     if (class(phy) != "phylo")
         stop('object "phy" is not of class "phylo"')
 {
     if (class(phy) != "phylo")
         stop('object "phy" is not of class "phylo"')

+    phy <- reorder(phy)

     Ntip <- length(phy$tip.label)
     NEWROOT <- ROOT <- Ntip + 1
     Nnode <- phy$Nnode
     Ntip <- length(phy$tip.label)
     NEWROOT <- ROOT <- Ntip + 1
     Nnode <- phy$Nnode


@@ -145,7 +146,10 @@ drop.tip <-
     TIPS <- phy$edge[, 2] <= Ntip
     ## keep the ordering so no need to reorder tip.label:
     phy$edge[TIPS, 2] <- rank(phy$edge[TIPS, 2])
     TIPS <- phy$edge[, 2] <= Ntip
     ## keep the ordering so no need to reorder tip.label:
     phy$edge[TIPS, 2] <- rank(phy$edge[TIPS, 2])

-    Ntip <- length(phy$tip.label) # update Ntip
+    ## 3) update node.label if needed
+    if (!is.null(phy$node.label))
+        phy$node.label <- phy$node.label[sort(unique(phy$edge[, 1])) - Ntip]
+    Ntip <- length(phy$tip.label) # 4) update Ntip

 
     ## make new tip labels if necessary
     if (subtree || !trim.internal) {
 
     ## make new tip labels if necessary
     if (subtree || !trim.internal) {

diff --git a/R/makeNodeLabel.R b/R/makeNodeLabel.R
new file mode 100644 (file)
index 0000000..d679f5e
--- /dev/null
+++ b/R/makeNodeLabel.R
@@ -0,0 +1,61 @@
+## makeNodeLabel.R (2009-03-22)
+
+##   Makes Node Labels
+
+## Copyright 2009 Emmanuel Paradis
+
+## This file is part of the R-package `ape'.
+## See the file ../COPYING for licensing issues.
+
+makeNodeLabel <- function(phy, method = "number", prefix = "Node",
+                          nodeList = list(), ...)
+{
+    method <- sapply(method, match.arg, c("number", "md5sum", "user"),
+                     USE.NAMES = FALSE)
+
+    if ("number" %in% method)
+        phy$node.label <- paste(prefix, 1:phy$Nnode, sep = "")
+
+    if ("md5sum" %in% method) {
+        nl <- character(phy$Nnode)
+        pp <- prop.part(phy, check.labels = FALSE)
+        labs <- attr(pp, "labels")
+        fl <- tempfile()
+        for (i in seq_len(phy$Nnode)) {
+            cat(sort(labs[pp[[i]]]), sep = "\n", file = fl)
+            nl[i] <- tools::md5sum(fl)
+        }
+        unlink(fl)
+        phy$node.label <- nl
+    }
+
+    if ("user" %in% method) {
+        if (is.null(phy$node.label))
+            phy$node.label <- character(phy$Nnode)
+        nl <- names(nodeList)
+        if (is.null(nl)) stop("argument 'nodeList' has no names")
+        Ntip <- length(phy$tip.label)
+        seq.nod <- .Call("seq_root2tip", phy$edge, Ntip, phy$Nnode,
+                         PACKAGE = "ape")
+        ## a local version to avoid the above call many times:
+        .getMRCA <- function(seq.nod, tip) {
+            sn <- seq.nod[tip]
+            MRCA <- Ntip + 1
+            i <- 2
+            repeat {
+                x <- unique(unlist(lapply(sn, "[", i)))
+                if (length(x) != 1) break
+                MRCA <- x
+                i <- i + 1
+            }
+            MRCA
+        }
+        for (i in seq_along(nodeList)) {
+            tips <- sapply(nodeList[[i]], grep, phy$tip.label, ...,
+                           USE.NAMES = FALSE)
+            j <- .getMRCA(seq.nod, unique(unlist(tips)))
+            phy$node.label[j - Ntip] <- nl[i]
+        }
+    }
+    phy
+}

diff --git a/Thanks b/Thanks
index 8f93403e4748d0fc9b8329e11b75a3381835ba65..57e3cb259be8cb042aa554fa3149471ed787837d 100644 (file)
--- a/Thanks
+++ b/Thanks
@@ -6,8 +6,8 @@ comments, or bug reports: thanks to all of you!
 
 Significant bug fixes were provided by Cécile Ané, James Bullard,
 Éric Durand, Olivier François, Bret Larget, Nick Matzke,
 
 Significant bug fixes were provided by Cécile Ané, James Bullard,
 Éric Durand, Olivier François, Bret Larget, Nick Matzke,

-Elizabeth Purdom, Klaus Schliep, Li-San Wang, Yan Wong, and Peter
-Wragg. Contact me if I forgot someone.
+Elizabeth Purdom, Dan Rabosky, Klaus Schliep, Li-San Wang, Yan Wong,
+and Peter Wragg. Contact me if I forgot someone.

 
 Kurt Hornik, of the R Core Team, helped in several occasions to
 fix some problems and bugs.
 
 Kurt Hornik, of the R Core Team, helped in several occasions to
 fix some problems and bugs.

diff --git a/man/ace.Rd b/man/ace.Rd
index 85cfb210bbbea57af3f7f0786d6b4f4028ea8e89..73934dfb303498aacbdcaf7a2fe78bc536f7b532 100644 (file)
--- a/man/ace.Rd
+++ b/man/ace.Rd
@@ -41,7 +41,7 @@ ace(x, phy, type = "continuous", method = "ML", CI = TRUE,
   \item{k}{a numeric value giving the penalty per estimated parameter;
     the default is \code{k = 2} which is the classical Akaike
     information criterion.}
   \item{k}{a numeric value giving the penalty per estimated parameter;
     the default is \code{k = 2} which is the classical Akaike
     information criterion.}

-  \item{...}{further arguments passed to or from other methods.}
+  \item{\dots}{further arguments passed to or from other methods.}

 }
 \description{
   This function estimates ancestral character states, and the associated
 }
 \description{
   This function estimates ancestral character states, and the associated

diff --git a/man/makeNodeLabel.Rd b/man/makeNodeLabel.Rd
new file mode 100644 (file)
index 0000000..501f405
--- /dev/null
+++ b/man/makeNodeLabel.Rd
@@ -0,0 +1,69 @@
+\name{makeNodeLabel}
+\alias{makeNodeLabel}
+\title{Makes Node Labels}
+\description{
+  This function makes node labels in a tree in a flexible way.
+}
+
+\usage{
+makeNodeLabel(phy, method = "number", prefix = "Node", nodeList = list(), ...)
+}
+\arguments{
+  \item{phy}{an object of class \code{"phylo"}.}
+  \item{method}{a character string giving the method used to create the
+    labels. Three choices are possible: \code{"number"} (the default),
+    \code{"md5sum"}, and \code{"user"}, or any unambiguous abbreviation
+    of these.}
+  \item{prefix}{the prefix used if \code{method = "number"}.}
+  \item{nodeList}{a named list specifying how nodes are names if
+    \code{method = "user"} (see details and examples).}
+   \item{\dots}{further arguments passed to \code{grep}.}
+}
+\details{
+  The three methods are described below:
+
+  \itemize{
+    \item{``number''}{The labels are created with 1, 2, \dots suffixed
+      with the argument \code{prefix}; thus the default is to have
+      Node1, Node2, \dots Set \code{prefix = ""} to have only numbers.}
+    \item{``md5sum''}{For each node, the labels of the tips descendant
+      from this node are extracted, sorted alphabetically, and written
+      into a temporary file, then the md5sum of this file is extracted
+      and used as label. This results in a 32-character string which is
+      unique (even accross trees) for a given set of tip labels.}
+    \item{``user''}{the argument \code{nodeList} must be a list with
+      names, the latter will be used as node labels. For each element of
+      \code{nodeList}, the tip labels of the tree are searched for
+      patterns present in this element: this is done using
+      \code{\link[base]{grep}}. Then the most recent common ancestor of
+      the matching tips is given the corresponding names as labels. This
+      is repeated for each element of \code{nodeList}.}
+  }
+
+  The method \code{"user"} can be used in combination with either of the
+  two others (see examples). Note that this method only modifies the
+  specified node labels (so that if the other nodes have already labels
+  they are not modified) while the two others change all labels.
+}
+
+\author{Emmanuel Paradis \email{Emmanuel.Paradis@mpl.ird.fr}}
+\seealso{
+  \code{\link{makeLabel}}, \code{\link[base]{grep}}
+}
+\examples{
+tr <-
+"((Pan_paniscus,Pan_troglodytes),((Homo_sapiens,Hom_erectus),Homo_abilis));"
+tr <- read.tree(text = tr)
+tr <- makeNodeLabel(tr, "u", nodeList = list(Pan = "Pan", Homo = "Homo"))
+plot(tr, show.node.label = TRUE)
+### does not erase the previous node labels:
+tr <- makeNodeLabel(tr, "u", nodeList = list(Hominid = c("Pan","Homo")))
+plot(tr, show.node.label = TRUE)
+### the two previous commands could be combined:
+L <- list(Pan = "Pan", Homo = "Homo", Hominid = c("Pan","Homo"))
+tr <- makeNodeLabel(tr, "u", nodeList = L)
+### combining different methods:
+tr <- makeNodeLabel(tr, c("n", "u"), prefix = "#", nodeList = list(Hominid = c("Pan","Homo")))
+plot(tr, show.node.label = TRUE)
+}
+\keyword{manip}