From: paradis <paradis@6e262413-ae40-0410-9e79-b911bd7a66b7>
Date: Thu, 6 May 2010 16:48:18 +0000 (+0000)
Subject: 2nd try!!
X-Git-Url: https://git.donarmstrong.com/?p=ape.git;a=commitdiff_plain;h=06c3113db74a7cfa54c15a6f18163cd9b2c1f6db

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 430c210..0e80212 100644
--- 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 8fdb5c1..aa5fd01 100644
--- 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 f728ee7..9ea6590 100644
--- 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 1eb770a..779db02 100644
--- 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 642591e..8cc48e1 100644
--- 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))
         }
     }