-## chronopl.R (2008-03-26)
+## 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
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
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
}
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) {
if (length(node)) {
chr <- chronopl(tr, lambda, age.min, age.max, node,
S, tol, FALSE, eval.max, iter.max, ...)
+ tmp <-
+ if (Nnode(chr) == Nnode(ophy)) BT else BT[-(ophy$edge[j, 1] - n)]
D2[i] <- sum((tmp - branching.times(chr))^2 / tmp)
} else D2[i] <- 0
}