WIP: [FEATURE] Add the LSH fast layout algorithm.

.gitignore

@@ -32,6 +32,7 @@
 
 # build dir
 /build/
+/build_debug/
 
 # editor configuration
 /.vscode/

include/chopper/configuration.hpp

@@ -22,6 +22,9 @@ namespace chopper
 
 struct configuration
 {
+    //!\brief Whether to use the fast layout algorithm instead of the default one.
+    bool fast_layout{false};
+
     /*!\name General Configuration
      * \{
      */
@@ -77,6 +80,23 @@ struct configuration
     mutable seqan::hibf::concurrent_timer union_estimation_timer{};
     mutable seqan::hibf::concurrent_timer rearrangement_timer{};
     mutable seqan::hibf::concurrent_timer dp_algorithm_timer{};
+    /*!\brief Fast layout: time spent in LSH clustering (`lsh_in_seconds` in the timing output).
+     *
+     * Summed over all partitionings, including concurrent ones, so it can exceed the wall-clock time.
+     */
+    mutable seqan::hibf::concurrent_timer lsh_algorithm_timer{};
+    /*!\brief Fast layout: time spent assigning clusters to partitions by similarity (`search_best_p_in_seconds` in
+     *        the timing output).
+     *
+     * Summed over all partitionings, including concurrent ones, so it can exceed the wall-clock time.
+     */
+    mutable seqan::hibf::concurrent_timer search_partition_algorithm_timer{};
+    /*!\brief Fast layout: time for partitioning the top level (`initial_partition_timer_in_seconds` in the timing
+     *        output).
+     */
+    mutable seqan::hibf::concurrent_timer initial_partition_timer{};
+    //!\brief Fast layout: time for laying out all lower levels (`small_layouts_timer_in_seconds` in the timing output).
+    mutable seqan::hibf::concurrent_timer small_layouts_timer{};
 
     void read_from(std::istream & stream);
 

include/chopper/layout/determine_split_bins.hpp

@@ -0,0 +1,53 @@
+// --------------------------------------------------------------------------------------------------
+// Copyright (c) 2006-2023, Knut Reinert & Freie Universität Berlin
+// Copyright (c) 2016-2023, Knut Reinert & MPI für molekulare Genetik
+// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License
+// shipped with this file and also available at: https://github.com/seqan/chopper/blob/main/LICENSE.md
+// --------------------------------------------------------------------------------------------------
+
+/*!\file
+ * \brief Provides chopper::layout::determine_split_bins.
+ * \author Svenja Mehringer <svenja.mehringer AT fu-berlin.de>
+ */
+
+#pragma once
+
+#include <cstddef>
+#include <utility>
+#include <vector>
+
+#include <chopper/configuration.hpp>
+
+namespace chopper::layout
+{
+
+/*!\brief Distributes the largest user bins as split bins over technical bins at the back of `partitions`.
+ * \param[in]     config             The configuration (uses `maximum_fpr` and `number_of_hash_functions` of
+ *                                   `config.hibf_config`).
+ * \param[in]     positions          User bin indices, sorted by descending cardinality. The first `num_user_bins`
+ *                                   are split.
+ * \param[in]     cardinalities      The cardinality of each user bin, indexed by the values in `positions`.
+ * \param[in]     num_technical_bins The maximum number of technical bins for the split user bins. Must be at least
+ *                                   `num_user_bins`.
+ * \param[in]     num_user_bins      The number of user bins to split.
+ * \param[in,out] partitions         The technical bins. Split user bin indices are appended to the last technical
+ *                                   bins, one index per technical bin. `positions[num_user_bins - 1]` gets the last
+ *                                   technical bins, `positions[0]` the first of them. The technical bins of each user
+ *                                   bin are consecutive.
+ * \returns A pair of
+ *          1. the number of technical bins used, which may be less than `num_technical_bins`, and
+ *          2. the largest FPR-corrected cardinality per technical bin.
+ *          `{0, 0}` if `num_user_bins == 0`.
+ *
+ * A dynamic programming algorithm assigns each user bin a number of technical bins, such that the largest
+ * FPR-corrected cardinality per technical bin, `ceil(cardinality * fpr_correction[k] / k)` for a user bin in `k`
+ * technical bins, is minimal. If several numbers of technical bins give the same minimum, the smallest is used.
+ */
+std::pair<size_t, size_t> determine_split_bins(chopper::configuration const & config,
+                                               std::vector<size_t> const & positions,
+                                               std::vector<size_t> const & cardinalities,
+                                               size_t const num_technical_bins,
+                                               size_t const num_user_bins,
+                                               std::vector<std::vector<size_t>> & partitions);
+
+} // namespace chopper::layout

include/chopper/layout/execute.hpp

@@ -13,12 +13,32 @@
 #include <chopper/configuration.hpp>
 
 #include <hibf/sketch/hyperloglog.hpp>
+#include <hibf/sketch/minhashes.hpp>
 
 namespace chopper::layout
 {
 
+/*!\brief Computes the layout for the given user bins and writes it to `config.output_filename`.
+ * \param[in,out] config           The configuration. `config.hibf_config` is validated and completed
+ *                                 (`validate_and_set_defaults`), and the timers are updated.
+ * \param[in]     filenames        The file names of each user bin. They are written to the layout file.
+ * \param[in]     sketches         The HyperLogLog sketch of each user bin.
+ * \param[in]     minHash_sketches The MinHash sketches of each user bin. Only used, and then required for every user
+ *                                 bin, if `config.fast_layout` is set. May be empty otherwise.
+ * \returns 0.
+ * \throws std::invalid_argument If both `config.determine_best_tmax` and `config.fast_layout` are set.
+ *
+ * The layout is computed with
+ * - determine_best_number_of_technical_bins if `config.determine_best_tmax` is set,
+ * - fast_layout if `config.fast_layout` is set,
+ * - the DP layout of the HIBF library (seqan::hibf::layout::compute_layout) otherwise.
+ *
+ * Unless `config.determine_best_tmax` is set, `config.output_verbose_statistics` prints statistics of the layout to
+ * `std::cout`.
+ */
 int execute(chopper::configuration & config,
             std::vector<std::vector<std::string>> const & filenames,
-            std::vector<seqan::hibf::sketch::hyperloglog> const & sketches);
+            std::vector<seqan::hibf::sketch::hyperloglog> const & sketches,
+            std::vector<seqan::hibf::sketch::minhashes> const & minHash_sketches);
 
 } // namespace chopper::layout

include/chopper/layout/fast_layout.hpp

@@ -0,0 +1,50 @@
+// ---------------------------------------------------------------------------------------------------
+// Copyright (c) 2006-2023, Knut Reinert & Freie Universität Berlin
+// Copyright (c) 2016-2023, Knut Reinert & MPI für molekulare Genetik
+// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License
+// shipped with this file and also available at: https://github.com/seqan/chopper/blob/main/LICENSE.md
+// ---------------------------------------------------------------------------------------------------
+
+#pragma once
+
+#include <vector>
+
+#include <chopper/configuration.hpp>
+
+#include <hibf/layout/layout.hpp>
+#include <hibf/sketch/hyperloglog.hpp>
+#include <hibf/sketch/minhashes.hpp>
+
+namespace chopper::layout
+{
+
+/*!\brief Computes an HIBF layout with the fast (LSH/similarity-based) layout algorithm.
+ * \param[in]  config           The configuration. Uses `hibf_config` (`tmax`, `number_of_user_bins`, FPR settings)
+ *                              and the fast-layout timers.
+ * \param[in]  positions        The global indices of all user bins. Must be a permutation of
+ *                              `[0, config.hibf_config.number_of_user_bins)`.
+ * \param[in]  cardinalities    The cardinality of each user bin, indexed by global user bin index.
+ * \param[in]  sketches         The HyperLogLog sketch of each user bin, indexed by global user bin index.
+ * \param[in]  minHash_sketches The MinHash tables of each user bin, indexed by global user bin index.
+ * \param[out] hibf_layout      The resulting layout. Expected to be empty on entry. `user_bins` is resized to
+ *                              `number_of_user_bins`, so `user_bins[i].idx == i`.
+ *
+ * 1. **Top level:** partition_user_bins distributes all user bins onto `tmax` technical bins. The user bins are
+ *    initialised in `hibf_layout`: merged bins get `previous_TB_indices = {t}`; split and single bins get
+ *    `storage_TB_id` and their number of consecutive technical bins. `top_level_max_bin_id` is set to the technical
+ *    bin with the largest FPR-corrected size.
+ * 2. **Lower levels:** Merged bins are processed in parallel (OpenMP `taskloop`). Depending on
+ *    do_I_need_a_fast_layout, each is laid out recursively with the fast layout or with the regular DP layout,
+ *    and the result is grafted into `hibf_layout`.
+ * 3. `max_bins` is sorted by level (path length), then lexicographically by path.
+ *
+ * \throws std::logic_error In debug builds, if the layout's user bins are not a permutation of `positions`.
+ */
+void fast_layout(chopper::configuration const & config,
+                 std::vector<size_t> const & positions,
+                 std::vector<size_t> const & cardinalities,
+                 std::vector<seqan::hibf::sketch::hyperloglog> const & sketches,
+                 std::vector<seqan::hibf::sketch::minhashes> const & minHash_sketches,
+                 seqan::hibf::layout::layout & hibf_layout);
+
+} // namespace chopper::layout

include/chopper/layout/fast_layout_cluster.hpp

@@ -0,0 +1,187 @@
+// ---------------------------------------------------------------------------------------------------
+// Copyright (c) 2006-2023, Knut Reinert & Freie Universität Berlin
+// Copyright (c) 2016-2023, Knut Reinert & MPI für molekulare Genetik
+// This file may be used, modified and/or redistributed under the terms of the 3-clause BSD-License
+// shipped with this file and also available at: https://github.com/seqan/chopper/blob/main/LICENSE.md
+// ---------------------------------------------------------------------------------------------------
+
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <functional>
+#include <optional>
+#include <vector>
+
+namespace chopper::layout
+{
+
+/*!\brief A cluster of user bins, used by the LSH clustering of the fast layout.
+ *
+ * In a vector of clusters in which the cluster at position `i` was created with id `i`, the cluster at position `i`
+ * keeps `id() == i` and is either
+ * - **valid**: it has not been moved and contains at least one user bin, or
+ * - **moved**: its user bins were moved into another cluster (move_to), it is empty, and moved_to_cluster_id() is the
+ *   id of that cluster. That cluster may itself have been moved, so moves can form a chain.
+ *   LSH_find_representative_cluster follows the chain to the valid cluster.
+ *
+ * Note that `is_valid(i)` is true in both states: it checks that `id() == i` and that the cluster is in one of them.
+ */
+struct Cluster
+{
+private:
+    //!\brief The id of the cluster.
+    size_t representative_id{};
+
+    //!\brief The user bins contained in this cluster.
+    std::vector<size_t> user_bins{};
+
+    //!\brief The id of the cluster that this cluster's user bins were moved to, if they were moved.
+    std::optional<size_t> moved_id{std::nullopt};
+
+public:
+    /*!\name Constructors, destructor and assignment
+     * \{
+     */
+    Cluster() = default;                            //!< Defaulted.
+    Cluster(Cluster const &) = default;             //!< Defaulted.
+    Cluster(Cluster &&) = default;                  //!< Defaulted.
+    Cluster & operator=(Cluster const &) = default; //!< Defaulted.
+    Cluster & operator=(Cluster &&) = default;      //!< Defaulted.
+    ~Cluster() = default;                           //!< Defaulted.
+
+    /*!\brief Creates a cluster with id `id` that contains the single user bin `user_bins_id`.
+     * \param[in] id           The id of the cluster.
+     * \param[in] user_bins_id The user bin the cluster contains.
+     */
+    Cluster(size_t const id, size_t const user_bins_id) : representative_id{id}, user_bins({user_bins_id})
+    {}
+
+    /*!\brief Creates a cluster with id `id` that contains the single user bin `id`.
+     * \param[in] id The id of the cluster and the user bin it contains.
+     */
+    explicit Cluster(size_t const id) : Cluster{id, id}
+    {}
+    //!\}
+
+    //!\brief Returns the id of the cluster.
+    size_t id() const
+    {
+        return representative_id;
+    }
+
+    //!\brief Returns the user bins contained in the cluster.
+    std::vector<size_t> const & contained_user_bins() const
+    {
+        return user_bins;
+    }
+
+    //!\brief Returns whether the user bins of this cluster were moved to another cluster (see move_to).
+    bool has_been_moved() const
+    {
+        return moved_id.has_value();
+    }
+
+    //!\brief Returns whether the cluster contains no user bins.
+    bool empty() const
+    {
+        return user_bins.empty();
+    }
+
+    //!\brief Returns the number of user bins in the cluster.
+    size_t size() const
+    {
+        return user_bins.size();
+    }
+
+    //!\brief Removes the last user bin from the cluster and returns it. The cluster must not be empty.
+    size_t pop_back()
+    {
+        size_t last = user_bins.back();
+        user_bins.pop_back();
+        return last;
+    }
+
+    /*!\brief Adds a user bin to the cluster.
+     * \param[in] user_bin The user bin to add.
+     */
+    void add_user_bin(size_t const user_bin)
+    {
+        user_bins.push_back(user_bin);
+    }
+
+    /*!\brief Checks that the cluster at position `id` is either valid or moved, see Cluster.
+     * \param[in] id The position of the cluster, i.e., the id it must have.
+     * \returns `true` if `id() == id`, and the cluster is either not moved and non-empty, or moved and empty.
+     */
+    bool is_valid(size_t const id) const
+    {
+        bool const ids_equal = representative_id == id;
+        bool const properly_moved = has_been_moved() && empty();
+        bool const not_moved = !has_been_moved() && !empty();
+
+        return ids_equal && (properly_moved || not_moved);
+    }
+
+    //!\brief Returns the id of the cluster that this cluster's user bins were moved to. The cluster must be moved.
+    size_t moved_to_cluster_id() const
+    {
+        assert(moved_id.has_value());
+        assert(is_valid(representative_id));
+        return moved_id.value();
+    }
+
+    /*!\brief Moves all user bins of this cluster to `target_cluster` and marks this cluster as moved.
+     * \param[in,out] target_cluster The cluster to move the user bins to. Must be a different cluster.
+     *
+     * Afterwards, this cluster is empty, its memory is released, and moved_to_cluster_id() is `target_cluster.id()`.
+     */
+    void move_to(Cluster & target_cluster)
+    {
+        auto & target = target_cluster.user_bins;
+        auto & source = this->user_bins;
+        target.insert(target.end(), source.cbegin(), source.cend());
+        source = std::vector<size_t>{}; // .clear() AND release memory
+
+        moved_id = target_cluster.id();
+    }
+
+    /*!\brief Sorts the user bins of the cluster by descending cardinality.
+     * \param[in] cardinalities The cardinality of each user bin, indexed by user bin.
+     */
+    void sort_by_cardinality(std::vector<size_t> const & cardinalities)
+    {
+        std::ranges::sort(user_bins,
+                          std::ranges::greater{},
+                          [&cardinalities](size_t const i)
+                          {
+                              return cardinalities[i];
+                          });
+    }
+};
+
+/*!\brief Returns the position of the representative cluster of `clusters[current_id]`.
+ * \param[in] clusters   The clusters. The cluster at position `i` must have id `i`, see Cluster.
+ * \param[in] current_id The position of the cluster to start from.
+ * \returns The position of the representative cluster, i.e., the valid cluster that holds the user bins of
+ *          `clusters[current_id]` now.
+ *
+ * Follows the chain of moves, starting at `clusters[current_id]`. See Cluster for valid and moved clusters.
+ */
+inline size_t LSH_find_representative_cluster(std::vector<Cluster> const & clusters, size_t current_id)
+{
+    std::reference_wrapper<Cluster const> representative = clusters[current_id];
+
+    assert(representative.get().is_valid(current_id));
+
+    while (representative.get().has_been_moved())
+    {
+        current_id = representative.get().moved_to_cluster_id();
+        representative = clusters[current_id]; // replace by next cluster
+        assert(representative.get().is_valid(current_id));
+    }
+
+    return current_id;
+}
+
+} // namespace chopper::layout