- ##ini.time <- node.depth(phy)[-TIPS] - 1
- ini.time <- node.depth(phy) - 1
-
- ## first, rescale all times with respect to
- ## the age of the 1st element of `node':
- ratio <- age.min[1]/ini.time[node[1]]
- ini.time <- ini.time*ratio
-
- 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) 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]
+#############################################################################
+### This bit sets 'ini.time' and should result in no negative branch lengths
+
+ seq.nod <- .Call("seq_root2tip", phy$edge, n, phy$Nnode, PACKAGE = "ape")
+
+ ini.time <- age
+ ini.time[ROOT:(n + m)] <- NA
+ ini.time[node] <- if (is.null(age.max)) age.min else (age.min + age.max) / 2
+
+ ## if no age given for the root, find one approximately:
+ if (is.na(ini.time[ROOT]))
+ ini.time[ROOT] <- if (is.null(age.max)) 3 * max(age.min) else 3 * max(age.max)
+
+ ISnotNA.ALL <- unlist(lapply(seq.nod, function(x) sum(!is.na(ini.time[x]))))
+ o <- order(ISnotNA.ALL, decreasing = TRUE)
+
+ for (y in seq.nod[o]) {
+ ISNA <- is.na(ini.time[y])
+ if (any(ISNA)) {
+ i <- 2L # we know the 1st value is not NA, so we start at the 2nd one
+ while (i <= length(y)) {
+ if (ISNA[i]) { # we stop at the next NA
+ j <- i + 1L
+ while (ISNA[j]) j <- j + 1L # look for the next non-NA
+ nb.val <- j - i
+ by <- (ini.time[y[i - 1L]] - ini.time[y[j]]) / (nb.val + 1)
+ ini.time[y[i:(j - 1L)]] <- ini.time[y[i - 1L]] - by * seq_len(nb.val)
+ i <- j + 1L
+ } else i <- i + 1L