Large tree support: SplitList (#241) & splits_to_edge() (#242)

Author: ms609

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: TreeTools
 Title: Create, Modify and Analyse Phylogenetic Trees
-Version: 2.0.0.9000
+Version: 2.0.0.9001
 Authors@R: c(
   person("Martin R.", 'Smith', role = c("aut", "cre", "cph"), 
          email = "martin.smith@durham.ac.uk",

NEWS.md

@@ -1,3 +1,6 @@
+# TreeTools 2.0.0.9001 (development) #
+- Remove hard limit on tree size in `SplitList`.
+
 # TreeTools 2.0.0.9000 (development) #
 - `MatrixToPhyDat()` gains `tipLabels` parameter.
 - Document return value for `J1Index()`.

inst/include/TreeTools/SplitList.h

@@ -2,23 +2,25 @@
 #define _TREETOOLS_SPLITLIST_H
 
 #include <Rcpp/Lightest>
-
 #include <stdexcept> /* for errors */
+#include <vector>    /* for heap allocation */
+#include <algorithm> /* for std::fill */
 
 #include "assert.h" /* for ASSERT */
-#include "types.h" /* for int16 */
+#include "types.h" /* for int16, int32 */
 
 using splitbit = uint_fast64_t;
 
 #define R_BIN_SIZE int16(8)
 #define SL_BIN_SIZE int16(64)
 #define SL_MAX_BINS int16(32)
-/* 64*32 is about the largest size for which two SplitList objects reliably fit
- * on the stack (as required in TreeDist; supporting more leaves would mean
- * refactoring to run on the heap (and, trivially, converting int16 to int32
- * for split*bin implicit calculation in state[split][bin]?) */
-#define SL_MAX_TIPS (SL_BIN_SIZE * SL_MAX_BINS) // 32 * 64 = 2048
-#define SL_MAX_SPLITS (SL_MAX_TIPS - 3) /* no slower than a power of two */
+
+/* * Stack allocation limits (Legacy support for speed)
+ * Trees smaller than this will use stack arrays.
+ * Trees larger will trigger heap allocation.
+ */
+#define SL_MAX_TIPS (SL_BIN_SIZE * SL_MAX_BINS) // 2048
+#define SL_MAX_SPLITS (SL_MAX_TIPS - 3) 
 
 #define INLASTBIN(n, size) int16((size) - int16((size) - int16((n) % (size))) % (size))
 #define INSUBBIN(bin, offset)                                  \
@@ -38,31 +40,29 @@ namespace TreeTools {
 
 #if __cplusplus >= 202002L
 #include <bit> // C++20 header for std::popcount
-  
-  inline int16 count_bits(splitbit x) {
-    return static_cast<int16>(std::popcount(x));
+  inline int32 count_bits(splitbit x) {
+    return static_cast<int32>(std::popcount(x));
   }
-  
   // Option 2: Fallback for C++17 and older
 #else
 #if defined(__GNUC__) || defined(__clang__)
   // GCC and Clang support __builtin_popcountll for long long
-  inline int16 count_bits(splitbit x) {
-    return static_cast<int16>(__builtin_popcountll(x));
+  inline int32 count_bits(splitbit x) {
+    return static_cast<int32>(__builtin_popcountll(x));
   }
 #elif defined(_MSC_VER)
 #include <intrin.h>
-  inline int16 count_bits(splitbit x) {
-    return static_cast<int16>(__popcnt64(x));
+  inline int32 count_bits(splitbit x) {
+    return static_cast<int32>(__popcnt64(x));
   }
 #else
   // A slower, but safe and highly portable fallback for all other compilers
   // This is a last resort if no built-in is available.
-  inline int16_t count_bits(splitbit x) {
-    int16_t count = 0;
+  inline int32_t count_bits(splitbit x) {
+    int32_t count = 0;
     while (x != 0) {
       x &= (x - 1);
-      count++;
+      ++count;
     }
     return count;
   }
@@ -72,45 +72,73 @@ namespace TreeTools {
 
   class SplitList {
   public:
-    int16 n_splits, n_bins;
-    int16 in_split[SL_MAX_SPLITS];
-    splitbit state[SL_MAX_SPLITS][SL_MAX_BINS];
+    int32 n_splits;
+    int32 n_bins;
+    int32* in_split;
+    splitbit** state;
+  
+  private:
+    /* STACK STORAGE (Fast path for small trees) */
+    int32 stack_in_split[SL_MAX_SPLITS];
+    splitbit stack_state[SL_MAX_SPLITS][SL_MAX_BINS];
+    splitbit* stack_rows[SL_MAX_SPLITS];
+    
+    /* HEAP STORAGE (Large trees) */
+    std::vector<int32> heap_in_split;
+    std::vector<splitbit> heap_data;      
+    std::vector<splitbit*> heap_rows;     
+    
+  public:
     SplitList(const Rcpp::RawMatrix &x) {
-      if (double(x.rows()) > double(std::numeric_limits<int16>::max())) {
-        Rcpp::stop("This many splits cannot be supported. "
-                   "Please contact the TreeTools maintainer if "
-                   "you need to use more!");
-      }
-      if (double(x.cols()) > double(std::numeric_limits<int16>::max())) {
-        Rcpp::stop("This many leaves cannot be supported. "
-                     "Please contact the TreeTools maintainer if "
-                     "you need to use more!");
+      
+      const double n_rows = static_cast<double>(x.rows());
+      
+      /* Check limits */
+      if (n_rows > static_cast<double>(std::numeric_limits<int32>::max())) {
+        Rcpp::stop("Too many splits (exceeds int32 limit).");                   // #nocov
       }
 
-      n_splits = int16(x.rows());
+      n_splits = int32(x.rows());
       ASSERT(n_splits >= 0);
 
-      const int16 n_input_bins = int16(x.cols());
-
-      n_bins = int16(n_input_bins + R_BIN_SIZE - 1) / input_bins_per_bin;
-
-      if (n_bins > SL_MAX_BINS) {
-        Rcpp::stop("This many leaves cannot be supported. "
-                   "Please contact the TreeTools maintainer if "
-                   "you need to use more!");
-      }
+      const int32 n_input_bins = int32(x.cols());
+      ASSERT(n_input_bins > 0);
+      n_bins = int32(n_input_bins + R_BIN_SIZE - 1) / input_bins_per_bin;
+      
+      bool use_heap = (n_splits > SL_MAX_SPLITS) || (n_bins > SL_MAX_BINS);
 
-      for (int16 split = 0; split < n_splits; ++split) {
-        in_split[split] = 0;
+      if (use_heap) {
+        heap_in_split.resize(n_splits, 0);
+        in_split = heap_in_split.data();
+        
+        size_t total_elements = static_cast<size_t>(n_splits) *
+          static_cast<size_t>(n_bins);
+        heap_data.resize(total_elements); 
+        
+        heap_rows.resize(n_splits);
+        for (int32 i = 0; i < n_splits; ++i) {
+          heap_rows[i] = &heap_data[i * n_bins];
+        }
+        state = heap_rows.data();
+        
+      } else {
+        in_split = stack_in_split;
+        
+        for (int32 i = 0; i < n_splits; ++i) {
+          stack_rows[i] = stack_state[i];
+          in_split[i] = 0;
+        }
+        state = stack_rows;
       }
 
-      for (int16 bin = 0; bin < n_bins - 1; ++bin) {
-        const int16 bin_offset = bin * input_bins_per_bin;
+      
+      for (int32 bin = 0; bin < n_bins - 1; ++bin) {
+        const int32 bin_offset = bin * input_bins_per_bin;
 
-        for (int16 split = 0; split < n_splits; ++split) {
+        for (int32 split = 0; split < n_splits; ++split) {
           splitbit combined = splitbit(x(split, bin_offset));
 
-          for (int16 input_bin = 1; input_bin < input_bins_per_bin; ++input_bin) {
+          for (int32 input_bin = 1; input_bin < input_bins_per_bin; ++input_bin) {
             combined |= splitbit(x(split, bin_offset + input_bin)) <<
               (R_BIN_SIZE * input_bin);
           }
@@ -120,19 +148,26 @@ namespace TreeTools {
         }
       }
 
-      const int16 last_bin = n_bins - 1;
-      const int16 raggedy_bins = INLASTBIN(n_input_bins, R_BIN_SIZE);
+      const int32 last_bin = n_bins - 1;
+      const int32 raggedy_bins = INLASTBIN(n_input_bins, R_BIN_SIZE);
 
-      for (int16 split = 0; split < n_splits; ++split) {
+      for (int32 split = 0; split < n_splits; ++split) {
         state[split][last_bin] = INSUBBIN(last_bin, 0);
 
-        for (int16 input_bin = 1; input_bin < raggedy_bins; ++input_bin) {
+        for (int32 input_bin = 1; input_bin < raggedy_bins; ++input_bin) {
           state[split][last_bin] += INBIN(input_bin, last_bin);
         }
 
         in_split[split] += count_bits(state[split][last_bin]);
       }
     }
+    
+    // Default destructor handles vector cleanup automatically
+    ~SplitList() = default;
+    
+    // Disable copy/move to prevent pointer invalidation issues
+    SplitList(const SplitList&) = delete;
+    SplitList& operator=(const SplitList&) = delete;
   };
 }
 

man/SplitInformation.Rd

@@ -43,7 +43,7 @@ Rcpp::RawMatrix cpp_edge_to_splits(const Rcpp::IntegerMatrix& edge,
   const uintx n_edge = edge.rows();
   if (n_edge + 1 >= NOT_TRIVIAL) {
     Rcpp::stop("Too many edges in tree for edge_to_splits: "                    // # nocov
-                 "Contact maintainer for advice");                   // # nocov
+               "Contact maintainer for advice");                                // # nocov
   }
 
   if (nTip[0] < 1) {

man/figures/Stemwardness.png

@@ -7,8 +7,8 @@ using namespace Rcpp;
 using namespace TreeTools;
 
 inline void insert_ancestor(const int16 tip, const int16 *next_node,
-                            std::array<int16, SL_MAX_TIPS + SL_MAX_SPLITS>& parent,
-                            std::array<int16, SL_MAX_TIPS>& patriarch) {
+                            int16* parent,
+                            int16* patriarch) {
   if (patriarch[tip]) {
     parent[patriarch[tip]] = *next_node;
   } else {
@@ -19,9 +19,6 @@ inline void insert_ancestor(const int16 tip, const int16 *next_node,
 
 // [[Rcpp::export]]
 IntegerMatrix splits_to_edge(const RawMatrix splits, const IntegerVector nTip) {
-  if (double(nTip[0]) > double(std::numeric_limits<int16>::max())) {
-    Rcpp::stop("This many tips are not (yet) supported.");
-  }
   const int16 n_tip = int16(nTip[0]);
   if (splits.nrow() == 0) {
     IntegerMatrix ret(n_tip, 2);
@@ -32,12 +29,38 @@ IntegerMatrix splits_to_edge(const RawMatrix splits, const IntegerVector nTip) {
     return ret;
   }
   const SplitList x(splits);
-  alignas(64) std::array<int16, SL_MAX_TIPS + SL_MAX_SPLITS> parent{};
-  alignas(64) std::array<int16, SL_MAX_TIPS> patriarch{};
-
-  std::array<int16, SL_MAX_SPLITS> split_order;
-  std::iota(split_order.begin(), split_order.begin() + x.n_splits, 0);
-  std::sort(split_order.begin(), split_order.begin() + x.n_splits,
+  
+  // Decide whether to use stack or heap allocation based on tree size
+  const bool use_heap = (n_tip > SL_MAX_TIPS) || (x.n_splits > SL_MAX_SPLITS);
+  
+  // Stack allocation for small trees (fast path)
+  alignas(64) std::array<int16, SL_MAX_TIPS + SL_MAX_SPLITS> stack_parent{};
+  alignas(64) std::array<int16, SL_MAX_TIPS> stack_patriarch{};
+  
+  // Heap allocation for large trees
+  std::vector<int16> heap_parent;
+  std::vector<int16> heap_patriarch;
+  
+  // Pointers to active storage
+  int16* parent;
+  int16* patriarch;
+  
+  if (use_heap) {
+    const size_t parent_size = static_cast<size_t>(n_tip) + 
+                               static_cast<size_t>(x.n_splits);
+    heap_parent.resize(parent_size, 0);
+    heap_patriarch.resize(n_tip, 0);
+    parent = heap_parent.data();
+    patriarch = heap_patriarch.data();
+  } else {
+    parent = stack_parent.data();
+    patriarch = stack_patriarch.data();
+  }
+  
+  // Allocate split_order appropriately
+  std::vector<int16> split_order(x.n_splits);
+  std::iota(split_order.begin(), split_order.end(), 0);
+  std::sort(split_order.begin(), split_order.end(),
             [&in_split = x.in_split](int16 a, int16 b) {
               return in_split[a] > in_split[b];
             });

src/splits.cpp

@@ -512,11 +512,15 @@ test_that("Split combination", {
   #TODO: Fully test splits with large (> 8 tip) trees
 })
 
-test_that("as.phylo.Splits() fails gracefully", {
-  expect_error(
-    as.phylo(as.Splits(BalancedTree(3000))),
-    "many leaves cannot be supported"
-  )
+test_that("as.phylo.Splits() supports large trees", {
+  tree3000 <- BalancedTree(3000)
+  expect_no_error(splits3000 <- as.Splits(tree3000))
+  result <- as.phylo(splits3000)
+  
+  # Verify it's a valid phylo object
+  expect_s3_class(result, "phylo")
+  expect_equal(length(result$tip.label), 3000)
+  expect_equal(NTip(result), 3000)
 })
 
 test_that("as.phylo.Splits()", {