Feature: PHIBF

include/chopper/configuration.hpp

@@ -10,6 +10,9 @@
 #include <cinttypes>
 #include <filesystem>
 #include <iosfwd>
+#include <istream>
+#include <ostream>
+#include <string>
 
 #include <cereal/cereal.hpp>
 
@@ -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

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.

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,16 +175,19 @@ 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
  *          `clusters[current_id]` now.
  *
  * Follows the chain of moves, starting at `clusters[current_id]`. See Cluster for valid and moved clusters.
  */
-inline size_t LSH_find_representative_cluster(std::vector<Cluster> const & clusters, size_t current_id)
+template <typename cluster_t>
+size_t LSH_find_representative_cluster(std::vector<cluster_t> const & clusters, size_t current_id)
 {
-    std::reference_wrapper<Cluster const> representative = clusters[current_id];
+    std::reference_wrapper<cluster_t const> representative = clusters[current_id];
 
     assert(representative.get().is_valid(current_id));
 

include/chopper/layout/input.hpp

@@ -21,12 +21,12 @@ namespace chopper::layout
 
 std::vector<std::vector<std::string>> 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<std::vector<std::vector<std::string>>, configuration, seqan::hibf::layout::layout>
 read_layout_file(std::istream & stream);

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 <svenja.mehringer AT fu-berlin.de>
+ */
+
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cstddef>
+#include <cstdint>
+#include <numeric>
+#include <ranges>
+#include <vector>
+
+#include <hibf/contrib/robin_hood.hpp>
+#include <hibf/sketch/minhashes.hpp>
+
+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<uint64_t> 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 <typename cluster_t>
+auto LSH_fill_hashtable(std::vector<cluster_t> const & clusters,
+                        std::vector<seqan::hibf::sketch::minhashes> const & minHash_sketches,
+                        size_t const current_sketch_index,
+                        size_t const current_number_of_sketch_hashes)
+{
+    robin_hood::unordered_flat_map<uint64_t, std::vector<size_t>> 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<decltype(current.contained_user_bins())>;
+
+        if constexpr (std::ranges::range<std::ranges::range_value_t<user_bins_t>>) // 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