Let --aux-len work correctly

This includes some refactoring. In particular, instead of aux_len, we now
use a main_hash_mask and instead of right-shifting by aux_len, we use the
mask to get only the relevant part of the hash.

Author: marcelm

CMakeLists.txt

@@ -132,6 +132,7 @@ add_executable(test-strobealign
   tests/test_aligner.cpp
   tests/test_cigar.cpp
   tests/test_randstrobes.cpp
+  tests/test_indexparameters.cpp
 )
 target_link_libraries(test-strobealign salib)
 target_include_directories(test-strobealign PUBLIC src/ ext/ ${PROJECT_BINARY_DIR})

src/arguments.hpp

@@ -26,7 +26,7 @@ struct SeedingArguments {
             "results with non default values.", {'s'}}
         , bits{parser, "INT", "No. of top bits of hash to use as bucket indices (8-31)"
             "[determined from reference size]", {'b'}}
-        , aux_len{parser, "INT", "No. of bits to use from secondary strobe hash [24]", {"aux-len"}}
+        , aux_len{parser, "INT", "No. of bits to use from secondary strobe hash [17]", {"aux-len"}}
     {
     }
     args::ArgumentParser& parser;

src/cmdline.hpp

@@ -52,7 +52,7 @@ struct CommandLineOptions {
     int u { 7 };
     int s { 16 };
     int c { 8 };
-    int aux_len{26};
+    int aux_len{17};
 
     // Alignment
     int A { 2 };

src/dumpstrobes.cpp

@@ -101,7 +101,7 @@ int run_dumpstrobes(int argc, char **argv) {
     }
 
     // Seeding
-    int r{150}, k{20}, s{16}, c{8}, l{1}, u{7}, aux_len{26};
+    int r{150}, k{20}, s{16}, c{8}, l{1}, u{7}, aux_len{17};
     int max_seed_len{};
 
     bool k_set{false}, s_set{false}, c_set{false}, max_seed_len_set{false}, l_set{false}, u_set{false};

src/index.hpp

@@ -59,25 +59,23 @@ struct StrobemerIndex {
 
     // Find first entry that matches the given key
     size_t find_full(randstrobe_hash_t key) const {
-        return find(key, 0);
+        return find(key, RANDSTROBE_HASH_MASK);
     }
 
     /*
      * Find the first entry that matches the main hash (ignoring the aux_len
      * least significant bits)
      */
     size_t find_partial(randstrobe_hash_t key) const {
-        return find(key, parameters.randstrobe.aux_len);
+        return find(key, parameters.randstrobe.main_hash_mask);
     }
 
     /*
-     * Find first entry whose hash matches the given key, but ignore the
-     * b least significant bits
+     * Find first entry whose hash matches the given key. Mask both key and
+     * entry by hash_mask.
      */
-    size_t find(randstrobe_hash_t key, uint8_t b) const {
-        const unsigned int aux_len = b;
-        randstrobe_hash_t key_prefix = key >> aux_len;
-
+    size_t find(randstrobe_hash_t key, uint64_t hash_mask) const {
+        randstrobe_hash_t masked_key = key & hash_mask;
         constexpr int MAX_LINEAR_SEARCH = 4;
         const unsigned int top_N = key >> (64 - bits);
         bucket_index_t position_start = randstrobe_start_indices[top_N];
@@ -88,19 +86,20 @@ struct StrobemerIndex {
 
         if (position_end - position_start < MAX_LINEAR_SEARCH) {
             for ( ; position_start < position_end; ++position_start) {
-                if (randstrobes[position_start].hash() >> aux_len == key_prefix) return position_start;
-                if (randstrobes[position_start].hash() >> aux_len > key_prefix) return end();
+                if ((randstrobes[position_start].hash() & hash_mask) == masked_key) return position_start;
+                if ((randstrobes[position_start].hash() & hash_mask) > masked_key) return end();
             }
             return end();
         }
-        auto cmp = [&aux_len](const RefRandstrobe lhs, const RefRandstrobe rhs) {
-            return (lhs.hash() >> aux_len) < (rhs.hash() >> aux_len); };
+        auto cmp = [&hash_mask](const RefRandstrobe lhs, const RefRandstrobe rhs) {
+            return (lhs.hash() & hash_mask) < (rhs.hash() & hash_mask);
+        };
 
         auto pos = std::lower_bound(randstrobes.begin() + position_start,
                                     randstrobes.begin() + position_end,
                                     RefRandstrobe{key, 0, 0, 0, 0},
                                     cmp);
-        if (pos->hash() >> aux_len == key_prefix) return pos - randstrobes.begin();
+        if ((pos->hash() & hash_mask) == masked_key) return pos - randstrobes.begin();
         return end();
     }
 
@@ -114,7 +113,7 @@ struct StrobemerIndex {
 
     randstrobe_hash_t get_main_hash(bucket_index_t position) const {
         if (position < randstrobes.size()) {
-            return randstrobes[position].hash() >> parameters.randstrobe.aux_len;
+            return randstrobes[position].hash() & parameters.randstrobe.main_hash_mask;
         } else {
             return end();
         }
@@ -129,8 +128,7 @@ struct StrobemerIndex {
     }
 
     bool is_partial_filtered(bucket_index_t position) const {
-        const unsigned int shift = parameters.randstrobe.aux_len;
-        return (get_hash(position) >> shift) == (get_hash(position + partial_filter_cutoff) >> shift);
+        return get_main_hash(position) == get_main_hash(position + partial_filter_cutoff);
     }
 
     unsigned int get_strobe1_position(bucket_index_t position) const {
@@ -163,14 +161,14 @@ struct StrobemerIndex {
     }
 
     unsigned int get_count_full(bucket_index_t position) const {
-        return get_count(position, 0);
+        return get_count(position, RANDSTROBE_HASH_MASK);
     }
 
     unsigned int get_count_partial(bucket_index_t position) const {
-        return get_count(position, parameters.randstrobe.aux_len);
+        return get_count(position, parameters.randstrobe.main_hash_mask);
     }
 
-    unsigned int get_count(bucket_index_t position, uint8_t b) const {
+    unsigned int get_count(bucket_index_t position, uint64_t hash_mask) const {
         // For 95% of cases, the result will be small and a brute force search
         // is the best option. Once, we go over MAX_LINEAR_SEARCH, though, we
         // use a binary search to get the next position
@@ -182,27 +180,25 @@ struct StrobemerIndex {
         // seed with the given hash to yield the number of seeds with this hash.
 
         constexpr unsigned int MAX_LINEAR_SEARCH = 8;
-        const unsigned int aux_len = b;
 
         const auto key = randstrobes[position].hash();
-        randstrobe_hash_t key_prefix = key >> aux_len;
+        randstrobe_hash_t masked_key = key & hash_mask;
 
         const unsigned int top_N = key >> (64 - bits);
         bucket_index_t position_end = randstrobe_start_indices[top_N + 1];
         uint64_t count = 1;
 
         if (position_end - position < MAX_LINEAR_SEARCH) {
             for (bucket_index_t position_start = position + 1; position_start < position_end; ++position_start) {
-                if (randstrobes[position_start].hash() >> aux_len == key_prefix){
+                if ((randstrobes[position_start].hash() & hash_mask) == masked_key) {
                     count += 1;
-                }
-                else{
+                } else {
                     break;
                 }
             }
             return count;
         }
-        auto cmp = [&aux_len](const RefRandstrobe lhs, const RefRandstrobe rhs) {return (lhs.hash() >> aux_len) < (rhs.hash() >> aux_len); };
+        auto cmp = [&hash_mask](const RefRandstrobe lhs, const RefRandstrobe rhs) {return (lhs.hash() & hash_mask) < (rhs.hash() & hash_mask); };
 
         auto pos = std::upper_bound(randstrobes.begin() + position,
                                     randstrobes.begin() + position_end,
@@ -223,8 +219,8 @@ struct StrobemerIndex {
         return bits;
     }
 
-    int get_aux_len() const {
-        return parameters.randstrobe.aux_len;
+    uint64_t get_main_hash_mask() const {
+        return parameters.randstrobe.main_hash_mask;
     }
 
 private:

src/indexparameters.cpp

@@ -17,7 +17,7 @@ bool RandstrobeParameters::operator==(const RandstrobeParameters& other) const {
         && this->max_dist == other.max_dist
         && this->w_min == other.w_min
         && this->w_max == other.w_max
-        && this->aux_len == other.aux_len;
+        && this->main_hash_mask == other.main_hash_mask;
 }
 
 /* Pre-defined index parameters that work well for a certain
@@ -78,7 +78,7 @@ IndexParameters IndexParameters::from_read_length(int read_length, int k, int s,
     }
     int q = std::pow(2, c == DEFAULT ? default_c : c) - 1;
     if (aux_len == DEFAULT) {
-        aux_len = 26;
+        aux_len = 17;
     }
 
     return IndexParameters(canonical_read_length, k, s, l, u, q, max_dist, aux_len);
@@ -92,7 +92,7 @@ void IndexParameters::write(std::ostream& os) const {
     write_int_to_ostream(os, randstrobe.w_max);
     write_int_to_ostream(os, randstrobe.q);
     write_int_to_ostream(os, randstrobe.max_dist);
-    write_int_to_ostream(os, randstrobe.aux_len);
+    write_uint64_to_ostream(os, randstrobe.main_hash_mask);
 }
 
 IndexParameters IndexParameters::read(std::istream& is) {
@@ -105,8 +105,8 @@ IndexParameters IndexParameters::read(std::istream& is) {
     uint32_t w_max = read_int_from_istream(is);
     uint64_t q = read_int_from_istream(is);
     int max_dist = read_int_from_istream(is);
-    uint32_t aux_len = read_int_from_istream(is);
-    const RandstrobeParameters randstrobe_parameters{q, max_dist, w_min, w_max, aux_len};
+    uint64_t main_hash_mask = read_uint64_from_istream(is);
+    const RandstrobeParameters randstrobe_parameters{q, max_dist, w_min, w_max, main_hash_mask};
 
     return IndexParameters(canonical_read_length, syncmer_parameters, randstrobe_parameters);
 }
@@ -132,6 +132,26 @@ std::string IndexParameters::filename_extension() const {
     return sstream.str();
 }
 
+std::ostream& operator<<(std::ostream& os, const SyncmerParameters& parameters) {
+    os << "SyncmerParameters("
+        << "k=" << parameters.k
+        << ", s=" << parameters.s
+        << ", t_syncmer=" << parameters.t_syncmer
+        << ")";
+    return os;
+}
+
+std::ostream& operator<<(std::ostream& os, const RandstrobeParameters& parameters) {
+    os << "RandstrobeParameters("
+        << "q=" << parameters.q
+        << ", max_dist=" << parameters.max_dist
+        << ", w_min=" << parameters.w_min
+        << ", w_max=" << parameters.w_max
+        << ", main_hash_mask=0x" << std::hex << parameters.main_hash_mask << std::dec
+        << ")";
+    return os;
+}
+
 std::ostream& operator<<(std::ostream& os, const IndexParameters& parameters) {
     os << "IndexParameters("
         << "r=" << parameters.canonical_read_length
@@ -142,6 +162,7 @@ std::ostream& operator<<(std::ostream& os, const IndexParameters& parameters) {
         << ", max_dist=" << parameters.randstrobe.max_dist
         << ", w_min=" << parameters.randstrobe.w_min
         << ", w_max=" << parameters.randstrobe.w_max
+        << ", main_hash_mask=0x" << std::hex << parameters.randstrobe.main_hash_mask << std::dec
         << ")";
     return os;
 }

src/indexparameters.hpp

@@ -41,22 +41,21 @@ struct RandstrobeParameters {
     const int max_dist;
     const unsigned w_min;
     const unsigned w_max;
-    const unsigned aux_len;
+    const uint64_t main_hash_mask;
 
-    RandstrobeParameters(uint64_t q, int max_dist, unsigned w_min, unsigned w_max, unsigned aux_len)
+    RandstrobeParameters(uint64_t q, int max_dist, unsigned w_min, unsigned w_max, uint64_t main_hash_mask)
         : q(q)
         , max_dist(max_dist)
         , w_min(w_min)
         , w_max(w_max)
-        , aux_len(aux_len)
+        , main_hash_mask(main_hash_mask)
     {
-        verify();
     }
 
     bool operator==(const RandstrobeParameters& other) const;
 
 private:
-    void verify() const {
+    void verify(unsigned aux_len) const {
         if (max_dist > 255) {
             throw BadParameter("maximum seed length (-m <max_dist>) is larger than 255");
         }
@@ -78,11 +77,12 @@ class IndexParameters {
 
     static const int DEFAULT = std::numeric_limits<int>::min();
 
-    IndexParameters(size_t canonical_read_length, int k, int s, int l, int u, int q, int max_dist, int aux_len)
+    IndexParameters(size_t canonical_read_length, int k, int s, int l, int u, uint64_t q, int max_dist, int aux_len)
         : canonical_read_length(canonical_read_length)
         , syncmer(k, s)
-        , randstrobe(q, max_dist, std::max(0, k / (k - s + 1) + l), k / (k - s + 1) + u, aux_len)
+        , randstrobe(q, max_dist, std::max(0, k / (k - s + 1) + l), k / (k - s + 1) + u, ~0ul << (9 + aux_len))
     {
+        verify(aux_len);
     }
 
     IndexParameters(size_t canonical_read_length, SyncmerParameters syncmer, RandstrobeParameters randstrobe)
@@ -92,6 +92,12 @@ class IndexParameters {
     {
     }
 
+    void verify(unsigned aux_len) const {
+        if (aux_len > 27) {
+            throw BadParameter("aux_len must be less than 28");
+        }
+    }
+
     static IndexParameters from_read_length(
         int read_length, int k = DEFAULT, int s = DEFAULT, int l = DEFAULT, int u = DEFAULT, int c = DEFAULT, int max_seed_len = DEFAULT, int aux_len = DEFAULT);
     static IndexParameters read(std::istream& os);
@@ -101,6 +107,8 @@ class IndexParameters {
     bool operator!=(const IndexParameters& other) const { return !(*this == other); }
 };
 
+std::ostream& operator<<(std::ostream& os, const SyncmerParameters& parameters);
+std::ostream& operator<<(std::ostream& os, const RandstrobeParameters& parameters);
 std::ostream& operator<<(std::ostream& os, const IndexParameters& parameters);
 
 #endif

src/io.cpp

@@ -11,3 +11,15 @@ int32_t read_int_from_istream(std::istream& is) {
     is.read(reinterpret_cast<char*>(&val), sizeof(val));
     return val;
 }
+
+void write_uint64_to_ostream(std::ostream& os, uint64_t value) {
+    uint64_t val;
+    val = value;
+    os.write(reinterpret_cast<const char*>(&val), sizeof(val));
+}
+
+uint64_t read_uint64_from_istream(std::istream& is) {
+    uint64_t val;
+    is.read(reinterpret_cast<char*>(&val), sizeof(val));
+    return val;
+}

src/io.hpp

@@ -8,6 +8,9 @@
 void write_int_to_ostream(std::ostream& os, int32_t value);
 int32_t read_int_from_istream(std::istream& is);
 
+void write_uint64_to_ostream(std::ostream& os, uint64_t value);
+uint64_t read_uint64_from_istream(std::istream& is);
+
 /* Write a vector to an output stream, preceded by its length */
 template <typename T>
 void write_vector(std::ostream& os, const std::vector<T>& v) {

src/main.cpp

@@ -230,7 +230,8 @@ int run_strobealign(int argc, char **argv) {
         throw InvalidFasta("Too many reference sequences. Current maximum is " + std::to_string(RefRandstrobe::max_number_of_references));
     }
 
-    logger.debug() << "Auxiliary hash length: " << index_parameters.randstrobe.aux_len << "\n";
+    logger.debug() << "Auxiliary hash length: " << opt.aux_len << "\n";
+    logger.debug() << "Base hash mask: " << std::hex << index_parameters.randstrobe.main_hash_mask << std::dec << '\n';
     logger.info() << "Using multi-context seeds: " << (map_param.use_mcs ? "yes" : "no") << '\n';
     StrobemerIndex index(references, index_parameters, opt.bits);
     if (opt.use_index) {

src/nam.cpp

@@ -234,7 +234,7 @@ std::tuple<float, int, std::vector<Nam>> find_nams(
     int total_hits = 0;
     for (const auto &q : query_randstrobes) {
         size_t position = index.find_full(q.hash);
-        if (position != index.end()){
+        if (position != index.end()) {
             total_hits++;
             if (index.is_filtered(position)) {
                 continue;
@@ -243,7 +243,7 @@ std::tuple<float, int, std::vector<Nam>> find_nams(
             add_to_matches_map_full(matches_map[q.is_reverse], q.start, q.end, index, position);
         }
         else if (use_mcs) {
-            PartialHit ph{q.hash >> index.get_aux_len(), q.partial_start, q.is_reverse};
+            PartialHit ph{q.hash & index.get_main_hash_mask(), q.partial_start, q.is_reverse};
             if (std::find(partial_queried.begin(), partial_queried.end(), ph) != partial_queried.end()) {
                 // already queried
                 continue;
@@ -312,7 +312,7 @@ std::pair<int, std::vector<Nam>> find_nams_rescue(
             }
         }
         else if (use_mcs) {
-            PartialHit ph = {qr.hash >> index.get_aux_len(), qr.partial_start, qr.is_reverse};
+            PartialHit ph = {qr.hash & index.get_main_hash_mask(), qr.partial_start, qr.is_reverse};
             if (std::find(partial_queried.begin(), partial_queried.end(), ph) != partial_queried.end()) {
                 // already queried
                 continue;