Merge pull request #479 from ksahlin/revert-basehash

marcelm · web-flow · commit c0695b042640 · 2025-01-29T23:11:52.000+01:00
Revert PR #472 "Always pick first syncmer as base"
diff --git a/CHANGES.md b/CHANGES.md
@@ -2,9 +2,9 @@
 
 ## development version
 
-* #476: Significantly improve speed and accuracy by enabling by default a new
-  variant of multi-context seeds: When no regular seeds - which consist
-  of two strobes - can be found for the entire query, strobealign attempts to find
+* #476: Improve accuracy by enabling (by default) a variant of multi-context
+  seeds: When no regular seeds - which consist of two strobes - can be found
+  for the *entire* query, strobealign now always attempts to find
   single-strobe ("partial") seeds.
   The `--mcs` option is still available for now. It is a bit slower, but
   slightly more accurate.
diff --git a/src/index.cpp b/src/index.cpp
@@ -309,7 +309,8 @@ void StrobemerIndex::assign_randstrobes(size_t ref_index, size_t offset) {
                 randstrobe.hash,
                 randstrobe.strobe1_pos,
                 static_cast<uint32_t>(ref_index),
-                static_cast<uint8_t>(randstrobe.strobe2_pos - randstrobe.strobe1_pos)
+                static_cast<uint8_t>(randstrobe.strobe2_pos - randstrobe.strobe1_pos),
+                randstrobe.first_strobe_is_main
             };
         }
         chunk.clear();
diff --git a/src/index.hpp b/src/index.hpp
@@ -100,7 +100,7 @@ struct StrobemerIndex {
 
         auto pos = std::lower_bound(randstrobes.begin() + position_start,
                                     randstrobes.begin() + position_end,
-                                    RefRandstrobe{key, 0, 0, 0},
+                                    RefRandstrobe{key, 0, 0, 0, 0},
                                     cmp);
         if ((pos->hash() & hash_mask) == masked_key) return pos - randstrobes.begin();
         return end();
@@ -122,6 +122,10 @@ struct StrobemerIndex {
         }
     }
 
+    bool first_strobe_is_main(bucket_index_t position) const {
+        return randstrobes[position].first_strobe_is_main();
+    }
+    
     bool is_filtered(bucket_index_t position) const {
         return get_hash(position) == get_hash(position + filter_cutoff);
     }
@@ -141,6 +145,9 @@ struct StrobemerIndex {
     std::pair<int, int> strobe_extent_partial(bucket_index_t position) const {
         // Construct the match from the strobe that was selected as the main part of the hash
         int ref_start = get_strobe1_position(position);
+        if (!first_strobe_is_main(position)) {
+            ref_start += strobe2_offset(position);
+        }
         return {ref_start, ref_start + k()};
     }
 
@@ -198,7 +205,7 @@ struct StrobemerIndex {
 
         auto pos = std::upper_bound(randstrobes.begin() + position,
                                     randstrobes.begin() + position_end,
-                                    RefRandstrobe{key, 0, 0, 0},
+                                    RefRandstrobe{key, 0, 0, 0, 0},
                                     cmp);
         return (pos - randstrobes.begin() - 1) - position + 1;
     }
diff --git a/src/nam.cpp b/src/nam.cpp
@@ -13,6 +13,19 @@ bool operator==(const Match& lhs, const Match& rhs) {
     return (lhs.query_start == rhs.query_start) && (lhs.query_end == rhs.query_end) && (lhs.ref_start == rhs.ref_start) && (lhs.ref_end == rhs.ref_end);
 }
 
+/*
+ * A partial hit is a hit where not the full randstrobe hash could be found in
+ * the index but only the "main" hash (only the first aux_len bits).
+ */
+struct PartialHit {
+    randstrobe_hash_t hash;
+    unsigned int start;  // position in strobemer index
+    bool is_reverse;
+    bool operator==(const PartialHit& rhs) const {
+        return (hash == rhs.hash) && (start == rhs.start) && (is_reverse == rhs.is_reverse);
+    }
+};
+
 inline void add_to_matches_map_full(
     robin_hood::unordered_map<unsigned int, std::vector<Match>>& matches_map,
     int query_start,
@@ -34,6 +47,25 @@ inline void add_to_matches_map_full(
     }
 }
 
+/*
+This function often produces the same Match multiple times.
+This happens when multiple randstrobes involve a syncmer that has a low hash
+value (and which therefore gets chosen as main hash multiple times).
+
+For example, assume we have syncmers starting at positions 14, 21, 39, 51
+and randstrobes like this:
+
+14, 21, 39, 51
+14------39
+    21--39
+        39--51
+
+If 39 now has the lowest hash value, it will be chosen as the main hash for all
+three randstrobes. Thus, when looking up that main hash, all three randstrobes
+are found. After the start coordinate of the randstrobe is converted to the
+start coordinate of the single syncmer (using strobe_extent_partial()), they all
+result in a hit that starts at 39 (and has length k).
+*/
 inline void add_to_matches_map_partial(
     robin_hood::unordered_map<unsigned int, std::vector<Match>>& matches_map,
     int query_start,
@@ -191,6 +223,10 @@ std::tuple<float, int, std::vector<Nam>> find_nams(
     const StrobemerIndex& index,
     bool use_mcs
 ) {
+    std::vector<PartialHit> partial_queried; // TODO: is a small set more efficient than linear search in a small vector?
+    if (use_mcs) {
+        partial_queried.reserve(10);
+    }
     std::array<robin_hood::unordered_map<unsigned int, std::vector<Match>>, 2> matches_map;
     matches_map[0].reserve(100);
     matches_map[1].reserve(100);
@@ -207,15 +243,22 @@ std::tuple<float, int, std::vector<Nam>> find_nams(
             add_to_matches_map_full(matches_map[q.is_revcomp], q.start, q.end, index, position);
         }
         else if (use_mcs) {
+            PartialHit ph{q.hash & index.get_main_hash_mask(), q.partial_start, q.is_revcomp};
+            if (std::find(partial_queried.begin(), partial_queried.end(), ph) != partial_queried.end()) {
+                // already queried
+                continue;
+            }
             size_t partial_pos = index.find_partial(q.hash);
             if (partial_pos != index.end()) {
                 total_hits++;
                 if (index.is_partial_filtered(partial_pos)) {
+                    partial_queried.push_back(ph);
                     continue;
                 }
                 nr_good_hits++;
-                add_to_matches_map_partial(matches_map[q.is_revcomp], q.start, q.start + index.k(), index, partial_pos);
+                add_to_matches_map_partial(matches_map[q.is_revcomp], q.partial_start, q.partial_end, index, partial_pos);
             }
+            partial_queried.push_back(ph);
         }
     }
 
@@ -263,6 +306,8 @@ std::pair<int, std::vector<Nam>> find_nams_rescue(
                 < std::tie(rhs.count, rhs.query_start, rhs.query_end);
         }
     };
+    std::vector<PartialHit> partial_queried;  // TODO: is a small set more efficient than linear search in a small vector?
+    partial_queried.reserve(10);
     std::array<robin_hood::unordered_map<unsigned int, std::vector<Match>>, 2> matches_map;
     std::array<std::vector<RescueHit>, 2> hits;
     matches_map[0].reserve(100);
@@ -278,12 +323,18 @@ std::pair<int, std::vector<Nam>> find_nams_rescue(
             hits[qr.is_revcomp].push_back(rh);
         }
         else if (use_mcs) {
+            PartialHit ph = {qr.hash & index.get_main_hash_mask(), qr.partial_start, qr.is_revcomp};
+            if (std::find(partial_queried.begin(), partial_queried.end(), ph) != partial_queried.end()) {
+                // already queried
+                continue;
+            }
             size_t partial_pos = index.find_partial(qr.hash);
             if (partial_pos != index.end()) {
                 unsigned int partial_count = index.get_count_partial(partial_pos);
-                RescueHit rh{partial_pos, partial_count, qr.start, qr.start + index.k(), true};
+                RescueHit rh{partial_pos, partial_count, qr.partial_start, qr.partial_end, true};
                 hits[qr.is_revcomp].push_back(rh);
             }
+            partial_queried.push_back(ph);
         }
     }
 
diff --git a/src/randstrobes.cpp b/src/randstrobes.cpp
@@ -49,7 +49,8 @@ static inline syncmer_hash_t syncmer_smer_hash(uint64_t packed) {
  * This function combines two individual syncmer hashes into a single hash
  * for the randstrobe.
  *
- * One syncmer is designated as the "main", the other is the "auxiliary".
+ * The syncmer with the smaller hash is designated as the "main", the other is
+ * the "auxiliary".
  * The combined hash is obtained by setting the top bits to the bits of
  * the main hash and the bottom bits to the bits of the auxiliary
  * hash. Since entries in the index are sorted by randstrobe hash, this allows
@@ -58,6 +59,10 @@ static inline syncmer_hash_t syncmer_smer_hash(uint64_t packed) {
 static inline randstrobe_hash_t randstrobe_hash(
     syncmer_hash_t hash1, syncmer_hash_t hash2, randstrobe_hash_t main_hash_mask
 ) {
+    // Make the function symmetric
+    if (hash1 > hash2) {
+        std::swap(hash1, hash2);
+    }
     return ((hash1 & main_hash_mask) | (hash2 & ~main_hash_mask)) & RANDSTROBE_HASH_MASK;
 }
 
@@ -141,7 +146,7 @@ std::vector<Syncmer> canonical_syncmers(
 
 std::ostream& operator<<(std::ostream& os, const Randstrobe& randstrobe) {
     os << "Randstrobe(hash=" << randstrobe.hash << ", strobe1_pos=" << randstrobe.strobe1_pos << ", strobe2_pos="
-       << randstrobe.strobe2_pos << ")";
+       << randstrobe.strobe2_pos << ", first_strobe_is_main=" << randstrobe.first_strobe_is_main << ")";
     return os;
 }
 
@@ -156,10 +161,12 @@ std::ostream& operator<<(std::ostream& os, const QueryRandstrobe& randstrobe) {
 }
 
 Randstrobe make_randstrobe(Syncmer strobe1, Syncmer strobe2, randstrobe_hash_t main_hash_mask) {
+    bool first_strobe_is_main = strobe1.hash < strobe2.hash;
     return Randstrobe{
         randstrobe_hash(strobe1.hash, strobe2.hash, main_hash_mask),
         static_cast<uint32_t>(strobe1.position),
-        static_cast<uint32_t>(strobe2.position)
+        static_cast<uint32_t>(strobe2.position),
+        first_strobe_is_main
     };
 }
 
@@ -239,9 +246,11 @@ std::vector<QueryRandstrobe> randstrobes_query(const std::string_view seq, const
     RandstrobeIterator randstrobe_fwd_iter{syncmers, parameters.randstrobe};
     while (randstrobe_fwd_iter.has_next()) {
         auto randstrobe = randstrobe_fwd_iter.next();
+        const unsigned int partial_start = randstrobe.first_strobe_is_main ? randstrobe.strobe1_pos : randstrobe.strobe2_pos;
         randstrobes.push_back(
             QueryRandstrobe {
-                randstrobe.hash, randstrobe.strobe1_pos, randstrobe.strobe2_pos + parameters.syncmer.k, false
+                randstrobe.hash, randstrobe.strobe1_pos, randstrobe.strobe2_pos + parameters.syncmer.k,
+                partial_start, partial_start + parameters.syncmer.k,  false
             }
         );
     }
@@ -262,9 +271,12 @@ std::vector<QueryRandstrobe> randstrobes_query(const std::string_view seq, const
     RandstrobeIterator randstrobe_rc_iter{syncmers, parameters.randstrobe};
     while (randstrobe_rc_iter.has_next()) {
         auto randstrobe = randstrobe_rc_iter.next();
+        bool first_strobe_is_main = randstrobe.first_strobe_is_main;
+        const unsigned int partial_start = first_strobe_is_main ? randstrobe.strobe1_pos : randstrobe.strobe2_pos;
         randstrobes.push_back(
             QueryRandstrobe {
-                randstrobe.hash, randstrobe.strobe1_pos, randstrobe.strobe2_pos + parameters.syncmer.k, true
+                randstrobe.hash, randstrobe.strobe1_pos, randstrobe.strobe2_pos + parameters.syncmer.k,
+                partial_start, partial_start + parameters.syncmer.k, true
             }
         );
     }
diff --git a/src/randstrobes.hpp b/src/randstrobes.hpp
@@ -17,20 +17,20 @@
 using syncmer_hash_t = uint64_t;
 using randstrobe_hash_t = uint64_t;
 
-static constexpr uint64_t RANDSTROBE_HASH_MASK = 0xFFFFFFFFFFFFFF00;
+static constexpr uint64_t RANDSTROBE_HASH_MASK = 0xFFFFFFFFFFFFFE00;
 
 struct RefRandstrobe {
 private:
-    // packed representation of hash and offset
+    // packed representation of hash, offset and first_strobe_is_main
     randstrobe_hash_t m_hash_offset_flag;
     uint32_t m_position;
     uint32_t m_ref_index;
 
 public:
     RefRandstrobe() : m_hash_offset_flag(0), m_position(0), m_ref_index(0) { }
 
-    RefRandstrobe(randstrobe_hash_t hash, uint32_t position, uint32_t ref_index, uint8_t offset)
-        : m_hash_offset_flag((hash & RANDSTROBE_HASH_MASK) | offset)
+    RefRandstrobe(randstrobe_hash_t hash, uint32_t position, uint32_t ref_index, uint8_t offset, bool first_strobe_is_main)
+        : m_hash_offset_flag((hash & RANDSTROBE_HASH_MASK) | (offset << 1) | first_strobe_is_main)
         , m_position(position)
         , m_ref_index(ref_index)
     { }
@@ -45,12 +45,16 @@ struct RefRandstrobe {
         return lhs < rhs;
     }
 
+    bool first_strobe_is_main() const {
+        return m_hash_offset_flag & 1;
+    }
+
     unsigned reference_index() const {
         return m_ref_index;
     }
 
     unsigned strobe2_offset() const {
-        return m_hash_offset_flag & 0xff;
+        return (m_hash_offset_flag >> 1) & 0xff;
     }
 
     randstrobe_hash_t hash() const {
@@ -68,6 +72,11 @@ struct QueryRandstrobe {
     randstrobe_hash_t hash;
     unsigned int start;
     unsigned int end;
+    /* Start and end of the main syncmer (relevant if the randstrobe couldn’t
+     * be found in the index and we fall back to a partial hit)
+     */
+    unsigned int partial_start;
+    unsigned int partial_end;
     bool is_revcomp;
 };
 
@@ -79,6 +88,7 @@ struct Randstrobe {
     randstrobe_hash_t hash;
     unsigned int strobe1_pos;
     unsigned int strobe2_pos;
+    bool first_strobe_is_main;
 
     bool operator==(const Randstrobe& other) const {
         return hash == other.hash && strobe1_pos == other.strobe1_pos && strobe2_pos == other.strobe2_pos;
@@ -176,7 +186,7 @@ class RandstrobeGenerator {
     { }
 
     Randstrobe next();
-    Randstrobe end() const { return Randstrobe{0, 0, 0}; }
+    Randstrobe end() const { return Randstrobe{0, 0, 0, false}; }
 
 private:
     SyncmerIterator syncmer_iterator;
diff --git a/tests/baseline-commit.txt b/tests/baseline-commit.txt
@@ -1 +1 @@
-acc4cffe5ac2c4db266c58d00b7b6462c6b4189c
+b63002e80213afd97517d47162ba5e51f2bed56e
diff --git a/tests/phix.mcs.se.paf b/tests/phix.mcs.se.paf
@@ -1,11 +1,11 @@
-SRR1377138.32	301	2	293	+	NC_001422.1	5386	1434	1725	49	291	255
-SRR1377138.33	301	2	297	+	NC_001422.1	5386	3818	4113	51	295	255
-SRR1377138.34	301	33	299	-	NC_001422.1	5386	844	1110	43	266	255
-SRR1377138.35	301	5	298	-	NC_001422.1	5386	4041	4334	49	293	255
-SRR1377138.36	301	3	301	+	NC_001422.1	5386	4997	5295	53	298	255
-SRR1377138.37	301	2	299	-	NC_001422.1	5386	794	1091	43	297	255
-SRR1377138.38	301	32	284	-	NC_001422.1	5386	4971	5223	44	252	255
-SRR1377138.39/1	301	1	295	-	NC_001422.1	5386	1791	2085	50	294	255
-SRR1377138.40	301	4	293	-	NC_001422.1	5386	3020	3309	48	289	255
-rescuable.42	301	4	293	-	NC_001422.1	5386	3020	3309	48	289	255
-not.rescuable	301	4	293	-	NC_001422.1	5386	3020	3309	48	289	255
+SRR1377138.32	301	2	293	+	NC_001422.1	5386	1434	1725	47	291	255
+SRR1377138.33	301	2	297	+	NC_001422.1	5386	3818	4113	47	295	255
+SRR1377138.34	301	33	299	-	NC_001422.1	5386	844	1110	40	266	255
+SRR1377138.35	301	5	298	-	NC_001422.1	5386	4041	4334	47	293	255
+SRR1377138.36	301	3	301	+	NC_001422.1	5386	4997	5295	51	298	255
+SRR1377138.37	301	2	299	-	NC_001422.1	5386	794	1091	42	297	255
+SRR1377138.38	301	32	284	-	NC_001422.1	5386	4971	5223	42	252	255
+SRR1377138.39/1	301	1	295	-	NC_001422.1	5386	1791	2085	48	294	255
+SRR1377138.40	301	4	293	-	NC_001422.1	5386	3020	3309	44	289	255
+rescuable.42	301	4	293	-	NC_001422.1	5386	3020	3309	44	289	255
+not.rescuable	301	4	293	-	NC_001422.1	5386	3020	3309	44	289	255
diff --git a/tests/test_randstrobes.cpp b/tests/test_randstrobes.cpp
@@ -39,27 +39,15 @@ TEST_CASE("RefRandstrobe constructor") {
     randstrobe_hash_t hash = 0x1234567890ABCDEF & RANDSTROBE_HASH_MASK;
     uint32_t position = ~0u;
     uint32_t ref_index = RefRandstrobe::max_number_of_references - 1;
-    SUBCASE("one") {
-        uint8_t offset = 255;
-        RefRandstrobe rr{hash, position, ref_index, offset};
-
-        CHECK(rr.hash() == hash);
-        CHECK(rr.position() == position);
-        CHECK(rr.reference_index() == ref_index);
-        CHECK(rr.strobe2_offset() == offset);
-    }
-
-    SUBCASE("two") {
-        uint8_t offset = 0;
-        RefRandstrobe rr{hash, position, ref_index, offset};
-
-        CHECK(rr.hash() == hash);
-        CHECK(rr.position() == position);
-        CHECK(rr.reference_index() == ref_index);
-        CHECK(rr.strobe2_offset() == offset);
-    }
-
-
+    uint8_t offset = 255;
+    bool first_strobe_is_main = true;
+    RefRandstrobe rr{hash, position, ref_index, offset, first_strobe_is_main};
+
+    CHECK(rr.hash() == hash);
+    CHECK(rr.position() == position);
+    CHECK(rr.reference_index() == ref_index);
+    CHECK(rr.strobe2_offset() == offset);
+    CHECK(rr.first_strobe_is_main() == first_strobe_is_main);
 }
 
 TEST_CASE("SyncmerIterator") {

Original file line number	Diff line number	Diff line change
`@@ -1 +1 @@`
`1`		`-acc4cffe5ac2c4db266c58d00b7b6462c6b4189c`
	`1`	`+b63002e80213afd97517d47162ba5e51f2bed56e`