From 3296f9714915f00f70a43f3f3d460e1b408178ec Mon Sep 17 00:00:00 2001 From: RevBayes analysis <1695515+ms609@users.noreply.github.com> Date: Fri, 8 Aug 2025 10:51:00 +0100 Subject: [PATCH 1/8] Replace insertion_sort entirely --- inst/include/TreeTools/renumber_tree.h | 63 ++++---------------------- 1 file changed, 9 insertions(+), 54 deletions(-) diff --git a/inst/include/TreeTools/renumber_tree.h b/inst/include/TreeTools/renumber_tree.h index 5d8274959..ff56d5601 100644 --- a/inst/include/TreeTools/renumber_tree.h +++ b/inst/include/TreeTools/renumber_tree.h @@ -12,58 +12,7 @@ #include "types.h" namespace TreeTools { - inline void swap(int32 *a, int32 *b) { - const int32 temp = *a; - *a = *b; - *b = temp; - } - - inline void insertion_sort_by_smallest(int32* arr, const int32 arr_len, - const int32* sort_by) { - ASSERT(arr_len > 0); - switch (arr_len) { - // case 0: return; - case 1: return; - case 2: - if (sort_by[arr[0]] > sort_by[arr[1]]) { - swap(&arr[0], &arr[1]); - } - return; - } - - for (int32 i = 1; i != arr_len; ++i) { - const int32 - tmp = arr[i], - key = sort_by[tmp] - ; - int32 j = i; - while (j && sort_by[arr[j - 1]] > key) { - arr[j] = arr[j - 1]; - --j; - } - arr[j] = tmp; - } - } - - inline void insertion_sort_by_smallest(std::vector& arr, - const int32 arr_len, - const int32* sort_by) { - ASSERT(arr_len > 0); - switch (arr_len) { - // case 0: - case 1: return; - case 2: - if (sort_by[arr[0]] > sort_by[arr[1]]) { - std::swap(arr[0], arr[1]); - } - return; - default: - std::sort(arr.begin(), arr.end(), [&sort_by](int32 a, int32 b) { - return sort_by[a] < sort_by[b]; - }); - } - } - + struct Frame { int32 node; int32 parent_label; @@ -263,7 +212,10 @@ namespace TreeTools { for (int32 node = n_tip + 1; node < node_limit; ++node) { int32* node_children = children_data + children_start_idx[node]; - insertion_sort_by_smallest(node_children, n_children[node], smallest_desc); + std::sort(node_children, node_children + n_children[node], + [&smallest_desc](int32 a, int32 b) { + return smallest_desc[a] < smallest_desc[b]; + }); } int32 next_label = n_tip + 2; @@ -363,7 +315,10 @@ namespace TreeTools { for (int32 node = n_tip + 1; node < node_limit; ++node) { int32* node_children = children_data + children_start_idx[node]; - insertion_sort_by_smallest(node_children, n_children[node], smallest_desc); + std::sort(node_children, node_children + n_children[node], + [&smallest_desc](int32 a, int32 b) { + return smallest_desc[a] < smallest_desc[b]; + }); } int32 next_label = n_tip + 2; From 4e540084941f6c201f9c8f9e8b2d8c4c49e47dae Mon Sep 17 00:00:00 2001 From: RevBayes analysis <1695515+ms609@users.noreply.github.com> Date: Fri, 8 Aug 2025 11:18:13 +0100 Subject: [PATCH 2/8] Templating underway --- inst/include/TreeTools/renumber_tree.h | 498 +++++++++---------------- inst/include/TreeTools/root_tree.h | 22 +- 2 files changed, 196 insertions(+), 324 deletions(-) diff --git a/inst/include/TreeTools/renumber_tree.h b/inst/include/TreeTools/renumber_tree.h index ff56d5601..67cb0896f 100644 --- a/inst/include/TreeTools/renumber_tree.h +++ b/inst/include/TreeTools/renumber_tree.h @@ -12,341 +12,198 @@ #include "types.h" namespace TreeTools { + +// We'll use this a sentinel type to handle the unweighted case +struct NoWeights {}; +// Used to conditionally create a type +struct DummyDoubleVector {}; + +// Struct to encapsulate memory allocation and pointers +struct TreeData { + int32_t n_edge; + int32_t node_limit; - struct Frame { - int32 node; - int32 parent_label; - int32 child_index; // which child to process next - int32 child_count; - const int32* node_children; - }; + std::vector memory_block; + + int32_t* parent_of; + int32_t* n_children; + int32_t* smallest_desc; + int32_t* children_start_idx; + int32_t* children_data; + + TreeData(int32_t num_edges) + : n_edge(num_edges), + node_limit(num_edges + 2), + memory_block( + node_limit + // parent_of + node_limit + // n_children + node_limit + // smallest_desc + node_limit + // children_start_idx + n_edge, // children_data + 0) + { + auto it = memory_block.begin(); + parent_of = &*it; it += node_limit; + n_children = &*it; it += node_limit; + smallest_desc = &*it; it += node_limit; + children_start_idx = &*it; it += node_limit; + children_data = &*it; + } +}; - inline void add_child_edges( - const int32 root_node, - const int32 root_label, - const int32* children_data, - const int32* children_start_idx, - const int32* n_children, - Rcpp::IntegerMatrix& ret, - int32* next_edge, - int32* next_label) { - - std::stack stack; - - // Initialize with root node - { - int32 child_count = n_children[root_node]; - assert(child_count > 0); - stack.push(Frame{root_node, root_label, 0, child_count, children_data + children_start_idx[root_node]}); +// The core templated function that handles all logic +template +RetType preorder_core( + const Rcpp::IntegerVector& parent, + const Rcpp::IntegerVector& child, + const W& weights) +{ + const int32_t n_edge = parent.length(); + if (R_xlen_t(2LL + child.length() + 2LL + child.length()) > R_xlen_t(INT_FAST32_MAX)) { + Rcpp::stop("Too many edges in tree: Contact 'TreeTools' maintainer for support."); + } + + // Check if parent and child lengths match + if (child.length() != n_edge) { + Rcpp::stop("Length of parent and child must match"); + } + + TreeData data(n_edge); + int32_t root_node = n_edge * 2; + int32_t n_tip = 0; + + // Store edge weights if provided + std::conditional_t, + DummyDoubleVector, std::vector> wt_above_storage; + const std::vector* wt_above_ptr = nullptr; + + if constexpr (!std::is_same_v) { + if (weights.length() != n_edge) { + Rcpp::stop("weights must match number of edges"); } - - while (!stack.empty()) { - Frame& top = stack.top(); - - if (top.child_index == top.child_count) { - // All children processed for this node, pop stack - stack.pop(); - continue; - } - - int32 child = top.node_children[top.child_index]; - - // Write edge info - ret(*next_edge, 0) = top.parent_label; - - if (n_children[child] == 0) { - ret(*next_edge, 1) = child; - ++(*next_edge); - ++top.child_index; - } else { - int32 child_label = *next_label; - ret(*next_edge, 1) = child_label; - ++(*next_label); - ++(*next_edge); - ++top.child_index; - - // Push child frame on stack to process its children next - int32 child_count = n_children[child]; - const int32* child_children = children_data + children_start_idx[child]; - stack.push(Frame{child, child_label, 0, child_count, child_children}); - } + wt_above_storage.resize(data.node_limit); + wt_above_ptr = &wt_above_storage; + } + + for (int32_t i = 0; i < n_edge; ++i) { + const int32_t child_i = child[i]; + const int32_t parent_i = parent[i]; + data.parent_of[child_i] = parent_i; + ++data.n_children[parent_i]; + if constexpr (!std::is_same_v) { + wt_above_storage[child_i] = weights[i]; } } - - inline void add_child_edges( - const int32 root_node, - const int32 root_label, - const int32* children_data, - const int32* children_start_idx, - const int32* n_children, - const std::vector& wt_above, - Rcpp::IntegerMatrix& ret, - Rcpp::NumericVector& weight, - int32* next_edge, - int32* next_label) { - - std::stack stack; - - // Initialize with root node - { - int32 child_count = n_children[root_node]; - assert(child_count > 0); - stack.push(Frame{root_node, root_label, 0, child_count, children_data + children_start_idx[root_node]}); + + int32_t current_idx = 0; + for (int32_t i = 1; i < data.node_limit; i++) { + if (!data.parent_of[i]) { + root_node = i; } - - while (!stack.empty()) { - Frame& top = stack.top(); + if (!data.n_children[i]) { + ++n_tip; + } else { + data.children_start_idx[i] = current_idx; + current_idx += data.n_children[i]; + } + } + + for (int32_t tip = 1; tip < n_tip + 1; ++tip) { + data.smallest_desc[tip] = tip; + int32_t parent_node = data.parent_of[tip]; + while (parent_node && !data.smallest_desc[parent_node]) { + data.smallest_desc[parent_node] = tip; + parent_node = data.parent_of[parent_node]; + } + } + + std::fill(data.n_children, data.n_children + data.node_limit, 0); + for (int32_t i = 0; i < n_edge; ++i) { + int32_t p = parent[i]; + int32_t insert_pos = data.children_start_idx[p] + data.n_children[p]; + data.children_data[insert_pos] = child[i]; + ++data.n_children[p]; + } + + for (int32_t node = n_tip + 1; node < data.node_limit; ++node) { + int32_t* node_children = data.children_data + data.children_start_idx[node]; + std::sort(node_children, node_children + data.n_children[node], + [&](int32_t a, int32_t b) { + return data.smallest_desc[a] < data.smallest_desc[b]; + }); + } + + // Now for the traversal and output generation + int32_t next_edge = 0; + Rcpp::IntegerMatrix ret_edges(n_edge, 2); + + std::conditional_t, + DummyDoubleVector, Rcpp::NumericVector> ret_weights; + + if constexpr (!std::is_same_v) { + ret_weights = Rcpp::NumericVector(n_edge); + } + + std::function add_child_edges = + [&](int32_t node, int32_t current_node_label) { + ret_edges(next_edge, 0) = current_node_label; + ret_edges(next_edge, 1) = node; - if (top.child_index == top.child_count) { - // All children processed for this node, pop stack - stack.pop(); - continue; + if constexpr (!std::is_same_v) { + ret_weights[next_edge] = (*wt_above_ptr)[node]; } - int32 child = top.node_children[top.child_index]; - - // Write edge info - ret(*next_edge, 0) = top.parent_label; - weight[*next_edge] = wt_above[child]; + ++next_edge; - if (n_children[child] == 0) { - ret(*next_edge, 1) = child; - ++(*next_edge); - ++top.child_index; - } else { - int32 child_label = *next_label; - ret(*next_edge, 1) = child_label; - ++(*next_label); - ++(*next_edge); - ++top.child_index; - - // Push child frame on stack to process its children next - int32 child_count = n_children[child]; - const int32* child_children = children_data + children_start_idx[child]; - stack.push(Frame{child, child_label, 0, child_count, child_children}); + const int32_t* node_children = data.children_data + data.children_start_idx[node]; + for (int32_t i = 0; i < data.n_children[node]; ++i) { + add_child_edges(node_children[i], node); } - } - } - - // [[Rcpp::export]] - inline Rcpp::IntegerMatrix preorder_edges_and_nodes( - const Rcpp::IntegerVector parent, - const Rcpp::IntegerVector child) - { - if (R_xlen_t(2LL + child.length() + 2LL + child.length()) > R_xlen_t(INT_FAST32_MAX)) { - Rcpp::stop("Too many edges in tree: " // #nocov - "Contact 'TreeTools' maintainer for support."); // #nocov - } + }; + + add_child_edges(root_node, n_tip + 1); - ASSERT(parent.length() < INT_FAST32_MAX - 2); - const int32 n_edge = int32(parent.length()); - if (child.length() != n_edge) { - Rcpp::stop("Length of parent and child must match"); + if constexpr (std::is_same_v) { + return ret_edges; + } else { + return std::make_pair(ret_edges, ret_weights); } - const int32 max_node = n_edge + 1; - assert(max_node == *std::max_element(parent.begin(), parent.end())); - const int32 node_limit = max_node + 1; +} + +// === PUBLIC EXPORTED FUNCTIONS === + +// [[Rcpp::export]] +inline Rcpp::IntegerMatrix preorder_edges_and_nodes( + const Rcpp::IntegerVector parent, + const Rcpp::IntegerVector child) +{ + return preorder_core(parent, child, NoWeights{}); +} + +// [[Rcpp::export]] +inline Rcpp::List preorder_weighted( + const Rcpp::IntegerVector& parent, + const Rcpp::IntegerVector& child, + const Rcpp::DoubleVector& weight) +{ + std::pair result = + preorder_core>(parent, child, weight); + + Rcpp::List ret = Rcpp::List::create( + Rcpp::Named("edge") = result.first, + Rcpp::Named("edge.length") = result.second + ); + + return ret; +} + - int32 next_edge = 0; - int32 root_node = n_edge * 2; /* Initialize with too-big value */ - int32 n_tip = 0; - - // Single large allocation instead of many small ones - const size_t total_ints_needed = - node_limit + // parent_of - node_limit + // n_children - node_limit + // smallest_desc - node_limit + // children_start_idx - n_edge; // children_data - - std::vector memory_block(total_ints_needed, 0); - - auto it = memory_block.begin(); - int32* parent_of = &*it; - it += node_limit; - int32* n_children = &*it; - it += node_limit; - int32* smallest_desc = &*it; - it += node_limit; - int32* children_start_idx = &*it; - it += node_limit; - int32* children_data = &*it; - - for (int32 i = n_edge; i--; ) { - const int32 parent_i = parent[i]; - parent_of[child[i]] = parent_i; - ++n_children[parent_i]; - } - - int32 current_idx = 0; - for (int32 i = 1; i < node_limit; i++) { - if (!parent_of[i]) { - root_node = i; - } - if (!n_children[i]) { - ++n_tip; - } else { - children_start_idx[i] = current_idx; - current_idx += n_children[i]; - } - } - for (int32 tip = 1; tip < n_tip + 1; ++tip) { - smallest_desc[tip] = tip; - int32 parent = parent_of[tip]; - while (!smallest_desc[parent]) { - smallest_desc[parent] = tip; - parent = parent_of[parent]; - } - } - - // Reset n_children - use as insertion counter - std::fill(n_children, n_children + node_limit, 0); - for (int32 i = 0; i < n_edge; ++i) { - int32 p = parent[i]; - int32 insert_pos = children_start_idx[p] + n_children[p]; - children_data[insert_pos] = child[i]; - ++n_children[p]; - } - for (int32 node = n_tip + 1; node < node_limit; ++node) { - int32* node_children = children_data + children_start_idx[node]; - std::sort(node_children, node_children + n_children[node], - [&smallest_desc](int32 a, int32 b) { - return smallest_desc[a] < smallest_desc[b]; - }); - } - - int32 next_label = n_tip + 2; - Rcpp::IntegerMatrix ret(n_edge, 2); - add_child_edges(root_node, n_tip + 1, children_data, children_start_idx, - n_children, ret, &next_edge, &next_label); - return ret; - } - inline std::pair preorder_weighted_pair( - const Rcpp::IntegerVector& parent, - const Rcpp::IntegerVector& child, - const Rcpp::DoubleVector& weight) - { - if (R_xlen_t(2LL + child.length() + 2LL + child.length()) > - R_xlen_t(INT_FAST32_MAX)) { - Rcpp::stop("Too many edges in tree: " // #nocov - "Contact 'TreeTools' maintainer for support."); // #nocov - } - - ASSERT(parent.length() < INT_FAST32_MAX - 2); - const int32 n_edge = int32(parent.length()); - const int32 node_limit = n_edge + 2; - - if (child.length() != n_edge) { - Rcpp::stop("Length of parent and child must match"); - } - if (weight.length() != n_edge) { - Rcpp::stop("weights must match number of edges"); - } - - int32 next_edge = 0; - int32 root_node = n_edge * 2; /* Initialize with too-big value */ - int32 n_tip = 0; - - const size_t total_ints_needed = - node_limit + // parent_of - node_limit + // n_children - node_limit + // smallest_desc - node_limit + // children_start_idx - n_edge; // children_data - - std::vector memory_block(total_ints_needed, 0); - - auto it = memory_block.begin(); - int32* parent_of = &*it; - it += node_limit; - int32* n_children = &*it; - it += node_limit; - int32* smallest_desc = &*it; - it += node_limit; - int32* children_start_idx = &*it; - it += node_limit; - int32* children_data = &*it; - - std::vector wt_above(node_limit); - - for (int32 i = 0; i < n_edge; ++i) { - const int32 child_i = child[i]; - const int32 parent_i = parent[i]; - wt_above[child_i] = weight[i]; - parent_of[child_i] = parent_i; - ++n_children[parent_i]; - } - - int32 current_idx = 0; - for (int32 i = 1; i < node_limit; i++) { - if (!parent_of[i]) { - root_node = i; - } - if (!n_children[i]) { - ++n_tip; - } else { - children_start_idx[i] = current_idx; - current_idx += n_children[i]; - } - } - - for (int32 tip = 1; tip < n_tip + 1; ++tip) { - smallest_desc[tip] = tip; - int32 parent = parent_of[tip]; - while (!smallest_desc[parent]) { - smallest_desc[parent] = tip; - parent = parent_of[parent]; - } - } - - // Reset n_children - use as insertion counter - std::fill(n_children, n_children + node_limit, 0); - for (int32 i = 0; i < n_edge; ++i) { - int32 p = parent[i]; - int32 insert_pos = children_start_idx[p] + n_children[p]; - children_data[insert_pos] = child[i]; - ++n_children[p]; - } - - for (int32 node = n_tip + 1; node < node_limit; ++node) { - int32* node_children = children_data + children_start_idx[node]; - std::sort(node_children, node_children + n_children[node], - [&smallest_desc](int32 a, int32 b) { - return smallest_desc[a] < smallest_desc[b]; - }); - } - - int32 next_label = n_tip + 2; - Rcpp::IntegerMatrix ret(n_edge, 2); - Rcpp::NumericVector ret_wt(n_edge); - add_child_edges(root_node, n_tip + 1, children_data, children_start_idx, - n_children, wt_above, ret, ret_wt, &next_edge, &next_label); - - return std::make_pair(ret, ret_wt); - } - inline Rcpp::List preorder_weighted( - const Rcpp::IntegerVector& parent, - const Rcpp::IntegerVector& child, - const Rcpp::DoubleVector& weight) - { - // Call the core function to get the pair - std::pair result = - preorder_weighted_pair(parent, child, weight); - - // Manually create an Rcpp::List and populate it - Rcpp::List ret = Rcpp::List::create( - Rcpp::Named("edge") = result.first, - Rcpp::Named("edge.length") = result.second - ); - - return ret; - } inline int32 get_subtree_size(int32 node, int32 *subtree_size, int32 *n_children, int32 **children_of, @@ -359,6 +216,17 @@ namespace TreeTools { } return subtree_size[node]; } + + + + + + + + + + + template struct SmallBuffer { diff --git a/inst/include/TreeTools/root_tree.h b/inst/include/TreeTools/root_tree.h index 464964e4a..2790fe84b 100644 --- a/inst/include/TreeTools/root_tree.h +++ b/inst/include/TreeTools/root_tree.h @@ -13,11 +13,11 @@ namespace TreeTools { const Rcpp::IntegerVector parent, const Rcpp::IntegerVector child); - extern inline std::pair - preorder_weighted_pair( + template + extern inline RetType preorder_core( const Rcpp::IntegerVector& parent, const Rcpp::IntegerVector& child, - const Rcpp::DoubleVector& weight); + const W& weights); // edge must be BINARY // edge must be in preorder @@ -110,10 +110,14 @@ namespace TreeTools { const bool weighted = phy.containsElementNamed("edge.length"); if (weighted) { - std::tie(edge, weight) = preorder_weighted_pair( - edge(Rcpp::_, 0), - edge(Rcpp::_, 1), - phy["edge.length"] + + Rcpp::IntegerVector parent_col(edge(Rcpp::_, 0)); + Rcpp::IntegerVector child_col(edge(Rcpp::_, 1)); + Rcpp::NumericVector edge_len = phy["edge.length"]; + std::tie(edge, weight) = preorder_core>( + parent_col, + child_col, + edge_len ); } else { edge = preorder_edges_and_nodes(edge(Rcpp::_, 0), edge(Rcpp::_, 1)); @@ -176,7 +180,7 @@ namespace TreeTools { new_edge(root_edges[spare_edge], 1) = outgroup; if (weighted) { Rcpp::List preorder_res; - auto [edge, edge_weight] = preorder_weighted_pair(new_edge(Rcpp::_, 0), + auto [edge, edge_weight] = preorder_core(new_edge(Rcpp::_, 0), new_edge(Rcpp::_, 1), weight); ret["edge"] = edge; @@ -217,7 +221,7 @@ namespace TreeTools { ret["Nnode"] = n_node + 1; if (weighted) { Rcpp::List preorder_res; - auto [edge, weight] = preorder_weighted_pair( + auto [edge, weight] = preorder_core( new_edge(Rcpp::_, 0), new_edge(Rcpp::_, 1), new_wt); From 5bd1890663d0913a2073c055568ed25ab96eb638 Mon Sep 17 00:00:00 2001 From: RevBayes analysis <1695515+ms609@users.noreply.github.com> Date: Fri, 8 Aug 2025 11:19:01 +0100 Subject: [PATCH 3/8] #include --- inst/include/TreeTools/renumber_tree.h | 1 + 1 file changed, 1 insertion(+) diff --git a/inst/include/TreeTools/renumber_tree.h b/inst/include/TreeTools/renumber_tree.h index 67cb0896f..183649823 100644 --- a/inst/include/TreeTools/renumber_tree.h +++ b/inst/include/TreeTools/renumber_tree.h @@ -7,6 +7,7 @@ #include /* for stack */ #include /* for errors */ #include /* for pair */ +#include #include #include "assert.h" /* for ASSERT */ #include "types.h" From 97266056d52bf10ea0380b7ae1dfeb64474f93e9 Mon Sep 17 00:00:00 2001 From: RevBayes analysis <1695515+ms609@users.noreply.github.com> Date: Fri, 8 Aug 2025 11:44:30 +0100 Subject: [PATCH 4/8] Almost working Just need checks? --- inst/include/TreeTools/renumber_tree.h | 84 +++++++++++++++++--------- 1 file changed, 54 insertions(+), 30 deletions(-) diff --git a/inst/include/TreeTools/renumber_tree.h b/inst/include/TreeTools/renumber_tree.h index 183649823..d8f06b429 100644 --- a/inst/include/TreeTools/renumber_tree.h +++ b/inst/include/TreeTools/renumber_tree.h @@ -19,7 +19,6 @@ struct NoWeights {}; // Used to conditionally create a type struct DummyDoubleVector {}; -// Struct to encapsulate memory allocation and pointers struct TreeData { int32_t n_edge; int32_t node_limit; @@ -52,7 +51,14 @@ struct TreeData { } }; -// The core templated function that handles all logic +struct Frame { + int32_t node; + int32_t parent_label; + int32_t child_index; + int32_t child_count; + const int32_t* node_children; +}; + template RetType preorder_core( const Rcpp::IntegerVector& parent, @@ -64,7 +70,6 @@ RetType preorder_core( Rcpp::stop("Too many edges in tree: Contact 'TreeTools' maintainer for support."); } - // Check if parent and child lengths match if (child.length() != n_edge) { Rcpp::stop("Length of parent and child must match"); } @@ -73,9 +78,7 @@ RetType preorder_core( int32_t root_node = n_edge * 2; int32_t n_tip = 0; - // Store edge weights if provided - std::conditional_t, - DummyDoubleVector, std::vector> wt_above_storage; + std::conditional_t, DummyDoubleVector, std::vector> wt_above_storage; const std::vector* wt_above_ptr = nullptr; if constexpr (!std::is_same_v) { @@ -134,41 +137,63 @@ RetType preorder_core( }); } - // Now for the traversal and output generation int32_t next_edge = 0; + int32_t next_label = n_tip + 2; + Rcpp::IntegerMatrix ret_edges(n_edge, 2); - std::conditional_t, - DummyDoubleVector, Rcpp::NumericVector> ret_weights; + std::conditional_t, DummyDoubleVector, Rcpp::NumericVector> ret_weights; if constexpr (!std::is_same_v) { ret_weights = Rcpp::NumericVector(n_edge); } - std::function add_child_edges = - [&](int32_t node, int32_t current_node_label) { - ret_edges(next_edge, 0) = current_node_label; - ret_edges(next_edge, 1) = node; - - if constexpr (!std::is_same_v) { - ret_weights[next_edge] = (*wt_above_ptr)[node]; - } - - ++next_edge; - - const int32_t* node_children = data.children_data + data.children_start_idx[node]; - for (int32_t i = 0; i < data.n_children[node]; ++i) { - add_child_edges(node_children[i], node); - } - }; + std::stack stack; + + // Initialize with root node children + { + int32_t child_count = data.n_children[root_node]; + if (child_count > 0) { + stack.push(Frame{root_node, n_tip + 1, 0, child_count, data.children_data + data.children_start_idx[root_node]}); + } + } + + while (!stack.empty()) { + Frame& top = stack.top(); + + if (top.child_index == top.child_count) { + stack.pop(); + continue; + } - add_child_edges(root_node, n_tip + 1); + int32_t child_node = top.node_children[top.child_index]; - if constexpr (std::is_same_v) { - return ret_edges; + ret_edges(next_edge, 0) = top.parent_label; + if constexpr (!std::is_same_v) { + ret_weights[next_edge] = (*wt_above_ptr)[child_node]; + } + + if (data.n_children[child_node] == 0) { + ret_edges(next_edge, 1) = child_node; + ++next_edge; + ++top.child_index; } else { - return std::make_pair(ret_edges, ret_weights); + int32_t child_label = next_label++; + ret_edges(next_edge, 1) = child_label; + ++next_edge; + ++top.child_index; + + int32_t child_count = data.n_children[child_node]; + const int32_t* child_children = data.children_data + data.children_start_idx[child_node]; + stack.push(Frame{child_node, child_label, 0, child_count, child_children}); } + } + + if constexpr (std::is_same_v) { + return ret_edges; + } else { + return std::make_pair(ret_edges, ret_weights); + } } // === PUBLIC EXPORTED FUNCTIONS === @@ -205,7 +230,6 @@ inline Rcpp::List preorder_weighted( - inline int32 get_subtree_size(int32 node, int32 *subtree_size, int32 *n_children, int32 **children_of, int32 n_edge) { From a3106efd667c0251a11599a4e70a5bf9683ecec0 Mon Sep 17 00:00:00 2001 From: RevBayes analysis <1695515+ms609@users.noreply.github.com> Date: Fri, 8 Aug 2025 11:50:05 +0100 Subject: [PATCH 5/8] root_label --- inst/include/TreeTools/renumber_tree.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/inst/include/TreeTools/renumber_tree.h b/inst/include/TreeTools/renumber_tree.h index d8f06b429..243d32f7a 100644 --- a/inst/include/TreeTools/renumber_tree.h +++ b/inst/include/TreeTools/renumber_tree.h @@ -149,12 +149,13 @@ RetType preorder_core( } std::stack stack; + int32_t root_label = n_tip + 1; // Initialize with root node children { int32_t child_count = data.n_children[root_node]; if (child_count > 0) { - stack.push(Frame{root_node, n_tip + 1, 0, child_count, data.children_data + data.children_start_idx[root_node]}); + stack.push(Frame{root_node, root_label, 0, child_count, data.children_data + data.children_start_idx[root_node]}); } } @@ -195,7 +196,6 @@ RetType preorder_core( return std::make_pair(ret_edges, ret_weights); } } - // === PUBLIC EXPORTED FUNCTIONS === // [[Rcpp::export]] From 7595b4965908aa58d22a57f8999b5836ad29f5d8 Mon Sep 17 00:00:00 2001 From: RevBayes analysis <1695515+ms609@users.noreply.github.com> Date: Fri, 8 Aug 2025 12:39:02 +0100 Subject: [PATCH 6/8] complete merge --- inst/include/TreeTools/renumber_tree.h | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/inst/include/TreeTools/renumber_tree.h b/inst/include/TreeTools/renumber_tree.h index e27bc8005..e509c1dba 100644 --- a/inst/include/TreeTools/renumber_tree.h +++ b/inst/include/TreeTools/renumber_tree.h @@ -128,8 +128,7 @@ RetType preorder_core( data.children_data[insert_pos] = child[i]; ++data.n_children[p]; } -HEAD - + for (int32_t node = n_tip + 1; node < data.node_limit; ++node) { int32_t* node_children = data.children_data + data.children_start_idx[node]; std::sort(node_children, node_children + data.n_children[node], From 0e8b7069fff305344239ea8b2dd5138d5ab0bc0b Mon Sep 17 00:00:00 2001 From: RevBayes analysis <1695515+ms609@users.noreply.github.com> Date: Fri, 8 Aug 2025 12:40:33 +0100 Subject: [PATCH 7/8] preorder_core --- inst/include/TreeTools/root_tree.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/inst/include/TreeTools/root_tree.h b/inst/include/TreeTools/root_tree.h index c7b7a5cc0..2dd7a27ec 100644 --- a/inst/include/TreeTools/root_tree.h +++ b/inst/include/TreeTools/root_tree.h @@ -181,7 +181,7 @@ namespace TreeTools { new_child[root_edges[spare_edge]] = outgroup; if (weighted) { - std::tie(edge, weight) = preorder_core(new_parent, new_child, weight); + std::tie(edge, weight) = preorder_core>(new_parent, new_child, weight); ret["edge"] = edge; ret["edge.length"] = weight; } else { @@ -211,7 +211,7 @@ namespace TreeTools { Rcpp::NumericVector new_wt(n_edge + 1); std::copy(weight.begin(), weight.end(), new_wt.begin()); new_wt[n_edge] = 0; - std::tie(edge, weight) = preorder_core(new_parent, new_child, new_wt); + std::tie(edge, weight) = preorder_core>(new_parent, new_child, new_wt); ret["edge"] = edge; ret["edge.length"] = weight; } else { From 2d90ef3fc71012c4171188e6fb66ead337331907 Mon Sep 17 00:00:00 2001 From: RevBayes analysis <1695515+ms609@users.noreply.github.com> Date: Fri, 8 Aug 2025 13:12:53 +0100 Subject: [PATCH 8/8] Mark inline --- inst/include/TreeTools/renumber_tree.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inst/include/TreeTools/renumber_tree.h b/inst/include/TreeTools/renumber_tree.h index e509c1dba..fe1c0f92b 100644 --- a/inst/include/TreeTools/renumber_tree.h +++ b/inst/include/TreeTools/renumber_tree.h @@ -60,7 +60,7 @@ struct Frame { }; template -RetType preorder_core( +inline RetType preorder_core( const Rcpp::IntegerVector& parent, const Rcpp::IntegerVector& child, const W& weights)