From 7597f19aa847653bd743bc528bc010d06aec7db6 Mon Sep 17 00:00:00 2001 From: smehringer Date: Mon, 21 Jul 2025 12:30:43 +0200 Subject: [PATCH 01/35] [FEATURE] Add the LSH fast layout algorithm. --- .gitignore | 1 + include/chopper/configuration.hpp | 7 + .../chopper/layout/determine_split_bins.hpp | 29 + include/chopper/layout/execute.hpp | 4 +- include/chopper/layout/fast_layout.hpp | 50 ++ .../chopper/layout/fast_layout_cluster.hpp | 133 +++ .../fast_layout_find_bins_to_be_split.hpp | 89 +++ .../chopper/layout/partition_user_bins.hpp | 59 ++ src/chopper_layout.cpp | 19 +- src/layout/CMakeLists.txt | 12 +- src/layout/determine_split_bins.cpp | 153 ++++ src/layout/execute.cpp | 43 +- src/layout/fast_layout.cpp | 457 +++++++++++ src/layout/partition_user_bins.cpp | 755 ++++++++++++++++++ src/set_up_parser.cpp | 6 + test/api/layout/CMakeLists.txt | 4 + test/api/layout/determine_split_bins_test.cpp | 58 ++ test/api/layout/execute_layout_test.cpp | 398 ++++++++- .../layout/execute_with_estimation_test.cpp | 396 ++++++++- test/api/layout/fast_layout_cluster_test.cpp | 72 ++ ...fast_layout_find_bins_to_be_split_test.cpp | 66 ++ test/api/layout/hibf_statistics_test.cpp | 70 +- test/api/layout/partition_user_bins_test.cpp | 178 +++++ test/cli/cli_chopper_pipeline_test.cpp | 105 +++ test/cli/cli_output_sketches.cpp | 2 +- 25 files changed, 3124 insertions(+), 42 deletions(-) create mode 100644 include/chopper/layout/determine_split_bins.hpp create mode 100644 include/chopper/layout/fast_layout.hpp create mode 100644 include/chopper/layout/fast_layout_cluster.hpp create mode 100644 include/chopper/layout/fast_layout_find_bins_to_be_split.hpp create mode 100644 include/chopper/layout/partition_user_bins.hpp create mode 100644 src/layout/determine_split_bins.cpp create mode 100644 src/layout/fast_layout.cpp create mode 100644 src/layout/partition_user_bins.cpp create mode 100644 test/api/layout/determine_split_bins_test.cpp create mode 100644 test/api/layout/fast_layout_cluster_test.cpp create mode 100644 test/api/layout/fast_layout_find_bins_to_be_split_test.cpp create mode 100644 test/api/layout/partition_user_bins_test.cpp diff --git a/.gitignore b/.gitignore index 7dfea74f..90a24ee0 100644 --- a/.gitignore +++ b/.gitignore @@ -32,6 +32,7 @@ # build dir /build/ +/build_debug/ # editor configuration /.vscode/ diff --git a/include/chopper/configuration.hpp b/include/chopper/configuration.hpp index bae58feb..7a1bfad0 100644 --- a/include/chopper/configuration.hpp +++ b/include/chopper/configuration.hpp @@ -22,6 +22,9 @@ namespace chopper struct configuration { + //!\brief Whether to use the fast layout algorithm instead of the default one. + bool fast_layout{false}; + /*!\name General Configuration * \{ */ @@ -77,6 +80,10 @@ struct configuration mutable seqan::hibf::concurrent_timer union_estimation_timer{}; mutable seqan::hibf::concurrent_timer rearrangement_timer{}; mutable seqan::hibf::concurrent_timer dp_algorithm_timer{}; + mutable seqan::hibf::concurrent_timer lsh_algorithm_timer{}; + mutable seqan::hibf::concurrent_timer search_partition_algorithm_timer{}; + mutable seqan::hibf::concurrent_timer intital_partition_timer{}; + mutable seqan::hibf::concurrent_timer small_layouts_timer{}; void read_from(std::istream & stream); diff --git a/include/chopper/layout/determine_split_bins.hpp b/include/chopper/layout/determine_split_bins.hpp new file mode 100644 index 00000000..82b5ac0a --- /dev/null +++ b/include/chopper/layout/determine_split_bins.hpp @@ -0,0 +1,29 @@ +// -------------------------------------------------------------------------------------------------- +// Copyright (c) 2006-2023, Knut Reinert & Freie Universität Berlin +// Copyright (c) 2016-2023, Knut Reinert & MPI für molekulare Genetik +// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License +// shipped with this file and also available at: https://github.com/seqan/chopper/blob/main/LICENSE.md +// -------------------------------------------------------------------------------------------------- + +/*!\file + * \brief Provides chopper::determine_split_bins. + * \author Svenja Mehringer + */ + +#pragma once + +#include + +#include + +namespace chopper::layout +{ + +std::pair determine_split_bins(chopper::configuration const & config, + std::vector const & positions, + std::vector const & cardinalities, + size_t const num_technical_bins, + size_t const num_user_bins, + std::vector> & partitions); + +} diff --git a/include/chopper/layout/execute.hpp b/include/chopper/layout/execute.hpp index 52bfc3b7..00c4fa2e 100644 --- a/include/chopper/layout/execute.hpp +++ b/include/chopper/layout/execute.hpp @@ -13,12 +13,14 @@ #include #include +#include namespace chopper::layout { int execute(chopper::configuration & config, std::vector> const & filenames, - std::vector const & sketches); + std::vector const & sketches, + std::vector const & minHash_sketches); } // namespace chopper::layout diff --git a/include/chopper/layout/fast_layout.hpp b/include/chopper/layout/fast_layout.hpp new file mode 100644 index 00000000..3541771d --- /dev/null +++ b/include/chopper/layout/fast_layout.hpp @@ -0,0 +1,50 @@ +// --------------------------------------------------------------------------------------------------- +// Copyright (c) 2006-2023, Knut Reinert & Freie Universität Berlin +// Copyright (c) 2016-2023, Knut Reinert & MPI für molekulare Genetik +// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License +// shipped with this file and also available at: https://github.com/seqan/chopper/blob/main/LICENSE.md +// --------------------------------------------------------------------------------------------------- + +#pragma once + +#include + +#include + +#include +#include +#include + +namespace chopper::layout +{ + +/*!\brief Computes an HIBF layout with the fast (LSH/similarity-based) layout algorithm. + * \param[in] config The configuration. Uses `hibf_config` (`tmax`, `number_of_user_bins`, FPR settings) + * and the fast-layout timers. + * \param[in] positions The global indices of all user bins. Must be a permutation of + * `[0, config.hibf_config.number_of_user_bins)`. + * \param[in] cardinalities The cardinality of each user bin, indexed by global user bin index. + * \param[in] sketches The HyperLogLog sketch of each user bin, indexed by global user bin index. + * \param[in] minHash_sketches The MinHash tables of each user bin, indexed by global user bin index. + * \param[out] hibf_layout The resulting layout. Expected to be empty on entry. `user_bins` is resized to + * `number_of_user_bins`, so `user_bins[i].idx == i`. + * + * 1. **Top level:** partition_user_bins distributes all user bins onto `tmax` technical bins. The user bins are + * initialised in `hibf_layout`: merged bins get `previous_TB_indices = {t}`; split and single bins get + * `storage_TB_id` and their number of consecutive technical bins. `top_level_max_bin_id` is set to the technical + * bin with the largest FPR-corrected size. + * 2. **Lower levels:** Merged bins are processed in parallel (OpenMP `taskloop`). Depending on + * do_I_need_a_fast_layout, each is laid out recursively with the fast layout or with the regular DP layout, + * and the result is grafted into `hibf_layout`. + * 3. `max_bins` is sorted by level (path length), then lexicographically by path. + * + * \throws std::logic_error In debug builds, if the layout's user bins are not a permutation of `positions`. + */ +void fast_layout(chopper::configuration const & config, + std::vector const & positions, + std::vector const & cardinalities, + std::vector const & sketches, + std::vector const & minHash_sketches, + seqan::hibf::layout::layout & hibf_layout); + +} // namespace chopper::layout diff --git a/include/chopper/layout/fast_layout_cluster.hpp b/include/chopper/layout/fast_layout_cluster.hpp new file mode 100644 index 00000000..f02a14e9 --- /dev/null +++ b/include/chopper/layout/fast_layout_cluster.hpp @@ -0,0 +1,133 @@ +// --------------------------------------------------------------------------------------------------- +// Copyright (c) 2006-2023, Knut Reinert & Freie Universität Berlin +// Copyright (c) 2016-2023, Knut Reinert & MPI für molekulare Genetik +// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License +// shipped with this file and also available at: https://github.com/seqan/chopper/blob/main/LICENSE.md +// --------------------------------------------------------------------------------------------------- + +#pragma once + +#include +#include +#include + +namespace chopper::layout +{ + +/*\brief foo + */ +struct Cluster +{ +protected: + size_t representative_id{}; // representative id of the cluster; identifier; + + std::vector user_bins{}; // the user bins contained in thus cluster + + std::optional moved_id{std::nullopt}; // where this Clusters user bins where moved to + +public: + Cluster() = default; + Cluster(Cluster const &) = default; + Cluster(Cluster &&) = default; + Cluster & operator=(Cluster const &) = default; + Cluster & operator=(Cluster &&) = default; + ~Cluster() = default; + + Cluster(size_t const id, size_t const user_bins_id) : representative_id{id}, user_bins({user_bins_id}) + {} + + Cluster(size_t const id) : Cluster{id, id} + {} + + size_t id() const + { + return representative_id; + } + + std::vector const & contained_user_bins() const + { + return user_bins; + } + + bool has_been_moved() const + { + return moved_id.has_value(); + } + + bool empty() const + { + return user_bins.empty(); + } + + size_t size() const + { + return user_bins.size(); + } + + size_t pop_back() + { + size_t last = user_bins.back(); + user_bins.pop_back(); + return last; + } + + void add_user_bin(size_t user_bin) + { + user_bins.push_back(user_bin); + } + + bool is_valid(size_t id) const + { + bool const ids_equal = representative_id == id; + bool const properly_moved = has_been_moved() && empty(); + bool const not_moved = !has_been_moved() && !empty(); + + return ids_equal && (properly_moved || not_moved); + } + + size_t moved_to_cluster_id() const + { + assert(moved_id.has_value()); + assert(is_valid(representative_id)); + return moved_id.value(); + } + + void move_to(Cluster & target_cluster) + { + target_cluster.user_bins.insert(target_cluster.user_bins.end(), this->user_bins.begin(), this->user_bins.end()); + this->user_bins.clear(); + moved_id = target_cluster.id(); + } + + void sort_by_cardinality(std::vector const & cardinalities) + { + std::ranges::sort(user_bins, + [&cardinalities](auto const & v1, auto const & v2) + { + return cardinalities[v1] > cardinalities[v2]; + }); + } +}; + +// A valid cluster is one that hasn't been moved but actually contains user bins +// A valid cluster at position i is identified by the following equality: cluster[i].size() >= 1 && cluster[i][0] == i +// A moved cluster is one that has been joined and thereby moved to another cluster +// A moved cluster i is identified by the following: cluster[i].size() == 1 && cluster[i][0] != i +// returns position of the representative cluster +inline size_t LSH_find_representative_cluster(std::vector const & clusters, size_t current_id) +{ + std::reference_wrapper representative = clusters[current_id]; + + assert(representative.get().is_valid(current_id)); + + while (representative.get().has_been_moved()) + { + current_id = representative.get().moved_to_cluster_id(); + representative = clusters[current_id]; // replace by next cluster + assert(representative.get().is_valid(current_id)); + } + + return current_id; +} + +} // namespace chopper::layout \ No newline at end of file diff --git a/include/chopper/layout/fast_layout_find_bins_to_be_split.hpp b/include/chopper/layout/fast_layout_find_bins_to_be_split.hpp new file mode 100644 index 00000000..6891a544 --- /dev/null +++ b/include/chopper/layout/fast_layout_find_bins_to_be_split.hpp @@ -0,0 +1,89 @@ +// -------------------------------------------------------------------------------------------------- +// Copyright (c) 2006-2023, Knut Reinert & Freie Universität Berlin +// Copyright (c) 2016-2023, Knut Reinert & MPI für molekulare Genetik +// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License +// shipped with this file and also available at: https://github.com/seqan/chopper/blob/main/LICENSE.md +// -------------------------------------------------------------------------------------------------- + +/*!\file + * \brief Provides chopper::layout::find_bins_to_be_split. + * \author Svenja Mehringer + */ + +#pragma once + +#include +#include +#include +#include +#include + +namespace chopper::layout +{ + +/*!\brief Determines which user bins are split and how many technical bins the split user bins occupy. + * \param[in] sorted_positions User bin indices, sorted by descending cardinality + * (see seqan::hibf::sketch::toolbox::sort_by_cardinalities). + * \param[in] cardinalities The cardinality of each user bin, indexed by the values in `sorted_positions`. + * \param[in] threshold The initial cardinality threshold. User bins with a cardinality greater than + * `threshold` are split. Must be greater than 0. + * \param[in] max_bins The maximum number of technical bins that split user bins may occupy. + * \returns A pair of + * 1. the number of user bins to split, i.e., the length of the prefix of `sorted_positions`, and + * 2. the number of technical bins for these split user bins, which is at most `max_bins`. + * + * The user bins with a cardinality greater than `threshold` form a prefix of `sorted_positions`. The number of + * technical bins they need is `sum / threshold` (rounded down), where `sum` is the sum of their cardinalities. + * If this exceeds `max_bins`, `threshold` is multiplied by `max(1.01, sum / (threshold * max_bins))` and the prefix + * is computed again, until the split user bins fit into `max_bins` technical bins. + * + * The number of split user bins is clamped to the number of technical bins, so each split user bin gets at least + * one technical bin. + */ +inline std::pair find_bins_to_be_split(std::vector const & sorted_positions, + std::vector const & cardinalities, + size_t threshold, + size_t const max_bins) +{ + assert(threshold > 0); + assert(max_bins > 0); + + size_t idx{0}; + size_t sum{0}; + // update idx and sum + auto find_idx_and_sum = [&]() + { + while (idx < sorted_positions.size() && cardinalities[sorted_positions[idx]] > threshold) + { + sum += cardinalities[sorted_positions[idx]]; + ++idx; + } + }; + find_idx_and_sum(); + + // SInce the threshold is the expected size of each technical bin, the number of split bins needed i: + size_t number_of_split_bins = sum / threshold; + + // If number_of_split_bins is more than the available technical bins, the threshold must be adjusted + while (number_of_split_bins > max_bins) + { + // Since number_of_split_bins is too large iwth the given threshold + // we need to increase the threshold such that less user bins are targetted for splitting. + // The new threshold is therefore adjusted to the new expected average technical bin size of + // distributing the split bin content evenly to `max_bin` technical bins. + // max(threshold, ...) to avoid an endless loop. + threshold = std::max(threshold + 1, static_cast(sum) / max_bins); + + // update idx, sum and number_of_split_bins + idx = 0; + sum = 0; + find_idx_and_sum(); + number_of_split_bins = sum / threshold; + } + + assert(idx <= number_of_split_bins); // there should never be more user bins than available split bins + + return {idx, number_of_split_bins}; +} + +} // namespace chopper::layout diff --git a/include/chopper/layout/partition_user_bins.hpp b/include/chopper/layout/partition_user_bins.hpp new file mode 100644 index 00000000..0eb9b258 --- /dev/null +++ b/include/chopper/layout/partition_user_bins.hpp @@ -0,0 +1,59 @@ +// -------------------------------------------------------------------------------------------------- +// Copyright (c) 2006-2023, Knut Reinert & Freie Universität Berlin +// Copyright (c) 2016-2023, Knut Reinert & MPI für molekulare Genetik +// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License +// shipped with this file and also available at: https://github.com/seqan/chopper/blob/main/LICENSE.md +// -------------------------------------------------------------------------------------------------- + +/*!\file + * \brief Provides chopper::layout::partition_user_bins. + * \author Svenja Mehringer + */ + +#pragma once + +#include +#include + +#include + +#include +#include + +namespace chopper::layout +{ + +/*!\brief Distributes user bins onto `tmax` technical bins (partitions) of a single IBF. + * \param[in] config The configuration. Uses `tmax`, `sketch_bits`, `maximum_fpr`, `relaxed_fpr` and + * `number_of_hash_functions` of `config.hibf_config`, and the layout timers. + * \param[in] positions The indices of the user bins to distribute (need not be sorted). + * \param[in] cardinalities The cardinality of each user bin, indexed by the values in `positions`. + * \param[in] sketches The HyperLogLog sketch of each user bin, indexed by the values in `positions`. + * \param[in] minHash_sketches The MinHash tables of each user bin, indexed by the values in `positions`. + * Each table needs at least 3 sketches of at least 5 hashes each. + * \param[out] partitions Must have size `tmax` on entry. On return, `partitions[i]` holds the user bin + * indices assigned to technical bin `i`. + * + * User bins are sorted by descending cardinality and handled in two groups: + * + * 1. **Split bins**: The largest user bins (cardinality above a threshold, see find_bins_to_be_split) are + * distributed over the technical bins at the *back* of `partitions` by determine_split_bins. At least one + * technical bin is left for merged bins. If fewer user bins remain than technical bins would be left, more + * technical bins are given to the split bins. + * 2. **Merged bins**: The remaining user bins go into the technical bins at the *front* of `partitions`. They are + * clustered by similarity with MinHash LSH and assigned by HyperLogLog union estimates (see lsh_sim_approach). + * The target size per merged technical bin is `max(max_split_size / relaxed_fpr_correction, split_threshold)`. + * + * The first `split_threshold` is `ceil(joint_estimate / tmax)`, a lower bound that holds only if all user bins are + * identical. If there are merged bins, the threshold is recalibrated **once**, based on the ratio of + * `max_merged_size * relaxed_fpr_correction` to `max_split_size` (averaged with `max_merged_size` if nothing was + * split). `partitions` is then cleared and both steps run again with the new threshold. + */ +void partition_user_bins(chopper::configuration const & config, + std::vector const & positions, + std::vector const & cardinalities, + std::vector const & sketches, + std::vector const & minHash_sketches, + std::vector> & partitions); + +} // namespace chopper::layout diff --git a/src/chopper_layout.cpp b/src/chopper_layout.cpp index 5fda4ddb..9df893a3 100644 --- a/src/chopper_layout.cpp +++ b/src/chopper_layout.cpp @@ -93,6 +93,7 @@ int chopper_layout(chopper::configuration & config, sharg::parser & parser) std::vector> filenames{}; std::vector sketches{}; + std::vector minHash_sketches{}; if (input_is_a_sketch_file) { @@ -106,6 +107,7 @@ int chopper_layout(chopper::configuration & config, sharg::parser & parser) filenames = std::move(sin.filenames); // No need to call check_filenames because the files are not read. sketches = std::move(sin.hll_sketches); + minHash_sketches = std::move(sin.minHash_sketches); validate_configuration(parser, config, sin.chopper_config); } else @@ -128,17 +130,18 @@ int chopper_layout(chopper::configuration & config, sharg::parser & parser) if (!input_is_a_sketch_file) { config.compute_sketches_timer.start(); - seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches); + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); config.compute_sketches_timer.stop(); } - exit_code |= chopper::layout::execute(config, filenames, sketches); + exit_code |= chopper::layout::execute(config, filenames, sketches, minHash_sketches); if (!config.disable_sketch_output) { chopper::sketch::sketch_file sout{.chopper_config = config, .filenames = std::move(filenames), - .hll_sketches = std::move(sketches)}; + .hll_sketches = std::move(sketches), + .minHash_sketches = std::move(minHash_sketches)}; std::ofstream os{config.sketch_directory, std::ios::binary}; cereal::BinaryOutputArchive oarchive{os}; oarchive(sout); @@ -151,11 +154,19 @@ int chopper_layout(chopper::configuration & config, sharg::parser & parser) output_stream << "sketching_in_seconds\t" << "layouting_in_seconds\t" << "union_estimation_in_seconds\t" - << "rearrangement_in_seconds\n"; + << "rearrangement_in_seconds\t" + << "lsh_in_seconds\t" + << "intital_partition_timer_in_seconds\t" + << "small_layouts_timer_in_seconds\t" + << "search_best_p_in_seconds\n"; output_stream << config.compute_sketches_timer.in_seconds() << '\t'; output_stream << config.dp_algorithm_timer.in_seconds() << '\t'; output_stream << config.union_estimation_timer.in_seconds() << '\t'; output_stream << config.rearrangement_timer.in_seconds() << '\t'; + output_stream << config.lsh_algorithm_timer.in_seconds() << '\t'; + output_stream << config.intital_partition_timer.in_seconds() << '\t'; + output_stream << config.small_layouts_timer.in_seconds() << '\t'; + output_stream << config.search_partition_algorithm_timer.in_seconds() << '\n'; } return exit_code; diff --git a/src/layout/CMakeLists.txt b/src/layout/CMakeLists.txt index 915c83d7..7d2f98b5 100644 --- a/src/layout/CMakeLists.txt +++ b/src/layout/CMakeLists.txt @@ -4,8 +4,16 @@ if (TARGET chopper::layout) return () endif () -add_library (chopper_layout STATIC determine_best_number_of_technical_bins.cpp execute.cpp hibf_statistics.cpp - ibf_query_cost.cpp input.cpp output.cpp +add_library (chopper_layout STATIC + determine_best_number_of_technical_bins.cpp + determine_split_bins.cpp + execute.cpp + fast_layout.cpp + hibf_statistics.cpp + ibf_query_cost.cpp + input.cpp + output.cpp + partition_user_bins.cpp ) target_link_libraries (chopper_layout PUBLIC chopper::shared) add_library (chopper::layout ALIAS chopper_layout) diff --git a/src/layout/determine_split_bins.cpp b/src/layout/determine_split_bins.cpp new file mode 100644 index 00000000..fb4ba653 --- /dev/null +++ b/src/layout/determine_split_bins.cpp @@ -0,0 +1,153 @@ +// --------------------------------------------------------------------------------------------------- +// Copyright (c) 2006-2023, Knut Reinert & Freie Universität Berlin +// Copyright (c) 2016-2023, Knut Reinert & MPI für molekulare Genetik +// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License +// shipped with this file and also available at: https://github.com/seqan/chopper/blob/main/LICENSE.md +// --------------------------------------------------------------------------------------------------- + +#include // for allocator, string +#include // for vector + +#include +#include + +#include +#include // for data_store +#include + +namespace chopper::layout +{ + +std::pair determine_split_bins(chopper::configuration const & config, + std::vector const & positions, + std::vector const & cardinalities, + size_t const num_technical_bins, + size_t const num_user_bins, + std::vector> & partitions) +{ + if (num_user_bins == 0) + return {0, 0}; + + assert(num_technical_bins > 0u); + assert(num_user_bins > 0u); + assert(num_user_bins <= num_technical_bins); + + auto const fpr_correction = + seqan::hibf::layout::compute_fpr_correction({.fpr = config.hibf_config.maximum_fpr, // + .hash_count = config.hibf_config.number_of_hash_functions, + .t_max = num_technical_bins}); + + std::vector> matrix(num_technical_bins); // rows + for (auto & v : matrix) + v.resize(num_user_bins, std::numeric_limits::max()); // columns + + std::vector> trace(num_technical_bins); // rows + for (auto & v : trace) + v.resize(num_user_bins, std::numeric_limits::max()); // columns + + size_t const extra_bins = num_technical_bins - num_user_bins + 1; + + // initialize first column (first row is initialized with inf) + double const ub_cardinality = static_cast(cardinalities[positions[0]]); + for (size_t i = 0; i < extra_bins; ++i) + { + size_t const corrected_ub_cardinality = static_cast(ub_cardinality * fpr_correction[i + 1]); + matrix[i][0] = seqan::hibf::divide_and_ceil(corrected_ub_cardinality, i + 1u); + } + + // we must iterate column wise + for (size_t j = 1; j < num_user_bins; ++j) + { + double const ub_cardinality = static_cast(cardinalities[positions[j]]); + + for (size_t i = j; i < j + extra_bins; ++i) + { + size_t minimum{std::numeric_limits::max()}; + + for (size_t i_prime = j - 1; i_prime < i; ++i_prime) + { + size_t const corrected_ub_cardinality = + static_cast(ub_cardinality * fpr_correction[(i - i_prime)]); + size_t score = std::max(seqan::hibf::divide_and_ceil(corrected_ub_cardinality, i - i_prime), + matrix[i_prime][j - 1]); + + // std::cout << "j:" << j << " i:" << i << " i':" << i_prime << " score:" << score << std::endl; + + minimum = (score < minimum) ? (trace[i][j] = i_prime, score) : minimum; + } + + matrix[i][j] = minimum; + } + } + + // seqan::hibf::layout::print_matrix(matrix, num_technical_bins, num_user_bins, std::numeric_limits::max()); + //seqan::hibf::layout::print_matrix(trace, num_technical_bins, num_user_bins, std::numeric_limits::max()); + + // backtracking + // first, in the last column, find the row with minimum score (it can happen that the last rows are equally good) + // the less rows, the better + size_t trace_i{num_technical_bins - 1}; + size_t best_score{std::numeric_limits::max()}; + for (size_t best_i{0}; best_i < num_technical_bins; ++best_i) + { + if (matrix[best_i][num_user_bins - 1] < best_score) + { + best_score = matrix[best_i][num_user_bins - 1]; + trace_i = best_i; + } + } + size_t const number_of_split_tbs{trace_i + 1}; // trace_i is the position. So +1 for the number + + // now that we found the best trace_i start usual backtracking + size_t trace_j = num_user_bins - 1; + + // size_t max_id{}; + size_t max_size{}; + + size_t bin_id{}; + + while (trace_j > 0) + { + size_t next_i = trace[trace_i][trace_j]; + size_t const number_of_bins = (trace_i - next_i); + size_t const cardinality = cardinalities[positions[trace_j]]; + size_t const corrected_cardinality = static_cast(cardinality * fpr_correction[number_of_bins]); + size_t const cardinality_per_bin = seqan::hibf::divide_and_ceil(corrected_cardinality, number_of_bins); + + if (cardinality_per_bin > max_size) + { + // max_id = bin_id; + max_size = cardinality_per_bin; + } + + for (size_t splits{0}; splits < number_of_bins; ++splits) + { + partitions[partitions.size() - 1 - bin_id].push_back(positions[trace_j]); + ++bin_id; + } + + trace_i = trace[trace_i][trace_j]; + --trace_j; + } + ++trace_i; // because we want the length not the index. Now trace_i == number_of_bins + size_t const cardinality = cardinalities[positions[0]]; + size_t const corrected_cardinality = static_cast(cardinality * fpr_correction[trace_i]); + // NOLINTNEXTLINE(clang-analyzer-core.DivideZero) + size_t const cardinality_per_bin = seqan::hibf::divide_and_ceil(corrected_cardinality, trace_i); + + if (cardinality_per_bin > max_size) + { + // max_id = bin_id; + max_size = cardinality_per_bin; + } + + for (size_t splits{0}; splits < trace_i; ++splits) + { + partitions[partitions.size() - 1 - bin_id].push_back(positions[0]); + ++bin_id; + } + + return {number_of_split_tbs, max_size}; +} + +} // namespace chopper::layout diff --git a/src/layout/execute.cpp b/src/layout/execute.cpp index 6630447e..8fd9278f 100644 --- a/src/layout/execute.cpp +++ b/src/layout/execute.cpp @@ -9,23 +9,19 @@ #include #include #include -#include #include #include #include #include -#include -#include #include #include #include +#include #include #include -#include #include -#include #include #include // for estimate_kmer_counts #include @@ -35,33 +31,48 @@ namespace chopper::layout int execute(chopper::configuration & config, std::vector> const & filenames, - std::vector const & sketches) + std::vector const & sketches, + std::vector const & minHash_sketches) { config.hibf_config.validate_and_set_defaults(); + std::vector cardinalities; + seqan::hibf::sketch::estimate_kmer_counts(sketches, cardinalities); + seqan::hibf::layout::layout hibf_layout; - std::vector kmer_counts; - seqan::hibf::sketch::estimate_kmer_counts(sketches, kmer_counts); if (config.determine_best_tmax) { - hibf_layout = determine_best_number_of_technical_bins(config, kmer_counts, sketches); + // ToDo, what about determine_best_tmax iwth fast layout? + hibf_layout = determine_best_number_of_technical_bins(config, cardinalities, sketches); } else { config.dp_algorithm_timer.start(); - hibf_layout = seqan::hibf::layout::compute_layout(config.hibf_config, - kmer_counts, - sketches, - seqan::hibf::iota_vector(sketches.size()), - config.union_estimation_timer, - config.rearrangement_timer); + if (config.fast_layout) + { + fast_layout(config, + seqan::hibf::iota_vector(sketches.size()), + cardinalities, + sketches, + minHash_sketches, + hibf_layout); + } + else + { + hibf_layout = seqan::hibf::layout::compute_layout(config.hibf_config, + cardinalities, + sketches, + seqan::hibf::iota_vector(sketches.size()), + config.union_estimation_timer, + config.rearrangement_timer); + } config.dp_algorithm_timer.stop(); if (config.output_verbose_statistics) { size_t dummy{}; - chopper::layout::hibf_statistics global_stats{config, sketches, kmer_counts}; + chopper::layout::hibf_statistics global_stats{config, sketches, cardinalities}; global_stats.hibf_layout = hibf_layout; global_stats.print_header_to(std::cout); global_stats.print_summary_to(dummy, std::cout); diff --git a/src/layout/fast_layout.cpp b/src/layout/fast_layout.cpp new file mode 100644 index 00000000..0d661433 --- /dev/null +++ b/src/layout/fast_layout.cpp @@ -0,0 +1,457 @@ +// --------------------------------------------------------------------------------------------------- +// Copyright (c) 2006-2023, Knut Reinert & Freie Universität Berlin +// Copyright (c) 2016-2023, Knut Reinert & MPI für molekulare Genetik +// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License +// shipped with this file and also available at: https://github.com/seqan/chopper/blob/main/LICENSE.md +// --------------------------------------------------------------------------------------------------- + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +namespace chopper::layout +{ + +/*!\brief Computes the layout of a subset of user bins with the regular (DP-based) HIBF layout algorithm. + * \param[in] config The configuration (uses `hibf_config`). + * \param[in] positions The global indices of the user bins to lay out. + * \param[in] cardinalities The cardinality of each user bin, indexed by global user bin index. + * \param[in] sketches The HyperLogLog sketch of each user bin, indexed by global user bin index. + * \returns A layout whose top level is the IBF of the given subset. `user_bins[i].idx` are global user bin indices; + * `previous_TB_indices` and `max_bins` are relative to this subset's top-level IBF. + * + * The union estimation and rearrangement timers are local and discarded. + */ +seqan::hibf::layout::layout general_layout(chopper::configuration const & config, + std::vector positions, + std::vector const & cardinalities, + std::vector const & sketches) +{ + seqan::hibf::layout::layout hibf_layout; + + seqan::hibf::concurrent_timer union_estimation_timer{}; + seqan::hibf::concurrent_timer rearrangement_timer{}; + seqan::hibf::concurrent_timer dp_algorithm_timer{}; + + dp_algorithm_timer.start(); + hibf_layout = seqan::hibf::layout::compute_layout(config.hibf_config, + cardinalities, + sketches, + std::move(positions), + union_estimation_timer, + rearrangement_timer); + dp_algorithm_timer.stop(); + + return hibf_layout; +} + +/*!\brief Decides whether the lower-level IBF for a merged bin is laid out by the fast layout or by general_layout. + * \param[in] config The configuration (uses `hibf_config.tmax`). + * \param[in] positions The global indices of the user bins in the merged bin. + * \param[in] cardinalities The cardinality of each user bin, indexed by global user bin index. + * \returns `true` if the fast layout should be used: + * - `false` if there are fewer than `64 * tmax` user bins. With few user bins per technical bin, the greedy + * fast layout merges only a handful of bins at a time, and the resulting heavy splitting on lower levels + * raises the FPR correction. + * - `true` if there are more than 500'000 user bins. The DP layout would take more than about half a day. + * - otherwise `true` only if no user bin is larger than `sum_of_cardinalities / tmax`, i.e., no user bin + * needs splitting and a merge-only layout suffices. + */ +bool do_I_need_a_fast_layout(chopper::configuration const & config, + std::vector const & positions, + std::vector const & cardinalities) +{ + // the fast layout heuristic would greedily merge even if merging only 2 bins at a time + // merging only little number of bins is highly disadvantegous for lower levels because few bins + // will be heavily split and this will raise the fpr correction for split bins + // Thus, if the average number of user bins per technical bin is less then 64, we should not fast layout + if (positions.size() < (64 * config.hibf_config.tmax)) + return false; + + if (positions.size() > 500'000) // layout takes more than half a day (should this be a user option?) + return true; + + size_t largest_size{0}; + size_t sum_of_cardinalities{0}; + + for (size_t const i : positions) + { + sum_of_cardinalities += cardinalities[i]; + largest_size = std::max(largest_size, cardinalities[i]); + } + + size_t const cardinality_per_tb = sum_of_cardinalities / config.hibf_config.tmax; + + bool const largest_user_bin_might_be_split = largest_size > cardinality_per_tb; + + // if no splitting is needed, its worth it to use a fast-merge-only algorithm + if (!largest_user_bin_might_be_split) + return true; + + return false; +} + +/*!\brief Records a lower-level IBF, computed by partition_user_bins, in the layout. + * \param[in] config The configuration (uses the FPR settings of `hibf_config`). + * \param[in,out] hibf_layout The global layout. Its `user_bins` must be indexed by global user bin index + * (`user_bins[i].idx == i`), as initialised by fast_layout. + * \param[in] partitions The technical bins of the new IBF: `partitions[t]` holds the global user bin indices + * in technical bin `t`. + * \param[in] sketches The HyperLogLog sketch of each user bin, indexed by global user bin index. + * \param[in] previous The merged-bin path identifying the new IBF. Every user bin in `partitions` must + * currently have exactly this path as `previous_TB_indices`. + * + * - **Merged bin** (more than one user bin): `t` is appended to each user bin's `previous_TB_indices`. Their + * `storage_TB_id` is set later, when the next level down is recorded. + * - **Split or single bin** (one user bin): `storage_TB_id` is set to `t`, and `number_of_technical_bins` to the + * number of consecutive technical bins holding the same user bin. partition_user_bins places split bins + * contiguously. + * - **Empty technical bins** are skipped. + * + * The technical bin with the largest FPR-corrected size (relaxed correction for merged bins, split correction for + * split bins) is appended to `hibf_layout.max_bins` as `(previous, max_bin_id)`. + * + * Not thread-safe; callers serialise it with `omp critical`. + */ +void add_level_to_layout(chopper::configuration const & config, + seqan::hibf::layout::layout & hibf_layout, + std::vector> const & partitions, + std::vector const & sketches, + std::vector const & previous) +{ + size_t max_bin_id{0}; + size_t max_size{0}; + + auto const split_fpr_correction = + seqan::hibf::layout::compute_fpr_correction({.fpr = config.hibf_config.maximum_fpr, // + .hash_count = config.hibf_config.number_of_hash_functions, + .t_max = partitions.size()}); + + double const relaxed_fpr_correction = seqan::hibf::layout::compute_relaxed_fpr_correction( + {.fpr = config.hibf_config.maximum_fpr, // + .relaxed_fpr = config.hibf_config.relaxed_fpr, + .hash_count = config.hibf_config.number_of_hash_functions}); + + // we assume here that the user bins have been sorted by user bin id such that pos = idx + for (size_t partition_idx{0}; partition_idx < partitions.size(); ++partition_idx) + { + auto const & partition = partitions[partition_idx]; + + if (partition.size() > 1) // merged bin + { + seqan::hibf::sketch::hyperloglog current_sketch{sketches[0]}; // ensure same bit size + current_sketch.reset(); + + for (size_t const user_bin_id : partition) + { + assert(hibf_layout.user_bins[user_bin_id].idx == user_bin_id); + auto & current_user_bin = hibf_layout.user_bins[user_bin_id]; + + // update + assert(previous == current_user_bin.previous_TB_indices); + current_user_bin.previous_TB_indices.push_back(partition_idx); + current_sketch.merge(sketches[user_bin_id]); + } + + // update max_bin_id, max_size + size_t const current_size = current_sketch.estimate() * relaxed_fpr_correction; + if (current_size > max_size) + { + max_bin_id = partition_idx; + max_size = current_size; + } + } + else if (partition.size() == 0) // should not happen.. dge case? + { + continue; + } + else // single or split bin (partition.size() == 1) + { + auto & current_user_bin = hibf_layout.user_bins[partitions[partition_idx][0]]; + assert(current_user_bin.idx == partitions[partition_idx][0]); + current_user_bin.storage_TB_id = partition_idx; + current_user_bin.number_of_technical_bins = 1; // initialise to 1 + + while (partition_idx + 1 < partitions.size() && partitions[partition_idx].size() == 1 + && partitions[partition_idx + 1].size() == 1 + && partitions[partition_idx][0] == partitions[partition_idx + 1][0]) + { + ++current_user_bin.number_of_technical_bins; + ++partition_idx; + } + + // update max_bin_id, max_size + size_t const current_size = sketches[current_user_bin.idx].estimate() + * split_fpr_correction[current_user_bin.number_of_technical_bins]; + if (current_size > max_size) + { + max_bin_id = current_user_bin.storage_TB_id; + max_size = current_size; + } + } + } + + hibf_layout.max_bins.emplace_back(previous, max_bin_id); // add lower level meta information +} + +/*!\brief Grafts a layout computed for a merged bin (see general_layout) into the global layout. + * \param[in,out] child_layout The layout of the merged bin's subtree. Its `max_bins` are modified (prefixed with + * `new_previous`) and copied. + * \param[in,out] hibf_layout The global layout. Its `user_bins` must be indexed by global user bin index. + * \param[in] new_previous The path of the merged bin that the child layout's top level is attached to. + * + * - The child's top-level IBF is added to `hibf_layout.max_bins` as `(new_previous, child.top_level_max_bin_id)`. + * - Every other child max bin is added with `new_previous` prepended to its path. + * - For every child user bin, the child's relative path is appended to the global user bin's `previous_TB_indices` + * (which already equals `new_previous`), and `storage_TB_id` and `number_of_technical_bins` are copied. + * + * Not thread-safe; callers serialise it with `omp critical`. + */ +void update_layout_from_child_layout(seqan::hibf::layout::layout & child_layout, + seqan::hibf::layout::layout & hibf_layout, + std::vector const & new_previous) +{ + hibf_layout.max_bins.emplace_back(new_previous, child_layout.top_level_max_bin_id); + + for (auto & max_bin : child_layout.max_bins) + { + max_bin.previous_TB_indices.insert(max_bin.previous_TB_indices.begin(), + new_previous.begin(), + new_previous.end()); + hibf_layout.max_bins.push_back(max_bin); + } + + for (auto const & user_bin : child_layout.user_bins) + { + auto & actual_user_bin = hibf_layout.user_bins[user_bin.idx]; + + actual_user_bin.previous_TB_indices.insert(actual_user_bin.previous_TB_indices.end(), + user_bin.previous_TB_indices.begin(), + user_bin.previous_TB_indices.end()); + actual_user_bin.number_of_technical_bins = user_bin.number_of_technical_bins; + actual_user_bin.storage_TB_id = user_bin.storage_TB_id; + } +} + +/*!\brief Lays out the lower-level IBF of one merged bin with the fast layout, and recursively lays out its children. + * \param[in] config The configuration. + * \param[in] positions The global indices of the user bins in the merged bin. + * \param[in] cardinalities The cardinality of each user bin, indexed by global user bin index. + * \param[in] sketches The HyperLogLog sketch of each user bin, indexed by global user bin index. + * \param[in] minHash_sketches The MinHash tables of each user bin, indexed by global user bin index. + * \param[in,out] hibf_layout The global layout, updated in place. + * \param[in] previous The merged-bin path of this IBF. + * + * Partitions `positions` into `tmax` technical bins (partition_user_bins) and records them (add_level_to_layout). + * Each resulting merged bin is then handled like in fast_layout: another fast-layout recursion or a general_layout, + * decided by do_I_need_a_fast_layout. The recursion runs sequentially within the calling OpenMP task. Only the + * updates to `hibf_layout` are in critical sections. + */ +void fast_layout_recursion(chopper::configuration const & config, + std::vector const & positions, + std::vector const & cardinalities, + std::vector const & sketches, + std::vector const & minHash_sketches, + seqan::hibf::layout::layout & hibf_layout, + std::vector const & previous) +{ + std::vector> tmax_partitions(config.hibf_config.tmax); + + // here we assume that we want to start with a fast layout + partition_user_bins(config, positions, cardinalities, sketches, minHash_sketches, tmax_partitions); + +#pragma omp critical + { + add_level_to_layout(config, hibf_layout, tmax_partitions, sketches, previous); + } + + for (size_t partition_idx = 0; partition_idx < tmax_partitions.size(); ++partition_idx) + { + auto const & partition = tmax_partitions[partition_idx]; + auto const new_previous = [&]() + { + auto cpy{previous}; + cpy.push_back(partition_idx); + return cpy; + }(); + + if (partition.empty() || partition.size() == 1) // nothing to merge + continue; + + if (do_I_need_a_fast_layout(config, partition, cardinalities)) + { + fast_layout_recursion(config, + partition, + cardinalities, + sketches, + minHash_sketches, + hibf_layout, + new_previous); // recurse fast_layout + } + else + { + auto child_layout = general_layout(config, partition, cardinalities, sketches); + +#pragma omp critical + { + update_layout_from_child_layout(child_layout, hibf_layout, new_previous); + } + } + } +} + +void fast_layout(chopper::configuration const & config, + std::vector const & positions, + std::vector const & cardinalities, + std::vector const & sketches, + std::vector const & minHash_sketches, + seqan::hibf::layout::layout & hibf_layout) +{ + std::vector> tmax_partitions(config.hibf_config.tmax); + + // here we assume that we want to start with a fast layout + config.intital_partition_timer.start(); + partition_user_bins(config, positions, cardinalities, sketches, minHash_sketches, tmax_partitions); + config.intital_partition_timer.stop(); + + auto const split_fpr_correction = + seqan::hibf::layout::compute_fpr_correction({.fpr = config.hibf_config.maximum_fpr, // + .hash_count = config.hibf_config.number_of_hash_functions, + .t_max = config.hibf_config.tmax}); + + double const relaxed_fpr_correction = seqan::hibf::layout::compute_relaxed_fpr_correction( + {.fpr = config.hibf_config.maximum_fpr, // + .relaxed_fpr = config.hibf_config.relaxed_fpr, + .hash_count = config.hibf_config.number_of_hash_functions}); + + size_t max_bin_id{0}; + size_t max_size{0}; + hibf_layout.user_bins.resize(config.hibf_config.number_of_user_bins); + + // initialise user bins in layout + for (size_t partition_idx = 0; partition_idx < tmax_partitions.size(); ++partition_idx) + { + if (tmax_partitions[partition_idx].size() > 1) // merged bin + { + seqan::hibf::sketch::hyperloglog current_sketch{sketches[0]}; // ensure same bit size + current_sketch.reset(); + + for (size_t const user_bin_id : tmax_partitions[partition_idx]) + { + hibf_layout.user_bins[user_bin_id] = {.previous_TB_indices = {partition_idx}, + .storage_TB_id = 0 /*not determiend yet*/, + .number_of_technical_bins = 1 /*not determiend yet*/, + .idx = user_bin_id}; + current_sketch.merge(sketches[user_bin_id]); + } + + // update max_bin_id, max_size + size_t const current_size = current_sketch.estimate() * relaxed_fpr_correction; + if (current_size > max_size) + { + max_bin_id = partition_idx; + max_size = current_size; + } + } + else if (tmax_partitions[partition_idx].size() == 0) // should not happen. Edge case? + { + continue; + } + else // single or split bin (tmax_partitions[partition_idx].size() == 1) + { + assert(tmax_partitions[partition_idx].size() == 1); + size_t const user_bin_id = tmax_partitions[partition_idx][0]; + hibf_layout.user_bins[user_bin_id] = {.previous_TB_indices = {}, + .storage_TB_id = partition_idx, + .number_of_technical_bins = 1 /*determiend below*/, + .idx = user_bin_id}; + + while (partition_idx + 1 < tmax_partitions.size() && tmax_partitions[partition_idx].size() == 1 + && tmax_partitions[partition_idx + 1].size() == 1 + && tmax_partitions[partition_idx][0] == tmax_partitions[partition_idx + 1][0]) + { + ++hibf_layout.user_bins[user_bin_id].number_of_technical_bins; + ++partition_idx; + } + + // update max_bin_id, max_size + size_t const current_size = + sketches[user_bin_id].estimate() + * split_fpr_correction[hibf_layout.user_bins[user_bin_id].number_of_technical_bins]; + if (current_size > max_size) + { + max_bin_id = hibf_layout.user_bins[user_bin_id].storage_TB_id; + max_size = current_size; + } + } + } + + hibf_layout.top_level_max_bin_id = max_bin_id; + + config.small_layouts_timer.start(); +#pragma omp parallel +#pragma omp single + { +#pragma omp taskloop + for (size_t partition_idx = 0; partition_idx < tmax_partitions.size(); ++partition_idx) + { + auto const & partition = tmax_partitions[partition_idx]; + + if (partition.empty() || partition.size() == 1) // nothing to merge + continue; + + if (do_I_need_a_fast_layout(config, partition, cardinalities)) + { + fast_layout_recursion(config, + partition, + cardinalities, + sketches, + minHash_sketches, + hibf_layout, + {partition_idx}); // recurse fast_layout + } + else + { + auto small_layout = general_layout(config, partition, cardinalities, sketches); + +#pragma omp critical + { + update_layout_from_child_layout(small_layout, hibf_layout, std::vector{partition_idx}); + } + } + } + } + config.small_layouts_timer.stop(); + + // sort records ascending by the number of bin indices (corresponds to the IBF levels) + // GCOVR_EXCL_START + std::ranges::sort(hibf_layout.max_bins, + [](auto const & r, auto const & l) + { + if (r.previous_TB_indices.size() == l.previous_TB_indices.size()) + return std::ranges::lexicographical_compare(r.previous_TB_indices, l.previous_TB_indices); + else + return r.previous_TB_indices.size() < l.previous_TB_indices.size(); + }); + // GCOVR_EXCL_STOP + +#ifndef NDEBUG + // sanity check in debug + std::vector layout_user_bins{}; + for (auto & user_bin : hibf_layout.user_bins) + layout_user_bins.push_back(user_bin.idx); + if (!std::ranges::is_permutation(layout_user_bins, positions)) + throw std::logic_error{"Not all/Wrong user bins have been assigned to the layout!"}; +#endif +} + +} // namespace chopper::layout \ No newline at end of file diff --git a/src/layout/partition_user_bins.cpp b/src/layout/partition_user_bins.cpp new file mode 100644 index 00000000..16bd7b19 --- /dev/null +++ b/src/layout/partition_user_bins.cpp @@ -0,0 +1,755 @@ +// -------------------------------------------------------------------------------------------------- +// Copyright (c) 2006-2023, Knut Reinert & Freie Universität Berlin +// Copyright (c) 2016-2023, Knut Reinert & MPI für molekulare Genetik +// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License +// shipped with this file and also available at: https://github.com/seqan/chopper/blob/main/LICENSE.md +// -------------------------------------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace chopper::layout +{ + +/*!\brief Combines the first `number_of_hashes_to_consider` MinHash values of a sketch into a single LSH key. + * \param[in] sketch A single MinHash sketch (one row of seqan::hibf::sketch::minhashes::table). + * \param[in] number_of_hashes_to_consider The number of leading hashes to combine (LSH parameter r). + * Must be `<= sketch.size()`. + * \returns The sum of the first `number_of_hashes_to_consider` hashes, with unsigned wrap-around. + * + * This is the AND step of the LSH AND-OR scheme: two user bins get the same key only if all r hashes agree, + * apart from sum collisions. + */ +uint64_t lsh_hash_the_sketch(std::vector const & sketch, size_t const number_of_hashes_to_consider) +{ + assert(number_of_hashes_to_consider <= sketch.size()); + return std::reduce(sketch.begin(), sketch.begin() + number_of_hashes_to_consider); +} + +/*!\brief Builds the LSH collision table of the current clusters for one LSH band. + * \param[in] clusters The current clusters. Clusters that were moved are skipped. + * \param[in] minHash_sketches The MinHash tables of all user bins, indexed by global user bin index. + * \param[in] current_sketch_index The sketch (row of the MinHash table) to use (LSH band index). + * \param[in] current_number_of_sketch_hashes The number of hashes combined per key (LSH parameter r). + * \returns A map from LSH key to the sorted, unique ids of the representative clusters that produced the key. + * + * Each user bin in a valid cluster adds its key, and the cluster's id is stored under that key. A multi-member + * cluster can therefore appear under several keys, so clusters that share a key with *any* member collide. + */ +auto LSH_fill_hashtable(std::vector const & clusters, + std::vector const & minHash_sketches, + size_t const current_sketch_index, + size_t const current_number_of_sketch_hashes) +{ + robin_hood::unordered_flat_map> table; + + [[maybe_unused]] size_t processed_user_bins{0}; // only for sanity check + + for (size_t pos = 0; pos < clusters.size(); ++pos) + { + auto const & current = clusters[pos]; + assert(current.is_valid(pos)); + + if (current.has_been_moved()) // cluster has been moved somewhere else, don't process + continue; + + for (size_t const user_bin_idx : current.contained_user_bins()) + { + ++processed_user_bins; + uint64_t const key = lsh_hash_the_sketch(minHash_sketches[user_bin_idx].table[current_sketch_index], + current_number_of_sketch_hashes); + table[key].push_back(current.id()); // insert representative for all user bins + } + } + assert(processed_user_bins == clusters.size()); // all user bins should've been processed by one of the clusters + + // uniquify list. Since I am inserting representative_idx's into the table, the same number can + // be inserted into multiple splots, and multiple times in the same slot. + for (auto & [key, list] : table) + { + std::ranges::sort(list); + auto const ret = std::ranges::unique(list); + list.erase(ret.begin(), ret.end()); + } + + return table; +} + +/*!\brief Clusters very similar user bins by iterative MinHash LSH. + * \param[in] minHash_sketches The MinHash tables of all user bins, indexed by global user bin index. + * \param[in] positions The global indices of the user bins to cluster. + * \param[in] cardinalities The cardinality of each user bin, indexed by global user bin index. + * \param[in] sketches The HyperLogLog sketch of each user bin, indexed by global user bin index. + * \param[in] average_technical_bin_size The threshold. Stop clustering as soon as a cluster's estimated cardinality + * reaches this. + * \param[in] config The configuration (uses `hibf_config.sketch_bits`). + * \returns One Cluster per entry in `positions`. Cluster `i` is created with id `i` (a local index) and holds user bin + * `positions[i]` (a global index). After merging, a cluster is either valid (`id() == i`, contains >= 1 + * user bins) or moved (empty, and `moved_to_cluster_id()` points to the cluster it was merged into, + * possibly through a chain of moves; resolve with LSH_find_representative_cluster). + * + * In each round, one LSH band (`current_sketch_index`) is used to build a collision table (LSH_fill_hashtable), and + * all clusters in a bucket are merged into the representative of the first one (OR step). The merged cluster's + * cardinality is re-estimated from the union of its HyperLogLog sketches. Rounds stop after + * `number_of_max_minHash_sketches` (b = 3) bands or when the largest cluster cardinality reaches + * `average_technical_bin_size`. The number of hashes per key stays at `minHash_sketch_size` (r = 5) in all rounds. + * b and r were chosen by experiment. + */ +std::vector very_similar_LSH_clustering(std::vector const & minHash_sketches, + std::vector const & positions, + std::vector const & cardinalities, + std::vector const & sketches, + size_t const average_technical_bin_size, + chopper::configuration const & config) +{ + // The following two parameters are experimentally derived: + size_t const number_of_max_minHash_sketches{3}; // LSH ADD+OR parameter b + size_t const minHash_sketch_size{5}; // LSH ADD+OR parameter r + size_t const number_of_user_bins{positions.size()}; + seqan::hibf::sketch::hyperloglog const empty_sketch{config.hibf_config.sketch_bits}; + + assert(!minHash_sketches.empty()); + assert(minHash_sketches[0].table.size() >= number_of_max_minHash_sketches); + assert(minHash_sketches[0].table[0].size() >= minHash_sketch_size); + assert(number_of_user_bins <= minHash_sketches.size()); + + // initialise clusters with a signle user bin per cluster. + // clusters are either + // 1) of size 1; containing an id != position where the id points to the cluster it has been moved to + // e.g. cluster[Y] = {Z} (Y has been moved into Z, so Z could look likes this cluster[Z] = {Z, Y}) + // 2) of size >= 1; with the first entry beging id == position (a valid cluster) + // e.g. cluster[X] = {X} // valid singleton + // e.g. cluster[X] = {X, a, b, c, ...} // valid cluster with more joined entries + // The clusters could me moved recursively, s.t. + // cluster[A] = {B} + // cluster[B] = {C} + // cluster[C] = {C, A, B} // is valid cluster since cluster[C][0] == C; contains A and B + std::vector clusters; + clusters.reserve(number_of_user_bins); + + std::vector current_cluster_cardinality(number_of_user_bins); + std::vector current_cluster_sketches(number_of_user_bins, empty_sketch); + size_t current_max_cluster_size{0}; + size_t current_sketch_index{0}; + + for (size_t pos = 0; pos < number_of_user_bins; ++pos) + { + clusters.emplace_back(pos, positions[pos]); + current_cluster_cardinality[pos] = cardinalities[positions[pos]]; + current_cluster_sketches[pos] = sketches[positions[pos]]; + current_max_cluster_size = std::max(current_max_cluster_size, cardinalities[positions[pos]]); + } + + auto get_cluster = [&clusters](size_t const current) -> Cluster & + { + size_t const cluster_id = LSH_find_representative_cluster(clusters, current); + return clusters[cluster_id]; + }; + + // refine clusters + while (current_max_cluster_size < average_technical_bin_size + && current_sketch_index < number_of_max_minHash_sketches) + { + // fill LSH collision hashtable + robin_hood::unordered_flat_map> table = + LSH_fill_hashtable(clusters, minHash_sketches, current_sketch_index, minHash_sketch_size); + + // read out LSH collision hashtable + // for each present key, if the list contains more than one cluster, we merge everything contained in the list + // into the first cluster, since those clusters collide and should be joined in the same bucket + for (auto & [key, list] : table) + { + assert(!list.empty()); + + if (list.size() <= 1) // nothing to do here + continue; + + // Now combine all clusters into the first. + + // 1) find the representative cluster to merge everything else into + // It can happen, that the representative has already been joined with another cluster + // e.g. + // [key1] = {0,11} // then clusters[11] is merged into clusters[0] + // [key2] = {11,13} // now I want to merge clusters[13] into clusters[11] but the latter has been moved + auto & representative_cluster = get_cluster(list[0]); + assert(representative_cluster.id() == clusters[representative_cluster.id()].id()); + + auto & representative_cluster_sketch = current_cluster_sketches[representative_cluster.id()]; + + for (size_t const current : std::views::drop(list, 1)) + { + // For every other entry in the list, it can happen that I already joined that list with another + // e.g. + // [key1] = {0,11} // then clusters[11] is merged into clusters[0] + // [key2] = {0, 2, 11} // now I want to do it again + auto & next_cluster = get_cluster(current); + + if (next_cluster.id() == representative_cluster.id()) // already joined + continue; + + next_cluster.move_to(representative_cluster); // otherwise join next_cluster into representative_cluster + assert(next_cluster.empty()); + assert(next_cluster.has_been_moved()); + assert(representative_cluster.size() > 1); // there should be at least two user bins now + + representative_cluster_sketch.merge(current_cluster_sketches[next_cluster.id()]); + } + + current_cluster_cardinality[representative_cluster.id()] = representative_cluster_sketch.estimate(); + } + current_max_cluster_size = std::ranges::max(current_cluster_cardinality); + + ++current_sketch_index; + } + + return clusters; +} + +/*!\brief Orders the clusters so that lsh_sim_approach can take the leading ones as partition seeds. + * \param[in,out] clusters The clusters returned by very_similar_LSH_clustering. Reordered in place. + * \param[in] cardinalities The cardinality of each user bin, indexed by global user bin index. + * \param[in] config The configuration (uses `hibf_config.tmax`). + * \throws std::runtime_error if an empty cluster ends up before a non-empty one (sanity check). + * + * 1. The user bins inside each cluster are sorted by descending cardinality, so `contained_user_bins().front()` is + * the largest. + * 2. The first `tmax` positions receive the largest clusters by number of user bins, in descending order. Ties are + * broken by the cardinality of the largest user bin. + * 3. The remaining clusters are sorted by the cardinality of their largest user bin, in descending order. Empty + * (moved) clusters go last. + * + * After this, a cluster's position no longer matches its id(), so moved-to links (and + * LSH_find_representative_cluster) must not be used anymore. + */ +void post_process_clusters(std::vector & clusters, + std::vector const & cardinalities, + chopper::configuration const & config) +{ + // clusters are done. Start post processing + // since post processing involves re-ordering the clusters, the moved_to_cluster_id value of a cluster will not + // refer to the position of the cluster in the `clusters` vecto anymore but the cluster with the resprive id() + // would neet to be found + for (size_t pos = 0; pos < clusters.size(); ++pos) + { + assert(clusters[pos].is_valid(pos)); + clusters[pos].sort_by_cardinality(cardinalities); + } + + // push largest p clusters to the front + std::ranges::partial_sort(clusters, + std::ranges::next(clusters.begin(), config.hibf_config.tmax, clusters.end()), + [&cardinalities](auto const & v1, auto const & v2) + { + // Note: If v2 is empty, so is v1. + if (v1.size() == v2.size() && !v2.empty()) + return cardinalities[v1.contained_user_bins().front()] + > cardinalities[v2.contained_user_bins().front()]; + + return v1.size() > v2.size(); + }); + + // after filling up the partitions with the biggest clusters, sort the clusters by cardinality of the biggest ub + // s.t. that euqally sizes ub are assigned after each other and the small stuff is added at last. + // the largest ub is already at the start because of former sorting. + std::ranges::sort(std::ranges::next(clusters.begin(), config.hibf_config.tmax, clusters.end()), + clusters.end(), + [&cardinalities](auto const & v1, auto const & v2) + { + if (v1.empty()) + return false; // v1 can never be larger than v2 then + + if (v2.empty()) // and v1 is not, since the first if would catch + return true; + + return cardinalities[v1.contained_user_bins().front()] + > cardinalities[v2.contained_user_bins().front()]; + }); + + assert(clusters.size() < 2 || clusters[0].size() >= clusters[1].size()); // sanity check + +#ifndef NDEBUG + for (size_t cidx = 1; cidx < clusters.size(); ++cidx) + { + // once empty - always empty; all empty clusters should be at the end + if (clusters[cidx - 1].empty() && !clusters[cidx].empty()) + throw std::runtime_error{"sorting did not work"}; + } +#endif +} + +/*!\brief Assigns a whole cluster of user bins to the partition where adding it costs the least. + * \param[in] config The configuration (uses `hibf_config.tmax` and `sketch_bits`). + * \param[in] number_of_partitions Only partitions `[0, number_of_partitions)` are considered. + * \param[in,out] corrected_estimate_per_part The current target cardinality per partition. Raised to the chosen + * partition's new estimate if that is larger. + * \param[in] cluster The global user bin indices to assign together. + * \param[in] cardinalities The cardinality of each user bin, indexed by global user bin index. + * \param[in] sketches The HyperLogLog sketch of each user bin, indexed by global user bin index. + * \param[in,out] positions The user bins per partition. `cluster` is appended to the chosen one. + * \param[in,out] partition_sketches The union sketch per partition. Updated for the chosen partition. + * \param[in,out] max_partition_cardinality The largest user bin cardinality per partition. Updated. + * \param[in,out] min_partition_cardinality The smallest user bin cardinality per partition. Updated. + * \returns `true`. Throws std::runtime_error if no partition was selected, i.e., if `number_of_partitions == 0`. + * + * For each partition p, the cost of adding the cluster is the sum of + * - `union - |p|`: the new k-mers the cluster adds to p (HyperLogLog estimates), + * - `tmax * max(0, union - corrected_estimate_per_part)`: growth of the IBF bin size beyond the current target, and + * - a lower-level penalty depending on `max_card`, the largest user bin cardinality in `cluster`: + * - p already holds more than `tmax` user bins (lower level exists): `max_card * log_tmax(#UBs after adding)`, + * an estimate of how often the content is stored again on lower levels; + * - adding the cluster pushes p above `tmax` user bins (new lower level): `min(min_p, max_card) * tmax`; + * - otherwise: `max_p - max_card` if the cluster is smaller than every user bin in p (wasted space), or + * `(max_card - max_p) * tmax` if it is larger than every user bin in p (IBF grows), else 0. + * + * The partition with the smallest cost is chosen. A partition with zero cost always replaces the current best, so if + * several have zero cost, the last one wins. + */ +bool find_best_partition(chopper::configuration const & config, + size_t const number_of_partitions, + size_t & corrected_estimate_per_part, + std::vector const & cluster, + std::vector const & cardinalities, + std::vector const & sketches, + std::vector> & positions, + std::vector & partition_sketches, + std::vector & max_partition_cardinality, + std::vector & min_partition_cardinality) +{ + seqan::hibf::sketch::hyperloglog const current_sketch = [&sketches, &cluster, &config]() + { + seqan::hibf::sketch::hyperloglog result{config.hibf_config.sketch_bits}; + + for (size_t const user_bin_idx : cluster) + result.merge(sketches[user_bin_idx]); + + return result; + }(); + + size_t const max_card = [&cardinalities, &cluster]() + { + size_t max{0}; + + for (size_t const user_bin_idx : cluster) + max = std::max(max, cardinalities[user_bin_idx]); + + return max; + }(); + + // Search best partition fit by similarity. Similarity here is defined as: + // "whose (<-partition) effective text size is subsumed most by the current user bin". Or in other words: + // "which partition has the largest intersection with user bin b compared to its own (partition) size." + size_t smallest_change{std::numeric_limits::max()}; + size_t best_p{0}; + bool best_p_found{false}; + + auto penalty_lower_level = [&](size_t const additional_number_of_user_bins, size_t const p) -> size_t + { + assert(positions[p].size() != 0); // partitions should be initialised beforehand + size_t const min = min_partition_cardinality[p]; + size_t const max = max_partition_cardinality[p]; + + if (positions[p].size() > config.hibf_config.tmax) // already a third level + { + // if there must already be another lower level because the current merged bin contains more than tmax + // user bins, then the current user bin is very likely stored multiple times. Therefore, the penalty is set + // to the cardinality of the current user bin times the number of levels, e.g. the number of times this user + // bin needs to be stored additionally + size_t const num_ubs_in_merged_bin{positions[p].size() + additional_number_of_user_bins}; + double const levels = std::log(num_ubs_in_merged_bin) / std::log(config.hibf_config.tmax); + return static_cast(max_card * levels); + } + else if (positions[p].size() + additional_number_of_user_bins > config.hibf_config.tmax) // now a third level + { + // if the current merged bin contains exactly tmax UBS, adding otherone must + // result in another lower level. Most likely, the smallest user bin will end up on the lower level + // therefore the penalty is set to 'min * tmax' + // of course, there could also be a third level with a lower number of user bins, but this is hard to + // estimate. + size_t const penalty = std::min(min, max_card) * config.hibf_config.tmax; + return penalty; + } + else // positions[p].size() + additional_number_of_user_bins <= tmax + { + // if the new user bin is smaller than all other already contained user bins + // the waste of space is high if stored in a single technical bin + if (max_card < min) + return (max - max_card); + // if the new user bin is bigger than all other already contained user bins, the IBF size increases + else if (max_card > max) + return (max_card - max) * config.hibf_config.tmax; + // else, end if-else-block and zero is returned + } + + return 0u; + }; + + for (size_t p = 0; p < number_of_partitions; ++p) + { + seqan::hibf::sketch::hyperloglog union_sketch = current_sketch; + union_sketch.merge(partition_sketches[p]); + size_t const union_estimate = union_sketch.estimate(); + size_t const current_partition_size = partition_sketches[p].estimate(); + + assert(union_estimate >= current_partition_size); + size_t const penalty_current_bin = union_estimate - current_partition_size; + size_t const penalty_current_ibf = + config.hibf_config.tmax + * ((union_estimate <= corrected_estimate_per_part) ? 0u : union_estimate - corrected_estimate_per_part); + size_t const change = penalty_current_bin + penalty_current_ibf + penalty_lower_level(cluster.size(), p); + + if (change == 0 || /* If there is no penalty at all, this is a best fit even if the partition is "full"*/ + (smallest_change > change)) + { + smallest_change = change; + best_p = p; + best_p_found = true; + } + } + + if (!best_p_found) + throw std::runtime_error{"currently there are no safety measures if a partition is not found"}; + + // now that we know which partition fits best (`best_p`), add those indices to it + for (size_t const user_bin_idx : cluster) + { + positions[best_p].push_back(user_bin_idx); + max_partition_cardinality[best_p] = std::max(max_partition_cardinality[best_p], cardinalities[user_bin_idx]); + min_partition_cardinality[best_p] = std::min(min_partition_cardinality[best_p], cardinalities[user_bin_idx]); + } + partition_sketches[best_p].merge(current_sketch); + corrected_estimate_per_part = std::max(corrected_estimate_per_part, partition_sketches[best_p].estimate()); + + return true; +} + +/*!\brief Distributes the merged-bin candidates onto `number_of_remaining_tbs` partitions by LSH clustering and + * similarity-based assignment. + * \param[in] config The configuration (uses `hibf_config` and the LSH/search timers). + * \param[in] sorted_positions2 The global indices of the user bins to distribute, sorted by descending + * cardinality (the remainder after split bins were removed). + * \param[in] cardinalities The cardinality of each user bin, indexed by global user bin index. + * \param[in] sketches The HyperLogLog sketch of each user bin, indexed by global user bin index. + * \param[in] minHash_sketches The MinHash tables of each user bin, indexed by global user bin index. + * \param[in,out] partitions Receives the assignment in `partitions[0, number_of_remaining_tbs)`. + * Must have at least `number_of_remaining_tbs` entries. + * \param[in] number_of_remaining_tbs The number of technical bins available for merged bins. + * \param[in] technical_bin_size_threshold The target cardinality per technical bin. Stops LSH clustering and + * triggers spill-over while seeding partitions. + * \param[in] sum_of_cardinalities The sum of the cardinalities of *all* user bins of this IBF. + * \returns The largest estimated cardinality of any of the `number_of_remaining_tbs` partitions. + * + * 1. **Cluster:** very_similar_LSH_clustering followed by post_process_clusters. + * 2. **Ensure enough clusters:** If there are fewer non-empty clusters than partitions, user bins are moved out of + * the last cluster with more than one user bin, into their own clusters, until there are enough clusters. + * 3. **Seed:** Partition p receives the next cluster in order. A cluster is spread over consecutive partitions + * whenever `tmax` user bins have been placed or the partition's estimate exceeds `technical_bin_size_threshold`. + * A cluster with more than `tmax` user bins whose total cardinality exceeds `0.05 * sum_of_cardinalities / tmax` + * is placed in full. Any other cluster places at most `tmax` user bins. User bins not placed, or left over + * because the partitions ran out, go to a list of remaining clusters. + * 4. **Assign the rest:** Each remaining cluster, plus all clusters that were not used as seeds, is assigned as a + * whole by find_best_partition. The per-partition target starts at `technical_bin_size_threshold` and only grows. + */ +size_t lsh_sim_approach(chopper::configuration const & config, + std::vector const & sorted_positions2, + std::vector const & cardinalities, + std::vector const & sketches, + std::vector const & minHash_sketches, + std::vector> & partitions, + size_t const number_of_remaining_tbs, + size_t const technical_bin_size_threshold, + size_t const sum_of_cardinalities) +{ + uint8_t const sketch_bits{config.hibf_config.sketch_bits}; + std::vector partition_sketches(number_of_remaining_tbs, + seqan::hibf::sketch::hyperloglog(sketch_bits)); + + std::vector max_partition_cardinality(number_of_remaining_tbs, 0u); + std::vector min_partition_cardinality(number_of_remaining_tbs, std::numeric_limits::max()); + + // initial partitioning using locality sensitive hashing (LSH) + config.lsh_algorithm_timer.start(); + std::vector clusters = very_similar_LSH_clustering(minHash_sketches, + sorted_positions2, + cardinalities, + sketches, + technical_bin_size_threshold, + config); + post_process_clusters(clusters, cardinalities, config); + config.lsh_algorithm_timer.stop(); + + // There must be more non-empty clusters than technical bins + size_t number_of_remaining_clusters = std::ranges::count_if(clusters, + [](auto const & c) + { + return !c.empty(); + }); + clusters.resize(std::max(number_of_remaining_tbs, clusters.size())); + while (number_of_remaining_tbs > number_of_remaining_clusters) + { + // If there are not enough remaining clusters, we need to split up clusters. This is rather an edge case. + auto empty_cluster_it = clusters.begin() + number_of_remaining_clusters; + // We can only split up clusters of size > 1. Start with smaller ones to the back of `clusters` + auto cluster_it = std::ranges::find_if(clusters.rbegin(), + clusters.rend(), + [](auto const & c) + { + return c.size() > 1; + }); + assert(cluster_it != clusters.rend()); + for (; cluster_it->size() != 1 && number_of_remaining_tbs > number_of_remaining_clusters;) + { + assert(empty_cluster_it->empty()); + empty_cluster_it->add_user_bin(cluster_it->pop_back()); // not super efficient but fine + ++number_of_remaining_clusters; + ++empty_cluster_it; + } + } + assert(number_of_remaining_tbs <= static_cast(std::ranges::count_if(clusters, + [](auto const & c) + { + return !c.empty(); + }))); + + std::vector> remaining_clusters{}; + + // initialise partitions with the first p largest clusters (post_processing sorts by size) + size_t cidx{0}; // current cluster index + for (size_t p = 0; p < number_of_remaining_tbs; ++p) + { + assert(!clusters[cidx].empty()); + auto const & cluster = clusters[cidx].contained_user_bins(); + bool split_cluster = false; + + if (cluster.size() > config.hibf_config.tmax) + { + size_t card{0}; + for (size_t uidx = 0; uidx < cluster.size(); ++uidx) + card += cardinalities[cluster[uidx]]; + + if (card > 0.05 * sum_of_cardinalities / config.hibf_config.tmax) + split_cluster = true; + } + + size_t end = (split_cluster) ? cluster.size() : std::min(cluster.size(), config.hibf_config.tmax); + for (size_t uidx = 0; uidx < end; ++uidx) + { + size_t const user_bin_idx = cluster[uidx]; + // if a single cluster already exceeds the cardinality_per_part, + // then the remaining user bins of the cluster must spill over into the next partition + if ((uidx != 0 && (uidx % config.hibf_config.tmax == 0)) + || partition_sketches[p].estimate() > technical_bin_size_threshold) + { + ++p; + + if (p >= number_of_remaining_tbs) + { + split_cluster = true; + end = uidx; + break; + } + } + + partition_sketches[p].merge(sketches[user_bin_idx]); + partitions[p].push_back(user_bin_idx); + max_partition_cardinality[p] = std::max(max_partition_cardinality[p], cardinalities[user_bin_idx]); + min_partition_cardinality[p] = std::min(min_partition_cardinality[p], cardinalities[user_bin_idx]); + } + + if (split_cluster) + { + std::vector remainder(cluster.begin() + end, cluster.end()); + remaining_clusters.insert(remaining_clusters.end(), remainder); + } + + ++cidx; + } + + for (size_t i = cidx; i < clusters.size(); ++i) + { + if (clusters[i].empty()) + break; + + remaining_clusters.insert(remaining_clusters.end(), clusters[i].contained_user_bins()); + } + + // assign the rest by similarity + size_t merged_threshold{technical_bin_size_threshold}; + for (size_t ridx = 0; ridx < remaining_clusters.size(); ++ridx) + { + auto const & cluster = remaining_clusters[ridx]; + + config.search_partition_algorithm_timer.start(); + find_best_partition(config, + number_of_remaining_tbs, + merged_threshold, + cluster, + cardinalities, + sketches, + partitions, + partition_sketches, + max_partition_cardinality, + min_partition_cardinality); + config.search_partition_algorithm_timer.stop(); + } + + // compute actual max size + size_t max_size{0}; + for (auto const & sketch : partition_sketches) + max_size = std::max(max_size, (size_t)sketch.estimate()); + + return max_size; +} + +void partition_user_bins(chopper::configuration const & config, + std::vector const & positions, + std::vector const & cardinalities, + std::vector const & sketches, + std::vector const & minHash_sketches, + std::vector> & partitions) +{ + // all approaches need sorted positions + std::vector const sorted_positions = [&positions, &cardinalities]() + { + std::vector ps(positions.begin(), positions.end()); + seqan::hibf::sketch::toolbox::sort_by_cardinalities(cardinalities, ps); + return ps; + }(); + + auto const [sum_of_cardinalities, joint_estimate] = [&]() + { + size_t sum{0}; + seqan::hibf::sketch::hyperloglog sketch{config.hibf_config.sketch_bits}; + + for (size_t const pos : positions) + { + sum += cardinalities[pos]; + sketch.merge(sketches[pos]); + } + + return std::tuple{sum, sketch.estimate()}; + }(); + + double const relaxed_fpr_correction = seqan::hibf::layout::compute_relaxed_fpr_correction( + {.fpr = config.hibf_config.maximum_fpr, // + .relaxed_fpr = config.hibf_config.relaxed_fpr, + .hash_count = config.hibf_config.number_of_hash_functions}); + + size_t idx{0}; // start in sorted positions + size_t number_of_split_tbs{0}; + size_t number_of_merged_tbs{config.hibf_config.tmax}; + size_t max_split_size{0}; + size_t max_merged_size{0}; + + // Start with the lower bound on the split threshold: All user bin content is the same and when split/merged + // one technical bin is never more than joint_estimate/tmax per bin. + // (unrealistic but a good starting point. threshold will be revised in a second iteration) + size_t split_threshold = seqan::hibf::divide_and_ceil(joint_estimate, config.hibf_config.tmax); + + auto partition_split_bins = [&]() + { + size_t number_of_potential_split_bins{0}; // determined in find_bins_to_be_split + size_t max_tbs{config.hibf_config.tmax - 1}; // leave one bin for merging + + std::tie(idx, number_of_potential_split_bins) = + find_bins_to_be_split(sorted_positions, cardinalities, split_threshold, max_tbs); + + // if (remaining UBs < remaining TBs for merging) + // then assign more split bins, such that there are at least as manu UBs as TBs for merging + if (sorted_positions.size() - idx < config.hibf_config.tmax - number_of_potential_split_bins) + number_of_potential_split_bins += + (config.hibf_config.tmax - number_of_potential_split_bins) - (sorted_positions.size() - idx); + + std::tie(number_of_split_tbs, max_split_size) = + chopper::layout::determine_split_bins(config, + sorted_positions, + cardinalities, + number_of_potential_split_bins, + idx, + partitions); + number_of_merged_tbs = config.hibf_config.tmax - number_of_split_tbs; + }; + + auto partition_merged_bins = [&]() + { + // determine number of split bins + std::vector const sorted_positions2(sorted_positions.begin() + idx, sorted_positions.end()); + + // distribute the rest to merged bins + size_t const corrected_max_split_size = max_split_size / relaxed_fpr_correction; + size_t const merged_threshold = std::max(corrected_max_split_size, split_threshold); + + max_merged_size = lsh_sim_approach(config, + sorted_positions2, + cardinalities, + sketches, + minHash_sketches, + partitions, + number_of_merged_tbs, + merged_threshold, + sum_of_cardinalities); + }; + + partition_split_bins(); + + // All user bins can be assigned as split bins (idx == sorted_positions.size()). + // In that case there are no remaining bins to distribute via partition_merged_bins. + // And no reconfiguration of the threshold needs to be done since splitting is done with an "optimal" DP + if (idx < sorted_positions.size()) + { + partition_merged_bins(); + + int64_t const difference = + static_cast(max_merged_size * relaxed_fpr_correction) - static_cast(max_split_size); + + if (number_of_split_tbs == 0) + split_threshold = (split_threshold + max_merged_size) / 2; // increase threshold + else if (difference > 0) // need more merged bins -> increase threshold + split_threshold = static_cast(split_threshold) + * ((static_cast(max_merged_size) * relaxed_fpr_correction) + / static_cast(max_split_size)); + else // need more split bins -> decrease threshold + split_threshold = std::max(1.0, + static_cast(split_threshold) + * ((static_cast(max_merged_size) * relaxed_fpr_correction) + / static_cast(max_split_size))); + + // reset result + partitions.clear(); + partitions.resize(config.hibf_config.tmax); + idx = 0; + number_of_split_tbs = 0; + number_of_merged_tbs = config.hibf_config.tmax; + max_split_size = 0; + max_merged_size = 0; + + partition_split_bins(); + // All user bins can be assigned as split bins (idx == sorted_positions.size()). + // In that case there are no remaining bins to distribute via partition_merged_bins. + if (idx < sorted_positions.size()) + partition_merged_bins(); + } +} + +} // namespace chopper::layout diff --git a/src/set_up_parser.cpp b/src/set_up_parser.cpp index 08c36b4f..e34f6a2a 100644 --- a/src/set_up_parser.cpp +++ b/src/set_up_parser.cpp @@ -245,6 +245,12 @@ void set_up_parser(sharg::parser & parser, configuration & config) "ignored and has no effect.", .advanced = true}); + parser.add_flag(config.fast_layout, + sharg::config{.short_id = '\0', + .long_id = "fast-layout", + .description = "Uses the fast layout algorithm instead of the default one.", + .advanced = false}); + parser.add_flag( config.output_verbose_statistics, sharg::config{.short_id = '\0', diff --git a/test/api/layout/CMakeLists.txt b/test/api/layout/CMakeLists.txt index 9cdb4076..92101d16 100644 --- a/test/api/layout/CMakeLists.txt +++ b/test/api/layout/CMakeLists.txt @@ -26,3 +26,7 @@ if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU" endif () add_api_test (user_bin_io_test.cpp) add_api_test (input_test.cpp) +add_api_test (determine_split_bins_test.cpp) +add_api_test (fast_layout_cluster_test.cpp) +add_api_test (fast_layout_find_bins_to_be_split_test.cpp) +add_api_test (partition_user_bins_test.cpp) diff --git a/test/api/layout/determine_split_bins_test.cpp b/test/api/layout/determine_split_bins_test.cpp new file mode 100644 index 00000000..f061fd2d --- /dev/null +++ b/test/api/layout/determine_split_bins_test.cpp @@ -0,0 +1,58 @@ +#include // for Test, TestInfo, EXPECT_EQ, Message, TEST, TestPartResult + +#include // for allocator, string +#include // for vector + +#include +#include + +TEST(simple_split, first) +{ + chopper::configuration config{}; + + std::vector cardinalities{}; + for (size_t i{0}; i < 20; ++i) + cardinalities.push_back(2400); + cardinalities.push_back(200); + for (size_t i{0}; i < 84; ++i) + cardinalities.push_back(50); + + std::vector positions(cardinalities.size()); + std::iota(positions.begin(), positions.end(), 0u); + + size_t num_technical_bins{63}; + size_t num_user_bins{21}; + + std::vector> partitions(64); + + auto const [num_splits, max_size] = chopper::layout::determine_split_bins(config, + positions, + cardinalities, + num_technical_bins, + num_user_bins, + partitions); + + EXPECT_EQ(num_splits, 61); + EXPECT_EQ(max_size, 1452); + + size_t ub_idx{20}; + size_t tb_idx{partitions.size() - 1}; + + // TBs: 0, 1, 2, 3, ... + // UBs: 20,19,19,19,18,18,18,17,17,17,16,16,16,15,15,15,14,14,14,13,13,13,12,12,12,11,11,11,10,10,10,9,9,9,8,8,8,7,7,7,6,6,6,5,5,5,4,4,4,3,3,3,2,2,2,1,1,1,0,0,0, + + EXPECT_EQ(partitions[tb_idx][0], ub_idx); // single bin with card = 200 + --tb_idx; + --ub_idx; + + for (; tb_idx > partitions.size() - num_splits - 1; --tb_idx) + { + for (size_t three = 0; three < 3; ++three) + { + EXPECT_EQ(partitions[tb_idx][0], ub_idx); + --tb_idx; + } + ++tb_idx; // one too much in loop before + --ub_idx; + } +} \ No newline at end of file diff --git a/test/api/layout/execute_layout_test.cpp b/test/api/layout/execute_layout_test.cpp index 412dd27b..41948467 100644 --- a/test/api/layout/execute_layout_test.cpp +++ b/test/api/layout/execute_layout_test.cpp @@ -46,9 +46,10 @@ TEST(execute_test, few_ubs) filenames{{"seq0a", "seq0b"}, {"seq1"}, {"seq2"}, {"seq3"}, {"seq4"}, {"seq5"}, {"seq6"}, {"seq7"}}; std::vector sketches; - seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches); + std::vector minHash_sketches{}; + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); - chopper::layout::execute(config, filenames, sketches); + chopper::layout::execute(config, filenames, sketches, minHash_sketches); std::string const expected_file{"@CHOPPER_USER_BINS\n" "@0 seq0a seq0b\n" @@ -120,6 +121,106 @@ TEST(execute_test, few_ubs) EXPECT_EQ(actual_file, expected_file) << actual_file; } +TEST(execute_test, few_ubs_fast_layout) +{ + seqan3::test::tmp_directory tmp_dir{}; + std::filesystem::path const layout_file{tmp_dir.path() / "layout.tsv"}; + + auto simulated_input = [&](size_t const num, seqan::hibf::insert_iterator it) + { + size_t const desired_kmer_count = (num == 1) ? 1760 : 940; + for (auto hash : std::views::iota(0u, desired_kmer_count)) + it = hash; + }; + + chopper::configuration config{}; + config.fast_layout = true; + config.hibf_config.input_fn = simulated_input; + config.hibf_config.number_of_user_bins = 8; + config.hibf_config.tmax = 64; + config.output_filename = layout_file; + config.disable_sketch_output = true; + config.hibf_config.disable_estimate_union = true; // also disables rearrangement + + std::vector> + filenames{{"seq0a", "seq0b"}, {"seq1"}, {"seq2"}, {"seq3"}, {"seq4"}, {"seq5"}, {"seq6"}, {"seq7"}}; + + std::vector sketches; + std::vector minHash_sketches{}; + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); + + chopper::layout::execute(config, filenames, sketches, minHash_sketches); + + std::string const expected_file{"@CHOPPER_USER_BINS\n" + "@0 seq0a seq0b\n" + "@1 seq1\n" + "@2 seq2\n" + "@3 seq3\n" + "@4 seq4\n" + "@5 seq5\n" + "@6 seq6\n" + "@7 seq7\n" + "@CHOPPER_USER_BINS_END\n" + "@CHOPPER_CONFIG\n" + "@{\n" + "@ \"chopper_config\": {\n" + "@ \"version\": 2,\n" + "@ \"data_file\": {\n" + "@ \"value0\": \"\"\n" + "@ },\n" + "@ \"debug\": false,\n" + "@ \"sketch_directory\": {\n" + "@ \"value0\": \"\"\n" + "@ },\n" + "@ \"k\": 19,\n" + "@ \"window_size\": 19,\n" + "@ \"disable_sketch_output\": true,\n" + "@ \"precomputed_files\": false,\n" + "@ \"output_filename\": {\n" + "@ \"value0\": \"" + + layout_file.string() + + "\"\n" + "@ },\n" + "@ \"determine_best_tmax\": false,\n" + "@ \"force_all_binnings\": false\n" + "@ }\n" + "@}\n" + "@CHOPPER_CONFIG_END\n" + "@HIBF_CONFIG\n" + "@{\n" + "@ \"hibf_config\": {\n" + "@ \"version\": 3,\n" + "@ \"number_of_user_bins\": 8,\n" + "@ \"number_of_hash_functions\": 2,\n" + "@ \"maximum_fpr\": 0.05,\n" + "@ \"relaxed_fpr\": 0.3,\n" + "@ \"threads\": 1,\n" + "@ \"sketch_bits\": 12,\n" + "@ \"tmax\": 64,\n" + "@ \"empty_bin_fraction\": 0.0,\n" + "@ \"track_occupancy\": false,\n" + "@ \"alpha\": 1.2,\n" + "@ \"max_rearrangement_ratio\": 0.5,\n" + "@ \"disable_estimate_union\": true,\n" + "@ \"disable_rearrangement\": true\n" + "@ }\n" + "@}\n" + "@HIBF_CONFIG_END\n" + "#TOP_LEVEL_IBF fullest_technical_bin_idx:4\n" + "#USER_BIN_IDX\tTECHNICAL_BIN_INDICES\tNUMBER_OF_TECHNICAL_BINS\n" + "0\t29\t5\n" + "1\t4\t25\n" + "2\t34\t5\n" + "3\t39\t5\n" + "4\t44\t5\n" + "5\t49\t5\n" + "6\t54\t5\n" + "7\t59\t5\n"}; + std::string const actual_file{string_from_file(layout_file)}; + + EXPECT_EQ(actual_file, expected_file) << actual_file; +} + TEST(execute_test, set_default_tmax) { seqan3::test::tmp_directory tmp_dir{}; @@ -143,9 +244,10 @@ TEST(execute_test, set_default_tmax) filenames{{"seq0"}, {"seq1"}, {"seq2"}, {"seq3"}, {"seq4"}, {"seq5"}, {"seq6"}, {"seq7"}}; std::vector sketches; - seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches); + std::vector minHash_sketches{}; + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); - chopper::layout::execute(config, filenames, sketches); + chopper::layout::execute(config, filenames, sketches, minHash_sketches); EXPECT_EQ(config.hibf_config.tmax, 64u); } @@ -179,9 +281,10 @@ TEST(execute_test, many_ubs) config.hibf_config.disable_estimate_union = true; // also disables rearrangement std::vector sketches; - seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches); + std::vector minHash_sketches{}; + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); - chopper::layout::execute(config, many_filenames, sketches); + chopper::layout::execute(config, many_filenames, sketches, minHash_sketches); std::string const expected_file{"@CHOPPER_USER_BINS\n" "@0 seq0\n" @@ -438,3 +541,286 @@ TEST(execute_test, many_ubs) EXPECT_EQ(actual_file, expected_file) << actual_file << std::endl; } + +TEST(execute_test, many_ubs_fast_layout) +{ + seqan3::test::tmp_directory tmp_dir{}; + std::filesystem::path const layout_file{tmp_dir.path() / "layout.tsv"}; + + std::vector> many_filenames; + + for (size_t i{0}; i < 96u; ++i) + many_filenames.push_back({seqan3::detail::to_string("seq", i)}); + + // Creates sizes of the following series + // [801,802,...,820,922,923,...,941,1043,1044,...,1062,1164,1165,...,1183,1285,1286,...,1300] + // See also https://godbolt.org/z/9517eaaaG + auto simulated_input = [&](size_t const num, seqan::hibf::insert_iterator it) + { + size_t const desired_kmer_count = 101 * ((num + 20) / 20) + num + 700; + for (auto hash : std::views::iota(0u, desired_kmer_count)) + it = hash; + }; + + chopper::configuration config{}; + config.fast_layout = true; + config.output_filename = layout_file; + config.disable_sketch_output = true; + config.hibf_config.tmax = 64; + config.hibf_config.input_fn = simulated_input; + config.hibf_config.number_of_user_bins = many_filenames.size(); + config.hibf_config.disable_estimate_union = true; // also disables rearrangement + + std::vector sketches; + std::vector minHash_sketches{}; + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); + + chopper::layout::execute(config, many_filenames, sketches, minHash_sketches); + + std::string const expected_file{"@CHOPPER_USER_BINS\n" + "@0 seq0\n" + "@1 seq1\n" + "@2 seq2\n" + "@3 seq3\n" + "@4 seq4\n" + "@5 seq5\n" + "@6 seq6\n" + "@7 seq7\n" + "@8 seq8\n" + "@9 seq9\n" + "@10 seq10\n" + "@11 seq11\n" + "@12 seq12\n" + "@13 seq13\n" + "@14 seq14\n" + "@15 seq15\n" + "@16 seq16\n" + "@17 seq17\n" + "@18 seq18\n" + "@19 seq19\n" + "@20 seq20\n" + "@21 seq21\n" + "@22 seq22\n" + "@23 seq23\n" + "@24 seq24\n" + "@25 seq25\n" + "@26 seq26\n" + "@27 seq27\n" + "@28 seq28\n" + "@29 seq29\n" + "@30 seq30\n" + "@31 seq31\n" + "@32 seq32\n" + "@33 seq33\n" + "@34 seq34\n" + "@35 seq35\n" + "@36 seq36\n" + "@37 seq37\n" + "@38 seq38\n" + "@39 seq39\n" + "@40 seq40\n" + "@41 seq41\n" + "@42 seq42\n" + "@43 seq43\n" + "@44 seq44\n" + "@45 seq45\n" + "@46 seq46\n" + "@47 seq47\n" + "@48 seq48\n" + "@49 seq49\n" + "@50 seq50\n" + "@51 seq51\n" + "@52 seq52\n" + "@53 seq53\n" + "@54 seq54\n" + "@55 seq55\n" + "@56 seq56\n" + "@57 seq57\n" + "@58 seq58\n" + "@59 seq59\n" + "@60 seq60\n" + "@61 seq61\n" + "@62 seq62\n" + "@63 seq63\n" + "@64 seq64\n" + "@65 seq65\n" + "@66 seq66\n" + "@67 seq67\n" + "@68 seq68\n" + "@69 seq69\n" + "@70 seq70\n" + "@71 seq71\n" + "@72 seq72\n" + "@73 seq73\n" + "@74 seq74\n" + "@75 seq75\n" + "@76 seq76\n" + "@77 seq77\n" + "@78 seq78\n" + "@79 seq79\n" + "@80 seq80\n" + "@81 seq81\n" + "@82 seq82\n" + "@83 seq83\n" + "@84 seq84\n" + "@85 seq85\n" + "@86 seq86\n" + "@87 seq87\n" + "@88 seq88\n" + "@89 seq89\n" + "@90 seq90\n" + "@91 seq91\n" + "@92 seq92\n" + "@93 seq93\n" + "@94 seq94\n" + "@95 seq95\n" + "@CHOPPER_USER_BINS_END\n" + "@CHOPPER_CONFIG\n" + "@{\n" + "@ \"chopper_config\": {\n" + "@ \"version\": 2,\n" + "@ \"data_file\": {\n" + "@ \"value0\": \"\"\n" + "@ },\n" + "@ \"debug\": false,\n" + "@ \"sketch_directory\": {\n" + "@ \"value0\": \"\"\n" + "@ },\n" + "@ \"k\": 19,\n" + "@ \"window_size\": 19,\n" + "@ \"disable_sketch_output\": true,\n" + "@ \"precomputed_files\": false,\n" + "@ \"output_filename\": {\n" + "@ \"value0\": \"" + + layout_file.string() + + "\"\n" + "@ },\n" + "@ \"determine_best_tmax\": false,\n" + "@ \"force_all_binnings\": false\n" + "@ }\n" + "@}\n" + "@CHOPPER_CONFIG_END\n" + "@HIBF_CONFIG\n" + "@{\n" + "@ \"hibf_config\": {\n" + "@ \"version\": 3,\n" + "@ \"number_of_user_bins\": 96,\n" + "@ \"number_of_hash_functions\": 2,\n" + "@ \"maximum_fpr\": 0.05,\n" + "@ \"relaxed_fpr\": 0.3,\n" + "@ \"threads\": 1,\n" + "@ \"sketch_bits\": 12,\n" + "@ \"tmax\": 64,\n" + "@ \"empty_bin_fraction\": 0.0,\n" + "@ \"track_occupancy\": false,\n" + "@ \"alpha\": 1.2,\n" + "@ \"max_rearrangement_ratio\": 0.5,\n" + "@ \"disable_estimate_union\": true,\n" + "@ \"disable_rearrangement\": true\n" + "@ }\n" + "@}\n" + "@HIBF_CONFIG_END\n" + "#TOP_LEVEL_IBF fullest_technical_bin_idx:0\n" + "#LOWER_LEVEL_IBF_63 fullest_technical_bin_idx:62\n" + "#LOWER_LEVEL_IBF_63;53 fullest_technical_bin_idx:0\n" + "#USER_BIN_IDX\tTECHNICAL_BIN_INDICES\tNUMBER_OF_TECHNICAL_BINS\n" + "0\t63;0\t1;3\n" + "1\t63;3\t1;2\n" + "2\t63;5\t1;2\n" + "3\t63;7\t1;2\n" + "4\t63;9\t1;2\n" + "5\t63;11\t1;2\n" + "6\t63;13\t1;2\n" + "7\t63;15\t1;2\n" + "8\t63;17\t1;2\n" + "9\t63;19\t1;2\n" + "10\t63;21\t1;2\n" + "11\t63;23\t1;2\n" + "12\t63;25\t1;2\n" + "13\t63;27\t1;2\n" + "14\t63;29\t1;2\n" + "15\t63;31\t1;2\n" + "16\t63;33\t1;2\n" + "17\t63;35\t1;2\n" + "18\t63;37\t1;2\n" + "19\t63;39\t1;2\n" + "20\t63;41\t1;2\n" + "21\t63;43\t1;2\n" + "22\t63;45\t1;2\n" + "23\t63;47\t1;2\n" + "24\t63;49\t1;2\n" + "25\t63;51\t1;2\n" + "26\t63;53;32\t1;1;32\n" + "27\t63;53;0\t1;1;32\n" + "28\t63;54\t1;2\n" + "29\t63;56\t1;2\n" + "30\t63;58\t1;2\n" + "31\t63;60\t1;2\n" + "32\t63;62\t1;2\n" + "33\t62\t1\n" + "34\t61\t1\n" + "35\t60\t1\n" + "36\t59\t1\n" + "37\t58\t1\n" + "38\t57\t1\n" + "39\t56\t1\n" + "40\t55\t1\n" + "41\t54\t1\n" + "42\t53\t1\n" + "43\t52\t1\n" + "44\t51\t1\n" + "45\t50\t1\n" + "46\t49\t1\n" + "47\t48\t1\n" + "48\t47\t1\n" + "49\t46\t1\n" + "50\t45\t1\n" + "51\t44\t1\n" + "52\t43\t1\n" + "53\t42\t1\n" + "54\t41\t1\n" + "55\t40\t1\n" + "56\t39\t1\n" + "57\t38\t1\n" + "58\t37\t1\n" + "59\t36\t1\n" + "60\t35\t1\n" + "61\t34\t1\n" + "62\t33\t1\n" + "63\t32\t1\n" + "64\t31\t1\n" + "65\t30\t1\n" + "66\t29\t1\n" + "67\t28\t1\n" + "68\t27\t1\n" + "69\t26\t1\n" + "70\t25\t1\n" + "71\t24\t1\n" + "72\t23\t1\n" + "73\t22\t1\n" + "74\t21\t1\n" + "75\t20\t1\n" + "76\t19\t1\n" + "77\t18\t1\n" + "78\t17\t1\n" + "79\t16\t1\n" + "80\t15\t1\n" + "81\t14\t1\n" + "82\t13\t1\n" + "83\t12\t1\n" + "84\t11\t1\n" + "85\t10\t1\n" + "86\t9\t1\n" + "87\t8\t1\n" + "88\t7\t1\n" + "89\t6\t1\n" + "90\t5\t1\n" + "91\t4\t1\n" + "92\t3\t1\n" + "93\t2\t1\n" + "94\t1\t1\n" + "95\t0\t1\n"}; + std::string const actual_file{string_from_file(layout_file)}; + + EXPECT_EQ(actual_file, expected_file) << actual_file << std::endl; +} diff --git a/test/api/layout/execute_with_estimation_test.cpp b/test/api/layout/execute_with_estimation_test.cpp index 7ac2db6a..518a2a60 100644 --- a/test/api/layout/execute_with_estimation_test.cpp +++ b/test/api/layout/execute_with_estimation_test.cpp @@ -53,9 +53,68 @@ TEST(execute_estimation_test, few_ubs) filenames{{"seq0"}, {"seq1"}, {"seq2"}, {"seq3"}, {"seq4"}, {"seq5"}, {"seq6"}, {"seq7"}}; std::vector sketches; - seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches); + std::vector minHash_sketches{}; + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); - chopper::layout::execute(config, filenames, sketches); + chopper::layout::execute(config, filenames, sketches, minHash_sketches); + + ASSERT_TRUE(std::filesystem::exists(stats_file)); + + std::string const written_file{string_from_file(stats_file)}; + + EXPECT_EQ(written_file, + R"expected_cout(## ### Parameters ### +## number of user bins = 8 +## number of hash functions = 2 +## maximum false positive rate = 0.05 +## relaxed false positive rate = 0.3 +## ### Notation ### +## X-IBF = An IBF with X number of bins. +## X-HIBF = An HIBF with tmax = X, e.g a maximum of X technical bins on each level. +## ### Column Description ### +## tmax : The maximum number of technical bin on each level +## c_tmax : The technical extra cost of querying an tmax-IBF, compared to 64-IBF +## l_tmax : The estimated query cost for an tmax-HIBF, compared to an 64-HIBF +## m_tmax : The estimated memory consumption for an tmax-HIBF, compared to an 64-HIBF +## (l*m)_tmax : Computed by l_tmax * m_tmax +## size : The expected total size of an tmax-HIBF +# tmax c_tmax l_tmax m_tmax (l*m)_tmax size +64 1.00 1.00 1.00 1.00 38.2KiB +# Best t_max (regarding expected query runtime): 64 +)expected_cout"); +} + +TEST(execute_estimation_test, few_ubs_fast_layout) +{ + seqan3::test::tmp_directory tmp_dir{}; + std::filesystem::path const layout_file{tmp_dir.path() / "layout.tsv"}; + std::filesystem::path const stats_file{layout_file.string() + ".stats"}; + + auto simulated_input = [&](size_t const num, seqan::hibf::insert_iterator it) + { + size_t const desired_kmer_count = (num == 1) ? 1700 : 1200; + for (auto hash : std::views::iota(0u, desired_kmer_count)) + it = hash; + }; + + chopper::configuration config{}; + config.fast_layout = true; + config.hibf_config.tmax = 64; + config.hibf_config.input_fn = simulated_input; + config.hibf_config.number_of_user_bins = 8; + config.determine_best_tmax = true; + config.disable_sketch_output = true; + config.output_filename = layout_file; + config.hibf_config.disable_estimate_union = true; // also disables rearrangement + + std::vector> + filenames{{"seq0"}, {"seq1"}, {"seq2"}, {"seq3"}, {"seq4"}, {"seq5"}, {"seq6"}, {"seq7"}}; + + std::vector sketches; + std::vector minHash_sketches{}; + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); + + chopper::layout::execute(config, filenames, sketches, minHash_sketches); ASSERT_TRUE(std::filesystem::exists(stats_file)); @@ -114,9 +173,328 @@ TEST(execute_estimation_test, many_ubs) config.hibf_config.disable_estimate_union = true; // also disables rearrangement std::vector sketches; - seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches); + std::vector minHash_sketches{}; + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); + + chopper::layout::execute(config, many_filenames, sketches, minHash_sketches); + + ASSERT_TRUE(std::filesystem::exists(stats_file)); + + std::string const written_file{string_from_file(stats_file)}; + + EXPECT_EQ(written_file, + R"expected_cout(## ### Parameters ### +## number of user bins = 96 +## number of hash functions = 2 +## maximum false positive rate = 0.05 +## relaxed false positive rate = 0.3 +## ### Notation ### +## X-IBF = An IBF with X number of bins. +## X-HIBF = An HIBF with tmax = X, e.g a maximum of X technical bins on each level. +## ### Column Description ### +## tmax : The maximum number of technical bin on each level +## c_tmax : The technical extra cost of querying an tmax-IBF, compared to 64-IBF +## l_tmax : The estimated query cost for an tmax-HIBF, compared to an 64-HIBF +## m_tmax : The estimated memory consumption for an tmax-HIBF, compared to an 64-HIBF +## (l*m)_tmax : Computed by l_tmax * m_tmax +## size : The expected total size of an tmax-HIBF +# tmax c_tmax l_tmax m_tmax (l*m)_tmax size +64 1.00 1.36 1.00 1.36 269.0KiB +128 1.22 1.42 1.00 1.42 269.8KiB +# Best t_max (regarding expected query runtime): 64 +)expected_cout"); + + std::string const expected_file{"@CHOPPER_USER_BINS\n" + "@0 seq0\n" + "@1 seq1\n" + "@2 seq2\n" + "@3 seq3\n" + "@4 seq4\n" + "@5 seq5\n" + "@6 seq6\n" + "@7 seq7\n" + "@8 seq8\n" + "@9 seq9\n" + "@10 seq10\n" + "@11 seq11\n" + "@12 seq12\n" + "@13 seq13\n" + "@14 seq14\n" + "@15 seq15\n" + "@16 seq16\n" + "@17 seq17\n" + "@18 seq18\n" + "@19 seq19\n" + "@20 seq20\n" + "@21 seq21\n" + "@22 seq22\n" + "@23 seq23\n" + "@24 seq24\n" + "@25 seq25\n" + "@26 seq26\n" + "@27 seq27\n" + "@28 seq28\n" + "@29 seq29\n" + "@30 seq30\n" + "@31 seq31\n" + "@32 seq32\n" + "@33 seq33\n" + "@34 seq34\n" + "@35 seq35\n" + "@36 seq36\n" + "@37 seq37\n" + "@38 seq38\n" + "@39 seq39\n" + "@40 seq40\n" + "@41 seq41\n" + "@42 seq42\n" + "@43 seq43\n" + "@44 seq44\n" + "@45 seq45\n" + "@46 seq46\n" + "@47 seq47\n" + "@48 seq48\n" + "@49 seq49\n" + "@50 seq50\n" + "@51 seq51\n" + "@52 seq52\n" + "@53 seq53\n" + "@54 seq54\n" + "@55 seq55\n" + "@56 seq56\n" + "@57 seq57\n" + "@58 seq58\n" + "@59 seq59\n" + "@60 seq60\n" + "@61 seq61\n" + "@62 seq62\n" + "@63 seq63\n" + "@64 seq64\n" + "@65 seq65\n" + "@66 seq66\n" + "@67 seq67\n" + "@68 seq68\n" + "@69 seq69\n" + "@70 seq70\n" + "@71 seq71\n" + "@72 seq72\n" + "@73 seq73\n" + "@74 seq74\n" + "@75 seq75\n" + "@76 seq76\n" + "@77 seq77\n" + "@78 seq78\n" + "@79 seq79\n" + "@80 seq80\n" + "@81 seq81\n" + "@82 seq82\n" + "@83 seq83\n" + "@84 seq84\n" + "@85 seq85\n" + "@86 seq86\n" + "@87 seq87\n" + "@88 seq88\n" + "@89 seq89\n" + "@90 seq90\n" + "@91 seq91\n" + "@92 seq92\n" + "@93 seq93\n" + "@94 seq94\n" + "@95 seq95\n" + "@CHOPPER_USER_BINS_END\n" + "@CHOPPER_CONFIG\n" + "@{\n" + "@ \"chopper_config\": {\n" + "@ \"version\": 2,\n" + "@ \"data_file\": {\n" + "@ \"value0\": \"\"\n" + "@ },\n" + "@ \"debug\": false,\n" + "@ \"sketch_directory\": {\n" + "@ \"value0\": \"\"\n" + "@ },\n" + "@ \"k\": 19,\n" + "@ \"window_size\": 19,\n" + "@ \"disable_sketch_output\": true,\n" + "@ \"precomputed_files\": false,\n" + "@ \"output_filename\": {\n" + "@ \"value0\": \"" + + layout_file.string() + + "\"\n" + "@ },\n" + "@ \"determine_best_tmax\": true,\n" + "@ \"force_all_binnings\": false\n" + "@ }\n" + "@}\n" + "@CHOPPER_CONFIG_END\n" + "@HIBF_CONFIG\n" + "@{\n" + "@ \"hibf_config\": {\n" + "@ \"version\": 3,\n" + "@ \"number_of_user_bins\": 96,\n" + "@ \"number_of_hash_functions\": 2,\n" + "@ \"maximum_fpr\": 0.05,\n" + "@ \"relaxed_fpr\": 0.3,\n" + "@ \"threads\": 1,\n" + "@ \"sketch_bits\": 12,\n" + "@ \"tmax\": 64,\n" + "@ \"empty_bin_fraction\": 0.0,\n" + "@ \"track_occupancy\": false,\n" + "@ \"alpha\": 1.2,\n" + "@ \"max_rearrangement_ratio\": 0.5,\n" + "@ \"disable_estimate_union\": true,\n" + "@ \"disable_rearrangement\": true\n" + "@ }\n" + "@}\n" + "@HIBF_CONFIG_END\n" + "#TOP_LEVEL_IBF fullest_technical_bin_idx:63\n" + "#LOWER_LEVEL_IBF_0 fullest_technical_bin_idx:0\n" + "#LOWER_LEVEL_IBF_1 fullest_technical_bin_idx:52\n" + "#LOWER_LEVEL_IBF_2 fullest_technical_bin_idx:52\n" + "#LOWER_LEVEL_IBF_3 fullest_technical_bin_idx:52\n" + "#LOWER_LEVEL_IBF_4 fullest_technical_bin_idx:31\n" + "#LOWER_LEVEL_IBF_5 fullest_technical_bin_idx:0\n" + "#LOWER_LEVEL_IBF_6 fullest_technical_bin_idx:0\n" + "#LOWER_LEVEL_IBF_7 fullest_technical_bin_idx:0\n" + "#LOWER_LEVEL_IBF_8 fullest_technical_bin_idx:0\n" + "#LOWER_LEVEL_IBF_9 fullest_technical_bin_idx:0\n" + "#USER_BIN_IDX\tTECHNICAL_BIN_INDICES\tNUMBER_OF_TECHNICAL_BINS\n" + "1\t0;0\t1;32\n" + "0\t0;32\t1;32\n" + "6\t1;0\t1;13\n" + "5\t1;13\t1;13\n" + "4\t1;26\t1;13\n" + "3\t1;39\t1;13\n" + "2\t1;52\t1;12\n" + "11\t2;0\t1;13\n" + "10\t2;13\t1;13\n" + "9\t2;26\t1;13\n" + "8\t2;39\t1;13\n" + "7\t2;52\t1;12\n" + "16\t3;0\t1;13\n" + "15\t3;13\t1;13\n" + "14\t3;26\t1;13\n" + "13\t3;39\t1;13\n" + "12\t3;52\t1;12\n" + "21\t4;0\t1;16\n" + "20\t4;16\t1;15\n" + "19\t4;31\t1;11\n" + "18\t4;42\t1;11\n" + "17\t4;53\t1;11\n" + "25\t5;0\t1;16\n" + "24\t5;16\t1;16\n" + "23\t5;32\t1;16\n" + "22\t5;48\t1;16\n" + "29\t6;0\t1;16\n" + "28\t6;16\t1;16\n" + "27\t6;32\t1;16\n" + "26\t6;48\t1;16\n" + "33\t7;0\t1;16\n" + "32\t7;16\t1;16\n" + "31\t7;32\t1;16\n" + "30\t7;48\t1;16\n" + "37\t8;0\t1;16\n" + "36\t8;16\t1;16\n" + "35\t8;32\t1;16\n" + "34\t8;48\t1;16\n" + "41\t9;0\t1;18\n" + "40\t9;18\t1;18\n" + "39\t9;36\t1;14\n" + "38\t9;50\t1;14\n" + "42\t10\t1\n" + "43\t11\t1\n" + "44\t12\t1\n" + "45\t13\t1\n" + "46\t14\t1\n" + "47\t15\t1\n" + "48\t16\t1\n" + "49\t17\t1\n" + "50\t18\t1\n" + "51\t19\t1\n" + "52\t20\t1\n" + "53\t21\t1\n" + "54\t22\t1\n" + "55\t23\t1\n" + "56\t24\t1\n" + "57\t25\t1\n" + "58\t26\t1\n" + "59\t27\t1\n" + "60\t28\t1\n" + "61\t29\t1\n" + "62\t30\t1\n" + "63\t31\t1\n" + "64\t32\t1\n" + "65\t33\t1\n" + "66\t34\t1\n" + "67\t35\t1\n" + "68\t36\t1\n" + "69\t37\t1\n" + "70\t38\t1\n" + "71\t39\t1\n" + "72\t40\t1\n" + "73\t41\t1\n" + "74\t42\t1\n" + "75\t43\t1\n" + "76\t44\t1\n" + "77\t45\t1\n" + "78\t46\t1\n" + "79\t47\t1\n" + "80\t48\t1\n" + "81\t49\t1\n" + "82\t50\t1\n" + "83\t51\t1\n" + "84\t52\t1\n" + "85\t53\t1\n" + "86\t54\t1\n" + "87\t55\t1\n" + "88\t56\t1\n" + "89\t57\t1\n" + "90\t58\t1\n" + "91\t59\t1\n" + "92\t60\t1\n" + "93\t61\t1\n" + "94\t62\t1\n" + "95\t63\t1\n"}; + std::string const actual_file{string_from_file(layout_file)}; + EXPECT_EQ(actual_file, expected_file) << actual_file; +} + +TEST(execute_estimation_test, many_ubs_fast_layout) +{ + seqan3::test::tmp_directory tmp_dir{}; + std::filesystem::path const layout_file{tmp_dir.path() / "layout.tsv"}; + std::filesystem::path const stats_file{layout_file.string() + ".stats"}; + + std::vector> many_filenames; + + for (size_t i{0}; i < 96u; ++i) + many_filenames.push_back({seqan3::detail::to_string("seq", i)}); + + // Creates sizes of the following series + // [801,802,...,820,922,923,...,941,1043,1044,...,1062,1164,1165,...,1183,1285,1286,...,1300] + // See also https://godbolt.org/z/9517eaaaG + auto simulated_input = [&](size_t const num, seqan::hibf::insert_iterator it) + { + size_t const desired_kmer_count = 101 * ((num + 20) / 20) + num + 700; + for (auto hash : std::views::iota(0u, desired_kmer_count)) + it = hash; + }; + + chopper::configuration config{}; + config.fast_layout = true; + config.determine_best_tmax = true; + config.output_filename = layout_file; + config.disable_sketch_output = true; + config.hibf_config.tmax = 1024; + config.hibf_config.input_fn = simulated_input; + config.hibf_config.number_of_user_bins = many_filenames.size(); + config.hibf_config.disable_estimate_union = true; // also disables rearrangement + + std::vector sketches; + std::vector minHash_sketches{}; + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); - chopper::layout::execute(config, many_filenames, sketches); + chopper::layout::execute(config, many_filenames, sketches, minHash_sketches); ASSERT_TRUE(std::filesystem::exists(stats_file)); @@ -429,9 +807,10 @@ TEST(execute_estimation_test, many_ubs_force_all) config.hibf_config.disable_estimate_union = true; // also disables rearrangement std::vector sketches; - seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches); + std::vector minHash_sketches{}; + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); - chopper::layout::execute(config, many_filenames, sketches); + chopper::layout::execute(config, many_filenames, sketches, minHash_sketches); ASSERT_TRUE(std::filesystem::exists(stats_file)); @@ -519,9 +898,10 @@ TEST(execute_estimation_test, with_rearrangement) config.hibf_config.number_of_user_bins = filenames.size(); std::vector sketches; - seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches); + std::vector minHash_sketches{}; + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); - chopper::layout::execute(config, filenames, sketches); + chopper::layout::execute(config, filenames, sketches, minHash_sketches); ASSERT_TRUE(std::filesystem::exists(stats_file)); diff --git a/test/api/layout/fast_layout_cluster_test.cpp b/test/api/layout/fast_layout_cluster_test.cpp new file mode 100644 index 00000000..928eb708 --- /dev/null +++ b/test/api/layout/fast_layout_cluster_test.cpp @@ -0,0 +1,72 @@ +#include // for Test, TestInfo, EXPECT_EQ, Message, TEST, TestPartResult + +#include // for size_t +#include // for operator<<, char_traits, basic_ostream, basic_stringstream, strings... +#include // for allocator, string +#include // for operator<< +#include // for vector + +#include + +TEST(Cluster_test, ctor_from_id) +{ + size_t user_bin_idx{5}; + chopper::layout::Cluster const cluster{user_bin_idx}; + + EXPECT_EQ(cluster.id(), user_bin_idx); + EXPECT_FALSE(cluster.empty()); + EXPECT_EQ(cluster.size(), 1u); + ASSERT_EQ(cluster.contained_user_bins().size(), 1u); + EXPECT_EQ(cluster.contained_user_bins()[0], user_bin_idx); + EXPECT_TRUE(cluster.is_valid(user_bin_idx)); +} + +TEST(Cluster_test, move_to) +{ + size_t user_bin_idx1{5}; + size_t user_bin_idx2{7}; + chopper::layout::Cluster cluster1{user_bin_idx1}; + chopper::layout::Cluster cluster2{user_bin_idx2}; + + EXPECT_TRUE(cluster1.is_valid(user_bin_idx1)); + EXPECT_TRUE(cluster2.is_valid(user_bin_idx2)); + + cluster2.move_to(cluster1); + + // cluster1 now contains user bins 5 and 7 + EXPECT_EQ(cluster1.size(), 2u); + ASSERT_EQ(cluster1.contained_user_bins().size(), 2u); + EXPECT_EQ(cluster1.contained_user_bins()[0], user_bin_idx1); + EXPECT_EQ(cluster1.contained_user_bins()[1], user_bin_idx2); + + // cluster 2 is empty + EXPECT_TRUE(cluster2.has_been_moved()); + EXPECT_TRUE(cluster2.empty()); + EXPECT_EQ(cluster2.size(), 0u); + EXPECT_EQ(cluster2.contained_user_bins().size(), 0u); + EXPECT_EQ(cluster2.moved_to_cluster_id(), cluster1.id()); + + // both should still be valid + EXPECT_TRUE(cluster1.is_valid(user_bin_idx1)); + EXPECT_TRUE(cluster2.is_valid(user_bin_idx2)); +} + +TEST(LSH_find_representative_cluster_test, cluster_one_move) +{ + std::vector clusters{chopper::layout::Cluster{0}, chopper::layout::Cluster{1}}; + clusters[1].move_to(clusters[0]); + + EXPECT_EQ(chopper::layout::LSH_find_representative_cluster(clusters, clusters[1].id()), clusters[0].id()); +} + +TEST(LSH_find_representative_cluster_test, cluster_two_moves) +{ + std::vector clusters{chopper::layout::Cluster{0}, + chopper::layout::Cluster{1}, + chopper::layout::Cluster{2}}; + clusters[2].move_to(clusters[1]); + clusters[1].move_to(clusters[0]); + + EXPECT_EQ(chopper::layout::LSH_find_representative_cluster(clusters, clusters[1].id()), clusters[0].id()); + EXPECT_EQ(chopper::layout::LSH_find_representative_cluster(clusters, clusters[2].id()), clusters[0].id()); +} diff --git a/test/api/layout/fast_layout_find_bins_to_be_split_test.cpp b/test/api/layout/fast_layout_find_bins_to_be_split_test.cpp new file mode 100644 index 00000000..a96db447 --- /dev/null +++ b/test/api/layout/fast_layout_find_bins_to_be_split_test.cpp @@ -0,0 +1,66 @@ +// -------------------------------------------------------------------------------------------------- +// Copyright (c) 2006-2023, Knut Reinert & Freie Universität Berlin +// Copyright (c) 2016-2023, Knut Reinert & MPI für molekulare Genetik +// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License +// shipped with this file and also available at: https://github.com/seqan/chopper/blob/main/LICENSE.md +// -------------------------------------------------------------------------------------------------- + +#include // for Test, TestInfo, EXPECT_EXIT, TEST + +#include // for exit +#include // for iota +#include // for alarm +#include // for vector + +#include + +TEST(find_bins_to_be_split_test, dissertation_mehringer_tiny_example) +{ + std::vector const cardinalities{120, 60, 30, 30}; + std::vector const sorted_positions{0, 1, 2, 3}; + + // t_s = number of split techincal bins + auto const [idx, t_s] = chopper::layout::find_bins_to_be_split(sorted_positions, cardinalities, /*t=*/45, 3); + EXPECT_EQ(t_s, 2); // two split bins + EXPECT_EQ(idx, 1); // points to the end of the range of indices to split -> only user bin 0 (cardinality 120) +} + +TEST(find_bins_to_be_split_test, all_bins_split) +{ + std::vector const cardinalities(60, 1100); // 60 times a value of 1100 + std::vector sorted_positions(cardinalities.size()); + std::iota(sorted_positions.begin(), sorted_positions.end(), 0); + + // t_s = number of split techincal bins + auto const [idx, t_s] = chopper::layout::find_bins_to_be_split(sorted_positions, cardinalities, /*t=*/1000, 63); + EXPECT_EQ(t_s, 63); // 63 split bins + EXPECT_EQ(idx, 60); // all user bins shall be split +} + +TEST(find_bins_to_be_split_test, small_number_edge_case) +{ + // small number edge case: the `std::max(threshold + 1, ... )` catches here + // but increasing by 1 results in none of the bins split. They will be "merged" then + // which is not optimal but fine. The merging algorithm will distribute each user bin + // into one technical bin and 4 technical bins will be empty. With these small numbers + // that is fine. + std::vector const cardinalities(60, 11); // 60 times a value of 11 + std::vector sorted_positions(cardinalities.size()); + std::iota(sorted_positions.begin(), sorted_positions.end(), 0); + + // t_s = number of split techincal bins + auto const [idx, t_s] = chopper::layout::find_bins_to_be_split(sorted_positions, cardinalities, /*t=*/10, 63); + EXPECT_EQ(t_s, 0); // 63 split bins + EXPECT_EQ(idx, 0); // no user bins shall be split +} + +TEST(find_bins_to_be_split_test, threshold_zero_death_in_debug) +{ +#ifdef NDEBUG + GTEST_SKIP() << "Debug-only test"; +#endif + std::vector const cardinalities{10, 20, 30}; + std::vector const sorted_positions{0, 1, 2}; + + EXPECT_DEATH(chopper::layout::find_bins_to_be_split(sorted_positions, cardinalities, /*t=*/0, 2), "threshold > 0"); +} diff --git a/test/api/layout/hibf_statistics_test.cpp b/test/api/layout/hibf_statistics_test.cpp index 26c1dcb6..5ce53a9e 100644 --- a/test/api/layout/hibf_statistics_test.cpp +++ b/test/api/layout/hibf_statistics_test.cpp @@ -133,11 +133,12 @@ TEST(execute_test, chopper_layout_statistics) .disable_estimate_union = true /* also disable rearrangement */}}; std::vector sketches; - seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches); + std::vector minHash_sketches{}; + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); testing::internal::CaptureStdout(); testing::internal::CaptureStderr(); - chopper::layout::execute(config, many_filenames, sketches); + chopper::layout::execute(config, many_filenames, sketches, minHash_sketches); std::string layout_result_stdout = testing::internal::GetCapturedStdout(); std::string layout_result_stderr = testing::internal::GetCapturedStderr(); @@ -161,6 +162,66 @@ TEST(execute_test, chopper_layout_statistics) EXPECT_EQ(layout_result_stderr, std::string{}); } +TEST(execute_test, chopper_layout_statistics_fast_layout) +{ + seqan3::test::tmp_directory tmp_dir{}; + std::filesystem::path const layout_file{tmp_dir.path() / "layout.tsv"}; + + std::vector> many_filenames; + + for (size_t i{0}; i < 96u; ++i) + many_filenames.push_back({seqan3::detail::to_string("seq", i)}); + + // Creates sizes of the following series + // [801,802,...,820,922,923,...,941,1043,1044,...,1062,1164,1165,...,1183,1285,1286,...,1300] + // See also https://godbolt.org/z/9517eaaaG + auto simulated_input = [&](size_t const num, seqan::hibf::insert_iterator it) + { + size_t const desired_kmer_count = 101 * ((num + 20) / 20) + num + 700; + for (auto hash : std::views::iota(0u, desired_kmer_count)) + it = hash; + }; + + chopper::configuration config{.fast_layout = true, + .data_file = "not needed", + .output_filename = layout_file.c_str(), + .disable_sketch_output = true, + .output_verbose_statistics = true, + .hibf_config = {.input_fn = simulated_input, + .number_of_user_bins = many_filenames.size(), + .tmax = 64, + .disable_estimate_union = true /* also disable rearrangement */}}; + + std::vector sketches; + std::vector minHash_sketches{}; + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); + + testing::internal::CaptureStdout(); + testing::internal::CaptureStderr(); + chopper::layout::execute(config, many_filenames, sketches, minHash_sketches); + std::string layout_result_stdout = testing::internal::GetCapturedStdout(); + std::string layout_result_stderr = testing::internal::GetCapturedStderr(); + + std::string expected_cout = + R"expected_cout(## ### Notation ### +## X-IBF = An IBF with X number of bins. +## X-HIBF = An HIBF with tmax = X, e.g a maximum of X technical bins on each level. +## ### Column Description ### +## tmax : The maximum number of technical bin on each level +## c_tmax : The technical extra cost of querying an tmax-IBF, compared to 64-IBF +## l_tmax : The estimated query cost for an tmax-HIBF, compared to an 64-HIBF +## m_tmax : The estimated memory consumption for an tmax-HIBF, compared to an 64-HIBF +## (l*m)_tmax : Computed by l_tmax * m_tmax +## size : The expected total size of an tmax-HIBF +## uncorr_size : The expected size of an tmax-HIBF without FPR correction +# tmax c_tmax l_tmax m_tmax (l*m)_tmax size uncorr_size level num_ibfs level_size level_size_no_corr total_num_tbs avg_num_tbs split_tb_percentage max_split_tb avg_split_tb max_factor avg_factor +64 1.00 1.24 1.00 1.24 672.6KiB 2.0MiB :0:1:2 :1:1:1 :624.8KiB:47.8KiB:0Bytes :1.9MiB:130.0KiB:0Bytes :64:64:0 :64:64:0 :98.44:98.44:-nan :1:3:- :1.00:2.03:- :1.00:1.81:- :1.00:1.47:- +)expected_cout"; + + EXPECT_EQ(layout_result_stdout, expected_cout) << layout_result_stdout; + EXPECT_EQ(layout_result_stderr, std::string{}); +} + TEST(execute_test, chopper_layout_statistics_determine_best_bins) { seqan3::test::tmp_directory tmp_dir{}; @@ -190,9 +251,10 @@ TEST(execute_test, chopper_layout_statistics_determine_best_bins) .disable_estimate_union = true /* also disable rearrangement */}}; std::vector sketches; - seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches); + std::vector minHash_sketches{}; + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); - chopper::layout::execute(config, filenames, sketches); + chopper::layout::execute(config, filenames, sketches, minHash_sketches); std::string expected_cout = R"expected_cout(## ### Parameters ### diff --git a/test/api/layout/partition_user_bins_test.cpp b/test/api/layout/partition_user_bins_test.cpp new file mode 100644 index 00000000..161a7f99 --- /dev/null +++ b/test/api/layout/partition_user_bins_test.cpp @@ -0,0 +1,178 @@ +// -------------------------------------------------------------------------------------------------- +// Copyright (c) 2006-2023, Knut Reinert & Freie Universität Berlin +// Copyright (c) 2016-2023, Knut Reinert & MPI für molekulare Genetik +// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License +// shipped with this file and also available at: https://github.com/seqan/chopper/blob/main/LICENSE.md +// -------------------------------------------------------------------------------------------------- + +#include // for Test, TestInfo, EXPECT_EQ, TEST + +#include // for size_t +#include // for uint64_t +#include // for exit +#include // for iota +#include // for to_string +#include // for vector + +#include +#include + +#include +#include +#include + +namespace +{ + +// splitmix64 finaliser: turns consecutive integers into well-distributed hashes. HyperLogLog and the MinHash buckets +// (hash & 15, each needs 40 values) both assume uniformly distributed hashes, which plain consecutive integers are not. +uint64_t scramble(uint64_t x) +{ + x += 0x9e3779b97f4a7c15ULL; + x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9ULL; + x = (x ^ (x >> 27)) * 0x94d049bb133111ebULL; + return x ^ (x >> 31); +} + +// User bin `i` consists of `kmer_counts[i]` distinct hashes drawn from content `content_ids[i]`. +// User bins with the same content id share their first min(kmer_counts) hashes, all others are disjoint. +// If `content_ids` is empty, every user bin has its own content. +std::vector> run_partition_user_bins(std::vector const & kmer_counts, + size_t const tmax, + std::vector content_ids = {}) +{ + if (content_ids.empty()) + { + content_ids.resize(kmer_counts.size()); + std::iota(content_ids.begin(), content_ids.end(), 0u); + } + + chopper::configuration config{}; + config.hibf_config.tmax = tmax; + config.hibf_config.number_of_user_bins = kmer_counts.size(); + config.hibf_config.input_fn = [&kmer_counts, &content_ids](size_t const ub, seqan::hibf::insert_iterator it) + { + uint64_t const offset = static_cast(content_ids[ub]) << 32; + for (uint64_t i = 0; i < kmer_counts[ub]; ++i) + it = scramble(offset + i); + }; + + std::vector sketches{}; + std::vector minHash_sketches{}; + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); + + std::vector cardinalities(kmer_counts.size()); + for (size_t i = 0; i < sketches.size(); ++i) + cardinalities[i] = sketches[i].estimate(); + + std::vector positions(kmer_counts.size()); + std::iota(positions.begin(), positions.end(), 0u); + + std::vector> partitions(tmax); + chopper::layout::partition_user_bins(config, positions, cardinalities, sketches, minHash_sketches, partitions); + + return partitions; +} + +} // namespace + +// A user bin is split if it is assigned to more than one technical bin. +// A technical bin is merged if it contains more than one user bin. + +TEST(partition_user_bins_test, only_split_bins) +{ + // 4 large user bins for 8 technical bins: every user bin exceeds the split threshold. + auto const partitions = run_partition_user_bins(std::vector(4, 10'000), /*tmax*/ 8); + + // each user bin is split into two technical bins + std::vector> expected_partitions{{3}, {3}, {1}, {1}, {2}, {2}, {0}, {0}}; + + ASSERT_EQ(partitions.size(), expected_partitions.size()); + for (size_t tb = 0; tb < partitions.size(); ++tb) + EXPECT_EQ(partitions[tb], expected_partitions[tb]) << "technical bin " << tb << " is not correct"; +} + +TEST(partition_user_bins_test, initial_split_threshold_small) +{ + std::vector kmer_counts(2, 200'000); + kmer_counts.resize(12, 1000); + std::vector content_ids{0, 0, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11}; // user bins 0 and 1 are identical + + auto const partitions = run_partition_user_bins(kmer_counts, /*tmax*/ 7, content_ids); + + // ToDo this partitioning is not perfect + std::vector> expected_partitions{{8, 4, 6, 10, 9, 3}, {5, 7}, {11, 2}, {0}, {0}, {1}, {1}}; + + ASSERT_EQ(partitions.size(), expected_partitions.size()); + for (size_t tb = 0; tb < partitions.size(); ++tb) + EXPECT_EQ(partitions[tb], expected_partitions[tb]) << "technical bin " << tb << " is not correct"; +} + +TEST(partition_user_bins_test, only_merged_bins) +{ + std::vector kmer_counts(200, 2'000); + // 20 small user bins for 2 technical bins: no use + auto const partitions = run_partition_user_bins(kmer_counts, 2); + + // since all user bins have the same cardinality, both partitions should contain equally many + // SInce not cardinalities but HLL sketches are used, there is some noise + // Thats why there are not exactly equal + ASSERT_EQ(partitions.size(), 2); + EXPECT_EQ(partitions[0].size(), 102); + EXPECT_EQ(partitions[1].size(), 98); +} + +TEST(partition_user_bins_test, another_edge_case) +{ + // one very large user bin and one very small one. + std::vector const kmer_counts{400'000, 2'000}; + + auto const partitions = run_partition_user_bins(kmer_counts, 64); + + std::vector> expected_partitions{{8, 4, 6, 10, 9, 3}, {5, 7}, {11, 2}, {0}, {0}, {1}, {1}}; + + ASSERT_EQ(partitions.size(), 64); + + ASSERT_EQ(partitions[0].size(), 1); + EXPECT_EQ(partitions[0][0], 1); + for (size_t tb = 1; tb < partitions.size(); ++tb) + { + ASSERT_EQ(partitions[tb].size(), 1); + EXPECT_EQ(partitions[tb][0], 0) << "technical bin " << tb << " is not correct"; + } +} + +// KI tests that failed and caught some bugs that are now fixed: + +TEST(partition_user_bins_test, one_large_cluster) +{ + // This case simulates that there are more available merged bins than clusters + // with only one large cluster + std::vector kmer_counts{50'635, 16'021, 1'734, 305'923, 1'844, 7'087}; + std::vector content_ids{0, 0, 0, 0, 0, 0}; + auto const partitions = run_partition_user_bins(kmer_counts, 64, content_ids); + + ASSERT_EQ(partitions.size(), 64); +} + +TEST(partition_user_bins_test, one_large_one_small_cluster) +{ + // This case simulates that there are more available merged bins than clusters + // with two clusters, of which one only has size 1 + std::vector kmer_counts{4'169, 14'011, 4'311, 2'016, 239'260, 16'111}; + std::vector content_ids{0, 0, 0, 0, 1, 0}; + auto const partitions = run_partition_user_bins(kmer_counts, 64, content_ids); + + ASSERT_EQ(partitions.size(), 64); +} + +TEST(partition_user_bins_test, several_clusters) +{ + // This case simulates that there are more available merged bins than clusters + // with several clusters + std::vector kmer_counts{5'829, 1'875, 8'912, 254'423, 6'445, 188'171, 7'621, 7'127, 7'579, 11'915, 344'000}; + std::vector content_ids{1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1}; + auto const partitions = run_partition_user_bins(kmer_counts, 64, content_ids); + + ASSERT_EQ(partitions.size(), 64); +} diff --git a/test/cli/cli_chopper_pipeline_test.cpp b/test/cli/cli_chopper_pipeline_test.cpp index 308299ce..1ca7a53d 100644 --- a/test/cli/cli_chopper_pipeline_test.cpp +++ b/test/cli/cli_chopper_pipeline_test.cpp @@ -222,3 +222,108 @@ TEST_F(cli_test, chopper_layout2) std::string const actual_file{string_from_file(binning_filename)}; EXPECT_EQ(actual_file, expected_file); } + +TEST_F(cli_test, chopper_layout2_fast_layout) +{ + std::string const seq1_filename = data("seq1.fa"); + std::string const seq2_filename = data("seq2.fa"); + std::string const seq3_filename = data("seq3.fa"); + std::string const seq4_filename = data("small.fa"); + seqan3::test::tmp_directory tmp_dir{}; + std::filesystem::path const taxa_filename{tmp_dir.path() / "data.tsv"}; + std::filesystem::path const binning_filename{tmp_dir.path() / "output.binning"}; + + // we need to have filenames from the user + { + std::ofstream fout{taxa_filename}; + fout << seq1_filename << '\n' << seq2_filename << '\n' << seq3_filename << '\n' << seq4_filename << '\n'; + } + + cli_test_result result = execute_app("chopper", + "--threads", + "2", + "--sketch-bits", + "12", + "--fast-layout", + "--input", + taxa_filename.c_str(), + "--tmax", + "64", + "--output", + binning_filename.c_str()); + + EXPECT_EQ(result.exit_code, 0); + EXPECT_EQ(result.out, std::string{}); + EXPECT_EQ(result.err, std::string{}); + + std::string const expected_file{"@CHOPPER_USER_BINS\n" + "@0 " + + seq1_filename + + "\n" + "@1 " + + seq2_filename + + "\n" + "@2 " + + seq3_filename + + "\n" + "@3 " + + seq4_filename + + "\n" + "@CHOPPER_USER_BINS_END\n" + "@CHOPPER_CONFIG\n" + "@{\n" + "@ \"chopper_config\": {\n" + "@ \"version\": 2,\n" + "@ \"data_file\": {\n" + "@ \"value0\": \"" + + taxa_filename.string() + + "\"\n" + "@ },\n" + "@ \"debug\": false,\n" + "@ \"sketch_directory\": {\n" + "@ \"value0\": \"\"\n" + "@ },\n" + "@ \"k\": 19,\n" + "@ \"window_size\": 19,\n" + "@ \"disable_sketch_output\": true,\n" + "@ \"precomputed_files\": false,\n" + "@ \"output_filename\": {\n" + "@ \"value0\": \"" + + binning_filename.string() + + "\"\n" + "@ },\n" + "@ \"determine_best_tmax\": false,\n" + "@ \"force_all_binnings\": false\n" + "@ }\n" + "@}\n" + "@CHOPPER_CONFIG_END\n" + "@HIBF_CONFIG\n" + "@{\n" + "@ \"hibf_config\": {\n" + "@ \"version\": 3,\n" + "@ \"number_of_user_bins\": 4,\n" + "@ \"number_of_hash_functions\": 2,\n" + "@ \"maximum_fpr\": 0.05,\n" + "@ \"relaxed_fpr\": 0.3,\n" + "@ \"threads\": 2,\n" + "@ \"sketch_bits\": 12,\n" + "@ \"tmax\": 64,\n" + "@ \"empty_bin_fraction\": 0.0,\n" + "@ \"track_occupancy\": false,\n" + "@ \"alpha\": 1.2,\n" + "@ \"max_rearrangement_ratio\": 0.5,\n" + "@ \"disable_estimate_union\": false,\n" + "@ \"disable_rearrangement\": false\n" + "@ }\n" + "@}\n" + "@HIBF_CONFIG_END\n" + "#TOP_LEVEL_IBF fullest_technical_bin_idx:0\n" + "#USER_BIN_IDX\tTECHNICAL_BIN_INDICES\tNUMBER_OF_TECHNICAL_BINS\n" + "0\t56\t8\n" + "1\t39\t9\n" + "2\t48\t8\n" + "3\t0\t39\n"}; + + std::string const actual_file{string_from_file(binning_filename)}; + EXPECT_EQ(actual_file, expected_file); +} diff --git a/test/cli/cli_output_sketches.cpp b/test/cli/cli_output_sketches.cpp index 6d0662db..2d47e019 100644 --- a/test/cli/cli_output_sketches.cpp +++ b/test/cli/cli_output_sketches.cpp @@ -96,7 +96,7 @@ TEST_F(cli_test, chopper_layout) EXPECT_EQ(sin.filenames.size(), 3); EXPECT_EQ(sin.hll_sketches.size(), 3); - EXPECT_EQ(sin.minHash_sketches.size(), 0); // currently, no minhash sketches are needed in chopper layout + EXPECT_EQ(sin.minHash_sketches.size(), 3); // currently, no minhash sketches are needed in chopper layout EXPECT_EQ(sin.filenames[0][0], data("seq1.fa").string()); EXPECT_EQ(sin.filenames[1][0], data("seq2.fa").string()); From 23b60fc83cfe1ddd58f38031f047d2701e600e04 Mon Sep 17 00:00:00 2001 From: Svenja Mehringer Date: Tue, 29 Sep 2026 12:30:30 +0200 Subject: [PATCH 02/35] [HEADER TEST] Add . --- include/chopper/layout/fast_layout_cluster.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/include/chopper/layout/fast_layout_cluster.hpp b/include/chopper/layout/fast_layout_cluster.hpp index f02a14e9..f6654af9 100644 --- a/include/chopper/layout/fast_layout_cluster.hpp +++ b/include/chopper/layout/fast_layout_cluster.hpp @@ -8,6 +8,7 @@ #pragma once #include +#include #include #include From 0ee07d285a439952f0fba36820733c4fbf84c667 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 18:32:58 +0200 Subject: [PATCH 03/35] fix: missing include --- include/chopper/layout/fast_layout_cluster.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/chopper/layout/fast_layout_cluster.hpp b/include/chopper/layout/fast_layout_cluster.hpp index f6654af9..e08566b4 100644 --- a/include/chopper/layout/fast_layout_cluster.hpp +++ b/include/chopper/layout/fast_layout_cluster.hpp @@ -9,6 +9,7 @@ #include #include +#include #include #include @@ -131,4 +132,4 @@ inline size_t LSH_find_representative_cluster(std::vector const & clust return current_id; } -} // namespace chopper::layout \ No newline at end of file +} // namespace chopper::layout From 73b705f65e2a1c21e888759e52ea754ddff50ee0 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 18:37:16 +0200 Subject: [PATCH 04/35] fix: partition test clang --- test/api/layout/partition_user_bins_test.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/api/layout/partition_user_bins_test.cpp b/test/api/layout/partition_user_bins_test.cpp index 161a7f99..1db95fd0 100644 --- a/test/api/layout/partition_user_bins_test.cpp +++ b/test/api/layout/partition_user_bins_test.cpp @@ -118,8 +118,13 @@ TEST(partition_user_bins_test, only_merged_bins) // SInce not cardinalities but HLL sketches are used, there is some noise // Thats why there are not exactly equal ASSERT_EQ(partitions.size(), 2); +#ifdef _LIBCPP_VERSION + EXPECT_EQ(partitions[0].size(), 99); + EXPECT_EQ(partitions[1].size(), 101); +#else EXPECT_EQ(partitions[0].size(), 102); EXPECT_EQ(partitions[1].size(), 98); +#endif } TEST(partition_user_bins_test, another_edge_case) From a4370c38c929b7169a0f37636266598408dae47f Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 18:55:55 +0200 Subject: [PATCH 05/35] fix: just hardoce NaN string 0/0.0 is nan on gcc, but -nan on clang... --- src/layout/hibf_statistics.cpp | 2 ++ test/api/layout/hibf_statistics_test.cpp | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/layout/hibf_statistics.cpp b/src/layout/hibf_statistics.cpp index df86775c..f966896d 100644 --- a/src/layout/hibf_statistics.cpp +++ b/src/layout/hibf_statistics.cpp @@ -127,6 +127,8 @@ void hibf_statistics::print_summary_to(size_t & t_max_64_memory, std::ostream & // go through each level and collect and output the statistics auto to_string_with_precision = [](auto num) { + if (std::isnan(num)) + return "NaN"; std::stringstream ss; ss << std::fixed << std::setprecision(2) << num; return ss.str(); diff --git a/test/api/layout/hibf_statistics_test.cpp b/test/api/layout/hibf_statistics_test.cpp index 5ce53a9e..92bf9745 100644 --- a/test/api/layout/hibf_statistics_test.cpp +++ b/test/api/layout/hibf_statistics_test.cpp @@ -215,7 +215,7 @@ TEST(execute_test, chopper_layout_statistics_fast_layout) ## size : The expected total size of an tmax-HIBF ## uncorr_size : The expected size of an tmax-HIBF without FPR correction # tmax c_tmax l_tmax m_tmax (l*m)_tmax size uncorr_size level num_ibfs level_size level_size_no_corr total_num_tbs avg_num_tbs split_tb_percentage max_split_tb avg_split_tb max_factor avg_factor -64 1.00 1.24 1.00 1.24 672.6KiB 2.0MiB :0:1:2 :1:1:1 :624.8KiB:47.8KiB:0Bytes :1.9MiB:130.0KiB:0Bytes :64:64:0 :64:64:0 :98.44:98.44:-nan :1:3:- :1.00:2.03:- :1.00:1.81:- :1.00:1.47:- +64 1.00 1.24 1.00 1.24 672.6KiB 2.0MiB :0:1:2 :1:1:1 :624.8KiB:47.8KiB:0Bytes :1.9MiB:130.0KiB:0Bytes :64:64:0 :64:64:0 :98.44:98.44:NaN :1:3:- :1.00:2.03:- :1.00:1.81:- :1.00:1.47:- )expected_cout"; EXPECT_EQ(layout_result_stdout, expected_cout) << layout_result_stdout; From 1d937682ae973468cf8adfae71db78e7b9d549ea Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 19:19:49 +0200 Subject: [PATCH 06/35] fix: lambda return deduction --- src/layout/hibf_statistics.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/layout/hibf_statistics.cpp b/src/layout/hibf_statistics.cpp index f966896d..ca9a9679 100644 --- a/src/layout/hibf_statistics.cpp +++ b/src/layout/hibf_statistics.cpp @@ -125,7 +125,7 @@ void hibf_statistics::print_summary_to(size_t & t_max_64_memory, std::ostream & size_t total_size_no_corr{}; // go through each level and collect and output the statistics - auto to_string_with_precision = [](auto num) + auto to_string_with_precision = [](auto num) -> std::string { if (std::isnan(num)) return "NaN"; From 40e9fa11abd2ab727bce5d5dbf5552e80f29fcb0 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 20:10:19 +0200 Subject: [PATCH 07/35] fix: bogus warning --- test/api/layout/partition_user_bins_test.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/test/api/layout/partition_user_bins_test.cpp b/test/api/layout/partition_user_bins_test.cpp index 1db95fd0..c56df2e8 100644 --- a/test/api/layout/partition_user_bins_test.cpp +++ b/test/api/layout/partition_user_bins_test.cpp @@ -43,7 +43,14 @@ std::vector> run_partition_user_bins(std::vector con { if (content_ids.empty()) { +#if CHOPPER_WORKAROUND_GCC_BOGUS_MEMMOV +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Warray-bounds=" +#endif // CHOPPER_WORKAROUND_GCC_BOGUS_MEMMOV content_ids.resize(kmer_counts.size()); +#if CHOPPER_WORKAROUND_GCC_BOGUS_MEMMOV +# pragma GCC diagnostic pop +#endif // CHOPPER_WORKAROUND_GCC_BOGUS_MEMMOV std::iota(content_ids.begin(), content_ids.end(), 0u); } From bc9fdee52f508d425e98e21eb0c33ba7417624c8 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 20:12:26 +0200 Subject: [PATCH 08/35] refactor: const and ranges --- .../chopper/layout/fast_layout_cluster.hpp | 20 +++++--- src/layout/fast_layout.cpp | 11 ++-- src/layout/partition_user_bins.cpp | 50 +++++++------------ 3 files changed, 38 insertions(+), 43 deletions(-) diff --git a/include/chopper/layout/fast_layout_cluster.hpp b/include/chopper/layout/fast_layout_cluster.hpp index e08566b4..808e46c0 100644 --- a/include/chopper/layout/fast_layout_cluster.hpp +++ b/include/chopper/layout/fast_layout_cluster.hpp @@ -73,12 +73,12 @@ struct Cluster return last; } - void add_user_bin(size_t user_bin) + void add_user_bin(size_t const user_bin) { user_bins.push_back(user_bin); } - bool is_valid(size_t id) const + bool is_valid(size_t const id) const { bool const ids_equal = representative_id == id; bool const properly_moved = has_been_moved() && empty(); @@ -96,17 +96,25 @@ struct Cluster void move_to(Cluster & target_cluster) { - target_cluster.user_bins.insert(target_cluster.user_bins.end(), this->user_bins.begin(), this->user_bins.end()); - this->user_bins.clear(); + auto & target = target_cluster.user_bins; + auto & source = this->user_bins; +#if __cpp_lib_containers_ranges + target.append_range(source); +#else + target.insert(target.end(), source.cbegin(), source.cend()); +#endif + source = std::vector{}; // .clear() AND release memory + moved_id = target_cluster.id(); } void sort_by_cardinality(std::vector const & cardinalities) { std::ranges::sort(user_bins, - [&cardinalities](auto const & v1, auto const & v2) + std::ranges::greater{}, + [&cardinalities](size_t const i) { - return cardinalities[v1] > cardinalities[v2]; + return cardinalities[i]; }); } }; diff --git a/src/layout/fast_layout.cpp b/src/layout/fast_layout.cpp index 0d661433..89a06f7b 100644 --- a/src/layout/fast_layout.cpp +++ b/src/layout/fast_layout.cpp @@ -8,7 +8,9 @@ #include #include #include +#include #include +#include #include #include @@ -435,12 +437,11 @@ void fast_layout(chopper::configuration const & config, // sort records ascending by the number of bin indices (corresponds to the IBF levels) // GCOVR_EXCL_START std::ranges::sort(hibf_layout.max_bins, - [](auto const & r, auto const & l) + std::ranges::less{}, + [](auto const & mb) { - if (r.previous_TB_indices.size() == l.previous_TB_indices.size()) - return std::ranges::lexicographical_compare(r.previous_TB_indices, l.previous_TB_indices); - else - return r.previous_TB_indices.size() < l.previous_TB_indices.size(); + // std::cref: compare the vector by reference instead of copying it + return std::make_tuple(mb.previous_TB_indices.size(), std::cref(mb.previous_TB_indices)); }); // GCOVR_EXCL_STOP diff --git a/src/layout/partition_user_bins.cpp b/src/layout/partition_user_bins.cpp index 16bd7b19..87f8d810 100644 --- a/src/layout/partition_user_bins.cpp +++ b/src/layout/partition_user_bins.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -90,8 +91,8 @@ auto LSH_fill_hashtable(std::vector const & clusters, for (auto & [key, list] : table) { std::ranges::sort(list); - auto const ret = std::ranges::unique(list); - list.erase(ret.begin(), ret.end()); + auto const [first, last] = std::ranges::unique(list); + list.erase(first, last); } return table; @@ -231,7 +232,6 @@ std::vector very_similar_LSH_clustering(std::vector & clusters, clusters[pos].sort_by_cardinality(cardinalities); } + // The user bins are sorted by cardinality, so the first one is the largest. + // Non-empty clusters have cardinality > 0, hence empty clusters are sorted last. + auto const largest_user_bin_cardinality = [&cardinalities](Cluster const & c) + { + return c.empty() ? size_t{} : cardinalities[c.contained_user_bins().front()]; + }; + // push largest p clusters to the front std::ranges::partial_sort(clusters, std::ranges::next(clusters.begin(), config.hibf_config.tmax, clusters.end()), - [&cardinalities](auto const & v1, auto const & v2) + std::ranges::greater{}, + [&largest_user_bin_cardinality](Cluster const & c) { - // Note: If v2 is empty, so is v1. - if (v1.size() == v2.size() && !v2.empty()) - return cardinalities[v1.contained_user_bins().front()] - > cardinalities[v2.contained_user_bins().front()]; - - return v1.size() > v2.size(); + return std::tuple{c.size(), largest_user_bin_cardinality(c)}; }); // after filling up the partitions with the biggest clusters, sort the clusters by cardinality of the biggest ub // s.t. that euqally sizes ub are assigned after each other and the small stuff is added at last. // the largest ub is already at the start because of former sorting. - std::ranges::sort(std::ranges::next(clusters.begin(), config.hibf_config.tmax, clusters.end()), - clusters.end(), - [&cardinalities](auto const & v1, auto const & v2) - { - if (v1.empty()) - return false; // v1 can never be larger than v2 then - - if (v2.empty()) // and v1 is not, since the first if would catch - return true; - - return cardinalities[v1.contained_user_bins().front()] - > cardinalities[v2.contained_user_bins().front()]; - }); + std::ranges::sort(clusters | std::views::drop(config.hibf_config.tmax), + std::ranges::greater{}, + largest_user_bin_cardinality); assert(clusters.size() < 2 || clusters[0].size() >= clusters[1].size()); // sanity check - -#ifndef NDEBUG - for (size_t cidx = 1; cidx < clusters.size(); ++cidx) - { - // once empty - always empty; all empty clusters should be at the end - if (clusters[cidx - 1].empty() && !clusters[cidx].empty()) - throw std::runtime_error{"sorting did not work"}; - } -#endif + // once empty - always empty; all empty clusters should be at the end + assert(std::ranges::is_partitioned(clusters, std::not_fn(&Cluster::empty))); } /*!\brief Assigns a whole cluster of user bins to the partition where adding it costs the least. From f9e95a50ea3ea9bdbc8113850194bf20bc684df3 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 20:44:55 +0200 Subject: [PATCH 09/35] fix: bogus only on gcc 16 --- include/chopper/workarounds.hpp | 10 ++++++++++ test/api/layout/partition_user_bins_test.cpp | 8 ++++---- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/include/chopper/workarounds.hpp b/include/chopper/workarounds.hpp index c8d98a13..bbc2b39a 100644 --- a/include/chopper/workarounds.hpp +++ b/include/chopper/workarounds.hpp @@ -33,3 +33,13 @@ # define CHOPPER_WORKAROUND_GCC_BOGUS_MEMMOV 0 # endif #endif + +/*!\brief Workaround bogus memmov errors in GCC 16. (Warray-bounds) + */ +#ifndef CHOPPER_WORKAROUND_GCC_BOGUS_ARRAY +# if CHOPPER_COMPILER_IS_GCC && (__GNUC__ == 16) +# define CHOPPER_WORKAROUND_GCC_BOGUS_ARRAY 1 +# else +# define CHOPPER_WORKAROUND_GCC_BOGUS_ARRAY 0 +# endif +#endif diff --git a/test/api/layout/partition_user_bins_test.cpp b/test/api/layout/partition_user_bins_test.cpp index c56df2e8..e5bb7bc5 100644 --- a/test/api/layout/partition_user_bins_test.cpp +++ b/test/api/layout/partition_user_bins_test.cpp @@ -43,14 +43,14 @@ std::vector> run_partition_user_bins(std::vector con { if (content_ids.empty()) { -#if CHOPPER_WORKAROUND_GCC_BOGUS_MEMMOV +#if CHOPPER_WORKAROUND_GCC_BOGUS_ARRAY # pragma GCC diagnostic push # pragma GCC diagnostic ignored "-Warray-bounds=" -#endif // CHOPPER_WORKAROUND_GCC_BOGUS_MEMMOV +#endif // CHOPPER_WORKAROUND_GCC_BOGUS_ARRAY content_ids.resize(kmer_counts.size()); -#if CHOPPER_WORKAROUND_GCC_BOGUS_MEMMOV +#if CHOPPER_WORKAROUND_GCC_BOGUS_ARRAY # pragma GCC diagnostic pop -#endif // CHOPPER_WORKAROUND_GCC_BOGUS_MEMMOV +#endif // CHOPPER_WORKAROUND_GCC_BOGUS_ARRAY std::iota(content_ids.begin(), content_ids.end(), 0u); } From d7e186987ba4414e15ebe97e7f69e511c6bbd78e Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 20:56:28 +0200 Subject: [PATCH 10/35] fix: do not drop user bins of large seed clusters A seed cluster with more than tmax user bins whose cardinality is at most 0.05 * sum_of_cardinalities / tmax places only tmax user bins. The rest was only queued for assignment if split_cluster was set, so these user bins ended up in no partition. Queue the unplaced user bins whenever there are any. This also stops queueing empty remainders. Co-Authored-By: Claude Opus 5.5 --- src/layout/partition_user_bins.cpp | 8 ++--- test/api/layout/partition_user_bins_test.cpp | 34 ++++++++++++++++---- 2 files changed, 31 insertions(+), 11 deletions(-) diff --git a/src/layout/partition_user_bins.cpp b/src/layout/partition_user_bins.cpp index 87f8d810..f9802460 100644 --- a/src/layout/partition_user_bins.cpp +++ b/src/layout/partition_user_bins.cpp @@ -563,11 +563,9 @@ size_t lsh_sim_approach(chopper::configuration const & config, min_partition_cardinality[p] = std::min(min_partition_cardinality[p], cardinalities[user_bin_idx]); } - if (split_cluster) - { - std::vector remainder(cluster.begin() + end, cluster.end()); - remaining_clusters.insert(remaining_clusters.end(), remainder); - } + // User bins that were not placed, either because at most tmax are placed or because the partitions ran out. + if (end < cluster.size()) + remaining_clusters.emplace_back(cluster.begin() + end, cluster.end()); ++cidx; } diff --git a/test/api/layout/partition_user_bins_test.cpp b/test/api/layout/partition_user_bins_test.cpp index e5bb7bc5..da32bc1f 100644 --- a/test/api/layout/partition_user_bins_test.cpp +++ b/test/api/layout/partition_user_bins_test.cpp @@ -7,12 +7,13 @@ #include // for Test, TestInfo, EXPECT_EQ, TEST -#include // for size_t -#include // for uint64_t -#include // for exit -#include // for iota -#include // for to_string -#include // for vector +#include // for sort, unique +#include // for size_t +#include // for uint64_t +#include // for exit +#include // for iota +#include // for to_string +#include // for vector #include #include @@ -188,3 +189,24 @@ TEST(partition_user_bins_test, several_clusters) ASSERT_EQ(partitions.size(), 64); } + +TEST(partition_user_bins_test, cluster_larger_than_tmax_with_small_cardinality) +{ + // User bins 2 to 6 are identical and form one cluster of 5 > tmax user bins. Its cardinality is below + // 0.05 * sum_of_cardinalities / tmax, so only tmax user bins seed a partition. The others must still be assigned. + // User bins 7 to 9 provide enough clusters, so the large cluster is not broken up beforehand. + std::vector const kmer_counts{300'000, 300'000, 1'000, 1'000, 1'000, 1'000, 1'000, 1'000, 1'000, 1'000}; + std::vector const content_ids{0, 1, 2, 2, 2, 2, 2, 3, 4, 5}; + auto const partitions = run_partition_user_bins(kmer_counts, /*tmax*/ 4, content_ids); + + std::vector assigned_user_bins{}; + for (auto const & partition : partitions) + assigned_user_bins.insert(assigned_user_bins.end(), partition.begin(), partition.end()); + std::ranges::sort(assigned_user_bins); + auto const [first, last] = std::ranges::unique(assigned_user_bins); + assigned_user_bins.erase(first, last); + + std::vector expected_user_bins(kmer_counts.size()); + std::iota(expected_user_bins.begin(), expected_user_bins.end(), 0u); + EXPECT_EQ(assigned_user_bins, expected_user_bins); +} From 2a63a7e5fed83fe587afe0732218f3599e8260e2 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 20:57:57 +0200 Subject: [PATCH 11/35] fix: compute minhash sketches only for the fast layout hibf's compute_sketches throws if a user bin has too few k-mers to fill all MinHash sketches (about 640). Since chopper layout always computed them, the default layout failed for such user bins, which worked before. Compute MinHash sketches only with --fast-layout. Sketch files written without --fast-layout contain no MinHash sketches again. Co-Authored-By: Claude Opus 5.5 --- src/chopper_layout.cpp | 7 ++++- test/cli/cli_chopper_basic_test.cpp | 44 +++++++++++++++++++++++++++++ test/cli/cli_output_sketches.cpp | 43 +++++++++++++++++++++++++++- 3 files changed, 92 insertions(+), 2 deletions(-) diff --git a/src/chopper_layout.cpp b/src/chopper_layout.cpp index 9df893a3..580ad61e 100644 --- a/src/chopper_layout.cpp +++ b/src/chopper_layout.cpp @@ -130,7 +130,12 @@ int chopper_layout(chopper::configuration & config, sharg::parser & parser) if (!input_is_a_sketch_file) { config.compute_sketches_timer.start(); - seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); + // Only the fast layout needs MinHash sketches. Computing them requires enough k-mers per user bin and throws + // otherwise, so the default layout must not compute them. + if (config.fast_layout) + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); + else + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches); config.compute_sketches_timer.stop(); } diff --git a/test/cli/cli_chopper_basic_test.cpp b/test/cli/cli_chopper_basic_test.cpp index 80d31251..6efb7a89 100644 --- a/test/cli/cli_chopper_basic_test.cpp +++ b/test/cli/cli_chopper_basic_test.cpp @@ -94,3 +94,47 @@ TEST_F(cli_test, chopper_cmd_kmer_bigger_than_window) EXPECT_EQ(result.out, std::string{}); EXPECT_EQ(result.err, std::string{"[ERROR] The k-mer size cannot be bigger than the window size.\n"}); } + +TEST_F(cli_test, chopper_user_bin_with_few_kmers) +{ + seqan3::test::tmp_directory tmp_dir{}; + std::filesystem::path const input_filename{tmp_dir.path() / "data.filenames"}; + std::filesystem::path const short_filename{tmp_dir.path() / "short.fa"}; + std::filesystem::path const layout_filename{tmp_dir.path() / "output.binning"}; + + // 200 bases give 182 k-mers, too few to fill the MinHash sketches that only the fast layout needs. + { + std::ofstream fout{short_filename}; + fout << ">short\n"; + for (size_t i = 0; i < 200; ++i) + fout << "ACGT"[(i * 7 + i / 3) % 4]; + fout << '\n'; + } + + { + std::ofstream fout{input_filename}; + fout << data("seq1.fa").string() << '\n' << short_filename.string() << '\n'; + } + + cli_test_result result = + execute_app("chopper", "--input", input_filename.c_str(), "--tmax", "64", "--output", layout_filename.c_str()); + + EXPECT_EQ(result.exit_code, 0); + EXPECT_EQ(result.out, std::string{}); + EXPECT_EQ(result.err, std::string{}); + EXPECT_TRUE(std::filesystem::exists(layout_filename)); + + // The fast layout needs MinHash sketches and fails. + result = execute_app("chopper", + "--fast-layout", + "--input", + input_filename.c_str(), + "--tmax", + "64", + "--output", + layout_filename.c_str()); + + EXPECT_NE(result.exit_code, 0); + EXPECT_EQ(result.out, std::string{}); + EXPECT_TRUE(result.err.starts_with("[ERROR] Not enough kmers")) << result.err; +} diff --git a/test/cli/cli_output_sketches.cpp b/test/cli/cli_output_sketches.cpp index 2d47e019..571408db 100644 --- a/test/cli/cli_output_sketches.cpp +++ b/test/cli/cli_output_sketches.cpp @@ -96,9 +96,50 @@ TEST_F(cli_test, chopper_layout) EXPECT_EQ(sin.filenames.size(), 3); EXPECT_EQ(sin.hll_sketches.size(), 3); - EXPECT_EQ(sin.minHash_sketches.size(), 3); // currently, no minhash sketches are needed in chopper layout + EXPECT_EQ(sin.minHash_sketches.size(), 0); // MinHash sketches are only computed for --fast-layout EXPECT_EQ(sin.filenames[0][0], data("seq1.fa").string()); EXPECT_EQ(sin.filenames[1][0], data("seq2.fa").string()); EXPECT_EQ(sin.filenames[2][0], data("seq3.fa").string()); } + +TEST_F(cli_test, chopper_layout_fast_layout) +{ + seqan3::test::tmp_directory tmp_dir{}; + std::filesystem::path const input_filename{tmp_dir.path() / "data.filenames"}; + std::filesystem::path const layout_filename{tmp_dir.path() / "output.binning"}; + std::filesystem::path const sketches_filename{tmp_dir.path() / "out.sketches"}; + + { + std::ofstream fout{input_filename}; + fout << data("seq1.fa").string() << '\n' + << data("seq2.fa").string() << '\n' + << data("seq3.fa").string() << '\n'; + } + + cli_test_result result = execute_app("chopper", + "--fast-layout", + "--input", + input_filename.c_str(), + "--tmax", + "64", + "--output-sketches-to", + sketches_filename.c_str(), + "--output", + layout_filename.c_str()); + + EXPECT_EQ(result.exit_code, 0); + EXPECT_EQ(result.out, std::string{}); + EXPECT_EQ(result.err, std::string{}); + + ASSERT_TRUE(std::filesystem::exists(sketches_filename)); + + chopper::sketch::sketch_file sin{}; + + std::ifstream is{sketches_filename}; + cereal::BinaryInputArchive iarchive{is}; + iarchive(sin); + + EXPECT_EQ(sin.hll_sketches.size(), 3); + EXPECT_EQ(sin.minHash_sketches.size(), 3); // the fast layout needs MinHash sketches +} From 48cb208485e24731d5e23139f32bf06bbe8753c6 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 20:59:09 +0200 Subject: [PATCH 12/35] fix: reject --fast-layout for sketch files without minhash sketches Sketch files written without --fast-layout, including all sketch files written before the fast layout existed, contain no MinHash sketches. The fast layout then failed an assertion in Debug and read out of bounds in Release. Co-Authored-By: Claude Opus 5.5 --- src/chopper_layout.cpp | 4 ++++ test/cli/cli_chopper_layout_from_sketch_file.cpp | 15 +++++++++++++++ 2 files changed, 19 insertions(+) diff --git a/src/chopper_layout.cpp b/src/chopper_layout.cpp index 580ad61e..a46905ed 100644 --- a/src/chopper_layout.cpp +++ b/src/chopper_layout.cpp @@ -109,6 +109,10 @@ int chopper_layout(chopper::configuration & config, sharg::parser & parser) sketches = std::move(sin.hll_sketches); minHash_sketches = std::move(sin.minHash_sketches); validate_configuration(parser, config, sin.chopper_config); + + if (config.fast_layout && minHash_sketches.size() != sketches.size()) + throw sharg::parser_error{"The sketch file does not contain MinHash sketches, which --fast-layout needs. " + "Create the sketch file with --fast-layout."}; } else { diff --git a/test/cli/cli_chopper_layout_from_sketch_file.cpp b/test/cli/cli_chopper_layout_from_sketch_file.cpp index 6e2e6288..4d80a1db 100644 --- a/test/cli/cli_chopper_layout_from_sketch_file.cpp +++ b/test/cli/cli_chopper_layout_from_sketch_file.cpp @@ -166,4 +166,19 @@ TEST_F(cli_test, chopper_layout_from_sketch_file) EXPECT_NE(result3.exit_code, 0); EXPECT_EQ(result3.out, std::string{}); EXPECT_EQ(result3.err, std::string{"[ERROR] You cannot set --sketch-bits when using a sketch file as input.\n"}); + + // The sketch file contains no MinHash sketches. + cli_test_result result4 = execute_app("chopper", + "--threads 2", + "--fast-layout", + "--input", + input_filename.c_str(), + "--tmax 64", + "--output", + binning_filename.c_str()); + EXPECT_NE(result4.exit_code, 0); + EXPECT_EQ(result4.out, std::string{}); + EXPECT_EQ(result4.err, + std::string{"[ERROR] The sketch file does not contain MinHash sketches, which --fast-layout needs. " + "Create the sketch file with --fast-layout.\n"}); } From 15274984b5df96a9f2fea2346c3578bb6582c76a Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 20:59:54 +0200 Subject: [PATCH 13/35] fix: reject --fast-layout together with --determine-best-tmax determine_best_number_of_technical_bins always uses the DP layout, so --fast-layout was silently ignored. Reject the combination right after parsing, before any sketching. Co-Authored-By: Claude Opus 5.5 --- src/chopper_layout.cpp | 3 +++ src/layout/execute.cpp | 2 +- test/cli/cli_chopper_basic_test.cpp | 15 +++++++++++++++ 3 files changed, 19 insertions(+), 1 deletion(-) diff --git a/src/chopper_layout.cpp b/src/chopper_layout.cpp index a46905ed..c193211b 100644 --- a/src/chopper_layout.cpp +++ b/src/chopper_layout.cpp @@ -78,6 +78,9 @@ int chopper_layout(chopper::configuration & config, sharg::parser & parser) else if (config.k > config.window_size) throw sharg::parser_error{"The k-mer size cannot be bigger than the window size."}; + if (config.fast_layout && config.determine_best_tmax) + throw sharg::parser_error{"You cannot combine --fast-layout with --determine-best-tmax."}; + auto has_sketch_file_extension = [](std::filesystem::path const & path) { return path.string().ends_with(".sketch") || path.string().ends_with(".sketches"); diff --git a/src/layout/execute.cpp b/src/layout/execute.cpp index 8fd9278f..b3d1c70b 100644 --- a/src/layout/execute.cpp +++ b/src/layout/execute.cpp @@ -43,7 +43,7 @@ int execute(chopper::configuration & config, if (config.determine_best_tmax) { - // ToDo, what about determine_best_tmax iwth fast layout? + // Always uses the DP layout; config.fast_layout is ignored. chopper_layout rejects the combination. hibf_layout = determine_best_number_of_technical_bins(config, cardinalities, sketches); } else diff --git a/test/cli/cli_chopper_basic_test.cpp b/test/cli/cli_chopper_basic_test.cpp index 6efb7a89..3f0093d7 100644 --- a/test/cli/cli_chopper_basic_test.cpp +++ b/test/cli/cli_chopper_basic_test.cpp @@ -138,3 +138,18 @@ TEST_F(cli_test, chopper_user_bin_with_few_kmers) EXPECT_EQ(result.out, std::string{}); EXPECT_TRUE(result.err.starts_with("[ERROR] Not enough kmers")) << result.err; } + +TEST_F(cli_test, chopper_fast_layout_with_determine_best_tmax) +{ + cli_test_result result = execute_app("chopper", + "--fast-layout", + "--determine-best-tmax", + "--input", + data("seq1.fa").c_str(), + "--output", + "output.binning"); + + EXPECT_NE(result.exit_code, 0); + EXPECT_EQ(result.out, std::string{}); + EXPECT_EQ(result.err, std::string{"[ERROR] You cannot combine --fast-layout with --determine-best-tmax.\n"}); +} From bb1aa72e24d568b8d3366541142038c3ed8ff550 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 21:03:53 +0200 Subject: [PATCH 14/35] fix: use the configured number of threads in the fast layout The parallel region had no num_threads clause, so it used all hardware threads (or OMP_NUM_THREADS) and ignored --threads. Co-Authored-By: Claude Opus 5.5 --- src/layout/fast_layout.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/layout/fast_layout.cpp b/src/layout/fast_layout.cpp index 89a06f7b..bd6bc936 100644 --- a/src/layout/fast_layout.cpp +++ b/src/layout/fast_layout.cpp @@ -400,7 +400,7 @@ void fast_layout(chopper::configuration const & config, hibf_layout.top_level_max_bin_id = max_bin_id; config.small_layouts_timer.start(); -#pragma omp parallel +#pragma omp parallel num_threads(config.hibf_config.threads) #pragma omp single { #pragma omp taskloop From ace02e1f82adba66a1ee3352fc29b9332e6883e1 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 21:05:59 +0200 Subject: [PATCH 15/35] fix: data race on the lsh and search timers lsh_sim_approach runs concurrently for different merged bins, but concurrent_timer::start() and stop() write shared, non-atomic time points; only operator+=() is thread-safe. Besides wrong timings, stop()'s assertion stop_point >= start_point failed in Debug when another thread restarted the timer in between (12 of 40 runs of the new fast layout recursion test with 8 threads). Time with local serial_timers and add them to the configuration's timers. Co-Authored-By: Claude Opus 5.5 --- src/layout/partition_user_bins.cpp | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/layout/partition_user_bins.cpp b/src/layout/partition_user_bins.cpp index f9802460..9c808270 100644 --- a/src/layout/partition_user_bins.cpp +++ b/src/layout/partition_user_bins.cpp @@ -29,6 +29,7 @@ #include #include #include +#include #include namespace chopper::layout @@ -474,8 +475,13 @@ size_t lsh_sim_approach(chopper::configuration const & config, std::vector max_partition_cardinality(number_of_remaining_tbs, 0u); std::vector min_partition_cardinality(number_of_remaining_tbs, std::numeric_limits::max()); + // lsh_sim_approach runs concurrently for different merged bins. concurrent_timer::start() and stop() are not + // thread-safe, only operator+=() is. Hence, time locally and add the result to the configuration's timers. + seqan::hibf::serial_timer lsh_algorithm_timer{}; + seqan::hibf::serial_timer search_partition_algorithm_timer{}; + // initial partitioning using locality sensitive hashing (LSH) - config.lsh_algorithm_timer.start(); + lsh_algorithm_timer.start(); std::vector clusters = very_similar_LSH_clustering(minHash_sketches, sorted_positions2, cardinalities, @@ -483,7 +489,8 @@ size_t lsh_sim_approach(chopper::configuration const & config, technical_bin_size_threshold, config); post_process_clusters(clusters, cardinalities, config); - config.lsh_algorithm_timer.stop(); + lsh_algorithm_timer.stop(); + config.lsh_algorithm_timer += lsh_algorithm_timer; // There must be more non-empty clusters than technical bins size_t number_of_remaining_clusters = std::ranges::count_if(clusters, @@ -584,7 +591,7 @@ size_t lsh_sim_approach(chopper::configuration const & config, { auto const & cluster = remaining_clusters[ridx]; - config.search_partition_algorithm_timer.start(); + search_partition_algorithm_timer.start(); find_best_partition(config, number_of_remaining_tbs, merged_threshold, @@ -595,8 +602,9 @@ size_t lsh_sim_approach(chopper::configuration const & config, partition_sketches, max_partition_cardinality, min_partition_cardinality); - config.search_partition_algorithm_timer.stop(); + search_partition_algorithm_timer.stop(); } + config.search_partition_algorithm_timer += search_partition_algorithm_timer; // compute actual max size size_t max_size{0}; From 6bcde8a58b30e8d9ea00824528fe31aec4d773a5 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 21:06:08 +0200 Subject: [PATCH 16/35] test: cover the recursive fast layout No test reached fast_layout_recursion and add_level_to_layout: a merged bin needs at least 64 * tmax user bins for that, and the largest fast layout test had 96 user bins with tmax = 64. With tmax = 4 and 2000 user bins, each of the 4 top-level merged bins (about 500 user bins) is laid out recursively, and the DP layout handles the level below. The layout is checked structurally: every user bin is stored exactly once, merged and stored technical bins do not overlap, and there is one max bin entry per lower-level IBF. Co-Authored-By: Claude Opus 5.5 --- test/api/layout/CMakeLists.txt | 1 + test/api/layout/fast_layout_test.cpp | 146 +++++++++++++++++++++++++++ 2 files changed, 147 insertions(+) create mode 100644 test/api/layout/fast_layout_test.cpp diff --git a/test/api/layout/CMakeLists.txt b/test/api/layout/CMakeLists.txt index 92101d16..907851aa 100644 --- a/test/api/layout/CMakeLists.txt +++ b/test/api/layout/CMakeLists.txt @@ -29,4 +29,5 @@ add_api_test (input_test.cpp) add_api_test (determine_split_bins_test.cpp) add_api_test (fast_layout_cluster_test.cpp) add_api_test (fast_layout_find_bins_to_be_split_test.cpp) +add_api_test (fast_layout_test.cpp) add_api_test (partition_user_bins_test.cpp) diff --git a/test/api/layout/fast_layout_test.cpp b/test/api/layout/fast_layout_test.cpp new file mode 100644 index 00000000..1e91d041 --- /dev/null +++ b/test/api/layout/fast_layout_test.cpp @@ -0,0 +1,146 @@ +// -------------------------------------------------------------------------------------------------- +// Copyright (c) 2006-2023, Knut Reinert & Freie Universität Berlin +// Copyright (c) 2016-2023, Knut Reinert & MPI für molekulare Genetik +// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License +// shipped with this file and also available at: https://github.com/seqan/chopper/blob/main/LICENSE.md +// -------------------------------------------------------------------------------------------------- + +#include // for Test, TestInfo, EXPECT_EQ, TEST + +#include // for max, ranges::max +#include // for size_t +#include // for uint64_t +#include // for map +#include // for iota +#include // for set +#include // for vector + +#include +#include + +#include +#include +#include +#include +#include + +namespace +{ + +// splitmix64 finaliser, see partition_user_bins_test.cpp. +uint64_t scramble(uint64_t x) +{ + x += 0x9e3779b97f4a7c15ULL; + x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9ULL; + x = (x ^ (x >> 27)) * 0x94d049bb133111ebULL; + return x ^ (x >> 31); +} + +// Checks that `hibf_layout` describes a consistent HIBF over `number_of_user_bins` user bins: +// * user_bins[i].idx == i +// * In every IBF (identified by its path of merged bin indices), each technical bin is either a merged bin (it is on the +// path of some user bin) or holds exactly one user bin. +// * max_bins contains exactly one entry per lower-level IBF. +// * All technical bin indices are below next_multiple_of_64(tmax): The lowest levels of the DP layout use +// next_multiple_of_64(#user bins) technical bins, which is only bounded by tmax if tmax is a multiple of 64. +void check_layout(seqan::hibf::layout::layout const & hibf_layout, size_t const number_of_user_bins, size_t const tmax) +{ + ASSERT_EQ(hibf_layout.user_bins.size(), number_of_user_bins); + + size_t const max_technical_bins{seqan::hibf::next_multiple_of_64(tmax)}; + + constexpr size_t merged_bin{static_cast(-1)}; + constexpr size_t empty_bin{static_cast(-2)}; + std::map, std::vector> ibfs{}; // path -> content of each technical bin + + auto get_ibf = [&](std::vector const & path) -> std::vector & + { + return ibfs.try_emplace(path, max_technical_bins, empty_bin).first->second; + }; + + for (size_t ub = 0; ub < number_of_user_bins; ++ub) + { + auto const & user_bin = hibf_layout.user_bins[ub]; + ASSERT_EQ(user_bin.idx, ub); + + std::vector path{}; + for (size_t const tb : user_bin.previous_TB_indices) + { + ASSERT_LT(tb, max_technical_bins); + auto & bin = get_ibf(path)[tb]; + ASSERT_TRUE(bin == empty_bin || bin == merged_bin) << "user bin " << ub << " passes through a filled bin"; + bin = merged_bin; + path.push_back(tb); + } + + ASSERT_GE(user_bin.number_of_technical_bins, 1u); + ASSERT_LE(user_bin.storage_TB_id + user_bin.number_of_technical_bins, max_technical_bins); + auto & ibf = get_ibf(path); + for (size_t tb = user_bin.storage_TB_id; tb < user_bin.storage_TB_id + user_bin.number_of_technical_bins; ++tb) + { + ASSERT_EQ(ibf[tb], empty_bin) << "user bin " << ub << " is stored in an occupied bin"; + ibf[tb] = ub; + } + } + + EXPECT_LT(hibf_layout.top_level_max_bin_id, tmax); // the top level is laid out by the fast layout + + std::set> lower_level_ibfs{}; + for (auto const & [path, bins] : ibfs) + if (!path.empty()) + lower_level_ibfs.insert(path); + + std::set> max_bin_ibfs{}; + for (auto const & max_bin : hibf_layout.max_bins) + { + EXPECT_LT(max_bin.id, max_technical_bins); + EXPECT_TRUE(max_bin_ibfs.insert(max_bin.previous_TB_indices).second) << "duplicate max bin entry"; + } + + EXPECT_EQ(max_bin_ibfs, lower_level_ibfs); +} + +} // namespace + +TEST(fast_layout_test, recursion) +{ + // With tmax = 4, a merged bin is laid out recursively with the fast layout if it contains at least 64 * 4 = 256 + // user bins of similar size. 2000 user bins in 4 top-level technical bins give about 500 per merged bin. + size_t const tmax{4}; + size_t const number_of_user_bins{2000}; + + chopper::configuration config{}; + config.fast_layout = true; + config.hibf_config.tmax = tmax; + config.hibf_config.threads = 2; + config.hibf_config.number_of_user_bins = number_of_user_bins; + config.hibf_config.disable_estimate_union = true; // also disables rearrangement + config.hibf_config.input_fn = [](size_t const ub, seqan::hibf::insert_iterator it) + { + uint64_t const offset = static_cast(ub) << 32; + for (uint64_t i = 0; i < 3'000 + ub % 7 * 100; ++i) + it = scramble(offset + i); + }; + + std::vector sketches{}; + std::vector minHash_sketches{}; + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); + + std::vector cardinalities(number_of_user_bins); + for (size_t i = 0; i < number_of_user_bins; ++i) + cardinalities[i] = sketches[i].estimate(); + + std::vector positions(number_of_user_bins); + std::iota(positions.begin(), positions.end(), 0u); + + seqan::hibf::layout::layout hibf_layout{}; + chopper::layout::fast_layout(config, positions, cardinalities, sketches, minHash_sketches, hibf_layout); + + check_layout(hibf_layout, number_of_user_bins, tmax); + + // Recursive layouts produce at least three levels. + size_t max_depth{0}; + for (auto const & user_bin : hibf_layout.user_bins) + max_depth = std::max(max_depth, user_bin.previous_TB_indices.size()); + EXPECT_GE(max_depth, 2u); +} From 3442ee0a44f8c3ac122532feab7a16debdd5e6a2 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 21:06:55 +0200 Subject: [PATCH 17/35] fix: underflow of the partition growth penalty find_best_partition computed union_estimate - current_partition_size and asserted that the union estimate is not smaller. HyperLogLog estimates are not monotonic under merging: at the switch from linear counting to the raw estimate, a superset can be estimated smaller (4 of 800000 steps in a probe with 12 sketch bits, by up to 51). Debug builds then aborted and Release builds wrapped to a huge penalty that excluded the partition. Clamp the penalty to 0. Co-Authored-By: Claude Opus 5.5 --- src/layout/partition_user_bins.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/layout/partition_user_bins.cpp b/src/layout/partition_user_bins.cpp index 9c808270..d6b63f13 100644 --- a/src/layout/partition_user_bins.cpp +++ b/src/layout/partition_user_bins.cpp @@ -399,8 +399,9 @@ bool find_best_partition(chopper::configuration const & config, size_t const union_estimate = union_sketch.estimate(); size_t const current_partition_size = partition_sketches[p].estimate(); - assert(union_estimate >= current_partition_size); - size_t const penalty_current_bin = union_estimate - current_partition_size; + // HyperLogLog estimates are not monotonic under merging: Where the estimate switches from linear counting to + // the raw estimate, the union can be estimated smaller than the partition. Clamp to 0 instead of underflowing. + size_t const penalty_current_bin = union_estimate - std::min(union_estimate, current_partition_size); size_t const penalty_current_ibf = config.hibf_config.tmax * ((union_estimate <= corrected_estimate_per_part) ? 0u : union_estimate - corrected_estimate_per_part); From 8527d87f580072c8c63f6c012d8f8aeab3811960 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 21:07:09 +0200 Subject: [PATCH 18/35] docs: describe the threshold update of find_bins_to_be_split correctly The documentation claimed that the threshold is multiplied by max(1.01, sum / (threshold * max_bins)) and that the number of split user bins is clamped. The code sets the threshold to max(threshold + 1, sum / max_bins), and the number of split user bins cannot exceed the number of technical bins, which is only asserted. Co-Authored-By: Claude Opus 5.5 --- .../chopper/layout/fast_layout_find_bins_to_be_split.hpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/include/chopper/layout/fast_layout_find_bins_to_be_split.hpp b/include/chopper/layout/fast_layout_find_bins_to_be_split.hpp index 6891a544..04ba0d22 100644 --- a/include/chopper/layout/fast_layout_find_bins_to_be_split.hpp +++ b/include/chopper/layout/fast_layout_find_bins_to_be_split.hpp @@ -34,11 +34,12 @@ namespace chopper::layout * * The user bins with a cardinality greater than `threshold` form a prefix of `sorted_positions`. The number of * technical bins they need is `sum / threshold` (rounded down), where `sum` is the sum of their cardinalities. - * If this exceeds `max_bins`, `threshold` is multiplied by `max(1.01, sum / (threshold * max_bins))` and the prefix - * is computed again, until the split user bins fit into `max_bins` technical bins. + * If this exceeds `max_bins`, `threshold` is set to `max(threshold + 1, sum / max_bins)` (rounded down) and the + * prefix is computed again, until the split user bins fit into `max_bins` technical bins. * - * The number of split user bins is clamped to the number of technical bins, so each split user bin gets at least - * one technical bin. + * The number of split user bins never exceeds the number of technical bins, so each split user bin gets at least one + * technical bin: Each split user bin has a cardinality greater than `threshold`, hence `sum / threshold` is at least + * the number of split user bins. This is asserted, not enforced. */ inline std::pair find_bins_to_be_split(std::vector const & sorted_positions, std::vector const & cardinalities, From 65d2727a3af621e35be56ee360f0b387bb605453 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 21:09:05 +0200 Subject: [PATCH 19/35] fix: sort empty clusters last regardless of cardinality post_process_clusters sorted the clusters after the first tmax by the cardinality of their largest user bin, with 0 for empty clusters. This put empty clusters last only if every non-empty cluster has a cardinality estimate above 0. Otherwise, a non-empty cluster could end up among the empty ones. Debug builds then failed the is_partitioned assertion, and in Release builds lsh_sim_approach stopped at the first empty cluster, so the user bins of such clusters were not assigned. Sort by (non-empty, cardinality) instead. For all other clusters, the order is unchanged. Co-Authored-By: Claude Opus 5.5 --- src/layout/partition_user_bins.cpp | 8 +++- test/api/layout/partition_user_bins_test.cpp | 43 +++++++++++++++----- 2 files changed, 39 insertions(+), 12 deletions(-) diff --git a/src/layout/partition_user_bins.cpp b/src/layout/partition_user_bins.cpp index d6b63f13..bb87804b 100644 --- a/src/layout/partition_user_bins.cpp +++ b/src/layout/partition_user_bins.cpp @@ -259,7 +259,6 @@ void post_process_clusters(std::vector & clusters, } // The user bins are sorted by cardinality, so the first one is the largest. - // Non-empty clusters have cardinality > 0, hence empty clusters are sorted last. auto const largest_user_bin_cardinality = [&cardinalities](Cluster const & c) { return c.empty() ? size_t{} : cardinalities[c.contained_user_bins().front()]; @@ -277,9 +276,14 @@ void post_process_clusters(std::vector & clusters, // after filling up the partitions with the biggest clusters, sort the clusters by cardinality of the biggest ub // s.t. that euqally sizes ub are assigned after each other and the small stuff is added at last. // the largest ub is already at the start because of former sorting. + // Empty clusters are sorted last explicitly. Their cardinality key 0 does not suffice, because a non-empty cluster + // can have a cardinality estimate of 0, too. std::ranges::sort(clusters | std::views::drop(config.hibf_config.tmax), std::ranges::greater{}, - largest_user_bin_cardinality); + [&largest_user_bin_cardinality](Cluster const & c) + { + return std::tuple{!c.empty(), largest_user_bin_cardinality(c)}; + }); assert(clusters.size() < 2 || clusters[0].size() >= clusters[1].size()); // sanity check // once empty - always empty; all empty clusters should be at the end diff --git a/test/api/layout/partition_user_bins_test.cpp b/test/api/layout/partition_user_bins_test.cpp index da32bc1f..25301551 100644 --- a/test/api/layout/partition_user_bins_test.cpp +++ b/test/api/layout/partition_user_bins_test.cpp @@ -38,9 +38,11 @@ uint64_t scramble(uint64_t x) // User bin `i` consists of `kmer_counts[i]` distinct hashes drawn from content `content_ids[i]`. // User bins with the same content id share their first min(kmer_counts) hashes, all others are disjoint. // If `content_ids` is empty, every user bin has its own content. +// The cardinalities of the user bins in `zero_cardinality_user_bins` are set to 0, regardless of their content. std::vector> run_partition_user_bins(std::vector const & kmer_counts, size_t const tmax, - std::vector content_ids = {}) + std::vector content_ids = {}, + std::vector const & zero_cardinality_user_bins = {}) { if (content_ids.empty()) { @@ -72,6 +74,8 @@ std::vector> run_partition_user_bins(std::vector con std::vector cardinalities(kmer_counts.size()); for (size_t i = 0; i < sketches.size(); ++i) cardinalities[i] = sketches[i].estimate(); + for (size_t const ub : zero_cardinality_user_bins) + cardinalities[ub] = 0u; std::vector positions(kmer_counts.size()); std::iota(positions.begin(), positions.end(), 0u); @@ -82,6 +86,22 @@ std::vector> run_partition_user_bins(std::vector con return partitions; } +// Expects that every user bin in [0, number_of_user_bins) is assigned to at least one partition. +void expect_all_user_bins_assigned(std::vector> const & partitions, + size_t const number_of_user_bins) +{ + std::vector assigned_user_bins{}; + for (auto const & partition : partitions) + assigned_user_bins.insert(assigned_user_bins.end(), partition.begin(), partition.end()); + std::ranges::sort(assigned_user_bins); + auto const [first, last] = std::ranges::unique(assigned_user_bins); + assigned_user_bins.erase(first, last); + + std::vector expected_user_bins(number_of_user_bins); + std::iota(expected_user_bins.begin(), expected_user_bins.end(), 0u); + EXPECT_EQ(assigned_user_bins, expected_user_bins); +} + } // namespace // A user bin is split if it is assigned to more than one technical bin. @@ -199,14 +219,17 @@ TEST(partition_user_bins_test, cluster_larger_than_tmax_with_small_cardinality) std::vector const content_ids{0, 1, 2, 2, 2, 2, 2, 3, 4, 5}; auto const partitions = run_partition_user_bins(kmer_counts, /*tmax*/ 4, content_ids); - std::vector assigned_user_bins{}; - for (auto const & partition : partitions) - assigned_user_bins.insert(assigned_user_bins.end(), partition.begin(), partition.end()); - std::ranges::sort(assigned_user_bins); - auto const [first, last] = std::ranges::unique(assigned_user_bins); - assigned_user_bins.erase(first, last); + expect_all_user_bins_assigned(partitions, kmer_counts.size()); +} - std::vector expected_user_bins(kmer_counts.size()); - std::iota(expected_user_bins.begin(), expected_user_bins.end(), 0u); - EXPECT_EQ(assigned_user_bins, expected_user_bins); +TEST(partition_user_bins_test, user_bin_with_zero_cardinality) +{ + // User bins 0 to 11 form four clusters of three identical user bins, which leaves eight empty (moved) clusters. + // User bin 14 has an estimated cardinality of 0. Its non-empty cluster must still be sorted before the empty + // clusters, which also have the key 0, or it is not assigned. + std::vector const kmer_counts(15, 3'000); + std::vector const content_ids{0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 5, 6}; + auto const partitions = run_partition_user_bins(kmer_counts, /*tmax*/ 4, content_ids, {14}); + + expect_all_user_bins_assigned(partitions, kmer_counts.size()); } From fa4a1a4dc4d5d5e466d921362734e61bcc22fa98 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 21:37:03 +0200 Subject: [PATCH 20/35] fix: include workarounds.hpp where its macros are used, warn on undefined macros partition_user_bins_test.cpp and check_filenames.cpp use the CHOPPER_WORKAROUND_GCC_BOGUS_* macros without including chopper/workarounds.hpp. The undefined macros evaluate to 0 in #if, so the diagnostic pragmas were never active. This is why the GCC 16 -Warray-bounds workaround in partition_user_bins_test.cpp had no effect. -Wundef turns such mistakes into errors. For it, the feature test in fast_layout_cluster.hpp uses #ifdef, because __cpp_lib_containers_ranges is not defined before C++23 or with older standard libraries. Co-Authored-By: Claude Opus 5.5 --- include/chopper/layout/fast_layout_cluster.hpp | 2 +- src/CMakeLists.txt | 2 +- src/sketch/check_filenames.cpp | 1 + test/api/layout/partition_user_bins_test.cpp | 1 + 4 files changed, 4 insertions(+), 2 deletions(-) diff --git a/include/chopper/layout/fast_layout_cluster.hpp b/include/chopper/layout/fast_layout_cluster.hpp index 808e46c0..b50b6bad 100644 --- a/include/chopper/layout/fast_layout_cluster.hpp +++ b/include/chopper/layout/fast_layout_cluster.hpp @@ -98,7 +98,7 @@ struct Cluster { auto & target = target_cluster.user_bins; auto & source = this->user_bins; -#if __cpp_lib_containers_ranges +#ifdef __cpp_lib_containers_ranges target.append_range(source); #else target.insert(target.end(), source.cbegin(), source.cend()); diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 75637058..7fdcf51f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -14,7 +14,7 @@ endif () add_library (chopper_interface INTERFACE) target_link_libraries (chopper_interface INTERFACE sharg::sharg seqan3::seqan3 seqan::hibf) target_include_directories (chopper_interface INTERFACE "${PROJECT_SOURCE_DIR}/include") -target_compile_options (chopper_interface INTERFACE "-pedantic" "-Wall" "-Wextra") +target_compile_options (chopper_interface INTERFACE "-pedantic" "-Wall" "-Wextra" "-Wundef") add_library (chopper::interface ALIAS chopper_interface) option (CHOPPER_WITH_WERROR "Report compiler warnings as errors." ON) diff --git a/src/sketch/check_filenames.cpp b/src/sketch/check_filenames.cpp index 2bf48d44..df7074e3 100644 --- a/src/sketch/check_filenames.cpp +++ b/src/sketch/check_filenames.cpp @@ -16,6 +16,7 @@ #include #include +#include namespace chopper::sketch { diff --git a/test/api/layout/partition_user_bins_test.cpp b/test/api/layout/partition_user_bins_test.cpp index 25301551..11a81b99 100644 --- a/test/api/layout/partition_user_bins_test.cpp +++ b/test/api/layout/partition_user_bins_test.cpp @@ -17,6 +17,7 @@ #include #include +#include #include #include From 717fc00d1608da36c7bc8e75aa1f3c2886946c9f Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 22:09:20 +0200 Subject: [PATCH 21/35] refactor: record the top level with add_level_to_layout fast_layout duplicated add_level_to_layout to record the top-level IBF. add_level_to_layout now sets the idx of the user bins it records and returns the maximum technical bin. fast_layout stores it as top_level_max_bin_id, fast_layout_recursion appends it to max_bins. The layouts are unchanged (compared on a recursive and a mixed input). Co-Authored-By: Claude Opus 5.5 --- src/layout/fast_layout.cpp | 111 +++++++------------------------------ 1 file changed, 21 insertions(+), 90 deletions(-) diff --git a/src/layout/fast_layout.cpp b/src/layout/fast_layout.cpp index bd6bc936..dce9810a 100644 --- a/src/layout/fast_layout.cpp +++ b/src/layout/fast_layout.cpp @@ -102,15 +102,17 @@ bool do_I_need_a_fast_layout(chopper::configuration const & config, return false; } -/*!\brief Records a lower-level IBF, computed by partition_user_bins, in the layout. +/*!\brief Records an IBF, computed by partition_user_bins, in the layout. * \param[in] config The configuration (uses the FPR settings of `hibf_config`). - * \param[in,out] hibf_layout The global layout. Its `user_bins` must be indexed by global user bin index - * (`user_bins[i].idx == i`), as initialised by fast_layout. + * \param[in,out] hibf_layout The global layout. Its `user_bins` must be indexed by global user bin index, i.e., + * `user_bins[i]` describes user bin `i`. The `idx` of each recorded user bin is set. * \param[in] partitions The technical bins of the new IBF: `partitions[t]` holds the global user bin indices * in technical bin `t`. * \param[in] sketches The HyperLogLog sketch of each user bin, indexed by global user bin index. - * \param[in] previous The merged-bin path identifying the new IBF. Every user bin in `partitions` must - * currently have exactly this path as `previous_TB_indices`. + * \param[in] previous The merged-bin path identifying the new IBF, empty for the top level. Every user bin in + * `partitions` must currently have exactly this path as `previous_TB_indices`. + * \returns The id of the technical bin with the largest FPR-corrected size (relaxed correction for merged bins, split + * correction for split bins). * * - **Merged bin** (more than one user bin): `t` is appended to each user bin's `previous_TB_indices`. Their * `storage_TB_id` is set later, when the next level down is recorded. @@ -119,16 +121,13 @@ bool do_I_need_a_fast_layout(chopper::configuration const & config, * contiguously. * - **Empty technical bins** are skipped. * - * The technical bin with the largest FPR-corrected size (relaxed correction for merged bins, split correction for - * split bins) is appended to `hibf_layout.max_bins` as `(previous, max_bin_id)`. - * - * Not thread-safe; callers serialise it with `omp critical`. + * Not thread-safe; fast_layout_recursion serialises it with `omp critical`. */ -void add_level_to_layout(chopper::configuration const & config, - seqan::hibf::layout::layout & hibf_layout, - std::vector> const & partitions, - std::vector const & sketches, - std::vector const & previous) +size_t add_level_to_layout(chopper::configuration const & config, + seqan::hibf::layout::layout & hibf_layout, + std::vector> const & partitions, + std::vector const & sketches, + [[maybe_unused]] std::vector const & previous) { size_t max_bin_id{0}; size_t max_size{0}; @@ -155,11 +154,11 @@ void add_level_to_layout(chopper::configuration const & config, for (size_t const user_bin_id : partition) { - assert(hibf_layout.user_bins[user_bin_id].idx == user_bin_id); auto & current_user_bin = hibf_layout.user_bins[user_bin_id]; // update assert(previous == current_user_bin.previous_TB_indices); + current_user_bin.idx = user_bin_id; current_user_bin.previous_TB_indices.push_back(partition_idx); current_sketch.merge(sketches[user_bin_id]); } @@ -179,7 +178,8 @@ void add_level_to_layout(chopper::configuration const & config, else // single or split bin (partition.size() == 1) { auto & current_user_bin = hibf_layout.user_bins[partitions[partition_idx][0]]; - assert(current_user_bin.idx == partitions[partition_idx][0]); + assert(previous == current_user_bin.previous_TB_indices); + current_user_bin.idx = partitions[partition_idx][0]; current_user_bin.storage_TB_id = partition_idx; current_user_bin.number_of_technical_bins = 1; // initialise to 1 @@ -202,7 +202,7 @@ void add_level_to_layout(chopper::configuration const & config, } } - hibf_layout.max_bins.emplace_back(previous, max_bin_id); // add lower level meta information + return max_bin_id; } /*!\brief Grafts a layout computed for a merged bin (see general_layout) into the global layout. @@ -273,7 +273,8 @@ void fast_layout_recursion(chopper::configuration const & config, #pragma omp critical { - add_level_to_layout(config, hibf_layout, tmax_partitions, sketches, previous); + size_t const max_bin_id = add_level_to_layout(config, hibf_layout, tmax_partitions, sketches, previous); + hibf_layout.max_bins.emplace_back(previous, max_bin_id); // add lower level meta information } for (size_t partition_idx = 0; partition_idx < tmax_partitions.size(); ++partition_idx) @@ -325,79 +326,9 @@ void fast_layout(chopper::configuration const & config, partition_user_bins(config, positions, cardinalities, sketches, minHash_sketches, tmax_partitions); config.intital_partition_timer.stop(); - auto const split_fpr_correction = - seqan::hibf::layout::compute_fpr_correction({.fpr = config.hibf_config.maximum_fpr, // - .hash_count = config.hibf_config.number_of_hash_functions, - .t_max = config.hibf_config.tmax}); - - double const relaxed_fpr_correction = seqan::hibf::layout::compute_relaxed_fpr_correction( - {.fpr = config.hibf_config.maximum_fpr, // - .relaxed_fpr = config.hibf_config.relaxed_fpr, - .hash_count = config.hibf_config.number_of_hash_functions}); - - size_t max_bin_id{0}; - size_t max_size{0}; - hibf_layout.user_bins.resize(config.hibf_config.number_of_user_bins); - // initialise user bins in layout - for (size_t partition_idx = 0; partition_idx < tmax_partitions.size(); ++partition_idx) - { - if (tmax_partitions[partition_idx].size() > 1) // merged bin - { - seqan::hibf::sketch::hyperloglog current_sketch{sketches[0]}; // ensure same bit size - current_sketch.reset(); - - for (size_t const user_bin_id : tmax_partitions[partition_idx]) - { - hibf_layout.user_bins[user_bin_id] = {.previous_TB_indices = {partition_idx}, - .storage_TB_id = 0 /*not determiend yet*/, - .number_of_technical_bins = 1 /*not determiend yet*/, - .idx = user_bin_id}; - current_sketch.merge(sketches[user_bin_id]); - } - - // update max_bin_id, max_size - size_t const current_size = current_sketch.estimate() * relaxed_fpr_correction; - if (current_size > max_size) - { - max_bin_id = partition_idx; - max_size = current_size; - } - } - else if (tmax_partitions[partition_idx].size() == 0) // should not happen. Edge case? - { - continue; - } - else // single or split bin (tmax_partitions[partition_idx].size() == 1) - { - assert(tmax_partitions[partition_idx].size() == 1); - size_t const user_bin_id = tmax_partitions[partition_idx][0]; - hibf_layout.user_bins[user_bin_id] = {.previous_TB_indices = {}, - .storage_TB_id = partition_idx, - .number_of_technical_bins = 1 /*determiend below*/, - .idx = user_bin_id}; - - while (partition_idx + 1 < tmax_partitions.size() && tmax_partitions[partition_idx].size() == 1 - && tmax_partitions[partition_idx + 1].size() == 1 - && tmax_partitions[partition_idx][0] == tmax_partitions[partition_idx + 1][0]) - { - ++hibf_layout.user_bins[user_bin_id].number_of_technical_bins; - ++partition_idx; - } - - // update max_bin_id, max_size - size_t const current_size = - sketches[user_bin_id].estimate() - * split_fpr_correction[hibf_layout.user_bins[user_bin_id].number_of_technical_bins]; - if (current_size > max_size) - { - max_bin_id = hibf_layout.user_bins[user_bin_id].storage_TB_id; - max_size = current_size; - } - } - } - - hibf_layout.top_level_max_bin_id = max_bin_id; + hibf_layout.user_bins.resize(config.hibf_config.number_of_user_bins); + hibf_layout.top_level_max_bin_id = add_level_to_layout(config, hibf_layout, tmax_partitions, sketches, {}); config.small_layouts_timer.start(); #pragma omp parallel num_threads(config.hibf_config.threads) From 2c923407a3e1feb57a052bc52e3d93b7c1ad7132 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 22:10:53 +0200 Subject: [PATCH 22/35] refactor: find_best_partition returns void find_best_partition always returned true and threw only if no partition was selected. With at least one partition, partition 0 is always selected unless its cost is SIZE_MAX, and lsh_sim_approach always passes at least one partition. Assert that instead of tracking best_p_found. Co-Authored-By: Claude Opus 5.5 --- src/layout/partition_user_bins.cpp | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/src/layout/partition_user_bins.cpp b/src/layout/partition_user_bins.cpp index bb87804b..d3821fd9 100644 --- a/src/layout/partition_user_bins.cpp +++ b/src/layout/partition_user_bins.cpp @@ -16,7 +16,6 @@ #include #include #include -#include #include #include #include @@ -292,7 +291,8 @@ void post_process_clusters(std::vector & clusters, /*!\brief Assigns a whole cluster of user bins to the partition where adding it costs the least. * \param[in] config The configuration (uses `hibf_config.tmax` and `sketch_bits`). - * \param[in] number_of_partitions Only partitions `[0, number_of_partitions)` are considered. + * \param[in] number_of_partitions Only partitions `[0, number_of_partitions)` are considered. Must be + * greater than 0. * \param[in,out] corrected_estimate_per_part The current target cardinality per partition. Raised to the chosen * partition's new estimate if that is larger. * \param[in] cluster The global user bin indices to assign together. @@ -302,7 +302,6 @@ void post_process_clusters(std::vector & clusters, * \param[in,out] partition_sketches The union sketch per partition. Updated for the chosen partition. * \param[in,out] max_partition_cardinality The largest user bin cardinality per partition. Updated. * \param[in,out] min_partition_cardinality The smallest user bin cardinality per partition. Updated. - * \returns `true`. Throws std::runtime_error if no partition was selected, i.e., if `number_of_partitions == 0`. * * For each partition p, the cost of adding the cluster is the sum of * - `union - |p|`: the new k-mers the cluster adds to p (HyperLogLog estimates), @@ -317,7 +316,7 @@ void post_process_clusters(std::vector & clusters, * The partition with the smallest cost is chosen. A partition with zero cost always replaces the current best, so if * several have zero cost, the last one wins. */ -bool find_best_partition(chopper::configuration const & config, +void find_best_partition(chopper::configuration const & config, size_t const number_of_partitions, size_t & corrected_estimate_per_part, std::vector const & cluster, @@ -328,6 +327,8 @@ bool find_best_partition(chopper::configuration const & config, std::vector & max_partition_cardinality, std::vector & min_partition_cardinality) { + assert(number_of_partitions > 0u); + seqan::hibf::sketch::hyperloglog const current_sketch = [&sketches, &cluster, &config]() { seqan::hibf::sketch::hyperloglog result{config.hibf_config.sketch_bits}; @@ -353,7 +354,6 @@ bool find_best_partition(chopper::configuration const & config, // "which partition has the largest intersection with user bin b compared to its own (partition) size." size_t smallest_change{std::numeric_limits::max()}; size_t best_p{0}; - bool best_p_found{false}; auto penalty_lower_level = [&](size_t const additional_number_of_user_bins, size_t const p) -> size_t { @@ -416,13 +416,9 @@ bool find_best_partition(chopper::configuration const & config, { smallest_change = change; best_p = p; - best_p_found = true; } } - if (!best_p_found) - throw std::runtime_error{"currently there are no safety measures if a partition is not found"}; - // now that we know which partition fits best (`best_p`), add those indices to it for (size_t const user_bin_idx : cluster) { @@ -432,8 +428,6 @@ bool find_best_partition(chopper::configuration const & config, } partition_sketches[best_p].merge(current_sketch); corrected_estimate_per_part = std::max(corrected_estimate_per_part, partition_sketches[best_p].estimate()); - - return true; } /*!\brief Distributes the merged-bin candidates onto `number_of_remaining_tbs` partitions by LSH clustering and From 2e16669586d4d7c160769e0691bd81c028917d82 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 22:14:15 +0200 Subject: [PATCH 23/35] refactor: clean up Cluster * Document Cluster instead of the "foo" placeholder, including that is_valid accepts both valid and moved clusters. * Make the members private; nothing derives from Cluster. * Make the single-argument constructor explicit. * Replace the comments above LSH_find_representative_cluster and in very_similar_LSH_clustering, which described an older representation (cluster[i][0] == i). Co-Authored-By: Claude Opus 5.5 --- .../chopper/layout/fast_layout_cluster.hpp | 22 ++++++++++++------- src/layout/partition_user_bins.cpp | 15 ++++--------- 2 files changed, 18 insertions(+), 19 deletions(-) diff --git a/include/chopper/layout/fast_layout_cluster.hpp b/include/chopper/layout/fast_layout_cluster.hpp index b50b6bad..ce6e0c13 100644 --- a/include/chopper/layout/fast_layout_cluster.hpp +++ b/include/chopper/layout/fast_layout_cluster.hpp @@ -16,11 +16,20 @@ namespace chopper::layout { -/*\brief foo +/*!\brief A cluster of user bins, used by the LSH clustering of the fast layout. + * + * In a vector of clusters in which the cluster at position `i` was created with id `i`, the cluster at position `i` + * keeps `id() == i` and is either + * - **valid**: it has not been moved and contains at least one user bin, or + * - **moved**: its user bins were moved into another cluster (move_to), it is empty, and moved_to_cluster_id() is the + * id of that cluster. That cluster may itself have been moved, so moves can form a chain. + * LSH_find_representative_cluster follows the chain to the valid cluster. + * + * Note that `is_valid(i)` is true in both states: it checks that `id() == i` and that the cluster is in one of them. */ struct Cluster { -protected: +private: size_t representative_id{}; // representative id of the cluster; identifier; std::vector user_bins{}; // the user bins contained in thus cluster @@ -38,7 +47,7 @@ struct Cluster Cluster(size_t const id, size_t const user_bins_id) : representative_id{id}, user_bins({user_bins_id}) {} - Cluster(size_t const id) : Cluster{id, id} + explicit Cluster(size_t const id) : Cluster{id, id} {} size_t id() const @@ -119,11 +128,8 @@ struct Cluster } }; -// A valid cluster is one that hasn't been moved but actually contains user bins -// A valid cluster at position i is identified by the following equality: cluster[i].size() >= 1 && cluster[i][0] == i -// A moved cluster is one that has been joined and thereby moved to another cluster -// A moved cluster i is identified by the following: cluster[i].size() == 1 && cluster[i][0] != i -// returns position of the representative cluster +// Follows the chain of moves, starting at clusters[current_id], and returns the position of the representative +// cluster, i.e., the valid cluster that holds the user bins now. See Cluster for valid and moved clusters. inline size_t LSH_find_representative_cluster(std::vector const & clusters, size_t current_id) { std::reference_wrapper representative = clusters[current_id]; diff --git a/src/layout/partition_user_bins.cpp b/src/layout/partition_user_bins.cpp index d3821fd9..4535591c 100644 --- a/src/layout/partition_user_bins.cpp +++ b/src/layout/partition_user_bins.cpp @@ -136,17 +136,10 @@ std::vector very_similar_LSH_clustering(std::vector= minHash_sketch_size); assert(number_of_user_bins <= minHash_sketches.size()); - // initialise clusters with a signle user bin per cluster. - // clusters are either - // 1) of size 1; containing an id != position where the id points to the cluster it has been moved to - // e.g. cluster[Y] = {Z} (Y has been moved into Z, so Z could look likes this cluster[Z] = {Z, Y}) - // 2) of size >= 1; with the first entry beging id == position (a valid cluster) - // e.g. cluster[X] = {X} // valid singleton - // e.g. cluster[X] = {X, a, b, c, ...} // valid cluster with more joined entries - // The clusters could me moved recursively, s.t. - // cluster[A] = {B} - // cluster[B] = {C} - // cluster[C] = {C, A, B} // is valid cluster since cluster[C][0] == C; contains A and B + // Initialise one cluster per user bin. Merging moves the user bins of a cluster into another cluster, e.g., if + // clusters[A] is moved into clusters[B] and clusters[B] into clusters[C], clusters[A] and clusters[B] are empty + // and point to B and C, respectively, while clusters[C] holds the user bins of A, B and C. + // See Cluster for valid and moved clusters. std::vector clusters; clusters.reserve(number_of_user_bins); From 821d265ada8822016c03013a425b659fe980d738 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 22:15:29 +0200 Subject: [PATCH 24/35] refactor: use insert in Cluster::move_to unconditionally append_range and insert do the same for a std::vector, but each toolchain only compiled and tested one of the two branches. Keep insert. Co-Authored-By: Claude Opus 5.5 --- include/chopper/layout/fast_layout_cluster.hpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/include/chopper/layout/fast_layout_cluster.hpp b/include/chopper/layout/fast_layout_cluster.hpp index ce6e0c13..ab731164 100644 --- a/include/chopper/layout/fast_layout_cluster.hpp +++ b/include/chopper/layout/fast_layout_cluster.hpp @@ -107,11 +107,7 @@ struct Cluster { auto & target = target_cluster.user_bins; auto & source = this->user_bins; -#ifdef __cpp_lib_containers_ranges - target.append_range(source); -#else target.insert(target.end(), source.cbegin(), source.cend()); -#endif source = std::vector{}; // .clear() AND release memory moved_id = target_cluster.id(); From f588a36395811dfb10daedbbe30ada476422c059 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 22:17:28 +0200 Subject: [PATCH 25/35] refactor: tidy up determine_split_bins * Document determine_split_bins. * Remove commented-out debug output and the unused max_id. * Remove the unused print_matrix.hpp include, include what is used instead of , and fix the namespace in the file brief. Co-Authored-By: Claude Opus 5.5 --- .../chopper/layout/determine_split_bins.hpp | 28 +++++++++++++++++-- src/layout/determine_split_bins.cpp | 17 ++++------- 2 files changed, 32 insertions(+), 13 deletions(-) diff --git a/include/chopper/layout/determine_split_bins.hpp b/include/chopper/layout/determine_split_bins.hpp index 82b5ac0a..d0b173d6 100644 --- a/include/chopper/layout/determine_split_bins.hpp +++ b/include/chopper/layout/determine_split_bins.hpp @@ -6,12 +6,14 @@ // -------------------------------------------------------------------------------------------------- /*!\file - * \brief Provides chopper::determine_split_bins. + * \brief Provides chopper::layout::determine_split_bins. * \author Svenja Mehringer */ #pragma once +#include +#include #include #include @@ -19,6 +21,28 @@ namespace chopper::layout { +/*!\brief Distributes the largest user bins as split bins over technical bins at the back of `partitions`. + * \param[in] config The configuration (uses `maximum_fpr` and `number_of_hash_functions` of + * `config.hibf_config`). + * \param[in] positions User bin indices, sorted by descending cardinality. The first `num_user_bins` + * are split. + * \param[in] cardinalities The cardinality of each user bin, indexed by the values in `positions`. + * \param[in] num_technical_bins The maximum number of technical bins for the split user bins. Must be at least + * `num_user_bins`. + * \param[in] num_user_bins The number of user bins to split. + * \param[in,out] partitions The technical bins. Split user bin indices are appended to the last technical + * bins, one index per technical bin. `positions[num_user_bins - 1]` gets the last + * technical bins, `positions[0]` the first of them. The technical bins of each user + * bin are consecutive. + * \returns A pair of + * 1. the number of technical bins used, which may be less than `num_technical_bins`, and + * 2. the largest FPR-corrected cardinality per technical bin. + * `{0, 0}` if `num_user_bins == 0`. + * + * A dynamic programming algorithm assigns each user bin a number of technical bins, such that the largest + * FPR-corrected cardinality per technical bin, `ceil(cardinality * fpr_correction[k] / k)` for a user bin in `k` + * technical bins, is minimal. If several numbers of technical bins give the same minimum, the smallest is used. + */ std::pair determine_split_bins(chopper::configuration const & config, std::vector const & positions, std::vector const & cardinalities, @@ -26,4 +50,4 @@ std::pair determine_split_bins(chopper::configuration const & co size_t const num_user_bins, std::vector> & partitions); -} +} // namespace chopper::layout diff --git a/src/layout/determine_split_bins.cpp b/src/layout/determine_split_bins.cpp index fb4ba653..c4e85f3d 100644 --- a/src/layout/determine_split_bins.cpp +++ b/src/layout/determine_split_bins.cpp @@ -5,14 +5,17 @@ // shipped with this file and also available at: https://github.com/seqan/chopper/blob/main/LICENSE.md // --------------------------------------------------------------------------------------------------- -#include // for allocator, string -#include // for vector +#include +#include +#include +#include +#include +#include #include #include #include -#include // for data_store #include namespace chopper::layout @@ -71,8 +74,6 @@ std::pair determine_split_bins(chopper::configuration const & co size_t score = std::max(seqan::hibf::divide_and_ceil(corrected_ub_cardinality, i - i_prime), matrix[i_prime][j - 1]); - // std::cout << "j:" << j << " i:" << i << " i':" << i_prime << " score:" << score << std::endl; - minimum = (score < minimum) ? (trace[i][j] = i_prime, score) : minimum; } @@ -80,9 +81,6 @@ std::pair determine_split_bins(chopper::configuration const & co } } - // seqan::hibf::layout::print_matrix(matrix, num_technical_bins, num_user_bins, std::numeric_limits::max()); - //seqan::hibf::layout::print_matrix(trace, num_technical_bins, num_user_bins, std::numeric_limits::max()); - // backtracking // first, in the last column, find the row with minimum score (it can happen that the last rows are equally good) // the less rows, the better @@ -101,7 +99,6 @@ std::pair determine_split_bins(chopper::configuration const & co // now that we found the best trace_i start usual backtracking size_t trace_j = num_user_bins - 1; - // size_t max_id{}; size_t max_size{}; size_t bin_id{}; @@ -116,7 +113,6 @@ std::pair determine_split_bins(chopper::configuration const & co if (cardinality_per_bin > max_size) { - // max_id = bin_id; max_size = cardinality_per_bin; } @@ -137,7 +133,6 @@ std::pair determine_split_bins(chopper::configuration const & co if (cardinality_per_bin > max_size) { - // max_id = bin_id; max_size = cardinality_per_bin; } From bfd8fad0b6550a2774b7d039918e10e7d8fcd605 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 22:18:41 +0200 Subject: [PATCH 26/35] refactor: remove the unused timer in general_layout general_layout started and stopped a local dp_algorithm_timer that was never read. Return the result of compute_layout directly. Co-Authored-By: Claude Opus 5.5 --- src/layout/fast_layout.cpp | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/src/layout/fast_layout.cpp b/src/layout/fast_layout.cpp index dce9810a..b02c1cdd 100644 --- a/src/layout/fast_layout.cpp +++ b/src/layout/fast_layout.cpp @@ -38,22 +38,15 @@ seqan::hibf::layout::layout general_layout(chopper::configuration const & config std::vector const & cardinalities, std::vector const & sketches) { - seqan::hibf::layout::layout hibf_layout; - seqan::hibf::concurrent_timer union_estimation_timer{}; seqan::hibf::concurrent_timer rearrangement_timer{}; - seqan::hibf::concurrent_timer dp_algorithm_timer{}; - - dp_algorithm_timer.start(); - hibf_layout = seqan::hibf::layout::compute_layout(config.hibf_config, - cardinalities, - sketches, - std::move(positions), - union_estimation_timer, - rearrangement_timer); - dp_algorithm_timer.stop(); - - return hibf_layout; + + return seqan::hibf::layout::compute_layout(config.hibf_config, + cardinalities, + sketches, + std::move(positions), + union_estimation_timer, + rearrangement_timer); } /*!\brief Decides whether the lower-level IBF for a merged bin is laid out by the fast layout or by general_layout. From 1828a12e17b39551f9d914a415ace05cb1744bdb Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 22:19:49 +0200 Subject: [PATCH 27/35] refactor: use push_back for the remaining clusters remaining_clusters.insert(remaining_clusters.end(), x) is push_back(x). Co-Authored-By: Claude Opus 5.5 --- src/layout/partition_user_bins.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/layout/partition_user_bins.cpp b/src/layout/partition_user_bins.cpp index 4535591c..9183af04 100644 --- a/src/layout/partition_user_bins.cpp +++ b/src/layout/partition_user_bins.cpp @@ -574,7 +574,7 @@ size_t lsh_sim_approach(chopper::configuration const & config, if (clusters[i].empty()) break; - remaining_clusters.insert(remaining_clusters.end(), clusters[i].contained_user_bins()); + remaining_clusters.push_back(clusters[i].contained_user_bins()); } // assign the rest by similarity From d216d8aa558e6c4e046e461340e4f18415471490 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 22:21:07 +0200 Subject: [PATCH 28/35] refactor: rename the positions parameter of find_best_partition The parameter holds the user bins per partition, and the caller passes partitions. Name it partitions. Co-Authored-By: Claude Opus 5.5 --- src/layout/partition_user_bins.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/layout/partition_user_bins.cpp b/src/layout/partition_user_bins.cpp index 9183af04..552f64ad 100644 --- a/src/layout/partition_user_bins.cpp +++ b/src/layout/partition_user_bins.cpp @@ -291,7 +291,7 @@ void post_process_clusters(std::vector & clusters, * \param[in] cluster The global user bin indices to assign together. * \param[in] cardinalities The cardinality of each user bin, indexed by global user bin index. * \param[in] sketches The HyperLogLog sketch of each user bin, indexed by global user bin index. - * \param[in,out] positions The user bins per partition. `cluster` is appended to the chosen one. + * \param[in,out] partitions The user bins per partition. `cluster` is appended to the chosen one. * \param[in,out] partition_sketches The union sketch per partition. Updated for the chosen partition. * \param[in,out] max_partition_cardinality The largest user bin cardinality per partition. Updated. * \param[in,out] min_partition_cardinality The smallest user bin cardinality per partition. Updated. @@ -315,7 +315,7 @@ void find_best_partition(chopper::configuration const & config, std::vector const & cluster, std::vector const & cardinalities, std::vector const & sketches, - std::vector> & positions, + std::vector> & partitions, std::vector & partition_sketches, std::vector & max_partition_cardinality, std::vector & min_partition_cardinality) @@ -350,21 +350,21 @@ void find_best_partition(chopper::configuration const & config, auto penalty_lower_level = [&](size_t const additional_number_of_user_bins, size_t const p) -> size_t { - assert(positions[p].size() != 0); // partitions should be initialised beforehand + assert(partitions[p].size() != 0); // partitions should be initialised beforehand size_t const min = min_partition_cardinality[p]; size_t const max = max_partition_cardinality[p]; - if (positions[p].size() > config.hibf_config.tmax) // already a third level + if (partitions[p].size() > config.hibf_config.tmax) // already a third level { // if there must already be another lower level because the current merged bin contains more than tmax // user bins, then the current user bin is very likely stored multiple times. Therefore, the penalty is set // to the cardinality of the current user bin times the number of levels, e.g. the number of times this user // bin needs to be stored additionally - size_t const num_ubs_in_merged_bin{positions[p].size() + additional_number_of_user_bins}; + size_t const num_ubs_in_merged_bin{partitions[p].size() + additional_number_of_user_bins}; double const levels = std::log(num_ubs_in_merged_bin) / std::log(config.hibf_config.tmax); return static_cast(max_card * levels); } - else if (positions[p].size() + additional_number_of_user_bins > config.hibf_config.tmax) // now a third level + else if (partitions[p].size() + additional_number_of_user_bins > config.hibf_config.tmax) // now a third level { // if the current merged bin contains exactly tmax UBS, adding otherone must // result in another lower level. Most likely, the smallest user bin will end up on the lower level @@ -374,7 +374,7 @@ void find_best_partition(chopper::configuration const & config, size_t const penalty = std::min(min, max_card) * config.hibf_config.tmax; return penalty; } - else // positions[p].size() + additional_number_of_user_bins <= tmax + else // partitions[p].size() + additional_number_of_user_bins <= tmax { // if the new user bin is smaller than all other already contained user bins // the waste of space is high if stored in a single technical bin @@ -415,7 +415,7 @@ void find_best_partition(chopper::configuration const & config, // now that we know which partition fits best (`best_p`), add those indices to it for (size_t const user_bin_idx : cluster) { - positions[best_p].push_back(user_bin_idx); + partitions[best_p].push_back(user_bin_idx); max_partition_cardinality[best_p] = std::max(max_partition_cardinality[best_p], cardinalities[user_bin_idx]); min_partition_cardinality[best_p] = std::min(min_partition_cardinality[best_p], cardinalities[user_bin_idx]); } From 7eeee80ab933a2e654ffa9305de7e457f2213446 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 22:23:03 +0200 Subject: [PATCH 29/35] fix: typo in the initial partition timer Rename intital_partition_timer to initial_partition_timer. This also changes the column header in the --timing-output file to initial_partition_timer_in_seconds. Co-Authored-By: Claude Opus 5.5 --- include/chopper/configuration.hpp | 2 +- src/chopper_layout.cpp | 4 ++-- src/layout/fast_layout.cpp | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/include/chopper/configuration.hpp b/include/chopper/configuration.hpp index 7a1bfad0..c8d42e29 100644 --- a/include/chopper/configuration.hpp +++ b/include/chopper/configuration.hpp @@ -82,7 +82,7 @@ struct configuration mutable seqan::hibf::concurrent_timer dp_algorithm_timer{}; mutable seqan::hibf::concurrent_timer lsh_algorithm_timer{}; mutable seqan::hibf::concurrent_timer search_partition_algorithm_timer{}; - mutable seqan::hibf::concurrent_timer intital_partition_timer{}; + mutable seqan::hibf::concurrent_timer initial_partition_timer{}; mutable seqan::hibf::concurrent_timer small_layouts_timer{}; void read_from(std::istream & stream); diff --git a/src/chopper_layout.cpp b/src/chopper_layout.cpp index c193211b..28027989 100644 --- a/src/chopper_layout.cpp +++ b/src/chopper_layout.cpp @@ -168,7 +168,7 @@ int chopper_layout(chopper::configuration & config, sharg::parser & parser) << "union_estimation_in_seconds\t" << "rearrangement_in_seconds\t" << "lsh_in_seconds\t" - << "intital_partition_timer_in_seconds\t" + << "initial_partition_timer_in_seconds\t" << "small_layouts_timer_in_seconds\t" << "search_best_p_in_seconds\n"; output_stream << config.compute_sketches_timer.in_seconds() << '\t'; @@ -176,7 +176,7 @@ int chopper_layout(chopper::configuration & config, sharg::parser & parser) output_stream << config.union_estimation_timer.in_seconds() << '\t'; output_stream << config.rearrangement_timer.in_seconds() << '\t'; output_stream << config.lsh_algorithm_timer.in_seconds() << '\t'; - output_stream << config.intital_partition_timer.in_seconds() << '\t'; + output_stream << config.initial_partition_timer.in_seconds() << '\t'; output_stream << config.small_layouts_timer.in_seconds() << '\t'; output_stream << config.search_partition_algorithm_timer.in_seconds() << '\n'; } diff --git a/src/layout/fast_layout.cpp b/src/layout/fast_layout.cpp index b02c1cdd..1342fe85 100644 --- a/src/layout/fast_layout.cpp +++ b/src/layout/fast_layout.cpp @@ -315,9 +315,9 @@ void fast_layout(chopper::configuration const & config, std::vector> tmax_partitions(config.hibf_config.tmax); // here we assume that we want to start with a fast layout - config.intital_partition_timer.start(); + config.initial_partition_timer.start(); partition_user_bins(config, positions, cardinalities, sketches, minHash_sketches, tmax_partitions); - config.intital_partition_timer.stop(); + config.initial_partition_timer.stop(); // initialise user bins in layout hibf_layout.user_bins.resize(config.hibf_config.number_of_user_bins); From ddcb2bd8982abb7cccb966a35698273f4bbb8b53 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 22:52:42 +0200 Subject: [PATCH 30/35] fix: throw in execute for determine_best_tmax with fast_layout determine_best_number_of_technical_bins always uses the DP layout, so execute silently ignored fast_layout if determine_best_tmax was set. The command line already rejects the combination; the API now throws std::invalid_argument. The two execute_estimation_test cases with both options set expected the DP output. They are kept under #if 0 until the combination is supported, and same-named tests check the throw instead. Co-Authored-By: Claude Opus 5.5 --- src/layout/execute.cpp | 4 +- .../layout/execute_with_estimation_test.cpp | 74 +++++++++++++++++++ 2 files changed, 77 insertions(+), 1 deletion(-) diff --git a/src/layout/execute.cpp b/src/layout/execute.cpp index b3d1c70b..d18be871 100644 --- a/src/layout/execute.cpp +++ b/src/layout/execute.cpp @@ -34,6 +34,9 @@ int execute(chopper::configuration & config, std::vector const & sketches, std::vector const & minHash_sketches) { + if (config.determine_best_tmax && config.fast_layout) + throw std::invalid_argument{"determine_best_tmax is not supported with fast_layout."}; + config.hibf_config.validate_and_set_defaults(); std::vector cardinalities; @@ -43,7 +46,6 @@ int execute(chopper::configuration & config, if (config.determine_best_tmax) { - // Always uses the DP layout; config.fast_layout is ignored. chopper_layout rejects the combination. hibf_layout = determine_best_number_of_technical_bins(config, cardinalities, sketches); } else diff --git a/test/api/layout/execute_with_estimation_test.cpp b/test/api/layout/execute_with_estimation_test.cpp index 518a2a60..913ee19d 100644 --- a/test/api/layout/execute_with_estimation_test.cpp +++ b/test/api/layout/execute_with_estimation_test.cpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include #include @@ -84,6 +85,7 @@ TEST(execute_estimation_test, few_ubs) )expected_cout"); } +#if 0 // determine_best_tmax + fast_layout currently not supported TEST(execute_estimation_test, few_ubs_fast_layout) { seqan3::test::tmp_directory tmp_dir{}; @@ -141,6 +143,39 @@ TEST(execute_estimation_test, few_ubs_fast_layout) # Best t_max (regarding expected query runtime): 64 )expected_cout"); } +#else +TEST(execute_estimation_test, few_ubs_fast_layout) +{ + seqan3::test::tmp_directory tmp_dir{}; + std::filesystem::path const layout_file{tmp_dir.path() / "layout.tsv"}; + + auto simulated_input = [&](size_t const num, seqan::hibf::insert_iterator it) + { + size_t const desired_kmer_count = (num == 1) ? 1700 : 1200; + for (auto hash : std::views::iota(0u, desired_kmer_count)) + it = hash; + }; + + chopper::configuration config{}; + config.fast_layout = true; + config.hibf_config.tmax = 64; + config.hibf_config.input_fn = simulated_input; + config.hibf_config.number_of_user_bins = 8; + config.determine_best_tmax = true; + config.disable_sketch_output = true; + config.output_filename = layout_file; + config.hibf_config.disable_estimate_union = true; // also disables rearrangement + + std::vector> + filenames{{"seq0"}, {"seq1"}, {"seq2"}, {"seq3"}, {"seq4"}, {"seq5"}, {"seq6"}, {"seq7"}}; + + std::vector sketches; + std::vector minHash_sketches{}; + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); + + EXPECT_THROW(chopper::layout::execute(config, filenames, sketches, minHash_sketches), std::invalid_argument); +} +#endif TEST(execute_estimation_test, many_ubs) { @@ -459,6 +494,7 @@ TEST(execute_estimation_test, many_ubs) EXPECT_EQ(actual_file, expected_file) << actual_file; } +#if 0 // determine_best_tmax + fast_layout currently not supported TEST(execute_estimation_test, many_ubs_fast_layout) { seqan3::test::tmp_directory tmp_dir{}; @@ -776,6 +812,44 @@ TEST(execute_estimation_test, many_ubs_fast_layout) std::string const actual_file{string_from_file(layout_file)}; EXPECT_EQ(actual_file, expected_file) << actual_file; } +#else +TEST(execute_estimation_test, many_ubs_fast_layout) +{ + seqan3::test::tmp_directory tmp_dir{}; + std::filesystem::path const layout_file{tmp_dir.path() / "layout.tsv"}; + + std::vector> many_filenames; + + for (size_t i{0}; i < 96u; ++i) + many_filenames.push_back({seqan3::detail::to_string("seq", i)}); + + // Creates sizes of the following series + // [801,802,...,820,922,923,...,941,1043,1044,...,1062,1164,1165,...,1183,1285,1286,...,1300] + // See also https://godbolt.org/z/9517eaaaG + auto simulated_input = [&](size_t const num, seqan::hibf::insert_iterator it) + { + size_t const desired_kmer_count = 101 * ((num + 20) / 20) + num + 700; + for (auto hash : std::views::iota(0u, desired_kmer_count)) + it = hash; + }; + + chopper::configuration config{}; + config.fast_layout = true; + config.determine_best_tmax = true; + config.output_filename = layout_file; + config.disable_sketch_output = true; + config.hibf_config.tmax = 1024; + config.hibf_config.input_fn = simulated_input; + config.hibf_config.number_of_user_bins = many_filenames.size(); + config.hibf_config.disable_estimate_union = true; // also disables rearrangement + + std::vector sketches; + std::vector minHash_sketches{}; + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); + + EXPECT_THROW(chopper::layout::execute(config, many_filenames, sketches, minHash_sketches), std::invalid_argument); +} +#endif TEST(execute_estimation_test, many_ubs_force_all) { From d7a03d6fe8f104d43074b242ea40f21c87d97f64 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 22:56:31 +0200 Subject: [PATCH 31/35] docs: document the members of Cluster Co-Authored-By: Claude Opus 5.5 --- .../chopper/layout/fast_layout_cluster.hpp | 56 +++++++++++++++---- 1 file changed, 46 insertions(+), 10 deletions(-) diff --git a/include/chopper/layout/fast_layout_cluster.hpp b/include/chopper/layout/fast_layout_cluster.hpp index ab731164..a5cdb54f 100644 --- a/include/chopper/layout/fast_layout_cluster.hpp +++ b/include/chopper/layout/fast_layout_cluster.hpp @@ -30,51 +30,71 @@ namespace chopper::layout struct Cluster { private: - size_t representative_id{}; // representative id of the cluster; identifier; + //!\brief The id of the cluster. + size_t representative_id{}; - std::vector user_bins{}; // the user bins contained in thus cluster + //!\brief The user bins contained in this cluster. + std::vector user_bins{}; - std::optional moved_id{std::nullopt}; // where this Clusters user bins where moved to + //!\brief The id of the cluster that this cluster's user bins were moved to, if they were moved. + std::optional moved_id{std::nullopt}; public: - Cluster() = default; - Cluster(Cluster const &) = default; - Cluster(Cluster &&) = default; - Cluster & operator=(Cluster const &) = default; - Cluster & operator=(Cluster &&) = default; - ~Cluster() = default; - + /*!\name Constructors, destructor and assignment + * \{ + */ + Cluster() = default; //!< Defaulted. + Cluster(Cluster const &) = default; //!< Defaulted. + Cluster(Cluster &&) = default; //!< Defaulted. + Cluster & operator=(Cluster const &) = default; //!< Defaulted. + Cluster & operator=(Cluster &&) = default; //!< Defaulted. + ~Cluster() = default; //!< Defaulted. + + /*!\brief Creates a cluster with id `id` that contains the single user bin `user_bins_id`. + * \param[in] id The id of the cluster. + * \param[in] user_bins_id The user bin the cluster contains. + */ Cluster(size_t const id, size_t const user_bins_id) : representative_id{id}, user_bins({user_bins_id}) {} + /*!\brief Creates a cluster with id `id` that contains the single user bin `id`. + * \param[in] id The id of the cluster and the user bin it contains. + */ explicit Cluster(size_t const id) : Cluster{id, id} {} + //!\} + //!\brief Returns the id of the cluster. size_t id() const { return representative_id; } + //!\brief Returns the user bins contained in the cluster. std::vector const & contained_user_bins() const { return user_bins; } + //!\brief Returns whether the user bins of this cluster were moved to another cluster (see move_to). bool has_been_moved() const { return moved_id.has_value(); } + //!\brief Returns whether the cluster contains no user bins. bool empty() const { return user_bins.empty(); } + //!\brief Returns the number of user bins in the cluster. size_t size() const { return user_bins.size(); } + //!\brief Removes the last user bin from the cluster and returns it. The cluster must not be empty. size_t pop_back() { size_t last = user_bins.back(); @@ -82,11 +102,18 @@ struct Cluster return last; } + /*!\brief Adds a user bin to the cluster. + * \param[in] user_bin The user bin to add. + */ void add_user_bin(size_t const user_bin) { user_bins.push_back(user_bin); } + /*!\brief Checks that the cluster at position `id` is either valid or moved, see Cluster. + * \param[in] id The position of the cluster, i.e., the id it must have. + * \returns `true` if `id() == id`, and the cluster is either not moved and non-empty, or moved and empty. + */ bool is_valid(size_t const id) const { bool const ids_equal = representative_id == id; @@ -96,6 +123,7 @@ struct Cluster return ids_equal && (properly_moved || not_moved); } + //!\brief Returns the id of the cluster that this cluster's user bins were moved to. The cluster must be moved. size_t moved_to_cluster_id() const { assert(moved_id.has_value()); @@ -103,6 +131,11 @@ struct Cluster return moved_id.value(); } + /*!\brief Moves all user bins of this cluster to `target_cluster` and marks this cluster as moved. + * \param[in,out] target_cluster The cluster to move the user bins to. Must be a different cluster. + * + * Afterwards, this cluster is empty, its memory is released, and moved_to_cluster_id() is `target_cluster.id()`. + */ void move_to(Cluster & target_cluster) { auto & target = target_cluster.user_bins; @@ -113,6 +146,9 @@ struct Cluster moved_id = target_cluster.id(); } + /*!\brief Sorts the user bins of the cluster by descending cardinality. + * \param[in] cardinalities The cardinality of each user bin, indexed by user bin. + */ void sort_by_cardinality(std::vector const & cardinalities) { std::ranges::sort(user_bins, From 44da2b79d7f31c14ec801063c57c474e99c89eeb Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 22:57:36 +0200 Subject: [PATCH 32/35] docs: document LSH_find_representative_cluster with doxygen It had a plain comment, so it was missing from the documentation. Co-Authored-By: Claude Opus 5.5 --- include/chopper/layout/fast_layout_cluster.hpp | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/include/chopper/layout/fast_layout_cluster.hpp b/include/chopper/layout/fast_layout_cluster.hpp index a5cdb54f..b40c885f 100644 --- a/include/chopper/layout/fast_layout_cluster.hpp +++ b/include/chopper/layout/fast_layout_cluster.hpp @@ -160,8 +160,14 @@ struct Cluster } }; -// Follows the chain of moves, starting at clusters[current_id], and returns the position of the representative -// cluster, i.e., the valid cluster that holds the user bins now. See Cluster for valid and moved clusters. +/*!\brief Returns the position of the representative cluster of `clusters[current_id]`. + * \param[in] clusters The clusters. The cluster at position `i` must have id `i`, see Cluster. + * \param[in] current_id The position of the cluster to start from. + * \returns The position of the representative cluster, i.e., the valid cluster that holds the user bins of + * `clusters[current_id]` now. + * + * Follows the chain of moves, starting at `clusters[current_id]`. See Cluster for valid and moved clusters. + */ inline size_t LSH_find_representative_cluster(std::vector const & clusters, size_t current_id) { std::reference_wrapper representative = clusters[current_id]; From 27cc446bc9b2a11c0a69bc91b3ec61c814227d30 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 22:59:46 +0200 Subject: [PATCH 33/35] docs: document the fast layout timers Co-Authored-By: Claude Opus 5.5 --- include/chopper/configuration.hpp | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/include/chopper/configuration.hpp b/include/chopper/configuration.hpp index c8d42e29..8bb2e653 100644 --- a/include/chopper/configuration.hpp +++ b/include/chopper/configuration.hpp @@ -80,9 +80,22 @@ struct configuration mutable seqan::hibf::concurrent_timer union_estimation_timer{}; mutable seqan::hibf::concurrent_timer rearrangement_timer{}; mutable seqan::hibf::concurrent_timer dp_algorithm_timer{}; + /*!\brief Fast layout: time spent in LSH clustering (`lsh_in_seconds` in the timing output). + * + * Summed over all partitionings, including concurrent ones, so it can exceed the wall-clock time. + */ mutable seqan::hibf::concurrent_timer lsh_algorithm_timer{}; + /*!\brief Fast layout: time spent assigning clusters to partitions by similarity (`search_best_p_in_seconds` in + * the timing output). + * + * Summed over all partitionings, including concurrent ones, so it can exceed the wall-clock time. + */ mutable seqan::hibf::concurrent_timer search_partition_algorithm_timer{}; + /*!\brief Fast layout: time for partitioning the top level (`initial_partition_timer_in_seconds` in the timing + * output). + */ mutable seqan::hibf::concurrent_timer initial_partition_timer{}; + //!\brief Fast layout: time for laying out all lower levels (`small_layouts_timer_in_seconds` in the timing output). mutable seqan::hibf::concurrent_timer small_layouts_timer{}; void read_from(std::istream & stream); From e6fb149a6ef063d7dfef1ba5c806cc1af21ed1c2 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Tue, 29 Sep 2026 23:01:11 +0200 Subject: [PATCH 34/35] docs: document execute Co-Authored-By: Claude Opus 5.5 --- include/chopper/layout/execute.hpp | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/include/chopper/layout/execute.hpp b/include/chopper/layout/execute.hpp index 00c4fa2e..f3e4c161 100644 --- a/include/chopper/layout/execute.hpp +++ b/include/chopper/layout/execute.hpp @@ -18,6 +18,24 @@ namespace chopper::layout { +/*!\brief Computes the layout for the given user bins and writes it to `config.output_filename`. + * \param[in,out] config The configuration. `config.hibf_config` is validated and completed + * (`validate_and_set_defaults`), and the timers are updated. + * \param[in] filenames The file names of each user bin. They are written to the layout file. + * \param[in] sketches The HyperLogLog sketch of each user bin. + * \param[in] minHash_sketches The MinHash sketches of each user bin. Only used, and then required for every user + * bin, if `config.fast_layout` is set. May be empty otherwise. + * \returns 0. + * \throws std::invalid_argument If both `config.determine_best_tmax` and `config.fast_layout` are set. + * + * The layout is computed with + * - determine_best_number_of_technical_bins if `config.determine_best_tmax` is set, + * - fast_layout if `config.fast_layout` is set, + * - the DP layout of the HIBF library (seqan::hibf::layout::compute_layout) otherwise. + * + * Unless `config.determine_best_tmax` is set, `config.output_verbose_statistics` prints statistics of the layout to + * `std::cout`. + */ int execute(chopper::configuration & config, std::vector> const & filenames, std::vector const & sketches, From 3acabe6d9c3e6f643a132d6183296209755ab722 Mon Sep 17 00:00:00 2001 From: Enrico Seiler Date: Wed, 30 Sep 2026 00:00:12 +0200 Subject: [PATCH 35/35] test: add the measurement of the post_process_clusters seeding variants post_process_clusters sorts the first tmax clusters by size, but lsh_sim_approach seeds only tmax - number_of_split_tbs partitions. The alternative, sorting the first number_of_remaining_tbs clusters, was measured on 60 synthetic inputs (400 to 50000 user bins). There was no systematic difference in expected HIBF size or query cost, so chopper keeps the current order. The directory contains the write-up, the measurement code (driver, patch for the alternative, build, run and analysis scripts) and the results. Co-Authored-By: Claude Opus 5.5 --- .../post_process_clusters_variants/README.md | 118 +++++++++ .../post_process_clusters_variants/analyse.py | 76 ++++++ .../post_process_clusters_variants/build.sh | 48 ++++ .../measure.cpp | 230 ++++++++++++++++++ .../results_large.tsv | 25 ++ .../results_small.tsv | 37 +++ .../post_process_clusters_variants/run.sh | 39 +++ .../variant_b.patch | 52 ++++ 8 files changed, 625 insertions(+) create mode 100644 test/benchmark/benchmark_data/post_process_clusters_variants/README.md create mode 100755 test/benchmark/benchmark_data/post_process_clusters_variants/analyse.py create mode 100755 test/benchmark/benchmark_data/post_process_clusters_variants/build.sh create mode 100644 test/benchmark/benchmark_data/post_process_clusters_variants/measure.cpp create mode 100644 test/benchmark/benchmark_data/post_process_clusters_variants/results_large.tsv create mode 100644 test/benchmark/benchmark_data/post_process_clusters_variants/results_small.tsv create mode 100755 test/benchmark/benchmark_data/post_process_clusters_variants/run.sh create mode 100644 test/benchmark/benchmark_data/post_process_clusters_variants/variant_b.patch diff --git a/test/benchmark/benchmark_data/post_process_clusters_variants/README.md b/test/benchmark/benchmark_data/post_process_clusters_variants/README.md new file mode 100644 index 00000000..e33070f6 --- /dev/null +++ b/test/benchmark/benchmark_data/post_process_clusters_variants/README.md @@ -0,0 +1,118 @@ +## Seeding clusters in `post_process_clusters` + +Measurement behind the decision to keep how `post_process_clusters` (`src/layout/partition_user_bins.cpp`) orders the +clusters for the fast layout. + +### Question + +`lsh_sim_approach` seeds the merged technical bins with the first clusters, then assigns all other clusters greedily by +similarity (`find_best_partition`). Before that, `post_process_clusters` orders the clusters: + +1. The first `tmax` positions receive the clusters with the most user bins (partial sort). +2. The rest is sorted by the cardinality of each cluster's largest user bin, so that clusters with user bins of similar + size are assigned after each other and small ones last. + +If there are split bins, only `tmax - number_of_split_tbs` technical bins are left for merged bins and seeded. The +clusters at positions `[tmax - number_of_split_tbs, tmax)` are then assigned first, in order of their number of user +bins, not by cardinality as step 2 intends. + +* **Variant A** (chopper): partial sort of the first `tmax` clusters. +* **Variant B**: partial sort of the first `number_of_remaining_tbs = tmax - number_of_split_tbs` clusters. + +The variants differ only if there are split bins. Because the second calibration pass of `partition_user_bins` depends +on the first pass, a difference can remain even if the final partitioning has no split bins. + +**Decision: keep variant A.** No systematic difference was found (see [Results](#results)). + +### Method + +`measure.cpp` generates one synthetic input, computes its sketches once and lays it out with both variants. The variant +is switched at run time by `chopper_variant_b`, which `variant_b.patch` adds to a copy of `partition_user_bins.cpp`. +With `chopper_variant_b = false`, the patched copy produces the same layouts as chopper. + +**Inputs.** Each input imitates a collection of related genomes: + +* User bins belong to families. Family k-mer pools have log-normal sizes (`median`, `sigma`), clamped to `[min, max]`. +* A user bin contains a random 60 % to 98 % of its family's pool, plus `(1 - share) * pool + 500` unique k-mers. +* `families = n / family_divisor`. `family_divisor = 1` means every user bin has its own family. +* `tmax` is chosen as by the command line interface (`seqan::hibf::config::validate_and_set_defaults`). + +The inputs are deterministic. Running an input twice gives identical results. + +| Grid | n | family_divisor | sigma | seeds | median, min, max (k-mers) | tmax | inputs | +|---|---|---|---|---|---|---|---| +| small | 400, 1000 | 1, 4, 25 | 0.7, 1.4 | 3 | 8000, 5000, 400000 | 64 | 36 | +| large | 10000, 25000, 50000 | 5, 50 | 1.4, 2.0 | 3, 2, 1 | 4000, 2500, 2000000 | 128, 192, 256 | 24 | + +The minimum size is due to the MinHash sketches: computing them can fail for user bins with fewer than about 2000 +k-mers. + +**Metrics** (lower is better), per variant: + +* `topmax`: the largest FPR-corrected technical bin of the top-level partitioning. It determines the size of the + top-level IBF. +* `size`: expected total size of the HIBF in bytes (`hibf_statistics::total_hibf_size_in_byte`). +* `cost`: expected query cost of the HIBF (`hibf_statistics::expected_HIBF_query_cost`). + +Also recorded: `s`, the number of top-level technical bins of user bins that span at least two technical bins, and +`same_top`, whether both variants produce the same top-level partitioning. + +**Not covered.** A merged bin is laid out recursively with the fast layout only if it has at least `64 * tmax` user +bins. In the small grid this is impossible (`n < 64 * tmax`), so the lower levels were laid out by the DP algorithm of +the HIBF library and only the top-level partitioning differs between the variants. For the large grid, this was not +checked. All inputs are synthetic. + +### Results + +Ratios B/A over the inputs whose top-level partitionings differ. p-values are two-sided; the Wilcoxon signed-rank test +is on log(B/A). The per-input values are in `results_small.tsv` and `results_large.tsv`; `analyse.py` prints them. + +**Small grid.** The top-level partitionings differ in 21 of 36 inputs. + +| Metric | B better | B worse | equal | geo-mean B/A | range B/A | sign test p | Wilcoxon p | +|---|--:|--:|--:|--:|--:|--:|--:| +| largest top-level technical bin | 1 | 2 | 18 | 1.0004 | 0.979 - 1.031 | 1.00 | 0.75 | +| expected HIBF size | 7 | 10 | 4 | 0.9996 | 0.982 - 1.008 | 0.63 | 0.58 | +| expected query cost | 11 | 6 | 4 | 0.9988 | 0.971 - 1.015 | 0.33 | 0.52 | + +**Large grid.** The top-level partitionings differ in only 9 of 24 inputs, and in none with 50000 user bins: with many +small user bins, the split threshold (joint size / `tmax`) is high and few user bins are split. + +| Metric | B better | B worse | equal | geo-mean B/A | range B/A | sign test p | Wilcoxon p | +|---|--:|--:|--:|--:|--:|--:|--:| +| largest top-level technical bin | 4 | 2 | 3 | 0.9806 | 0.905 - 1.012 | 0.69 | 0.16 | +| expected HIBF size | 3 | 6 | 0 | 1.0062 | 0.997 - 1.034 | 0.51 | 0.25 | +| expected query cost | 3 | 5 | 1 | 1.0119 | 0.995 - 1.089 | 0.73 | 0.46 | + +The largest effect is a single input (n = 10000, 200 families, sigma = 2.0, seed 3, s = 76): variant B shrinks the +largest top-level technical bin by 6.5 %, but the HIBF grows by 3.5 % and the query cost by 8.9 %. It dominates the +geometric means of the large grid. + +Over both grids (30 inputs that differ), B is better in size for 10 inputs and worse for 16 (sign test p = 0.33), and +better in query cost for 14 and worse for 11 (p = 0.69). The differences go in both directions, are mostly below 1 %, +and no direction is significant. There is no measurable reason to change variant A. + +### Reproduce + +Requires CMake, a C++ compiler supported by chopper, `patch` and Python 3 (SciPy optional, for the Wilcoxon test). + +```bash +CXX=clang++-23 ./build.sh /tmp/ppcv # configures a chopper Release build, builds /tmp/ppcv/measure +cd /tmp/ppcv && /path/to/this/directory/run.sh ./measure # results_small.tsv, results_large.tsv (about 30 min) +/path/to/this/directory/analyse.py results_small.tsv results_large.tsv +``` + +`variant_b.patch` applies to `src/layout/partition_user_bins.cpp` as of the commit that added this directory. If the +file changes, the patch may need to be adapted. + +`results_small.tsv` was produced by `run.sh`. `results_large.tsv` comes from an earlier run of the same input generator +before it was moved into `measure.cpp`; the first 12 of its 24 rows were checked against `run.sh` and are identical. The +timing columns (`t_sketch`, `t_A`, `t_B`, in seconds) depend on the machine. + +### Environment + +* chopper `e6fb149` (`final_fast_layout` branch), hibf `7f252fd`, Release build (`-O3 -DNDEBUG`). +* Compiler: Debian clang version 23.1.2, libstdc++. +* OS: Linux 7.2.8-WSL2-STABLE x86_64. +* CPU: AMD Ryzen 9 9950X, 16 cores, 32 threads. Memory: 49 GiB. +* Threads: 8 (small grid), 32 (large grid). diff --git a/test/benchmark/benchmark_data/post_process_clusters_variants/analyse.py b/test/benchmark/benchmark_data/post_process_clusters_variants/analyse.py new file mode 100755 index 00000000..d5368cc3 --- /dev/null +++ b/test/benchmark/benchmark_data/post_process_clusters_variants/analyse.py @@ -0,0 +1,76 @@ +#!/usr/bin/env python3 +# --------------------------------------------------------------------------------------------------- +# Copyright (c) 2006-2023, Knut Reinert & Freie Universität Berlin +# Copyright (c) 2016-2023, Knut Reinert & MPI für molekulare Genetik +# This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License +# shipped with this file and also available at: https://github.com/seqan/chopper/blob/main/LICENSE.md +# --------------------------------------------------------------------------------------------------- +"""Summarises the output of run.sh. See README.md. + +Usage: analyse.py results_small.tsv [results_large.tsv ...] + +For each file, prints the ratio variant B / variant A per input and, over the inputs whose top-level partitions differ, +the number of inputs where B is better/worse, the geometric mean and range of B/A, a two-sided sign test and, if SciPy +is installed, a two-sided Wilcoxon signed-rank test on log(B/A). Lower is better for all metrics. +""" + +import csv +import math +import sys + +try: + from scipy.stats import wilcoxon +except ImportError: + wilcoxon = None + +METRICS = [("topmax", "largest top-level technical bin"), ("size", "expected HIBF size"), ("cost", "expected query cost")] + + +def sign_test(better: int, worse: int) -> float: + """Exact two-sided sign test.""" + n = better + worse + if n == 0: + return 1.0 + tail = sum(math.comb(n, i) for i in range(min(better, worse) + 1)) / 2**n + return min(1.0, 2 * tail) + + +def summarise(rows: list[dict], title: str) -> None: + print(f"## {title}\n") + print("| n | families | sigma | seed | tmax | s | same | topmax B/A | size B/A | cost B/A |") + print("|--:|--:|--:|--:|--:|--:|:-:|--:|--:|--:|") + for r in rows: + ratios = " | ".join(f"{float(r[m + '_B']) / float(r[m + '_A']):.4f}" for m, _ in METRICS) + same = "yes" if r["same_top"] == "1" else "no" + print(f"| {r['n']} | {r['families']} | {r['sigma']} | {r['seed']} | {r['tmax']} | {r['s_A']} | {same} | {ratios} |") + + differing = [r for r in rows if r["same_top"] == "0"] + print(f"\nInputs: {len(rows)}. Top-level partitions differ in {len(differing)}.\n") + if not differing: + return + print("| Metric | B better | B worse | equal | geo-mean B/A | range B/A | sign test p | Wilcoxon p |") + print("|---|--:|--:|--:|--:|--:|--:|--:|") + for metric, name in METRICS: + logs = [math.log(float(r[metric + "_B"]) / float(r[metric + "_A"])) for r in differing] + nonzero = [x for x in logs if abs(x) > 1e-12] + better = sum(x < 0 for x in nonzero) + worse = sum(x > 0 for x in nonzero) + geo_mean = math.exp(sum(logs) / len(logs)) + wilcoxon_p = f"{wilcoxon(nonzero).pvalue:.2f}" if wilcoxon and nonzero else "n/a" + print( + f"| {name} | {better} | {worse} | {len(logs) - len(nonzero)} | {geo_mean:.4f} | " + f"{math.exp(min(logs)):.3f} - {math.exp(max(logs)):.3f} | {sign_test(better, worse):.2f} | {wilcoxon_p} |" + ) + print() + + +def main() -> None: + if len(sys.argv) < 2: + sys.exit(__doc__) + for path in sys.argv[1:]: + with open(path, newline="") as f: + summarise(list(csv.DictReader(f, delimiter="\t")), path) + + +if __name__ == "__main__": + main() diff --git a/test/benchmark/benchmark_data/post_process_clusters_variants/build.sh b/test/benchmark/benchmark_data/post_process_clusters_variants/build.sh new file mode 100755 index 00000000..e7db5d7d --- /dev/null +++ b/test/benchmark/benchmark_data/post_process_clusters_variants/build.sh @@ -0,0 +1,48 @@ +#!/usr/bin/env bash +# Builds the `measure` binary. See README.md. +# +# Usage: CXX= ./build.sh +# +# Configures a chopper Release build in /chopper, builds chopper_layout, applies variant_b.patch to +# a copy of src/layout/partition_user_bins.cpp and compiles measure.cpp with the flags that chopper uses for +# partition_user_bins.cpp. + +set -Eeuo pipefail + +if [[ $# -ne 1 ]]; then + echo "Usage: CXX= $0 " >&2 + exit 1 +fi + +SCRIPT_DIR=$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd) +CHOPPER_DIR=$(cd "${SCRIPT_DIR}/../../../.." && pwd) +OUT_DIR=$(mkdir -p "$1" && cd "$1" && pwd) +BUILD_DIR="${OUT_DIR}/chopper" + +cmake -S "${CHOPPER_DIR}" -B "${BUILD_DIR}" -DCMAKE_BUILD_TYPE=Release -DCMAKE_EXPORT_COMPILE_COMMANDS=ON +cmake --build "${BUILD_DIR}" --target chopper_layout --parallel + +PATCHED="${OUT_DIR}/partition_user_bins_variants.cpp" +patch --output="${PATCHED}" "${CHOPPER_DIR}/src/layout/partition_user_bins.cpp" < "${SCRIPT_DIR}/variant_b.patch" + +# The compiler and flags of partition_user_bins.cpp, without -o , -c and the source file. +mapfile -t COMMAND < <(python3 - "${BUILD_DIR}/compile_commands.json" <<'PYTHON' +import json, shlex, sys +entry = next(e for e in json.load(open(sys.argv[1])) if e["file"].endswith("src/layout/partition_user_bins.cpp")) +args = shlex.split(entry["command"]) +result, skip = [], False +for arg in args: + if skip: + skip = False + elif arg == "-o": + skip = True + elif arg != "-c" and arg != entry["file"]: + result.append(arg) +print("\n".join(arg for arg in result if not arg.endswith("ccache"))) +PYTHON +) + +"${COMMAND[@]}" "${SCRIPT_DIR}/measure.cpp" "${PATCHED}" -o "${OUT_DIR}/measure" \ + "${BUILD_DIR}/lib/libchopper_layout.a" "${BUILD_DIR}/lib/libchopper_shared.a" "${BUILD_DIR}/lib/libhibf.a" + +echo "Built ${OUT_DIR}/measure" diff --git a/test/benchmark/benchmark_data/post_process_clusters_variants/measure.cpp b/test/benchmark/benchmark_data/post_process_clusters_variants/measure.cpp new file mode 100644 index 00000000..4807a649 --- /dev/null +++ b/test/benchmark/benchmark_data/post_process_clusters_variants/measure.cpp @@ -0,0 +1,230 @@ +// --------------------------------------------------------------------------------------------------- +// Copyright (c) 2006-2023, Knut Reinert & Freie Universität Berlin +// Copyright (c) 2016-2023, Knut Reinert & MPI für molekulare Genetik +// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License +// shipped with this file and also available at: https://github.com/seqan/chopper/blob/main/LICENSE.md +// --------------------------------------------------------------------------------------------------- + +// Compares two variants of post_process_clusters (src/layout/partition_user_bins.cpp) on one synthetic input. +// Variant A sorts the first tmax clusters by size (as in chopper), variant B the first number_of_remaining_tbs. +// Must be linked with a copy of partition_user_bins.cpp that has variant_b.patch applied. See README.md. +// +// Usage: measure +// Prints one tab-separated line; run.sh prints the header. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +// Defined in the patched partition_user_bins.cpp. +extern bool chopper_variant_b; + +namespace +{ + +// splitmix64 finaliser: turns consecutive integers into well-distributed hashes. +uint64_t scramble(uint64_t x) +{ + x += 0x9e3779b97f4a7c15ULL; + x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9ULL; + x = (x ^ (x >> 27)) * 0x94d049bb133111ebULL; + return x ^ (x >> 31); +} + +// A user bin consists of a random subset of its family's k-mer pool plus unique k-mers. +struct user_bin_spec +{ + uint64_t family; + uint64_t pool_size; + uint32_t keep_threshold; // Keep a pool k-mer if a hash of (user bin, k-mer) & 0xffff is below this. + uint64_t unique; +}; + +struct top_metrics +{ + size_t split_tbs{}; // Technical bins holding a user bin that spans at least two technical bins. + double max_corrected{}; // Largest FPR-corrected technical bin. + std::vector> partitions{}; +}; + +// Partitions the top level and computes its metrics. +top_metrics top_level(chopper::configuration const & config, + std::vector const & positions, + std::vector const & cardinalities, + std::vector const & sketches, + std::vector const & minhashes) +{ + top_metrics m{}; + m.partitions.resize(config.hibf_config.tmax); + chopper::layout::partition_user_bins(config, positions, cardinalities, sketches, minhashes, m.partitions); + + auto const split_corr = + seqan::hibf::layout::compute_fpr_correction({.fpr = config.hibf_config.maximum_fpr, + .hash_count = config.hibf_config.number_of_hash_functions, + .t_max = config.hibf_config.tmax}); + double const relaxed = seqan::hibf::layout::compute_relaxed_fpr_correction( + {.fpr = config.hibf_config.maximum_fpr, + .relaxed_fpr = config.hibf_config.relaxed_fpr, + .hash_count = config.hibf_config.number_of_hash_functions}); + + std::map occurrences{}; + for (auto const & p : m.partitions) + if (p.size() == 1) + ++occurrences[p[0]]; + + for (auto const & p : m.partitions) + { + if (p.empty()) + continue; + if (p.size() > 1) + { + seqan::hibf::sketch::hyperloglog u{config.hibf_config.sketch_bits}; + for (size_t ub : p) + u.merge(sketches[ub]); + m.max_corrected = std::max(m.max_corrected, u.estimate() * relaxed); + } + else + { + size_t const k = occurrences[p[0]]; + if (k > 1) + ++m.split_tbs; + m.max_corrected = std::max(m.max_corrected, cardinalities[p[0]] * split_corr[k] / k); + } + } + return m; +} + +} // namespace + +int main(int argc, char ** argv) +{ + if (argc != 9) + { + std::fprintf(stderr, "Usage: %s \n", argv[0]); + return 1; + } + size_t const n = std::stoul(argv[1]); + size_t const family_divisor = std::stoul(argv[2]); + double const sigma = std::stod(argv[3]); + uint64_t const seed = std::stoul(argv[4]); + double const median = std::stod(argv[5]); + uint64_t const min_pool = std::stoul(argv[6]); + uint64_t const max_pool = std::stoul(argv[7]); + size_t const threads = std::stoul(argv[8]); + + // Input generation. Family pool sizes are log-normal; each user bin keeps 60 % to 98 % of its family's pool. + size_t const families = std::max(1, n / family_divisor); + std::mt19937_64 rng{seed * 1000003 + n * 31 + family_divisor * 7 + static_cast(sigma * 10)}; + std::lognormal_distribution size_dist{std::log(median), sigma}; + std::uniform_real_distribution sim_dist{0.6, 0.98}; + + std::vector family_pool(families); + for (auto & b : family_pool) + b = std::clamp(static_cast(size_dist(rng)), min_pool, max_pool); + + std::vector specs(n); + for (size_t u = 0; u < n; ++u) + { + uint64_t const f = rng() % families; + double const sim = sim_dist(rng); + specs[u] = {f, + family_pool[f], + static_cast(sim * 65536), + static_cast((1.0 - sim) * family_pool[f]) + 500}; + } + + chopper::configuration config{}; + config.hibf_config.number_of_user_bins = n; + config.hibf_config.threads = threads; + config.hibf_config.input_fn = [&specs, seed](size_t const u, seqan::hibf::insert_iterator it) + { + auto const & sp = specs[u]; + for (uint64_t j = 0; j < sp.pool_size; ++j) + if ((scramble(seed ^ (static_cast(u) << 32) ^ j) & 0xffff) < sp.keep_threshold) + it = scramble((sp.family << 36) + j + (seed << 60)); + for (uint64_t i = 0; i < sp.unique; ++i) + it = scramble((1ULL << 63) | (static_cast(u) << 34) | i) ^ seed; + }; + config.hibf_config.validate_and_set_defaults(); // tmax as chosen by the command line interface + + auto const t0 = std::chrono::steady_clock::now(); + std::vector sketches{}; + std::vector minhashes{}; + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minhashes); + std::vector cardinalities{}; + seqan::hibf::sketch::estimate_kmer_counts(sketches, cardinalities); + std::vector positions(n); + std::iota(positions.begin(), positions.end(), 0u); + double const t_sketch = std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); + + // Both variants on the same sketches. + top_metrics top[2]; + size_t size[2]; + double cost[2]; + double t_layout[2]; + size_t depth[2]; + for (int v = 0; v < 2; ++v) + { + chopper_variant_b = (v == 1); + top[v] = top_level(config, positions, cardinalities, sketches, minhashes); + + auto const t1 = std::chrono::steady_clock::now(); + seqan::hibf::layout::layout hibf_layout{}; + chopper::layout::fast_layout(config, positions, cardinalities, sketches, minhashes, hibf_layout); + t_layout[v] = std::chrono::duration(std::chrono::steady_clock::now() - t1).count(); + + depth[v] = 0; + for (auto const & ub : hibf_layout.user_bins) + depth[v] = std::max(depth[v], ub.previous_TB_indices.size() + 1); + + chopper::layout::hibf_statistics stats{config, sketches, cardinalities}; + stats.hibf_layout = hibf_layout; + size[v] = stats.total_hibf_size_in_byte(); + cost[v] = stats.expected_HIBF_query_cost; + } + + std::printf("%zu\t%zu\t%.1f\t%lu\t%.0f\t%lu\t%lu\t" // n families sigma seed median min max + "%zu\t%zu\t%zu\t%d\t%.0f\t%.0f\t" // tmax s_A s_B same_top topmax_A topmax_B + "%zu\t%zu\t%.4f\t%.4f\t%zu\t%zu\t" // size_A size_B cost_A cost_B depth_A depth_B + "%.1f\t%.1f\t%.1f\n", // t_sketch t_A t_B + n, + families, + sigma, + static_cast(seed), + median, + static_cast(min_pool), + static_cast(max_pool), + config.hibf_config.tmax, + top[0].split_tbs, + top[1].split_tbs, + top[0].partitions == top[1].partitions ? 1 : 0, + top[0].max_corrected, + top[1].max_corrected, + size[0], + size[1], + cost[0], + cost[1], + depth[0], + depth[1], + t_sketch, + t_layout[0], + t_layout[1]); +} diff --git a/test/benchmark/benchmark_data/post_process_clusters_variants/results_large.tsv b/test/benchmark/benchmark_data/post_process_clusters_variants/results_large.tsv new file mode 100644 index 00000000..ce59255e --- /dev/null +++ b/test/benchmark/benchmark_data/post_process_clusters_variants/results_large.tsv @@ -0,0 +1,25 @@ +n families sigma seed median min max tmax s_A s_B same_top topmax_A topmax_B size_A size_B cost_A cost_B depth_A depth_B t_sketch t_A t_B +10000 2000 1.4 1 4000 2500 2000000 128 0 0 0 420047 420047 271816094 271100603 2.4677 2.4572 3 3 0.3 3.8 3.7 +10000 2000 1.4 2 4000 2500 2000000 128 14 14 0 829415 829415 291320205 290368535 2.4083 2.4100 3 3 0.3 3.5 3.5 +10000 2000 1.4 3 4000 2500 2000000 128 2 2 0 616943 616943 307271390 307927583 2.4125 2.4125 3 3 0.3 3.6 3.6 +10000 2000 2.0 1 4000 2500 2000000 128 30 28 0 2767223 2800990 925152663 930539390 2.1987 2.2005 3 3 1.6 3.8 3.4 +10000 2000 2.0 2 4000 2500 2000000 128 40 41 0 1542129 1523446 582433437 581100821 2.1422 2.1323 3 3 1.3 12.0 8.8 +10000 2000 2.0 3 4000 2500 2000000 128 21 21 0 1800985 1818161 637861156 638600416 2.3269 2.3254 3 3 0.9 2.9 2.8 +10000 200 1.4 1 4000 2500 2000000 128 0 0 1 720738 720738 288551044 288551044 2.5977 2.5977 3 3 0.2 3.5 3.4 +10000 200 1.4 2 4000 2500 2000000 128 0 0 1 881423 881423 401792167 401792167 2.6344 2.6344 3 3 0.4 3.6 3.5 +10000 200 1.4 3 4000 2500 2000000 128 0 0 1 329586 329586 218169453 218169453 2.6613 2.6613 3 3 0.2 3.3 3.4 +10000 200 2.0 1 4000 2500 2000000 128 0 0 1 4454216 4454216 1192907945 1192907945 2.5568 2.5568 3 3 1.4 1.8 1.8 +10000 200 2.0 2 4000 2500 2000000 128 0 0 1 1090071 1090071 493768078 493768078 2.6159 2.6159 3 3 0.5 3.9 3.8 +10000 200 2.0 3 4000 2500 2000000 128 76 50 0 2383945 2230087 607717359 628659531 1.9230 2.0936 3 3 1.2 55.1 25.8 +25000 5000 1.4 1 4000 2500 2000000 192 0 0 1 720587 720587 661830984 661830984 2.5481 2.5481 3 3 0.7 14.7 14.9 +25000 5000 1.4 2 4000 2500 2000000 192 0 0 1 947761 947761 732207517 732207517 2.4958 2.4958 3 3 0.7 15.2 15.0 +25000 5000 2.0 1 4000 2500 2000000 192 0 0 1 4324752 4324752 2443841642 2443841642 2.6147 2.6147 3 3 3.3 14.1 13.9 +25000 5000 2.0 2 4000 2500 2000000 192 0 0 1 4299249 4299249 2326629222 2326629222 2.6397 2.6397 3 3 3.1 13.8 13.5 +25000 500 1.4 1 4000 2500 2000000 192 0 0 1 1205504 1205504 779610460 779610460 2.9260 2.9260 3 3 0.9 10.7 10.8 +25000 500 1.4 2 4000 2500 2000000 192 0 0 1 552034 552034 630975299 630975299 2.8742 2.8742 3 3 0.7 12.9 12.8 +25000 500 2.0 1 4000 2500 2000000 192 0 0 0 3252194 3190264 1565217124 1572320131 2.7201 2.7281 3 3 1.9 11.1 11.0 +25000 500 2.0 2 4000 2500 2000000 192 0 0 0 2370357 2144858 1725187738 1753735178 2.7782 2.8523 3 3 2.7 9.5 9.6 +50000 10000 1.4 1 4000 2500 2000000 256 0 0 1 1594383 1594383 1425899527 1425899527 2.5915 2.5915 3 3 1.5 42.1 39.6 +50000 10000 2.0 1 4000 2500 2000000 256 0 0 1 9668454 9668454 5762428676 5762428676 2.6987 2.6987 3 3 10.2 43.3 40.1 +50000 1000 1.4 1 4000 2500 2000000 256 0 0 1 488061 488061 1012672276 1012672276 3.0655 3.0655 3 3 1.3 33.0 33.3 +50000 1000 2.0 1 4000 2500 2000000 256 0 0 1 2779252 2779252 3254538846 3254538846 3.0766 3.0766 3 3 6.8 29.6 28.1 diff --git a/test/benchmark/benchmark_data/post_process_clusters_variants/results_small.tsv b/test/benchmark/benchmark_data/post_process_clusters_variants/results_small.tsv new file mode 100644 index 00000000..6bc7e6b2 --- /dev/null +++ b/test/benchmark/benchmark_data/post_process_clusters_variants/results_small.tsv @@ -0,0 +1,37 @@ +n families sigma seed median min max tmax s_A s_B same_top topmax_A topmax_B size_A size_B cost_A cost_B depth_A depth_B t_sketch t_A t_B +400 400 0.7 1 8000 5000 400000 64 0 0 1 20513 20513 19538753 19538753 2.1400 2.1400 2 2 0.0 0.1 0.1 +400 400 0.7 2 8000 5000 400000 64 0 0 0 59047 59047 12476837 12476837 1.7585 1.7585 2 2 0.0 0.1 0.1 +400 400 0.7 3 8000 5000 400000 64 0 0 1 21208 21208 20645811 20645811 2.2080 2.2080 2 2 0.0 0.1 0.1 +400 400 1.4 1 8000 5000 400000 64 19 19 0 215006 215006 24490420 24490420 1.5639 1.5639 2 2 0.0 0.1 0.0 +400 400 1.4 2 8000 5000 400000 64 23 23 0 202249 202249 24183185 24182364 1.5796 1.5795 3 3 0.0 0.0 0.1 +400 400 1.4 3 8000 5000 400000 64 11 11 0 152358 152358 19770128 19781317 1.6009 1.6025 2 2 0.0 0.1 0.1 +400 100 0.7 1 8000 5000 400000 64 4 4 0 41890 41890 11153578 11147445 1.8526 1.8471 3 3 0.0 0.1 0.1 +400 100 0.7 2 8000 5000 400000 64 8 8 0 35401 35401 12205763 12212464 1.9279 1.9299 2 2 0.0 0.1 0.1 +400 100 0.7 3 8000 5000 400000 64 2 2 0 46372 46372 11965348 12062892 1.8903 1.8886 2 2 0.0 0.1 0.1 +400 100 1.4 1 8000 5000 400000 64 30 30 0 76097 76110 14148739 14026920 1.6598 1.6121 3 3 0.0 0.0 0.0 +400 100 1.4 2 8000 5000 400000 64 32 34 0 102797 100662 17182527 16871246 1.7394 1.7251 2 2 0.0 0.0 0.0 +400 100 1.4 3 8000 5000 400000 64 0 0 1 311248 311248 39206976 39206976 1.6588 1.6588 2 2 0.0 0.0 0.0 +400 16 0.7 1 8000 5000 400000 64 24 24 0 29225 29225 9424022 9424022 1.9496 1.9496 2 2 0.0 0.0 0.0 +400 16 0.7 2 8000 5000 400000 64 4 4 0 42445 42445 13838217 13839481 1.8408 1.8404 3 3 0.0 0.0 0.0 +400 16 0.7 3 8000 5000 400000 64 0 0 1 48342 48342 11340763 11340763 1.8439 1.8439 3 3 0.0 0.0 0.0 +400 16 1.4 1 8000 5000 400000 64 0 0 1 52587 52587 9618920 9618920 1.6255 1.6255 2 2 0.0 0.0 0.0 +400 16 1.4 2 8000 5000 400000 64 0 0 1 27026 27026 34801504 34801504 6.2004 6.2004 2 2 0.0 0.1 0.1 +400 16 1.4 3 8000 5000 400000 64 34 34 0 53198 53198 12780279 12780279 1.8130 1.8130 2 2 0.0 0.0 0.0 +1000 1000 0.7 1 8000 5000 400000 64 0 0 1 51153 51153 36289950 36289950 2.1063 2.1063 2 2 0.0 0.2 0.2 +1000 1000 0.7 2 8000 5000 400000 64 0 0 1 74764 74764 34098151 34098151 2.0457 2.0457 3 3 0.0 0.2 0.2 +1000 1000 0.7 3 8000 5000 400000 64 0 0 1 58397 58397 32223575 32223575 2.0631 2.0631 3 3 0.0 0.2 0.2 +1000 1000 1.4 1 8000 5000 400000 64 2 2 0 323714 323714 60208762 60173993 1.8570 1.8520 3 3 0.1 0.2 0.2 +1000 1000 1.4 2 8000 5000 400000 64 0 0 0 275146 275146 43660682 43722571 1.7981 1.7947 2 3 0.1 0.2 0.2 +1000 1000 1.4 3 8000 5000 400000 64 4 4 0 355200 355200 57518117 57550674 1.7281 1.7298 3 3 0.1 0.2 0.2 +1000 250 0.7 1 8000 5000 400000 64 0 0 1 75316 75316 28970614 28970614 2.1098 2.1098 3 3 0.0 0.2 0.2 +1000 250 0.7 2 8000 5000 400000 64 0 0 1 38631 38631 27975704 27975704 2.1763 2.1763 3 3 0.0 0.2 0.2 +1000 250 0.7 3 8000 5000 400000 64 0 0 1 55895 55895 31672322 31672322 2.1554 2.1554 3 3 0.0 0.2 0.2 +1000 250 1.4 1 8000 5000 400000 64 31 29 0 299699 308854 45709105 45474949 1.7036 1.6895 3 3 0.1 0.1 0.1 +1000 250 1.4 2 8000 5000 400000 64 4 4 0 211132 211132 52492676 52508986 1.9482 1.9661 3 3 0.1 0.1 0.1 +1000 250 1.4 3 8000 5000 400000 64 17 17 0 244577 244577 43940607 44074690 1.7999 1.8234 3 3 0.1 0.1 0.1 +1000 40 0.7 1 8000 5000 400000 64 0 0 1 50187 50187 35395618 35395618 2.4898 2.4898 3 3 0.1 0.2 0.2 +1000 40 0.7 2 8000 5000 400000 64 14 14 0 82291 82291 26686518 26685507 1.9903 1.9897 2 2 0.1 0.1 0.1 +1000 40 0.7 3 8000 5000 400000 64 14 14 0 74042 74042 22976814 23038071 1.9091 1.8903 3 3 0.1 0.1 0.1 +1000 40 1.4 1 8000 5000 400000 64 12 12 0 464279 464279 48581499 48908774 1.5701 1.5939 3 3 0.1 0.1 0.1 +1000 40 1.4 2 8000 5000 400000 64 0 0 1 366254 366254 94906654 94906654 2.0779 2.0779 2 2 0.1 0.1 0.1 +1000 40 1.4 3 8000 5000 400000 64 0 0 1 104841 104841 68497595 68497595 3.0255 3.0255 3 3 0.1 0.1 0.1 diff --git a/test/benchmark/benchmark_data/post_process_clusters_variants/run.sh b/test/benchmark/benchmark_data/post_process_clusters_variants/run.sh new file mode 100755 index 00000000..0d63ad10 --- /dev/null +++ b/test/benchmark/benchmark_data/post_process_clusters_variants/run.sh @@ -0,0 +1,39 @@ +#!/usr/bin/env bash +# Runs both input grids of README.md and writes results_small.tsv and results_large.tsv to the current directory. +# +# Usage: ./run.sh + +set -Eeuo pipefail + +if [[ $# -ne 1 ]]; then + echo "Usage: $0 " >&2 + exit 1 +fi + +MEASURE=$1 +HEADER="n\tfamilies\tsigma\tseed\tmedian\tmin\tmax\ttmax\ts_A\ts_B\tsame_top\ttopmax_A\ttopmax_B\tsize_A\tsize_B\tcost_A\tcost_B\tdepth_A\tdepth_B\tt_sketch\tt_A\tt_B" + +# Small grid: 400 and 1000 user bins, median 8000 k-mers, 3 seeds each. +printf '%b\n' "${HEADER}" > results_small.tsv +for n in 400 1000; do + for family_divisor in 1 4 25; do + for sigma in 0.7 1.4; do + for seed in 1 2 3; do + "${MEASURE}" "${n}" "${family_divisor}" "${sigma}" "${seed}" 8000 5000 400000 8 >> results_small.tsv + done + done + done +done + +# Large grid: 10000, 25000 and 50000 small user bins, median 4000 k-mers, 3, 2 and 1 seeds. +printf '%b\n' "${HEADER}" > results_large.tsv +for spec in "10000 3" "25000 2" "50000 1"; do + read -r n seeds <<< "${spec}" + for family_divisor in 5 50; do + for sigma in 1.4 2.0; do + for seed in $(seq 1 "${seeds}"); do + "${MEASURE}" "${n}" "${family_divisor}" "${sigma}" "${seed}" 4000 2500 2000000 32 >> results_large.tsv + done + done + done +done diff --git a/test/benchmark/benchmark_data/post_process_clusters_variants/variant_b.patch b/test/benchmark/benchmark_data/post_process_clusters_variants/variant_b.patch new file mode 100644 index 00000000..cc016c72 --- /dev/null +++ b/test/benchmark/benchmark_data/post_process_clusters_variants/variant_b.patch @@ -0,0 +1,52 @@ +--- a/src/layout/partition_user_bins.cpp ++++ b/src/layout/partition_user_bins.cpp +@@ -31,6 +31,9 @@ + #include + #include + ++// Measurement switch: false = variant A (as in chopper), true = variant B (seed count = number_of_remaining_tbs). ++bool chopper_variant_b = false; ++ + namespace chopper::layout + { + +@@ -238,7 +241,8 @@ + */ + void post_process_clusters(std::vector & clusters, + std::vector const & cardinalities, +- chopper::configuration const & config) ++ [[maybe_unused]] chopper::configuration const & config, ++ size_t const number_of_seeds) + { + // clusters are done. Start post processing + // since post processing involves re-ordering the clusters, the moved_to_cluster_id value of a cluster will not +@@ -258,7 +262,7 @@ + + // push largest p clusters to the front + std::ranges::partial_sort(clusters, +- std::ranges::next(clusters.begin(), config.hibf_config.tmax, clusters.end()), ++ std::ranges::next(clusters.begin(), number_of_seeds, clusters.end()), + std::ranges::greater{}, + [&largest_user_bin_cardinality](Cluster const & c) + { +@@ -270,7 +274,7 @@ + // the largest ub is already at the start because of former sorting. + // Empty clusters are sorted last explicitly. Their cardinality key 0 does not suffice, because a non-empty cluster + // can have a cardinality estimate of 0, too. +- std::ranges::sort(clusters | std::views::drop(config.hibf_config.tmax), ++ std::ranges::sort(clusters | std::views::drop(number_of_seeds), + std::ranges::greater{}, + [&largest_user_bin_cardinality](Cluster const & c) + { +@@ -480,7 +484,10 @@ + sketches, + technical_bin_size_threshold, + config); +- post_process_clusters(clusters, cardinalities, config); ++ post_process_clusters(clusters, ++ cardinalities, ++ config, ++ ::chopper_variant_b ? number_of_remaining_tbs : config.hibf_config.tmax); + lsh_algorithm_timer.stop(); + config.lsh_algorithm_timer += lsh_algorithm_timer; +