2nd try!!

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

diff --git a/ChangeLog b/ChangeLog
index 430c21093791aab51038f32d86ce4575df8f8698..0e8021260d28bed95db39a36b15ab75bd782feb2 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,10 +1,20 @@
                CHANGES IN APE VERSION 2.5-2
 
 
+NEW FEATURES
+
+    o There is now a print method for results from ace().
+
+
 BUG FIXES
 
     o as.phylo.hclust() used to multiply edge lengths by 2.
 
+    o A minor bug was fixed in rTraitDisc().
+
+    o ace() sometimes failed (parameter value was NaN and the optimisation
+      failed).
+
 
 DEPRECATED & DEFUNCT
 

diff --git a/DESCRIPTION b/DESCRIPTION
index 8fdb5c182fb0bc37f4d5c3631472c02a522bb0a7..aa5fd013f129c10eca02e679e26bd7fbc883ae6a 100644 (file)
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: ape
 Version: 2.5-2
-Date: 2010-04-21
+Date: 2010-05-06
 Title: Analyses of Phylogenetics and Evolution
 Author: Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong, Richard Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel, Gangolf Jobb, Christoph Heibl, Daniel Lawson, Vincent Lefort, Pierre Legendre, Jim Lemon, Yvonnick Noel, Johan Nylander, Rainer Opgen-Rhein, Korbinian Strimmer, Damien de Vienne
 Maintainer: Emmanuel Paradis <Emmanuel.Paradis@ird.fr>

diff --git a/R/ace.R b/R/ace.R
index f728ee7c6045f9478e21fa0edd1e12fafa5aa4fa..9ea659001688562e42f6e6f557b22720ad0830b2 100644 (file)
--- a/R/ace.R
+++ b/R/ace.R
@@ -1,4 +1,4 @@
-## ace.R (2010-02-03)
+## ace.R (2010-04-23)
 
 ##   Ancestral Character Estimation
 
@@ -157,9 +157,10 @@ as the number of categories in `x'")
         phy <- reorder(phy, "pruningwise")
 
         Q <- matrix(0, nl, nl)
-        ## from Rich FitzJohn:
-        comp <- numeric(nb.tip + nb.node) # Storage...
         dev <- function(p, output.liks = FALSE) {
+            if (any(is.nan(p)) || any(is.infinite(p))) return(1e50)
+            ## from Rich FitzJohn:
+            comp <- numeric(nb.tip + nb.node) # Storage...
             Q[] <- c(p, 0)[rate]
             diag(Q) <- -rowSums(Q)
             for (i  in seq(from = 1, by = 2, length.out = nb.node)) {
@@ -177,7 +178,7 @@ as the number of categories in `x'")
             -2 * sum(log(comp[-TIPS]))
         }
         out <- nlminb(rep(ip, length.out = np), function(p) dev(p),
-                      lower = rep(0, np), upper = rep(Inf, np))
+                      lower = rep(0, np), upper = rep(1e50, np))
         obj$loglik <- -out$objective/2
         obj$rates <- out$par
         oldwarn <- options("warn")
@@ -229,3 +230,24 @@ anova.ace <- function(object, ...)
     structure(table, heading = "Likelihood Ratio Test Table",
               class = c("anova", "data.frame"))
 }
+
+print.ace <- function(x, digits = 4, ...)
+{
+    cat("\n    Ancestral Character Reconstruction\n\n")
+    cat("Call: ")
+    print(x$call)
+    cat("\n")
+    cat("    Log-likelihood:", x$loglik, "\n\n")
+    cat("Rate index matrix:\n")
+    ratemat <- x$index.matrix
+    dimnames(ratemat)[1:2] <- dimnames(x$lik.anc)[2]
+    print(ratemat, na.print = ".")
+    cat("\n")
+    npar <- length(x$rates)
+    estim <- data.frame(1:npar, round(x$rates, digits), round(x$se, digits))
+    cat("Parameter estimates:\n")
+    names(estim) <- c("rate index", "estimate", "std-err")
+    print(estim, row.names = FALSE)
+    cat("\nScaled likelihoods at the root (type 'x$lik.anc' to get them for all nodes):\n")
+    print(x$lik.anc[1, ])
+}

diff --git a/R/dist.topo.R b/R/dist.topo.R
index 1eb770af6a12f15948b70388a322eb0f6fcd7992..779db02b2f579e83277afd38fd0d72aab3a2c260 100644 (file)
--- a/R/dist.topo.R
+++ b/R/dist.topo.R
@@ -1,4 +1,4 @@
-## dist.topo.R (2010-01-25)
+## dist.topo.R (2010-05-06)
 
 ##      Topological Distances, Tree Bipartitions,
 ##   Consensus Trees, and Bootstrapping Phylogenies
@@ -47,7 +47,6 @@ dist.topo <- function(x, y, method = "PH85")
         for (i in 2:q1) {
             for (j in 2:q2) {
                 if (identical(bp1[[i]], bp2[[j]]) | identical(bp1[[i]], bp2.comp[[j]])) {
-                    if (i == 19) browser()
                     dT <- dT + (x$edge.length[which(x$edge[, 2] == nx + i)] -
                                 y$edge.length[which(y$edge[, 2] == ny + j)])^2
                     found1 <- found2[j] <- TRUE

diff --git a/R/rTrait.R b/R/rTrait.R
index 642591e6e2523d2ac0da057b861f1831546f6621..8cc48e17c9e45296c082b33ca33f29e01c3f3296 100644 (file)
--- a/R/rTrait.R
+++ b/R/rTrait.R
@@ -1,4 +1,4 @@
-## rTrait.R (2010-02-03)
+## rTrait.R (2010-05-06)
 
 ##   Trait Evolution
 
@@ -58,10 +58,12 @@ rTraitDisc <-
         environment(model) <- environment() # to find 'k'
         for (i in N:1) x[des[i]] <- model(x[anc[i]], el[i])
     } else {
-        diag(Q) <- -rowSums(Q)
         freq <- rep(freq, each = k)
+        Q <- Q * freq
+        diag(Q) <- 0
+        diag(Q) <- -rowSums(Q)
         for (i in N:1) {
-            p <- matexpo(Q * freq * el[i])[x[anc[i]], ]
+            p <- matexpo(Q * el[i])[x[anc[i]], ]
             x[des[i]] <- .Internal(sample(k, size = 1, FALSE, prob = p))
         }
     }