improved ace()

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

diff --git a/ChangeLog b/ChangeLog
index 8ad799a009ac2d7c8469030d915b855691453cfa..5705d9abd1fbd4e8ced94b76af8e0ecb5510d468 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -22,6 +22,10 @@ BUG FIXES
 
     o A small bug was fixed in CDAM.global().
 
+    o ace() failed with large data sets. Thanks to Rich FitzJohn for
+      the fix. With other improvements, this function is now about 6
+      times faster.
+
 
 OTHER CHANGES
 

diff --git a/DESCRIPTION b/DESCRIPTION
index fac51bde1208d73aef86d6e7c4cdaf02bc9a1af0..c12694c389d07fb849f7c9f72d6e41311c418f48 100644 (file)
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: ape
 Version: 2.3-1
-Date: 2009-06-08
+Date: 2009-06-10
 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 b45a88e2009cc5ee252c2e1e1246c6db0162423f..91994bed1ffeb9a35fd2dfdc9a9bc20aa2d11384 100644 (file)
--- a/R/ace.R
+++ b/R/ace.R
@@ -1,4 +1,4 @@
-## ace.R (2009-05-10)
+## ace.R (2009-06-10)
 
 ##   Ancestral Character Estimation
 
@@ -149,10 +149,13 @@ as the number of categories in `x'")
         rate[rate == 0] <- np + 1 # to avoid 0's since we will use this an numeric indexing
 
         liks <- matrix(0, nb.tip + nb.node, nl)
-        for (i in 1:nb.tip) liks[i, x[i]] <- 1
+        TIPS <- 1:nb.tip
+        liks[cbind(TIPS, x)] <- 1
         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) {
             Q[] <- c(p, 0)[rate]
             diag(Q) <- -rowSums(Q)
@@ -161,14 +164,14 @@ as the number of categories in `x'")
                 anc <- phy$edge[i, 1]
                 des1 <- phy$edge[i, 2]
                 des2 <- phy$edge[j, 2]
-                tmp <- eigen(Q * phy$edge.length[i], symmetric = FALSE)
-                P1 <- tmp$vectors %*% diag(exp(tmp$values)) %*% solve(tmp$vectors)
-                tmp <- eigen(Q * phy$edge.length[j], symmetric = FALSE)
-                P2 <- tmp$vectors %*% diag(exp(tmp$values)) %*% solve(tmp$vectors)
-                liks[anc, ] <- P1 %*% liks[des1, ] * P2 %*% liks[des2, ]
+                v.l <- matexpo(Q * phy$edge.length[i]) %*% liks[des1, ]
+                v.r <- matexpo(Q * phy$edge.length[j]) %*% liks[des2, ]
+                v <- v.l * v.r
+                comp[anc] <- sum(v)
+                liks[anc, ] <- v/comp[anc]
             }
-            if (output.liks) return(liks[-(1:nb.tip), ])
-            - 2 * log(sum(liks[nb.tip + 1, ]))
+            if (output.liks) return(liks[-TIPS, ])
+            -2 * sum(log(comp[-TIPS]))
         }
         out <- nlminb(rep(ip, length.out = np), function(p) dev(p),
                       lower = rep(0, np), upper = rep(Inf, np))
@@ -184,10 +187,8 @@ as the number of categories in `x'")
         else obj$se <- sqrt(diag(solve(h)))
         obj$index.matrix <- index.matrix
         if (CI) {
-            lik.anc <- dev(obj$rates, TRUE)
-            lik.anc <- lik.anc / rowSums(lik.anc)
-            colnames(lik.anc) <- lvls
-            obj$lik.anc <- lik.anc
+            obj$lik.anc <- dev(obj$rates, TRUE)
+            colnames(obj$lik.anc) <- lvls
         }
     }
     obj$call <- match.call()

diff --git a/Thanks b/Thanks
index 57e3cb259be8cb042aa554fa3149471ed787837d..157d4cc8dcee296a0d44a45c1c84655e5ff0d0e0 100644 (file)
--- a/Thanks
+++ b/Thanks
@@ -5,7 +5,7 @@ Many users gave important feed-back with their encouragements,
 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,
+Éric Durand, Olivier François, Rich FitzJohn Bret Larget, Nick Matzke,
 Elizabeth Purdom, Dan Rabosky, Klaus Schliep, Li-San Wang, Yan Wong,
 and Peter Wragg. Contact me if I forgot someone.