NNIDist(): make lower bounds valid; adopt Fack theorem 3

AGENTS.md

@@ -44,7 +44,7 @@ when ready.
 | `spr.cpp` | 7 KB | SPR distance approximation |
 | `spr_lookup.cpp` | — | SPR lookup-table implementation |
 | `nni_distance.cpp` | 16 KB | NNI distance approximations; HybridBuffer allocation |
-| `li_diameters.h` | 30 KB | Precomputed NNI diameter lookup tables |
+| `li_diameters.h` | <1 KB | Exact NNI diameters for trees with ≤ 12 leaves |
 | `information.h` | 6 KB | log₂ / factorial lookup tables (max 8192); cached at startup |
 | `binary_entropy_counts.cpp` | — | Binary entropy calculations |
 | `day_1985.cpp` | 10 KB | Consensus tree information |

NEWS.md

@@ -2,6 +2,17 @@
 
 ## Bug fixes
 
+- `NNIDist()`: `best_lower` was not a lower bound. It added the diameter of
+  each region of unmatched edges, which bounds the distance from above; it
+  could exceed the true distance (#220). `lower` added a lower bound on the
+  diameter, which exceeded the number of unmatched edges in regions of 37 or
+  more edges. Both now count unmatched edges, the only known lower bound, and
+  are equal.
+
+- `NNIDist()` and `NNIDiameter()` use theorem 3 of Fack _et al._ (2002) for
+  `fack_upper` and `fackMax`, tightening `loose_upper` and `best_upper`.
+  All seven fields of `NNIDist()` are now documented.
+
 - 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

R/tree_distance_nni.R

@@ -31,21 +31,46 @@
 #'    \tab 15 \tab 18 \tab 21 \tab ? \cr
 #'  }
 #' 
+#' Upper bounds treat each region separately: transforming each region in turn
+#' is a valid path between the trees, so the sum of the diameters of the regions
+#' bounds the distance from above.
+#' A diameter is the distance between the _most_ different pair of trees with
+#' that many leaves, so it says nothing about how close a particular pair might
+#' be: the only known lower bound is the number of unmatched edges
+#' \insertCite{Li1996}{TreeDist}.
+#' An optimal path need not leave matched edges untouched
+#' \insertCite{Li1996}{TreeDist}, so even the exact distance between the trees
+#' that each region induces is not a lower bound.
+#'
 #' @template tree12Params
 #' 
 #' @return `NNIDist()` returns, for each pair of trees, a named vector
-#'  containing three integers:
-#' 
-#' - `lower` is a lower bound on the NNI distance, and corresponds
-#' to the RF distance between the trees. 
-#' 
+#'  containing seven integers:
+#'
+#' - `lower` is a lower bound on the NNI distance: the number of edges in
+#' one tree that do not occur in the other, i.e. half the Robinson-Foulds
+#' distance between the trees.
+#'
+#' - `best_lower` is the tightest available lower bound.
+#' No bound tighter than `lower` is known, so the two values are equal.
+#'
 #' - `tight_upper` is an upper bound on the distance, based on calculated
 #' maximum diameters for trees with < 13 leaves.  _NA_ is returned if trees are
 #' too different to employ this approach.
-#' 
+#'
+#' - `best_upper` is the tightest available upper bound: it uses the exact
+#' diameter for regions with < 13 leaves, and `loose_upper` for larger
+#' regions.
+#'
 #' - `loose_upper` is a looser upper bound on the distance, using _n_ log _n_ +
-#' O(_n_).
-#' 
+#' O(_n_): the smaller of `fack_upper` and `li_upper`.
+#'
+#' - `fack_upper` is an upper bound based on the maximum diameter given by
+#' theorem 3 of \insertCite{Fack2002;textual}{TreeDist}.
+#'
+#' - `li_upper` is an upper bound based on the sorting approach of
+#' \insertCite{Li1996;textual}{TreeDist}.
+#'
 #' 
 #' @references 
 #' \insertAllCited{}
@@ -108,9 +133,9 @@ NNIDist <- function(tree1, tree2 = tree1) {
 #' - `liMax`: Upper bound on diameter following Li _et al._ (1996), i.e. 
 #'   \deqn{n \log2{n} + \textrm{O}(n)}{n log&#8322; _n_ + O(_n_)};
 #'   
-#' - `fackMax`: Upper bound on diameter following Fack _et al_. (2002), i.e.
-#'   (\deqn{N - 2}{_N_ &minus; 2}) ceiling(\deqn{\log2{n}}{log&#8322; _N_}) 
-#'   + _N_;
+#' - `fackMax`: Upper bound on diameter following theorem 3 of
+#'   Fack _et al_. (2002), i.e. _N_ _c_ &minus; 2^_c_ + 1 + 2(_N_ &minus;
+#'   ceiling(log&#8322; (_N_ + 1))), where _c_ = ceiling(log&#8322; _N_);
 #'   
 #' - `max`: The smaller of `liMax` and `fackMax`;
 #'   
@@ -145,7 +170,8 @@ NNIDiameter.numeric <- function(tree) {
   fackMin <- ceiling(lfactorial(n) / 4L / log(2))
 
   liMax <- .SortingNumber(tree) + (2 * .DegenerateDistance(tree))
-  fackMax <- ((n - 2L) * ceiling(log2(n))) + n
+  lgN <- ceiling(log2(n))
+  fackMax <- (n * lgN) - (2 ^ lgN) + 1 + (2 * (n - ceiling(log2(n + 1))))
   cbind(liMin = liMin,
         fackMin = fackMin,
         min = pmax(liMin, fackMin),

src/nni_distance.cpp

@@ -99,8 +99,14 @@ static void initialize_cache() {
   for (int32_t i = 8192 + 1; i != 16384 + 1; i++) lg2_ceiling[i] = 14;
   for (int32_t i = 16384 + 1; i != 32768 + 1; i++) lg2_ceiling[i] = 15;
 
+  /* Fack et al. (2002), theorem 3: with N = n - 2 internal nodes,
+   * diam <= N ceil(lg N) - 2^ceil(lg N) + 1 + 2 (N - ceil(lg(N + 1))).
+   * Never looser than the bound in their theorem 1. */
   for (int32_t i = 4; i != NNI_MAX_TIPS + 1; i++) {
-    fack_lookup[i] = int32_t(((i - 2 - 2) * lg2_ceiling[i - 2]) + i - 2);
+    const int32_t n_node = i - 2;
+    const int32_t lg_n = lg2_ceiling[n_node];
+    fack_lookup[i] = n_node * lg_n - (int32_t(1) << lg_n) + 1 +
+      2 * (n_node - lg2_ceiling[n_node + 1]);
   }
   for (int32_t i = 4; i != NNI_MAX_TIPS + 1; i++) {
     const int32_t log_ceiling = lg2_ceiling[i];
@@ -138,14 +144,16 @@ inline void ensure_cache() {
 
 // Score subtree, add to score, and reset subtree size
 void update_score(const int32_t subtree_edges, int32_t *lower_bound,
-                  int32_t *closest_lower_bound,
                   int32_t *tight_bound, 
                   int32_t *closest_upper_bound,
                   int32_t *loose_bound, 
                   int32_t *li_bound, int32_t *fack_bound) {
   if (subtree_edges) {
     const int32_t subtree_tips = subtree_edges + 3;
-    const int32_t lower = min_diameter[subtree_edges];
+    // Each unmatched edge must undergo at least one NNI (Li et al. 1996).
+    // Diameters bound a single pair's distance from above only, so they
+    // cannot tighten this lower bound.
+    const int32_t lower = subtree_edges;
     const int32_t this_li = li[subtree_tips],
       fack = fack_lookup[subtree_tips],
       upper = this_li < fack ? this_li : fack;
@@ -157,13 +165,11 @@ void update_score(const int32_t subtree_edges, int32_t *lower_bound,
 
     if (subtree_tips <= N_EXACT) {
       const int32_t exact = exact_diameter[subtree_tips];
-      *closest_lower_bound += exact;
       if (*tight_bound != NA_INT32) {
         *tight_bound += exact;
       }
       *closest_upper_bound += exact;
     } else {
-      *closest_lower_bound += lower;
       *tight_bound = NA_INT32;
       *closest_upper_bound += upper;
     }
@@ -296,7 +302,6 @@ IntegerVector cpp_nni_distance(const IntegerMatrix& edge1,
   const int32_t n_edge = int32_t(edge1.nrow());
 
   int32_t lower_bound = 0;
-  int32_t best_lower_bound = 0;
   int32_t tight_score_bound = 0;
 
   int32_t loose_score_bound = 0;
@@ -379,8 +384,7 @@ IntegerVector cpp_nni_distance(const IntegerMatrix& edge1,
         unmatched_below[parent_i - node_0] += unmatched_below[child_i - node_0];
       } else {
         update_score(unmatched_below[child_i - node_0], &lower_bound,
-                     &best_lower_bound, &tight_score_bound,
-                     &best_upper_bound, &loose_score_bound,
+                     &tight_score_bound, &best_upper_bound, &loose_score_bound,
                      &li_score_bound, &fack_score_bound);
       }
     }
@@ -400,31 +404,31 @@ IntegerVector cpp_nni_distance(const IntegerMatrix& edge1,
       if (!matched_1[root_child_2 - node_0]) {
         update_score(int32_t(unmatched_below[root_node]
                              + unmatched_1
-                             + unmatched_2), &lower_bound, &best_lower_bound,
+                             + unmatched_2), &lower_bound,
                      &tight_score_bound, &best_upper_bound, &loose_score_bound,
                      &li_score_bound, &fack_score_bound);
       } else {
-        update_score(unmatched_1, &lower_bound, &best_lower_bound, 
+        update_score(unmatched_1, &lower_bound,
                      &tight_score_bound, &best_upper_bound, &loose_score_bound,
                      &li_score_bound, &fack_score_bound);
-        update_score(unmatched_2, &lower_bound, &best_lower_bound, 
+        update_score(unmatched_2, &lower_bound,
                      &tight_score_bound, &best_upper_bound, &loose_score_bound,
                      &li_score_bound, &fack_score_bound);
       }
     } else {
-      update_score(unmatched_1, &lower_bound, &best_lower_bound, 
-                   &tight_score_bound, &best_upper_bound, &loose_score_bound, 
+      update_score(unmatched_1, &lower_bound,
+                   &tight_score_bound, &best_upper_bound, &loose_score_bound,
                    &li_score_bound, &fack_score_bound);
     }
   } else {
-    update_score(unmatched_below[root_node], &lower_bound, &best_lower_bound,
+    update_score(unmatched_below[root_node], &lower_bound,
                  &tight_score_bound, &best_upper_bound, &loose_score_bound,
                  &li_score_bound, &fack_score_bound);
   }
 
   return IntegerVector::create(
     Named("lower")   = lower_bound,
-    _["best_lower"] = best_lower_bound,
+    _["best_lower"] = lower_bound, // No tighter lower bound is known
     _["tight_upper"] = tight_score_bound == NA_INT32 
                      ? NA_INTEGER : tight_score_bound,
     _["best_upper"] = best_upper_bound,

tests/testthat/test-tree_distance_nni.R

@@ -57,7 +57,12 @@ test_that("Simple NNI approximations", {
   edge1 <- Postorder(tree1$edge)
   edge2 <- Postorder(tree2$edge)
 
-  Fack <- function(n) ((n - 2) * ceiling(log2(n))) + n
+  # Fack et al. 2002, theorem 3; N = internal nodes = n_edges + 1
+  Fack <- function(n_edges) {
+    N <- n_edges + 1
+    lgN <- ceiling(log2(N))
+    N * lgN - 2 ^ lgN + 1 + 2 * (N - ceiling(log2(N + 1)))
+  }
 
   Sorting <- function(n) {
     lc <- ceiling(log2(n))
@@ -76,19 +81,22 @@ test_that("Simple NNI approximations", {
                   best_upper = 0L, loose_upper = 0L, fack_upper = 0L,
                   li_upper = 0L)
   oneUnmatched <- c(lower = 1L, best_lower = 1L, tight_upper = 1L,
-                    best_upper = 1L, loose_upper = 2L, fack_upper = 2L,
+                    best_upper = 1L, loose_upper = 1L, fack_upper = Fack(1),
                     li_upper = Li(1))
-  fiveUnmatched <- c(lower = 5L, best_lower = 10L, tight_upper = 10L,
-                     best_upper = 10L, loose_upper = 18L, fack_upper = 18L,
+  fiveUnmatched <- c(lower = 5L, best_lower = 5L, tight_upper = 10L,
+                     best_upper = 10L, loose_upper = 17L, fack_upper = Fack(5),
                      li_upper = Li(5))
 
   Test <- function(tree, expect) {
     expectation <- rep(NA_integer_, 7L)
     names(expectation) <- c("lower", "best_lower", "tight_upper", "best_upper",
                             "loose_upper", "fack_upper", "li_upper")
     expectation[names(expect)] <- expect
-    if (is.na(expectation["best_lower"]) && !is.na(expect["tight_upper"])) {
-      expectation[c("best_lower", "best_upper")] <- expect["tight_upper"]
+    if (is.na(expectation["best_lower"])) {
+      expectation["best_lower"] <- expect["lower"]
+    }
+    if (is.na(expectation["best_upper"]) && !is.na(expect["tight_upper"])) {
+      expectation["best_upper"] <- expect["tight_upper"]
     }
     if (is.na(expect["loose_upper"])) {
       expectation["loose_upper"] <- min(expect[c("fack_upper", "li_upper")])
@@ -130,11 +138,11 @@ test_that("Simple NNI approximations", {
 
   # One region of three unmatched edges
   Test(read.tree(text="(((a, e), (c, d)), ((b, f), (g, h)));"),
-       c(lower = 3L, tight_upper = 5L, fack_upper = 8L, li_upper = Li(3)))
+       c(lower = 3L, tight_upper = 5L, fack_upper = Fack(3), li_upper = Li(3)))
 
   # One region of four unmatched edges
   Test(ape::read.tree(text="(((a, e), (f, d)), ((b, c), (g, h)));"),
-       c(lower = 4L, tight_upper = 7L, fack_upper = 14L, li_upper = Li(4)))
+       c(lower = 4L, tight_upper = 7L, fack_upper = Fack(4), li_upper = Li(4)))
 
   # One region of five unmatched edges
   Test(ape::read.tree(text="(((a, e), (f, d)), ((b, g), (c, h)));"),
@@ -154,6 +162,57 @@ test_that("Simple NNI approximations", {
 
 })
 
+test_that("NNIDist() lower bounds are valid", {
+  # Exact NNI distance is 2 (ms609/TreeDist issue 220)
+  a <- read.tree(text = "(t1,(t2,(t4,t5)),(t3,t6));")
+  b <- read.tree(text = "(t6,(t1,(t4,t5)),(t3,t2));")
+  expect_equal(NNIDist(a, b)[c("lower", "best_lower")],
+               c(lower = 2L, best_lower = 2L))
+})
+
+.ExpectBoundsBracketExact <- function(nTip) {
+  # Exact distances by breadth-first search over all trees, which are an NNI
+  # apart iff they differ in a single split
+  trees <- as.phylo(seq_len(NUnrooted(nTip)) - 1, nTip)
+  adjacent <- as.matrix(RobinsonFoulds(trees)) == 2
+  exact <- matrix(NA_integer_, length(trees), length(trees))
+  for (origin in seq_along(trees)) {
+    dist <- rep(NA_integer_, length(trees))
+    dist[origin] <- 0L
+    frontier <- origin
+    step <- 0L
+    while (length(frontier)) {
+      step <- step + 1L
+      reached <- which(colSums(adjacent[frontier, , drop = FALSE]) > 0 &
+                         is.na(dist))
+      dist[reached] <- step
+      frontier <- reached
+    }
+    exact[origin, ] <- dist
+  }
+  testthat::expect_equal(max(exact), NNIDiameter(nTip)[, "exact"],
+                         ignore_attr = TRUE)
+
+  bounds <- lapply(NNIDist(trees), as.matrix)
+  exact <- exact[lower.tri(exact)]
+  .Bound <- function(field) bounds[[field]][lower.tri(bounds[[field]])]
+  testthat::expect_true(all(.Bound("lower") <= exact))
+  testthat::expect_equal(.Bound("best_lower"), .Bound("lower"))
+  for (upper in c("tight_upper", "best_upper", "loose_upper", "fack_upper",
+                  "li_upper")) {
+    testthat::expect_true(all(.Bound(upper) >= exact))
+  }
+}
+
+test_that("NNIDist() bounds bracket exact distances: 6 leaves", {
+  .ExpectBoundsBracketExact(6)
+})
+
+test_that("NNIDist() bounds bracket exact distances: 7 leaves", {
+  skip_on_cran() # c. 70 s
+  .ExpectBoundsBracketExact(7)
+})
+
 test_that("NNI with lists of trees", {
   tree1 <- BalancedTree(1:8)
   list1 <- list(tree1, PectinateTree(as.character(1:8)),
@@ -192,7 +251,10 @@ test_that("NNIDiameter() is sane", {
   FackMin <- function(n) ceiling(0.25 * lfactorial(n) / log(2))
   exacts <- c(0, 0, 0, 1, 3, 5, 7, 10, 12, 15, 18, 21)
   liMaxes <- c(0, 1, 3, 5, 8, 13, 16, 21, 25, 31, 37, 43, 47, 53, 59, 65)
-  FackMax <- function(n) n * ceiling(log2(n)) + n - (2 * ceiling(log2(n)))
+  FackMax <- function(n) {
+    lgN <- ceiling(log2(n))
+    n * lgN - 2 ^ lgN + 1 + 2 * (n - ceiling(log2(n + 1)))
+  }
   n <- 4:8
   expect_equal(NNIDiameter(n), cbind(
     liMin = n - 3L,