+ ind <- matrix(0L, N - Nbasal, 2)
+ ind[, 1] <- EDGES[-basal]
+ ind[, 2] <- match(e1[EDGES[-basal]], e2)
+
+ age <- numeric(n + m)
+
+#############################################################################
+### 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
+ }
+ }
+ }
+
+ real.edge.length <- ini.time[e1] - ini.time[e2]
+
+ if (any(real.edge.length <= 0))
+ stop("some initial branch lengths are zero or negative;
+ maybe you need to adjust the given dates -- see '?chronopl' for details")
+#############################################################################
+
+ ## because if (!is.null(age.max)), 'node' is modified,
+ ## so we copy it in case CV = TRUE:
+ node.bak <- node
+
+ ## `unknown.ages' will contain the index of the nodes of unknown age:
+ unknown.ages <- n + 1:m
+
+ ## define the bounds for the node ages:
+ lower <- rep(tol, length(unknown.ages))
+ upper <- rep(1/tol, length(unknown.ages))
+
+ if (!is.null(age.max)) { # are some nodes known within some intervals?
+ lower[node - n] <- age.min
+ upper[node - n] <- age.max
+ ## find nodes known within an interval:
+ interv <- which(age.min != age.max)
+ ## drop them from the 'node' since they will be estimated:
+ node <- node[-interv]
+ if (length(node)) age[node] <- age.min[-interv] # update 'age'
+ } else age[node] <- age.min
+
+ if (length(node)) {
+ unknown.ages <- unknown.ages[n - node] # 'n - node' is simplification for '-(node - n)'
+ lower <- lower[n - node]
+ upper <- upper[n - node]