diff --git a/NEWS.md b/NEWS.md index ed8657f4..4db05c58 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,6 +2,11 @@ ## Bug fixes +- Fix `MASTSize()` under-counting on some trees: the root of each tree was + never compared with the subtrees of the other, so a maximum agreement + subtree confined to one side of a root could be missed. Affects both + `rooted = TRUE` and `rooted = FALSE`, and thus `MASTInfo()`. + - C++ code no longer executes the x86 POPCNT instruction on CPUs that lack it, which crashed R with "illegal operation" on some virtual machines (ms609/TreeTools#287). diff --git a/src/mast.cpp b/src/mast.cpp index c26af5ff..8ecb002c 100644 --- a/src/mast.cpp +++ b/src/mast.cpp @@ -118,6 +118,20 @@ int cpp_mast (IntegerMatrix edge1, IntegerMatrix edge2, IntegerVector nTip) { t2_left, t2_right, n_tip, n_all_nodes); } } + // Root of each tree against every non-root node of the other. + // The root recurrence reads M[root1, child of root2] and + // M[child of root1, root2], so these must be filled first. + // edge2 is in postorder, so M[root1, child of node2] precedes M[root1, node2]. + for (int16 j = 0; j != n_edge; j++) { + fill_M(n_tip, edge2(j, 1), M, t1_descendantsof, + t2_descendantsof, t1_left, t1_right, + t2_left, t2_right, n_tip, n_all_nodes); + } + for (int16 i = 0; i != n_edge; i++) { + fill_M(edge1(i, 1), n_tip, M, t1_descendantsof, + t2_descendantsof, t1_left, t1_right, + t2_left, t2_right, n_tip, n_all_nodes); + } // Root node: fill_M(n_tip, n_tip, M, t1_descendantsof, t2_descendantsof, t1_left, t1_right, diff --git a/tests/testthat/test-mast.R b/tests/testthat/test-mast.R index 2e73f23d..6bb21b49 100644 --- a/tests/testthat/test-mast.R +++ b/tests/testthat/test-mast.R @@ -65,11 +65,12 @@ test_that("MAST size calculated correctly on small trees", { nTip = 5L)) } - ME(2L, 1) + # Values checked against brute-force rooted MAST + ME(3L, 1) ME(3L, 2) ME(4L, 3) - ME(2L, 4) - ME(2L, 5) + ME(3L, 4) + ME(3L, 5) ME(4L, 6) ME(4L, 7) ME(3L, 8) @@ -77,3 +78,55 @@ test_that("MAST size calculated correctly on small trees", { expect_equal(4L, MASTSize(t1, t2, rooted = FALSE)) }) + +test_that("MASTSize() matches brute force", { + # Reported under-count: dropping t1 makes the trees identical + tree1 <- ape::read.tree(text = "(t1,(t2,(t3,t4)),((t5,t6),(t7,t8)));") + tree2 <- ape::read.tree(text = "(t1,(t5,t6),((t2,(t3,t4)),(t7,t8)));") + expect_equal(MASTSize(tree1, tree2, rooted = FALSE), 7L) + + # Rooted: MAST lies within one child of the root of tree2 + rooted1 <- ape::read.tree(text = "((a,b),(c,(d,e)));") + rooted2 <- ape::read.tree(text = "(e,((a,b),(c,d)));") + expect_equal(MASTSize(rooted1, rooted2), 4L) + expect_equal(MASTSize(rooted2, rooted1), 4L) + + AddRootLeaf <- function(tree) { + ape::read.tree(text = paste0("(ROOT,", sub(";$", "", ape::write.tree(tree)), + ");")) + } + BruteMAST <- function(tree1, tree2, rooted) { + labels <- tree1[["tip.label"]] + for (k in rev(seq_along(labels))[-length(labels)]) { + for (keep in combn(labels, k, simplify = FALSE)) { + a <- KeepTip(tree1, keep) + b <- KeepTip(tree2, keep) + if (rooted) { + a <- AddRootLeaf(a) + b <- AddRootLeaf(b) + } else if (k < 4) { + return(k) + } + if (RobinsonFoulds(a, b) == 0) return(k) + } + } + 1L + } + RandomPair <- function(nTip) { + labels <- paste0("t", seq_len(nTip)) + Tidy <- function(tr) { + ape::read.tree(text = ape::write.tree(RenumberTips(tr, labels))) + } + list(Tidy(RandomTree(labels, root = TRUE)), + Tidy(RandomTree(labels, root = TRUE))) + } + + set.seed(4110) + for (i in 1:12) { + pair <- RandomPair(sample(5:7, 1)) + expect_equal(MASTSize(pair[[1]], pair[[2]]), + BruteMAST(pair[[1]], pair[[2]], rooted = TRUE)) + expect_equal(MASTSize(pair[[1]], pair[[2]], rooted = FALSE), + BruteMAST(pair[[1]], pair[[2]], rooted = FALSE)) + } +})