diff --git a/include/chopper/configuration.hpp b/include/chopper/configuration.hpp index 34ed74aa..5371251a 100644 --- a/include/chopper/configuration.hpp +++ b/include/chopper/configuration.hpp @@ -10,6 +10,9 @@ #include #include #include +#include +#include +#include #include @@ -20,6 +23,19 @@ namespace chopper { +//!\brief The partitioning approaches of the partitioned HIBF. Selected by `configuration::partitioning_approach`. +enum struct partitioning_scheme : uint8_t +{ + blocked, // 0 + sorted, // 1 + folded, // 2 + weighted_fold, // 3 + similarity, // 4 + lsh, // 5 + lsh_sim, // 6 + unknown // for error parsing +}; + struct configuration { //!\brief Whether to use the fast layout algorithm instead of the default one. @@ -50,6 +66,23 @@ struct configuration bool precomputed_files{false}; //!\} + /*!\name Partitioned HIBF configuration + * \{ + */ + //!\brief The number of partitions for the HIBF index. 1 computes a single HIBF layout. + size_t number_of_partitions{1}; + + //!\brief The partitioning approach. See chopper::partitioning_scheme. + partitioning_scheme partitioning_approach{}; + + /*!\brief Whether `hibf_config.tmax` was given by the user. Not serialised. + * + * If set, every partition of a partitioned HIBF is laid out with `hibf_config.tmax`. Otherwise, the tmax of each + * partition is chosen based on its number of user bins. + */ + bool tmax_is_set{false}; + //!\} + /*!\name Configuration of size estimates * \{ */ @@ -111,7 +144,8 @@ struct configuration void serialize(archive_t & archive) { // Version 3 added fast_layout. - uint32_t version{3}; + // Version 4 added the partitioned HIBF configuration. + uint32_t version{4}; archive(CEREAL_NVP(version)); archive(CEREAL_NVP(data_file)); @@ -126,10 +160,69 @@ struct configuration archive(CEREAL_NVP(determine_best_tmax)); archive(CEREAL_NVP(force_all_binnings)); - // Files written before version 3 do not contain this field. Reading it unconditionally would throw. if (version >= 3) archive(CEREAL_NVP(fast_layout)); + + if (version >= 4) + { + archive(CEREAL_NVP(number_of_partitions)); + archive(CEREAL_NVP(partitioning_approach)); + } } }; } // namespace chopper + +namespace chopper +{ + +// Make chopper::partitioning_scheme satisfy sharg::ostreamable. Found by argument-dependent lookup. +inline std::ostream & operator<<(std::ostream & output, chopper::partitioning_scheme const & scheme) +{ + if (scheme == chopper::partitioning_scheme::blocked) + output << "blocked"; + else if (scheme == chopper::partitioning_scheme::sorted) + output << "sorted"; + else if (scheme == chopper::partitioning_scheme::folded) + output << "folded"; + else if (scheme == chopper::partitioning_scheme::weighted_fold) + output << "weighted_fold"; + else if (scheme == chopper::partitioning_scheme::similarity) + output << "similarity"; + else if (scheme == chopper::partitioning_scheme::lsh) + output << "lsh"; + else if (scheme == chopper::partitioning_scheme::lsh_sim) + output << "lsh_sim"; + else + output << "unknown"; + + return output; +} + +// Make chopper::partitioning_scheme satisfy sharg::istreamable. Found by argument-dependent lookup. +inline std::istream & operator>>(std::istream & input, chopper::partitioning_scheme & scheme) +{ + std::string text; + input >> text; + + if (text == "blocked") + scheme = chopper::partitioning_scheme::blocked; + else if (text == "sorted") + scheme = chopper::partitioning_scheme::sorted; + else if (text == "folded") + scheme = chopper::partitioning_scheme::folded; + else if (text == "weighted_fold") + scheme = chopper::partitioning_scheme::weighted_fold; + else if (text == "similarity") + scheme = chopper::partitioning_scheme::similarity; + else if (text == "lsh") + scheme = chopper::partitioning_scheme::lsh; + else if (text == "lsh_sim") + scheme = chopper::partitioning_scheme::lsh_sim; + else + scheme = chopper::partitioning_scheme::unknown; + + return input; +} + +} // namespace chopper diff --git a/include/chopper/layout/execute.hpp b/include/chopper/layout/execute.hpp index c4d580c9..6376960b 100644 --- a/include/chopper/layout/execute.hpp +++ b/include/chopper/layout/execute.hpp @@ -24,12 +24,16 @@ namespace chopper::layout * \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. + * bin, if `config.fast_layout` is set or a partitioned HIBF with the `lsh` or `lsh_sim` + * partitioning approach is computed. May be empty otherwise. * \returns 0. - * \throws std::invalid_argument If both `config.determine_best_tmax` and `config.fast_layout` are set, or if the - * computed layout is invalid (see seqan::hibf::layout::layout::validate). + * \throws std::invalid_argument If... + * * ...both `config.determine_best_tmax` and `config.fast_layout` are set. + * * ...`config.number_of_partitions >= 2` and either of them is set. + * * ...the computed layout is invalid (see seqan::hibf::layout::layout::validate). * - * The layout is computed with + * If `config.number_of_partitions >= 2`, a partitioned HIBF layout is computed with partition_user_bins. + * Otherwise, 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. diff --git a/include/chopper/layout/fast_layout_cluster.hpp b/include/chopper/layout/fast_layout_cluster.hpp index b40c885f..614bc1f6 100644 --- a/include/chopper/layout/fast_layout_cluster.hpp +++ b/include/chopper/layout/fast_layout_cluster.hpp @@ -131,6 +131,20 @@ struct Cluster return moved_id.value(); } + /*!\brief Moves end range of user bins of this cluster to `target_cluster`. + * \param[in,out] target_cluster The cluster to move the user bins to. Must be a different cluster. + * \param[in] size The size of the part that should be moved, e.g. the number of user bins. + * + * No ids are changed. This happens in post_process_clusters. + */ + void move_part_to(Cluster & target_cluster, size_t size) + { + auto & target = target_cluster.user_bins; + auto & source = this->user_bins; + target.insert(target.end(), source.cbegin() + size, source.cend()); + source.erase(source.begin() + size, source.end()); + } + /*!\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. * @@ -161,6 +175,8 @@ struct Cluster }; /*!\brief Returns the position of the representative cluster of `clusters[current_id]`. + * \tparam cluster_t The cluster type. Cluster or a type with the same interface, e.g., the MultiCluster of the + * partitioned HIBF. * \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 @@ -168,9 +184,10 @@ struct Cluster * * 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) +template +size_t LSH_find_representative_cluster(std::vector const & clusters, size_t current_id) { - std::reference_wrapper representative = clusters[current_id]; + std::reference_wrapper representative = clusters[current_id]; assert(representative.get().is_valid(current_id)); diff --git a/include/chopper/layout/input.hpp b/include/chopper/layout/input.hpp index f99d07ab..47900df5 100644 --- a/include/chopper/layout/input.hpp +++ b/include/chopper/layout/input.hpp @@ -21,12 +21,12 @@ namespace chopper::layout std::vector> read_filenames_from(std::istream & stream); -/*!\brief Reads a layout file: the file names of the user bins, the chopper configuration and the layout. - * \param[in] stream The content of a layout file. - * \returns The file names, the configuration and the layout. - * \details - * The layout is not validated. Validate it with seqan::hibf::layout::layout::validate and the returned - * `configuration::hibf_config` before using it. +/*!\brief Reads a layout file that contains one or more (partitioned HIBF) layouts. +* \param[in] stream The content of a layout file. +* \returns The filenames, the chopper configuration, and the layouts in the order they appear in the file. +* \details +* The layout is not validated. Validate it with seqan::hibf::layout::layout::validate and the returned +* `configuration::hibf_config` before using it. */ std::tuple>, configuration, seqan::hibf::layout::layout> read_layout_file(std::istream & stream); diff --git a/include/chopper/layout/lsh.hpp b/include/chopper/layout/lsh.hpp new file mode 100644 index 00000000..972feb20 --- /dev/null +++ b/include/chopper/layout/lsh.hpp @@ -0,0 +1,112 @@ +// -------------------------------------------------------------------------------------------------- +// 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 the MinHash LSH primitives shared by the fast layout and the partitioned HIBF. + * \author Svenja Mehringer + */ + +#pragma once + +#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. + */ +inline 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. + * \tparam cluster_t The cluster type. Either chopper::layout::Cluster, whose + * `contained_user_bins()` is a range of user bins, or a type whose + * `contained_user_bins()` is a range of ranges of user bins (e.g., the + * MultiCluster of the partitioned HIBF). + * \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. + */ +template +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; + + auto insert = [&](size_t const user_bin_idx) + { + ++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 + }; + + using user_bins_t = std::remove_cvref_t; + + if constexpr (std::ranges::range>) // e.g. MultiCluster + { + for (auto const & similarity_cluster : current.contained_user_bins()) + for (size_t const user_bin_idx : similarity_cluster) + insert(user_bin_idx); + } + else + { + for (size_t const user_bin_idx : current.contained_user_bins()) + insert(user_bin_idx); + } + } + 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; +} + +} // namespace chopper::layout diff --git a/include/chopper/layout/multi_cluster.hpp b/include/chopper/layout/multi_cluster.hpp new file mode 100644 index 00000000..d88992d3 --- /dev/null +++ b/include/chopper/layout/multi_cluster.hpp @@ -0,0 +1,151 @@ +// --------------------------------------------------------------------------------------------------- +// 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::MultiCluster. + * \author Svenja Mehringer + */ + +#pragma once + +#include +#include +#include +#include +#include + +#include + +namespace chopper::layout +{ + +/*!\brief A cluster of clusters of user bins, used by the `lsh_sim` partitioning of the partitioned HIBF. + * + * A MultiCluster is created from a Cluster (the result of a first, very similar, LSH clustering). The second LSH + * clustering then joins MultiClusters, such that a MultiCluster holds several clusters of very similar user bins. + * + * The validity invariants are the same as for chopper::layout::Cluster, with `contained_user_bins()` being a vector of + * clusters instead of a vector of user bins. Hence, chopper::layout::LSH_find_representative_cluster and + * chopper::layout::LSH_fill_hashtable work on MultiClusters, too. + */ +struct MultiCluster +{ +private: + //!\brief The id of the cluster. + size_t representative_id{}; + + //!\brief The clusters of 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 + * \{ + */ + MultiCluster() = default; //!< Defaulted. + MultiCluster(MultiCluster const &) = default; //!< Defaulted. + MultiCluster(MultiCluster &&) = default; //!< Defaulted. + MultiCluster & operator=(MultiCluster const &) = default; //!< Defaulted. + MultiCluster & operator=(MultiCluster &&) = default; //!< Defaulted. + ~MultiCluster() = default; //!< Defaulted. + + /*!\brief Creates a MultiCluster from a Cluster. Not explicit, such that `{Cluster{0}}` converts. + * \param[in] clust The cluster. If it has been moved, the MultiCluster is moved to the same id and is empty. + * Otherwise, the MultiCluster contains the user bins of `clust` as its single cluster. + */ + MultiCluster(Cluster const & clust) : representative_id{clust.id()} + { + if (clust.has_been_moved()) + moved_id = clust.moved_to_cluster_id(); + else + user_bins.push_back(clust.contained_user_bins()); + } + //!\} + + //!\brief Returns the id of the cluster. + size_t id() const + { + return representative_id; + } + + //!\brief Returns the clusters of 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 clusters. + bool empty() const + { + return user_bins.empty(); + } + + //!\brief Returns the number of clusters in the cluster. + size_t size() const + { + return user_bins.size(); + } + + /*!\brief Checks that the cluster at position `id` is either valid or moved, see chopper::layout::Cluster. + * \param[in] id The position of the cluster, i.e., the id it must have. + */ + 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 clusters of this cluster to `target_cluster` and marks this cluster as moved. + * \param[in,out] target_cluster The cluster to move the clusters to. Must be a different cluster. + */ + void move_to(MultiCluster & target_cluster) + { + target_cluster.user_bins.insert(target_cluster.user_bins.end(), this->user_bins.begin(), this->user_bins.end()); + this->user_bins.clear(); + moved_id = target_cluster.id(); + } + + /*!\brief Sorts the user bins within each cluster by descending cardinality, and the clusters by descending size. + * \param[in] cardinalities The cardinality of each user bin, indexed by user bin. + */ + void sort_by_cardinality(std::vector const & cardinalities) + { + auto cmp = [&cardinalities](auto const & v1, auto const & v2) + { + return cardinalities[v2] < cardinalities[v1]; + }; + for (auto & user_bin_cluster : user_bins) + std::sort(user_bin_cluster.begin(), user_bin_cluster.end(), cmp); + + auto cmp_clusters = [](auto const & c1, auto const & c2) + { + return c2.size() < c1.size(); + }; + std::sort(user_bins.begin(), user_bins.end(), cmp_clusters); + } +}; + +} // namespace chopper::layout diff --git a/include/chopper/layout/output.hpp b/include/chopper/layout/output.hpp index ac183769..0dd8721d 100644 --- a/include/chopper/layout/output.hpp +++ b/include/chopper/layout/output.hpp @@ -14,6 +14,8 @@ namespace chopper::layout { -void write_user_bins_to(std::vector> const & filenames, std::ostream & stream); +void write_user_bins_to(std::vector> const & filenames, + std::vector const & mapping, + std::ostream & stream); } // 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..c7883ec3 --- /dev/null +++ b/include/chopper/layout/partition_user_bins.hpp @@ -0,0 +1,48 @@ +// --------------------------------------------------------------------------------------------------- +// 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 all user bins onto `config.number_of_partitions` partitions of a partitioned HIBF. + * \param[in] config The configuration. Uses `number_of_partitions`, `partitioning_approach`, and + * `hibf_config.{sketch_bits, number_of_user_bins, tmax}`. + * \param[in] cardinalities The cardinality of each user bin. + * \param[in] sketches The HyperLogLog sketch of each user bin. + * \param[in] minHash_sketches The MinHash tables of each user bin. Only used by the `lsh` and `lsh_sim` approaches. + * \param[out] partitions Must have size `config.number_of_partitions` on entry. On return, `partitions[i]` holds + * the user bin indices assigned to partition `i`. No partition is empty. + * \throws std::invalid_argument If `config.partitioning_approach` is not a chopper::partitioning_scheme, or if there are + * fewer user bins than partitions. + * \throws std::logic_error If not all user bins have been assigned to a partition. + * \throws std::runtime_error If a partition is empty. + * + * Each partition is laid out as a separate HIBF. In contrast to the fast layout's + * chopper::layout::partition_user_bins, user bins are never split across partitions. + */ +void partition_user_bins(chopper::configuration const & config, + std::vector const & cardinalities, + std::vector const & sketches, + std::vector const & minHash_sketches, + std::vector> & partitions); + +} // namespace chopper::layout diff --git a/src/chopper_layout.cpp b/src/chopper_layout.cpp index f675f8ec..9585c8e2 100644 --- a/src/chopper_layout.cpp +++ b/src/chopper_layout.cpp @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -81,6 +82,15 @@ int chopper_layout(chopper::configuration & config, sharg::parser & parser) if (config.fast_layout && config.determine_best_tmax) throw sharg::parser_error{"You cannot combine --fast-layout with --determine-best-tmax."}; + bool const partitioned = config.number_of_partitions >= 2u; + config.tmax_is_set = parser.is_option_set("tmax"); + + if (partitioned && config.fast_layout) + throw sharg::parser_error{"You cannot combine --fast-layout with --number-of-partitions."}; + + if (partitioned && config.determine_best_tmax) + throw sharg::parser_error{"You cannot combine --determine-best-tmax with --number-of-partitions."}; + auto has_sketch_file_extension = [](std::filesystem::path const & path) { return path.string().ends_with(".sketch") || path.string().ends_with(".sketches"); @@ -116,6 +126,14 @@ int chopper_layout(chopper::configuration & config, sharg::parser & parser) 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."}; + + bool const partitioning_needs_minhashes = config.partitioning_approach == chopper::partitioning_scheme::lsh + || config.partitioning_approach == chopper::partitioning_scheme::lsh_sim; + + if (partitioned && partitioning_needs_minhashes && minHash_sketches.size() != sketches.size()) + throw sharg::parser_error{"The sketch file does not contain MinHash sketches, which the chosen " + "--partitioning-approach needs. Create the sketch file with " + "--number-of-partitions."}; } else { @@ -129,6 +147,13 @@ int chopper_layout(chopper::configuration & config, sharg::parser & parser) chopper::sketch::check_filenames(filenames, config); } + if (partitioned && filenames.size() < config.number_of_partitions) + throw sharg::parser_error{sharg::detail::to_string("The number of partitions (", + config.number_of_partitions, + ") must not exceed the number of user bins (", + filenames.size(), + ").")}; + config.hibf_config.input_fn = chopper::input_functor{filenames, config.precomputed_files, config.k, config.window_size}; config.hibf_config.number_of_user_bins = filenames.size(); @@ -137,9 +162,9 @@ int chopper_layout(chopper::configuration & config, sharg::parser & parser) if (!input_is_a_sketch_file) { config.compute_sketches_timer.start(); - // 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) + // Only the fast layout and the partitioned HIBF need 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 || partitioned) seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); else seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches); diff --git a/src/layout/CMakeLists.txt b/src/layout/CMakeLists.txt index 6b9d3b47..9007d8cf 100644 --- a/src/layout/CMakeLists.txt +++ b/src/layout/CMakeLists.txt @@ -14,6 +14,7 @@ add_library (chopper_layout STATIC input.cpp lsh_distributed_ibf_layout.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/execute.cpp b/src/layout/execute.cpp index ea52985a..985a053f 100644 --- a/src/layout/execute.cpp +++ b/src/layout/execute.cpp @@ -20,6 +20,8 @@ #include #include #include +#include +#include #include #include @@ -34,62 +36,113 @@ int execute(chopper::configuration & config, std::vector const & sketches, std::vector const & minHash_sketches) { + bool const partitioned = config.number_of_partitions >= 2u; // 0 == unset == single HIBF, 1 == single HIBF + assert(config.number_of_partitions >= 1u); // 0 == fail, 1 == no partitioning, => 2 PHIBF + if (config.determine_best_tmax && config.fast_layout) throw std::invalid_argument{"determine_best_tmax is not supported with fast_layout."}; + if (partitioned && (config.determine_best_tmax)) + throw std::invalid_argument{"The partitioned HIBF (number_of_partitions >= 2) is not supported with " + "determine_best_tmax."}; + config.hibf_config.validate_and_set_defaults(); std::vector cardinalities; seqan::hibf::sketch::estimate_kmer_counts(sketches, cardinalities); + std::vector> partitions(config.number_of_partitions); + std::vector hibf_layouts(config.number_of_partitions); - seqan::hibf::layout::layout hibf_layout; + if (partitioned) + partition_user_bins(config, cardinalities, sketches, minHash_sketches, partitions); + else + partitions[0] = seqan::hibf::iota_vector(sketches.size()); if (config.determine_best_tmax) { - hibf_layout = determine_best_number_of_technical_bins(config, cardinalities, sketches); + hibf_layouts[0] = determine_best_number_of_technical_bins(config, cardinalities, sketches); } else { - config.dp_algorithm_timer.start(); - if (config.fast_layout) - { - fast_layout(config, - seqan::hibf::iota_vector(sketches.size()), - cardinalities, - sketches, - minHash_sketches, - hibf_layout); - } - else +#pragma omp parallel for schedule(dynamic) num_threads(config.hibf_config.threads) + for (size_t i = 0; i < config.number_of_partitions; ++i) { - 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); + // reset tmax to sqrt(n) where n is the number of user bins in the partition + auto config_copy = config; + + if (!config.tmax_is_set) + { + config_copy.hibf_config.tmax = + chopper::next_multiple_of_64(static_cast(std::ceil(std::sqrt(partitions[i].size())))); + } + + // use timers of config_copy for thread safe start() and stop() and add result to config afterwards + config_copy.dp_algorithm_timer.start(); + if (config.fast_layout) + { + fast_layout(config_copy, + std::move(partitions[i]), + cardinalities, + sketches, + minHash_sketches, + hibf_layouts[i]); + } + else + { + hibf_layouts[i] = seqan::hibf::layout::compute_layout(config_copy.hibf_config, + cardinalities, + sketches, + std::move(partitions[i]), + config_copy.union_estimation_timer, + config_copy.rearrangement_timer); + config.union_estimation_timer += config_copy.union_estimation_timer; + config.rearrangement_timer += config_copy.rearrangement_timer; + } + config_copy.dp_algorithm_timer.stop(); + + config.dp_algorithm_timer += config_copy.dp_algorithm_timer; } - config.dp_algorithm_timer.stop(); } - hibf_layout.validate(config.hibf_config); + // local and global indices + std::vector> mapping_per_partition(config.number_of_partitions); - // With determine_best_tmax, determine_best_number_of_technical_bins already wrote the statistics of every tmax to - // `.stats`. - if (!config.determine_best_tmax && config.output_verbose_statistics) + for (size_t i = 0; i < config.number_of_partitions; ++i) { - size_t dummy{}; - 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); - } + auto & layout = hibf_layouts[i]; + std::vector mapping = seqan::hibf::iota_vector(layout.user_bins.size()); - // brief Write the output to the layout file. - std::ofstream fout{config.output_filename}; - chopper::layout::write_user_bins_to(filenames, fout); - config.write_to(fout); - hibf_layout.write_to(fout); + if (config.number_of_partitions > 1) + { + // update global indicies to local ones within the partitionied lauyout + mapping.resize(layout.user_bins.size()); + for (size_t ubx{0}; ubx < layout.user_bins.size(); ++ubx) + { + mapping[ubx] = layout.user_bins[ubx].idx; // for current local store its global UB.idx + layout.user_bins[ubx].idx = ubx; // overwrite global index with local one + } + + // update hibf_config.number_of_user_bins + config.hibf_config.number_of_user_bins = hibf_layouts[i].user_bins.size(); + } + + hibf_layouts[i].validate(config.hibf_config); + + // With determine_best_tmax, statistics of every tmax have already been written to `.stats` + if (!config.determine_best_tmax && config.output_verbose_statistics) + { + size_t dummy{}; + chopper::layout::hibf_statistics global_stats{config, sketches, cardinalities}; + global_stats.hibf_layout = hibf_layouts[i]; + global_stats.print_header_to(std::cout); + global_stats.print_summary_to(dummy, std::cout); + } + + std::ofstream fout{config.output_filename, std::ios_base::app}; + chopper::layout::write_user_bins_to(filenames, mapping, fout); + config.write_to(fout); + layout.write_to(fout); + } return 0; } diff --git a/src/layout/input.cpp b/src/layout/input.cpp index 2b1a7459..bec0f519 100644 --- a/src/layout/input.cpp +++ b/src/layout/input.cpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include @@ -20,6 +21,7 @@ #include #include +#include namespace chopper::layout { @@ -70,10 +72,13 @@ std::tuple>, configuration, seqan::hibf::la read_layout_file(std::istream & stream) { std::vector> filenames = chopper::layout::read_filenames_from(stream); + chopper::configuration chopper_config; chopper_config.read_from(stream); + seqan::hibf::layout::layout hibf_layout{}; hibf_layout.read_from(stream); + return std::make_tuple(std::move(filenames), std::move(chopper_config), std::move(hibf_layout)); } diff --git a/src/layout/lsh_distributed_ibf_layout.cpp b/src/layout/lsh_distributed_ibf_layout.cpp index c4625e63..3b2837d3 100644 --- a/src/layout/lsh_distributed_ibf_layout.cpp +++ b/src/layout/lsh_distributed_ibf_layout.cpp @@ -14,7 +14,6 @@ #include #include #include -#include #include #include #include @@ -23,6 +22,7 @@ #include #include #include +#include #include #include @@ -34,70 +34,6 @@ 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. diff --git a/src/layout/output.cpp b/src/layout/output.cpp index 71dbd8b8..72a09d01 100644 --- a/src/layout/output.cpp +++ b/src/layout/output.cpp @@ -19,17 +19,19 @@ namespace chopper::layout { -void write_user_bins_to(std::vector> const & filenames, std::ostream & stream) +void write_user_bins_to(std::vector> const & filenames, + std::vector const & mapping, + std::ostream & stream) { stream << chopper::prefix::meta_chopper_user_bins_start << '\n'; size_t counter{}; - for (auto const & filenames_of_user_bin : filenames) + for (auto const & idx : mapping) { // the below will write lines like this: // @0 file1.fa file2.fa // @1 fileABC.fa stream << seqan::hibf::prefix::meta_header << counter++; - for (std::string const & filename : filenames_of_user_bin) + for (std::string const & filename : filenames[idx]) stream << ' ' << filename; stream << '\n'; } diff --git a/src/layout/partition_user_bins.cpp b/src/layout/partition_user_bins.cpp new file mode 100644 index 00000000..e47b8efd --- /dev/null +++ b/src/layout/partition_user_bins.cpp @@ -0,0 +1,1113 @@ +// --------------------------------------------------------------------------------------------------- +// 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 Clusters similar user bins by iterative MinHash LSH (ADD+OR scheme). + * \param[in] minHash_sketches The MinHash tables of all user bins. + * \param[in] cardinalities The cardinality of each user bin. + * \param[in] average_technical_bin_size Stop clustering as soon as the largest cluster's cardinality (sum of its user + * bins' cardinalities) reaches this. + * \param[in] number_of_bands LSH parameter b: the maximum number of rounds, one MinHash sketch per round. + * \param[in] number_of_hashes_per_band LSH parameter r: the number of hashes combined per LSH key. + */ +std::vector LSH_clustering(std::vector const & minHash_sketches, + std::vector const & cardinalities, + size_t const average_technical_bin_size, + size_t const number_of_bands, + size_t const number_of_hashes_per_band) +{ + + size_t const number_of_user_bins{cardinalities.size()}; + size_t const number_of_max_minHash_sketches{number_of_bands}; // LSH ADD+OR parameter b + size_t const minHash_sketche_size{number_of_hashes_per_band}; // LSH ADD+OR parameter r + + assert(number_of_user_bins == minHash_sketches.size()); + assert(!minHash_sketches.empty()); + assert(minHash_sketches[0].table.size() >= number_of_max_minHash_sketches); + assert(minHash_sketches[0].table[0].size() >= minHash_sketche_size); + + // initialise clusters with a signle user bin per cluster. + // clusters are either + // 1) of size 1; containing an id != position where the id points to the cluster it has been moved to + // e.g. cluster[Y] = {Z} (Y has been moved into Z, so Z could look likes this cluster[Z] = {Z, Y}) + // 2) of size >= 1; with the first entry beging id == position (a valid cluster) + // e.g. cluster[X] = {X} // valid singleton + // e.g. cluster[X] = {X, a, b, c, ...} // valid cluster with more joined entries + // The clusters could me moved recursively, s.t. + // cluster[A] = {B} + // cluster[B] = {C} + // cluster[C] = {C, A, B} // is valid cluster since cluster[C][0] == C; contains A and B + std::vector clusters; + clusters.reserve(number_of_user_bins); + + std::vector current_cluster_cardinality(number_of_user_bins); + size_t current_max_cluster_size{0}; + size_t current_number_of_sketch_hashes{minHash_sketche_size}; + size_t current_sketch_index{0}; + + for (size_t user_bin_idx = 0; user_bin_idx < number_of_user_bins; ++user_bin_idx) + { + clusters.emplace_back(user_bin_idx); + current_cluster_cardinality[user_bin_idx] = cardinalities[user_bin_idx]; + current_max_cluster_size = std::max(current_max_cluster_size, cardinalities[user_bin_idx]); + } + + // 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, current_number_of_sketch_hashes); + + // 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()); + + // LSH_fill_hashtable returns sorted lists without duplicates. + auto const begin = list.begin(); + auto const end = list.end(); + + 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 + size_t const representative_cluster_id = LSH_find_representative_cluster(clusters, *begin); + auto & representative_cluster = clusters[representative_cluster_id]; + assert(representative_cluster.is_valid(representative_cluster_id)); + assert(representative_cluster.id() == clusters[representative_cluster.id()].id()); + + for (auto current = begin + 1; current < end; ++current) + { + // 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 + size_t const next_cluster_id = LSH_find_representative_cluster(clusters, *current); + auto & next_cluster = clusters[next_cluster_id]; + + 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.size() == 0); + assert(next_cluster.has_been_moved()); + assert(representative_cluster.size() > 1); // there should be at least two user bins now + assert(representative_cluster.is_valid(representative_cluster_id)); // and it should still be valid + + current_cluster_cardinality[representative_cluster.id()] += + current_cluster_cardinality[next_cluster.id()]; + } + + current_max_cluster_size = *std::ranges::max_element(current_cluster_cardinality); + } + + ++current_sketch_index; + } + + return clusters; +} + +// At the end of this function I have clusters of clusters (Multi clusters). +// E.g. [1:[A,B], 2:[C,D,E,F], 3:[G]] +// Let A-G be (initial_)clusters of very similar user bins obtained by previous clustering +// Then 1,2, and 3 Are a MultiCluster each. The clusters within each MultiCluster are rather similar +// Now clusters from 1,2 and 3 respectively can be used to initialise a partition because between the MultiClusters +// there is not much similarity +std::vector +most_distant_LSH_partitioning(std::vector const & initial_clusters, + std::vector const & minHash_sketches, + chopper::configuration const & config) +{ + assert(!minHash_sketches.empty()); + assert(!minHash_sketches[0].table.empty()); + + size_t const number_of_user_bins{initial_clusters.size()}; // empty clusters are still left in the vector + size_t const number_of_max_minHash_sketches{minHash_sketches[0].table.size()}; // LSH ADD+OR parameter b + size_t current_number_of_sketch_hashes{5}; + + assert(number_of_user_bins == minHash_sketches.size()); + assert(minHash_sketches[0].table[0].size() >= current_number_of_sketch_hashes); + + size_t current_sketch_index{0}; + + size_t current_number_of_clusters = [&initial_clusters]() + { + size_t number{0}; + for (auto const & cluster : initial_clusters) + { + if (cluster.size()) + ++number; + } + return number; + }(); + + std::vector clusters; + clusters.reserve(number_of_user_bins); + + for (size_t user_bin_idx = 0; user_bin_idx < number_of_user_bins; ++user_bin_idx) + { + clusters.emplace_back(initial_clusters[user_bin_idx]); + assert(clusters[user_bin_idx].is_valid(user_bin_idx)); + } + + // refine clusters. I want to cluster as much as possible to initialise partitions with very distant clusters + while (current_number_of_clusters > (config.hibf_config.tmax * 2) + && 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, current_number_of_sketch_hashes); + + // 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()); + + // LSH_fill_hashtable returns sorted lists without duplicates. + auto const begin = list.begin(); + auto const end = list.end(); + + 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 + size_t const representative_cluster_id = LSH_find_representative_cluster(clusters, *begin); + auto & representative_cluster = clusters[representative_cluster_id]; + assert(representative_cluster.is_valid(representative_cluster_id)); + assert(representative_cluster.id() == clusters[representative_cluster.id()].id()); + + for (auto current = begin + 1; current < end; ++current) + { + // 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 + size_t const next_cluster_id = LSH_find_representative_cluster(clusters, *current); + auto & next_cluster = clusters[next_cluster_id]; + + 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.size() == 0); + assert(next_cluster.has_been_moved()); + assert(representative_cluster.size() > 1); // there should be at least two user bins now + assert(representative_cluster.is_valid(representative_cluster_id)); // and it should still be valid + + --current_number_of_clusters; + } + } + + ++current_sketch_index; + + // The LSH matching looser over time: Every 4 rounds, the number of hashes that must agree drops by one. + // More and more dissimilar clusters merge, until there are at most 2·tmax clusters or the rounds run out. + if (current_sketch_index % 4 == 0 && current_number_of_sketch_hashes > 1) + --current_number_of_sketch_hashes; + } + + return clusters; +} + +void post_process_partition_clusters(std::vector & clusters, + std::vector const & cardinalities, + chopper::configuration const & config) +{ + assert(clusters.size() == cardinalities.size()); + assert(cardinalities.size() > config.number_of_partitions); + + // 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); + } + + // Edge case: The sequences are all very similar and there are fewer clusters than partitions + // Solution: Split largest cluster until there are as many clusters as partitions + size_t num_non_empty_clusters = std::ranges::count_if(clusters, + [](auto const & c) + { + return !c.empty(); + }); + while (num_non_empty_clusters < config.number_of_partitions) + { + auto it_larges_cluster = std::max_element(clusters.begin(), + clusters.end(), + [](auto const & a, auto const & b) + { + return a.size() < b.size(); + }); + auto it_empty_cluster = std::ranges::find_if(clusters.begin(), + clusters.end(), + [](auto const & c) + { + return c.empty(); + }); + assert(it_larges_cluster != clusters.end()); // There is always a non-empty cluster + assert( + it_empty_cluster + != clusters + .end()); // There is always an empty cluster because there are more user bis than partitions when we arrive here + it_larges_cluster->move_part_to(*it_empty_cluster, it_larges_cluster->size() / 2); + ++num_non_empty_clusters; + } + + // push largest p clusters to the front + auto cluster_size_cmp = [&cardinalities](auto const & v1, auto const & v2) + { + if (v2.size() == v1.size() && !v2.empty()) + return cardinalities[v2.contained_user_bins()[0]] < cardinalities[v1.contained_user_bins()[0]]; + return v2.size() < v1.size(); + }; + std::partial_sort(clusters.begin(), + clusters.begin() + config.number_of_partitions, + clusters.end(), + cluster_size_cmp); + + // 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. + auto compare_cardinality_and_move_empty_clusters_to_the_end = [&cardinalities](auto const & v1, auto const & v2) + { + if (v1.empty()) + return false; // v1 can never be larger than v2 then + + if (v2.empty()) // and v1 is not, since the first if would catch + return true; + + return cardinalities[v2.contained_user_bins()[0]] < cardinalities[v1.contained_user_bins()[0]]; + }; + std::sort(clusters.begin() + config.number_of_partitions, + clusters.end(), + compare_cardinality_and_move_empty_clusters_to_the_end); + + assert(clusters.size() < 2 || clusters[0].size() >= clusters[1].size()); // sanity check + // sanity check. Note that a cluster's id() is not its largest user bin; contained_user_bins()[0] is. + assert(std::is_sorted(clusters.begin() + config.number_of_partitions, + clusters.end(), + compare_cardinality_and_move_empty_clusters_to_the_end)); + + // debug + for (size_t cidx = 1; cidx < clusters.size(); ++cidx) + { + if (clusters[cidx - 1].empty() + && !clusters[cidx].empty()) // once empty - always empty; all empty clusters should be at the end + throw std::runtime_error{"sorting did not work"}; + } + // debug +} + +void post_process_partition_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); + } + + auto number_of_contained_user_bins = [](MultiCluster const & mc) + { + size_t sum{0}; + for (auto const & user_bin_cluster : mc.contained_user_bins()) + sum += user_bin_cluster.size(); + return sum; + }; + + // A MultiCluster can only be created from a Cluster, so the split parts are rebuilt from their (sub)clusters. + auto make_multi_cluster = [](size_t const id, auto const first, auto const last) + { + assert(first != last); + auto to_cluster = [id](std::vector const & user_bin_cluster) + { + assert(!user_bin_cluster.empty()); + Cluster cluster{id, user_bin_cluster[0]}; + for (size_t i = 1; i < user_bin_cluster.size(); ++i) + cluster.add_user_bin(user_bin_cluster[i]); + return cluster; + }; + + MultiCluster result{to_cluster(*first)}; + for (auto it = first + 1; it != last; ++it) + { + MultiCluster tmp{to_cluster(*it)}; + tmp.move_to(result); + } + return result; + }; + + // Edge case: The sequences are all very similar and there are fewer clusters than partitions + // Solution: Split largest cluster (by number of user bins) until there are as many clusters as partitions. + // A cluster with several (sub)clusters is split between its (sub)clusters, s.t. very similar user bins stay together. + // A cluster with a single (sub)cluster is split between its user bins. + size_t num_non_empty_clusters = std::ranges::count_if(clusters, + [](auto const & c) + { + return !c.empty(); + }); + while (num_non_empty_clusters < config.number_of_partitions) + { + auto it_largest_cluster = + std::ranges::max_element(clusters, std::ranges::less{}, number_of_contained_user_bins); + auto it_empty_cluster = std::ranges::find_if(clusters, + [](auto const & c) + { + return c.empty(); + }); + assert( + it_empty_cluster + != clusters + .end()); // There is always an empty cluster because there are more user bins than partitions when we arrive here + assert(number_of_contained_user_bins(*it_largest_cluster) + >= 2); // pigeonhole: more user bins than non-empty clusters + + std::vector> parts = + it_largest_cluster->contained_user_bins(); // sorted by size, descending + if (parts.size() == 1) // split the user bins (sorted by cardinality, descending) of the single (sub)cluster + { + auto & user_bin_cluster = parts[0]; + size_t const half = user_bin_cluster.size() / 2; + std::vector second_half(user_bin_cluster.begin() + half, user_bin_cluster.end()); + user_bin_cluster.resize(half); + parts.push_back(std::move(second_half)); + } + + size_t const keep = parts.size() / 2; + *it_empty_cluster = make_multi_cluster(it_empty_cluster->id(), parts.begin() + keep, parts.end()); + *it_largest_cluster = make_multi_cluster(it_largest_cluster->id(), parts.begin(), parts.begin() + keep); + ++num_non_empty_clusters; + } + + // push largest p clusters to the front + auto cluster_size_cmp = [&cardinalities](auto const & mc1, auto const & mc2) + { + if (mc1.empty()) + return false; // mc1 can never be larger than mc2 then + if (mc2.empty()) // and mc1 is not, since the first if would catch + return true; + + auto const & largest_c1 = mc1.contained_user_bins()[0]; // first cluster is the largest because of sorting above + auto const & largest_c2 = mc2.contained_user_bins()[0]; // first cluster is the largest because of sorting above + + assert(!largest_c1.empty()); // should never be empty + assert(!largest_c2.empty()); // should never be empty + + if (largest_c2.size() == largest_c1.size()) + return cardinalities[largest_c2[0]] < cardinalities[largest_c1[0]]; + return largest_c2.size() < largest_c1.size(); + }; + std::partial_sort(clusters.begin(), + clusters.begin() + config.number_of_partitions, + clusters.end(), + cluster_size_cmp); + + auto move_empty_clusters_to_the_end = [](auto const & v1, auto const & v2) + { + return v2.size() < v1.size(); + }; + std::sort(clusters.begin() + config.number_of_partitions, clusters.end(), move_empty_clusters_to_the_end); + + // debug sanity check + for (size_t cidx = 1; cidx < clusters.size(); ++cidx) + { + if (clusters[cidx - 1].empty() + && !clusters[cidx].empty()) // once empty - always empty; all empty clusters should be at the end + throw std::runtime_error{"sorting did not work"}; + } + // debug +} + +std::vector +sim_dist_LSH_partitioning(std::vector const & minHash_sketches, + std::vector const & cardinalities, + size_t const average_technical_bin_size, + chopper::configuration const & config) +{ + std::vector clusters = LSH_clustering(minHash_sketches, + cardinalities, + average_technical_bin_size, + 3, // b + minHash_sketches[0].table[0].size()); // r + std::vector multi_clusters = most_distant_LSH_partitioning(clusters, minHash_sketches, config); + post_process_partition_clusters(multi_clusters, cardinalities, config); + return multi_clusters; +} + +void find_best_partition(chopper::configuration const & config, + size_t const cardinality_per_part, + std::vector const & cluster, + std::vector const & cardinalities, + std::vector const & sketches, + std::vector> & positions, + std::vector & partition_sketches, + std::vector & partition_cardinality) +{ + seqan::hibf::sketch::hyperloglog const current_sketch = [&sketches, &cluster, &config]() + { + seqan::hibf::sketch::hyperloglog result{config.hibf_config.sketch_bits}; + + for (size_t const user_bin_idx : cluster) + result.merge(sketches[user_bin_idx]); + + return result; + }(); + + size_t const current_cardinality = [&cardinalities, &cluster]() + { + size_t result{0}; + + for (size_t const user_bin_idx : cluster) + result += cardinalities[user_bin_idx]; + + return result; + }(); + + // 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." + double best_subsume_ratio{0.0}; + size_t best_p{0}; + bool best_p_found{false}; + + for (size_t p = 0; p < config.number_of_partitions; ++p) + { + seqan::hibf::sketch::hyperloglog tmp = current_sketch; + tmp.merge(partition_sketches[p]); + size_t const union_size = tmp.estimate(); + size_t const current_partition_size = partition_sketches[p].estimate(); + assert(union_size >= current_partition_size); // partition cannot shrink by adding kmers + size_t const change = union_size - current_partition_size; + // HyperLogLog noise can make `change` larger than the estimate of the cluster (e.g., for disjoint sets). + double const intersection = std::max(0.0, current_sketch.estimate() - static_cast(change)); + double const subsume_ratio = static_cast(intersection) / current_partition_size; + + if (subsume_ratio > best_subsume_ratio && partition_cardinality[p] < cardinality_per_part + && (current_cardinality + partition_cardinality[p]) < (cardinality_per_part * 1.2)) + { + best_subsume_ratio = subsume_ratio; + best_p = p; + best_p_found = true; + } + } + + if (!best_p_found) + { + // No partition satisfies the size limits + // (e.g. too few user bins per partition, or a single user bin larger than 1.2 · cardinality_per_part, + // or a user bin with no sequence content and cardinality 0) + + // If we are working on a cluster with several user bins, assign user bins individually + if (cluster.size() > 1) + { + for (size_t const user_bin_idx : cluster) + { + find_best_partition(config, + cardinality_per_part, + {user_bin_idx}, + cardinalities, + sketches, + positions, + partition_sketches, + partition_cardinality); + } + return; + } + else // For single user bins, fall back to the least filled partition. + { + best_p = std::ranges::min_element(partition_cardinality) - partition_cardinality.begin(); + } + } + + // now that we know which partition fits best (`best_p`), add those indices to it + for (size_t const user_bin_idx : cluster) + { + positions[best_p].push_back(user_bin_idx); + partition_cardinality[best_p] += cardinalities[user_bin_idx]; + } + partition_sketches[best_p].merge(current_sketch); +} + +void partition_user_bins(chopper::configuration const & config, + std::vector const & cardinalities, + std::vector const & sketches, + std::vector const & minHash_sketches, + std::vector> & partitions) +{ + if (config.partitioning_approach < partitioning_scheme::blocked + || config.partitioning_approach > partitioning_scheme::lsh_sim) + { + throw std::invalid_argument{"Unknown partitioning approach " + + std::to_string(static_cast(config.partitioning_approach)) + "."}; + } + + if (cardinalities.size() < config.number_of_partitions) + { + throw std::invalid_argument{"The number of partitions (" + std::to_string(config.number_of_partitions) + + ") must not exceed the number of user bins (" + + std::to_string(cardinalities.size()) + ")."}; + } + + // all approaches need sorted positions + std::vector const sorted_positions = [&cardinalities]() + { + std::vector ps; + ps.resize(cardinalities.size()); + std::iota(ps.begin(), ps.end(), 0); + seqan::hibf::sketch::toolbox::sort_by_cardinalities(cardinalities, ps); + return ps; + }(); + + size_t const sum_of_cardinalities = [&cardinalities]() + { + size_t sum{0}; + for (size_t const card : cardinalities) + sum += card; + return sum; + }(); + + if (config.partitioning_approach == partitioning_scheme::blocked) + { + size_t const u_bins_per_part = cardinalities.size() / config.number_of_partitions; + // for small edge cases, just use u_bins_per_part, + // otherwise, use next_multiple_of_64(ceil(sqrt(u_bins_per_part))) to have am ibf friendly block size. + // That way, we loop over the partitions, filling them block wise, again and again. + size_t const block_size = + (u_bins_per_part < 64) + ? u_bins_per_part + : chopper::next_multiple_of_64(static_cast(std::ceil(std::sqrt(u_bins_per_part)))); + + size_t current_part{0u}; + size_t current_block_count{0}; + + for (size_t const current_user_bin_id : sorted_positions) + { + partitions[current_part].push_back(current_user_bin_id); + ++current_block_count; + + if (current_block_count >= block_size) + { + current_block_count = 0; + ++current_part; + if (current_part == config.number_of_partitions) // we need to circle back to the first partition + current_part = 0; + } + } + } + else if (config.partitioning_approach == partitioning_scheme::sorted) + { + size_t const cardinality_per_part = sum_of_cardinalities / config.number_of_partitions; + + std::vector partition_cardinalities(config.number_of_partitions, 0); + size_t current_part{0}; + + for (size_t i = 0; i < sorted_positions.size(); ++i) + { + partitions[current_part].push_back(sorted_positions[i]); + partition_cardinalities[current_part] += cardinalities[sorted_positions[i]]; + + if (partition_cardinalities[current_part] >= cardinality_per_part /*partition is full*/ + || (sorted_positions.size() - i + <= config.number_of_partitions - current_part) /* small number edge case:not enough UBs left*/) + { + ++current_part; + // If cardinality_per_part was not chosen optimally + // we might need start at position 0 again. + // If all paritions are full, this will add one user bin at a time + // to each partition until no user bins remain. + if (current_part >= config.number_of_partitions) + current_part = 0; + } + } + } + else if (config.partitioning_approach == partitioning_scheme::folded) + { + size_t const cardinality_per_part_halved = + seqan::hibf::divide_and_ceil(sum_of_cardinalities, config.number_of_partitions * 2); + + size_t current_cardinality{0u}; + std::vector const parts = [&config]() + { + size_t const len{config.number_of_partitions}; + std::vector result(len * 2); + std::iota(result.begin(), result.begin() + len, 0); + std::copy(result.rbegin() + len, result.rend(), result.begin() + len); + return result; + }(); + size_t current_part{0}; + + for (size_t const current_user_bin_id : sorted_positions) + { + partitions[parts[current_part]].push_back(current_user_bin_id); + current_cardinality += cardinalities[current_user_bin_id]; + + if (current_cardinality >= cardinality_per_part_halved) + { + current_cardinality = 0; + // Remaining user bins (with cardinality 0) stay in the last part. + if (current_part + 1 < parts.size()) + ++current_part; + } + } + } + else if (config.partitioning_approach == partitioning_scheme::weighted_fold) + { + size_t const cardinality_per_part = + seqan::hibf::divide_and_ceil(sum_of_cardinalities, config.number_of_partitions); + size_t const u_bins_per_part = cardinalities.size() / config.number_of_partitions; + + // The unassigned user bins are sorted_positions[big, small_end). Large user bins are taken from the front, + // small ones from the back. + size_t big{0}; + size_t small_end{sorted_positions.size()}; + + // How far a partition with the given cardinality and number of user bins is off the targets. Smaller is better. + auto compute_score = [&](size_t const cardinality, size_t const number_of_user_bins) + { + double const weight = static_cast(cardinality) / cardinality_per_part; + double const amount = static_cast(number_of_user_bins) / u_bins_per_part; + return std::abs(1.0 - weight) + std::abs(1.0 - amount); + }; + + for (size_t current_part = 0; current_part + 1 < config.number_of_partitions; ++current_part) + { + std::vector large_bins; + std::vector small_bins; + size_t current_cardinality{0}; + + // first add all large bins that fit + while (current_cardinality < cardinality_per_part && big < small_end) + { + large_bins.push_back(sorted_positions[big]); + current_cardinality += cardinalities[sorted_positions[big]]; + ++big; + } + + double local_optimum = compute_score(current_cardinality, large_bins.size()); + + // then remove large bins and add small bins until a local optimum is reached + while (!large_bins.empty()) + { + // remove the last large bin ... + size_t trial_cardinality = current_cardinality - cardinalities[large_bins.back()]; + size_t trial_number_of_user_bins = large_bins.size() - 1 + small_bins.size(); + size_t trial_small_end = small_end; + double trial_score = compute_score(trial_cardinality, trial_number_of_user_bins); + + // ... and add small bins as long as this improves the score. The removed large bin is + // sorted_positions[big - 1] and must not be added as a small bin. + while (trial_small_end > big) + { + size_t const candidate = sorted_positions[trial_small_end - 1]; + double const score = + compute_score(trial_cardinality + cardinalities[candidate], trial_number_of_user_bins + 1); + + if (score >= trial_score) // smaller is better + break; + + trial_cardinality += cardinalities[candidate]; + ++trial_number_of_user_bins; + --trial_small_end; + trial_score = score; + } + + if (trial_number_of_user_bins == 0 || local_optimum < trial_score) // score would increase. Stop + break; + + // update + large_bins.pop_back(); + --big; + for (size_t pos = small_end; pos > trial_small_end; --pos) + small_bins.push_back(sorted_positions[pos - 1]); + small_end = trial_small_end; + current_cardinality = trial_cardinality; + local_optimum = trial_score; + } + + partitions[current_part] = std::move(large_bins); + partitions[current_part].insert(partitions[current_part].end(), small_bins.begin(), small_bins.end()); + } + + // remaining user bins go to last partition + for (; big < small_end; ++big) + partitions[config.number_of_partitions - 1].push_back(sorted_positions[big]); + } + else if (config.partitioning_approach == partitioning_scheme::similarity) + { + uint8_t const sketch_bits{config.hibf_config.sketch_bits}; + std::vector partition_sketches(config.number_of_partitions, + seqan::hibf::sketch::hyperloglog(sketch_bits)); + std::vector partition_cardinality(config.number_of_partitions, 0u); + + size_t const cardinality_per_part = + seqan::hibf::divide_and_ceil(sum_of_cardinalities, config.number_of_partitions); + size_t const block_size = [&]() + { + size_t const u_bins_per_part = cardinalities.size() / config.number_of_partitions; + return std::min(u_bins_per_part, + chopper::next_multiple_of_64(static_cast(std::ceil(std::sqrt(u_bins_per_part))))); + }(); + size_t const number_of_blocks = seqan::hibf::divide_and_ceil(cardinalities.size(), block_size); + + // don't move from largest to smallest but pick the next block to process randomly. + // this probably leads to more evenly distributed partitions (evenly in terms of number of user bins) + std::vector indices(number_of_blocks); + std::iota(indices.begin(), indices.end(), 0); + std::random_device shuffle_random_device; + std::mt19937 shuffle_engine(shuffle_random_device()); + std::shuffle(indices.begin(), indices.end(), shuffle_engine); + + // initialise partitions with the first random config.number_of_partitions blocks + assert(number_of_blocks >= config.number_of_partitions); + for (size_t p = 0; p < config.number_of_partitions; ++p) + { + size_t const i = indices[p]; + size_t const end = (i == (indices.size() - 1) && (cardinalities.size() % block_size != 0)) + ? (cardinalities.size() % block_size) + : block_size; + for (size_t x = 0; x < end; ++x) + { + assert(block_size * i + x < sorted_positions.size()); + partition_sketches[p].merge(sketches[sorted_positions[block_size * i + x]]); + partition_cardinality[p] += cardinalities[sorted_positions[block_size * i + x]]; + partitions[p].push_back(sorted_positions[block_size * i + x]); + } + } + + auto has_been_processed_in_init = [&config, &indices, &block_size](size_t const i) + { + bool result{false}; + for (size_t p = 0; p < config.number_of_partitions; ++p) + result |= ((indices[p] * block_size) == i); + return result; + }; + + // assign the rest by similarity, user in by user bin (single bins, no blocks) + for (size_t i = 0; i < sorted_positions.size(); ++i) + { + if (has_been_processed_in_init(i)) + { + i += (block_size - 1); // -1 because there will be an increment after continue + continue; + } + + find_best_partition(config, + cardinality_per_part, + {sorted_positions[i]}, + cardinalities, + sketches, + partitions, + partition_sketches, + partition_cardinality); + } + } + else if (config.partitioning_approach == partitioning_scheme::lsh) + { + uint8_t const sketch_bits{config.hibf_config.sketch_bits}; + std::vector partition_sketches(config.number_of_partitions, + seqan::hibf::sketch::hyperloglog(sketch_bits)); + std::vector partition_cardinality(config.number_of_partitions, 0u); + + size_t const cardinality_per_part = + seqan::hibf::divide_and_ceil(sum_of_cardinalities, config.number_of_partitions); + + // initial partitioning using locality sensitive hashing (LSH) + std::vector raw_clusters = LSH_clustering(minHash_sketches, + cardinalities, + cardinality_per_part, + seqan::hibf::sketch::minhashes::num_sketches, // b + seqan::hibf::sketch::minhashes::sketch_size); // r + post_process_partition_clusters(raw_clusters, cardinalities, config); + std::vector const & clusters = raw_clusters; // clusters do not change from this point on + + // initialise partitions with the first p largest clusters (post_processing sorts by size) + std::vector spilled_user_bins{}; // user bins that spilled over past the last partition + size_t cidx{0}; // current cluster index + for (size_t p = 0; p < config.number_of_partitions; ++p) + { + assert(!clusters[cidx].empty()); // post_process_cluster ensures enough non-empty clusters + + bool p_has_been_incremented{false}; + + for (size_t uidx = 0; uidx < clusters[cidx].size(); ++uidx) + { + size_t const user_bin_idx = clusters[cidx].contained_user_bins()[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 (partition_cardinality[p] > cardinality_per_part) + { + if (p + 1 == config.number_of_partitions) // No partition left to spill into + { + spilled_user_bins.push_back(user_bin_idx); // assigned after seeding + continue; + } + else + { + ++p; + p_has_been_incremented = true; + } + } + + partition_sketches[p].merge(sketches[user_bin_idx]); + partition_cardinality[p] += cardinalities[user_bin_idx]; + partitions[p].push_back(user_bin_idx); + } + + if (p_has_been_incremented && partition_cardinality[p] < cardinality_per_part) + { + // p will be incremeted again in the next iteration of this for loop + // since we have incremented p previously because of some user bins that spilled over in this currnet p + // we want to decrease p here, s.t. we stay at the same p. + // this ensures that the current p doesn't only have a small rest that just didn't fit into the + // previous partition + assert(p > 0); + --p; + } + + ++cidx; + } + + for (size_t const user_bin_idx : spilled_user_bins) + find_best_partition(config, + cardinality_per_part, + {user_bin_idx}, + cardinalities, + sketches, + partitions, + partition_sketches, + partition_cardinality); + + // assign the rest by similarity + assert(cardinalities.size() == clusters.size()); + for (; cidx < cardinalities.size(); ++cidx) + { + auto const & cluster = clusters[cidx]; + + if (cluster.empty()) // non valid clusters have been removed. Since list was sorted, we can abort + break; + + find_best_partition(config, + cardinality_per_part, + cluster.contained_user_bins(), + cardinalities, + sketches, + partitions, + partition_sketches, + partition_cardinality); + } + } + else if (config.partitioning_approach == partitioning_scheme::lsh_sim) + { + uint8_t const sketch_bits{config.hibf_config.sketch_bits}; + std::vector partition_sketches(config.number_of_partitions, + seqan::hibf::sketch::hyperloglog(sketch_bits)); + std::vector partition_cardinality(config.number_of_partitions, 0u); + + size_t const cardinality_per_part = + seqan::hibf::divide_and_ceil(sum_of_cardinalities, config.number_of_partitions); + + // initial partitioning using locality sensitive hashing (LSH) + std::vector const clusters = + sim_dist_LSH_partitioning(minHash_sketches, cardinalities, cardinality_per_part, config); + + // initialise partitions with the first p largest clusters (post processing sorts by size) + std::vector spilled_user_bins{}; // user bins that spilled over past the last partition + size_t cidx{0}; // current cluster index + for (size_t p = 0; p < config.number_of_partitions; ++p) + { + // There may be fewer (non-empty) clusters than partitions. The remaining partitions stay empty, which is + // reported at the end. + if (clusters[cidx].empty()) + break; + + bool p_has_been_incremented{false}; + auto const & cluster = clusters[cidx].contained_user_bins()[0]; + + for (size_t uidx = 0; uidx < cluster.size(); ++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 (partition_cardinality[p] > cardinality_per_part) + { + if (p + 1 == config.number_of_partitions) // No partition left to spill into + { + spilled_user_bins.push_back(user_bin_idx); // assigned after seeding + continue; + } + else + { + ++p; + p_has_been_incremented = true; + } + } + + partition_sketches[p].merge(sketches[user_bin_idx]); + partition_cardinality[p] += cardinalities[user_bin_idx]; + partitions[p].push_back(user_bin_idx); + } + + if (p_has_been_incremented && partition_cardinality[p] < cardinality_per_part) + { + // p will be incremeted again in the next iteration of this for loop + // since we have incremented p previously because of some user bins that spilled over in this currnet p + // we want to decrease p here, s.t. we stay at the same p. + // this ensures that the current p doesn't only have a small rest that just didn't fit into the + // previous partition + assert(p > 0); + --p; + } + + ++cidx; + } + + for (size_t const user_bin_idx : spilled_user_bins) + find_best_partition(config, + cardinality_per_part, + {user_bin_idx}, + cardinalities, + sketches, + partitions, + partition_sketches, + partition_cardinality); + + // assign the rest, cluster by cluster (because the clusters are very similar), by similarity + assert(cardinalities.size() == clusters.size()); + + std::vector> remaining_clusters{}; + + for (size_t i = 0; i < clusters.size(); ++i) + { + if (clusters[i].empty()) + break; + // if i < cidx the first cluster was already used in initialising above + auto start = clusters[i].contained_user_bins().begin() + ((i < cidx) ? 1 : 0); + remaining_clusters.insert(remaining_clusters.end(), start, clusters[i].contained_user_bins().end()); + } + auto cmp = [&cardinalities](auto const & c1, auto const & c2) + { + return cardinalities[c2[0]] < cardinalities[c1[0]]; + }; + std::sort(remaining_clusters.begin(), remaining_clusters.end(), cmp); + + for (auto const & cluster : remaining_clusters) + { + find_best_partition(config, + cardinality_per_part, + cluster, + cardinalities, + sketches, + partitions, + partition_sketches, + partition_cardinality); + } + } + +#ifndef NDEBUG + // sanity check in debug + std::vector user_bins{}; + for (auto const & partition : partitions) + { + user_bins.insert(user_bins.end(), partition.begin(), partition.end()); + } + if (!std::ranges::is_permutation(user_bins, seqan::hibf::iota_vector(sketches.size()))) + { + std::string str{"Not all user bins have been assigned to the "}; + str += std::to_string(partitions.size()); + str += " partitions! ("; + str += std::to_string(user_bins.size()); + str += "/"; + str += std::to_string(sketches.size()); + str += ")\n"; + for (auto const & p : partitions) + { + str += "["; + for (auto const h : p) + { + str += std::to_string(h); + str += ","; + } + str.back() = ']'; + str += '\n'; + } + + throw std::logic_error{str}; + } +#endif + + // Each partition is laid out as a separate HIBF, which needs at least one user bin. + for (size_t p = 0; p < partitions.size(); ++p) + { + if (partitions[p].empty()) + { + throw std::runtime_error{"Partition " + std::to_string(p) + " of " + std::to_string(partitions.size()) + + " is empty. Use fewer partitions or another partitioning approach."}; + } + } +} + +} // namespace chopper::layout diff --git a/src/set_up_parser.cpp b/src/set_up_parser.cpp index e34f6a2a..261ee7ca 100644 --- a/src/set_up_parser.cpp +++ b/src/set_up_parser.cpp @@ -86,6 +86,35 @@ void set_up_parser(sharg::parser & parser, configuration & config) .default_message = "", .validator = sharg::output_file_validator{}}); + parser.add_option( + config.number_of_partitions, + sharg::config{ + .short_id = '\0', + .long_id = "number-of-partitions", + .description = + "The number of partitions of the HIBF. If at least 2, the user bins are distributed onto this many " + "partitions, each of which is laid out as a separate HIBF. This reduces the size of each index and " + "thereby the peak memory usage of searching with the HIBF.", + .advanced = true}); + + sharg::value_list_validator partitioning_approach_validator{ + chopper::partitioning_scheme::blocked, + chopper::partitioning_scheme::sorted, + chopper::partitioning_scheme::folded, + chopper::partitioning_scheme::weighted_fold, + chopper::partitioning_scheme::similarity, + chopper::partitioning_scheme::lsh, + chopper::partitioning_scheme::lsh_sim}; + + parser.add_option(config.partitioning_approach, + sharg::config{.short_id = '\0', + .long_id = "partitioning-approach", + .description = "How the user bins are distributed onto the partitions of a " + "partitioned HIBF." + "This is only configurable for debugging.", + .advanced = true, + .validator = partitioning_approach_validator}); + parser.add_option( config.hibf_config.tmax, sharg::config{ diff --git a/src/util/display_layout/general.cpp b/src/util/display_layout/general.cpp index 17fcd7bd..e9cbe758 100644 --- a/src/util/display_layout/general.cpp +++ b/src/util/display_layout/general.cpp @@ -184,30 +184,16 @@ void process_and_write_records_to(std::vector & records, std::ostream & stream << std::flush; } -int execute(config const & cfg) +int execute(config const & cfg, + std::vector> const & filenames, + chopper::configuration const & chopper_config, + seqan::hibf::layout::layout & hibf_layout) { - std::ifstream layout_file{cfg.input}; - - if (!layout_file.good() || !layout_file.is_open()) - throw std::logic_error{"Could not open file " + cfg.input.string() + " for reading"}; - -// https://godbolt.org/z/PeKnxzjn1 -#if defined(__clang__) - auto tuple = chopper::layout::read_layout_file(layout_file); - // https://godbolt.org/z/WoWf55KPb - auto filenames = std::move(std::get<0>(tuple)); - auto chopper_config = std::move(std::get<1>(tuple)); - auto hibf_layout = std::move(std::get<2>(tuple)); -#else - auto [filenames, chopper_config, hibf_layout] = chopper::layout::read_layout_file(layout_file); -#endif auto const & hibf_config = chopper_config.hibf_config; // read_layout_file does not validate the layout. hibf_layout.validate(hibf_config); - layout_file.close(); - // multiplied to cardinality of a merged bin double const merged_correction = seqan::hibf::layout::compute_relaxed_fpr_correction( {.fpr = chopper_config.hibf_config.maximum_fpr, @@ -221,7 +207,9 @@ int execute(config const & cfg) .hash_count = chopper_config.hibf_config.number_of_hash_functions, .t_max = chopper_config.hibf_config.tmax}); - std::ofstream output_stream{cfg.output}; + // Append: A partitioned HIBF layout file contains multiple layouts, which are written one after another. + // The ouput path is validated with output_file_validator, so it must not exist before and we can safely append + std::ofstream output_stream{cfg.output, std::ios_base::app}; if (!output_stream.good() || !output_stream.is_open()) throw std::logic_error{"Could not open file " + cfg.output.string() + " for reading"}; @@ -387,5 +375,16 @@ int execute(config const & cfg) void execute_general(config const & cfg) { - execute(cfg); + std::ifstream layout_file{cfg.input}; + + if (!layout_file.good() || !layout_file.is_open()) + throw std::logic_error{"Could not open file " + cfg.input.string() + " for reading"}; + + while (layout_file.good()) + { + auto [filenames, chopper_config, hibf_layout] = chopper::layout::read_layout_file(layout_file); + execute(cfg, filenames, chopper_config, hibf_layout); + } + + layout_file.close(); } diff --git a/src/util/display_layout/sizes.cpp b/src/util/display_layout/sizes.cpp index e2e37d58..b901a544 100644 --- a/src/util/display_layout/sizes.cpp +++ b/src/util/display_layout/sizes.cpp @@ -282,58 +282,50 @@ void execute_general_stats(config const & cfg) if (!layout_file.good() || !layout_file.is_open()) throw std::logic_error{"Could not open file " + cfg.input.string() + " for reading"}; -// https://godbolt.org/z/PeKnxzjn1 -#if defined(__clang__) - auto tuple = chopper::layout::read_layout_file(layout_file); - // https://godbolt.org/z/WoWf55KPb - auto filenames = std::move(std::get<0>(tuple)); - auto chopper_config = std::move(std::get<1>(tuple)); - auto hibf_layout = std::move(std::get<2>(tuple)); -#else - auto [filenames, chopper_config, hibf_layout] = chopper::layout::read_layout_file(layout_file); -#endif - - // read_layout_file does not validate the layout. - hibf_layout.validate(chopper_config.hibf_config); - - // Prepare configs - chopper_config.hibf_config.threads = cfg.threads; - auto input_lambda = [&filenames, &chopper_config](size_t const user_bin_id, seqan::hibf::insert_iterator it) + // Output. A partitioned HIBF layout file contains multiple layouts, whose stats are written one after another. + std::ofstream output_stream{cfg.output}; + if (!output_stream.good() || !output_stream.is_open()) + throw std::logic_error{"Could not open file " + cfg.output.string() + " for reading"}; + + while (layout_file.good()) { - std::vector current_kmers; + auto [filenames, chopper_config, hibf_layout] = chopper::layout::read_layout_file(layout_file); + hibf_layout.validate(chopper_config.hibf_config); - for (std::string const & filename : filenames[user_bin_id]) - process_file(filename, current_kmers, chopper_config.k, chopper_config.window_size); + // Prepare configs + chopper_config.hibf_config.threads = cfg.threads; + auto input_lambda = [&filenames, &chopper_config](size_t const user_bin_id, seqan::hibf::insert_iterator it) + { + std::vector current_kmers; - for (auto const kmer : current_kmers) - it = kmer; - }; - chopper_config.hibf_config.input_fn = input_lambda; - auto const & hibf_config = chopper_config.hibf_config; + for (std::string const & filename : filenames[user_bin_id]) + process_file(filename, current_kmers, chopper_config.k, chopper_config.window_size); - // Prepare stats - size_t const number_of_ibfs = hibf_layout.max_bins.size() + 1u; - std::vector stats(number_of_ibfs); + for (auto const kmer : current_kmers) + it = kmer; + }; + chopper_config.hibf_config.input_fn = input_lambda; + auto const & hibf_config = chopper_config.hibf_config; - // Prepare data - seqan::hibf::build::build_data data{.config = hibf_config, .ibf_graph = {hibf_layout}}; - seqan::hibf::layout::graph::node const & root_node = data.ibf_graph.root; - size_t const t_max{root_node.number_of_technical_bins}; - data.fpr_correction = seqan::hibf::layout::compute_fpr_correction( - {.fpr = hibf_config.maximum_fpr, .hash_count = hibf_config.number_of_hash_functions, .t_max = t_max}); + // Prepare stats + size_t const number_of_ibfs = hibf_layout.max_bins.size() + 1u; + std::vector stats(number_of_ibfs); - // Get stats - hierarchical_stats(stats, root_node, data); + // Prepare data + seqan::hibf::build::build_data data{.config = hibf_config, .ibf_graph = {hibf_layout}}; + seqan::hibf::layout::graph::node const & root_node = data.ibf_graph.root; + size_t const t_max{root_node.number_of_technical_bins}; + data.fpr_correction = seqan::hibf::layout::compute_fpr_correction( + {.fpr = hibf_config.maximum_fpr, .hash_count = hibf_config.number_of_hash_functions, .t_max = t_max}); - // Get stats per level - per_level_stats const level_stats{stats}; + // Get stats + hierarchical_stats(stats, root_node, data); - // Output - std::ofstream output_stream{cfg.output}; - if (!output_stream.good() || !output_stream.is_open()) - throw std::logic_error{"Could not open file " + cfg.output.string() + " for reading"}; + // Get stats per level + per_level_stats const level_stats{stats}; - level_stats.print(output_stream, hibf_layout.user_bins.size()); + level_stats.print(output_stream, hibf_layout.user_bins.size()); + } } void execute_sizes(config const & cfg) diff --git a/test/api/config_test.cpp b/test/api/config_test.cpp index 352ee060..81695d38 100644 --- a/test/api/config_test.cpp +++ b/test/api/config_test.cpp @@ -19,6 +19,8 @@ chopper::configuration generate_config() config.window_size = 24; config.disable_sketch_output = true; config.precomputed_files = true; + config.number_of_partitions = 4; + config.partitioning_approach = chopper::partitioning_scheme::lsh_sim; config.output_filename = "file.layout"; config.determine_best_tmax = true; config.force_all_binnings = true; @@ -52,6 +54,8 @@ bool operator==(chopper::configuration const & lhs, chopper::configuration const lhs.window_size == rhs.window_size && // lhs.disable_sketch_output == rhs.disable_sketch_output && // lhs.precomputed_files == rhs.precomputed_files && // + lhs.number_of_partitions == rhs.number_of_partitions && // + lhs.partitioning_approach == rhs.partitioning_approach && // lhs.output_filename == rhs.output_filename && // lhs.determine_best_tmax == rhs.determine_best_tmax && // lhs.force_all_binnings == rhs.force_all_binnings && // @@ -74,7 +78,7 @@ bool operator==(chopper::configuration const & lhs, chopper::configuration const static constexpr std::string_view config_string_view{"@CHOPPER_CONFIG\n" "@{\n" "@ \"chopper_config\": {\n" - "@ \"version\": 3,\n" + "@ \"version\": 4,\n" "@ \"data_file\": {\n" "@ \"value0\": \"/path/to/data.file\"\n" "@ },\n" @@ -91,7 +95,9 @@ static constexpr std::string_view config_string_view{"@CHOPPER_CONFIG\n" "@ },\n" "@ \"determine_best_tmax\": true,\n" "@ \"force_all_binnings\": true,\n" - "@ \"fast_layout\": true\n" + "@ \"fast_layout\": true,\n" + "@ \"number_of_partitions\": 4,\n" + "@ \"partitioning_approach\": 6\n" "@ }\n" "@}\n" "@CHOPPER_CONFIG_END\n" @@ -165,8 +171,12 @@ TEST(config_test, read_from_version_2) config_string.replace(pos, from.size(), to); }; - replace("\"version\": 3,\n@ \"data_file\"", "\"version\": 2,\n@ \"data_file\""); - replace("\"force_all_binnings\": true,\n@ \"fast_layout\": true\n", "\"force_all_binnings\": true\n"); + replace("\"version\": 4,\n@ \"data_file\"", "\"version\": 2,\n@ \"data_file\""); + replace("\"force_all_binnings\": true,\n" + "@ \"fast_layout\": true,\n" + "@ \"number_of_partitions\": 4,\n" + "@ \"partitioning_approach\": 6\n", + "\"force_all_binnings\": true\n"); std::stringstream ss{config_string}; @@ -175,41 +185,40 @@ TEST(config_test, read_from_version_2) chopper::configuration expected{generate_config()}; expected.fast_layout = false; + expected.number_of_partitions = 1; + expected.partitioning_approach = chopper::partitioning_scheme::blocked; EXPECT_EQ(config, expected); } -// Easier to do in the config_test because of existing helper functions -TEST(input, read_layout_file) +// Configs written before version 4 do not contain the partitioned HIBF fields. +TEST(config_test, read_from_version_3) { - std::string config_string{"@CHOPPER_USER_BINS\n" - "@0 file1.fa\n" - "@1 file2.fa\n" - "@2 path/to/file3.fa\n" - "@3 file4.fastq\n" - "@CHOPPER_USER_BINS_END\n"}; - config_string += config_string_view; + std::string config_string{config_string_view}; + auto erase_line = [&config_string](std::string_view const line) { - seqan::hibf::layout::layout hibf_layout{}; - hibf_layout.top_level_max_bin_id = 42u; - std::stringstream config_out{}; - hibf_layout.write_to(config_out); - config_string += config_out.str(); - } + size_t const pos = config_string.find(line); + ASSERT_NE(pos, std::string::npos) << line; + config_string.erase(pos, line.size()); + }; + erase_line(",\n@ \"number_of_partitions\": 4"); + erase_line(",\n@ \"partitioning_approach\": 6"); + + std::string_view const version_4{"@ \"version\": 4,\n"}; + size_t const pos = config_string.find(version_4); + ASSERT_NE(pos, std::string::npos); + config_string.replace(pos, version_4.size(), "@ \"version\": 3,\n"); std::stringstream ss{config_string}; - auto [filenames, config, layout] = chopper::layout::read_layout_file(ss); + chopper::configuration config; + config.read_from(ss); - std::vector> const expected_filenames{{"file1.fa"}, - {"file2.fa"}, - {"path/to/file3.fa"}, - {"file4.fastq"}}; + // The fields are not read and keep their default values. + chopper::configuration expected{generate_config()}; + expected.number_of_partitions = 1; + expected.partitioning_approach = chopper::partitioning_scheme::blocked; - EXPECT_EQ(filenames, expected_filenames); - EXPECT_EQ(config, generate_config()); - EXPECT_EQ(layout.top_level_max_bin_id, 42u); - EXPECT_TRUE(layout.max_bins.empty()); - EXPECT_TRUE(layout.user_bins.empty()); + EXPECT_EQ(config, expected); } diff --git a/test/api/layout/CMakeLists.txt b/test/api/layout/CMakeLists.txt index fda6d2b4..10e94b45 100644 --- a/test/api/layout/CMakeLists.txt +++ b/test/api/layout/CMakeLists.txt @@ -31,3 +31,5 @@ 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 (lsh_distributed_ibf_layout_test.cpp) +add_api_test (phibf_multi_cluster_test.cpp) +add_api_test (phibf_test.cpp) diff --git a/test/api/layout/execute_layout_test.cpp b/test/api/layout/execute_layout_test.cpp index 187acd61..87508a42 100644 --- a/test/api/layout/execute_layout_test.cpp +++ b/test/api/layout/execute_layout_test.cpp @@ -67,7 +67,7 @@ TEST(execute_test, few_ubs) "@CHOPPER_CONFIG\n" "@{\n" "@ \"chopper_config\": {\n" - "@ \"version\": 3,\n" + "@ \"version\": 4,\n" "@ \"data_file\": {\n" "@ \"value0\": \"\"\n" "@ },\n" @@ -86,7 +86,9 @@ TEST(execute_test, few_ubs) "@ },\n" "@ \"determine_best_tmax\": false,\n" "@ \"force_all_binnings\": false,\n" - "@ \"fast_layout\": false\n" + "@ \"fast_layout\": false,\n" + "@ \"number_of_partitions\": 1,\n" + "@ \"partitioning_approach\": 0\n" "@ }\n" "@}\n" "@CHOPPER_CONFIG_END\n" @@ -178,7 +180,7 @@ TEST(execute_test, few_ubs_fast_layout) "@CHOPPER_CONFIG\n" "@{\n" "@ \"chopper_config\": {\n" - "@ \"version\": 3,\n" + "@ \"version\": 4,\n" "@ \"data_file\": {\n" "@ \"value0\": \"\"\n" "@ },\n" @@ -197,7 +199,9 @@ TEST(execute_test, few_ubs_fast_layout) "@ },\n" "@ \"determine_best_tmax\": false,\n" "@ \"force_all_binnings\": false,\n" - "@ \"fast_layout\": true\n" + "@ \"fast_layout\": true,\n" + "@ \"number_of_partitions\": 1,\n" + "@ \"partitioning_approach\": 0\n" "@ }\n" "@}\n" "@CHOPPER_CONFIG_END\n" @@ -408,7 +412,7 @@ TEST(execute_test, many_ubs) "@CHOPPER_CONFIG\n" "@{\n" "@ \"chopper_config\": {\n" - "@ \"version\": 3,\n" + "@ \"version\": 4,\n" "@ \"data_file\": {\n" "@ \"value0\": \"\"\n" "@ },\n" @@ -427,7 +431,9 @@ TEST(execute_test, many_ubs) "@ },\n" "@ \"determine_best_tmax\": false,\n" "@ \"force_all_binnings\": false,\n" - "@ \"fast_layout\": false\n" + "@ \"fast_layout\": false,\n" + "@ \"number_of_partitions\": 1,\n" + "@ \"partitioning_approach\": 0\n" "@ }\n" "@}\n" "@CHOPPER_CONFIG_END\n" @@ -703,7 +709,7 @@ TEST(execute_test, many_ubs_fast_layout) "@CHOPPER_CONFIG\n" "@{\n" "@ \"chopper_config\": {\n" - "@ \"version\": 3,\n" + "@ \"version\": 4,\n" "@ \"data_file\": {\n" "@ \"value0\": \"\"\n" "@ },\n" @@ -722,7 +728,9 @@ TEST(execute_test, many_ubs_fast_layout) "@ },\n" "@ \"determine_best_tmax\": false,\n" "@ \"force_all_binnings\": false,\n" - "@ \"fast_layout\": true\n" + "@ \"fast_layout\": true,\n" + "@ \"number_of_partitions\": 1,\n" + "@ \"partitioning_approach\": 0\n" "@ }\n" "@}\n" "@CHOPPER_CONFIG_END\n" diff --git a/test/api/layout/execute_with_estimation_test.cpp b/test/api/layout/execute_with_estimation_test.cpp index ef679d5f..ca46cf89 100644 --- a/test/api/layout/execute_with_estimation_test.cpp +++ b/test/api/layout/execute_with_estimation_test.cpp @@ -346,7 +346,7 @@ TEST(execute_estimation_test, many_ubs) "@CHOPPER_CONFIG\n" "@{\n" "@ \"chopper_config\": {\n" - "@ \"version\": 3,\n" + "@ \"version\": 4,\n" "@ \"data_file\": {\n" "@ \"value0\": \"\"\n" "@ },\n" @@ -365,7 +365,9 @@ TEST(execute_estimation_test, many_ubs) "@ },\n" "@ \"determine_best_tmax\": true,\n" "@ \"force_all_binnings\": false,\n" - "@ \"fast_layout\": false\n" + "@ \"fast_layout\": false,\n" + "@ \"number_of_partitions\": 1,\n" + "@ \"partitioning_approach\": 0\n" "@ }\n" "@}\n" "@CHOPPER_CONFIG_END\n" @@ -666,7 +668,7 @@ TEST(execute_estimation_test, many_ubs_fast_layout) "@CHOPPER_CONFIG\n" "@{\n" "@ \"chopper_config\": {\n" - "@ \"version\": 3,\n" + "@ \"version\": 4,\n" "@ \"data_file\": {\n" "@ \"value0\": \"\"\n" "@ },\n" @@ -685,7 +687,9 @@ TEST(execute_estimation_test, many_ubs_fast_layout) "@ },\n" "@ \"determine_best_tmax\": true,\n" "@ \"force_all_binnings\": false,\n" - "@ \"fast_layout\": true\n" + "@ \"fast_layout\": true,\n" + "@ \"number_of_partitions\": 1,\n" + "@ \"partitioning_approach\": 0\n" "@ }\n" "@}\n" "@CHOPPER_CONFIG_END\n" diff --git a/test/api/layout/input_test.cpp b/test/api/layout/input_test.cpp index afde1222..fd4a6ddb 100644 --- a/test/api/layout/input_test.cpp +++ b/test/api/layout/input_test.cpp @@ -6,81 +6,164 @@ #include -std::string const layout_header{"@CHOPPER_USER_BINS\n" - "@0 seq0.fa\n" - "@1 seq1.fa\n" - "@2 seq2.fa\n" - "@3 seq3.fa\n" - "@4 seq4.fa\n" - "@5 seq5.fa\n" - "@6 seq6.fa\n" - "@7 seq7.fa\n" - "@CHOPPER_USER_BINS_END\n" - "@CHOPPER_CONFIG\n" - "@{\n" - "@ \"chopper_config\": {\n" - "@ \"version\": 3,\n" - "@ \"data_file\": {\n" - "@ \"value0\": \"input.txt\"\n" - "@ },\n" - "@ \"debug\": false,\n" - "@ \"sketch_directory\": {\n" - "@ \"value0\": \"\"\n" - "@ },\n" - "@ \"k\": 15,\n" - "@ \"window_size\": 15,\n" - "@ \"disable_sketch_output\": true,\n" - "@ \"precomputed_files\": false,\n" - "@ \"output_filename\": {\n" - "@ \"value0\": \"foo.layout\"\n" - "@ },\n" - "@ \"determine_best_tmax\": false,\n" - "@ \"force_all_binnings\": false,\n" - "@ \"fast_layout\": false\n" - "@ }\n" - "@}\n" - "@CHOPPER_CONFIG_END\n" - "@HIBF_CONFIG\n" - "@{\n" - "@ \"hibf_config\": {\n" - "@ \"version\": 3,\n" - "@ \"number_of_user_bins\": 3,\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"}; - -TEST(layout_test, read_single_layout) +// A layout file consists of one section per (partitioned HIBF) layout. Each section has its own user bin filenames, +// chopper config, and layout. User bin indices within a section are local to that section. +std::string layout_section(std::vector const & filenames, + size_t const number_of_partitions, + std::string const & layout) { - // layout consists of three partitions, written one after the other - std::stringstream ss{layout_header + R"layout_file(#TOP_LEVEL_IBF fullest_technical_bin_idx:111 + std::string section{"@CHOPPER_USER_BINS\n"}; + for (size_t i = 0; i < filenames.size(); ++i) + section += "@" + std::to_string(i) + " " + filenames[i] + "\n"; + + section += "@CHOPPER_USER_BINS_END\n" + "@CHOPPER_CONFIG\n" + "@{\n" + "@ \"chopper_config\": {\n" + "@ \"version\": 4,\n" + "@ \"data_file\": {\n" + "@ \"value0\": \"input.txt\"\n" + "@ },\n" + "@ \"debug\": false,\n" + "@ \"sketch_directory\": {\n" + "@ \"value0\": \"\"\n" + "@ },\n" + "@ \"k\": 15,\n" + "@ \"window_size\": 15,\n" + "@ \"disable_sketch_output\": true,\n" + "@ \"precomputed_files\": false,\n" + "@ \"output_filename\": {\n" + "@ \"value0\": \"foo.layout\"\n" + "@ },\n" + "@ \"determine_best_tmax\": false,\n" + "@ \"force_all_binnings\": false,\n" + "@ \"fast_layout\": false,\n" + "@ \"number_of_partitions\": " + + std::to_string(number_of_partitions) + + ",\n" + "@ \"partitioning_approach\": 1\n" + "@ }\n" + "@}\n" + "@CHOPPER_CONFIG_END\n" + "@HIBF_CONFIG\n" + "@{\n" + "@ \"hibf_config\": {\n" + "@ \"version\": 3,\n" + "@ \"number_of_user_bins\": " + + std::to_string(filenames.size()) + + ",\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" + + layout; + return section; +} + +std::string const layout_with_three_user_bins{R"layout_file(#TOP_LEVEL_IBF fullest_technical_bin_idx:111 #LOWER_LEVEL_IBF_0 fullest_technical_bin_idx:0 #LOWER_LEVEL_IBF_2 fullest_technical_bin_idx:2 #LOWER_LEVEL_IBF_1;2;3;4 fullest_technical_bin_idx:22 #USER_BIN_IDX TECHNICAL_BIN_INDICES NUMBER_OF_TECHNICAL_BINS -7 0 1 -4 1;0 1;22 -5 1;2;3;4;22 1;1;1;1;21 +2 0 1 +0 1;0 1;22 +1 1;2;3;4;22 1;1;1;1;21 )layout_file"}; - auto [filenames, chopper_config, layout] = chopper::layout::read_layout_file(ss); +std::string const layout_with_two_user_bins{R"layout_file(#TOP_LEVEL_IBF fullest_technical_bin_idx:0 +#USER_BIN_IDX TECHNICAL_BIN_INDICES NUMBER_OF_TECHNICAL_BINS +1 0 40 +0 40 24 +)layout_file"}; +void check_layout_with_three_user_bins(seqan::hibf::layout::layout const & layout) +{ EXPECT_EQ(layout.top_level_max_bin_id, 111); + ASSERT_EQ(layout.max_bins.size(), 3u); EXPECT_EQ(layout.max_bins[0], (seqan::hibf::layout::layout::max_bin{{0}, 0})); EXPECT_EQ(layout.max_bins[1], (seqan::hibf::layout::layout::max_bin{{2}, 2})); EXPECT_EQ(layout.max_bins[2], (seqan::hibf::layout::layout::max_bin{{1, 2, 3, 4}, 22})); - EXPECT_EQ(layout.user_bins[0], (seqan::hibf::layout::layout::user_bin{std::vector{}, 0, 1, 7})); - EXPECT_EQ(layout.user_bins[1], (seqan::hibf::layout::layout::user_bin{std::vector{1}, 0, 22, 4})); - EXPECT_EQ(layout.user_bins[2], (seqan::hibf::layout::layout::user_bin{std::vector{1, 2, 3, 4}, 22, 21, 5})); + ASSERT_EQ(layout.user_bins.size(), 3u); + EXPECT_EQ(layout.user_bins[0], (seqan::hibf::layout::layout::user_bin{std::vector{}, 0, 1, 2})); + EXPECT_EQ(layout.user_bins[1], (seqan::hibf::layout::layout::user_bin{std::vector{1}, 0, 22, 0})); + EXPECT_EQ(layout.user_bins[2], (seqan::hibf::layout::layout::user_bin{std::vector{1, 2, 3, 4}, 22, 21, 1})); +} + +void check_layout_with_two_user_bins(seqan::hibf::layout::layout const & layout) +{ + EXPECT_EQ(layout.top_level_max_bin_id, 0); + EXPECT_TRUE(layout.max_bins.empty()); + ASSERT_EQ(layout.user_bins.size(), 2u); + EXPECT_EQ(layout.user_bins[0], (seqan::hibf::layout::layout::user_bin{std::vector{}, 0, 40, 1})); + EXPECT_EQ(layout.user_bins[1], (seqan::hibf::layout::layout::user_bin{std::vector{}, 40, 24, 0})); +} + +TEST(layout_test, read_single_layout) +{ + std::stringstream ss{layout_section({"seq0.fa", "seq1.fa", "seq2.fa"}, 1u, layout_with_three_user_bins)}; + + auto [filenames, chopper_config, layout] = chopper::layout::read_layout_file(ss); + + EXPECT_EQ(filenames, (std::vector>{{"seq0.fa"}, {"seq1.fa"}, {"seq2.fa"}})); + EXPECT_EQ(chopper_config.k, 15u); + EXPECT_EQ(chopper_config.number_of_partitions, 1u); + EXPECT_EQ(chopper_config.hibf_config.number_of_user_bins, 3u); + check_layout_with_three_user_bins(layout); + + // The whole file has been consumed. + EXPECT_FALSE(ss.good()); +} + +TEST(layout_test, read_from_partitioned_layout) +{ + // The layout consists of three partitions, each written as its own section, one after the other. + std::stringstream ss{layout_section({"seq0.fa", "seq1.fa", "seq2.fa"}, 3u, layout_with_three_user_bins) + + layout_section({"seq3.fa", "seq4.fa"}, 3u, layout_with_two_user_bins) + + layout_section({"seq5.fa", "seq6.fa", "seq7.fa"}, 3u, layout_with_three_user_bins)}; + + std::vector>> filenames{}; + std::vector configs{}; + std::vector layouts{}; + + while (ss.good()) + { + auto [section_filenames, section_config, section_layout] = chopper::layout::read_layout_file(ss); + filenames.push_back(std::move(section_filenames)); + configs.push_back(std::move(section_config)); + layouts.push_back(std::move(section_layout)); + } + + ASSERT_EQ(layouts.size(), 3u); + + // Each section has its own filenames. + EXPECT_EQ(filenames[0], (std::vector>{{"seq0.fa"}, {"seq1.fa"}, {"seq2.fa"}})); + EXPECT_EQ(filenames[1], (std::vector>{{"seq3.fa"}, {"seq4.fa"}})); + EXPECT_EQ(filenames[2], (std::vector>{{"seq5.fa"}, {"seq6.fa"}, {"seq7.fa"}})); + + // Each section has its own config, whose number of user bins is the one of its partition. + for (auto const & config : configs) + { + EXPECT_EQ(config.k, 15u); + EXPECT_EQ(config.number_of_partitions, 3u); + EXPECT_EQ(config.partitioning_approach, chopper::partitioning_scheme::sorted); + } + EXPECT_EQ(configs[0].hibf_config.number_of_user_bins, 3u); + EXPECT_EQ(configs[1].hibf_config.number_of_user_bins, 2u); + EXPECT_EQ(configs[2].hibf_config.number_of_user_bins, 3u); + + // Each layout keeps the user bin indices local to its section. + check_layout_with_three_user_bins(layouts[0]); + check_layout_with_two_user_bins(layouts[1]); + check_layout_with_three_user_bins(layouts[2]); } diff --git a/test/api/layout/phibf_multi_cluster_test.cpp b/test/api/layout/phibf_multi_cluster_test.cpp new file mode 100644 index 00000000..0e609e1d --- /dev/null +++ b/test/api/layout/phibf_multi_cluster_test.cpp @@ -0,0 +1,130 @@ +#include // for Test, TestInfo, EXPECT_EQ, Message, TEST, TestPartResult + +#include // for size_t +#include // for vector + +#include +#include + +using chopper::layout::Cluster; +using chopper::layout::MultiCluster; + +TEST(Multicluster_test, ctor_from_cluster) +{ + Cluster const cluster1{5}; + MultiCluster const multi_cluster1{cluster1}; + + EXPECT_EQ(multi_cluster1.id(), cluster1.id()); + EXPECT_FALSE(multi_cluster1.empty()); + EXPECT_EQ(multi_cluster1.size(), 1u); + ASSERT_EQ(multi_cluster1.contained_user_bins().size(), 1u); + ASSERT_EQ(multi_cluster1.contained_user_bins()[0].size(), 1u); + EXPECT_EQ(multi_cluster1.contained_user_bins()[0][0], 5u); + EXPECT_TRUE(multi_cluster1.is_valid(cluster1.id())); +} + +TEST(Multicluster_test, ctor_from_moved_cluster) +{ + Cluster cluster1{5}; + Cluster cluster2{7}; + cluster2.move_to(cluster1); + + MultiCluster const multi_cluster2{cluster2}; + + EXPECT_EQ(multi_cluster2.id(), cluster2.id()); + EXPECT_TRUE(multi_cluster2.empty()); + EXPECT_TRUE(multi_cluster2.has_been_moved()); + EXPECT_EQ(multi_cluster2.moved_to_cluster_id(), cluster2.moved_to_cluster_id()); + EXPECT_EQ(multi_cluster2.size(), 0u); + ASSERT_EQ(multi_cluster2.contained_user_bins().size(), 0u); + EXPECT_TRUE(multi_cluster2.is_valid(cluster2.id())); +} + +TEST(Multicluster_test, move_to) +{ + Cluster cluster1{5}; + Cluster cluster2{7}; + cluster2.move_to(cluster1); + ASSERT_EQ(cluster1.size(), 2u); + Cluster const cluster3{13}; + + MultiCluster multi_cluster1{cluster1}; + EXPECT_TRUE(multi_cluster1.is_valid(cluster1.id())); + EXPECT_EQ(multi_cluster1.size(), 1u); + EXPECT_EQ(multi_cluster1.contained_user_bins().size(), 1u); + EXPECT_EQ(multi_cluster1.contained_user_bins()[0].size(), 2u); + + MultiCluster multi_cluster3{cluster3}; + EXPECT_TRUE(multi_cluster3.is_valid(cluster3.id())); + EXPECT_EQ(multi_cluster3.size(), 1u); + + multi_cluster1.move_to(multi_cluster3); + + EXPECT_TRUE(multi_cluster1.is_valid(cluster1.id())); + EXPECT_TRUE(multi_cluster3.is_valid(cluster3.id())); + + // multi_cluster1 has been moved and is empty now + EXPECT_TRUE(multi_cluster1.empty()); + EXPECT_EQ(multi_cluster1.size(), 0u); + EXPECT_TRUE(multi_cluster1.has_been_moved()); + EXPECT_EQ(multi_cluster1.moved_to_cluster_id(), multi_cluster3.id()); + + // multi_cluster3 contains 2 clusters now, {13} and {5, 7} + EXPECT_FALSE(multi_cluster3.empty()); + EXPECT_EQ(multi_cluster3.size(), 2u); // two clusters + ASSERT_EQ(multi_cluster3.contained_user_bins().size(), 2u); + ASSERT_EQ(multi_cluster3.contained_user_bins()[0].size(), 1u); + ASSERT_EQ(multi_cluster3.contained_user_bins()[1].size(), 2u); + EXPECT_EQ(multi_cluster3.contained_user_bins()[0][0], cluster3.id()); + EXPECT_EQ(multi_cluster3.contained_user_bins()[1][0], cluster1.id()); + EXPECT_EQ(multi_cluster3.contained_user_bins()[1][1], cluster2.id()); +} + +TEST(Multicluster_test, sort_by_cardinality) +{ + std::vector const cardinalities{10, 50, 20, 40, 30}; + + Cluster cluster0{0}; + Cluster cluster1{1}; + Cluster cluster2{2}; + Cluster cluster3{3}; + cluster1.move_to(cluster0); // {0, 1} + cluster3.move_to(cluster2); // {2, 3} + Cluster cluster4{4}; + + MultiCluster multi_cluster0{cluster0}; // {{0, 1}} + MultiCluster multi_cluster4{cluster4}; // {{4}} + MultiCluster multi_cluster2{cluster2}; // {{2, 3}} + multi_cluster0.move_to(multi_cluster4); // {{4}, {0, 1}} + multi_cluster2.move_to(multi_cluster4); // {{4}, {0, 1}, {2, 3}} + + multi_cluster4.sort_by_cardinality(cardinalities); + + auto const & clusters = multi_cluster4.contained_user_bins(); + ASSERT_EQ(clusters.size(), 3u); + // inner clusters are sorted by descending cardinality, clusters by descending size + EXPECT_EQ(clusters[2], (std::vector{4})); + for (size_t i = 0; i < 2; ++i) + { + ASSERT_EQ(clusters[i].size(), 2u); + EXPECT_GE(cardinalities[clusters[i][0]], cardinalities[clusters[i][1]]); + } +} + +TEST(LSH_find_representative_cluster_test, multi_cluster_one_move) +{ + std::vector mclusters{{Cluster{0}}, {Cluster{1}}}; + mclusters[1].move_to(mclusters[0]); + + EXPECT_EQ(chopper::layout::LSH_find_representative_cluster(mclusters, mclusters[1].id()), mclusters[0].id()); +} + +TEST(LSH_find_representative_cluster_test, multi_cluster_two_moves) +{ + std::vector mclusters{{Cluster{0}}, {Cluster{1}}, {Cluster{2}}}; + mclusters[2].move_to(mclusters[1]); + mclusters[1].move_to(mclusters[0]); + + EXPECT_EQ(chopper::layout::LSH_find_representative_cluster(mclusters, mclusters[1].id()), mclusters[0].id()); + EXPECT_EQ(chopper::layout::LSH_find_representative_cluster(mclusters, mclusters[2].id()), mclusters[0].id()); +} diff --git a/test/api/layout/phibf_test.cpp b/test/api/layout/phibf_test.cpp new file mode 100644 index 00000000..2696a86b --- /dev/null +++ b/test/api/layout/phibf_test.cpp @@ -0,0 +1,481 @@ +// -------------------------------------------------------------------------------------------------- +// 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 +#include // for size_t +#include // for uint64_t +#include // for ifstream +#include // for iota +#include // for invalid_argument +#include // for string +#include // for vector + +#include +#include +#include +#include + +#include +#include +#include + +#include "../api_test.hpp" + +// Not declared in the header. Tested directly in phibf_regression_test.find_best_partition_intersection_underflow. +namespace chopper::layout +{ +void find_best_partition(chopper::configuration const & config, + size_t const cardinality_per_part, + std::vector const & cluster, + std::vector const & cardinalities, + std::vector const & sketches, + std::vector> & positions, + std::vector & partition_sketches, + std::vector & partition_cardinality); +} // namespace chopper::layout + +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); +} + +// `number_of_user_bins` user bins in families of 5 that mostly overlap, with sketches and the configuration for +// partition_user_bins. dist: "uniform" (3'000 to 23'000 hashes), "equal" (5'000 hashes), or "skewed" (two user bins +// with 400'000 hashes, the rest 3'000 to 3'600). +struct phibf_data +{ + std::vector kmer_counts{}; + chopper::configuration config{}; + std::vector sketches{}; + std::vector minHash_sketches{}; + std::vector cardinalities{}; + + phibf_data(chopper::partitioning_scheme const partitioning_approach, + size_t const number_of_user_bins, + size_t const number_of_partitions, + std::string const & dist = "uniform") + { + for (size_t i = 0; i < number_of_user_bins; ++i) + { + if (dist == "equal") + kmer_counts.push_back(5000); + else if (dist == "skewed") + kmer_counts.push_back((i < 2) ? 400'000 : 3000 + 100 * (i % 7)); + else + kmer_counts.push_back(3000 + (scramble(i) % 20000)); + } + + config.number_of_partitions = number_of_partitions; + config.partitioning_approach = partitioning_approach; + config.hibf_config.number_of_user_bins = number_of_user_bins; + config.hibf_config.tmax = 64; + config.hibf_config.input_fn = [this](size_t const ub, seqan::hibf::insert_iterator it) + { + uint64_t const offset = static_cast(ub / 5) << 32; + for (uint64_t i = 0; i < kmer_counts[ub]; ++i) + it = scramble(offset + i + (ub % 5) * 97); + }; + config.hibf_config.validate_and_set_defaults(); + + seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches); + + for (auto const & sketch : sketches) + cardinalities.push_back(sketch.estimate()); + } + + std::vector> partition() const + { + std::vector> partitions(config.number_of_partitions); + chopper::layout::partition_user_bins(config, cardinalities, sketches, minHash_sketches, partitions); + return partitions; + } + + std::vector layouts() + { + seqan3::test::tmp_directory tmp_dir{}; + config.output_filename = tmp_dir.path() / "dummy.layout"; + std::vector> filenames(kmer_counts.size(), {"ub.fa"}); + chopper::layout::execute(config, filenames, sketches, minHash_sketches); + + std::vector layouts{}; + std::ifstream layout_stream{config.output_filename}; + while (layout_stream.good()) + { + layouts.push_back(std::get<2>(chopper::layout::read_layout_file(layout_stream))); + } + return layouts; + } +}; + +std::vector const phibf_approaches{chopper::partitioning_scheme::blocked, + chopper::partitioning_scheme::sorted, + chopper::partitioning_scheme::folded, + chopper::partitioning_scheme::weighted_fold, + chopper::partitioning_scheme::similarity, + chopper::partitioning_scheme::lsh, + chopper::partitioning_scheme::lsh_sim}; + +// Expects that every user bin in [0, number_of_user_bins) is assigned to exactly one partition. +bool each_user_bin_assigned_once(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); + + std::vector expected_user_bins(number_of_user_bins); + std::iota(expected_user_bins.begin(), expected_user_bins.end(), 0u); + return assigned_user_bins == expected_user_bins; +} + +} // namespace + +TEST(phibf_partition_test, each_user_bin_is_assigned_to_exactly_one_partition) +{ + size_t const number_of_user_bins{240}; + size_t const number_of_partitions{4}; + + for (chopper::partitioning_scheme const approach : phibf_approaches) + { + phibf_data data{approach, number_of_user_bins, number_of_partitions}; + auto partitions = data.partition(); + ASSERT_EQ(partitions.size(), number_of_partitions); + EXPECT_TRUE(each_user_bin_assigned_once(partitions, number_of_user_bins)); + } +} + +TEST(phibf_partition_test, unknown_partitioning_approach) +{ + phibf_data data{chopper::partitioning_scheme::unknown, 16, 4}; + EXPECT_THROW(data.partition(), std::invalid_argument); +} + +TEST(phibf_execute_test, writes_one_layout_per_partition) +{ + seqan3::test::tmp_directory tmp_dir{}; + std::filesystem::path const layout_file{tmp_dir.path() / "phibf.layout"}; + + size_t const number_of_user_bins{240}; + size_t const number_of_partitions{3}; + phibf_data data{chopper::partitioning_scheme::lsh_sim, number_of_user_bins, number_of_partitions}; + data.config.output_filename = layout_file; + + std::vector> filenames{}; + for (size_t ub = 0; ub < number_of_user_bins; ++ub) + filenames.push_back({"ub" + std::to_string(ub) + ".fa"}); + + chopper::layout::execute(data.config, filenames, data.sketches, data.minHash_sketches); + + std::ifstream layout_stream{layout_file}; + + size_t total_read_number_of_user_bins{0}; + size_t number_of_layouts{0}; + std::vector> total_read_filenames{}; + while (layout_stream.good()) + { + auto [read_filenames, read_config, layout] = chopper::layout::read_layout_file(layout_stream); + total_read_number_of_user_bins += layout.user_bins.size(); + ++number_of_layouts; + total_read_filenames.insert(total_read_filenames.end(), read_filenames.begin(), read_filenames.end()); + + EXPECT_EQ(read_config.partitioning_approach, chopper::partitioning_scheme::lsh_sim); + EXPECT_EQ(read_config.number_of_partitions, number_of_partitions); + EXPECT_EQ(read_config.hibf_config.number_of_user_bins, layout.user_bins.size()); + + std::vector> user_bins_single_parition(1); // enables us to use each_user_bin_assigned_once + for (auto const & user_bin : layout.user_bins) + user_bins_single_parition.back().push_back(user_bin.idx); + EXPECT_TRUE(each_user_bin_assigned_once(user_bins_single_parition, layout.user_bins.size())); + } + EXPECT_EQ(total_read_number_of_user_bins, number_of_user_bins); + EXPECT_EQ(number_of_layouts, number_of_partitions); + std::ranges::sort(total_read_filenames); + std::ranges::sort(filenames); + EXPECT_EQ(total_read_filenames, filenames); +} + +TEST(phibf_invalid_input_test, determine_best_tmax_is_set) +{ + phibf_data data{chopper::partitioning_scheme::sorted, 240, 2}; + data.config.determine_best_tmax = true; + EXPECT_THROW(data.layouts(), std::invalid_argument); +} + +TEST(phibf_invalid_input_test, fewer_user_bins_than_partitions) +{ + for (chopper::partitioning_scheme const approach : phibf_approaches) + { + phibf_data const data{approach, 1 /*UB*/, 2 /*partitions*/, "uniform"}; + EXPECT_THROW(data.partition(), std::invalid_argument) << "approach " << approach; + } +} + +TEST(phibf_regression_test, small_number_of_user_bins_per_part) +{ + // If for small number of user bins the block size could not be optimally chosen (runding errors) + // The last partition ended up empty. This was fixed. + phibf_data const data{chopper::partitioning_scheme::blocked, 100, 16, "uniform"}; + EXPECT_TRUE(each_user_bin_assigned_once(data.partition(), 100)); +} + +TEST(phibf_regression_test, few_user_bins_with_skewed_sizes) +{ + // If few user bins with skewed sizes are the input, + // then cardinality_per_part is chosen too high and partitions at the end get no user bins. + // This has beenf fixed. + phibf_data const data{chopper::partitioning_scheme::sorted, 64, 8, "skewed"}; + EXPECT_TRUE(each_user_bin_assigned_once(data.partition(), 64)); +} + +TEST(phibf_regression_test, zero_cardinality_user_bins) +{ + using chopper::partitioning_scheme; + + for (chopper::partitioning_scheme const approach : phibf_approaches) + { + phibf_data data{approach, 5, 2, "equal"}; + data.cardinalities = {10, 10, 10, 10, 0}; + + auto const partitions = data.partition(); + EXPECT_TRUE(each_user_bin_assigned_once(partitions, 5)); + } +} + +TEST(phibf_regression_test, fewer_clusters_than_partitions) +{ + using chopper::partitioning_scheme; + + for (chopper::partitioning_scheme const approach : {partitioning_scheme::lsh, partitioning_scheme::lsh_sim}) + { + phibf_data const data{approach, 9, 8, "equal"}; + EXPECT_TRUE(each_user_bin_assigned_once(data.partition(), 9)); + } +} + +TEST(phibf_regression_test, user_bin_does_not_fit) +{ + // Those inputs made the respective partitioning approach fail, because one user bin could not be assigned + // To a partition within find_best_partition. For example single very large user bins or very skewed distributions. + // find_best_partition now has a fallback approach. + phibf_data const data_lsh{chopper::partitioning_scheme::lsh, 8, 5, "equal"}; + EXPECT_TRUE(each_user_bin_assigned_once(data_lsh.partition(), 8)); + + phibf_data const data_lsh_sim{chopper::partitioning_scheme::lsh_sim, 257, 8, "skewed"}; + EXPECT_TRUE(each_user_bin_assigned_once(data_lsh_sim.partition(), 257)); + + phibf_data const data_similarity{chopper::partitioning_scheme::similarity, 4, 3, "equal"}; + EXPECT_TRUE(each_user_bin_assigned_once(data_similarity.partition(), 4)); +} + +TEST(phibf_regression_test, weighted_fold_threshold) +{ + // About two user bins per partition: some partitions stay empty, which is reported. + // Formerly, u_bins_per_part was rounded up, causing some partitions to use up to many user bins. + // Now, it is rounded down and the following edge case passes: + phibf_data const data3{chopper::partitioning_scheme::weighted_fold, 33, 16, "uniform"}; + EXPECT_TRUE(each_user_bin_assigned_once(data3.partition(), 33)); +} + +TEST(phibf_regression_test, weighted_fold_score_edge_case) +{ + // 8 user bins, 2 partitions: 125 k-mers and 4 user bins per partition. + phibf_data data{chopper::partitioning_scheme::weighted_fold, 8, 2, "equal"}; + data.cardinalities = {100, 90, 10, 10, 10, 10, 10, 10}; + + // Partition 0 takes the large user bins 100 and 90 (score |1 - 190/125| + |1 - 2/4| = 1.02). + // Removing 90 and adding small user bins while the score improves: 110/2 (0.62), 120/3 (0.29), 130/4 (0.04); + // a fourth one would be worse (140/5: 0.37). 0.04 < 1.02, so this is accepted. + // Removing 100 as well cannot improve on that (at best 30/4: 0.76). The rest goes to partition 1. + auto const partitions = data.partition(); + ASSERT_EQ(partitions.size(), 2u); + + auto partition_cardinalities = [&data](std::vector const & partition) + { + std::vector result; + for (size_t const user_bin : partition) + result.push_back(data.cardinalities[user_bin]); + std::ranges::sort(result, std::ranges::greater{}); + return result; + }; + + EXPECT_EQ(partition_cardinalities(partitions[0]), (std::vector{100, 10, 10, 10})); + EXPECT_EQ(partition_cardinalities(partitions[1]), (std::vector{90, 10, 10, 10})); + EXPECT_TRUE(each_user_bin_assigned_once(partitions, 8)); +} + +TEST(phibf_execute_test, tmaxset_or_not_set) +{ + // The number of technical bins of the top-level IBF of a layout. + auto top_level_technical_bins = [](seqan::hibf::layout::layout const & layout) + { + size_t result{}; + for (auto const & user_bin : layout.user_bins) + { + if (user_bin.previous_TB_indices.empty()) // single or split bin on the top level + result = std::max(result, user_bin.storage_TB_id + user_bin.number_of_technical_bins); + else // merged bin on the top level + result = std::max(result, user_bin.previous_TB_indices[0] + 1); + } + return result; + }; + + size_t const number_of_partitions{3}; + + // tmax not set: about 80 user bins per partition, next_multiple_of_64(ceil(sqrt(80))) = 64. + { + phibf_data data{chopper::partitioning_scheme::sorted, 240, number_of_partitions}; + auto const layouts = data.layouts(); + ASSERT_EQ(layouts.size(), number_of_partitions); + for (auto const & layout : layouts) + EXPECT_EQ(top_level_technical_bins(layout), 64u); + } + + // tmax set: every partition uses it. + { + phibf_data data{chopper::partitioning_scheme::sorted, 240, number_of_partitions}; + data.config.hibf_config.tmax = 128; + data.config.tmax_is_set = true; + auto const layouts = data.layouts(); + ASSERT_EQ(layouts.size(), number_of_partitions); + for (auto const & layout : layouts) + EXPECT_EQ(top_level_technical_bins(layout), 128u); + } +} + +namespace +{ + +// A HyperLogLog sketch of the hashes scramble(offset), ..., scramble(offset + count - 1). +seqan::hibf::sketch::hyperloglog make_sketch(uint64_t const offset, size_t const count, uint8_t const sketch_bits) +{ + seqan::hibf::sketch::hyperloglog sketch{sketch_bits}; + for (uint64_t i = 0; i < count; ++i) + sketch.add(scramble(offset + i)); + return sketch; +} + +} // namespace + +TEST(phibf_regression_test, seeding_spills_past_last_partition) +{ + // lsh and lsh_sim initialise the partitions with the first clusters, which post_process_partition_clusters orders by + // their number of user bins, not by cardinality. A cluster that exceeds the target cardinality per partition spills + // over into the next partition. Former bug: `++p` was only guarded by an assertion, so a heavy cluster in the last + // seeded position moved p past the last partition (assertion in Debug, out-of-bounds write in Release). + // Now, user bins that would spill past the last partition are assigned by find_best_partition after seeding. + // Note: phibf_data cannot trigger this (no lsh/lsh_sim input with up to 128 user bins and 16 partitions does). + + // User bins 0-3 (cardinality 1) and 4-6 (cardinality 60) have identical MinHash tables within each group and form + // the clusters {0, 1, 2, 3} and {4, 5, 6}. The first has more user bins and seeds partition 0, the second seeds + // partition 1. The target cardinality per partition is ceil(184 / 2) = 92, which {4, 5, 6} exceeds after two of + // its user bins, so the third one cannot spill over into the next partition. + size_t const number_of_user_bins{7}; + std::vector const cardinalities{1, 1, 1, 1, 60, 60, 60}; + + chopper::configuration config{}; + config.number_of_partitions = 2; + config.hibf_config.number_of_user_bins = number_of_user_bins; + config.hibf_config.tmax = 64; + + auto minhashes_of_group = [](uint64_t const group) + { + using seqan::hibf::sketch::minhashes; + minhashes result{}; + result.table.assign(minhashes::num_sketches, std::vector(minhashes::sketch_size)); + uint64_t value{group << 32}; + for (auto & sketch : result.table) + for (auto & hash : sketch) + hash = scramble(value++); + return result; + }; + + std::vector sketches{}; + std::vector minHash_sketches{}; + for (size_t ub = 0; ub < number_of_user_bins; ++ub) + { + sketches.push_back(make_sketch(uint64_t{ub} << 32, cardinalities[ub], config.hibf_config.sketch_bits)); + minHash_sketches.push_back(minhashes_of_group(ub < 4 ? 1 : 2)); + } + + using chopper::partitioning_scheme; + + for (chopper::partitioning_scheme const approach : {partitioning_scheme::lsh, partitioning_scheme::lsh_sim}) + { + config.partitioning_approach = approach; + + std::vector> partitions(config.number_of_partitions); + chopper::layout::partition_user_bins(config, cardinalities, sketches, minHash_sketches, partitions); + EXPECT_TRUE(each_user_bin_assigned_once(partitions, number_of_user_bins)) << "approach " << approach; + } +} + +TEST(phibf_regression_test, find_best_partition_intersection_underflow) +{ + // find_best_partition estimated the intersection of a cluster and a partition as + // est(cluster) - (est(union) - est(partition)) in size_t. For disjoint sets, HyperLogLog noise makes the difference + // larger than est(cluster) about half of the time. The subtraction then wrapped around, and the partition got a huge + // similarity, so user bins went to a partition they share nothing with. + + chopper::configuration config{}; + config.number_of_partitions = 2; + uint8_t const sketch_bits{config.hibf_config.sketch_bits}; + + // Partition 0 holds 200'000 k-mers that no user bin shares. Partition 1 holds 100'000 k-mers, and each user bin is + // a disjoint subset of 1'000 of them. Each user bin must go to partition 1: there, the union is the partition + // itself, so the similarity is est(user bin) / est(partition 1). Partition 0 can score at most + // est(user bin) / est(partition 0), which is smaller. + uint64_t const offset_1{uint64_t{1} << 40}; + seqan::hibf::sketch::hyperloglog const partition_0 = make_sketch(0, 200'000, sketch_bits); + seqan::hibf::sketch::hyperloglog const partition_1 = make_sketch(offset_1, 100'000, sketch_bits); + + size_t const number_of_user_bins{50}; + std::vector sketches{}; + std::vector cardinalities{}; + for (size_t ub = 0; ub < number_of_user_bins; ++ub) + { + sketches.push_back(make_sketch(offset_1 + ub * 1'000, 1'000, sketch_bits)); + cardinalities.push_back(sketches.back().estimate()); + } + + size_t wrong_partition{0}; + for (size_t ub = 0; ub < number_of_user_bins; ++ub) + { + std::vector> partitions(config.number_of_partitions); + std::vector partition_sketches{partition_0, partition_1}; + std::vector partition_cardinality{0, 0}; // Both partitions have room for the user bin. + + chopper::layout::find_best_partition(config, + /*cardinality_per_part*/ 1'000'000, + {ub}, + cardinalities, + sketches, + partitions, + partition_sketches, + partition_cardinality); + wrong_partition += (partitions[1] != std::vector{ub}); + } + EXPECT_EQ(wrong_partition, 0u) << "of " << number_of_user_bins << " user bins"; +} + +TEST(phibf_regression_test, unassigned_number_of_user_bins) +{ + size_t const number_of_user_bins{50}; + phibf_data data{chopper::partitioning_scheme::lsh, number_of_user_bins, 2, "uniform"}; + data.config.hibf_config.number_of_user_bins = 0; // should not be used, instead use data sizes like cardinalities + + EXPECT_TRUE(each_user_bin_assigned_once(data.partition(), number_of_user_bins)); +} diff --git a/test/api/layout/user_bin_io_test.cpp b/test/api/layout/user_bin_io_test.cpp index 8e9fc97c..c8252298 100644 --- a/test/api/layout/user_bin_io_test.cpp +++ b/test/api/layout/user_bin_io_test.cpp @@ -7,15 +7,18 @@ #include #include -TEST(output, user_bins) +#include + +TEST(output, write_user_bins_to) { std::vector> const filenames{{"file1.fa", "fileB.fa"}, {"file2.fa"}, {"path/to/file3.fa"}, {"file4.fastq"}}; + auto mapping = seqan::hibf::iota_vector(filenames.size()); std::stringstream ss{}; - chopper::layout::write_user_bins_to(filenames, ss); + chopper::layout::write_user_bins_to(filenames, mapping, ss); std::string const expected{"@CHOPPER_USER_BINS\n" "@0 file1.fa fileB.fa\n" @@ -27,7 +30,7 @@ TEST(output, user_bins) EXPECT_EQ(ss.str(), expected); } -TEST(input, user_bins) +TEST(input, read_filenames_from) { std::stringstream ss{"@CHOPPER_USER_BINS\n" "@0 file1.fa fileB.fa\n" diff --git a/test/cli/CMakeLists.txt b/test/cli/CMakeLists.txt index c1b53a46..fd7b7809 100644 --- a/test/cli/CMakeLists.txt +++ b/test/cli/CMakeLists.txt @@ -33,3 +33,9 @@ target_use_datasources (util_display_layout_test FILES seq2.fa) target_use_datasources (util_display_layout_test FILES seq3.fa) target_use_datasources (util_display_layout_test FILES small.fa) add_dependencies (util_display_layout_test display_layout) + +add_cli_test (cli_chopper_phibf_test.cpp) +target_use_datasources (cli_chopper_phibf_test FILES seq1.fa) +target_use_datasources (cli_chopper_phibf_test FILES seq2.fa) +target_use_datasources (cli_chopper_phibf_test FILES seq3.fa) +target_use_datasources (cli_chopper_phibf_test FILES small.fa) diff --git a/test/cli/cli_chopper_layout_from_sketch_file.cpp b/test/cli/cli_chopper_layout_from_sketch_file.cpp index d3c45a9d..5ae645ce 100644 --- a/test/cli/cli_chopper_layout_from_sketch_file.cpp +++ b/test/cli/cli_chopper_layout_from_sketch_file.cpp @@ -81,7 +81,7 @@ TEST_F(cli_test, chopper_layout_from_sketch_file) "@CHOPPER_CONFIG\n" "@{\n" "@ \"chopper_config\": {\n" - "@ \"version\": 3,\n" + "@ \"version\": 4,\n" "@ \"data_file\": {\n" "@ \"value0\": \"" + input_filename.string() @@ -102,7 +102,9 @@ TEST_F(cli_test, chopper_layout_from_sketch_file) "@ },\n" "@ \"determine_best_tmax\": false,\n" "@ \"force_all_binnings\": false,\n" - "@ \"fast_layout\": false\n" + "@ \"fast_layout\": false,\n" + "@ \"number_of_partitions\": 1,\n" + "@ \"partitioning_approach\": 0\n" "@ }\n" "@}\n" "@CHOPPER_CONFIG_END\n" diff --git a/test/cli/cli_chopper_phibf_test.cpp b/test/cli/cli_chopper_phibf_test.cpp new file mode 100644 index 00000000..d5f89dee --- /dev/null +++ b/test/cli/cli_chopper_phibf_test.cpp @@ -0,0 +1,330 @@ +// --------------------------------------------------------------------------------------------------- +// 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 "../api/api_test.hpp" +#include "cli_test.hpp" + +namespace +{ + +std::filesystem::path write_input_file(std::filesystem::path const & directory) +{ + std::filesystem::path const input_filename{directory / "data.tsv"}; + std::ofstream fout{input_filename}; + fout << data("seq1.fa").string() << '\n' + << data("seq2.fa").string() << '\n' + << data("seq3.fa").string() << '\n' + << data("small.fa").string() << '\n'; + return input_filename; +} + +} // namespace + +TEST_F(cli_test, chopper_layout_phibf) +{ + seqan3::test::tmp_directory tmp_dir{}; + std::filesystem::path const input_filename{write_input_file(tmp_dir.path())}; + std::filesystem::path const layout_filename{tmp_dir.path() / "phibf.layout"}; + + cli_test_result result = execute_app("chopper", + "--input", + input_filename.c_str(), + "--number-of-partitions", + "2", + "--partitioning-approach", + "sorted", + "--output", + layout_filename.c_str()); + + EXPECT_EQ(result.exit_code, 0); + EXPECT_EQ(result.out, std::string{}); + EXPECT_EQ(result.err, std::string{}); + + std::ifstream layout_stream{layout_filename}; + + ASSERT_TRUE(layout_stream.good()); + // Layout 1 + { + auto [filenames, config, layout] = chopper::layout::read_layout_file(layout_stream); + EXPECT_EQ(config.number_of_partitions, 2u); + EXPECT_EQ(config.partitioning_approach, chopper::partitioning_scheme::sorted); + EXPECT_EQ(filenames, + (std::vector>{std::vector{data("seq2.fa")}, + std::vector{data("small.fa")}})); + + ASSERT_EQ(layout.user_bins.size(), 2); + EXPECT_EQ(layout.user_bins[0].idx, 0); + EXPECT_EQ(layout.user_bins[1].idx, 1); + } + + ASSERT_TRUE(layout_stream.good()); + // Layout 2 + { + auto [filenames, config, layout] = chopper::layout::read_layout_file(layout_stream); + EXPECT_EQ(config.number_of_partitions, 2u); + EXPECT_EQ(config.partitioning_approach, chopper::partitioning_scheme::sorted); + + ASSERT_EQ(layout.user_bins.size(), 2); + EXPECT_EQ(layout.user_bins[0].idx, 0); + EXPECT_EQ(layout.user_bins[1].idx, 1); + EXPECT_EQ(filenames, + (std::vector>{std::vector{data("seq1.fa")}, + std::vector{data("seq3.fa")}})); + } + + EXPECT_FALSE(layout_stream.good()); +} + +TEST_F(cli_test, chopper_layout_phibf_all_approaches) +{ + + std::vector const phibf_approaches{chopper::partitioning_scheme::blocked, + chopper::partitioning_scheme::sorted, + chopper::partitioning_scheme::folded, + chopper::partitioning_scheme::weighted_fold, + chopper::partitioning_scheme::similarity, + chopper::partitioning_scheme::lsh, + chopper::partitioning_scheme::lsh_sim}; + + for (chopper::partitioning_scheme const approach : phibf_approaches) + { + std::stringstream ss{}; + ss << approach; + seqan3::test::tmp_directory tmp_dir{}; + std::filesystem::path const input_filename{write_input_file(tmp_dir.path())}; + std::filesystem::path const layout_filename{tmp_dir.path() / "phibf.layout"}; + + cli_test_result result = execute_app("chopper", + "--input", + input_filename.c_str(), + "--number-of-partitions", + "2", + "--partitioning-approach", + ss.str(), + "--output", + layout_filename.c_str()); + + EXPECT_EQ(result.exit_code, 0); + EXPECT_EQ(result.out, std::string{}); + EXPECT_EQ(result.err, std::string{}); + } +} + +TEST_F(cli_test, chopper_layout_phibf_with_fast_layout) +{ + seqan3::test::tmp_directory tmp_dir{}; + std::filesystem::path const input_filename{write_input_file(tmp_dir.path())}; + std::filesystem::path const layout_filename{tmp_dir.path() / "phibf.layout"}; + + cli_test_result result = execute_app("chopper", + "--input", + input_filename.c_str(), + "--number-of-partitions", + "2", + "--fast-layout", + "--output", + layout_filename.c_str()); + + 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 --number-of-partitions.\n"}); +} + +TEST_F(cli_test, chopper_layout_phibf_with_determine_best_tmax) +{ + seqan3::test::tmp_directory tmp_dir{}; + std::filesystem::path const input_filename{write_input_file(tmp_dir.path())}; + std::filesystem::path const layout_filename{tmp_dir.path() / "phibf.layout"}; + + cli_test_result result = execute_app("chopper", + "--input", + input_filename.c_str(), + "--number-of-partitions", + "2", + "--determine-best-tmax", + "--output", + layout_filename.c_str()); + + EXPECT_NE(result.exit_code, 0); + EXPECT_EQ(result.out, std::string{}); + EXPECT_EQ(result.err, + std::string{"[ERROR] You cannot combine --determine-best-tmax with --number-of-partitions.\n"}); +} + +TEST_F(cli_test, chopper_layout_phibf_from_sketch_file_without_minhashes) +{ + seqan3::test::tmp_directory tmp_dir{}; + std::filesystem::path const input_filename{write_input_file(tmp_dir.path())}; + std::filesystem::path const sketch_filename{tmp_dir.path() / "default.sketch"}; + + // The default layout does not compute MinHash sketches. + { + std::filesystem::path const layout_filename{tmp_dir.path() / "default.layout"}; + cli_test_result result = execute_app("chopper", + "--input", + input_filename.c_str(), + "--output-sketches-to", + sketch_filename.c_str(), + "--output", + layout_filename.c_str()); + ASSERT_EQ(result.exit_code, 0); + } + + // lsh_sim needs MinHash sketches. + { + std::filesystem::path const layout_filename{tmp_dir.path() / "lsh_sim.layout"}; + cli_test_result result = execute_app("chopper", + "--input", + sketch_filename.c_str(), + "--number-of-partitions", + "2", + "--partitioning-approach", + "lsh_sim", + "--output", + layout_filename.c_str()); + + EXPECT_NE(result.exit_code, 0); + EXPECT_EQ(result.out, std::string{}); + EXPECT_EQ(result.err, + std::string{"[ERROR] The sketch file does not contain MinHash sketches, which the chosen " + "--partitioning-approach needs. Create the sketch file with --number-of-partitions.\n"}); + } + + // sorted does not need MinHash sketches. + { + std::filesystem::path const layout_filename{tmp_dir.path() / "sorted.layout"}; + cli_test_result result = execute_app("chopper", + "--input", + sketch_filename.c_str(), + "--number-of-partitions", + "2", + "--partitioning-approach", + "sorted", + "--output", + layout_filename.c_str()); + + EXPECT_EQ(result.exit_code, 0); + EXPECT_EQ(result.out, std::string{}); + EXPECT_EQ(result.err, std::string{}); + + std::ifstream layout_stream{layout_filename}; + ASSERT_TRUE(layout_stream.good()); + // Layout 1 + { + auto [filenames, config, layout] = chopper::layout::read_layout_file(layout_stream); + EXPECT_EQ(config.partitioning_approach, chopper::partitioning_scheme::sorted); + } + ASSERT_TRUE(layout_stream.good()); + // Layout 2 + { + auto [filenames, config, layout] = chopper::layout::read_layout_file(layout_stream); + EXPECT_EQ(config.partitioning_approach, chopper::partitioning_scheme::sorted); + } + EXPECT_FALSE(layout_stream.good()); + } +} + +TEST_F(cli_test, chopper_layout_phibf_unknown_partitioning_approach) +{ + seqan3::test::tmp_directory tmp_dir{}; + std::filesystem::path const input_filename{write_input_file(tmp_dir.path())}; + std::filesystem::path const layout_filename{tmp_dir.path() / "phibf.layout"}; + + cli_test_result result = execute_app("chopper", + "--input", + input_filename.c_str(), + "--number-of-partitions", + "2", + "--partitioning-approach", + "foo", + "--output", + layout_filename.c_str()); + + EXPECT_NE(result.exit_code, 0); + EXPECT_EQ(result.out, std::string{}); + EXPECT_EQ(result.err, + std::string{"[ERROR] Validation failed for option --partitioning-approach: Value unknown is not one of " + "[blocked, sorted, folded, weighted_fold, similarity, lsh, lsh_sim].\n"}); +} + +TEST_F(cli_test, chopper_layout_phibf_more_partitions_than_user_bins) +{ + seqan3::test::tmp_directory tmp_dir{}; + std::filesystem::path const input_filename{write_input_file(tmp_dir.path())}; + std::filesystem::path const layout_filename{tmp_dir.path() / "phibf.layout"}; + + cli_test_result result = execute_app("chopper", + "--input", + input_filename.c_str(), + "--number-of-partitions", + "5", + "--output", + layout_filename.c_str()); + + EXPECT_NE(result.exit_code, 0); + EXPECT_EQ(result.out, std::string{}); + EXPECT_EQ(result.err, + std::string{"[ERROR] The number of partitions (5) must not exceed the number of user bins (4).\n"}); +} + +TEST_F(cli_test, chopper_layout_phibf_tmax) +{ + seqan3::test::tmp_directory tmp_dir{}; + std::filesystem::path const input_filename{write_input_file(tmp_dir.path())}; + size_t const number_of_partitions{2}; + + for (bool const set_tmax : {false, true}) + { + std::filesystem::path const layout_filename{tmp_dir.path() / (set_tmax ? "tmax.layout" : "default.layout")}; + std::string const tmax_option{set_tmax ? "--tmax 128" : ""}; + + cli_test_result result = execute_app("chopper", + "--input", + input_filename.c_str(), + "--number-of-partitions", + std::to_string(number_of_partitions), + "--partitioning-approach", + "sorted", + tmax_option, + "--output", + layout_filename.c_str()); + ASSERT_EQ(result.exit_code, 0) << result.err; + + // Without --tmax, each partition (2 user bins) gets next_multiple_of_64(ceil(sqrt(2))) = 64. + size_t const expected_tmax = set_tmax ? 128u : 64u; + + std::ifstream layout_stream{layout_filename}; + size_t number_of_layouts{0}; + while (layout_stream.good()) + { + auto [filenames, config, layout] = chopper::layout::read_layout_file(layout_stream); + ++number_of_layouts; + + size_t tbs{}; // count tmax + for (auto const & user_bin : layout.user_bins) + tbs = std::max(tbs, + user_bin.previous_TB_indices.empty() + ? user_bin.storage_TB_id + user_bin.number_of_technical_bins + : user_bin.previous_TB_indices[0] + 1); + + EXPECT_EQ(tbs, expected_tmax); + } + EXPECT_EQ(number_of_partitions, number_of_layouts); + } +} diff --git a/test/cli/cli_chopper_pipeline_test.cpp b/test/cli/cli_chopper_pipeline_test.cpp index 53fad559..63bf63c2 100644 --- a/test/cli/cli_chopper_pipeline_test.cpp +++ b/test/cli/cli_chopper_pipeline_test.cpp @@ -64,7 +64,7 @@ TEST_F(cli_test, chopper_layout) "@CHOPPER_CONFIG\n" "@{\n" "@ \"chopper_config\": {\n" - "@ \"version\": 3,\n" + "@ \"version\": 4,\n" "@ \"data_file\": {\n" "@ \"value0\": \"" + taxa_filename.string() @@ -85,7 +85,9 @@ TEST_F(cli_test, chopper_layout) "@ },\n" "@ \"determine_best_tmax\": false,\n" "@ \"force_all_binnings\": false,\n" - "@ \"fast_layout\": false\n" + "@ \"fast_layout\": false,\n" + "@ \"number_of_partitions\": 1,\n" + "@ \"partitioning_approach\": 0\n" "@ }\n" "@}\n" "@CHOPPER_CONFIG_END\n" @@ -169,7 +171,7 @@ TEST_F(cli_test, chopper_layout2) "@CHOPPER_CONFIG\n" "@{\n" "@ \"chopper_config\": {\n" - "@ \"version\": 3,\n" + "@ \"version\": 4,\n" "@ \"data_file\": {\n" "@ \"value0\": \"" + taxa_filename.string() @@ -190,7 +192,9 @@ TEST_F(cli_test, chopper_layout2) "@ },\n" "@ \"determine_best_tmax\": false,\n" "@ \"force_all_binnings\": false,\n" - "@ \"fast_layout\": false\n" + "@ \"fast_layout\": false,\n" + "@ \"number_of_partitions\": 1,\n" + "@ \"partitioning_approach\": 0\n" "@ }\n" "@}\n" "@CHOPPER_CONFIG_END\n" @@ -275,7 +279,7 @@ TEST_F(cli_test, chopper_layout2_fast_layout) "@CHOPPER_CONFIG\n" "@{\n" "@ \"chopper_config\": {\n" - "@ \"version\": 3,\n" + "@ \"version\": 4,\n" "@ \"data_file\": {\n" "@ \"value0\": \"" + taxa_filename.string() @@ -296,7 +300,9 @@ TEST_F(cli_test, chopper_layout2_fast_layout) "@ },\n" "@ \"determine_best_tmax\": false,\n" "@ \"force_all_binnings\": false,\n" - "@ \"fast_layout\": true\n" + "@ \"fast_layout\": true,\n" + "@ \"number_of_partitions\": 1,\n" + "@ \"partitioning_approach\": 0\n" "@ }\n" "@}\n" "@CHOPPER_CONFIG_END\n" diff --git a/test/cli/util_display_layout_test.cpp b/test/cli/util_display_layout_test.cpp index 2a1fd934..5ad96ffd 100644 --- a/test/cli/util_display_layout_test.cpp +++ b/test/cli/util_display_layout_test.cpp @@ -57,6 +57,8 @@ std::string get_layout_with_correct_filenames(std::string_view const seq1_filena "@ \"window_size\": 15,\n" "@ \"disable_sketch_output\": true,\n" "@ \"precomputed_files\": false,\n" + "@ \"number_of_partitions\": 1,\n" + "@ \"partitioning_approach\": 0,\n" "@ \"output_filename\": {\n" "@ \"value0\": \"" + output_filename.data()