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..8bb2e653 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,23 @@ 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); diff --git a/include/chopper/layout/determine_split_bins.hpp b/include/chopper/layout/determine_split_bins.hpp new file mode 100644 index 00000000..d0b173d6 --- /dev/null +++ b/include/chopper/layout/determine_split_bins.hpp @@ -0,0 +1,53 @@ +// -------------------------------------------------------------------------------------------------- +// 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::determine_split_bins. + * \author Svenja Mehringer + */ + +#pragma once + +#include +#include +#include + +#include + +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, + size_t const num_technical_bins, + size_t const num_user_bins, + std::vector> & partitions); + +} // namespace chopper::layout diff --git a/include/chopper/layout/execute.hpp b/include/chopper/layout/execute.hpp index 52bfc3b7..f3e4c161 100644 --- a/include/chopper/layout/execute.hpp +++ b/include/chopper/layout/execute.hpp @@ -13,12 +13,32 @@ #include #include +#include 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); + 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..b40c885f --- /dev/null +++ b/include/chopper/layout/fast_layout_cluster.hpp @@ -0,0 +1,187 @@ +// --------------------------------------------------------------------------------------------------- +// 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 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 +{ +private: + //!\brief The id of the cluster. + size_t representative_id{}; + + //!\brief The user bins contained in this cluster. + std::vector user_bins{}; + + //!\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: + /*!\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(); + user_bins.pop_back(); + 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; + bool const properly_moved = has_been_moved() && empty(); + bool const not_moved = !has_been_moved() && !empty(); + + 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()); + assert(is_valid(representative_id)); + 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; + auto & source = this->user_bins; + target.insert(target.end(), source.cbegin(), source.cend()); + source = std::vector{}; // .clear() AND release memory + + 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, + std::ranges::greater{}, + [&cardinalities](size_t const i) + { + return cardinalities[i]; + }); + } +}; + +/*!\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]; + + 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 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..04ba0d22 --- /dev/null +++ b/include/chopper/layout/fast_layout_find_bins_to_be_split.hpp @@ -0,0 +1,90 @@ +// -------------------------------------------------------------------------------------------------- +// 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 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 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, + 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/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/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/chopper_layout.cpp b/src/chopper_layout.cpp index 5fda4ddb..28027989 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"); @@ -93,6 +96,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,7 +110,12 @@ 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); + + 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 { @@ -128,17 +137,23 @@ 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); + // 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(); } - 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 +166,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" + << "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'; 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.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'; } 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..c4e85f3d --- /dev/null +++ b/src/layout/determine_split_bins.cpp @@ -0,0 +1,148 @@ +// --------------------------------------------------------------------------------------------------- +// 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 +{ + +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]); + + minimum = (score < minimum) ? (trace[i][j] = i_prime, score) : minimum; + } + + matrix[i][j] = minimum; + } + } + + // 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_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_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_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..d18be871 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,50 @@ 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) { + 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; + 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); + 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..1342fe85 --- /dev/null +++ b/src/layout/fast_layout.cpp @@ -0,0 +1,382 @@ +// --------------------------------------------------------------------------------------------------- +// 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 + +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::concurrent_timer union_estimation_timer{}; + seqan::hibf::concurrent_timer rearrangement_timer{}; + + 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. + * \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 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, 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, 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. + * - **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. + * + * Not thread-safe; fast_layout_recursion serialises it with `omp critical`. + */ +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}; + + 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) + { + 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]); + } + + // 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(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 + + 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; + } + } + } + + return max_bin_id; +} + +/*!\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 + { + 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) + { + 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.initial_partition_timer.start(); + partition_user_bins(config, positions, cardinalities, sketches, minHash_sketches, tmax_partitions); + config.initial_partition_timer.stop(); + + // initialise user bins in layout + 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) +#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, + std::ranges::less{}, + [](auto const & mb) + { + // 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 + +#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/hibf_statistics.cpp b/src/layout/hibf_statistics.cpp index df86775c..ca9a9679 100644 --- a/src/layout/hibf_statistics.cpp +++ b/src/layout/hibf_statistics.cpp @@ -125,8 +125,10 @@ 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"; std::stringstream ss; ss << std::fixed << std::setprecision(2) << num; return ss.str(); diff --git a/src/layout/partition_user_bins.cpp b/src/layout/partition_user_bins.cpp new file mode 100644 index 00000000..552f64ad --- /dev/null +++ b/src/layout/partition_user_bins.cpp @@ -0,0 +1,739 @@ +// -------------------------------------------------------------------------------------------------- +// 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 +#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 [first, last] = std::ranges::unique(list); + list.erase(first, last); + } + + 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 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); + + 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`). + * + * 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); + } + + // The user bins are sorted by cardinality, so the first one is the largest. + 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()), + std::ranges::greater{}, + [&largest_user_bin_cardinality](Cluster const & c) + { + 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. + // 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](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 + 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. + * \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. 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. + * \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] 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. + * + * 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. + */ +void 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> & partitions, + std::vector & partition_sketches, + 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}; + + 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}; + + auto penalty_lower_level = [&](size_t const additional_number_of_user_bins, size_t const p) -> size_t + { + 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 (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{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 (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 + // 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 // 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 + 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(); + + // 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); + 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; + } + } + + // now that we know which partition fits best (`best_p`), add those indices to it + for (size_t const user_bin_idx : cluster) + { + 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]); + } + partition_sketches[best_p].merge(current_sketch); + corrected_estimate_per_part = std::max(corrected_estimate_per_part, partition_sketches[best_p].estimate()); +} + +/*!\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()); + + // 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) + 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); + 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, + [](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]); + } + + // 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; + } + + for (size_t i = cidx; i < clusters.size(); ++i) + { + if (clusters[i].empty()) + break; + + remaining_clusters.push_back(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]; + + 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); + search_partition_algorithm_timer.stop(); + } + config.search_partition_algorithm_timer += search_partition_algorithm_timer; + + // 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/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/CMakeLists.txt b/test/api/layout/CMakeLists.txt index 9cdb4076..907851aa 100644 --- a/test/api/layout/CMakeLists.txt +++ b/test/api/layout/CMakeLists.txt @@ -26,3 +26,8 @@ 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 (fast_layout_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..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 @@ -53,9 +54,10 @@ 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)); @@ -83,6 +85,98 @@ 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{}; + 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)); + + 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"); +} +#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) { seqan3::test::tmp_directory tmp_dir{}; @@ -114,9 +208,10 @@ 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); + chopper::layout::execute(config, many_filenames, sketches, minHash_sketches); ASSERT_TRUE(std::filesystem::exists(stats_file)); @@ -399,6 +494,363 @@ 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{}; + 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, 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; +} +#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) { seqan3::test::tmp_directory tmp_dir{}; @@ -429,9 +881,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 +972,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/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); +} diff --git a/test/api/layout/hibf_statistics_test.cpp b/test/api/layout/hibf_statistics_test.cpp index 26c1dcb6..92bf9745 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..11a81b99 --- /dev/null +++ b/test/api/layout/partition_user_bins_test.cpp @@ -0,0 +1,236 @@ +// -------------------------------------------------------------------------------------------------- +// 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 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 +#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. +// 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 const & zero_cardinality_user_bins = {}) +{ + if (content_ids.empty()) + { +#if CHOPPER_WORKAROUND_GCC_BOGUS_ARRAY +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Warray-bounds=" +#endif // CHOPPER_WORKAROUND_GCC_BOGUS_ARRAY + content_ids.resize(kmer_counts.size()); +#if CHOPPER_WORKAROUND_GCC_BOGUS_ARRAY +# pragma GCC diagnostic pop +#endif // CHOPPER_WORKAROUND_GCC_BOGUS_ARRAY + 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(); + 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); + + std::vector> partitions(tmax); + chopper::layout::partition_user_bins(config, positions, cardinalities, sketches, minHash_sketches, partitions); + + 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. +// 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); +#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) +{ + // 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); +} + +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); + + expect_all_user_bins_assigned(partitions, kmer_counts.size()); +} + +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()); +} 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; + diff --git a/test/cli/cli_chopper_basic_test.cpp b/test/cli/cli_chopper_basic_test.cpp index 80d31251..3f0093d7 100644 --- a/test/cli/cli_chopper_basic_test.cpp +++ b/test/cli/cli_chopper_basic_test.cpp @@ -94,3 +94,62 @@ 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; +} + +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"}); +} 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"}); } 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..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(), 0); // 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 +}