diff --git a/R/RcppExports.R b/R/RcppExports.R index 1abe44b6..1adeb2a1 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -57,6 +57,10 @@ cpp_msi_distance <- function(x, y, nTip) { .Call(`_TreeDist_cpp_msi_distance`, x, y, nTip) } +cpp_all_pairs_mci <- function(x, nTip) { + .Call(`_TreeDist_cpp_all_pairs_mci`, x, nTip) +} + cpp_mutual_clustering <- function(x, y, nTip) { .Call(`_TreeDist_cpp_mutual_clustering`, x, y, nTip) } diff --git a/R/tree_distance_info.R b/R/tree_distance_info.R index 9f7548e5..5ca371ab 100644 --- a/R/tree_distance_info.R +++ b/R/tree_distance_info.R @@ -297,8 +297,15 @@ ExpectedVariation <- function (tree1, tree2, samples = 1e+4) { #' @export MutualClusteringInfo <- function (tree1, tree2 = NULL, normalize = FALSE, reportMatching = FALSE, diag = TRUE) { - unnormalized <- CalculateTreeDistance(MutualClusteringInfoSplits, tree1, tree2, - reportMatching) + unnormalized <- if (is.null(tree2) && is.list(tree1) && + !inherits(tree1, 'phylo')) { + structure(cpp_all_pairs_mci(as.Splits(tree1, tree1[[1]]), NTip(tree1[[1]])), + class = "dist", + Size = length(tree1), diag = FALSE, upper = FALSE) + } else { + CalculateTreeDistance(MutualClusteringInfoSplits, tree1, tree2, + reportMatching) + } if (diag && is.null(tree2)) { unnormalized <- as.matrix(unnormalized) diag(unnormalized) <- ClusteringEntropy(tree1) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 9989fd52..2afb43d7 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -177,6 +177,18 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// cpp_all_pairs_mci +NumericVector cpp_all_pairs_mci(const List x, const IntegerVector nTip); +RcppExport SEXP _TreeDist_cpp_all_pairs_mci(SEXP xSEXP, SEXP nTipSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const List >::type x(xSEXP); + Rcpp::traits::input_parameter< const IntegerVector >::type nTip(nTipSEXP); + rcpp_result_gen = Rcpp::wrap(cpp_all_pairs_mci(x, nTip)); + return rcpp_result_gen; +END_RCPP +} // cpp_mutual_clustering List cpp_mutual_clustering(const RawMatrix x, const RawMatrix y, const IntegerVector nTip); RcppExport SEXP _TreeDist_cpp_mutual_clustering(SEXP xSEXP, SEXP ySEXP, SEXP nTipSEXP) { @@ -219,6 +231,7 @@ static const R_CallMethodDef CallEntries[] = { {"_TreeDist_cpp_matching_split_distance", (DL_FUNC) &_TreeDist_cpp_matching_split_distance, 3}, {"_TreeDist_cpp_jaccard_similarity", (DL_FUNC) &_TreeDist_cpp_jaccard_similarity, 5}, {"_TreeDist_cpp_msi_distance", (DL_FUNC) &_TreeDist_cpp_msi_distance, 3}, + {"_TreeDist_cpp_all_pairs_mci", (DL_FUNC) &_TreeDist_cpp_all_pairs_mci, 2}, {"_TreeDist_cpp_mutual_clustering", (DL_FUNC) &_TreeDist_cpp_mutual_clustering, 3}, {"_TreeDist_cpp_shared_phylo", (DL_FUNC) &_TreeDist_cpp_shared_phylo, 3}, {NULL, NULL, 0} diff --git a/src/SplitList.cpp b/src/SplitList.cpp index 2657ef67..f672723c 100644 --- a/src/SplitList.cpp +++ b/src/SplitList.cpp @@ -1,7 +1,7 @@ #include using namespace Rcpp; -#include "tree_distances.h" // Before SplitList.h, for int16 #include +#include "tree_distances.h" // Before SplitList.h, for int16 #include "SplitList.h" #define INLASTBIN(n, size) ((size) - ((size) - ((n) % (size))) % (size)) @@ -52,3 +52,131 @@ SplitList::SplitList(RawMatrix x) { } } } + +SplitList::SplitList(RawMatrix x, int16 n_tips) { + n_splits = x.rows(); + const int16 + n_input_bins = x.cols(), + input_bins_per_bin = BIN_SIZE / R_BIN_SIZE + ; + + n_bins = (n_input_bins + R_BIN_SIZE - 1) / input_bins_per_bin; + + if (n_bins < 1) throw std::invalid_argument("No leaves present."); + if (n_splits < 0) throw std::invalid_argument("Negative number of splits!?"); + if (n_bins > MAX_BINS) { + throw std::length_error("This many leaves cannot be supported. " // # nocov + "Please contact the TreeDist maintainer if you " // # nocov + "need to use more!"); // # nocov + } + // Rcout << "\n=== SplitList::SplitList(x, n_tips) ===\n"; + + for (int16 split = 0; split != n_splits; split++) { + int16 last_bin = n_bins - 1; + const int16 raggedy_bins = INLASTBIN(n_input_bins, R_BIN_SIZE); + + // Rcout << n_input_bins << " bins in; " << raggedy_bins << " raggedy bins\n"; + assert(raggedy_bins <= BIN_SIZE / R_BIN_SIZE); + assert(raggedy_bins > 0); + state[split][last_bin] = INSUBBIN(last_bin, 0); + // Rcout << " State[" << split << "][" << last_bin << "] = " << state[split][last_bin] << ".\n"; + for (int16 input_bin = 1; input_bin != raggedy_bins; input_bin++) { + // Rcout << " state [" << split << "][" << last_bin + // << "] increased from " << state[split][last_bin]; + state[split][last_bin] += INBIN(input_bin, last_bin); + // Rcout << " by " << INBIN(input_bin, last_bin) << "\n"; + } + + const bool invert = state[split][last_bin] & splitbit(1); + if (invert) { + const int16 last_bin_tips = INLASTBIN(n_tips, BIN_SIZE); + // Rcout << "Last bin tips: " << last_bin_tips << ".\n"; + const splitbit last_mask = last_bin_tips == sizeof(splitbit) * 8 ? + splitbit(0) - 1 : + (splitbit(1) << INLASTBIN(n_tips, BIN_SIZE)) - 1; + // Rcout << "Last mask: " << last_mask << ".\n"; + + state[split][last_bin] = ~state[split][last_bin] & last_mask; + } + + in_split[split] = count_bits(state[split][last_bin]); + + for (int16 bin = 0; bin != n_bins - 1; bin++) { + + // Rcout << "Split " << split << ", bin << " << bin << ".\n"; + state[split][bin] = INSUBBIN(bin, 0); + + for (int16 input_bin = 1; input_bin != input_bins_per_bin; input_bin++) { + // Rcout << " Adding " << INBIN(input_bin, bin) << " = " + // << (splitbit (x(split, (bin * input_bins_per_bin) + input_bin))) << " << " + // << (R_BIN_SIZE * input_bin) << " to state [" << split << "][" + // << bin << "], was " << state[split][bin] << "\n"; + + state[split][bin] += INBIN(input_bin, bin); + } + + if (invert) { + state[split][bin] = ~state[split][bin]; + } + + in_split[split] += count_bits(state[split][bin]); + } + } +} + +void SplitList::swap(int16 a, int16 b) { + for (int16 bin = n_bins; bin--; ) { + const splitbit tmp = state[a][bin]; + state[a][bin] = state[b][bin]; + state[b][bin] = tmp; + } + const int16 tmp = in_split[a]; + in_split[a] = in_split[b]; + in_split[b] = tmp; +} + +bool SplitList::less_than(int16 a, int16 b) { + for (int16 bin = 0; bin != n_bins; ++bin) { + if (state[a][bin] < state[b][bin]) { + return true; + } + } + return false; +} + +bool SplitList::greater_than(int16 a, int16 b) { + for (int16 bin = 0; bin != n_bins; ++bin) { + if (state[a][bin] > state[b][bin]) { + return true; + } + } + return false; +} + +int16 SplitList::partition(int16 lo, int16 hi) { + const int16 pivot = (hi + lo) / 2; + for (int16 i = lo - 1, j = hi + 1; ; ) { + do { + ++i; + } while (less_than(i, pivot)); + do { + --j; + } while (greater_than(j, pivot)); + if (i >= j) { + return j; + } + swap(i, j); + } +} + +void SplitList::quicksort(int16 lo, int16 hi) { + if (lo < hi) { + const int16 p = partition(lo, hi); + quicksort(lo, p); + quicksort(p + 1, hi); + } +} + +void SplitList::quicksort() { + quicksort(0, n_splits - 1); +} diff --git a/src/SplitList.h b/src/SplitList.h index 771d7531..42fb7a6a 100644 --- a/src/SplitList.h +++ b/src/SplitList.h @@ -3,16 +3,26 @@ #include #include +#include "ints.h" +#include "tree_distances.h" + using namespace Rcpp; const int16 R_BIN_SIZE = 8; class SplitList { + int16 partition(int16 lo, int16 hi); + void swap(int16 a, int16 b); + bool less_than(int16 a, int16 b); + bool greater_than(int16 a, int16 b); public: int16 n_splits, n_bins; int16 in_split[MAX_SPLITS]; splitbit state[MAX_SPLITS][MAX_BINS]; - SplitList(RawMatrix); + SplitList(RawMatrix x); + SplitList(RawMatrix x, int16 n_tips); + void quicksort(); + void quicksort(int16 lo, int16 hi); }; #endif diff --git a/src/SplitRoster.cpp b/src/SplitRoster.cpp new file mode 100644 index 00000000..671b4ebb --- /dev/null +++ b/src/SplitRoster.cpp @@ -0,0 +1,300 @@ +#include +#include "ints.h" +#include "tree_distances.h" +#include "SplitList.h" +#include "SplitRoster.h" +using namespace Rcpp; + +SplitRoster::SplitRoster(const List x, const IntegerVector nTip) { + n_tips = nTip[0]; + n_trees = x.length(); + splits.reserve(n_trees); + + for (int32 i = 0; i != n_trees; ++i) { + const RawMatrix xi = x[i]; + splits.emplace_back(xi, n_tips); + splits[i].quicksort(); + } + n_bins = splits[0].n_bins; + + const int16 max_splits = n_tips - 3; + + roster_tree = std::make_unique(n_trees * max_splits); + roster_split = std::make_unique(n_trees * max_splits); + in_split = std::make_unique(n_trees * max_splits); + out_split = std::make_unique(n_trees * max_splits); + roster_hits = std::make_unique(n_trees * max_splits); + index.reserve(n_trees); + + // Populate roster using k-way merge with tournament tree + const int32 + tournament_games = n_trees - 1, + tournament_nodes = n_trees + tournament_games + ; + + auto winners = std::make_unique(tournament_nodes); + auto losers = std::make_unique(tournament_nodes); + auto which_split = std::make_unique(n_trees); + + + // Populate leaves of tournament tree + for (int32 i = 0; i != n_trees; ++i) { + which_split[i] = splits[i].n_splits - 1; + winners[i + tournament_games] = i; + } + for (int32 i = n_trees + tournament_games; i != tournament_nodes; ++i) { + which_split[i] = -1; + winners[i] = -1; + } + + // Initial games + for (int32 i = tournament_games; i--; ) { + play_game(&i, winners, losers, which_split); + } + + roster_len = 0; + roster_tree[0] = winners[0]; + roster_split[0] = which_split[winners[0]]; + in_split[0] = splits[winners[0]].in_split[roster_split[0]]; + out_split[0] = n_tips - in_split[0]; + roster_hits[0] = 1; + index[roster_tree[0]][roster_split[0]] = 0; + + for (; ; ) { + const int32 winner = winners[0]; + --(which_split[winner]); + + int32 i = winner + tournament_games; + do { + i = (i - 1) / 2; + play_game(&i, winners, losers, which_split); + } while (i); + + if (which_split[winners[0]] < 0) break; + + push(winners[0], which_split); + } + + ++roster_len; // Point to first position after array + score.reserve(roster_len * roster_len); +} + +inline bool SplitRoster::splits_equal( + const splitbit (&a)[MAX_SPLITS][MAX_BINS], const int16 split_a, + const splitbit (&b)[MAX_SPLITS][MAX_BINS], const int16 split_b) { + for (int16 bin = n_bins; bin--; ) { + if (a[split_a][bin] != b[split_b][bin]) { + return false; + } + } + return true; +} + +inline bool SplitRoster::game_result( + const splitbit (&a)[MAX_SPLITS][MAX_BINS], const int16 split_a, + const splitbit (&b)[MAX_SPLITS][MAX_BINS], const int16 split_b) { + for (int16 bin = 0; bin != n_bins; ++bin) { + if (a[split_a][bin] > b[split_b][bin]) { + return true; + } + } + return false; +} + +inline void SplitRoster::play_game( + const int32 *node, + std::unique_ptr &winners, + std::unique_ptr &losers, + std::unique_ptr &which_split) { + const int32 + child1 = *node * 2 + 1, + child2 = child1 + 1, + tree1 = winners[child1], + tree2 = winners[child2] + ; + + bool child1_greater; + if (tree1 < 0 || which_split[tree1] < 0) { + child1_greater = false; + } else if (tree2 < 0 || which_split[tree2] < 0) { + child1_greater = true; + } else { + assert(tree1 >= 0); + assert(tree2 >= 0); + child1_greater = game_result( + splits[tree1].state, which_split[tree1], + splits[tree2].state, which_split[tree2]); + } + + if (child1_greater) { + winners[*node] = tree1; + losers[*node] = tree2; + } else { + winners[*node] = tree2; + losers[*node] = tree1; + } +} + +inline void SplitRoster::push( + const int32 tree, + std::unique_ptr &which_split) { + + const int16 new_split = which_split[tree]; + if (splits_equal(splits[tree].state, + new_split, + splits[roster_tree[roster_len]].state, + roster_split[roster_len])) { + ++roster_hits[roster_len]; + } else { + ++roster_len; + roster_tree[roster_len] = tree; + roster_split[roster_len] = new_split; + in_split[roster_len] = splits[tree].in_split[new_split]; + out_split[roster_len] = n_tips - in_split[roster_len]; + roster_hits[roster_len] = 1; + } + index[tree][new_split] = roster_len; +} + +#define SCORE(a, b) ((a) > (b) ? score[(b) * roster_len + (a)] : score[(a) * roster_len + (b)]) + +#define SPLIT(i) splits[roster_tree[(i)]].state[roster_split[(i)]] + +void SplitRoster::mutual_clustering() { + const cost max_score = BIG; + + for (int16 ai = 0; ai != roster_len; ++ai) { + const int16 + na = in_split[ai], + nA = out_split[ai] + ; + assert(na + nA == n_tips); + + SCORE(ai, ai) = max_score - + cost(max_score * (ic_matching(na, nA, n_tips) / n_tips)); + assert(SCORE(ai, ai) >= 0); + assert(SCORE(ai, ai) <= max_score); + + for (int16 bi = ai + 1; bi != roster_len; ++bi) { + + // x divides tips into a|A; y divides tips into b|B + int16 a_and_b = 0; + for (int16 bin = n_bins; bin--; ) { + a_and_b += count_bits(SPLIT(ai)[bin] & SPLIT(bi)[bin]); + } + assert(a_and_b < n_tips - 1); + + const int16 + nb = in_split[bi], + nB = out_split[bi], + a_and_B = na - a_and_b, + A_and_b = nb - a_and_b, + A_and_B = nA - A_and_b + ; + assert(nb + nB == n_tips); + assert(a_and_B >= 0); + assert(a_and_b >= 0); + assert(A_and_B >= 0); + assert(A_and_b >= 0); + + if (a_and_b == A_and_b && + a_and_b == a_and_B && + a_and_b == A_and_B) { + SCORE(ai, bi) = max_score; // Don't risk rounding error + } else { + SCORE(ai, bi) = max_score - + // Division by n_tips converts n(A&B) to P(A&B) for each ic_element + cost(max_score * (( + // 0 < Sum of IC_elements <= n_tips + ic_element(a_and_b, na, nb, n_tips) + + ic_element(a_and_B, na, nB, n_tips) + + ic_element(A_and_b, nA, nb, n_tips) + + ic_element(A_and_B, nA, nB, n_tips) + ) / n_tips) + ); + // Rcout << " Splits: \n " + // << SPLIT(ai)[0] <<", " << SPLIT(bi)[0]<<"\n"; + // Rcout << " Elements: \n" + // << " ___|_" << nb << "_|_" << nB << "_|_\n" + // << " " << na << " | " << a_and_b << " | " << a_and_B << " | \n" + // << " " << nA << " | " << A_and_b << " | " << A_and_B << " | " + // << " = " << n_tips << ";\n\n"; + } + assert(SCORE(ai, bi) >= 0); + assert(SCORE(ai, bi) <= max_score); + } + } +} + +NumericVector SplitRoster::score_pairs() { + + int32 + i = 0, + entry = 0, + n_scores = n_trees * (n_trees - 1) / 2; + ; + NumericVector ret(n_scores); + const cost max_score = BIG; + const int16 max_splits = n_tips - 3; + cost** lap_score = new cost*[max_splits]; + for (int16 li = max_splits; li--; ) lap_score[li] = new cost[max_splits]; + lap_col *rowsol = new lap_col[max_splits]; + lap_row *colsol = new lap_row[max_splits]; + cost *u = new cost[max_splits], *v = new cost[max_splits]; + + + for (i = 0; i != n_trees - 1; ++i) { + const int16 i_splits = splits[i].n_splits; + // Calculate tree's similarity to self + // for (int16 sp = 0; sp != i_splits; ++sp) { + // const int32 sp_i = index[i][sp]; + // ret[entry] += SCORE(sp_i, sp_i); + // } + // entry++; + for (int32 j = i + 1; j != n_trees; ++j) { + const int16 + j_splits = splits[j].n_splits, + most_splits = i_splits > j_splits ? i_splits : j_splits + ; + + + // TODO declare lap_score once at maximum size n_tips - 3 and retain. + cost** lap_score = new cost*[most_splits]; + for (int16 li = most_splits; li--; ) lap_score[li] = new cost[most_splits]; + lap_col *rowsol = new lap_col[most_splits]; + lap_row *colsol = new lap_row[most_splits]; + cost *u = new cost[most_splits], *v = new cost[most_splits]; + + + + for (int16 sp_i = i_splits; sp_i--; ) { + for (int16 sp_j = j_splits; sp_j--; ) { + lap_score[sp_i][sp_j] = SCORE(index[i][sp_i], index[j][sp_j]); + } + } + + for (int16 sp_i = i_splits; sp_i < most_splits; ++sp_i) { + for (int16 sp_j = 0; sp_j != most_splits; ++sp_j) { + lap_score[sp_i][sp_j] = max_score; + } + } + for (int16 sp_j = j_splits; sp_j < most_splits; ++sp_j) { + for (int16 sp_i = 0; sp_i != most_splits; ++sp_i) { + lap_score[sp_i][sp_j] = max_score; + } + } + + ret[entry] = double( + (max_score * most_splits) - + lap(most_splits, lap_score, rowsol, colsol, u, v) + ) / max_score; + + + entry++; + } + } + for (int16 li = max_splits; li--; ) delete[] lap_score[li]; + delete[] colsol; delete[] u; delete[] v; delete[] lap_score; + delete[] rowsol; + return ret; +} diff --git a/src/SplitRoster.h b/src/SplitRoster.h new file mode 100644 index 00000000..1bf3106e --- /dev/null +++ b/src/SplitRoster.h @@ -0,0 +1,51 @@ +#ifndef _TREEDIST_SPLITROSTER_HPP +#define _TREEDIST_SPLITROSTER_HPP + +#include +#include +#include +#include +#include "ints.h" +#include "tree_distances.h" +#include "SplitList.h" + +using namespace Rcpp; + +class SplitRoster { + int16 n_bins; + int32 roster_len; + + std::vector splits; + std::vector score; + std::unique_ptr roster_tree; + std::unique_ptr roster_split; + std::unique_ptr in_split; + std::unique_ptr out_split; + std::unique_ptr roster_hits; + std::vector > index; + + bool splits_equal( + const splitbit (&a)[MAX_SPLITS][MAX_BINS], const int16 split_a, + const splitbit (&b)[MAX_SPLITS][MAX_BINS], const int16 split_b); + bool game_result( + const splitbit (&a)[MAX_SPLITS][MAX_BINS], const int16 split_a, + const splitbit (&b)[MAX_SPLITS][MAX_BINS], const int16 split_b); + void play_game( + const int32 *node, + std::unique_ptr &winners, + std::unique_ptr &losers, + std::unique_ptr &which_split); + void push( + const int32 tree, + std::unique_ptr &which_split); + +public: + int16 n_tips, n_splits; + int32 n_trees; + + SplitRoster(const List x, const IntegerVector nTip); + void mutual_clustering(); + NumericVector score_pairs(); +}; + +#endif diff --git a/src/information.h b/src/information.h index 1e169e54..4bda1bdb 100644 --- a/src/information.h +++ b/src/information.h @@ -1,5 +1,5 @@ -#ifndef _TREEDIST_INFO_H -#define _TREEDIST_INFO_H +#ifndef _TREEDIST_INFO_HPP +#define _TREEDIST_INFO_HPP #include /* for log2() */ #include "ints.h" /* for int16 */ @@ -72,4 +72,4 @@ inline double split_clust_info (const int16 n_in, const int16 *n_tip, ); } -#endif \ No newline at end of file +#endif diff --git a/src/ints.h b/src/ints.h index ea7c6e62..878eee60 100644 --- a/src/ints.h +++ b/src/ints.h @@ -1,5 +1,5 @@ -#ifndef _TREEDIST_INT_H -#define _TREEDIST_INT_H +#ifndef _TREEDIST_INT_HPP +#define _TREEDIST_INT_HPP #include diff --git a/src/timsort.h b/src/timsort.h new file mode 100644 index 00000000..c8cf438d --- /dev/null +++ b/src/timsort.h @@ -0,0 +1,811 @@ +/* + * C++ implementation of timsort + * + * ported from Python's and OpenJDK's: + * - http://svn.python.org/projects/python/trunk/Objects/listobject.c + * - http://cr.openjdk.java.net/~martin/webrevs/openjdk7/timsort/raw_files/new/src/share/classes/java/util/TimSort.java + * + * Copyright (c) 2011 Fuji, Goro (gfx) . + * Copyright (c) 2019-2021 Morwenn. + * Copyright (c) 2021 Igor Kushnir . + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to + * deal in the Software without restriction, including without limitation the + * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or + * sell copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS + * IN THE SOFTWARE. + */ + +#ifndef GFX_TIMSORT_HPP +#define GFX_TIMSORT_HPP + +#include +#include +#include +#include +#include + +// Semantic versioning macros + +#define GFX_TIMSORT_VERSION_MAJOR 2 +#define GFX_TIMSORT_VERSION_MINOR 1 +#define GFX_TIMSORT_VERSION_PATCH 0 + +// Diagnostic selection macros + +#if defined(GFX_TIMSORT_ENABLE_ASSERT) || defined(GFX_TIMSORT_ENABLE_AUDIT) +# include +#endif + +#ifdef GFX_TIMSORT_ENABLE_ASSERT +# define GFX_TIMSORT_ASSERT(expr) assert(expr) +#else +# define GFX_TIMSORT_ASSERT(expr) ((void)0) +#endif + +#ifdef GFX_TIMSORT_ENABLE_AUDIT +# define GFX_TIMSORT_AUDIT(expr) assert(expr) +#else +# define GFX_TIMSORT_AUDIT(expr) ((void)0) +#endif + +#ifdef GFX_TIMSORT_ENABLE_LOG +# include +# define GFX_TIMSORT_LOG(expr) (std::clog << "# " << __func__ << ": " << expr << std::endl) +#else +# define GFX_TIMSORT_LOG(expr) ((void)0) +#endif + + +namespace gfx { + +// --------------------------------------- +// Implementation details +// --------------------------------------- + +namespace detail { + +// Equivalent to C++20 std::identity +struct identity { + template + constexpr T&& operator()(T&& value) const noexcept + { + return std::forward(value); + } +}; + +// Merge a predicate and a projection function +template +struct projection_compare { + projection_compare(Compare comp, Projection proj) : compare(comp), projection(proj) { + } + + template + bool operator()(T &&lhs, U &&rhs) { +#ifdef __cpp_lib_invoke + return static_cast(std::invoke(compare, + std::invoke(projection, std::forward(lhs)), + std::invoke(projection, std::forward(rhs)) + )); +#else + return static_cast(compare( + projection(std::forward(lhs)), + projection(std::forward(rhs)) + )); +#endif + } + + Compare compare; + Projection projection; +}; + +template struct run { + typedef typename std::iterator_traits::difference_type diff_t; + + Iterator base; + diff_t len; + + run(Iterator b, diff_t l) : base(b), len(l) { + } +}; + +template class TimSort { + typedef RandomAccessIterator iter_t; + typedef typename std::iterator_traits::value_type value_t; + typedef typename std::iterator_traits::reference ref_t; + typedef typename std::iterator_traits::difference_type diff_t; + + static const int MIN_MERGE = 32; + static const int MIN_GALLOP = 7; + + int minGallop_; // default to MIN_GALLOP + + std::vector tmp_; // temp storage for merges + typedef typename std::vector::iterator tmp_iter_t; + + std::vector > pending_; + + static void binarySort(iter_t const lo, iter_t const hi, iter_t start, Compare compare) { + GFX_TIMSORT_ASSERT(lo <= start); + GFX_TIMSORT_ASSERT(start <= hi); + if (start == lo) { + ++start; + } + for (; start < hi; ++start) { + GFX_TIMSORT_ASSERT(lo <= start); + value_t pivot = std::move(*start); + + iter_t const pos = std::upper_bound(lo, start, pivot, compare); + for (iter_t p = start; p > pos; --p) { + *p = std::move(*(p - 1)); + } + *pos = std::move(pivot); + } + } + + static diff_t countRunAndMakeAscending(iter_t const lo, iter_t const hi, Compare compare) { + GFX_TIMSORT_ASSERT(lo < hi); + + iter_t runHi = lo + 1; + if (runHi == hi) { + return 1; + } + + if (compare(*runHi, *lo)) { // decreasing + do { + ++runHi; + } while (runHi < hi && compare(*runHi, *(runHi - 1))); + std::reverse(lo, runHi); + } else { // non-decreasing + do { + ++runHi; + } while (runHi < hi && !compare(*runHi, *(runHi - 1))); + } + + return runHi - lo; + } + + static diff_t minRunLength(diff_t n) { + GFX_TIMSORT_ASSERT(n >= 0); + + diff_t r = 0; + while (n >= 2 * MIN_MERGE) { + r |= (n & 1); + n >>= 1; + } + return n + r; + } + + TimSort() : minGallop_(MIN_GALLOP) { + } + + // Silence GCC -Winline warning + ~TimSort() {} + + void pushRun(iter_t const runBase, diff_t const runLen) { + pending_.push_back(run(runBase, runLen)); + } + + void mergeCollapse(Compare compare) { + while (pending_.size() > 1) { + diff_t n = pending_.size() - 2; + + if ((n > 0 && pending_[n - 1].len <= pending_[n].len + pending_[n + 1].len) || + (n > 1 && pending_[n - 2].len <= pending_[n - 1].len + pending_[n].len)) { + if (pending_[n - 1].len < pending_[n + 1].len) { + --n; + } + mergeAt(n, compare); + } else if (pending_[n].len <= pending_[n + 1].len) { + mergeAt(n, compare); + } else { + break; + } + } + } + + void mergeForceCollapse(Compare compare) { + while (pending_.size() > 1) { + diff_t n = pending_.size() - 2; + + if (n > 0 && pending_[n - 1].len < pending_[n + 1].len) { + --n; + } + mergeAt(n, compare); + } + } + + void mergeAt(diff_t const i, Compare compare) { + diff_t const stackSize = pending_.size(); + GFX_TIMSORT_ASSERT(stackSize >= 2); + GFX_TIMSORT_ASSERT(i >= 0); + GFX_TIMSORT_ASSERT(i == stackSize - 2 || i == stackSize - 3); + + iter_t base1 = pending_[i].base; + diff_t len1 = pending_[i].len; + iter_t base2 = pending_[i + 1].base; + diff_t len2 = pending_[i + 1].len; + + pending_[i].len = len1 + len2; + + if (i == stackSize - 3) { + pending_[i + 1] = pending_[i + 2]; + } + + pending_.pop_back(); + + mergeConsecutiveRuns(base1, len1, base2, len2, std::move(compare)); + } + + void mergeConsecutiveRuns(iter_t base1, diff_t len1, iter_t base2, diff_t len2, Compare compare) { + GFX_TIMSORT_ASSERT(len1 > 0); + GFX_TIMSORT_ASSERT(len2 > 0); + GFX_TIMSORT_ASSERT(base1 + len1 == base2); + + diff_t const k = gallopRight(*base2, base1, len1, 0, compare); + GFX_TIMSORT_ASSERT(k >= 0); + + base1 += k; + len1 -= k; + + if (len1 == 0) { + return; + } + + len2 = gallopLeft(*(base1 + (len1 - 1)), base2, len2, len2 - 1, compare); + GFX_TIMSORT_ASSERT(len2 >= 0); + if (len2 == 0) { + return; + } + + if (len1 <= len2) { + mergeLo(base1, len1, base2, len2, compare); + } else { + mergeHi(base1, len1, base2, len2, compare); + } + } + + template + static diff_t gallopLeft(ref_t key, Iter const base, diff_t const len, + diff_t const hint, Compare compare) { + GFX_TIMSORT_ASSERT(len > 0); + GFX_TIMSORT_ASSERT(hint >= 0); + GFX_TIMSORT_ASSERT(hint < len); + + diff_t lastOfs = 0; + diff_t ofs = 1; + + if (compare(*(base + hint), key)) { + diff_t const maxOfs = len - hint; + while (ofs < maxOfs && compare(*(base + (hint + ofs)), key)) { + lastOfs = ofs; + ofs = (ofs << 1) + 1; + + if (ofs <= 0) { // int overflow + ofs = maxOfs; + } + } + if (ofs > maxOfs) { + ofs = maxOfs; + } + + lastOfs += hint; + ofs += hint; + } else { + diff_t const maxOfs = hint + 1; + while (ofs < maxOfs && !compare(*(base + (hint - ofs)), key)) { + lastOfs = ofs; + ofs = (ofs << 1) + 1; + + if (ofs <= 0) { + ofs = maxOfs; + } + } + if (ofs > maxOfs) { + ofs = maxOfs; + } + + diff_t const tmp = lastOfs; + lastOfs = hint - ofs; + ofs = hint - tmp; + } + GFX_TIMSORT_ASSERT(-1 <= lastOfs); + GFX_TIMSORT_ASSERT(lastOfs < ofs); + GFX_TIMSORT_ASSERT(ofs <= len); + + return std::lower_bound(base + (lastOfs + 1), base + ofs, key, compare) - base; + } + + template + static diff_t gallopRight(ref_t key, Iter const base, diff_t const len, + diff_t const hint, Compare compare) { + GFX_TIMSORT_ASSERT(len > 0); + GFX_TIMSORT_ASSERT(hint >= 0); + GFX_TIMSORT_ASSERT(hint < len); + + diff_t ofs = 1; + diff_t lastOfs = 0; + + if (compare(key, *(base + hint))) { + diff_t const maxOfs = hint + 1; + while (ofs < maxOfs && compare(key, *(base + (hint - ofs)))) { + lastOfs = ofs; + ofs = (ofs << 1) + 1; + + if (ofs <= 0) { + ofs = maxOfs; + } + } + if (ofs > maxOfs) { + ofs = maxOfs; + } + + diff_t const tmp = lastOfs; + lastOfs = hint - ofs; + ofs = hint - tmp; + } else { + diff_t const maxOfs = len - hint; + while (ofs < maxOfs && !compare(key, *(base + (hint + ofs)))) { + lastOfs = ofs; + ofs = (ofs << 1) + 1; + + if (ofs <= 0) { // int overflow + ofs = maxOfs; + } + } + if (ofs > maxOfs) { + ofs = maxOfs; + } + + lastOfs += hint; + ofs += hint; + } + GFX_TIMSORT_ASSERT(-1 <= lastOfs); + GFX_TIMSORT_ASSERT(lastOfs < ofs); + GFX_TIMSORT_ASSERT(ofs <= len); + + return std::upper_bound(base + (lastOfs + 1), base + ofs, key, compare) - base; + } + + static void rotateLeft(iter_t first, iter_t last) + { + value_t tmp = std::move(*first); + iter_t last_1 = std::move(first + 1, last, first); + *last_1 = std::move(tmp); + } + + static void rotateRight(iter_t first, iter_t last) + { + iter_t last_1 = last - 1; + value_t tmp = std::move(*last_1); + std::move_backward(first, last_1, last); + *first = std::move(tmp); + } + + + void mergeLo(iter_t const base1, diff_t len1, iter_t const base2, diff_t len2, Compare compare) { + GFX_TIMSORT_ASSERT(len1 > 0); + GFX_TIMSORT_ASSERT(len2 > 0); + GFX_TIMSORT_ASSERT(base1 + len1 == base2); + + if (len1 == 1) { + return rotateLeft(base1, base2 + len2); + } + if (len2 == 1) { + return rotateRight(base1, base2 + len2); + } + + copy_to_tmp(base1, len1); + + tmp_iter_t cursor1 = tmp_.begin(); + iter_t cursor2 = base2; + iter_t dest = base1; + + *dest = std::move(*cursor2); + ++cursor2; + ++dest; + --len2; + + int minGallop(minGallop_); + + // outer: + while (true) { + diff_t count1 = 0; + diff_t count2 = 0; + + do { + GFX_TIMSORT_ASSERT(len1 > 1); + GFX_TIMSORT_ASSERT(len2 > 0); + + if (compare(*cursor2, *cursor1)) { + *dest = std::move(*cursor2); + ++cursor2; + ++dest; + ++count2; + count1 = 0; + if (--len2 == 0) { + goto epilogue; + } + } else { + *dest = std::move(*cursor1); + ++cursor1; + ++dest; + ++count1; + count2 = 0; + if (--len1 == 1) { + goto epilogue; + } + } + } while ((count1 | count2) < minGallop); + + do { + GFX_TIMSORT_ASSERT(len1 > 1); + GFX_TIMSORT_ASSERT(len2 > 0); + + count1 = gallopRight(*cursor2, cursor1, len1, 0, compare); + if (count1 != 0) { + std::move_backward(cursor1, cursor1 + count1, dest + count1); + dest += count1; + cursor1 += count1; + len1 -= count1; + + if (len1 <= 1) { + goto epilogue; + } + } + *dest = std::move(*cursor2); + ++cursor2; + ++dest; + if (--len2 == 0) { + goto epilogue; + } + + count2 = gallopLeft(*cursor1, cursor2, len2, 0, compare); + if (count2 != 0) { + std::move(cursor2, cursor2 + count2, dest); + dest += count2; + cursor2 += count2; + len2 -= count2; + if (len2 == 0) { + goto epilogue; + } + } + *dest = std::move(*cursor1); + ++cursor1; + ++dest; + if (--len1 == 1) { + goto epilogue; + } + + --minGallop; + } while ((count1 >= MIN_GALLOP) | (count2 >= MIN_GALLOP)); + + if (minGallop < 0) { + minGallop = 0; + } + minGallop += 2; + } // end of "outer" loop + + epilogue: // merge what is left from either cursor1 or cursor2 + + minGallop_ = (std::min)(minGallop, 1); + + if (len1 == 1) { + GFX_TIMSORT_ASSERT(len2 > 0); + std::move(cursor2, cursor2 + len2, dest); + *(dest + len2) = std::move(*cursor1); + } else { + GFX_TIMSORT_ASSERT(len1 != 0 && "Comparison function violates its general contract"); + GFX_TIMSORT_ASSERT(len2 == 0); + GFX_TIMSORT_ASSERT(len1 > 1); + std::move(cursor1, cursor1 + len1, dest); + } + } + + void mergeHi(iter_t const base1, diff_t len1, iter_t const base2, diff_t len2, Compare compare) { + GFX_TIMSORT_ASSERT(len1 > 0); + GFX_TIMSORT_ASSERT(len2 > 0); + GFX_TIMSORT_ASSERT(base1 + len1 == base2); + + if (len1 == 1) { + return rotateLeft(base1, base2 + len2); + } + if (len2 == 1) { + return rotateRight(base1, base2 + len2); + } + + copy_to_tmp(base2, len2); + + iter_t cursor1 = base1 + len1; + tmp_iter_t cursor2 = tmp_.begin() + (len2 - 1); + iter_t dest = base2 + (len2 - 1); + + *dest = std::move(*(--cursor1)); + --dest; + --len1; + + int minGallop(minGallop_); + + // outer: + while (true) { + diff_t count1 = 0; + diff_t count2 = 0; + + // The next loop is a hot path of the algorithm, so we decrement + // eagerly the cursor so that it always points directly to the value + // to compare, but we have to implement some trickier logic to make + // sure that it points to the next value again by the end of said loop + --cursor1; + + do { + GFX_TIMSORT_ASSERT(len1 > 0); + GFX_TIMSORT_ASSERT(len2 > 1); + + if (compare(*cursor2, *cursor1)) { + *dest = std::move(*cursor1); + --dest; + ++count1; + count2 = 0; + if (--len1 == 0) { + goto epilogue; + } + --cursor1; + } else { + *dest = std::move(*cursor2); + --cursor2; + --dest; + ++count2; + count1 = 0; + if (--len2 == 1) { + ++cursor1; // See comment before the loop + goto epilogue; + } + } + } while ((count1 | count2) < minGallop); + ++cursor1; // See comment before the loop + + do { + GFX_TIMSORT_ASSERT(len1 > 0); + GFX_TIMSORT_ASSERT(len2 > 1); + + count1 = len1 - gallopRight(*cursor2, base1, len1, len1 - 1, compare); + if (count1 != 0) { + dest -= count1; + cursor1 -= count1; + len1 -= count1; + std::move_backward(cursor1, cursor1 + count1, dest + (1 + count1)); + + if (len1 == 0) { + goto epilogue; + } + } + *dest = std::move(*cursor2); + --cursor2; + --dest; + if (--len2 == 1) { + goto epilogue; + } + + count2 = len2 - gallopLeft(*(cursor1 - 1), tmp_.begin(), len2, len2 - 1, compare); + if (count2 != 0) { + dest -= count2; + cursor2 -= count2; + len2 -= count2; + std::move(cursor2 + 1, cursor2 + (1 + count2), dest + 1); + if (len2 <= 1) { + goto epilogue; + } + } + *dest = std::move(*(--cursor1)); + --dest; + if (--len1 == 0) { + goto epilogue; + } + + --minGallop; + } while ((count1 >= MIN_GALLOP) | (count2 >= MIN_GALLOP)); + + if (minGallop < 0) { + minGallop = 0; + } + minGallop += 2; + } // end of "outer" loop + + epilogue: // merge what is left from either cursor1 or cursor2 + + minGallop_ = (std::min)(minGallop, 1); + + if (len2 == 1) { + GFX_TIMSORT_ASSERT(len1 > 0); + dest -= len1; + std::move_backward(cursor1 - len1, cursor1, dest + (1 + len1)); + *dest = std::move(*cursor2); + } else { + GFX_TIMSORT_ASSERT(len2 != 0 && "Comparison function violates its general contract"); + GFX_TIMSORT_ASSERT(len1 == 0); + GFX_TIMSORT_ASSERT(len2 > 1); + std::move(tmp_.begin(), tmp_.begin() + len2, dest - (len2 - 1)); + } + } + + void copy_to_tmp(iter_t const begin, diff_t len) { + tmp_.assign(std::make_move_iterator(begin), + std::make_move_iterator(begin + len)); + } + +public: + + static void merge(iter_t const lo, iter_t const mid, iter_t const hi, Compare compare) { + GFX_TIMSORT_ASSERT(lo <= mid); + GFX_TIMSORT_ASSERT(mid <= hi); + + if (lo == mid || mid == hi) { + return; // nothing to do + } + + TimSort ts; + ts.mergeConsecutiveRuns(lo, mid - lo, mid, hi - mid, std::move(compare)); + + GFX_TIMSORT_LOG("1st size: " << (mid - lo) << "; 2nd size: " << (hi - mid) + << "; tmp_.size(): " << ts.tmp_.size()); + } + + static void sort(iter_t const lo, iter_t const hi, Compare compare) { + GFX_TIMSORT_ASSERT(lo <= hi); + + diff_t nRemaining = (hi - lo); + if (nRemaining < 2) { + return; // nothing to do + } + + if (nRemaining < MIN_MERGE) { + diff_t const initRunLen = countRunAndMakeAscending(lo, hi, compare); + GFX_TIMSORT_LOG("initRunLen: " << initRunLen); + binarySort(lo, hi, lo + initRunLen, compare); + return; + } + + TimSort ts; + diff_t const minRun = minRunLength(nRemaining); + iter_t cur = lo; + do { + diff_t runLen = countRunAndMakeAscending(cur, hi, compare); + + if (runLen < minRun) { + diff_t const force = (std::min)(nRemaining, minRun); + binarySort(cur, cur + force, cur + runLen, compare); + runLen = force; + } + + ts.pushRun(cur, runLen); + ts.mergeCollapse(compare); + + cur += runLen; + nRemaining -= runLen; + } while (nRemaining != 0); + + GFX_TIMSORT_ASSERT(cur == hi); + ts.mergeForceCollapse(compare); + GFX_TIMSORT_ASSERT(ts.pending_.size() == 1); + + GFX_TIMSORT_LOG("size: " << (hi - lo) << " tmp_.size(): " << ts.tmp_.size() + << " pending_.size(): " << ts.pending_.size()); + } +}; + +} // namespace detail + + +// --------------------------------------- +// Public interface implementation +// --------------------------------------- + +/** + * Stably merges two consecutive sorted ranges [first, middle) and [middle, last) into one + * sorted range [first, last) with a comparison function and a projection function. + */ +template +void timmerge(RandomAccessIterator first, RandomAccessIterator middle, + RandomAccessIterator last, Compare compare, Projection projection) { + typedef detail::projection_compare compare_t; + compare_t comp(std::move(compare), std::move(projection)); + GFX_TIMSORT_AUDIT(std::is_sorted(first, middle, comp) && "Precondition"); + GFX_TIMSORT_AUDIT(std::is_sorted(middle, last, comp) && "Precondition"); + detail::TimSort::merge(first, middle, last, comp); + GFX_TIMSORT_AUDIT(std::is_sorted(first, last, comp) && "Postcondition"); +} + +/** + * Same as std::inplace_merge(first, middle, last, compare). + */ +template +void timmerge(RandomAccessIterator first, RandomAccessIterator middle, + RandomAccessIterator last, Compare compare) { + gfx::timmerge(first, middle, last, compare, detail::identity()); +} + +/** + * Same as std::inplace_merge(first, middle, last). + */ +template +void timmerge(RandomAccessIterator first, RandomAccessIterator middle, + RandomAccessIterator last) { + typedef typename std::iterator_traits::value_type value_type; + gfx::timmerge(first, middle, last, std::less(), detail::identity()); +} + +/** + * Stably sorts a range with a comparison function and a projection function. + */ +template +void timsort(RandomAccessIterator const first, RandomAccessIterator const last, + Compare compare, Projection projection) { + typedef detail::projection_compare compare_t; + compare_t comp(std::move(compare), std::move(projection)); + detail::TimSort::sort(first, last, comp); + GFX_TIMSORT_AUDIT(std::is_sorted(first, last, comp) && "Postcondition"); +} + +/** + * Same as std::stable_sort(first, last, compare). + */ +template +void timsort(RandomAccessIterator const first, RandomAccessIterator const last, Compare compare) { + gfx::timsort(first, last, compare, detail::identity()); +} + +/** + * Same as std::stable_sort(first, last). + */ +template +void timsort(RandomAccessIterator const first, RandomAccessIterator const last) { + typedef typename std::iterator_traits::value_type value_type; + gfx::timsort(first, last, std::less(), detail::identity()); +} + +/** + * Stably sorts a range with a comparison function and a projection function. + */ +template +void timsort(RandomAccessRange &range, Compare compare, Projection projection) { + gfx::timsort(std::begin(range), std::end(range), compare, projection); +} + +/** + * Same as std::stable_sort(std::begin(range), std::end(range), compare). + */ +template +void timsort(RandomAccessRange &range, Compare compare) { + gfx::timsort(std::begin(range), std::end(range), compare); +} + +/** + * Same as std::stable_sort(std::begin(range), std::end(range)). + */ +template +void timsort(RandomAccessRange &range) { + gfx::timsort(std::begin(range), std::end(range)); +} + +} // namespace gfx + +#undef GFX_TIMSORT_ENABLE_ASSERT +#undef GFX_TIMSORT_ASSERT +#undef GFX_TIMSORT_ENABLE_AUDIT +#undef GFX_TIMSORT_AUDIT +#undef GFX_TIMSORT_ENABLE_LOG +#undef GFX_TIMSORT_LOG + +#endif // GFX_TIMSORT_HPP \ No newline at end of file diff --git a/src/tree_distances.cpp b/src/tree_distances.cpp index 19767fac..f9da15e4 100644 --- a/src/tree_distances.cpp +++ b/src/tree_distances.cpp @@ -3,6 +3,7 @@ #include #include "tree_distances.h" #include "SplitList.h" +#include "SplitRoster.h" using namespace Rcpp; @@ -380,6 +381,13 @@ List cpp_msi_distance (const RawMatrix x, const RawMatrix y, } +// [[Rcpp::export]] +NumericVector cpp_all_pairs_mci(const List x, const IntegerVector nTip) { + SplitRoster roster(x, nTip); + roster.mutual_clustering(); + return roster.score_pairs(); +} + // [[Rcpp::export]] List cpp_mutual_clustering (const RawMatrix x, const RawMatrix y, const IntegerVector nTip) { @@ -492,7 +500,7 @@ List cpp_mutual_clustering (const RawMatrix x, const RawMatrix y, } } - const double lap_score = + const double lap_score = double((max_score * lap_dim) - lap(lap_dim, score, rowsol, colsol, u, v)) / max_score; NumericVector final_score = @@ -635,4 +643,4 @@ List cpp_shared_phylo (const RawMatrix x, const RawMatrix y, _["matching"] = final_matching); } - \ No newline at end of file + diff --git a/src/tree_distances.h b/src/tree_distances.h index 543c848c..4210d7a9 100644 --- a/src/tree_distances.h +++ b/src/tree_distances.h @@ -1,3 +1,6 @@ +#ifndef _TREEDIST_TREE_DISTANCES_H +#define _TREEDIST_TREE_DISTANCES_H + #include #include #include "ints.h" @@ -14,18 +17,18 @@ typedef int16 lap_col; /*************** CONSTANTS *******************/ -const int16 BIN_SIZE = 64, - MAX_BINS = 32, - MAX_TIPS = BIN_SIZE * MAX_BINS, - MAX_SPLITS = MAX_TIPS; /* -3, but quicker if a power of two? */ +#define BIN_SIZE 64 +#define MAX_BINS 32 +#define MAX_TIPS (BIN_SIZE * MAX_BINS) +#define MAX_SPLITS MAX_TIPS /* -3, but quicker if a power of two? */ -const splitbit ALL_ONES = (std::numeric_limits::max)(); +#define ALL_ONES splitbit((std::numeric_limits::max)()) /* For a reason I've not estabilshed, shrinking BIG is necessary to avoid * an infinite loop in lap. */ -const cost BIG = ((std::numeric_limits::max)() / MAX_SPLITS); +#define BIG cost((std::numeric_limits::max)() / MAX_SPLITS) -const splitbit right16bits = 65535U; +#define right16bits splitbit(65535U) const splitbit powers_of_two[64] = { 0x1, 0x2, 0x4, 0x8, 0x10, 0x20, 0x40, 0x80, @@ -45,7 +48,7 @@ const splitbit powers_of_two[64] = { 0x1000000000000000, 0x2000000000000000, 0x4000000000000000, 0x8000000000000000 }; -const cost ROUND_PRECISION = 2048 * 2048; +#define ROUND_PRECISION cost(2048 * 2048) /***** Constants requiring initialization *****/ @@ -77,3 +80,5 @@ extern double extern List cpp_robinson_foulds_distance (RawMatrix x, RawMatrix y, IntegerVector nTip); + +#endif