[Feature] Too few kmers in sketch

include/hibf/sketch/minhashes.hpp

@@ -35,11 +35,30 @@ namespace seqan::hibf::sketch
 struct minhashes
 {
     static constexpr uint64_t register_id_mask{15}; /// ...00001111
-    static constexpr size_t num_sketches{16};
-    static constexpr size_t sketch_size{40};
+    static constexpr size_t num_sketches{16};       // lsh parameter b
+    static constexpr size_t sketch_size{40};        // lsh parameter r
 
-    //!\brief A table of sketches. For LSH we need multiple sketches, stored in a table.
-    std::vector<std::vector<uint64_t>> table{}; // Each element (vector<uint64_t>) is a minhash.
+    struct minhash_sketch
+    {
+        std::vector<uint64_t> hashes{};
+
+        bool is_valid() const
+        {
+            return hashes.size() == sketch_size;
+        }
+
+    private:
+        friend class cereal::access;
+
+        template <typename archive_t>
+        void serialize(archive_t & archive)
+        {
+            archive(CEREAL_NVP(hashes));
+        }
+    };
+
+    //!\brief A table of minhash sketches. For LSH we need multiple sketches, stored in a table.
+    std::vector<minhash_sketch> table{};
 
     /*!\name Constructors, destructor and assignment
      * \{

src/sketch/compute_sketches.cpp

@@ -3,15 +3,12 @@
 // SPDX-License-Identifier: BSD-3-Clause
 
 #include <algorithm>  // for __fn, sort_heap, all_of, min
-#include <atomic>     // for atomic_flag
 #include <cassert>    // for assert
 #include <cstddef>    // for size_t
 #include <cstdint>    // for uint64_t, uint8_t
 #include <functional> // for equal_to, function
 #include <limits>     // for numeric_limits
 #include <span>       // for span
-#include <stdexcept>  // for runtime_error
-#include <string>     // for operator+, allocator, char_traits, to_string, basic_string
 #include <utility>    // for move
 #include <vector>     // for vector
 
@@ -42,41 +39,6 @@ void compute_sketches(config const & config, std::vector<sketch::hyperloglog> &
     }
 }
 
-/*!\brief Encapsulates handling of too few kmers to compute minHash sketches.
- * \details
- * OMP does not allow just throwing in a thread; it requires that the throw is handled in the same thread.
- * In order to properly throw the exeception in the main thread, we need to set a flag.
- * This flag then tells the other threads to stop working.
- * Afterwards, the main thread can throw the exception.
- */
-struct too_few_kmers_handler
-{
-    std::atomic_flag flag{};
-    size_t available{};
-
-    inline bool test() noexcept
-    {
-        return flag.test();
-    }
-
-    inline void set(size_t const new_available) noexcept
-    {
-        // Sets the flag to true and returns previous value.
-        // If the flag was already set, another thread encountered this block at the same time.
-        // This basically acts as a mutex for setting available.
-        if (!flag.test_and_set())
-            this->available = new_available;
-    }
-
-    inline void check_and_throw()
-    {
-        if (test())
-            throw std::runtime_error{"Not enough kmers (" + std::to_string(available) + ") to get "
-                                     + std::to_string(minhashes::num_sketches * minhashes::sketch_size)
-                                     + " hashes for all minHash sketches."};
-    }
-};
-
 // minhash_sketches data structure:
 // Vector L1 : number of user bins
 // Vector L2 : number_of_max_minhash_sketches (LSH ADD+OR parameter b)
@@ -92,17 +54,11 @@ void compute_sketches(config const & config,
     // compute hll_sketches
     robin_hood::unordered_flat_set<uint64_t> kmers{};
 
-    too_few_kmers_handler too_few_kmers{};
-
     std::vector<uint64_t> heap{};
 
 #pragma omp parallel for schedule(dynamic) num_threads(config.threads) private(kmers) private(heap)
     for (size_t i = 0; i < config.number_of_user_bins; ++i)
     {
-        // Skip work if we already know that we have too few kmers.
-        if (too_few_kmers.test())
-            continue;
-
         seqan::hibf::sketch::hyperloglog hll_sketch{config.sketch_bits};
 
         kmers.clear();
@@ -124,13 +80,11 @@ void compute_sketches(config const & config,
 
         // In case the former heap_size wasn't sufficient to fill the minhash table, we need a larger heap.
         // Once sorted, it begins with the values of the previous heap, which are already in minhash_sketch.
-        while (!minhash_sketch.is_valid() && !too_few_kmers.test())
+        while (!minhash_sketch.is_valid())
         {
+            // All kmers are in the heap: too few to fill the minhash table, which stays invalid.
             if (heap_size == kmers.size())
-            {
-                too_few_kmers.set(kmers.size());
                 break;
-            }
 
             size_t const previous_heap_size = heap_size;
             heap_size = std::min(2u * heap_size, kmers.size());
@@ -153,8 +107,6 @@ void compute_sketches(config const & config,
         hll_sketches[i] = std::move(hll_sketch);
         minhash_sketches[i] = std::move(minhash_sketch);
     }
-
-    too_few_kmers.check_and_throw();
 }
 
 } // namespace seqan::hibf::sketch

src/sketch/minhashes.cpp

@@ -21,24 +21,19 @@ minhashes::minhashes(std::vector<uint64_t> const & smallest_values)
     table.resize(num_sketches);
 
     for (auto & elem : table)
-        elem.reserve(sketch_size);
+        elem.hashes.reserve(sketch_size);
 
     for (uint64_t const hash : smallest_values)
     {
-        auto & hash_table = table[hash & register_id_mask];
+        auto & hash_table = table[hash & register_id_mask].hashes;
         if (hash_table.size() < sketch_size)
             hash_table.push_back(hash >> 4);
     }
 }
 
 bool minhashes::is_valid() const
 {
-    return table.size() == num_sketches
-        && std::ranges::all_of(table,
-                               [](auto const & minHash_sketch)
-                               {
-                                   return minHash_sketch.size() == sketch_size;
-                               });
+    return table.size() == num_sketches && std::ranges::all_of(table, &minhash_sketch::is_valid);
 }
 
 void minhashes::fill_incomplete_sketches(std::span<uint64_t> const & more_smallest_values)
@@ -47,7 +42,7 @@ void minhashes::fill_incomplete_sketches(std::span<uint64_t> const & more_smalle
 
     for (uint64_t const hash : more_smallest_values)
     {
-        auto & hash_table = table[hash & register_id_mask];
+        auto & hash_table = table[hash & register_id_mask].hashes;
         assert(std::ranges::find(hash_table, hash >> 4) == hash_table.end()); // hashes should be unique
         if (hash_table.size() < sketch_size)
             hash_table.push_back(hash >> 4);

test/unit/hibf/sketch/compute_sketches_test.cpp

@@ -22,7 +22,6 @@
 #include <hibf/sketch/minhashes.hpp>        // for minhashes
 #include <hibf/sketch/wyhash.hpp>           // for wyhash
 #include <hibf/test/expect_range_eq.hpp>    // for expect_range_eq, EXPECT_RANGE_EQ
-#include <hibf/test/expect_throw_msg.hpp>   // for EXPECT_THROW_MSG
 
 static constexpr uint64_t max_hash = std::numeric_limits<uint64_t>::max();
 
@@ -93,9 +92,9 @@ class compute_sketches_test : public ::testing::Test
             size_t sid{};
             for (auto & sketch : minhash.table)
             {
-                ASSERT_EQ(sketch.size(), seqan::hibf::sketch::minhashes::sketch_size)
+                ASSERT_EQ(sketch.hashes.size(), seqan::hibf::sketch::minhashes::sketch_size)
                     << "ub_id: " << ub_id << "sid: " << sid;
-                EXPECT_RANGE_EQ(sketch, expected.table[sid]) << "ub_id: " << ub_id << "sid: " << sid;
+                EXPECT_RANGE_EQ(sketch.hashes, expected.table[sid].hashes) << "ub_id: " << ub_id << "sid: " << sid;
                 ++sid;
             }
             ++ub_id;
@@ -124,13 +123,18 @@ TEST_F(compute_sketches_test, too_few_hashes)
     this->config.sketch_bits = 12;
     this->config.input_fn = [&](size_t const, seqan::hibf::insert_iterator it)
     {
-        for (size_t i = 0; i < 100; ++i)
+        for (size_t i = 0; i < 10; ++i)
             it = i;
     };
 
-    EXPECT_THROW(
-        seqan::hibf::sketch::compute_sketches(this->config, this->hyperloglog_sketches, this->minhash_sketches),
-        std::runtime_error);
+    seqan::hibf::sketch::compute_sketches(this->config, this->hyperloglog_sketches, this->minhash_sketches);
+
+    ASSERT_EQ(this->minhash_sketches.size(), 1); // only one user bin
+
+    for (auto const & sketch : this->minhash_sketches[0].table)
+    {
+        ASSERT_FALSE(sketch.is_valid()) << "It by chance is valid";
+    }
 }
 
 TEST_F(compute_sketches_test, minHash_refill)
@@ -159,17 +163,20 @@ TEST_F(compute_sketches_test, minHash_refill)
 
 TEST_F(compute_sketches_test, too_few_hashes_in_one_register)
 {
-    // 960 kmers, but only 31 hashes in register 15, so sketch 15 cannot be filled.
+    // 960 kmers, but only 31 hashes in register 15, so sketch 15 cannot be filled and stays incomplete.
     // With fewer than 1000 kmers, the old code kept max() placeholders in the heap. These filled sketch 15 with
-    // max() >> 4, and it did not throw. If the throw is removed, check that sketch 15 holds only the 31 hashes instead.
+    // max() >> 4.
     this->config.number_of_user_bins = 1;
     this->config.input_fn = [](size_t const, seqan::hibf::insert_iterator it)
     {
         sparse_register_15(990u, max_hash / 2u, it);
     };
 
-    EXPECT_THROW_MSG(
-        seqan::hibf::sketch::compute_sketches(this->config, this->hyperloglog_sketches, this->minhash_sketches),
-        std::runtime_error,
-        "Not enough kmers (960) to get 640 hashes for all minHash sketches.");
+    seqan::hibf::sketch::compute_sketches(this->config, this->hyperloglog_sketches, this->minhash_sketches);
+
+    ASSERT_EQ(this->minhash_sketches.size(), 1u);
+    auto const & table = this->minhash_sketches[0].table;
+    ASSERT_EQ(table.size(), seqan::hibf::sketch::minhashes::num_sketches);
+    EXPECT_FALSE(this->minhash_sketches[0].is_valid());
+    EXPECT_EQ(table[15].hashes.size(), 31u); // no max() >> 4 placeholders
 }

test/unit/hibf/sketch/minhashes_test.cpp

@@ -76,8 +76,9 @@ TEST(minhashes_test, ctor_sorted_list)
     size_t sid{};
     for (auto & sketch : sketches.table)
     {
-        ASSERT_EQ(sketch.size(), seqan::hibf::sketch::minhashes::sketch_size) << "sid: " << sid;
-        EXPECT_RANGE_EQ(sketch, std::views::iota(0u, seqan::hibf::sketch::minhashes::sketch_size)) << "sid: " << sid;
+        ASSERT_EQ(sketch.hashes.size(), seqan::hibf::sketch::minhashes::sketch_size) << "sid: " << sid;
+        EXPECT_RANGE_EQ(sketch.hashes, std::views::iota(0u, seqan::hibf::sketch::minhashes::sketch_size))
+            << "sid: " << sid;
         ++sid;
     }
 }
@@ -103,8 +104,9 @@ TEST(minhashes_test, fill_incomplete_sketches)
     size_t sid{};
     for (auto & sketch : sketches.table)
     {
-        ASSERT_EQ(sketch.size(), seqan::hibf::sketch::minhashes::sketch_size) << "sid: " << sid;
-        EXPECT_RANGE_EQ(sketch, std::views::iota(0u, seqan::hibf::sketch::minhashes::sketch_size)) << "sid: " << sid;
+        ASSERT_EQ(sketch.hashes.size(), seqan::hibf::sketch::minhashes::sketch_size) << "sid: " << sid;
+        EXPECT_RANGE_EQ(sketch.hashes, std::views::iota(0u, seqan::hibf::sketch::minhashes::sketch_size))
+            << "sid: " << sid;
         ++sid;
     }
 }