#' Count resolved quartets
#'
#' Counts how many quartets are resolved or unresolved in a given tree,
#' following \insertCite{Brodal2013;textual}{Quartet}.
#'
#' Trees with more than 477 leaves risk encountering integer overflow errors,
#' as the number of quartets is larger than can be stored in R's signed
#' 32-bit integer representation. If warnings are thrown, check subsequent
#' calculations for errors.
#'
#' @template treeParam
#' @param countTriplets Logical; if `TRUE`, the function will return the number
#' of triplets instead of the number of quartets.
#'
#' @return `ResolvedQuartets()` returns a vector of length two, listing the
#' number of quartets (or triplets) that are `[1]` resolved; `[2]` unresolved
#' in the specified tree.
#'
#' @template MRS
#'
#' @family quartet counting functions
#'
#' @examples
#' data(sq_trees)
#'
#' ResolvedTriplets(sq_trees$collapse_some)
#' # Equivalent to:
#' ResolvedQuartets(sq_trees$collapse_some, countTriplets = TRUE)
#'
#' vapply(sq_trees, ResolvedQuartets, integer(2))
#'
#'
#' @references \insertAllCited{}
#'
#' @importFrom TreeTools Preorder
#' @export
ResolvedQuartets <- function (tree, countTriplets = FALSE) {
.CheckSize(tree)
tree <- Preorder(tree)
nTip <- length(tree$tip.label)
nNode <- tree$Nnode
edge <- tree$edge
parent <- edge[, 1]
child <- edge[, 2]
children <- lapply(nTip + seq_len(nNode), function (node) {
child[parent == node]
})
# Algebraic terms follow Brodal et al. (2013)
n <- rep(1, nTip + nNode) # Will be overwritten
unresolvedTripletsRootedHere <-
unresolvedQuartetsRootedHere <- integer(nNode)
for (node in rev(seq_len(nNode)) + nTip) {
nodeChildren <- children[[node - nTip]]
n[node] <- sum(n[nodeChildren])
# s = subtree size
s_vi <- n[nodeChildren[1L]]
# p = pairs; t = triplets; q = quartets
p_vi <- t_vi <- q_vi <- 0L
for (i in seq_along(nodeChildren)[-1]) {
n_vi <- n[nodeChildren[i]]
# Order is important: we need to use previous values in each sum
q_vi <- q_vi + (n_vi * t_vi)
t_vi <- t_vi + (n_vi * p_vi)
p_vi <- p_vi + (n_vi * s_vi)
s_vi <- s_vi + n_vi
}
unresolvedTripletsRootedHere[node - nTip] <- t_vi
unresolvedQuartetsRootedHere[node - nTip] <- q_vi + t_vi * (nTip - s_vi)
}
unresolved <- ifelse(countTriplets, sum(unresolvedTripletsRootedHere),
sum(unresolvedQuartetsRootedHere))
resolved <- choose(nTip, ifelse(countTriplets, 3, 4)) - unresolved
if (any(c(resolved, unresolved) > .Machine$integer.max)) {
stop("Sorry: trees too large for integer representation")
} else if (resolved + unresolved > .Machine$integer.max) {
warning("Large numbers: integer overflow likely")
}
# Return:
as.integer(c(resolved, unresolved))
}
#' @describeIn ResolvedQuartets Convenience function to calculate the number of
#' resolved/unresolved triplets.
#' @export
ResolvedTriplets <- function (tree) ResolvedQuartets(tree = tree,
countTriplets = TRUE)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.