some modifs for ape 3.0-8

[ape.git] / R / ace.R

diff --git a/R/ace.R b/R/ace.R
index 55e32c35ebc1e3c2580ce0af1be4d04650c53de4..2b9230ec7c9d43e7063ea9add8f24e8cb87a88e3 100644 (file)
--- a/R/ace.R
+++ b/R/ace.R
@@ -1,8 +1,8 @@
-## ace.R (2012-06-25)
+## ace.R (2013-01-31)
 
 ##   Ancestral Character Estimation
 
-## Copyright 2005-2012 Emmanuel Paradis and Ben Bolker
+## Copyright 2005-2013 Emmanuel Paradis and Ben Bolker
 
 ## This file is part of the R-package `ape'.
 ## See the file ../COPYING for licensing issues.
@@ -21,7 +21,8 @@
 
 ace <- function(x, phy, type = "continuous", method = "ML", CI = TRUE,
                 model = if (type == "continuous") "BM" else "ER",
-                scaled = TRUE, kappa = 1, corStruct = NULL, ip = 0.1)
+                scaled = TRUE, kappa = 1, corStruct = NULL, ip = 0.1,
+                use.expm = FALSE)
 {
     if (!inherits(phy, "phylo"))
         stop('object "phy" is not of class "phylo"')
@@ -190,6 +191,8 @@ ace <- function(x, phy, type = "continuous", method = "ML", CI = TRUE,
         liks[cbind(TIPS, x)] <- 1
         phy <- reorder(phy, "pruningwise")
 
+        E <- if (use.expm) expm::expm else ape::matexpo
+
         Q <- matrix(0, nl, nl)
         dev <- function(p, output.liks = FALSE) {
             if (any(is.nan(p)) || any(is.infinite(p))) return(1e50)
@@ -202,8 +205,8 @@ ace <- function(x, phy, type = "continuous", method = "ML", CI = TRUE,
                 anc <- phy$edge[i, 1]
                 des1 <- phy$edge[i, 2]
                 des2 <- phy$edge[j, 2]
-                v.l <- matexpo(Q * phy$edge.length[i]) %*% liks[des1, ]
-                v.r <- matexpo(Q * phy$edge.length[j]) %*% liks[des2, ]
+                v.l <- E(Q * phy$edge.length[i]) %*% liks[des1, ]
+                v.r <- E(Q * phy$edge.length[j]) %*% liks[des2, ]
                 v <- v.l * v.r
                 comp[anc] <- sum(v)
                 liks[anc, ] <- v/comp[anc]