Merge branch 'transfer-consensus' into expose-lapjv

Author: ms609

src/day_1985.cpp

@@ -3,6 +3,8 @@ using namespace Rcpp;
 
 #include "tree_distances.h" /* includes <TreeTools/SplitList.h> */
 #include "information.h"
+#include <cassert>
+#include <TreeTools/assert.h>
 
 #include <TreeTools.h> /* for root_on_node() */
 #include <TreeTools/root_tree.h> /* for root_on_node() */
@@ -96,7 +98,7 @@ double calc_consensus_info(const List &trees, const LogicalVector &phylo,
 
   std::vector<ClusterTable> tables;
   if (std::size_t(n_trees) > tables.max_size()) {
-    Rcpp::stop("Not enough memory available to compute consensus of so many trees"); // LCOV_EXCL_LINE
+    ASSERT(false && "Not enough memory for consensus of so many trees"); // LCOV_EXCL_LINE
   }
 
   tables.reserve(n_trees);
@@ -321,11 +323,9 @@ IntegerVector robinson_foulds_all_pairs(const List& tables) {
 // [[Rcpp::export]]
 double consensus_info(const List trees, const LogicalVector phylo,
                       const NumericVector p) {
-  if (p[0] > 1 + 1e-15) { // epsilon catches floating point error
-    Rcpp::stop("p must be <= 1.0 in consensus_info()");
-  } else if (p[0] < 0.5) {
-    Rcpp::stop("p must be >= 0.5 in consensus_info()");
-  }
+  // Validated by R caller (ConsensusInfo checks p range)
+  ASSERT(p[0] <= 1 + 1e-15 && "p must be <= 1.0 in consensus_info()");
+  ASSERT(p[0] >= 0.5 && "p must be >= 0.5 in consensus_info()");
 
   // Peek at tree size to choose stack vs heap allocation for the work buffer
   ClusterTable temp_table(Rcpp::List(trees(0)));

src/hmi.cpp

@@ -139,9 +139,8 @@ double hierarchical_self_info(const std::vector<TreeDist::HNode>& nodes, size_t
 double HMI_xptr(SEXP ptr1, SEXP ptr2) {
   Rcpp::XPtr<TreeDist::HPart> hp1(ptr1);
   Rcpp::XPtr<TreeDist::HPart> hp2(ptr2);
-  if (hp1->nodes[hp1->root].n_tip != hp2->nodes[hp2->root].n_tip) {
-    Rcpp::stop("Trees must have the same number of leaves");
-  }
+  ASSERT(hp1->nodes[hp1->root].n_tip == hp2->nodes[hp2->root].n_tip
+         && "Trees must have the same number of leaves");
   return TreeDist::hierarchical_mutual_info(hp1->nodes, hp1->root,
                                             hp2->nodes, hp2->root).second;
 }
@@ -170,12 +169,8 @@ Rcpp::NumericVector EHMI_xptr(SEXP hp1_ptr, SEXP hp2_ptr,
                               double tolerance = 0.01,
                               int minResample = 36) {
 
-  if (minResample < 2) {
-    Rcpp::stop("Must perform at least one resampling");
-  }
-  if (tolerance < 1e-8) {
-    Rcpp::stop("Tolerance too low");
-  }
+  ASSERT(minResample >= 2 && "Must perform at least one resampling");
+  ASSERT(tolerance >= 1e-8 && "Tolerance too low");
 
   Rcpp::XPtr<TreeDist::HPart> hp1(hp1_ptr);
   Rcpp::XPtr<TreeDist::HPart> hp2(hp2_ptr);

src/hpart.cpp

@@ -126,7 +126,7 @@ size_t build_node_from_list(const RObject& node,
   if (Rf_isInteger(node) || Rf_isReal(node)) {
     const IntegerVector leaf_vec(node);
     if (leaf_vec.size() != 1) {
-      Rcpp::stop("List must only contain integers, not vectors of integers");   // #nocov
+      ASSERT(false && "List must only contain integers, not vectors of integers"); // #nocov
     }
     const int leaf_label = leaf_vec[0];         // 1-based R leaf label
     const size_t leaf_idx = leaf_label - 1;     // 0-based label for HNode
@@ -172,7 +172,7 @@ size_t build_node_from_list(const RObject& node,
   }
 
   // Invalid node type
-  Rcpp::stop("Invalid node type");                                              // #nocov
+  ASSERT(false && "Invalid node type");                                        // #nocov
 }
 
 

src/hpart_relabel.cpp

@@ -13,7 +13,7 @@ void recompute_bitsets_postorder(TreeDist::HPart &hpart, const size_t node_idx,
   if (node.children.empty()) {
     // Leaf node
     if (node.leaf_count != 1) {
-      Rcpp::stop("Leaf node has leaf_count != 1");                              // #nocov
+      ASSERT(node.leaf_count == 1 && "Leaf node has leaf_count != 1");         // #nocov
     }
     int new_index = mapping[node.label]; // mapping is 0-based
     node.label = new_index;

src/information.h

@@ -2,7 +2,7 @@
 #define _TREEDIST_INFO_H
 
 #include <Rcpp/Lightest>
-#include <cassert> /* for log2() */
+#include <TreeTools/assert.h>
 #include <cmath> /* for log2() */
 #include "ints.h" /* for int16 */
 
@@ -48,14 +48,12 @@ __attribute__((constructor))
 
 inline double split_phylo_info(const int32 n_in, const int32 *n_tip,
                                const double p) {
-  if (*n_tip > CT_MAX_LEAVES) {
-    Rcpp::stop("This many leaves are not yet supported.");
-  }
+  ASSERT(*n_tip <= CT_MAX_LEAVES && "This many leaves are not yet supported.");
   const int32 n_out = *n_tip - n_in;
-  assert(p > 0);
-  assert(p <= 1);
-  assert(n_in > 1);
-  assert(n_out > 1);
+  ASSERT(p > 0);
+  ASSERT(p <= 1);
+  ASSERT(n_in > 1);
+  ASSERT(n_out > 1);
   if (p == 1) {
     return (l2unrooted[*n_tip] - l2rooted[n_in] - l2rooted[n_out]);
   } else {
@@ -74,14 +72,12 @@ inline double split_phylo_info(const int32 n_in, const int32 *n_tip,
 
 inline double split_clust_info(const int32 n_in, const int32 *n_tip,
                                const double p) {
-  if (*n_tip > CT_MAX_LEAVES) {
-    Rcpp::stop("This many leaves are not yet supported.");
-  }
+  ASSERT(*n_tip <= CT_MAX_LEAVES && "This many leaves are not yet supported.");
   const int32 n_out = *n_tip - n_in;
-  assert(p > 0);
-  assert(p <= 1);
-  assert(n_in > 1);
-  assert(n_out > 1);
+  ASSERT(p > 0);
+  ASSERT(p <= 1);
+  ASSERT(n_in > 1);
+  ASSERT(n_out > 1);
   return -p * (
       (((n_in * l2[n_in]) + (n_out * l2[n_out])) / *n_tip) - l2[*n_tip]
   );

src/mast.cpp

@@ -1,5 +1,6 @@
 #include <cstdlib>
 #include <Rcpp/Lightest>
+#include <TreeTools/assert.h>
 #include "ints.h"
 using namespace Rcpp;
 
@@ -63,12 +64,8 @@ int cpp_mast (IntegerMatrix edge1, IntegerMatrix edge2, IntegerVector nTip) {
               n_internal = n_tip - 1,
               n_all_nodes = n_tip + n_internal,
               n_edge = edge1.nrow();
-  if (edge2.nrow() != n_edge) {
-    Rcpp::stop("Both trees must contain the same number of edges.");
-  }
-  if (n_tip > MAST_MAX_TIP) {
-    Rcpp::stop("Tree too large; please contact maintainer for advice.");
-  }
+  ASSERT(edge2.nrow() == n_edge && "Both trees must contain the same number of edges.");
+  ASSERT(n_tip <= MAST_MAX_TIP && "Tree too large; please contact maintainer for advice.");
 
   int16
     t1_left[MAST_MAX_NODE] = {}, t1_right[MAST_MAX_NODE] = {},

src/nni_distance.cpp

@@ -226,16 +226,7 @@ grf_match nni_rf_matching (
 
     ASSERT(n_splits > 0);
     ASSERT(n_tips > 3);
-    if (n_tips > NNI_MAX_TIPS) {
-      Rcpp::stop("Cannot calculate NNI distance for trees with so many tips.");
-    }
-    
-    // #nocov start
-    if (static_cast<int64_t>(n_splits) * static_cast<int64_t>(n_bins) >
-          static_cast<int64_t>(std::numeric_limits<int32_t>::max())) {
-      Rcpp::stop("Cannot calculate NNI distance for trees with so many splits.");
-    }
-    // #nocov end
+    ASSERT(n_tips <= NNI_MAX_TIPS); // Validated by exported wrapper
 
     const int32_t last_bin = n_bins - 1;
     const int32_t unset_tips = (n_tips % SL_BIN_SIZE) ? 
@@ -298,9 +289,7 @@ IntegerVector cpp_nni_distance(const IntegerMatrix& edge1,
                                const IntegerMatrix& edge2,
                                const IntegerVector& nTip) {
 
-  if (nTip[0] > NNI_MAX_TIPS) {
-    Rcpp::stop("Cannot calculate NNI distance for trees with so many tips.");
-  }
+  ASSERT(nTip[0] <= NNI_MAX_TIPS && "Cannot calculate NNI distance for trees with so many tips.");
   const int32_t n_tip = static_cast<int32_t>(nTip[0]);
   const int32_t node_0 = n_tip;
   const int32_t node_0_r = n_tip + 1;
@@ -315,10 +304,8 @@ IntegerVector cpp_nni_distance(const IntegerMatrix& edge1,
   int32_t fack_score_bound = 0;
   int32_t li_score_bound = 0;
 
-  if (n_edge != int32_t(edge2.nrow())) {
-    Rcpp::stop("Both trees must have the same number of edges. "
-             "Is one rooted and the other unrooted?");
-  }
+  ASSERT(n_edge == int32_t(edge2.nrow())
+         && "Both trees must have the same number of edges.");
 
   if (n_tip < 4) {
     return(IntegerVector::create(Named("lower") = 0,
@@ -348,7 +335,7 @@ IntegerVector cpp_nni_distance(const IntegerMatrix& edge1,
   const int32_t n_splits = n_distinct_edge - n_tip;
 
   if (n_splits < 1) {
-    Rcpp::stop("NNI distance is undefined for trees with no splits"); // #nocov
+    ASSERT(false && "NNI distance is undefined for trees with no splits"); // #nocov
   }
 
   std::unique_ptr<uint64_t[]> splits1(new uint64_t[n_splits * n_bin]);

src/spr.cpp

@@ -47,18 +47,11 @@ IntegerVector calc_mismatch_size(const RawMatrix& x, const RawMatrix& y) {
 
 // [[Rcpp::export]]
 IntegerVector mismatch_size (const RawMatrix& x, const RawMatrix& y) {
-  if (x.rows() != y.rows()) {
-    Rcpp::stop("`x` and `y` differ in number of splits.");
-  }
-  if (!x.hasAttribute("nTip")) {
-    Rcpp::stop("`x` lacks nTip attribute");
-  }
-  if (!y.hasAttribute("nTip")) {
-    Rcpp::stop("`y` lacks nTip attribute");
-  }
-  if (static_cast<int>(x.attr("nTip")) != static_cast<int>(y.attr("nTip"))) {
-    Rcpp::stop("`x` and `y` differ in `nTip`");
-  }
+  ASSERT(x.rows() == y.rows() && "`x` and `y` differ in number of splits.");
+  ASSERT(x.hasAttribute("nTip") && "`x` lacks nTip attribute");
+  ASSERT(y.hasAttribute("nTip") && "`y` lacks nTip attribute");
+  ASSERT(static_cast<int>(x.attr("nTip")) == static_cast<int>(y.attr("nTip"))
+         && "`x` and `y` differ in `nTip`");
   return calc_mismatch_size(x, y);
 }
 
@@ -102,25 +95,17 @@ IntegerVector calc_confusion(const RawMatrix &x, const RawMatrix &y) {
 
 // [[Rcpp::export]]
 IntegerVector confusion(const RawMatrix& x, const RawMatrix& y) {
-  if (x.rows() != y.rows()) {
-    Rcpp::stop("Input splits must contain same number of splits.");
-  }
-  if (!x.hasAttribute("nTip")) {
-    Rcpp::stop("`x` lacks nTip attribute");
-  }
-  if (!y.hasAttribute("nTip")) {
-    Rcpp::stop("`y` lacks nTip attribute");
-  }
-  if (static_cast<int>(x.attr("nTip")) != static_cast<int>(y.attr("nTip"))) {
-    Rcpp::stop("`x` and `y` differ in `nTip`");
-  }
+  ASSERT(x.rows() == y.rows() && "Input splits must contain same number of splits.");
+  ASSERT(x.hasAttribute("nTip") && "`x` lacks nTip attribute");
+  ASSERT(y.hasAttribute("nTip") && "`y` lacks nTip attribute");
+  ASSERT(static_cast<int>(x.attr("nTip")) == static_cast<int>(y.attr("nTip"))
+         && "`x` and `y` differ in `nTip`");
   return calc_confusion(x, y);
 }
 
 IntegerMatrix reverse (const IntegerMatrix x) {
-  if (double(x.nrow()) > double(std::numeric_limits<intx>::max())) {
-    Rcpp::stop("This many edges are not (yet) supported.");
-  }
+  ASSERT(double(x.nrow()) <= double(std::numeric_limits<intx>::max())
+         && "This many edges are not (yet) supported.");
   const intx n_edge = intx(x.nrow());
   ASSERT(n_edge % 2 == 0); // Tree is binary
   IntegerMatrix ret(n_edge, 2);

src/spr_lookup.cpp

@@ -223,7 +223,7 @@ int lookup_7(const SplitSet7& sp1, const SplitSet7& sp2) {
 
 inline SplitSet7 read_splits(const Rcpp::RawVector& r) {
   if (r.size() != 4)
-    Rcpp::stop("Expected a length-4 raw vector of splits");
+    ASSERT(false && "Expected a length-4 raw vector of splits");
 
   SplitSet7 sp{};
   for (int i = 0; i < 4; ++i) {

src/transfer_distance.cpp

@@ -12,6 +12,7 @@
  */
 
 #include <TreeTools/SplitList.h>
+#include <TreeTools/assert.h>
 #include <Rcpp/Lightest>
 #include <algorithm>
 #include <cstdint>
@@ -145,9 +146,8 @@ List cpp_transfer_dist(
     const RawMatrix& x, const RawMatrix& y,
     const IntegerVector& nTip
 ) {
-  if (x.cols() != y.cols()) {
-    Rcpp::stop("Input splits must address same number of tips.");
-  }
+  // Validated by R caller (.ValidateSplitArgs or as.Splits with shared tipLabels)
+  ASSERT(x.cols() == y.cols() && "Input splits must address same number of tips.");
 
   const int n_tip = nTip[0];
   SplitList sl_x(x);
@@ -204,9 +204,8 @@ List cpp_transfer_dist_scored(
     const RawMatrix& x, const RawMatrix& y,
     const IntegerVector& nTip, bool scale
 ) {
-  if (x.cols() != y.cols()) {
-    Rcpp::stop("Input splits must address same number of tips.");
-  }
+  // Validated by R caller (.ValidateSplitArgs or as.Splits with shared tipLabels)
+  ASSERT(x.cols() == y.cols() && "Input splits must address same number of tips.");
 
   const int n_tip = nTip[0];
   SplitList sl_x(x);