[FEATURE] Allow invalid minhash sketches in LSH by skipping those bins

src/layout/lsh_distributed_ibf_layout.cpp

@@ -73,13 +73,21 @@ auto LSH_fill_hashtable(std::vector<Cluster> const & clusters,
         auto const & current = clusters[pos];
         assert(current.is_valid(pos));
 
+        // user bins with invalid sketches (e.g. too foo kmers) are not considered for clustering
+        // They will be left as single clusters and merged by cardinality later in the algorithm
+        if (current.size() == 1 && !minHash_sketches[current.contained_user_bins()[0]].is_valid())
+        {
+            ++processed_user_bins;
+            continue;
+        }
+
         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],
+            uint64_t const key = lsh_hash_the_sketch(minHash_sketches[user_bin_idx].table[current_sketch_index].hashes,
                                                      current_number_of_sketch_hashes);
             table[key].push_back(current.id()); // insert representative for all user bins
         }
@@ -133,7 +141,6 @@ std::vector<Cluster> very_similar_LSH_clustering(std::vector<seqan::hibf::sketch
 
     assert(!minHash_sketches.empty());
     assert(minHash_sketches[0].table.size() >= number_of_max_minHash_sketches);
-    assert(minHash_sketches[0].table[0].size() >= minHash_sketch_size);
     assert(number_of_user_bins <= minHash_sketches.size());
 
     // Initialise one cluster per user bin. Merging moves the user bins of a cluster into another cluster, e.g., if

test/api/layout/fast_layout_test.cpp

@@ -79,3 +79,46 @@ TEST(fast_layout_test, recursion)
     EXPECT_NO_THROW(hibf_layout.validate(config.hibf_config, ignore_unexpected_empty_bins_and_throw));
     EXPECT_GE(hibf_layout.number_of_levels(), 2u);
 }
+
+TEST(fast_layout_test, too_few_kmers_for_minHash)
+{
+    // User bins 0-5 have about 20'000 k-mers. User bins 6 and 7 have 200 and 500 k-mers, which is too few to fill
+    // all minHash sketches (num_sketches * sketch_size = 640 hashes).
+    // They are skipped during LSH clustering and must still be assigned leter in the algorithm.
+    size_t const number_of_user_bins{8};
+    std::vector<size_t> const kmer_counts{20'000, 20'500, 19'500, 21'000, 19'000, 20'200, 200, 500};
+
+    chopper::configuration config{};
+    config.fast_layout = true;
+    config.hibf_config.tmax = 4;
+    config.hibf_config.threads = 2;
+    config.hibf_config.number_of_user_bins = number_of_user_bins;
+    config.hibf_config.disable_estimate_union = true; // also disables rearrangement
+    config.hibf_config.input_fn = [&kmer_counts](size_t const ub, seqan::hibf::insert_iterator it)
+    {
+        uint64_t const offset = static_cast<uint64_t>(ub) << 32;
+        for (uint64_t i = 0; i < kmer_counts[ub]; ++i)
+            it = scramble(offset + i);
+    };
+
+    std::vector<seqan::hibf::sketch::hyperloglog> sketches{};
+    std::vector<seqan::hibf::sketch::minhashes> minHash_sketches{};
+    seqan::hibf::sketch::compute_sketches(config.hibf_config, sketches, minHash_sketches);
+
+    ASSERT_EQ(minHash_sketches.size(), number_of_user_bins);
+    for (size_t ub = 0; ub < 6; ++ub)
+        ASSERT_TRUE(minHash_sketches[ub].is_valid()) << "ub: " << ub;
+    ASSERT_FALSE(minHash_sketches[6].is_valid());
+    ASSERT_FALSE(minHash_sketches[7].is_valid());
+
+    std::vector<size_t> cardinalities(number_of_user_bins);
+    for (size_t i = 0; i < number_of_user_bins; ++i)
+        cardinalities[i] = sketches[i].estimate();
+
+    std::vector<size_t> positions = seqan::hibf::iota_vector(number_of_user_bins);
+    seqan::hibf::layout::layout hibf_layout{};
+    chopper::layout::fast_layout(config, positions, cardinalities, sketches, minHash_sketches, hibf_layout);
+
+    // validate() also reports missing and duplicate user bins.
+    EXPECT_NO_THROW(hibf_layout.validate(config.hibf_config, ignore_unexpected_empty_bins_and_throw));
+}

test/cli/cli_chopper_basic_test.cpp

@@ -95,50 +95,6 @@ TEST_F(cli_test, chopper_cmd_kmer_bigger_than_window)
     EXPECT_EQ(result.err, std::string{"[ERROR] The k-mer size cannot be bigger than the window size.\n"});
 }
 
-TEST_F(cli_test, chopper_user_bin_with_few_kmers)
-{
-    seqan3::test::tmp_directory tmp_dir{};
-    std::filesystem::path const input_filename{tmp_dir.path() / "data.filenames"};
-    std::filesystem::path const short_filename{tmp_dir.path() / "short.fa"};
-    std::filesystem::path const layout_filename{tmp_dir.path() / "output.binning"};
-
-    // 200 bases give 182 k-mers, too few to fill the MinHash sketches that only the fast layout needs.
-    {
-        std::ofstream fout{short_filename};
-        fout << ">short\n";
-        for (size_t i = 0; i < 200; ++i)
-            fout << "ACGT"[(i * 7 + i / 3) % 4];
-        fout << '\n';
-    }
-
-    {
-        std::ofstream fout{input_filename};
-        fout << data("seq1.fa").string() << '\n' << short_filename.string() << '\n';
-    }
-
-    cli_test_result result =
-        execute_app("chopper", "--input", input_filename.c_str(), "--tmax", "64", "--output", layout_filename.c_str());
-
-    EXPECT_EQ(result.exit_code, 0);
-    EXPECT_EQ(result.out, std::string{});
-    EXPECT_EQ(result.err, std::string{});
-    EXPECT_TRUE(std::filesystem::exists(layout_filename));
-
-    // The fast layout needs MinHash sketches and fails.
-    result = execute_app("chopper",
-                         "--fast-layout",
-                         "--input",
-                         input_filename.c_str(),
-                         "--tmax",
-                         "64",
-                         "--output",
-                         layout_filename.c_str());
-
-    EXPECT_NE(result.exit_code, 0);
-    EXPECT_EQ(result.out, std::string{});
-    EXPECT_TRUE(result.err.starts_with("[ERROR] Not enough kmers")) << result.err;
-}
-
 TEST_F(cli_test, chopper_fast_layout_with_determine_best_tmax)
 {
     cli_test_result result = execute_app("chopper",