diff --git a/DESCRIPTION b/DESCRIPTION index 9049f813a..40740784c 100644 --- a/DESCRIPTION +++ b/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", diff --git a/NEWS.md b/NEWS.md index da4bc8bef..b0606256c 100644 --- a/NEWS.md +++ b/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()`. diff --git a/inst/include/TreeTools/SplitList.h b/inst/include/TreeTools/SplitList.h index f0a06bf69..e9217d562 100644 --- a/inst/include/TreeTools/SplitList.h +++ b/inst/include/TreeTools/SplitList.h @@ -2,23 +2,25 @@ #define _TREETOOLS_SPLITLIST_H #include - #include /* for errors */ +#include /* for heap allocation */ +#include /* 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 // C++20 header for std::popcount - - inline int16 count_bits(splitbit x) { - return static_cast(std::popcount(x)); + inline int32 count_bits(splitbit x) { + return static_cast(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(__builtin_popcountll(x)); + inline int32 count_bits(splitbit x) { + return static_cast(__builtin_popcountll(x)); } #elif defined(_MSC_VER) #include - inline int16 count_bits(splitbit x) { - return static_cast(__popcnt64(x)); + inline int32 count_bits(splitbit x) { + return static_cast(__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 heap_in_split; + std::vector heap_data; + std::vector heap_rows; + + public: SplitList(const Rcpp::RawMatrix &x) { - if (double(x.rows()) > double(std::numeric_limits::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::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(x.rows()); + + /* Check limits */ + if (n_rows > static_cast(std::numeric_limits::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(n_splits) * + static_cast(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; }; } diff --git a/man/SplitInformation.Rd b/man/SplitInformation.Rd index 010333f4e..db27e2894 100644 --- a/man/SplitInformation.Rd +++ b/man/SplitInformation.Rd @@ -23,7 +23,7 @@ into partitions of the specified sizes. \description{ Calculate the phylogenetic information content (\emph{sensu} \insertCite{Steel2006;nobrackets}{TreeTools}) of a split, which -reflects the probability that a uniformly selected random tree will contain# +reflects the probability that a uniformly selected random tree will contain the split: a split that is consistent with a smaller number of trees will have a higher information content. } diff --git a/man/figures/Stemwardness.png b/man/figures/Stemwardness.png index 3b49ff717..d2954f608 100644 Binary files a/man/figures/Stemwardness.png and b/man/figures/Stemwardness.png differ diff --git a/src/splits.cpp b/src/splits.cpp index 9cd54f1a6..096dc2d70 100644 --- a/src/splits.cpp +++ b/src/splits.cpp @@ -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) { diff --git a/src/splits_to_tree.cpp b/src/splits_to_tree.cpp index f0049d7b1..cbd3e5aba 100644 --- a/src/splits_to_tree.cpp +++ b/src/splits_to_tree.cpp @@ -7,8 +7,8 @@ using namespace Rcpp; using namespace TreeTools; inline void insert_ancestor(const int16 tip, const int16 *next_node, - std::array& parent, - std::array& 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::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 parent{}; - alignas(64) std::array patriarch{}; - - std::array 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 stack_parent{}; + alignas(64) std::array stack_patriarch{}; + + // Heap allocation for large trees + std::vector heap_parent; + std::vector heap_patriarch; + + // Pointers to active storage + int16* parent; + int16* patriarch; + + if (use_heap) { + const size_t parent_size = static_cast(n_tip) + + static_cast(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 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]; }); diff --git a/tests/testthat/test-Splits.R b/tests/testthat/test-Splits.R index 163077ec9..974a02e3e 100644 --- a/tests/testthat/test-Splits.R +++ b/tests/testthat/test-Splits.R @@ -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()", {