new files for ape 3.0

[ape.git] / R / chronopl.R

diff --git a/R/chronopl.R b/R/chronopl.R
index 6a0bd89be6064680f5112f7e889b85156ab140ef..0e72a3f65107c17b7de1d857bbcbe668fba5f93f 100644 (file)
--- a/R/chronopl.R
+++ b/R/chronopl.R
@@ -1,15 +1,16 @@
-## chronopl.R (2008-11-04)
+## chronopl.R (2011-07-04)
 
 ##   Molecular Dating With Penalized Likelihood
 
-## Copyright 2005-2008 Emmanuel Paradis
+## Copyright 2005-2011 Emmanuel Paradis
 
 ## This file is part of the R-package `ape'.
 ## See the file ../COPYING for licensing issues.
 
-chronopl <- function(phy, lambda, age.min = 1, age.max = NULL,
-                     node = "root", S = 1, tol = 1e-8,
-                     CV = FALSE, eval.max = 500, iter.max = 500, ...)
+chronopl <-
+    function(phy, lambda, age.min = 1, age.max = NULL,
+             node = "root", S = 1, tol = 1e-8,
+             CV = FALSE, eval.max = 500, iter.max = 500, ...)
 {
     n <- length(phy$tip.label)
     ROOT <- n + 1
@@ -79,32 +80,31 @@ chronopl <- function(phy, lambda, age.min = 1, age.max = NULL,
     ratio <- age.min[1]/ini.time[node[1]]
     ini.time <- ini.time*ratio
 
+    ## because if (!is.null(age.max)), 'node' is modified, so we copy it in case CV = TRUE:
+    node.bak <- node
+
     if (length(node) > 1) {
         ini.time[node] <- age.min
         real.edge.length <- ini.time[e[, 1]] - ini.time[e[, 2]]
         while (any(real.edge.length <= 0)) {
             for (i in EDGES) {
-                if (real.edge.length[i] <= 0) {
-                    if (e[i, 1] %in% node) {
-                        ini.time[e[i, 2]] <-
-                            ini.time[e[, 2]] - 2*real.edge.length[i]
-                        next
-                    }
-                    if (e[i, 2] %in% node) {
-                        ini.time[e[i, 1]] <-
-                            ini.time[e[, 1]] + 2*real.edge.length[i]
-                        next
-                    }
-                    ini.time[e[i, 2]] <-
-                        ini.time[e[, 2]] - real.edge.length[i]
-                    ini.time[e[i, 1]] <-
-                        ini.time[e[, 1]] + real.edge.length[i]
+                if (real.edge.length[i] > 0) next
+                if (e[i, 1] %in% node) {
+                    ini.time[e[i, 2]] <- ini.time[e[1, 2]] - 2 * real.edge.length[i]
+                    next
+                }
+                if (e[i, 2] %in% node) {
+                    ini.time[e[i, 1]] <- ini.time[e[1, 1]] + 2 * real.edge.length[i]
+                    next
                 }
+                ##browser()
+                ini.time[e[i, 2]] <- ini.time[e[1, 2]] - real.edge.length[i]
+                ini.time[e[i, 1]] <- ini.time[e[1, 1]] + real.edge.length[i]
             }
             real.edge.length <- ini.time[e[, 1]] - ini.time[e[, 2]]
+            ##print(min(real.edge.length))
         }
     }
-
     ## `unknown.ages' will contain the index of the nodes of unknown age:
     unknown.ages <- n + 1:m
 
@@ -201,7 +201,7 @@ chronopl <- function(phy, lambda, age.min = 1, age.max = NULL,
     if (CV) ophy <- phy
     phy$edge.length <- age[e[, 1]] - age[e[, 2]]
     if (CV) attr(phy, "D2") <-
-        chronopl.cv(ophy, lambda, age.min, age.max, node,
+        chronopl.cv(ophy, lambda, age.min, age.max, node.bak,
                     n, S, tol, eval.max, iter.max, ...)
     phy
 }
@@ -217,9 +217,8 @@ chronopl.cv <- function(ophy, lambda, age.min, age.max, nodes,
     BT <- branching.times(ophy)
     D2 <- numeric(n)
 
-    cat("  dropping tip")
     for (i in 1:n) {
-        cat(" ", i, sep = "")
+        cat("\r  dropping tip ", i, " / ", n, sep = "")
         tr <- drop.tip(ophy, i)
         j <- which(ophy$edge[, 2] == i)
         if (ophy$edge[j, 1] %in% nodes) {