Experiment: Asymmetric randstrobe hashes

src/aln.cpp

@@ -972,13 +972,18 @@ void align_or_map_paired(
         // Find NAMs
         Timer nam_timer;
         auto [nonrepetitive_fraction, nams] = find_nams(query_randstrobes, index);
+
         statistics.tot_find_nams += nam_timer.duration();
 
         if (map_param.rescue_level > 1) {
             Timer rescue_timer;
             if (nams.empty() || nonrepetitive_fraction < 0.7) {
                 nams = find_nams_rescue(query_randstrobes, index, map_param.rescue_cutoff);
                 details[is_revcomp].nam_rescue = true;
+                if (nams.empty()) {
+                    query_randstrobes = randstrobes_query_rescue(record.seq, index_parameters);
+                    std::tie(nonrepetitive_fraction, nams) = find_nams(query_randstrobes, index);
+                }
             }
             statistics.tot_time_rescue += rescue_timer.duration();
         }
@@ -1092,13 +1097,18 @@ void align_or_map_single(
     // Find NAMs
     Timer nam_timer;
     auto [nonrepetitive_fraction, nams] = find_nams(query_randstrobes, index);
+
     statistics.tot_find_nams += nam_timer.duration();
 
     if (map_param.rescue_level > 1) {
         Timer rescue_timer;
         if (nams.empty() || nonrepetitive_fraction < 0.7) {
             details.nam_rescue = true;
             nams = find_nams_rescue(query_randstrobes, index, map_param.rescue_cutoff);
+            if (nams.empty()) {
+                query_randstrobes = randstrobes_query_rescue(record.seq, index_parameters);
+                std::tie(nonrepetitive_fraction, nams) = find_nams(query_randstrobes, index);
+            }
         }
         statistics.tot_time_rescue += rescue_timer.duration();
     }

src/randstrobes.cpp

@@ -45,8 +45,10 @@ static inline syncmer_hash_t syncmer_smer_hash(uint64_t packed) {
     return xxh64(packed);
 }
 
-static inline randstrobe_hash_t randstrobe_hash(syncmer_hash_t hash1, syncmer_hash_t hash2) {
-    return hash1 + hash2;
+static inline randstrobe_hash_t randstrobe_hash(syncmer_hash_t syncmer1_hash, syncmer_hash_t syncmer2_hash) {
+    // alternative:
+    // (syncmer1_hash << 1) ^ syncmer2_hash
+    return 2 * syncmer1_hash - syncmer2_hash;
 }
 
 std::ostream& operator<<(std::ostream& os, const Syncmer& syncmer) {
@@ -251,3 +253,61 @@ QueryRandstrobeVector randstrobes_query(const std::string_view seq, const IndexP
     }
     return randstrobes;
 }
+
+/*
+ * Generate all possible strobes using a pre-computed vector of syncmers
+ */
+void randstrobes_rescue(const std::vector<Syncmer>& syncmers, std::vector<QueryRandstrobe>& randstrobes, bool is_rc, RandstrobeParameters parameters, int k) {
+    const unsigned w_min = parameters.w_min;
+    const unsigned w_max = parameters.w_max;
+    const unsigned int max_dist = parameters.max_dist;
+    if (w_min > w_max) {
+        throw std::invalid_argument("w_min is greater than w_max");
+    }
+    unsigned int strobe1_index = 0;
+
+    while (strobe1_index + w_min < syncmers.size()) {
+        unsigned int w_end = std::min(static_cast<size_t>(strobe1_index + w_max), syncmers.size() - 1);
+        auto strobe1 = syncmers[strobe1_index];
+        auto max_position = strobe1.position + max_dist;
+        unsigned int w_start = strobe1_index + w_min;
+        Syncmer strobe2 = strobe1;
+
+        for (auto i = w_start; i <= w_end && syncmers[i].position <= max_position; i++) {
+            strobe2 = syncmers[i];
+            randstrobes.push_back(
+                QueryRandstrobe{
+                    randstrobe_hash(strobe1.hash, strobe2.hash), static_cast<uint32_t>(strobe1.position), static_cast<uint32_t>(strobe2.position + k), is_rc
+                }
+            );
+        }
+        strobe1_index++;
+    }
+}
+
+/*
+ * Generate randstrobes for a query sequence and its reverse complement.
+ */
+QueryRandstrobeVector randstrobes_query_rescue(const std::string_view seq, const IndexParameters& parameters) {
+    QueryRandstrobeVector randstrobes;
+    if (seq.length() < parameters.randstrobe.w_max) {
+        return randstrobes;
+    }
+
+    // Generate syncmers for the forward sequence
+    auto syncmers = canonical_syncmers(seq, parameters.syncmer);
+    if (syncmers.empty()) {
+        return randstrobes;
+    }
+
+    randstrobes_rescue(syncmers, randstrobes, false, parameters.randstrobe, parameters.syncmer.k);
+
+    std::reverse(syncmers.begin(), syncmers.end());
+    for (size_t i = 0; i < syncmers.size(); i++) {
+        syncmers[i].position = seq.length() - syncmers[i].position - parameters.syncmer.k;
+    }
+
+    randstrobes_rescue(syncmers, randstrobes, true, parameters.randstrobe, parameters.syncmer.k);
+
+    return randstrobes;
+}

src/randstrobes.hpp

@@ -69,6 +69,7 @@ std::ostream& operator<<(std::ostream& os, const QueryRandstrobe& randstrobe);
 using QueryRandstrobeVector = std::vector<QueryRandstrobe>;
 
 QueryRandstrobeVector randstrobes_query(const std::string_view seq, const IndexParameters& parameters);
+QueryRandstrobeVector randstrobes_query_rescue(const std::string_view seq, const IndexParameters& parameters);
 
 struct Randstrobe {
     randstrobe_hash_t hash;