Modernize C++ in NNI distance

.Rbuildignore

@@ -31,5 +31,6 @@ revdep
 ^\.travis\.yml$
 ^\.zenodo\.json$
 ^\.covrignore$
+^\.vs.*$
 ^\.vscode.*$
 ^CRAN-SUBMISSION$

.github/workflows/benchmark.yml

@@ -2,7 +2,11 @@ name: Benchmark
 
 on:
   workflow_dispatch:
-  pull_request:
+  workflow_dispatch:
+  workflow_run:
+    workflows: ["R-CMD-check"]
+    types:
+      - completed
     paths:
       - "src/**"
       - "R/**"

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: TreeDist
 Type: Package
 Title: Calculate and Map Distances Between Phylogenetic Trees
-Version: 2.10.1.9002
+Version: 2.10.1.9003
 Authors@R: c(person("Martin R.", "Smith",
                     email = "martin.smith@durham.ac.uk", 
                     role = c("aut", "cre", "cph", "prg"), 

NEWS.md

@@ -1,7 +1,8 @@
-# TreeDist 2.10.1.9002 (development)
+# TreeDist 2.10.1.9003 (development)
 
 - `HierarchicalMutualInformation()` calculates the information shared between
   pairs of hierarchical partition structures <doi:10.1103/PhysRevE.92.062825>.
+- Support larger trees in NNI distance calculations.
 - Fix crash in `robinson_foulds_all_pairs()` and `RobinsonFoulds(list)`.
 - Fix bug in calculation of `MutualClusteringInfo()`: greedy optimization
   was not guaranteed to find globally optimal matching, causing distances to be

cran-comments.md

@@ -17,12 +17,10 @@
 ## Downstream dependencies
 
 ["revdepcheck"](https://github.com/ms609/TreeDist/actions/workflows/revdepcheck.yml)
-confirmed no changes to worse in the three downstream dependencies:
+confirmed no changes to worse in the four downstream dependencies:
+
+  'Rogue', 'rwty', 'TBRDist', and 'TreeSearch'
 
-  'Rogue', 'TBRDist', and 'TreeSearch'
-
-(all of which I maintain).
-
 
 ## R CMD check results
 There were no ERRORs or WARNINGs.

src/RcppExports.cpp

@@ -175,14 +175,14 @@ BEGIN_RCPP
 END_RCPP
 }
 // cpp_nni_distance
-IntegerVector cpp_nni_distance(const IntegerMatrix edge1, const IntegerMatrix edge2, const IntegerVector nTip);
+IntegerVector cpp_nni_distance(const IntegerMatrix& edge1, const IntegerMatrix& edge2, const IntegerVector& nTip);
 RcppExport SEXP _TreeDist_cpp_nni_distance(SEXP edge1SEXP, SEXP edge2SEXP, SEXP nTipSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
-    Rcpp::traits::input_parameter< const IntegerMatrix >::type edge1(edge1SEXP);
-    Rcpp::traits::input_parameter< const IntegerMatrix >::type edge2(edge2SEXP);
-    Rcpp::traits::input_parameter< const IntegerVector >::type nTip(nTipSEXP);
+    Rcpp::traits::input_parameter< const IntegerMatrix& >::type edge1(edge1SEXP);
+    Rcpp::traits::input_parameter< const IntegerMatrix& >::type edge2(edge2SEXP);
+    Rcpp::traits::input_parameter< const IntegerVector& >::type nTip(nTipSEXP);
     rcpp_result_gen = Rcpp::wrap(cpp_nni_distance(edge1, edge2, nTip));
     return rcpp_result_gen;
 END_RCPP

src/ints.h

@@ -5,11 +5,7 @@
 
 using int16 = int_fast16_t;
 using int32 = int_fast32_t;
-using uint16 = uint_fast16_t;
-using grf_match = std::vector<int> ;
-
-constexpr uint16 INT_16_MAX = INT_FAST16_MAX;
-constexpr uint16 UINT_16_MAX = UINT_FAST16_MAX;
-constexpr int16 NA_INT16 = -0x7fff;
+using uint32 = uint_fast32_t;
+using grf_match = std::vector<int>;
 
 #endif

src/tree_distances.h

@@ -388,7 +388,7 @@ namespace TreeDist {
 
   // Returns lg2_unrooted[x] - lg2_trees_matching_split(y, x - y)
   [[nodiscard]] inline double mmsi_pair_score(const int16 x, const int16 y) noexcept {
-    assert(SL_MAX_TIPS + 2 <= INT_16_MAX); // verify int16 ok
+    assert(SL_MAX_TIPS + 2 <= std::numeric_limits<int16>::max()); // verify int16 ok
 
     return lg2_unrooted[x] - (lg2_rooted[y] + lg2_rooted[x - y]);
   }
@@ -417,7 +417,7 @@ namespace TreeDist {
   }
 
 [[nodiscard]] inline double one_overlap(const int16 a, const int16 b, const int16 n) noexcept {
-    assert(SL_MAX_TIPS + 2 <= INT_16_MAX); // verify int16 ok
+    assert(SL_MAX_TIPS + 2 <= std::numeric_limits<int16>::max()); // verify int16 ok
     if (a == b) {
       return lg2_rooted[a] + lg2_rooted[n - a];
     } else if (a < b) {
@@ -428,7 +428,7 @@ namespace TreeDist {
   }
 
   [[nodiscard]] inline double one_overlap_notb(const int16 a, const int16 n_minus_b, const int16 n) noexcept {
-    assert(SL_MAX_TIPS + 2 <= INT_16_MAX); // verify int16 ok
+    assert(SL_MAX_TIPS + 2 <= std::numeric_limits<int16>::max()); // verify int16 ok
     const int16 b = n - n_minus_b;
     if (a == b) {
       return lg2_rooted[b] + lg2_rooted[n_minus_b];

tests/testthat/test-tree_distance_nni.R

@@ -17,7 +17,18 @@ test_that("NNIDist() handles exceptions", {
 })
 
 test_that("NNIDist() at NNI_MAX_TIPS", {
-  maxTips <- 2048
+  maxTips <- 32768
+  more <- maxTips + 1
+  expect_error(.NNIDistSingle(PectinateTree(more), BalancedTree(more), more),
+               "so many tips")
+  goingQuickly <- TRUE
+  skip_if(goingQuickly)
+
+  heapTips <- 16384 + 1
+  skip_if_not_installed("testthat", "3.2.2")
+  expect_no_error(.NNIDistSingle(PectinateTree(heapTips),
+                                 BalancedTree(heapTips), heapTips))
+
   n <- .NNIDistSingle(PectinateTree(maxTips), BalancedTree(maxTips),
                            maxTips)
   expect_gt(n[["best_upper"]], n[["best_lower"]])
@@ -27,9 +38,6 @@ test_that("NNIDist() at NNI_MAX_TIPS", {
   expect_gte(n[["loose_upper"]], n[["best_upper"]])
   expect_gte(n[["fack_upper"]], n[["best_upper"]])
   expect_gte(n[["li_upper"]], n[["best_upper"]])
-  more <- maxTips + 1
-  expect_error(.NNIDistSingle(PectinateTree(more), BalancedTree(more), more),
-               "so many tips")
 })
 
 test_that("Simple NNI approximations", {
@@ -64,7 +72,7 @@ test_that("Simple NNI approximations", {
                      best_upper = 10L, loose_upper = 18L, fack_upper = 18L,
                      li_upper = Li(5))
 
-  Test <- function(expect, tree) {
+  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")
@@ -76,56 +84,56 @@ test_that("Simple NNI approximations", {
       expectation["loose_upper"] <- min(expect[c("fack_upper", "li_upper")])
     }
 
-    expect_equal(expectation, NNIDist(tree1, tree))
+    expect_equal(NNIDist(tree1, tree), expectation)
     for (i in c(2L, 3L, 4L, 6L)) {
       tree1i <- RootOnNode(tree1, i)
       j <- 0
       for (t2 in unique(lapply(1:9, RootOnNode, tree = tree))) {
-        expect_equal(expectation, NNIDist(tree1i, t2))
+        expect_equal(NNIDist(tree1i, t2), expectation)
       }
     }
   }
 
-  expect_equal(allMatched, NNIDist(BalancedTree(2), PectinateTree(2)))
+  expect_equal(NNIDist(BalancedTree(2), PectinateTree(2)), allMatched)
 
-  expect_equal(oneUnmatched, cpp_nni_distance(edge1, edge2, NTip(tree1)))
-  Test(oneUnmatched, PectinateTree(nTip))
+  expect_equal(cpp_nni_distance(edge1, edge2, NTip(tree1)), oneUnmatched)
+  Test(PectinateTree(nTip), oneUnmatched)
 
   # Identical trees
   tree1 <- Postorder(read.tree(text = "(((a, b), (c, d)), ((e, f), (g, h)));"))
   tree2 <- Postorder(read.tree(text = "(((a, b), (d, c)), ((h, g), (f, e)));"))
-  Test(allMatched, tree1)
-  Test(allMatched, tree2)
+  Test(tree1, allMatched)
+  Test(tree2, allMatched)
 
   # Tree names
   output <- NNIDist(list(bal = tree1, pec = tree2), 
                     as.phylo(0:2, tipLabels = letters[1:8]))
   expect_equal(rownames(output), c("bal", "pec"))
 
   # Only root edge is different
-  Test(oneUnmatched, 
-       Postorder(ape::read.tree(text="(((a, b), (e, f)), ((c, d), (g, h)));")))
+  Test(Postorder(ape::read.tree(text="(((a, b), (e, f)), ((c, d), (g, h)));")),
+       oneUnmatched)
 
   # Two separate regions of difference one
-  Test(oneUnmatched * 2, 
-       read.tree(text="((((a, b), c), d), (e, (f, (g, h))));"))
+  Test(read.tree(text="((((a, b), c), d), (e, (f, (g, h))));"),
+       oneUnmatched * 2)
 
   # One region of three unmatched edges
-  Test(c(lower = 3L, tight_upper = 5L, fack_upper = 8L, li_upper = Li(3)),
-       read.tree(text="(((a, e), (c, d)), ((b, f), (g, h)));"))
+  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)))
 
   # One region of four unmatched edges
-  Test(c(lower = 4L, tight_upper = 7L, fack_upper = 14L, li_upper = Li(4)),
-       tree2 <- ape::read.tree(text="(((a, e), (f, d)), ((b, c), (g, h)));"))
+  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)))
 
   # One region of five unmatched edges
-  Test(fiveUnmatched,
-       ape::read.tree(text="(((a, e), (f, d)), ((b, g), (c, h)));"))
+  Test(ape::read.tree(text="(((a, e), (f, d)), ((b, g), (c, h)));"),
+       fiveUnmatched)
 
-  # Trees with different leaves at root
-  tree1 <- PectinateTree(1:8)
-  Test(fiveUnmatched, 
-       ape::read.tree(text = "(3, ((5, 6), (7, (1, (2, (4, 8))))));"))
+  # Trees with different leaves at root.
+  tree1 <- PectinateTree(1:8) # used in Test()
+  Test(ape::read.tree(text = "(3, ((5, 6), (7, (1, (2, (4, 8))))));"),
+       fiveUnmatched)
 
   # Too different for tight upper bound
   expect_true(is.na(NNIDist(BalancedTree(100), 
@@ -164,28 +172,27 @@ test_that("NNI with lists of trees", {
 test_that("NNIDiameter() is sane", {
 
   exacts <- NNIDiameter(3:12)
-  expect_equal(exacts, do.call(rbind, NNIDiameter(lapply(3:12, as.integer))))
+  expect_equal(do.call(rbind, NNIDiameter(lapply(3:12, as.integer))), exacts)
   expect_true(all(exacts[, "min"] <= exacts[, "exact"]))
   expect_true(all(exacts[, "max"] >= exacts[, "exact"]))
   expect_true(is.na(NNIDiameter(13)[, "exact"]))
   expect_true(is.na(NNIDiameter(1)[, "exact"]))
-  expect_equal(c(exact = 10L), NNIDiameter(BalancedTree(8))[, "exact"])
+  expect_equal(NNIDiameter(BalancedTree(8))[, "exact"], c(exact = 10L))
 
   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) n * ceiling(log2(n)) + n - (2 * ceiling(log2(n)))
   n <- 4:8
-  expect_equal(cbind(
+  expect_equal(NNIDiameter(n), cbind(
     liMin = n - 3L,
     fackMin = FackMin(n - 2L),
     min = pmax(n - 3L, FackMin(4:8 - 2L)),
     exact = exacts[n],
     liMax = liMaxes[n],
     fackMax = FackMax(n - 2L),
     max = pmin(liMaxes[n], FackMax(n - 2L))
-  ), NNIDiameter(n))
+  ))
 
-  expect_equal(NNIDiameter(c(6, 6)), NNIDiameter(as.phylo(0:1, 6)))
-
+  expect_equal(NNIDiameter(as.phylo(0:1, 6)), NNIDiameter(c(6, 6)))
 })