Fix MASTSize() under-counting (root row/column of MAST table never filled)

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).

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,

tests/testthat/test-mast.R

@@ -65,15 +65,68 @@ 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)
   ME(3L, 9)
 
   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))
+  }
+})