NNIDist(): best_lower is not a valid lower bound #220

NNIDist() can return a best_lower that is larger than the true NNI distance.

Reprex (TreeDist 2.15.0)

library(TreeDist)
a <- ape::read.tree(text = "(t1,(t2,(t4,t5)),(t3,t6));")
b <- ape::read.tree(text = "(t6,(t1,(t4,t5)),(t3,t2));")
NNIDist(a, b)
#>  lower  best_lower tight_upper  best_upper loose_upper  fack_upper    li_upper
#>      2           3           3           3           5           5           8

The true NNI distance is 2. An explicit path: a → (t2,(t1,(t4,t5)),(t3,t6)) (swap t1↔t2) → b (swap t2↔t6). An exhaustive BFS over the NNI graph also gives 2.

Suspected cause

In update_score(), each region of unmatched edges that is small enough has its exact diameter added to closest_lower_bound, which is returned as best_lower (L159-L160, L427). The diameter is the maximum distance between two trees with that many leaves, so it is an upper bound, not a lower one. Here, one region has 2 unmatched edges (5 leaves), and exact_diameter[5] = 3, which gives best_lower = 3. For regions too large for an exact value, min_diameter is added instead, so best_lower mixes an upper and a lower quantity. The existing test encodes this behaviour (test-tree_distance_nni.R L81: 5 unmatched edges, best_lower = 10L).

Frequency

• Hyolith dataset (reported separately, not re-run here): best_lower > lower in 99.3% of consecutive tree pairs.
• Brute force, n = 6–8. The exact distance came from BFS over the NNI graph; random walks of k NNIs from random trees, 60 walks per k:
  • best_lower exceeded the exact distance in 0/60 walks at k = 1, 23/60 at k = 2, 30/60 at k = 3 and 42/60 at k = 4.
  • Over 60 random pairs at n = 6, it exceeded the exact distance in 44.
• Over all 300 cases, lower, tight_upper, best_upper, loose_upper, fack_upper and li_upper were always valid bounds.

Documentation

?NNIDist says a "named vector containing three integers" (R/tree_distance_nni.R L37) and documents only lower, tight_upper and loose_upper. The returned best_lower, best_upper, fack_upper and li_upper are undocumented.

Suggested direction

Make best_lower a genuine lower bound: add lower (not exact) for regions with an exact diameter, or drop the field. Then document all seven fields, and update the L81 expectation. A side note, not tested: min_diameter[k] exceeds k from k = 37. If it is a lower bound on the diameter (cf. fackMin), it is not a lower bound on one pair's distance, so lower may also be too high for very large regions.