From 6c53ae953be703e4c6e0430c0a31df372fc1c7b9 Mon Sep 17 00:00:00 2001 From: Giulio Ermanno Pibiri Date: Sat, 25 Jul 2026 16:18:06 +0000 Subject: [PATCH 01/14] one indexing modality: canonical minimizer at non-canonical density SSHash had two parsing modalities. The regular one picks, among the m-mers of a kmer, the one of smallest hash; the canonical one (--canonical) picks the smaller *value* between the best forward m-mer and the best reverse m-mer, so that a kmer and its reverse complement share a bucket and a lookup costs one probe instead of two. That second selection is made under an order unrelated to the hash, which is what made canonical mode expensive: comparing the outcomes of two independent minimizer processes under a third order costs a 4/3 factor in break density, hence ~1.28x the super-kmers and +0.5..0.9 bits/kmer. The 4/3 is the price of the order mismatch, not of canonicality. Keeping a single order recovers it. The minimizer is now the locus i minimizing h(kappa(i)), where kappa(i) = min(m-mer at i, rc(m-mer at i)) is the canonical m-mer at that locus, and the minimizer's value is kappa(i*): - kappa is invariant under reverse-complementation and the kappa-sequence of rc(x) is the reversal of that of x, so rc(x) selects the mirrored locus and literally the same m-mer -> same bucket, one probe; - the selection is still an argmin under one random order, so the density is the plain forward one, 2/(k-m+2) -> the super-kmer count of the regular modality. Since the canonical index is now no more expensive than the regular one and answers strictly more (it reports the orientation of the match), there is no trade-off left to expose: --canonical is gone and there is one build path, one lookup path, one streaming query. Measured on salmonella_enterica k=31 (m=13 / m=20): super-kmers bits/kmer build [musec] regular (old) 481390 / 501234 4.876 / 5.129 371711 / 630311 canonical (old) 619049 / 646108 5.557 / 5.799 497243 / 765721 new 482031 / 502192 4.946 / 5.128 402535 / 637544 and streaming query (best of 5, ms): 55/46/35 on the SRR fastq, 210/152/127 on ecoli1.fasta, 81/76/73 on salmonella_enterica.fasta for regular/canonical/new respectively. All three report identical positive, negative, invalid, search and extension counts. What changed: - util::compute_minimizer implements the scheme above; canonical_mmer and is_canonical are the two new primitives. Selecting on h(kappa) rather than on min(h(x), h(rc(x))) is equivalent for density and costs one hash per m-mer instead of two. - dna_uint_kmer_t::reverse_complement_mmer reverse-complements an m-mer known to fit in a word with a single crc64, regardless of how wide the kmer type is. alpha_kmer_t gets the identity fallback and a has_reverse_complement flag, so alphabets without a reverse complement (amino acids) reduce to exactly the previous forward minimizer. - minimizer_iterator computes the same thing incrementally; minimizer_iterator_rc is gone, and so is the per-kmer reverse complement in the build. - super-kmer parse: the minimizer value is now a function of the anchored locus, so the break test drops the value comparison. - lookup: one probe, and the two candidate kmer starts (j - pos_in_kmer and j - (k - m - pos_in_kmer)) with no case analysis, instead of up to two probes of two candidates each. lookup_regular/lookup_canonical collapse into lookup. - streaming_query: one anchor instead of two, one comparison for the early exit, no longer templated on the modality. - the index no longer stores the canonical flag: index version 5.1.1 -> 6.0.0, existing indexes need rebuilding. Ties in h(kappa) -- two loci of a window carrying the same canonical m-mer, which happens when x_i == x_j or x_i == rc(x_j) -- cannot be broken by position without breaking mirror-equivariance, which would send x and rc(x) to different buckets: a correctness failure, not a density one. They are broken in the frame of the canonical kmer, the same string for x and rc(x), i.e. leftmost when the kmer is canonical and rightmost otherwise. The rule fires on ~1e-5 of the windows at m=13 and makes the minimizer not strictly forward, which the parser and the lookup already tolerated. test/test_minimizer.cpp checks mirror-equivariance directly and the incremental iterator against the brute-force reference, over k in {5..63} and every m; small m makes 16.5% of the sampled kmers tie, so the tie-breaking rule is heavily exercised (removing it fails the test on 307513 kmers). Verified with --check (lookup, access, membership, negative lookup, navigational kmer and string queries, kmer and string iterators) on salmonella_enterica, salmonella_100, ecoli1, penicillium_chrysogenum and se.ust.k31/k47/k63, for m in {13,15,17,20,21,31}, single- and 4-threaded, plain and --weighted, in both the 64-bit and the 128-bit kmer builds, with assertions enabled -- which checks the incremental minimizer against the reference for every kmer of the input and every streaming lookup against the corresponding single lookup. --- CMakeLists.txt | 5 + README.md | 19 +-- include/builder/dictionary_builder.hpp | 2 - include/builder/util.hpp | 5 +- include/constants.hpp | 6 +- include/dictionary.hpp | 23 ++- include/kmer.hpp | 23 +++ include/minimizer_iterator.hpp | 163 +++++++++---------- include/spectrum_preserving_string_set.hpp | 115 ++++--------- include/streaming_query.hpp | 47 +----- include/util.hpp | 114 +++++++++++-- src/builder/build_sparse_and_skew_index.cpp | 11 +- src/builder/compute_minimizer_tuples.cpp | 23 +-- src/dictionary.cpp | 69 ++------ src/info.cpp | 1 - src/query.cpp | 31 +--- test/test_minimizer.cpp | 172 ++++++++++++++++++++ tools/build.cpp | 4 - tools/sshash.cpp | 1 - 19 files changed, 467 insertions(+), 367 deletions(-) create mode 100644 test/test_minimizer.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index a4c624c..3c4ed9c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -106,6 +106,11 @@ if(SSHASH_BUILD_EXECUTABLES) sshash_static ) + add_executable(test_minimizer test/test_minimizer.cpp) + target_link_libraries(test_minimizer + sshash_static + ) + add_executable(check test/check.cpp) target_link_libraries(check sshash_static diff --git a/README.md b/README.md index 77101b1..d1ab97c 100644 --- a/README.md +++ b/README.md @@ -204,20 +204,14 @@ if your queries are to be read from a (multi-line) FASTA file. ### Example 3 - ./sshash build -i ../data/unitigs_stitched/salmonella_100_k31_ust.fa.gz -k 31 -m 13 --canonical -o salmonella_100.canon.sshash - -This example builds a dictionary from the input file `../data/unitigs_stitched/salmonella_100_k31_ust.fa.gz` (same used in Example 2), with k = 31, m = 13, and with the canonical parsing modality (option `--canonical`). The dictionary is serialized on disk to the file `salmonella_100.canon.sshash`. - -The "canonical" version of the dictionary offers more speed for only a little space increase, especially under low-hit workloads -- when the majority of k-mers are not found in the dictionary. (For all details, refer to the paper.) - -Below a comparison between the dictionary built in Example 2 (not canonical) -and the one just built (Example 3, canonical). - ./sshash query -i salmonella_100.sshash -q ../data/queries/SRR5833294.10K.fastq.gz - ./sshash query -i salmonella_100.canon.sshash -q ../data/queries/SRR5833294.10K.fastq.gz +A k-mer and its reverse complement always share a bucket, so a lookup costs a +single bucket probe and reports, in `kmer_orientation`, which of the two was +found. This is the only indexing modality: there is no `--canonical` flag to +choose, and no space-versus-speed trade-off to make. -Both queries should originate the following report (reported here for reference): +The query above should originate the following report (reported here for reference): ==== query report: num_kmers = 460000 @@ -225,9 +219,6 @@ Both queries should originate the following report (reported here for reference) num_searches = 42/46 (91.3043%) num_extensions = 4/46 (8.69565%) -The canonical dictionary can be twice as fast as the regular dictionary -for low-hit workloads, even on this tiny example, for only +0.3 bits/k-mer. - ### Example 4 ./sshash permute -i ../data/unitigs_stitched/with_weights/ecoli_sakai.ust.k31.fa.gz -k 31 -o ecoli_sakai.permuted.fa diff --git a/include/builder/dictionary_builder.hpp b/include/builder/dictionary_builder.hpp index dea89a5..309343f 100644 --- a/include/builder/dictionary_builder.hpp +++ b/include/builder/dictionary_builder.hpp @@ -20,13 +20,11 @@ struct dictionary_builder // d.m_m = build_config.m; d.m_spss.k = build_config.k; d.m_spss.m = build_config.m; - d.m_canonical = build_config.canonical; d.m_hasher.seed(build_config.seed); build_stats.add("input_filename", filename.c_str()); build_stats.add("k", d.m_k); build_stats.add("m", d.m_m); - build_stats.add("canonical", d.m_canonical ? "true" : "false"); build_stats.add("seed", build_config.seed); build_stats.add("num_threads", build_config.num_threads); diff --git a/include/builder/util.hpp b/include/builder/util.hpp index b390436..2cdbe94 100644 --- a/include/builder/util.hpp +++ b/include/builder/util.hpp @@ -93,8 +93,9 @@ struct bucket_type { iterator end() const { return iterator(m_end); } /* - When a canonical index is built (option `--canonical`), - a minimizer offset can correspond to more than one super-kmer. + A minimizer offset can correspond to more than one super-kmer: the + minimizer is not strictly forward (see `util::compute_minimizer`), so a + locus can be abandoned and later re-selected. A super-kmer is uniquely identified by the couple (minimizer offset, position of minimizer in the first kmer of the super-kmer). These two components, together, give the diff --git a/include/constants.hpp b/include/constants.hpp index b252f6e..215901a 100644 --- a/include/constants.hpp +++ b/include/constants.hpp @@ -20,9 +20,9 @@ constexpr int forward_orientation = 1; constexpr int backward_orientation = -1; namespace current_version_number { -constexpr uint8_t x = 5; -constexpr uint8_t y = 1; -constexpr uint8_t z = 1; +constexpr uint8_t x = 6; +constexpr uint8_t y = 0; +constexpr uint8_t z = 0; } // namespace current_version_number } // namespace sshash::constants diff --git a/include/dictionary.hpp b/include/dictionary.hpp index a30b8c4..4eea0a1 100644 --- a/include/dictionary.hpp +++ b/include/dictionary.hpp @@ -22,8 +22,7 @@ struct dictionary // , m_num_kmers(0) , m_num_strings(0) , m_k(0) - , m_m(0) - , m_canonical(false) {} + , m_m(0) {} /* Build from input file. */ void build(std::string const& input_filename, build_configuration const& build_config); @@ -33,11 +32,16 @@ struct dictionary // uint64_t num_strings() const { return m_num_strings; } uint64_t k() const { return m_k; } uint64_t m() const { return m_m; } - bool canonical() const { return m_canonical; } bool weighted() const { return !m_weights.empty(); } hasher_type const& hasher() const { return m_hasher; } - /* Lookup queries. */ + /* + Lookup queries. A kmer and its reverse complement share a bucket, so a + lookup always costs a single probe and reports, via + `lookup_result::kmer_orientation`, which of the two was found. Pass + `check_reverse_complement = false` to restrict the answer to kmers + occurring in forward orientation. + */ lookup_result lookup(char const* string_kmer, bool check_reverse_complement = true) const; lookup_result lookup(Kmer uint_kmer, bool check_reverse_complement = true) const; @@ -75,7 +79,7 @@ struct dictionary // bool is_member(char const* string_kmer, bool check_reverse_complement = true) const; bool is_member(Kmer uint_kmer, bool check_reverse_complement = true) const; - template + template friend struct streaming_query; streaming_query_report // @@ -144,7 +148,6 @@ struct dictionary // visitor.visit(t.m_num_strings); visitor.visit(t.m_k); visitor.visit(t.m_m); - visitor.visit(t.m_canonical); visitor.visit(t.m_hasher); visitor.visit(t.m_spss); visitor.visit(t.m_ssi); @@ -156,7 +159,6 @@ struct dictionary // uint64_t m_num_strings; uint16_t m_k; uint16_t m_m; - bool m_canonical; hasher_type m_hasher; spectrum_preserving_string_set m_spss; @@ -164,12 +166,7 @@ struct dictionary // weights m_weights; - lookup_result lookup_regular(Kmer uint_kmer) const; - lookup_result lookup_regular(Kmer uint_kmer, minimizer_info mini_info) const; - - lookup_result lookup_canonical(Kmer uint_kmer) const; - lookup_result lookup_canonical(Kmer uint_kmer, Kmer uint_kmer_rc, - minimizer_info mini_info) const; + lookup_result lookup(Kmer uint_kmer, Kmer uint_kmer_rc, minimizer_info mini_info) const; void forward_neighbours(Kmer suffix, neighbourhood& res, bool check_reverse_complement) const; diff --git a/include/kmer.hpp b/include/kmer.hpp index 6fe1520..ea0bb9f 100644 --- a/include/kmer.hpp +++ b/include/kmer.hpp @@ -105,7 +105,14 @@ struct alpha_kmer_t : uint_kmer_t { static char uint64_to_char(uint64_t x) { return alphabet[x]; } // Revcompl only makes sense for DNA, fallback to noop otherwise + static constexpr bool has_reverse_complement = false; [[maybe_unused]] virtual void reverse_complement_inplace(uint64_t) {} + + /* Reverse complement of an m-mer packed in the low m*bits_per_char bits of a + word, for m <= max_m. Same fallback as above: the identity. */ + [[maybe_unused]] static uint64_t reverse_complement_mmer(uint64_t mmer, uint64_t) { + return mmer; + } [[maybe_unused]] static void compute_reverse_complement(char const* input, char* output, uint64_t size) { for (uint64_t i = 0; i != size; ++i) output[i] = input[i]; @@ -156,6 +163,8 @@ struct dna_uint_kmer_t : alpha_kmer_t { return res; } + static constexpr bool has_reverse_complement = true; + [[maybe_unused]] void reverse_complement_inplace(uint64_t k) override { assert(k <= max_k); dna_uint_kmer_t rev(0); @@ -164,6 +173,20 @@ struct dna_uint_kmer_t : alpha_kmer_t { *this = rev; } + /* + Reverse complement of an m-mer packed in the low m*bits_per_char bits of a + single word. This is `reverse_complement_inplace` specialized to a value + that is known to fit in 64 bits, hence a single crc64 regardless of how + wide the kmer type is. Used to canonicalize the m-mers of a kmer, which is + done once per character of the input during construction and once per + character of the query at lookup time. + */ + [[maybe_unused]] static uint64_t reverse_complement_mmer(uint64_t mmer, uint64_t m) { + assert(m <= max_m); + assert(m * bits_per_char < 64); + return crc64(mmer) >> (64 - m * bits_per_char); + } + #ifdef SSHASH_USE_TRADITIONAL_NUCLEOTIDE_ENCODING /* char decimal binary diff --git a/include/minimizer_iterator.hpp b/include/minimizer_iterator.hpp index b53adc8..82c18c9 100644 --- a/include/minimizer_iterator.hpp +++ b/include/minimizer_iterator.hpp @@ -1,11 +1,22 @@ #pragma once #include "kmer.hpp" +#include "util.hpp" namespace sshash { /* - "Re-scan" method. + "Re-scan" method: computes the minimizer of each kmer of a sequence, + sliding one character at a time. + + See `util::compute_minimizer` for the definition of the minimizer; this + iterator computes exactly the same thing incrementally, which the assertion + at the end of `next` checks. + + The only extra state compared to a plain forward minimizer is `m_num_mins`, + the number of loci of the current window attaining the minimum hash: a tie + cannot be broken by position alone without breaking mirror-equivariance, so + when there is one we have to look at the kmer's own orientation. */ template struct minimizer_iterator { @@ -30,6 +41,7 @@ struct minimizer_iterator { void reset() { m_min_pos_in_kmer = 0; m_min_position = m_position - 1; + m_num_mins = 0; } minimizer_info next(Kmer kmer) { @@ -41,130 +53,105 @@ struct minimizer_iterator { m_position += 1; Kmer mmer = kmer; mmer.drop_chars(m_k - m_m); - uint64_t hash = m_hasher.hash(uint64_t(mmer)); + uint64_t value = util::canonical_mmer(uint64_t(mmer), m_m); + uint64_t hash = m_hasher.hash(value); if (hash < m_min_hash) { m_min_hash = hash; - m_min_value = uint64_t(mmer); + m_min_value = value; m_min_position = m_position; m_min_pos_in_kmer = m_k - m_m; + m_num_mins = 1; } else { + /* Only the leftmost minimum can leave the window without a + re-scan, so the count stays valid across the slide. */ assert(m_min_pos_in_kmer > 0); m_min_pos_in_kmer -= 1; + if (hash == m_min_hash) m_num_mins += 1; } } - assert(minimizer_info(m_min_value, m_min_pos_in_kmer) == + minimizer_info mini_info(m_min_value, m_min_position, m_min_pos_in_kmer); + if (m_num_mins > 1) break_tie(kmer, mini_info); + + assert(minimizer_info(mini_info.minimizer, mini_info.pos_in_kmer) == util::compute_minimizer(kmer, m_k, m_m, m_hasher)); - return {m_min_value, m_min_position, m_min_pos_in_kmer}; + return mini_info; } private: uint64_t m_k, m_m; uint64_t m_position, m_min_pos_in_kmer; uint64_t m_min_value, m_min_position, m_min_hash; + uint64_t m_num_mins; hasher_type m_hasher; void rescan(Kmer kmer) { - m_min_hash = constants::invalid_uint64; - m_min_value = constants::invalid_uint64; - m_min_pos_in_kmer = 0; - uint64_t begin = m_position; - for (uint64_t i = 0; i != m_k - m_m + 1; ++i, ++m_position) { + const uint64_t begin = m_position; + + /* first locus, peeled off the loop: see `util::compute_minimizer` */ + { Kmer mmer = kmer; kmer.drop_char(); mmer.take_chars(m_m); - uint64_t hash = m_hasher.hash(uint64_t(mmer)); - if (hash < m_min_hash) { // leftmost - m_min_hash = hash; - m_min_value = uint64_t(mmer); - m_min_pos_in_kmer = i; - } + m_min_value = util::canonical_mmer(uint64_t(mmer), m_m); + m_min_hash = m_hasher.hash(m_min_value); + m_min_pos_in_kmer = 0; + m_num_mins = 1; + ++m_position; } - m_position -= 1; - m_min_position = begin + m_min_pos_in_kmer; - } -}; - -/* - "Re-scan" method. -*/ -template -struct minimizer_iterator_rc { - minimizer_iterator_rc() {} - - minimizer_iterator_rc(uint64_t k, uint64_t m, hasher_type const& hasher, uint64_t position = 0) - : m_k(k) - , m_m(m) - , m_min_value(constants::invalid_uint64) - , m_min_hash(constants::invalid_uint64) - , m_hasher(hasher) // - { - assert(k > 0 and m <= k); - set_position(position); - } - - void set_position(uint64_t position) { - m_position = position; - reset(); - } - void reset() { - m_min_pos_in_kmer = m_k - m_m; - m_min_position = m_position - 1; - } - - minimizer_info next(Kmer kmer) { - if (m_min_pos_in_kmer == m_k - m_m) { - /* min leaves the window: re-scan to compute the new min */ - m_position = m_min_position + 1; - rescan(kmer); - } else { - m_position += 1; + for (uint64_t i = 1; i != m_k - m_m + 1; ++i, ++m_position) { Kmer mmer = kmer; + kmer.drop_char(); mmer.take_chars(m_m); - uint64_t hash = m_hasher.hash(uint64_t(mmer)); - if (hash <= m_min_hash) { + uint64_t value = util::canonical_mmer(uint64_t(mmer), m_m); + uint64_t hash = m_hasher.hash(value); + if (hash < m_min_hash) { // leftmost m_min_hash = hash; - m_min_value = uint64_t(mmer); - m_min_position = m_position; - m_min_pos_in_kmer = 0; - } else { - m_min_pos_in_kmer += 1; - assert(m_min_pos_in_kmer <= m_k - m_m); + m_min_value = value; + m_min_pos_in_kmer = i; + m_num_mins = 1; + } else if (hash == m_min_hash) { + m_num_mins += 1; } } - assert(minimizer_info(m_min_value, m_min_pos_in_kmer) == - util::compute_minimizer(kmer, m_k, m_m, m_hasher)); - - return {m_min_value, m_min_position, m_min_pos_in_kmer}; + m_position -= 1; + m_min_position = begin + m_min_pos_in_kmer; } -private: - uint64_t m_k, m_m; - uint64_t m_position, m_min_pos_in_kmer; - uint64_t m_min_value, m_min_position, m_min_hash; - hasher_type m_hasher; - - void rescan(Kmer kmer) { - m_min_hash = constants::invalid_uint64; - m_min_value = constants::invalid_uint64; - m_min_pos_in_kmer = 0; - uint64_t begin = m_position; - for (int64_t i = m_k - m_m; i >= 0; --i, ++m_position) { - Kmer mmer = kmer; - mmer.drop_chars(i); + /* + Two or more loci of the window attain the minimum hash, and the leftmost + of them is the one currently held. Leftmost is the right answer when the + kmer is canonical; otherwise the canonical frame is rc(kmer), whose + leftmost tied locus is this window's rightmost one. + + This is the only place where the whole kmer, rather than just its m-mers, + has to be reverse-complemented. It runs on ~1e-5 of the windows. + */ + void break_tie(Kmer kmer, minimizer_info& mini_info) const { + assert(m_num_mins > 1); + if (util::is_canonical(kmer, m_k)) return; + + const uint64_t window_begin = m_min_position - m_min_pos_in_kmer; + uint64_t pos_in_kmer = m_min_pos_in_kmer; + uint64_t value = m_min_value; + + Kmer window = kmer; + window.drop_chars(m_min_pos_in_kmer + 1); + for (uint64_t i = m_min_pos_in_kmer + 1; i != m_k - m_m + 1; ++i) { + Kmer mmer = window; mmer.take_chars(m_m); - uint64_t hash = m_hasher.hash(uint64_t(mmer)); - if (hash <= m_min_hash) { // rightmost - m_min_hash = hash; - m_min_value = uint64_t(mmer); - m_min_pos_in_kmer = i; + uint64_t v = util::canonical_mmer(uint64_t(mmer), m_m); + if (m_hasher.hash(v) == m_min_hash) { // rightmost + pos_in_kmer = i; + value = v; } + window.drop_char(); } - m_position -= 1; - m_min_position = begin + (m_k - m_min_pos_in_kmer - m_m); + + mini_info = minimizer_info(value, window_begin + pos_in_kmer, pos_in_kmer); } }; diff --git a/include/spectrum_preserving_string_set.hpp b/include/spectrum_preserving_string_set.hpp index feca86b..5b88aae 100644 --- a/include/spectrum_preserving_string_set.hpp +++ b/include/spectrum_preserving_string_set.hpp @@ -27,9 +27,9 @@ struct spectrum_preserving_string_set // } template - lookup_result lookup_regular(Iterator it, // - const Kmer kmer, // - const minimizer_info mini_info) const // + lookup_result lookup(Iterator it, // + const Kmer kmer, const Kmer kmer_rc, // + const minimizer_info mini_info) const // { const uint64_t size = it.size(); assert(size > 0); @@ -43,13 +43,16 @@ struct spectrum_preserving_string_set // v[i] = strings_offsets.decode(minimizer_offset); } - /* check minimizer first */ + /* Check minimizer first. The minimizer's value is the canonical m-mer at + the anchored locus, so the m-mer stored at that offset is either the + minimizer itself or its reverse complement. */ if (uint64_t read_mmer = uint64_t( util::read_kmer_at(strings, m, Kmer::bits_per_char * v[0].absolute_offset)); - read_mmer != mini_info.minimizer) // + read_mmer != mini_info.minimizer and + util::canonical_mmer(read_mmer, m) != mini_info.minimizer) // { /* - The function `lookup_regular` determines if the minimizer is found at the + This function determines if the minimizer is found at the offset `Kmer::bits_per_char * p.absolute_offset`, not whether the minimizer does not appear at all. In fact, it can happen that the minimizer appear but not at the specified offset, so it would be wrong to set `res.minimizer_found` @@ -66,46 +69,7 @@ struct spectrum_preserving_string_set // lookup_result res; for (uint64_t i = 0; i != size; ++i) { - if (_lookup_regular(res, v[i], kmer, mini_info)) return res; - } - - return lookup_result(); - } - - template - lookup_result lookup_canonical(Iterator it, // - const Kmer kmer, const Kmer kmer_rc, // - const minimizer_info mini_info) const // - { - const uint64_t size = it.size(); - assert(size > 0); - - static thread_local // - std::array - v; - - for (uint64_t i = 0; i != size; ++i, ++it) { - uint64_t minimizer_offset = *it; - v[i] = strings_offsets.decode(minimizer_offset); - } - - /* check minimizer first */ - if (uint64_t read_mmer = uint64_t( - util::read_kmer_at(strings, m, Kmer::bits_per_char * v[0].absolute_offset)); - read_mmer != mini_info.minimizer) // - { - Kmer tmp = mini_info.minimizer; - tmp.reverse_complement_inplace(m); - uint64_t minimizer_rc = uint64_t(tmp); - if (read_mmer != minimizer_rc) { - /* Same note as for the function `lookup_regular`. */ - return lookup_result(it.bucket_type() != bucket_t::HEAVYLOAD ? false : true); - } - } - - lookup_result res; - for (uint64_t i = 0; i != size; ++i) { - if (_lookup_canonical(res, v[i], kmer, kmer_rc, mini_info)) return res; + if (_lookup(res, v[i], kmer, kmer_rc, mini_info)) return res; } return lookup_result(); @@ -210,47 +174,34 @@ struct spectrum_preserving_string_set // visitor.visit(t.strings); } - bool _lookup_regular(lookup_result& res, // - typename Offsets::decoded_offset p, // - const Kmer kmer, // - const minimizer_info mini_info) const // + /* + The minimizer is anchored at locus `pos_in_kmer` of the kmer and, by + mirror-equivariance, at locus `k - m - pos_in_kmer` of its reverse + complement. So a minimizer occurrence at offset j is the anchor of a kmer + starting either at j - pos_in_kmer or at j - (k - m - pos_in_kmer): + always exactly two candidates, with no case analysis. They lie within + 2k-m characters of each other, hence usually on the same cache line. + */ + bool _lookup(lookup_result& res, // + typename Offsets::decoded_offset p, // + const Kmer kmer, // + const Kmer kmer_rc, // + const minimizer_info mini_info) const // { - if (p.absolute_offset < mini_info.pos_in_kmer) return false; - - res.kmer_offset = p.absolute_offset - mini_info.pos_in_kmer; - - if (kmer != util::read_kmer_at(strings, k, Kmer::bits_per_char * res.kmer_offset)) { - return false; - } - - if (res.kmer_offset >= res.string_begin and res.kmer_offset < res.string_end - k + 1) { - res.kmer_id = res.kmer_offset - res.string_id * (k - 1); // absolute kmer id - res.kmer_id_in_string = res.kmer_offset - res.string_begin; // relative kmer id + if constexpr (Kmer::has_reverse_complement) { + if (_lookup_at(res, p, kmer, kmer_rc, mini_info.pos_in_kmer)) return true; + return _lookup_at(res, p, kmer, kmer_rc, k - m - mini_info.pos_in_kmer); } else { - strings_offsets.offset_to_id(res, p, k); + /* no reverse complement: the two candidates would coincide */ + return _lookup_at(res, p, kmer, kmer_rc, mini_info.pos_in_kmer); } - - if (res.kmer_offset < res.string_end - k + 1) return true; - return false; - } - - bool _lookup_canonical(lookup_result& res, // - typename Offsets::decoded_offset p, // - const Kmer kmer, // - const Kmer kmer_rc, // - const minimizer_info mini_info) const // - { - uint64_t pos_in_kmer = mini_info.pos_in_kmer; - if (__lookup_canonical(res, p, kmer, kmer_rc, pos_in_kmer)) return true; - pos_in_kmer = k - m - mini_info.pos_in_kmer; - return __lookup_canonical(res, p, kmer, kmer_rc, pos_in_kmer); } - bool __lookup_canonical(lookup_result& res, // - typename Offsets::decoded_offset p, // - const Kmer kmer, // - const Kmer kmer_rc, // - const uint64_t pos_in_kmer) const // + bool _lookup_at(lookup_result& res, // + typename Offsets::decoded_offset p, // + const Kmer kmer, // + const Kmer kmer_rc, // + const uint64_t pos_in_kmer) const // { if (p.absolute_offset < pos_in_kmer) return false; diff --git a/include/streaming_query.hpp b/include/streaming_query.hpp index e5b8bd7..c039578 100644 --- a/include/streaming_query.hpp +++ b/include/streaming_query.hpp @@ -6,7 +6,7 @@ namespace sshash { -template +template struct streaming_query // { using kmer_t = typename Dict::kmer_type; @@ -22,11 +22,8 @@ struct streaming_query // , m_m(dict->m_m) , m_minimizer_it(dict->m_k, dict->m_m, dict->m_hasher) - , m_minimizer_it_rc(dict->m_k, dict->m_m, dict->m_hasher) , m_curr_mini_info() , m_prev_mini_info() - , m_curr_mini_info_rc() - , m_prev_mini_info_rc() , m_it(dict->m_spss.strings, m_k) , m_remaining_string_bases(0) @@ -36,21 +33,13 @@ struct streaming_query // , m_num_invalid(0) , m_num_negative(0) - { - if (canonical != m_dict->m_canonical) { - std::stringstream ss; - ss << "dict.canonical() = " << (m_dict->canonical() ? "true" : "false") - << " but required " << (canonical ? "true" : "false"); - throw std::runtime_error(ss.str()); - } - } + {} void reset() { m_start = true; m_remaining_string_bases = 0; m_res = lookup_result(); m_minimizer_it.reset(); - m_minimizer_it_rc.reset(); } lookup_result lookup(char const* kmer) // @@ -80,7 +69,6 @@ struct streaming_query // } m_curr_mini_info = m_minimizer_it.next(m_kmer); - m_curr_mini_info_rc = m_minimizer_it_rc.next(m_kmer_rc); /* 3. compute result */ if (m_remaining_string_bases == 0) { @@ -101,7 +89,6 @@ struct streaming_query // /* 4. update state */ m_prev_mini_info = m_curr_mini_info; - m_prev_mini_info_rc = m_curr_mini_info_rc; m_start = false; assert(equal_lookup_result(m_dict->lookup(kmer), m_res)); @@ -127,9 +114,7 @@ struct streaming_query // /* minimizer state */ minimizer_iterator m_minimizer_it; - minimizer_iterator_rc m_minimizer_it_rc; minimizer_info m_curr_mini_info, m_prev_mini_info; - minimizer_info m_curr_mini_info_rc, m_prev_mini_info_rc; /* string state */ kmer_iterator m_it; @@ -147,37 +132,15 @@ struct streaming_query // /* if minimizer does not change and previous minimizer was not found, surely any kmer having the same minimizer cannot be found as well */ - if (m_curr_mini_info.minimizer == m_prev_mini_info.minimizer and // - m_curr_mini_info_rc.minimizer == m_prev_mini_info_rc.minimizer and // - m_res.minimizer_found == false) // + if (m_curr_mini_info.minimizer == m_prev_mini_info.minimizer and // + m_res.minimizer_found == false) // { assert(m_res.kmer_id == constants::invalid_uint64); m_num_negative += 1; return; } - if constexpr (canonical) { - if (m_curr_mini_info.minimizer < m_curr_mini_info_rc.minimizer) { - m_res = m_dict->lookup_canonical(m_kmer, m_kmer_rc, m_curr_mini_info); - } else if (m_curr_mini_info_rc.minimizer < m_curr_mini_info.minimizer) { - m_res = m_dict->lookup_canonical(m_kmer, m_kmer_rc, m_curr_mini_info_rc); - } else { - m_res = m_dict->lookup_canonical(m_kmer, m_kmer_rc, m_curr_mini_info); - if (m_res.kmer_id == constants::invalid_uint64) { - m_res = m_dict->lookup_canonical(m_kmer, m_kmer_rc, m_curr_mini_info_rc); - } - } - } else { - m_res = m_dict->lookup_regular(m_kmer, m_curr_mini_info); - bool minimizer_found = m_res.minimizer_found; - if (m_res.kmer_id == constants::invalid_uint64) { - assert(m_res.kmer_orientation == constants::forward_orientation); - m_res = m_dict->lookup_regular(m_kmer_rc, m_curr_mini_info_rc); - m_res.kmer_orientation = constants::backward_orientation; - bool minimizer_rc_found = m_res.minimizer_found; - m_res.minimizer_found = minimizer_rc_found or minimizer_found; - } - } + m_res = m_dict->lookup(m_kmer, m_kmer_rc, m_curr_mini_info); if (m_res.kmer_id == constants::invalid_uint64) { m_num_negative += 1; diff --git a/include/util.hpp b/include/util.hpp index bf9bebd..06fc7f8 100644 --- a/include/util.hpp +++ b/include/util.hpp @@ -150,7 +150,6 @@ struct build_configuration { , lambda(constants::lambda) - , canonical(false) , weighted(false) , verbose(true) @@ -166,7 +165,6 @@ struct build_configuration { double lambda; // drive PTHash trade-off - bool canonical; bool weighted; bool verbose; @@ -179,7 +177,6 @@ struct build_configuration { << ", num_threads = " << num_threads // << ", ram_limit_in_GiB = " << ram_limit_in_GiB // << ", lambda = " << lambda // - << ", canonical = " << (canonical ? "true" : "false") // << ", weighted = " << (weighted ? "true" : "false") // << ", verbose = " << (verbose ? "true" : "false") // << ", tmp_dirname = '" << tmp_dirname << "'" << std::endl; // @@ -257,7 +254,69 @@ static kmer_t read_kmer_at(bits::bit_vector const& bv, const uint64_t k, const u } /* - This implements the random minimizer. + The canonical m-mer at a locus: the smaller of the m-mer and its reverse + complement, under the numeric order on the packed encoding. For alphabets + that have no reverse complement (e.g. amino acids) this is the identity. +*/ +template +inline uint64_t canonical_mmer(const uint64_t mmer, const uint64_t m) { + if constexpr (!kmer_t::has_reverse_complement) { + (void)m; + return mmer; + } else { + return std::min(mmer, kmer_t::reverse_complement_mmer(mmer, m)); + } +} + +/* + True if `kmer` is the canonical one of the pair (kmer, rc(kmer)). A kmer that + equals its own reverse complement (possible only for even k) is deemed + canonical, so that the answer is always well defined. +*/ +template +inline bool is_canonical(kmer_t kmer, const uint64_t k) { + if constexpr (!kmer_t::has_reverse_complement) { + (void)k; + return true; + } else { + kmer_t kmer_rc = kmer; + kmer_rc.reverse_complement_inplace(k); + return !(kmer_rc < kmer); + } +} + +/* + This implements the random minimizer of a kmer x: the locus i* in [0, k-m] + minimizing h(kappa(i)), where kappa(i) := min(x_i, rc(x_i)) is the canonical + m-mer at locus i. The minimizer's value is kappa(i*). + + Because kappa(i) is invariant under reverse-complementation and the sequence + (kappa(0), ..., kappa(k-m)) of rc(x) is the reversal of that of x, the locus + selected for rc(x) is the mirror image of the one selected for x, and the + value is literally the same m-mer. So x and rc(x) always land in the same + bucket -- a lookup costs a single probe -- while the minimizer remains a + minimizer of the window under a single random order, which is what keeps the + density (and hence the number of super-kmers) at that of the plain forward + minimizer, 2/(k-m+2). + + When the alphabet has no reverse complement, kappa is the identity and this + reduces exactly to the plain forward minimizer. + + Selecting on h(kappa(i)) rather than on min(h(x_i), h(rc(x_i))) -- the two + are interchangeable, since both are strand-symmetric and induce a uniformly + random order on the loci of a window whose 2(k-m+1) m-mers are distinct, so + both have density 2/(k-m+2) -- costs one hash per m-mer instead of two, the + same as the plain forward minimizer. + + Ties -- h(kappa(i)) == h(kappa(j)) for i != j, which happens when x_i == x_j + or x_i == rc(x_j) -- must be broken in a mirror-equivariant way, or x and + rc(x) would be sent to different buckets, which is a correctness failure and + not merely a density one. We break them in the frame of the canonical kmer + min(x, rc(x)), which is literally the same string for x and rc(x): that + amounts to taking the leftmost tied locus when x is canonical and the + rightmost one otherwise. The rule fires on ~1e-5 of the windows for k=31, + m=13, and makes the minimizer not strictly forward, which the parser and the + lookup already tolerate. */ template minimizer_info compute_minimizer(kmer_t kmer, const uint64_t k, const uint64_t m, @@ -265,21 +324,46 @@ minimizer_info compute_minimizer(kmer_t kmer, const uint64_t k, const uint64_t m { assert(m <= kmer_t::max_m); assert(m <= k); - uint64_t min_hash = constants::invalid_uint64; - kmer_t minimizer = kmer_t(-1); - uint64_t pos = 0; - for (uint64_t i = 0; i != k - m + 1; ++i) { - kmer_t mmer = kmer; + + /* The first locus is peeled off the loop so that `min_hash` starts out at a + real hash value: initializing it to invalid_uint64 would make an actual + hash of invalid_uint64 register as a tie rather than as the minimum. */ + kmer_t window = kmer; + kmer_t first = window; + first.take_chars(m); + uint64_t minimizer = canonical_mmer(uint64_t(first), m); + uint64_t min_hash = hasher.hash(minimizer); + uint64_t leftmost = 0; + uint64_t rightmost = 0; + bool tie = false; + window.drop_char(); + + for (uint64_t i = 1; i != k - m + 1; ++i) { + kmer_t mmer = window; mmer.take_chars(m); - uint64_t hash = hasher.hash(uint64_t(mmer)); - if (hash < min_hash) { + uint64_t value = canonical_mmer(uint64_t(mmer), m); + uint64_t hash = hasher.hash(value); + if (hash < min_hash) { // leftmost min_hash = hash; - minimizer = mmer; - pos = i; + minimizer = value; + leftmost = i; + rightmost = i; + tie = false; + } else if (hash == min_hash) { + rightmost = i; + tie = true; } - kmer.drop_char(); + window.drop_char(); + } + + if (tie and !is_canonical(kmer, k)) { + /* rc(kmer) is the canonical frame: mirror its leftmost tied locus */ + kmer.drop_chars(rightmost); + kmer.take_chars(m); + return {canonical_mmer(uint64_t(kmer), m), rightmost}; } - return {uint64_t(minimizer), pos}; + + return {minimizer, leftmost}; } } // namespace util diff --git a/src/builder/build_sparse_and_skew_index.cpp b/src/builder/build_sparse_and_skew_index.cpp index 7ed9886..724b354 100644 --- a/src/builder/build_sparse_and_skew_index.cpp +++ b/src/builder/build_sparse_and_skew_index.cpp @@ -459,12 +459,11 @@ void dictionary_builder::build_sparse_and_skew_index( d.m_spss.strings, k, Kmer::bits_per_char * starting_pos_of_super_kmer); for (uint64_t i = 0; i != mt.num_kmers_in_super_kmer; ++i) { auto kmer = it.get(); - if (build_config.canonical) { /* take the canonical kmer */ - auto kmer_rc = kmer; - kmer_rc.reverse_complement_inplace(k); - kmer = std::min(kmer, kmer_rc); - } - kmers.push_back(kmer); + /* take the canonical kmer: a kmer and its reverse complement + share a bucket, hence a skew-index partition */ + auto kmer_rc = kmer; + kmer_rc.reverse_complement_inplace(k); + kmers.push_back(std::min(kmer, kmer_rc)); positions_in_bucket.push_back(pos_in_bucket); it.next(); } diff --git a/src/builder/compute_minimizer_tuples.cpp b/src/builder/compute_minimizer_tuples.cpp index 94916d6..c9190ce 100644 --- a/src/builder/compute_minimizer_tuples.cpp +++ b/src/builder/compute_minimizer_tuples.cpp @@ -50,7 +50,6 @@ void dictionary_builder::compute_minimizer_tuples() // kmer_iterator kmer_it(strings_builder, k); hasher_type hasher(build_config.seed); minimizer_iterator minimizer_it(k, m, hasher); - minimizer_iterator_rc minimizer_it_rc(k, m, hasher); for (uint64_t i = index_begin; i < index_end; ++i) // { @@ -65,7 +64,6 @@ void dictionary_builder::compute_minimizer_tuples() // kmer_it.at(Kmer::bits_per_char * begin); minimizer_it.set_position(begin); - minimizer_it_rc.set_position(begin); for (uint64_t j = 0; j != sequence_len - k + 1; ++j) { auto uint_kmer = kmer_it.get(); @@ -73,18 +71,6 @@ void dictionary_builder::compute_minimizer_tuples() // assert(mini_info.pos_in_seq < end - m + 1); assert(mini_info.pos_in_kmer < k - m + 1); - if (build_config.canonical) { - auto uint_kmer_rc = uint_kmer; - uint_kmer_rc.reverse_complement_inplace(k); - auto mini_info_rc = minimizer_it_rc.next(uint_kmer_rc); - assert(mini_info_rc.pos_in_seq < end - m + 1); - assert(mini_info_rc.pos_in_kmer < k - m + 1); - if (mini_info_rc.minimizer < mini_info.minimizer) { - mini_info = mini_info_rc; - mini_info.pos_in_kmer = k - m - mini_info.pos_in_kmer; - } - } - mini_info.pos_in_seq = strings_offsets_builder.encode(mini_info.pos_in_seq, begin, i); @@ -92,8 +78,13 @@ void dictionary_builder::compute_minimizer_tuples() // prev_mini_info = mini_info; } - if (mini_info.minimizer != prev_mini_info.minimizer or - mini_info.pos_in_seq != prev_mini_info.pos_in_seq) // + /* + The minimizer's value is the canonical m-mer at the + anchored locus, hence a function of `pos_in_seq` alone: + comparing positions is enough to detect a super-kmer + break, no need to also compare values. + */ + if (mini_info.pos_in_seq != prev_mini_info.pos_in_seq) // { save(prev_mini_info, num_kmers_in_super_kmer); prev_mini_info = mini_info; diff --git a/src/dictionary.cpp b/src/dictionary.cpp index c1b43f7..b4d50ce 100644 --- a/src/dictionary.cpp +++ b/src/dictionary.cpp @@ -4,55 +4,22 @@ namespace sshash { -template -lookup_result dictionary::lookup_regular(const Kmer uint_kmer) const { - auto mini_info = util::compute_minimizer(uint_kmer, m_k, m_m, m_hasher); - return lookup_regular(uint_kmer, mini_info); -} - -template -lookup_result dictionary::lookup_regular(const Kmer uint_kmer, // - const minimizer_info mini_info) const // +/* + The minimizer of a kmer and of its reverse complement is the same m-mer, so a + single bucket probe resolves both orientations at once. This is the only + lookup path: there is no separate "canonical" modality. +*/ +template +lookup_result dictionary::lookup(const Kmer uint_kmer, // + const Kmer uint_kmer_rc, // + const minimizer_info mini_info) const // { assert(minimizer_info(mini_info.minimizer, mini_info.pos_in_kmer) == util::compute_minimizer(uint_kmer, m_k, m_m, m_hasher)); - auto it = m_ssi.lookup(uint_kmer, mini_info); - return m_spss.lookup_regular(it, uint_kmer, mini_info); -} - -template -lookup_result dictionary::lookup_canonical(Kmer uint_kmer) const // -{ - Kmer uint_kmer_rc = uint_kmer; - uint_kmer_rc.reverse_complement_inplace(m_k); - auto mini_info = util::compute_minimizer(uint_kmer, m_k, m_m, m_hasher); - auto mini_info_rc = util::compute_minimizer(uint_kmer_rc, m_k, m_m, m_hasher); - if (mini_info.minimizer < mini_info_rc.minimizer) { - return lookup_canonical(uint_kmer, uint_kmer_rc, mini_info); - } else if (mini_info_rc.minimizer < mini_info.minimizer) { - return lookup_canonical(uint_kmer, uint_kmer_rc, mini_info_rc); - } else { - auto res = lookup_canonical(uint_kmer, uint_kmer_rc, mini_info); - if (res.kmer_id == constants::invalid_uint64) { - res = lookup_canonical(uint_kmer, uint_kmer_rc, mini_info_rc); - } - return res; - } -} - -template -lookup_result dictionary::lookup_canonical(const Kmer uint_kmer, // - const Kmer uint_kmer_rc, // - const minimizer_info mini_info) const // -{ - assert(mini_info.minimizer == - std::min(util::compute_minimizer(uint_kmer, m_k, m_m, m_hasher).minimizer, - util::compute_minimizer(uint_kmer_rc, m_k, m_m, m_hasher).minimizer)); - const Kmer uint_kmer_canon = std::min(uint_kmer, uint_kmer_rc); auto it = m_ssi.lookup(uint_kmer_canon, mini_info); - return m_spss.lookup_canonical(it, uint_kmer, uint_kmer_rc, mini_info); + return m_spss.lookup(it, uint_kmer, uint_kmer_rc, mini_info); } template @@ -65,14 +32,12 @@ template lookup_result dictionary::lookup(Kmer uint_kmer, bool check_reverse_complement) const // { - if (m_canonical) return lookup_canonical(uint_kmer); - auto res = lookup_regular(uint_kmer); - assert(res.kmer_orientation == constants::forward_orientation); - if (check_reverse_complement and res.kmer_id == constants::invalid_uint64) { - Kmer uint_kmer_rc = uint_kmer; - uint_kmer_rc.reverse_complement_inplace(m_k); - res = lookup_regular(uint_kmer_rc); - res.kmer_orientation = constants::backward_orientation; + Kmer uint_kmer_rc = uint_kmer; + uint_kmer_rc.reverse_complement_inplace(m_k); + auto mini_info = util::compute_minimizer(uint_kmer, m_k, m_m, m_hasher); + auto res = lookup(uint_kmer, uint_kmer_rc, mini_info); + if (!check_reverse_complement and res.kmer_orientation == constants::backward_orientation) { // + return lookup_result(); } return res; } @@ -203,7 +168,7 @@ neighbourhood dictionary::string_neighbours( template uint64_t dictionary::num_bits() const { return 8 * (sizeof(m_vnum) + sizeof(m_num_kmers) + sizeof(m_num_strings) + sizeof(m_k) + - sizeof(m_m) + sizeof(m_canonical) + sizeof(m_hasher)) + + sizeof(m_m) + sizeof(m_hasher)) + m_spss.num_bits() + m_ssi.num_bits() + m_weights.num_bits(); } diff --git a/src/info.cpp b/src/info.cpp index cf214fc..a1c5334 100644 --- a/src/info.cpp +++ b/src/info.cpp @@ -55,7 +55,6 @@ void dictionary::print_info() const { std::cout << "k = " << k() << '\n'; std::cout << "num_minimizers = " << m_ssi.codewords.size() << std::endl; std::cout << "m = " << m() << '\n'; - std::cout << "canonical = " << (canonical() ? "true" : "false") << '\n'; std::cout << "weighted = " << (weighted() ? "true" : "false") << '\n'; print_space_breakdown(); } diff --git a/src/query.cpp b/src/query.cpp index e367eee..29f6b1e 100644 --- a/src/query.cpp +++ b/src/query.cpp @@ -121,8 +121,7 @@ dictionary::streaming_query_from_file(std::string const& filename bool multiline) const // { using dictionary_type = dictionary; - using regular_query = streaming_query; - using canonical_query = streaming_query; + using query_type = streaming_query; std::ifstream is(filename.c_str()); if (!is.good()) throw std::runtime_error("error in opening the file '" + filename + "'"); @@ -130,42 +129,22 @@ dictionary::streaming_query_from_file(std::string const& filename if (util::ends_with(filename, ".fa.gz") or util::ends_with(filename, ".fasta.gz")) { zip_istream zis(is); - if (canonical()) { - report = streaming_query_from_fasta_file(this, zis, - multiline); - } else { - report = streaming_query_from_fasta_file(this, zis, - multiline); - } + report = streaming_query_from_fasta_file(this, zis, multiline); } else if (util::ends_with(filename, ".fq.gz") or util::ends_with(filename, ".fastq.gz")) { if (multiline) { std::cout << "==> Warning: option 'multiline' is only valid for FASTA files, not FASTQ." << std::endl; } zip_istream zis(is); - if (canonical()) { - report = streaming_query_from_fastq_file(this, zis); - } else { - report = streaming_query_from_fastq_file(this, zis); - } + report = streaming_query_from_fastq_file(this, zis); } else if (util::ends_with(filename, ".fa") or util::ends_with(filename, ".fasta")) { - if (canonical()) { - report = streaming_query_from_fasta_file(this, is, - multiline); - } else { - report = streaming_query_from_fasta_file(this, is, - multiline); - } + report = streaming_query_from_fasta_file(this, is, multiline); } else if (util::ends_with(filename, ".fq") or util::ends_with(filename, ".fastq")) { if (multiline) { std::cout << "==> Warning: option 'multiline' is only valid for FASTA files, not FASTQ." << std::endl; } - if (canonical()) { - report = streaming_query_from_fastq_file(this, is); - } else { - report = streaming_query_from_fastq_file(this, is); - } + report = streaming_query_from_fastq_file(this, is); } else { std::cerr << "unsupported query file format" << std::endl; } diff --git a/test/test_minimizer.cpp b/test/test_minimizer.cpp new file mode 100644 index 0000000..81e2b49 --- /dev/null +++ b/test/test_minimizer.cpp @@ -0,0 +1,172 @@ +#include +#include +#include +#include + +#include "include/util.hpp" +#include "include/kmer.hpp" +#include "include/minimizer_iterator.hpp" + +using namespace sshash; + +static uint64_t num_failed = 0; +static uint64_t num_checked = 0; +static uint64_t num_ties = 0; + +static void fail(std::string const& what) { + if (num_failed < 10) std::cerr << "FAILED: " << what << std::endl; + ++num_failed; +} + +/* + The property the whole design rests on: the minimizer must be + mirror-equivariant. The m-mer selected for a kmer and for its reverse + complement must be literally the same string -- otherwise the two would land + in different buckets and a lookup by one orientation could not find the other + -- and the selected locus must be the mirror image. +*/ +template +static void check_equivariance(kmer_t kmer, uint64_t k, uint64_t m, hasher_type const& hasher) { + kmer_t kmer_rc = kmer; + kmer_rc.reverse_complement_inplace(k); + + auto a = util::compute_minimizer(kmer, k, m, hasher); + auto b = util::compute_minimizer(kmer_rc, k, m, hasher); + + ++num_checked; + + if (a.minimizer != b.minimizer) { + fail("minimizer of kmer '" + util::uint_kmer_to_string(kmer, k) + "' (m=" + + std::to_string(m) + ") is " + util::uint_minimizer_to_string(a.minimizer, m) + + " but that of its reverse complement is " + + util::uint_minimizer_to_string(b.minimizer, m)); + return; + } + if (a.pos_in_kmer + b.pos_in_kmer != k - m) { + fail("locus of kmer '" + util::uint_kmer_to_string(kmer, k) + + "' is " + std::to_string(a.pos_in_kmer) + " but that of its reverse complement is " + + std::to_string(b.pos_in_kmer) + " (should mirror to " + + std::to_string(k - m - a.pos_in_kmer) + ")"); + return; + } + + /* the reported value must really be the canonical m-mer at the reported locus */ + kmer_t anchored = kmer; + anchored.drop_chars(a.pos_in_kmer); + anchored.take_chars(m); + if (a.minimizer != util::canonical_mmer(uint64_t(anchored), m)) { + fail("minimizer value does not match the m-mer at the selected locus"); + } +} + +/* Count the loci attaining the minimum hash, to report how often ties fire. */ +template +static bool has_tie(kmer_t kmer, uint64_t k, uint64_t m, hasher_type const& hasher) { + uint64_t min_hash = constants::invalid_uint64; + uint64_t count = 0; + for (uint64_t i = 0; i != k - m + 1; ++i) { + kmer_t mmer = kmer; + mmer.take_chars(m); + uint64_t hash = hasher.hash(util::canonical_mmer(uint64_t(mmer), m)); + if (i == 0 or hash < min_hash) { + min_hash = hash; + count = 1; + } else if (hash == min_hash) { + ++count; + } + kmer.drop_char(); + } + return count > 1; +} + +template +static void test_kmers(uint64_t k, uint64_t m, uint64_t num_kmers, std::mt19937_64& gen) { + hasher_type hasher(1234); + std::string s(k, 0); + for (uint64_t n = 0; n != num_kmers; ++n) { + for (uint64_t i = 0; i != k; ++i) s[i] = "ACGT"[gen() & 3]; + kmer_t kmer = util::string_to_uint_kmer(s.data(), k); + if (has_tie(kmer, k, m, hasher)) ++num_ties; + check_equivariance(kmer, k, m, hasher); + } +} + +/* + The incremental "re-scan" iterator must agree with the brute-force reference + on every window of a sequence, including the windows where a tie fires and + the anchor consequently moves backwards. +*/ +template +static void test_iterator(uint64_t k, uint64_t m, uint64_t length, std::mt19937_64& gen) { + hasher_type hasher(1234); + std::string s(length, 0); + for (uint64_t i = 0; i != length; ++i) s[i] = "ACGT"[gen() & 3]; + + minimizer_iterator it(k, m, hasher); + it.set_position(0); + + for (uint64_t i = 0; i + k <= length; ++i) { + kmer_t kmer = util::string_to_uint_kmer(s.data() + i, k); + auto got = it.next(kmer); + auto expected = util::compute_minimizer(kmer, k, m, hasher); + ++num_checked; + if (got.minimizer != expected.minimizer or got.pos_in_kmer != expected.pos_in_kmer) { + fail("iterator disagrees with the reference at position " + std::to_string(i) + + " (k=" + std::to_string(k) + ", m=" + std::to_string(m) + ")"); + return; + } + /* pos_in_seq must be the absolute position of the selected locus */ + if (got.pos_in_seq != i + got.pos_in_kmer) { + fail("iterator reports pos_in_seq " + std::to_string(got.pos_in_seq) + " but expected " + + std::to_string(i + got.pos_in_kmer)); + return; + } + } +} + +int main() { + using kmer_t = dna_uint_kmer_t; + using wide_kmer_t = dna_uint_kmer_t<__uint128_t>; + + std::mt19937_64 gen(42); + + std::cout << "checking mirror-equivariance of the minimizer..." << std::endl; + for (uint64_t k : {5, 15, 21, 31}) { + for (uint64_t m = 1; m <= std::min(k, kmer_t::max_m); ++m) { + test_kmers(k, m, 20000, gen); + } + } + for (uint64_t k : {31, 47, 63}) { + for (uint64_t m : {1, 2, 3, 7, 13, 21, 31}) { + if (m > k) continue; + test_kmers(k, m, 20000, gen); + } + } + + /* + Small m makes the m-mer universe tiny, so a window very often contains + two loci carrying the same canonical m-mer: this is what stresses the + tie-breaking rule, which on real data with m=13 fires on ~1e-5 of the + windows only. + */ + std::cout << " " << num_ties << "/" << num_checked << " of the kmers checked had a tie" + << std::endl; + + std::cout << "checking the incremental iterator against the reference..." << std::endl; + for (uint64_t k : {5, 15, 31}) { + for (uint64_t m = 1; m <= std::min(k, kmer_t::max_m); ++m) { + test_iterator(k, m, 20000, gen); + } + } + for (uint64_t k : {47, 63}) { + for (uint64_t m : {2, 3, 13, 31}) { test_iterator(k, m, 20000, gen); } + } + + std::cout << "checked " << num_checked << " kmers" << std::endl; + if (num_failed != 0) { + std::cerr << num_failed << " CHECKS FAILED" << std::endl; + return 1; + } + std::cout << "EVERYTHING OK!" << std::endl; + return 0; +} diff --git a/tools/build.cpp b/tools/build.cpp index 6630386..977c4dc 100644 --- a/tools/build.cpp +++ b/tools/build.cpp @@ -39,9 +39,6 @@ int build(int argc, char** argv) { "RAM limit in GiB. Default value is " + std::to_string(constants::default_ram_limit_in_GiB) + " GiB.", "-g", false); - parser.add("canonical", - "This option results in a trade-off between index space and lookup time.", - "--canonical", false, true); parser.add("weighted", "Also store the weights in compressed format.", "--weighted", false, true); parser.add("check", "Check correctness after construction.", "--check", false, true); @@ -59,7 +56,6 @@ int build(int argc, char** argv) { if (parser.parsed("seed")) build_config.seed = parser.get("seed"); if (parser.parsed("lambda")) build_config.lambda = parser.get("lambda"); - build_config.canonical = parser.get("canonical"); build_config.weighted = parser.get("weighted"); build_config.verbose = parser.get("verbose"); if (parser.parsed("tmp_dirname")) { diff --git a/tools/sshash.cpp b/tools/sshash.cpp index 8138a93..5d14509 100644 --- a/tools/sshash.cpp +++ b/tools/sshash.cpp @@ -53,7 +53,6 @@ int bench(int argc, char** argv) { perf_stats.add("index_filename", index_filename.c_str()); perf_stats.add("k", dict.k()); perf_stats.add("m", dict.m()); - perf_stats.add("canonical", dict.canonical() ? "true" : "false"); perf_test_lookup_access(dict, perf_stats); if (dict.weighted()) perf_test_lookup_weight(dict, perf_stats); From 236b8993b8cc098cc075776627546875e024463f Mon Sep 17 00:00:00 2001 From: jermp Date: Sat, 25 Jul 2026 22:46:41 +0200 Subject: [PATCH 02/14] minor --- include/util.hpp | 84 +++++++++++++++++++++++----------------------- src/dictionary.cpp | 5 --- 2 files changed, 42 insertions(+), 47 deletions(-) diff --git a/include/util.hpp b/include/util.hpp index 06fc7f8..8b83684 100644 --- a/include/util.hpp +++ b/include/util.hpp @@ -74,10 +74,10 @@ inline std::ostream& operator<<(std::ostream& os, lookup_result const& res) { return os; } -template +template struct neighbourhood { - std::array forward; - std::array backward; + std::array forward; + std::array backward; }; struct minimizer_info { @@ -201,55 +201,55 @@ static inline uint64_t get_seed_for_hash_function(build_configuration const& bui return std::equal(pattern.begin(), pattern.end(), str.end() - pattern.size()); } -template -[[maybe_unused]] static kmer_t string_to_uint_kmer(char const* str, uint64_t k) { - assert(k <= kmer_t::max_k); - kmer_t x = 0; - for (uint64_t i = 0; i != k; ++i) x.set(i, kmer_t::char_to_uint(str[i])); +template +[[maybe_unused]] static Kmer string_to_uint_kmer(char const* str, uint64_t k) { + assert(k <= Kmer::max_k); + Kmer x = 0; + for (uint64_t i = 0; i != k; ++i) x.set(i, Kmer::char_to_uint(str[i])); return x; } -template -static void uint_kmer_to_string(kmer_t x, char* str, uint64_t k) { - assert(k <= kmer_t::max_k); - for (uint64_t i = 0; i != k; ++i) str[i] = kmer_t::uint64_to_char(x.pop_char()); +template +static void uint_kmer_to_string(Kmer x, char* str, uint64_t k) { + assert(k <= Kmer::max_k); + for (uint64_t i = 0; i != k; ++i) str[i] = Kmer::uint64_to_char(x.pop_char()); } -template -[[maybe_unused]] static std::string uint_kmer_to_string(kmer_t x, uint64_t k) { - assert(k <= kmer_t::max_k); +template +[[maybe_unused]] static std::string uint_kmer_to_string(Kmer x, uint64_t k) { + assert(k <= Kmer::max_k); std::string str; str.resize(k); uint_kmer_to_string(x, str.data(), k); return str; } -template +template [[maybe_unused]] static std::string uint_minimizer_to_string(uint64_t minimizer, uint64_t m) { - assert(m <= kmer_t::max_m); + assert(m <= Kmer::max_m); std::string str; str.resize(m); - kmer_t x = minimizer; + Kmer x = minimizer; uint_kmer_to_string(x, str.data(), m); return str; } -template +template [[maybe_unused]] static bool is_valid(char const* str, uint64_t size) { for (uint64_t i = 0; i != size; ++i) { - if (!kmer_t::is_valid(str[i])) return false; + if (!Kmer::is_valid(str[i])) return false; } return true; } -template -static kmer_t read_kmer_at(bits::bit_vector const& bv, const uint64_t k, const uint64_t pos) { - static_assert(kmer_t::uint_kmer_bits % 64 == 0); - kmer_t kmer = 0; - for (int i = kmer_t::uint_kmer_bits - 64; i >= 0; i -= 64) { +template +static Kmer read_kmer_at(bits::bit_vector const& bv, const uint64_t k, const uint64_t pos) { + static_assert(Kmer::uint_kmer_bits % 64 == 0); + Kmer kmer = 0; + for (int i = Kmer::uint_kmer_bits - 64; i >= 0; i -= 64) { if (pos + i < bv.num_bits()) kmer.append64(bv.get_word64(pos + i)); } - kmer.take(kmer_t::bits_per_char * k); + kmer.take(Kmer::bits_per_char * k); return kmer; } @@ -258,13 +258,13 @@ static kmer_t read_kmer_at(bits::bit_vector const& bv, const uint64_t k, const u complement, under the numeric order on the packed encoding. For alphabets that have no reverse complement (e.g. amino acids) this is the identity. */ -template +template inline uint64_t canonical_mmer(const uint64_t mmer, const uint64_t m) { - if constexpr (!kmer_t::has_reverse_complement) { + if constexpr (!Kmer::has_reverse_complement) { (void)m; return mmer; } else { - return std::min(mmer, kmer_t::reverse_complement_mmer(mmer, m)); + return std::min(mmer, Kmer::reverse_complement_mmer(mmer, m)); } } @@ -273,13 +273,13 @@ inline uint64_t canonical_mmer(const uint64_t mmer, const uint64_t m) { equals its own reverse complement (possible only for even k) is deemed canonical, so that the answer is always well defined. */ -template -inline bool is_canonical(kmer_t kmer, const uint64_t k) { - if constexpr (!kmer_t::has_reverse_complement) { +template +inline bool is_canonical(Kmer kmer, const uint64_t k) { + if constexpr (!Kmer::has_reverse_complement) { (void)k; return true; } else { - kmer_t kmer_rc = kmer; + Kmer kmer_rc = kmer; kmer_rc.reverse_complement_inplace(k); return !(kmer_rc < kmer); } @@ -318,20 +318,20 @@ inline bool is_canonical(kmer_t kmer, const uint64_t k) { m=13, and makes the minimizer not strictly forward, which the parser and the lookup already tolerate. */ -template -minimizer_info compute_minimizer(kmer_t kmer, const uint64_t k, const uint64_t m, +template +minimizer_info compute_minimizer(Kmer kmer, const uint64_t k, const uint64_t m, hasher_type const& hasher) // { - assert(m <= kmer_t::max_m); + assert(m <= Kmer::max_m); assert(m <= k); /* The first locus is peeled off the loop so that `min_hash` starts out at a real hash value: initializing it to invalid_uint64 would make an actual hash of invalid_uint64 register as a tie rather than as the minimum. */ - kmer_t window = kmer; - kmer_t first = window; + Kmer window = kmer; + Kmer first = window; first.take_chars(m); - uint64_t minimizer = canonical_mmer(uint64_t(first), m); + uint64_t minimizer = canonical_mmer(uint64_t(first), m); uint64_t min_hash = hasher.hash(minimizer); uint64_t leftmost = 0; uint64_t rightmost = 0; @@ -339,9 +339,9 @@ minimizer_info compute_minimizer(kmer_t kmer, const uint64_t k, const uint64_t m window.drop_char(); for (uint64_t i = 1; i != k - m + 1; ++i) { - kmer_t mmer = window; + Kmer mmer = window; mmer.take_chars(m); - uint64_t value = canonical_mmer(uint64_t(mmer), m); + uint64_t value = canonical_mmer(uint64_t(mmer), m); uint64_t hash = hasher.hash(value); if (hash < min_hash) { // leftmost min_hash = hash; @@ -360,7 +360,7 @@ minimizer_info compute_minimizer(kmer_t kmer, const uint64_t k, const uint64_t m /* rc(kmer) is the canonical frame: mirror its leftmost tied locus */ kmer.drop_chars(rightmost); kmer.take_chars(m); - return {canonical_mmer(uint64_t(kmer), m), rightmost}; + return {canonical_mmer(uint64_t(kmer), m), rightmost}; } return {minimizer, leftmost}; diff --git a/src/dictionary.cpp b/src/dictionary.cpp index b4d50ce..83ce986 100644 --- a/src/dictionary.cpp +++ b/src/dictionary.cpp @@ -4,11 +4,6 @@ namespace sshash { -/* - The minimizer of a kmer and of its reverse complement is the same m-mer, so a - single bucket probe resolves both orientations at once. This is the only - lookup path: there is no separate "canonical" modality. -*/ template lookup_result dictionary::lookup(const Kmer uint_kmer, // const Kmer uint_kmer_rc, // From ba4172d9458b6f8f2569dbf4b21a4a017db990eb Mon Sep 17 00:00:00 2001 From: jermp Date: Wed, 29 Jul 2026 14:09:06 +0200 Subject: [PATCH 03/14] minor --- include/sparse_and_skew_index.hpp | 7 +++++-- include/util.hpp | 12 ++---------- src/dictionary.cpp | 23 ++++++++++------------- 3 files changed, 17 insertions(+), 25 deletions(-) diff --git a/include/sparse_and_skew_index.hpp b/include/sparse_and_skew_index.hpp index 4509558..96f73fc 100644 --- a/include/sparse_and_skew_index.hpp +++ b/include/sparse_and_skew_index.hpp @@ -109,7 +109,9 @@ struct sparse_and_skew_index // bits::compact_vector::iterator m_it; }; - bucket_iterator lookup(const Kmer uint_kmer, const minimizer_info mini_info) const // + bucket_iterator lookup(const Kmer uint_kmer, // + const Kmer uint_kmer_rc, // + const minimizer_info mini_info) const // { uint64_t code = codewords.lookup(mini_info.minimizer); @@ -132,7 +134,8 @@ struct sparse_and_skew_index // } assert(status == bucket_t::HEAVYLOAD); // minimizer is part of the skew index - uint64_t offset = ski.lookup(uint_kmer, code); + const Kmer uint_kmer_canon = std::min(uint_kmer, uint_kmer_rc); + uint64_t offset = ski.lookup(uint_kmer_canon, code); return {this, offset, 1, bucket_t::HEAVYLOAD}; } diff --git a/include/util.hpp b/include/util.hpp index 8b83684..af2415a 100644 --- a/include/util.hpp +++ b/include/util.hpp @@ -254,7 +254,7 @@ static Kmer read_kmer_at(bits::bit_vector const& bv, const uint64_t k, const uin } /* - The canonical m-mer at a locus: the smaller of the m-mer and its reverse + The canonical m-mer: the smaller of the m-mer and its reverse complement, under the numeric order on the packed encoding. For alphabets that have no reverse complement (e.g. amino acids) this is the identity. */ @@ -302,21 +302,13 @@ inline bool is_canonical(Kmer kmer, const uint64_t k) { When the alphabet has no reverse complement, kappa is the identity and this reduces exactly to the plain forward minimizer. - Selecting on h(kappa(i)) rather than on min(h(x_i), h(rc(x_i))) -- the two - are interchangeable, since both are strand-symmetric and induce a uniformly - random order on the loci of a window whose 2(k-m+1) m-mers are distinct, so - both have density 2/(k-m+2) -- costs one hash per m-mer instead of two, the - same as the plain forward minimizer. - Ties -- h(kappa(i)) == h(kappa(j)) for i != j, which happens when x_i == x_j or x_i == rc(x_j) -- must be broken in a mirror-equivariant way, or x and rc(x) would be sent to different buckets, which is a correctness failure and not merely a density one. We break them in the frame of the canonical kmer min(x, rc(x)), which is literally the same string for x and rc(x): that amounts to taking the leftmost tied locus when x is canonical and the - rightmost one otherwise. The rule fires on ~1e-5 of the windows for k=31, - m=13, and makes the minimizer not strictly forward, which the parser and the - lookup already tolerate. + rightmost one otherwise. */ template minimizer_info compute_minimizer(Kmer kmer, const uint64_t k, const uint64_t m, diff --git a/src/dictionary.cpp b/src/dictionary.cpp index 83ce986..748d0c4 100644 --- a/src/dictionary.cpp +++ b/src/dictionary.cpp @@ -4,19 +4,6 @@ namespace sshash { -template -lookup_result dictionary::lookup(const Kmer uint_kmer, // - const Kmer uint_kmer_rc, // - const minimizer_info mini_info) const // -{ - assert(minimizer_info(mini_info.minimizer, mini_info.pos_in_kmer) == - util::compute_minimizer(uint_kmer, m_k, m_m, m_hasher)); - - const Kmer uint_kmer_canon = std::min(uint_kmer, uint_kmer_rc); - auto it = m_ssi.lookup(uint_kmer_canon, mini_info); - return m_spss.lookup(it, uint_kmer, uint_kmer_rc, mini_info); -} - template lookup_result dictionary::lookup(char const* string_kmer, bool check_reverse_complement) const { @@ -36,6 +23,16 @@ lookup_result dictionary::lookup(Kmer uint_kmer, } return res; } +template +lookup_result dictionary::lookup(const Kmer uint_kmer, // + const Kmer uint_kmer_rc, // + const minimizer_info mini_info) const // +{ + assert(minimizer_info(mini_info.minimizer, mini_info.pos_in_kmer) == + util::compute_minimizer(uint_kmer, m_k, m_m, m_hasher)); + auto it = m_ssi.lookup(uint_kmer, uint_kmer_rc, mini_info); + return m_spss.lookup(it, uint_kmer, uint_kmer_rc, mini_info); +} template bool dictionary::is_member(char const* string_kmer, From c8481bccec6825b18d8fecf3d5a8d085bba43da0 Mon Sep 17 00:00:00 2001 From: Giulio Ermanno Pibiri Date: Mon, 10 Aug 2026 10:37:47 +0000 Subject: [PATCH 04/14] reverse-complement the kmer once, not every m-mer Canonicalizing the m-mers of a kmer called reverse_complement_mmer once per locus: k-m+1 crc64's per kmer, each reverse-complementing m symbols. Since rc(x_i) = rc(x)_{k-m-i}, the reverse complements of all the loci are windows of the single reverse complement of the whole kmer, so one suffices and each rc(x_i) is a shift out of it. New primitive util::canonical_mmer_at(window, kmer_rc, k, m, i), used by both the stateless compute_minimizer and the sliding iterator. When the kmer fits in a word -- always for the 64-bit kmer type, and for k <= 32 with the wider one -- the shift is a plain 64-bit one; otherwise it is a shift of the wide type. That distinction matters: extracting from a 128-bit value costs a few instructions more than a crc64 on an already-truncated 64-bit m-mer, so without the narrow path the wide kmer type at small k *regressed* by up to 23%, and with it the same case improves by 7-22%. Each of the two callers gets the reverse complement for free: - dictionary::lookup already computes it, one line above, to tell the two orientations apart; it is now passed in. compute_minimizer keeps a four-argument overload that computes it for callers that have none. - minimizer_iterator now carries it and slides it one character at a time (the new character's complement in at the front, the oldest out at the back) exactly as streaming_query already does for its own copy, so it costs a handful of instructions per kmer rather than a crc64. Like the window minimum, this assumes consecutive kmers, which set_position and reset restart; an assertion checks the slid value against a recomputed one on every kmer. It also removes the last per-locus reverse complement in break_tie and the separate is_canonical call, both of which now read the maintained value. Measured, salmonella_100 k=31, interleaved runs (min / median of 7): step 2, the minimizer iterator m=13 -10.9% / -12.8% m=20 -8.0% / -11.4% total build time m=13 -9.1% / -4.3% m=20 -1.8% / -0.2% streaming query, low hit -2.5% / -2.1% streaming query, 100% hit -3.7% / -5.4% and for compute_minimizer in isolation, per kmer: 64-bit kmer, k=31, m=11..24 -13% .. -30% 128-bit kmer, k=63, m=11..24 -14% .. -17% 128-bit kmer, k=31, m=11..24 -7% .. -22% The random-lookup path shows no measurable end-to-end change: the minimizer is ~85ns of a ~500ns lookup and the run-to-run spread on this machine is larger than the ~4% expected, so only the isolated figures above are meaningful there. The build and streaming paths, which go through the iterator, do move consistently. Index output is unchanged: same minimizer for every kmer, hence the same 482031 super-kmers and 4.94628 bits/kmer on salmonella_enterica k=31 m=13. Verified with --check on salmonella_enterica (m in {11,13,20,31}), ecoli1 (4 threads), a --weighted build, and se.ust.k47/k63, with assertions enabled -- which check the identity rc(x_i) = rc(x)_{k-m-i} at every locus of every kmer, the slid reverse complement against a recomputed one, and the iterator against the reference minimizer. --- include/kmer.hpp | 16 +++++++ include/minimizer_iterator.hpp | 71 +++++++++++++++++++-------- include/util.hpp | 88 ++++++++++++++++++++++++++++++---- src/dictionary.cpp | 6 ++- 4 files changed, 150 insertions(+), 31 deletions(-) diff --git a/include/kmer.hpp b/include/kmer.hpp index ea0bb9f..4db3d00 100644 --- a/include/kmer.hpp +++ b/include/kmer.hpp @@ -113,6 +113,9 @@ struct alpha_kmer_t : uint_kmer_t { [[maybe_unused]] static uint64_t reverse_complement_mmer(uint64_t mmer, uint64_t) { return mmer; } + + /* Complement of a single packed character. Same fallback: the identity. */ + [[maybe_unused]] static uint64_t complement_char(uint64_t c) { return c; } [[maybe_unused]] static void compute_reverse_complement(char const* input, char* output, uint64_t size) { for (uint64_t i = 0; i != size; ++i) output[i] = input[i]; @@ -187,6 +190,19 @@ struct dna_uint_kmer_t : alpha_kmer_t { return crc64(mmer) >> (64 - m * bits_per_char); } + /* + Complement of a single packed character: the per-character half of what + crc64 does to a whole word. Lets a reverse complement be slid one + character at a time instead of recomputed. + */ + [[maybe_unused]] static uint64_t complement_char(uint64_t c) { +#ifdef SSHASH_USE_TRADITIONAL_NUCLEOTIDE_ENCODING + return c ^ 3; // A<->T is 00<->11, C<->G is 01<->10 +#else + return c ^ 2; // A<->T is 00<->10, C<->G is 01<->11 +#endif + } + #ifdef SSHASH_USE_TRADITIONAL_NUCLEOTIDE_ENCODING /* char decimal binary diff --git a/include/minimizer_iterator.hpp b/include/minimizer_iterator.hpp index 82c18c9..37a81fa 100644 --- a/include/minimizer_iterator.hpp +++ b/include/minimizer_iterator.hpp @@ -13,10 +13,17 @@ namespace sshash { iterator computes exactly the same thing incrementally, which the assertion at the end of `next` checks. - The only extra state compared to a plain forward minimizer is `m_num_mins`, - the number of loci of the current window attaining the minimum hash: a tie + The extra state compared to a plain forward minimizer is `m_num_mins`, the + number of loci of the current window attaining the minimum hash -- a tie cannot be broken by position alone without breaking mirror-equivariance, so - when there is one we have to look at the kmer's own orientation. + when there is one we have to look at the kmer's own orientation -- and + `m_kmer_rc`, the reverse complement of the current kmer. + + `m_kmer_rc` serves every locus at once, via rc(x_i) = rc(x)_{k-m-i}, so no + m-mer is ever reverse-complemented on its own; and it is itself slid one + character at a time rather than recomputed. Like the window minimum, it + assumes the kmers handed to `next` are consecutive, which `set_position` and + `reset` restart. */ template struct minimizer_iterator { @@ -42,18 +49,21 @@ struct minimizer_iterator { m_min_pos_in_kmer = 0; m_min_position = m_position - 1; m_num_mins = 0; + m_fresh = true; } minimizer_info next(Kmer kmer) { + slide_reverse_complement(kmer); + if (m_min_pos_in_kmer == 0) { /* min leaves the window: re-scan to compute the new min */ m_position = m_min_position + 1; rescan(kmer); } else { m_position += 1; - Kmer mmer = kmer; - mmer.drop_chars(m_k - m_m); - uint64_t value = util::canonical_mmer(uint64_t(mmer), m_m); + Kmer window = kmer; + window.drop_chars(m_k - m_m); + uint64_t value = util::canonical_mmer_at(window, m_kmer_rc, m_k, m_m, m_k - m_m); uint64_t hash = m_hasher.hash(value); if (hash < m_min_hash) { m_min_hash = hash; @@ -74,7 +84,7 @@ struct minimizer_iterator { if (m_num_mins > 1) break_tie(kmer, mini_info); assert(minimizer_info(mini_info.minimizer, mini_info.pos_in_kmer) == - util::compute_minimizer(kmer, m_k, m_m, m_hasher)); + util::compute_minimizer(kmer, m_kmer_rc, m_k, m_m, m_hasher)); return mini_info; } @@ -84,17 +94,45 @@ struct minimizer_iterator { uint64_t m_position, m_min_pos_in_kmer; uint64_t m_min_value, m_min_position, m_min_hash; uint64_t m_num_mins; + Kmer m_kmer_rc; + bool m_fresh; hasher_type m_hasher; + /* + rc(x) for the kmer just handed in. After a restart it is computed + outright; otherwise the kmer has slid by one character, so its reverse + complement has too: the new character's complement enters at the front + and the oldest one falls off the back. + */ + void slide_reverse_complement(Kmer const& kmer) { + if constexpr (!Kmer::has_reverse_complement) { + (void)kmer; + return; + } else { + if (m_fresh) { + m_kmer_rc = kmer; + m_kmer_rc.reverse_complement_inplace(m_k); + m_fresh = false; + } else { + m_kmer_rc.pad_char(); + m_kmer_rc.set(0, Kmer::complement_char(kmer.at(m_k - 1))); + m_kmer_rc.take(m_k * Kmer::bits_per_char); + } + assert([&] { + Kmer expected = kmer; + expected.reverse_complement_inplace(m_k); + return expected == m_kmer_rc; + }()); + } + } + void rescan(Kmer kmer) { const uint64_t begin = m_position; /* first locus, peeled off the loop: see `util::compute_minimizer` */ { - Kmer mmer = kmer; + m_min_value = util::canonical_mmer_at(kmer, m_kmer_rc, m_k, m_m, 0); kmer.drop_char(); - mmer.take_chars(m_m); - m_min_value = util::canonical_mmer(uint64_t(mmer), m_m); m_min_hash = m_hasher.hash(m_min_value); m_min_pos_in_kmer = 0; m_num_mins = 1; @@ -102,10 +140,8 @@ struct minimizer_iterator { } for (uint64_t i = 1; i != m_k - m_m + 1; ++i, ++m_position) { - Kmer mmer = kmer; + uint64_t value = util::canonical_mmer_at(kmer, m_kmer_rc, m_k, m_m, i); kmer.drop_char(); - mmer.take_chars(m_m); - uint64_t value = util::canonical_mmer(uint64_t(mmer), m_m); uint64_t hash = m_hasher.hash(value); if (hash < m_min_hash) { // leftmost m_min_hash = hash; @@ -127,12 +163,11 @@ struct minimizer_iterator { kmer is canonical; otherwise the canonical frame is rc(kmer), whose leftmost tied locus is this window's rightmost one. - This is the only place where the whole kmer, rather than just its m-mers, - has to be reverse-complemented. It runs on ~1e-5 of the windows. + It runs on ~1e-5 of the windows. */ void break_tie(Kmer kmer, minimizer_info& mini_info) const { assert(m_num_mins > 1); - if (util::is_canonical(kmer, m_k)) return; + if (!(m_kmer_rc < kmer)) return; // the kmer is already the canonical frame const uint64_t window_begin = m_min_position - m_min_pos_in_kmer; uint64_t pos_in_kmer = m_min_pos_in_kmer; @@ -141,9 +176,7 @@ struct minimizer_iterator { Kmer window = kmer; window.drop_chars(m_min_pos_in_kmer + 1); for (uint64_t i = m_min_pos_in_kmer + 1; i != m_k - m_m + 1; ++i) { - Kmer mmer = window; - mmer.take_chars(m_m); - uint64_t v = util::canonical_mmer(uint64_t(mmer), m_m); + uint64_t v = util::canonical_mmer_at(window, m_kmer_rc, m_k, m_m, i); if (m_hasher.hash(v) == m_min_hash) { // rightmost pos_in_kmer = i; value = v; diff --git a/include/util.hpp b/include/util.hpp index af2415a..9561dd9 100644 --- a/include/util.hpp +++ b/include/util.hpp @@ -268,6 +268,48 @@ inline uint64_t canonical_mmer(const uint64_t mmer, const uint64_t m) { } } +/* + The canonical m-mer at locus i of a kmer x, given rc(x). + + Uses rc(x_i) = rc(x)_{k-m-i}: the reverse complements of all k-m+1 loci are + windows of the single reverse complement of the whole kmer, so the kmer is + reverse-complemented once instead of once per locus. `window` is x shifted so + that locus i sits at its low end -- the forward sliding window the callers + already carry along. + + When the whole kmer fits in a word -- always so for the 64-bit kmer type, and + for k <= 32 with the wider one -- the reverse window is extracted with a plain + 64-bit shift; otherwise with a shift of the wide type, which costs a few + instructions more but still beats reverse-complementing every locus. +*/ +template +inline uint64_t canonical_mmer_at(Kmer window, Kmer const& kmer_rc, const uint64_t k, + const uint64_t m, const uint64_t i) // +{ + window.take_chars(m); + if constexpr (!Kmer::has_reverse_complement) { + (void)kmer_rc; + (void)k; + (void)i; + return uint64_t(window); + } else { + constexpr uint64_t b = Kmer::bits_per_char; + assert(i + m <= k); + uint64_t rc; + if (k * b <= 64) { + const uint64_t mask = m * b == 64 ? ~uint64_t(0) : (uint64_t(1) << (m * b)) - 1; + rc = (uint64_t(kmer_rc) >> (b * (k - m - i))) & mask; + } else { + Kmer mmer_rc = kmer_rc; + mmer_rc.drop_chars(k - m - i); + mmer_rc.take_chars(m); + rc = uint64_t(mmer_rc); + } + assert(rc == Kmer::reverse_complement_mmer(uint64_t(window), m)); + return std::min(uint64_t(window), rc); + } +} + /* True if `kmer` is the canonical one of the pair (kmer, rc(kmer)). A kmer that equals its own reverse complement (possible only for even k) is deemed @@ -309,21 +351,41 @@ inline bool is_canonical(Kmer kmer, const uint64_t k) { min(x, rc(x)), which is literally the same string for x and rc(x): that amounts to taking the leftmost tied locus when x is canonical and the rightmost one otherwise. + + Since rc(x_i) = rc(x)_{k-m-i}, the reverse complements of all k-m+1 loci are + just windows of the single reverse complement of the whole kmer: `kmer_rc` is + reverse-complemented once and each rc(x_i) is read out of it with a shift, + instead of reverse-complementing every locus separately. The caller usually + has rc(x) already -- a lookup needs it anyway to recognise which orientation + it found -- in which case the reverse complementation is free. */ template -minimizer_info compute_minimizer(Kmer kmer, const uint64_t k, const uint64_t m, +minimizer_info compute_minimizer(Kmer kmer, Kmer const& kmer_rc, const uint64_t k, const uint64_t m, hasher_type const& hasher) // { assert(m <= Kmer::max_m); assert(m <= k); + /* + The canonical m-mer at locus i. The forward window slides left to right, + so it is carried along; the reverse one runs right to left over kmer_rc, + so it is extracted with a shift rather than slid. + + When the whole kmer fits in a word -- always so for the 64-bit kmer type, + and for k <= 32 with the wider one -- that shift is done on a plain + uint64_t. Otherwise it has to be a shift of the wide type, which costs a + few instructions more but is still cheaper than reverse-complementing + every locus separately. + */ + auto kappa = [&](Kmer const& window, const uint64_t i) { + return canonical_mmer_at(window, kmer_rc, k, m, i); + }; + /* The first locus is peeled off the loop so that `min_hash` starts out at a real hash value: initializing it to invalid_uint64 would make an actual hash of invalid_uint64 register as a tie rather than as the minimum. */ Kmer window = kmer; - Kmer first = window; - first.take_chars(m); - uint64_t minimizer = canonical_mmer(uint64_t(first), m); + uint64_t minimizer = kappa(window, 0); uint64_t min_hash = hasher.hash(minimizer); uint64_t leftmost = 0; uint64_t rightmost = 0; @@ -331,9 +393,7 @@ minimizer_info compute_minimizer(Kmer kmer, const uint64_t k, const uint64_t m, window.drop_char(); for (uint64_t i = 1; i != k - m + 1; ++i) { - Kmer mmer = window; - mmer.take_chars(m); - uint64_t value = canonical_mmer(uint64_t(mmer), m); + uint64_t value = kappa(window, i); uint64_t hash = hasher.hash(value); if (hash < min_hash) { // leftmost min_hash = hash; @@ -348,16 +408,24 @@ minimizer_info compute_minimizer(Kmer kmer, const uint64_t k, const uint64_t m, window.drop_char(); } - if (tie and !is_canonical(kmer, k)) { + if (tie and kmer_rc < kmer) { /* rc(kmer) is the canonical frame: mirror its leftmost tied locus */ kmer.drop_chars(rightmost); - kmer.take_chars(m); - return {canonical_mmer(uint64_t(kmer), m), rightmost}; + return {kappa(kmer, rightmost), rightmost}; } return {minimizer, leftmost}; } +template +minimizer_info compute_minimizer(Kmer kmer, const uint64_t k, const uint64_t m, + hasher_type const& hasher) // +{ + Kmer kmer_rc = kmer; + kmer_rc.reverse_complement_inplace(k); + return compute_minimizer(kmer, kmer_rc, k, m, hasher); +} + } // namespace util struct buffered_lines_iterator { diff --git a/src/dictionary.cpp b/src/dictionary.cpp index 748d0c4..ae5fa5d 100644 --- a/src/dictionary.cpp +++ b/src/dictionary.cpp @@ -16,7 +16,9 @@ lookup_result dictionary::lookup(Kmer uint_kmer, { Kmer uint_kmer_rc = uint_kmer; uint_kmer_rc.reverse_complement_inplace(m_k); - auto mini_info = util::compute_minimizer(uint_kmer, m_k, m_m, m_hasher); + /* the reverse complement is needed anyway to tell the two orientations + apart, so computing the minimizer from it costs nothing extra */ + auto mini_info = util::compute_minimizer(uint_kmer, uint_kmer_rc, m_k, m_m, m_hasher); auto res = lookup(uint_kmer, uint_kmer_rc, mini_info); if (!check_reverse_complement and res.kmer_orientation == constants::backward_orientation) { // return lookup_result(); @@ -29,7 +31,7 @@ lookup_result dictionary::lookup(const Kmer uint_kmer, const minimizer_info mini_info) const // { assert(minimizer_info(mini_info.minimizer, mini_info.pos_in_kmer) == - util::compute_minimizer(uint_kmer, m_k, m_m, m_hasher)); + util::compute_minimizer(uint_kmer, uint_kmer_rc, m_k, m_m, m_hasher)); auto it = m_ssi.lookup(uint_kmer, uint_kmer_rc, mini_info); return m_spss.lookup(it, uint_kmer, uint_kmer_rc, mini_info); } From f57aa7dfde885f5f88fca15110dd2d04f63f9146 Mon Sep 17 00:00:00 2001 From: Giulio Ermanno Pibiri Date: Mon, 10 Aug 2026 22:46:43 +0000 Subject: [PATCH 05/14] drop the incremental reverse-complement slide from minimizer_iterator The previous commit had the iterator carry the kmer's reverse complement as state and slide it one character at a time, on the assumption that this beat recomputing it. It does, but by much less than it costs to read: the whole mechanism -- Kmer::complement_char with its two encodings, the m_fresh flag, slide_reverse_complement, an assertion re-deriving the value to guard the "kmers must be consecutive" contract it introduced -- bought 3 percentage points on one build step. Measured on salmonella_100 k=31 m=13, step 2, min of 12 interleaved runs: before the optimization 386937 with the slide 352011 -9.0% recomputing rc per kmer 364386 -5.8% so the slide is worth about 3% of step 2, which is under 1% of total build time, for roughly thirty lines and a stateful invariant. Recomputing wins on simplicity: rc(x) is now a local in next(), passed down to rescan and break_tie, and minimizer_iterator ends up four code lines shorter than it was before the optimization rather than twenty-four longer. (An isolated microbenchmark of the iterator had suggested recomputing was the faster of the two. It was measuring the wrong thing: it rebuilt each kmer with string_to_uint_kmer per call, which dominated. The build numbers above are what counts.) Also trims two bits of incidental bulk in util.hpp: the (void) casts in canonical_mmer_at become [[maybe_unused]] parameters, and a comment block above the kappa lambda that restated what canonical_mmer_at already documents is gone. The net cost of the optimization is now +27 lines of real code, all of it canonical_mmer_at and the four-argument compute_minimizer overload, against +97 before this commit. Index output unchanged. Verified with --check on salmonella_enterica (m in {11,13,20,31}), ecoli1 (4 threads), --weighted, and se.ust.k47/k63, with assertions enabled. --- include/kmer.hpp | 15 ------- include/minimizer_iterator.hpp | 71 +++++++++------------------------- include/util.hpp | 19 ++------- 3 files changed, 22 insertions(+), 83 deletions(-) diff --git a/include/kmer.hpp b/include/kmer.hpp index 4db3d00..3dae600 100644 --- a/include/kmer.hpp +++ b/include/kmer.hpp @@ -114,8 +114,6 @@ struct alpha_kmer_t : uint_kmer_t { return mmer; } - /* Complement of a single packed character. Same fallback: the identity. */ - [[maybe_unused]] static uint64_t complement_char(uint64_t c) { return c; } [[maybe_unused]] static void compute_reverse_complement(char const* input, char* output, uint64_t size) { for (uint64_t i = 0; i != size; ++i) output[i] = input[i]; @@ -190,19 +188,6 @@ struct dna_uint_kmer_t : alpha_kmer_t { return crc64(mmer) >> (64 - m * bits_per_char); } - /* - Complement of a single packed character: the per-character half of what - crc64 does to a whole word. Lets a reverse complement be slid one - character at a time instead of recomputed. - */ - [[maybe_unused]] static uint64_t complement_char(uint64_t c) { -#ifdef SSHASH_USE_TRADITIONAL_NUCLEOTIDE_ENCODING - return c ^ 3; // A<->T is 00<->11, C<->G is 01<->10 -#else - return c ^ 2; // A<->T is 00<->10, C<->G is 01<->11 -#endif - } - #ifdef SSHASH_USE_TRADITIONAL_NUCLEOTIDE_ENCODING /* char decimal binary diff --git a/include/minimizer_iterator.hpp b/include/minimizer_iterator.hpp index 37a81fa..cf570cc 100644 --- a/include/minimizer_iterator.hpp +++ b/include/minimizer_iterator.hpp @@ -13,17 +13,14 @@ namespace sshash { iterator computes exactly the same thing incrementally, which the assertion at the end of `next` checks. - The extra state compared to a plain forward minimizer is `m_num_mins`, the - number of loci of the current window attaining the minimum hash -- a tie + The only extra state compared to a plain forward minimizer is `m_num_mins`, + the number of loci of the current window attaining the minimum hash: a tie cannot be broken by position alone without breaking mirror-equivariance, so - when there is one we have to look at the kmer's own orientation -- and - `m_kmer_rc`, the reverse complement of the current kmer. - - `m_kmer_rc` serves every locus at once, via rc(x_i) = rc(x)_{k-m-i}, so no - m-mer is ever reverse-complemented on its own; and it is itself slid one - character at a time rather than recomputed. Like the window minimum, it - assumes the kmers handed to `next` are consecutive, which `set_position` and - `reset` restart. + when there is one we have to look at the kmer's own orientation. + + The kmer is reverse-complemented once per call and that single value serves + every locus, via rc(x_i) = rc(x)_{k-m-i}, so no m-mer is ever + reverse-complemented on its own. */ template struct minimizer_iterator { @@ -49,21 +46,21 @@ struct minimizer_iterator { m_min_pos_in_kmer = 0; m_min_position = m_position - 1; m_num_mins = 0; - m_fresh = true; } minimizer_info next(Kmer kmer) { - slide_reverse_complement(kmer); + Kmer kmer_rc = kmer; + kmer_rc.reverse_complement_inplace(m_k); if (m_min_pos_in_kmer == 0) { /* min leaves the window: re-scan to compute the new min */ m_position = m_min_position + 1; - rescan(kmer); + rescan(kmer, kmer_rc); } else { m_position += 1; Kmer window = kmer; window.drop_chars(m_k - m_m); - uint64_t value = util::canonical_mmer_at(window, m_kmer_rc, m_k, m_m, m_k - m_m); + uint64_t value = util::canonical_mmer_at(window, kmer_rc, m_k, m_m, m_k - m_m); uint64_t hash = m_hasher.hash(value); if (hash < m_min_hash) { m_min_hash = hash; @@ -81,10 +78,10 @@ struct minimizer_iterator { } minimizer_info mini_info(m_min_value, m_min_position, m_min_pos_in_kmer); - if (m_num_mins > 1) break_tie(kmer, mini_info); + if (m_num_mins > 1) break_tie(kmer, kmer_rc, mini_info); assert(minimizer_info(mini_info.minimizer, mini_info.pos_in_kmer) == - util::compute_minimizer(kmer, m_kmer_rc, m_k, m_m, m_hasher)); + util::compute_minimizer(kmer, kmer_rc, m_k, m_m, m_hasher)); return mini_info; } @@ -94,44 +91,14 @@ struct minimizer_iterator { uint64_t m_position, m_min_pos_in_kmer; uint64_t m_min_value, m_min_position, m_min_hash; uint64_t m_num_mins; - Kmer m_kmer_rc; - bool m_fresh; hasher_type m_hasher; - /* - rc(x) for the kmer just handed in. After a restart it is computed - outright; otherwise the kmer has slid by one character, so its reverse - complement has too: the new character's complement enters at the front - and the oldest one falls off the back. - */ - void slide_reverse_complement(Kmer const& kmer) { - if constexpr (!Kmer::has_reverse_complement) { - (void)kmer; - return; - } else { - if (m_fresh) { - m_kmer_rc = kmer; - m_kmer_rc.reverse_complement_inplace(m_k); - m_fresh = false; - } else { - m_kmer_rc.pad_char(); - m_kmer_rc.set(0, Kmer::complement_char(kmer.at(m_k - 1))); - m_kmer_rc.take(m_k * Kmer::bits_per_char); - } - assert([&] { - Kmer expected = kmer; - expected.reverse_complement_inplace(m_k); - return expected == m_kmer_rc; - }()); - } - } - - void rescan(Kmer kmer) { + void rescan(Kmer kmer, Kmer const& kmer_rc) { const uint64_t begin = m_position; /* first locus, peeled off the loop: see `util::compute_minimizer` */ { - m_min_value = util::canonical_mmer_at(kmer, m_kmer_rc, m_k, m_m, 0); + m_min_value = util::canonical_mmer_at(kmer, kmer_rc, m_k, m_m, 0); kmer.drop_char(); m_min_hash = m_hasher.hash(m_min_value); m_min_pos_in_kmer = 0; @@ -140,7 +107,7 @@ struct minimizer_iterator { } for (uint64_t i = 1; i != m_k - m_m + 1; ++i, ++m_position) { - uint64_t value = util::canonical_mmer_at(kmer, m_kmer_rc, m_k, m_m, i); + uint64_t value = util::canonical_mmer_at(kmer, kmer_rc, m_k, m_m, i); kmer.drop_char(); uint64_t hash = m_hasher.hash(value); if (hash < m_min_hash) { // leftmost @@ -165,9 +132,9 @@ struct minimizer_iterator { It runs on ~1e-5 of the windows. */ - void break_tie(Kmer kmer, minimizer_info& mini_info) const { + void break_tie(Kmer kmer, Kmer const& kmer_rc, minimizer_info& mini_info) const { assert(m_num_mins > 1); - if (!(m_kmer_rc < kmer)) return; // the kmer is already the canonical frame + if (!(kmer_rc < kmer)) return; // the kmer is already the canonical frame const uint64_t window_begin = m_min_position - m_min_pos_in_kmer; uint64_t pos_in_kmer = m_min_pos_in_kmer; @@ -176,7 +143,7 @@ struct minimizer_iterator { Kmer window = kmer; window.drop_chars(m_min_pos_in_kmer + 1); for (uint64_t i = m_min_pos_in_kmer + 1; i != m_k - m_m + 1; ++i) { - uint64_t v = util::canonical_mmer_at(window, m_kmer_rc, m_k, m_m, i); + uint64_t v = util::canonical_mmer_at(window, kmer_rc, m_k, m_m, i); if (m_hasher.hash(v) == m_min_hash) { // rightmost pos_in_kmer = i; value = v; diff --git a/include/util.hpp b/include/util.hpp index 9561dd9..6adac4d 100644 --- a/include/util.hpp +++ b/include/util.hpp @@ -283,14 +283,12 @@ inline uint64_t canonical_mmer(const uint64_t mmer, const uint64_t m) { instructions more but still beats reverse-complementing every locus. */ template -inline uint64_t canonical_mmer_at(Kmer window, Kmer const& kmer_rc, const uint64_t k, - const uint64_t m, const uint64_t i) // +inline uint64_t canonical_mmer_at(Kmer window, [[maybe_unused]] Kmer const& kmer_rc, + [[maybe_unused]] const uint64_t k, const uint64_t m, + [[maybe_unused]] const uint64_t i) // { window.take_chars(m); if constexpr (!Kmer::has_reverse_complement) { - (void)kmer_rc; - (void)k; - (void)i; return uint64_t(window); } else { constexpr uint64_t b = Kmer::bits_per_char; @@ -366,17 +364,6 @@ minimizer_info compute_minimizer(Kmer kmer, Kmer const& kmer_rc, const uint64_t assert(m <= Kmer::max_m); assert(m <= k); - /* - The canonical m-mer at locus i. The forward window slides left to right, - so it is carried along; the reverse one runs right to left over kmer_rc, - so it is extracted with a shift rather than slid. - - When the whole kmer fits in a word -- always so for the 64-bit kmer type, - and for k <= 32 with the wider one -- that shift is done on a plain - uint64_t. Otherwise it has to be a shift of the wide type, which costs a - few instructions more but is still cheaper than reverse-complementing - every locus separately. - */ auto kappa = [&](Kmer const& window, const uint64_t i) { return canonical_mmer_at(window, kmer_rc, k, m, i); }; From b1d0706825e34e852f668c844469ab2bcd4d9497 Mon Sep 17 00:00:00 2001 From: Giulio Ermanno Pibiri Date: Tue, 18 Aug 2026 19:12:01 +0000 Subject: [PATCH 06/14] centre-closest tie-break (Proposition 24): the scheme is now forward Replace the tie-break of the canonical minimizer with the one of [Cologni and Pibiri, "Canonical Schemes and Minimizers", Proposition 24]: among the loci tied at the minimum, take the one closest to the centre (k-m)/2 of the window; when two are equally close -- they are then mirror images i and k-m-i -- take the smaller index if x <= rc(x) and the larger otherwise. The previous rule (leftmost in the frame of the canonical kmer) was mirror-equivariant but not forward: the anchor could move backwards on a tie, re-opening a super-kmer at a position already paid for. The centre rule is both. Equivariance holds because the mirror fixes the centre, so distances to it are strand-independent; forwardness because the distance of a fixed locus to the centre is monotone as the window slides right, so a locus at least as close as another stays at least as close, and the anchor never decreases (Proposition 24 of the paper). No leftmost rule could have worked: the mirror reverses the order of the tied loci, so "leftmost" names different loci on the two strands (Proposition 23). Confirmed: the number of super-kmers now equals the number of sampled positions, exactly, on every dataset and every m tried -- before (gap) after (gap) salmonella_enterica m=11 435321 (+897) 434313 (+0) m=13 482031 (+220) 481761 (+0) m=17 598144 (+34) 598098 (+0) m=20 734279 (+11) 734265 (+0) ecoli1 m=11..20 gaps +693..+20 all +0 salmonella_100 m=11..20 gaps +2528..+78 all +0 se.ust.k63 (128b) m=21 -- 223765 = 223765 and test/test_minimizer.cpp now asserts forwardness directly: the sampled position never decreases as the iterator slides. Restoring the old rule fails that check immediately (first violation at k=5), so the test has teeth. Mirror-equivariance is checked as before, over k in {5..63} and every valid m, tie rate 16.5% of sampled kmers at small m. The number of sampled positions itself also drops slightly (e.g. 434424 -> 434313 at m=11): the centre anchor is more stable than the frame-dependent one, so fewer distinct loci are ever selected. Performance: the hot loop is unchanged -- same comparisons, same updates; only the tie path (rate ~ sigma^{-m/2}/w, about 1e-4 of windows at m=13) differs. One codegen trap mattered: with the tie scan written inline, gcc let it bloat the body of dictionary::lookup and positive lookups regressed by 15% even though the scan executes almost never (and compute_minimizer in isolation had gotten *faster*, 59ns vs 79ns at m=13 -- the old tie path kept more state live). Outlining it as a noinline resolve_tie helper behind __builtin_expect restores parity: step 2 of the build measures +0.2%, and lookups are within this machine's noise floor (an A/A test of the same binary against itself shows +-4% on minima, and old-vs-new sits inside that: pos +5.2%, neg -1.1% on minima over 15 interleaved pairs). The sampled positions of tied windows differ from those of version 6 indexes, so an old index queried with this code would miss its tied kmers: index version 6.0.0 -> 7.0.0, existing indexes need rebuilding. Verified with --check (all six suites) on salmonella_enterica (m in {11,13,20,31}), ecoli1 (4 threads), a --weighted build, and se.ust.k47/k63, with assertions enabled. --- include/builder/util.hpp | 10 ++--- include/constants.hpp | 2 +- include/minimizer_iterator.hpp | 44 +++++++++++++-------- include/util.hpp | 72 ++++++++++++++++++++++++++++------ test/test_minimizer.cpp | 27 ++++++++++--- 5 files changed, 116 insertions(+), 39 deletions(-) diff --git a/include/builder/util.hpp b/include/builder/util.hpp index 2cdbe94..55de9b8 100644 --- a/include/builder/util.hpp +++ b/include/builder/util.hpp @@ -93,16 +93,16 @@ struct bucket_type { iterator end() const { return iterator(m_end); } /* - A minimizer offset can correspond to more than one super-kmer: the - minimizer is not strictly forward (see `util::compute_minimizer`), so a - locus can be abandoned and later re-selected. A super-kmer is uniquely identified by the couple (minimizer offset, position of minimizer in the first kmer of the super-kmer). These two components, together, give the starting position of a super-kmer in the sequence. - So the method size() returns the number of minimizer - positions which is <= the number of superkmers. + The minimizer is forward (see `util::compute_minimizer`), so a locus is + never abandoned and later re-selected, and the number of super-kmers + equals the number of minimizer positions. The code below does not rely + on that and stays correct for a non-forward scheme too, where size() -- + the number of minimizer positions -- is < the number of super-kmers. */ uint64_t num_super_kmers() const { return m_num_super_kmers; } uint64_t size() const { return m_num_minimizer_positions; } diff --git a/include/constants.hpp b/include/constants.hpp index 215901a..b5f214d 100644 --- a/include/constants.hpp +++ b/include/constants.hpp @@ -20,7 +20,7 @@ constexpr int forward_orientation = 1; constexpr int backward_orientation = -1; namespace current_version_number { -constexpr uint8_t x = 6; +constexpr uint8_t x = 7; constexpr uint8_t y = 0; constexpr uint8_t z = 0; } // namespace current_version_number diff --git a/include/minimizer_iterator.hpp b/include/minimizer_iterator.hpp index cf570cc..56f26f0 100644 --- a/include/minimizer_iterator.hpp +++ b/include/minimizer_iterator.hpp @@ -125,33 +125,45 @@ struct minimizer_iterator { } /* - Two or more loci of the window attain the minimum hash, and the leftmost - of them is the one currently held. Leftmost is the right answer when the - kmer is canonical; otherwise the canonical frame is rc(kmer), whose - leftmost tied locus is this window's rightmost one. - - It runs on ~1e-5 of the windows. + Two or more loci of the window attain the minimum hash (all carrying the + same class, so only the sampled position is at stake). Apply the + centre-closest tie-break of `util::compute_minimizer`: take the tied + locus closest to the window centre, and between two equally close ones + -- mirror images i and k-m-i -- the smaller index if kmer <= rc(kmer), + the larger otherwise. Forward and mirror-equivariant. It runs on ~1e-4 + of the windows, so the rescan below costs nothing overall. + + Tied loci cannot precede the leftmost minimum, so the scan starts there. */ void break_tie(Kmer kmer, Kmer const& kmer_rc, minimizer_info& mini_info) const { assert(m_num_mins > 1); - if (!(kmer_rc < kmer)) return; // the kmer is already the canonical frame const uint64_t window_begin = m_min_position - m_min_pos_in_kmer; - uint64_t pos_in_kmer = m_min_pos_in_kmer; - uint64_t value = m_min_value; + const uint64_t two_c = m_k - m_m; + uint64_t best_dist = constants::invalid_uint64; + uint64_t lo = m_min_pos_in_kmer; + uint64_t hi = m_min_pos_in_kmer; Kmer window = kmer; - window.drop_chars(m_min_pos_in_kmer + 1); - for (uint64_t i = m_min_pos_in_kmer + 1; i != m_k - m_m + 1; ++i) { - uint64_t v = util::canonical_mmer_at(window, kmer_rc, m_k, m_m, i); - if (m_hasher.hash(v) == m_min_hash) { // rightmost - pos_in_kmer = i; - value = v; + window.drop_chars(m_min_pos_in_kmer); + for (uint64_t i = m_min_pos_in_kmer; i != m_k - m_m + 1; ++i) { + const uint64_t v = util::canonical_mmer_at(window, kmer_rc, m_k, m_m, i); + if (m_hasher.hash(v) == m_min_hash) { + assert(v == m_min_value); + const uint64_t dist = 2 * i > two_c ? 2 * i - two_c : two_c - 2 * i; + if (dist < best_dist) { + best_dist = dist; + lo = i; + hi = i; + } else if (dist == best_dist) { + hi = i; + } } window.drop_char(); } - mini_info = minimizer_info(value, window_begin + pos_in_kmer, pos_in_kmer); + const uint64_t chosen = (lo == hi or !(kmer_rc < kmer)) ? lo : hi; + mini_info = minimizer_info(m_min_value, window_begin + chosen, chosen); } }; diff --git a/include/util.hpp b/include/util.hpp index 6adac4d..9d4b18c 100644 --- a/include/util.hpp +++ b/include/util.hpp @@ -342,13 +342,23 @@ inline bool is_canonical(Kmer kmer, const uint64_t k) { When the alphabet has no reverse complement, kappa is the identity and this reduces exactly to the plain forward minimizer. - Ties -- h(kappa(i)) == h(kappa(j)) for i != j, which happens when x_i == x_j - or x_i == rc(x_j) -- must be broken in a mirror-equivariant way, or x and + Ties -- h(kappa(i)) == h(kappa(j)) for i != j, which happens when x_i = x_j + or x_i = rc(x_j) -- must be broken in a mirror-equivariant way, or x and rc(x) would be sent to different buckets, which is a correctness failure and - not merely a density one. We break them in the frame of the canonical kmer - min(x, rc(x)), which is literally the same string for x and rc(x): that - amounts to taking the leftmost tied locus when x is canonical and the - rightmost one otherwise. + not merely a density one. No leftmost rule can be equivariant: the mirror + reverses the order of the tied loci, so "leftmost" names different loci on + the two strands. We use the tie-break of [Cologni and Pibiri, "Canonical + Schemes and Minimizers", Proposition 24]: among the tied loci, take the one + closest to the centre (k-m)/2 of the window -- the one reference point both + strands agree on -- and when two are equally close (they are then mirror + images i and k-m-i), take the smaller index if x <= rc(x) and the larger + otherwise. This rule is mirror-equivariant (for odd k, where x != rc(x) + always) and, unlike the previous frame-of-the-canonical-kmer rule, it is + also *forward*: the distance of a fixed locus to the centre is monotone as + the window slides right, so a locus at least as close as another stays so, + and the anchor never moves backwards. Forwardness makes the number of + super-kmers equal the number of sampled positions -- nothing is ever + abandoned and re-opened. Since rc(x_i) = rc(x)_{k-m-i}, the reverse complements of all k-m+1 loci are just windows of the single reverse complement of the whole kmer: `kmer_rc` is @@ -357,6 +367,49 @@ inline bool is_canonical(Kmer kmer, const uint64_t k) { has rc(x) already -- a lookup needs it anyway to recognise which orientation it found -- in which case the reverse complementation is free. */ +/* + A tie at the minimum (rare, rate ~ sigma^{-m/2}): all tied loci carry the + same class, so only the sampled position is at stake, not the value. Take + the tied locus closest to the centre of the window; the doubled distance + |2i - (k-m)| avoids the half-integer centre of an odd window. + + Deliberately not inlined: this runs on ~1e-4 of the windows, and letting its + loop inline into the lookup path costs measurably more than the call. +*/ +template +__attribute__((noinline)) minimizer_info resolve_tie(Kmer const& kmer, Kmer const& kmer_rc, + const uint64_t k, const uint64_t m, + hasher_type const& hasher, + const uint64_t min_hash, + const uint64_t minimizer, + const uint64_t leftmost, + const uint64_t rightmost) // +{ + const uint64_t two_c = k - m; + uint64_t best_dist = constants::invalid_uint64; + uint64_t lo = leftmost; + uint64_t hi = leftmost; + Kmer window = kmer; + window.drop_chars(leftmost); + for (uint64_t i = leftmost; i <= rightmost; ++i) { + if (hasher.hash(canonical_mmer_at(window, kmer_rc, k, m, i)) == min_hash) { + assert(canonical_mmer_at(window, kmer_rc, k, m, i) == minimizer); + const uint64_t dist = 2 * i > two_c ? 2 * i - two_c : two_c - 2 * i; + if (dist < best_dist) { + best_dist = dist; + lo = i; + hi = i; + } else if (dist == best_dist) { + hi = i; + } + } + window.drop_char(); + } + /* lo and hi are the two equally-closest loci (mirror images), or one locus */ + const uint64_t chosen = (lo == hi or !(kmer_rc < kmer)) ? lo : hi; + return {minimizer, chosen}; +} + template minimizer_info compute_minimizer(Kmer kmer, Kmer const& kmer_rc, const uint64_t k, const uint64_t m, hasher_type const& hasher) // @@ -395,12 +448,9 @@ minimizer_info compute_minimizer(Kmer kmer, Kmer const& kmer_rc, const uint64_t window.drop_char(); } - if (tie and kmer_rc < kmer) { - /* rc(kmer) is the canonical frame: mirror its leftmost tied locus */ - kmer.drop_chars(rightmost); - return {kappa(kmer, rightmost), rightmost}; + if (__builtin_expect(tie, false)) { + return resolve_tie(kmer, kmer_rc, k, m, hasher, min_hash, minimizer, leftmost, rightmost); } - return {minimizer, leftmost}; } diff --git a/test/test_minimizer.cpp b/test/test_minimizer.cpp index 81e2b49..484bc0d 100644 --- a/test/test_minimizer.cpp +++ b/test/test_minimizer.cpp @@ -43,8 +43,8 @@ static void check_equivariance(kmer_t kmer, uint64_t k, uint64_t m, hasher_type return; } if (a.pos_in_kmer + b.pos_in_kmer != k - m) { - fail("locus of kmer '" + util::uint_kmer_to_string(kmer, k) + - "' is " + std::to_string(a.pos_in_kmer) + " but that of its reverse complement is " + + fail("locus of kmer '" + util::uint_kmer_to_string(kmer, k) + "' is " + + std::to_string(a.pos_in_kmer) + " but that of its reverse complement is " + std::to_string(b.pos_in_kmer) + " (should mirror to " + std::to_string(k - m - a.pos_in_kmer) + ")"); return; @@ -93,8 +93,14 @@ static void test_kmers(uint64_t k, uint64_t m, uint64_t num_kmers, std::mt19937_ /* The incremental "re-scan" iterator must agree with the brute-force reference - on every window of a sequence, including the windows where a tie fires and - the anchor consequently moves backwards. + on every window of a sequence, including the windows where a tie fires. + + The scheme must also be *forward*: the sampled position must never decrease + as the window slides, so that the number of super-kmers equals the number of + sampled positions. The centre-closest tie-break is what guarantees this + (leftmost-if-canonical / rightmost-otherwise, the previous rule, is + mirror-equivariant too but moves the anchor backwards on ~1e-5 of the + windows). */ template static void test_iterator(uint64_t k, uint64_t m, uint64_t length, std::mt19937_64& gen) { @@ -105,6 +111,7 @@ static void test_iterator(uint64_t k, uint64_t m, uint64_t length, std::mt19937_ minimizer_iterator it(k, m, hasher); it.set_position(0); + uint64_t prev_pos_in_seq = 0; for (uint64_t i = 0; i + k <= length; ++i) { kmer_t kmer = util::string_to_uint_kmer(s.data() + i, k); auto got = it.next(kmer); @@ -117,10 +124,18 @@ static void test_iterator(uint64_t k, uint64_t m, uint64_t length, std::mt19937_ } /* pos_in_seq must be the absolute position of the selected locus */ if (got.pos_in_seq != i + got.pos_in_kmer) { - fail("iterator reports pos_in_seq " + std::to_string(got.pos_in_seq) + " but expected " + - std::to_string(i + got.pos_in_kmer)); + fail("iterator reports pos_in_seq " + std::to_string(got.pos_in_seq) + + " but expected " + std::to_string(i + got.pos_in_kmer)); return; } + /* forwardness: the sampled position never decreases */ + if (got.pos_in_seq < prev_pos_in_seq) { + fail("scheme is not forward: sampled position " + std::to_string(got.pos_in_seq) + + " after " + std::to_string(prev_pos_in_seq) + " (k=" + std::to_string(k) + + ", m=" + std::to_string(m) + ")"); + return; + } + prev_pos_in_seq = got.pos_in_seq; } } From 2cafe284dcf2eb31d92e52ab5a15c35b2aeed4e2 Mon Sep 17 00:00:00 2001 From: Giulio Ermanno Pibiri Date: Wed, 19 Aug 2026 20:03:28 +0000 Subject: [PATCH 07/14] drop num_super_kmers: one tuple per minimizer position, checked at merge With a forward scheme, the number of super-kmers equals the number of minimizer positions by construction: a tuple is a maximal run of consecutive kmers anchored at the same pos_in_seq, forwardness means a position is never re-selected after being abandoned (so one run per position), and positions are absolute offsets in the global strings (so no cross-sequence collisions). Index space was proportional to num_minimizer_positions all along -- mid/heavy load buckets store one entry per distinct position, the skew index one per kmer, the codewords one per minimizer -- so num_super_kmers only ever named the record count of the tuples file, which now coincides with the position count. Keeping both was redundancy. What remains of the quantity is its one useful residue: a forwardness check. `minimizers_tuples::merge` now verifies, while counting, that no (minimizer, pos_in_seq) pair occurs twice among the sorted tuples, and throws otherwise -- everything downstream sizes the index by the tuple count, so a violation must stop the build rather than corrupt it. The check is one comparison per record inside a disk-bound loop. Verified to fire: rebuilding with the old non-forward tie-break aborts with "the minimizer scheme is not forward". Simplifications that follow: - bucket_type no longer scans its tuples to count distinct positions: size() is end - begin, O(1) instead of O(bucket) at every construction (it is constructed per bucket in each of the three passes over the merged file). The distinctness scan survives as a debug assert. - minimizers_tuples loses m_num_super_kmers and its accessor; step 5 reads the tuples-file record count from num_minimizer_positions(). - step 7.1: the singleton branch reads its one tuple directly; the mid/heavy write loops push one position per tuple, with the codeword set once from the first tuple, instead of dedup-guarded loops; the transient `tuples` buffer is reserved from the two position counters instead of a third super-kmer counter. - step 7.2: pos_in_bucket increments per tuple, no prev_pos_in_seq. - compute_minimizer_tuples: the buffer-merging branch in save() is replaced by an assert -- it was dead code even for the old scheme, since consecutive saves always differ in pos_in_seq by the break test. Net -32 lines. Behaviour-preserving: the produced index is byte-identical to the parent commit's on salmonella_enterica and salmonella_100 (m=13 and m=20) and on ecoli1 with 4 threads, which exercises the multi-file merge path. Full --check matrix passes with assertions enabled, including --weighted and se.ust.k47/k63. --- include/builder/dictionary_builder.hpp | 16 +-- include/builder/util.hpp | 107 ++++++++++---------- src/builder/build_sparse_and_skew_index.cpp | 74 +++++--------- src/builder/compute_minimizer_tuples.cpp | 11 +- 4 files changed, 88 insertions(+), 120 deletions(-) diff --git a/include/builder/dictionary_builder.hpp b/include/builder/dictionary_builder.hpp index 309343f..5f63d54 100644 --- a/include/builder/dictionary_builder.hpp +++ b/include/builder/dictionary_builder.hpp @@ -48,7 +48,6 @@ struct dictionary_builder // std::cout << "num_minimizers = " << minimizers.num_minimizers() << std::endl; std::cout << "num_minimizer_positions = " << minimizers.num_minimizer_positions() << std::endl; - std::cout << "num_super_kmers = " << minimizers.num_super_kmers() << std::endl; } do_step("step 4 (build mphf)", [&]() { build_mphf(d); }); @@ -141,21 +140,22 @@ struct dictionary_builder // } } - const uint64_t num_super_kmers = minimizers.num_super_kmers(); + /* one tuple per minimizer position (the scheme is forward), so this is + the number of records in the tuples file */ + const uint64_t num_tuples = minimizers.num_minimizer_positions(); const uint64_t buffer_size = num_files_to_merge == 1 - ? num_super_kmers + ? num_tuples : (RAM_available_in_bytes / (3 * sizeof(minimizer_tuple))); - const uint64_t num_blocks = (num_super_kmers + buffer_size - 1) / buffer_size; - assert(num_super_kmers > (num_blocks - 1) * buffer_size); + const uint64_t num_blocks = (num_tuples + buffer_size - 1) / buffer_size; + assert(num_tuples > (num_blocks - 1) * buffer_size); std::vector threads; threads.reserve(num_threads); std::vector buffer; for (uint64_t i = 0; i != num_blocks; ++i) { - const uint64_t n = (i == num_blocks - 1) - ? num_super_kmers - (num_blocks - 1) * buffer_size - : buffer_size; + const uint64_t n = + (i == num_blocks - 1) ? num_tuples - (num_blocks - 1) * buffer_size : buffer_size; buffer.resize(n); input.read(reinterpret_cast(buffer.data()), buffer.size() * sizeof(minimizer_tuple)); diff --git a/include/builder/util.hpp b/include/builder/util.hpp index 55de9b8..add5e00 100644 --- a/include/builder/util.hpp +++ b/include/builder/util.hpp @@ -58,23 +58,24 @@ inline std::ostream& operator<<(std::ostream& os, minimizer_tuple const& mt) { return os; } +/* + The bucket of a minimizer: its tuples, sorted by pos_in_seq. The minimizer + is forward (see `util::compute_minimizer`), so a locus is never abandoned + and later re-selected: each tuple carries a distinct pos_in_seq, and the + number of tuples (= super-kmers) equals the number of minimizer positions, + which is what size() returns. `minimizers_tuples::merge` checks this. +*/ struct bucket_type { bucket_type(minimizer_tuple const* begin, minimizer_tuple const* end) : m_begin(begin) - , m_end(end) - , m_num_super_kmers(std::distance(begin, end)) - , m_num_minimizer_positions(0) // + , m_end(end) // { - uint64_t prev_pos_in_seq = constants::invalid_uint64; - while (begin != end) { - uint64_t pos_in_seq = (*begin).pos_in_seq; - if (pos_in_seq != prev_pos_in_seq) { - ++m_num_minimizer_positions; - prev_pos_in_seq = pos_in_seq; - } - ++begin; - } - assert(m_num_minimizer_positions <= m_num_super_kmers); + assert([&] { /* one tuple per minimizer position: the scheme is forward */ + for (auto it = begin; it + 1 < end; ++it) { + if (it->pos_in_seq == (it + 1)->pos_in_seq) return false; + } + return true; + }()); } struct iterator { @@ -92,20 +93,7 @@ struct bucket_type { iterator begin() const { return iterator(m_begin); } iterator end() const { return iterator(m_end); } - /* - A super-kmer is uniquely identified by the couple - (minimizer offset, position of minimizer in the first kmer of the super-kmer). - These two components, together, give the - starting position of a super-kmer in the sequence. - - The minimizer is forward (see `util::compute_minimizer`), so a locus is - never abandoned and later re-selected, and the number of super-kmers - equals the number of minimizer positions. The code below does not rely - on that and stays correct for a non-forward scheme too, where size() -- - the number of minimizer positions -- is < the number of super-kmers. - */ - uint64_t num_super_kmers() const { return m_num_super_kmers; } - uint64_t size() const { return m_num_minimizer_positions; } + uint64_t size() const { return std::distance(m_begin, m_end); } minimizer_tuple const* begin_ptr() const { return m_begin; } minimizer_tuple const* end_ptr() const { return m_end; } @@ -113,8 +101,6 @@ struct bucket_type { private: minimizer_tuple const* m_begin; minimizer_tuple const* m_end; - uint64_t m_num_super_kmers; - uint64_t m_num_minimizer_positions; }; /* @@ -160,7 +146,6 @@ struct minimizers_tuples { minimizers_tuples(build_configuration const& build_config) : m_num_minimizers(0) , m_num_minimizer_positions(0) - , m_num_super_kmers(0) , m_run_identifier(pthash::clock_type::now().time_since_epoch().count()) , m_build_config(build_config) // { @@ -218,17 +203,9 @@ struct minimizers_tuples { assert(m_num_minimizers == 0); assert(m_num_minimizer_positions == 0); - assert(m_num_super_kmers == 0); mm::file_source input(get_minimizers_filename(), mm::advice::sequential); - for (minimizers_tuples_iterator it(input.data(), input.data() + input.size()); - it.has_next(); it.next()) // - { - auto bucket = it.bucket(); - m_num_minimizers += 1; - m_num_minimizer_positions += bucket.size(); - m_num_super_kmers += bucket.num_super_kmers(); - } + count_and_check_forward(input.data(), input.data() + input.size()); input.close(); return; } @@ -246,23 +223,15 @@ struct minimizers_tuples { m_num_minimizers = 0; m_num_minimizer_positions = 0; - m_num_super_kmers = 0; uint64_t prev_minimizer = constants::invalid_uint64; uint64_t prev_pos_in_seq = constants::invalid_uint64; while (fm_iterator.has_next()) { minimizer_tuple mt = *fm_iterator; - if (mt.minimizer != prev_minimizer) { - prev_minimizer = mt.minimizer; - ++m_num_minimizers; - ++m_num_minimizer_positions; - } else { - if (mt.pos_in_seq != prev_pos_in_seq) ++m_num_minimizer_positions; - } + count_and_check_forward_one(mt, prev_minimizer, prev_pos_in_seq); out.write(reinterpret_cast(&mt), sizeof(minimizer_tuple)); - prev_pos_in_seq = mt.pos_in_seq; - ++m_num_super_kmers; - if (m_build_config.verbose and m_num_super_kmers % 100'000'000 == 0) { - std::cout << "processed " << m_num_super_kmers << " minimizer tuples" << std::endl; + if (m_build_config.verbose and m_num_minimizer_positions % 100'000'000 == 0) { + std::cout << "processed " << m_num_minimizer_positions << " minimizer tuples" + << std::endl; } fm_iterator.next(); } @@ -279,8 +248,11 @@ struct minimizers_tuples { uint64_t num_files_to_merge() const { return m_num_files_to_merge; } uint64_t num_minimizers() const { return m_num_minimizers; } + + /* One tuple per minimizer position, the scheme being forward (checked by + `merge`), so this is also the number of super-kmers and the number of + records in the merged tuples file. */ uint64_t num_minimizer_positions() const { return m_num_minimizer_positions; } - uint64_t num_super_kmers() const { return m_num_super_kmers; } void remove_tmp_file() { std::remove(get_minimizers_filename().c_str()); } @@ -288,10 +260,39 @@ struct minimizers_tuples { std::atomic m_num_files_to_merge; uint64_t m_num_minimizers; uint64_t m_num_minimizer_positions; - uint64_t m_num_super_kmers; uint64_t m_run_identifier; build_configuration m_build_config; + /* + Count minimizers and minimizer positions over tuples sorted by + (minimizer, pos_in_seq), checking the forwardness requirement: the + scheme never re-selects an abandoned locus, so no (minimizer, + pos_in_seq) pair may appear twice. Everything downstream sizes the + index by the tuple count, so a violation must stop the build. + */ + void count_and_check_forward_one(minimizer_tuple const& mt, uint64_t& prev_minimizer, + uint64_t& prev_pos_in_seq) // + { + if (mt.minimizer != prev_minimizer) { + prev_minimizer = mt.minimizer; + ++m_num_minimizers; + } else if (mt.pos_in_seq == prev_pos_in_seq) { + throw std::runtime_error( + "the minimizer scheme is not forward: " + "a (minimizer, position) pair occurs in more than one super-kmer"); + } + prev_pos_in_seq = mt.pos_in_seq; + ++m_num_minimizer_positions; + } + + void count_and_check_forward(minimizer_tuple const* begin, minimizer_tuple const* end) { + uint64_t prev_minimizer = constants::invalid_uint64; + uint64_t prev_pos_in_seq = constants::invalid_uint64; + for (; begin != end; ++begin) { + count_and_check_forward_one(*begin, prev_minimizer, prev_pos_in_seq); + } + } + std::string get_tmp_output_filename(uint64_t id) const { std::stringstream filename; filename << m_build_config.tmp_dirname << "/sshash.tmp.run_" << m_run_identifier diff --git a/src/builder/build_sparse_and_skew_index.cpp b/src/builder/build_sparse_and_skew_index.cpp index 724b354..2ef51dc 100644 --- a/src/builder/build_sparse_and_skew_index.cpp +++ b/src/builder/build_sparse_and_skew_index.cpp @@ -21,7 +21,6 @@ void dictionary_builder::build_sparse_and_skew_index( uint64_t num_buckets_larger_than_1_not_in_skew_index = 0; uint64_t num_buckets_in_skew_index = 0; - uint64_t num_super_kmers_in_buckets_larger_than_1 = 0; uint64_t num_minimizer_positions_of_buckets_larger_than_1 = 0; uint64_t num_minimizer_positions_of_buckets_in_skew_index = 0; @@ -41,7 +40,6 @@ void dictionary_builder::build_sparse_and_skew_index( ++num_buckets_in_skew_index; num_minimizer_positions_of_buckets_in_skew_index += bucket_size; } - num_super_kmers_in_buckets_larger_than_1 += bucket.num_super_kmers(); } for (auto mt : bucket) { @@ -97,7 +95,8 @@ void dictionary_builder::build_sparse_and_skew_index( std::vector buckets; buckets.reserve(num_buckets_larger_than_1_not_in_skew_index + num_buckets_in_skew_index); std::vector tuples; // backed memory - tuples.reserve(num_super_kmers_in_buckets_larger_than_1); + tuples.reserve(num_minimizer_positions_of_buckets_larger_than_1 + + num_minimizer_positions_of_buckets_in_skew_index); // Second pass: collect buckets > 1 for sorting AND handle size-1 buckets for (minimizers_tuples_iterator it(input.data(), input.data() + input.size()); // @@ -108,20 +107,13 @@ void dictionary_builder::build_sparse_and_skew_index( const uint64_t bucket_size = bucket.size(); if (bucket_size == 1) { - // Handle size-1 buckets: encode directly into control codewords - uint64_t prev_pos_in_seq = constants::invalid_uint64; - for (auto mt : bucket) { - if (mt.pos_in_seq != prev_pos_in_seq) { - /* - For minimizers occurring once, store a (log(N)+1)-bit - code, as follows: |offset|0|, i.e., the LSB is 0. - */ - uint64_t code = mt.pos_in_seq << 1; // first LS bit encodes status code: 0 - assert(code < (uint64_t(1) << num_bits_for_control)); - control_codewords_builder.set(bucket_id, code); - prev_pos_in_seq = mt.pos_in_seq; - } - } + /* + For minimizers occurring once, store a (log(N)+1)-bit + code, as follows: |offset|0|, i.e., the LSB is 0. + */ + const uint64_t code = (*bucket.begin()).pos_in_seq << 1; // LSB encodes status: 0 + assert(code < (uint64_t(1) << num_bits_for_control)); + control_codewords_builder.set(bucket_id, code); } else { // Collect buckets > 1 for later processing minimizer_tuple const* begin = tuples.data() + tuples.size(); @@ -202,37 +194,21 @@ void dictionary_builder::build_sparse_and_skew_index( } if (curr_bucket_size <= min_size) { - uint64_t prev_pos_in_seq = constants::invalid_uint64; - for (auto mt : bucket) { - if (prev_pos_in_seq == constants::invalid_uint64) { // only once - uint64_t p = (list_id << constants::min_l) | (curr_bucket_size - 2); - uint64_t code = (p << 2) | 1; // first two LS bits encode status code: 01 - assert(code < (uint64_t(1) << num_bits_for_control)); - control_codewords_builder.set(mt.minimizer, code); - } - if (mt.pos_in_seq != prev_pos_in_seq) { - mid_load_buckets_builder.push_back(mt.pos_in_seq); - prev_pos_in_seq = mt.pos_in_seq; - mid_load_buckets_size += 1; - } - } + uint64_t p = (list_id << constants::min_l) | (curr_bucket_size - 2); + uint64_t code = (p << 2) | 1; // first two LS bits encode status code: 01 + assert(code < (uint64_t(1) << num_bits_for_control)); + control_codewords_builder.set(bucket.begin_ptr()->minimizer, code); + for (auto mt : bucket) mid_load_buckets_builder.push_back(mt.pos_in_seq); + mid_load_buckets_size += bucket_size; ++list_id; } else { - uint64_t prev_pos_in_seq = constants::invalid_uint64; - for (auto mt : bucket) { - if (prev_pos_in_seq == constants::invalid_uint64) { // only once - assert(partition_id < 8); - uint64_t p = (heavy_load_buckets_size << 3) | partition_id; - uint64_t code = (p << 2) | 3; // first two LS bits encode status code: 11 - assert(code < (uint64_t(1) << num_bits_for_control)); - control_codewords_builder.set(mt.minimizer, code); - } - if (mt.pos_in_seq != prev_pos_in_seq) { - heavy_load_buckets_builder.push_back(mt.pos_in_seq); - prev_pos_in_seq = mt.pos_in_seq; - heavy_load_buckets_size += 1; - } - } + assert(partition_id < 8); + uint64_t p = (heavy_load_buckets_size << 3) | partition_id; + uint64_t code = (p << 2) | 3; // first two LS bits encode status code: 11 + assert(code < (uint64_t(1) << num_bits_for_control)); + control_codewords_builder.set(bucket.begin_ptr()->minimizer, code); + for (auto mt : bucket) heavy_load_buckets_builder.push_back(mt.pos_in_seq); + heavy_load_buckets_size += bucket_size; } } @@ -443,13 +419,9 @@ void dictionary_builder::build_sparse_and_skew_index( assert(bucket.size() > lower and bucket.size() <= upper); uint64_t pos_in_bucket = -1; - uint64_t prev_pos_in_seq = constants::invalid_uint64; for (auto mt : bucket) // { - if (mt.pos_in_seq != prev_pos_in_seq) { - prev_pos_in_seq = mt.pos_in_seq; - ++pos_in_bucket; - } + ++pos_in_bucket; // one position per tuple: the scheme is forward assert(mt.pos_in_seq >= mt.pos_in_kmer); mt.pos_in_seq = d.m_spss.strings_offsets.decode(mt.pos_in_seq).absolute_offset; diff --git a/src/builder/compute_minimizer_tuples.cpp b/src/builder/compute_minimizer_tuples.cpp index c9190ce..4070c57 100644 --- a/src/builder/compute_minimizer_tuples.cpp +++ b/src/builder/compute_minimizer_tuples.cpp @@ -28,14 +28,9 @@ void dictionary_builder::compute_minimizer_tuples() // uint64_t num_kmers_in_super_kmer) // { assert(num_kmers_in_super_kmer <= k - m + 1 /* max num kmers in super-kmer */); - if (!buffer.empty() and // - buffer.back().minimizer == mini_info.minimizer and // - buffer.back().pos_in_seq == mini_info.pos_in_seq and // - buffer.back().pos_in_kmer == mini_info.pos_in_kmer) // - { - buffer.back().num_kmers_in_super_kmer += num_kmers_in_super_kmer; - return; - } + /* the scheme is forward: a minimizer position is never saved twice */ + assert(buffer.empty() or !(buffer.back().minimizer == mini_info.minimizer and + buffer.back().pos_in_seq == mini_info.pos_in_seq)); if (buffer.size() == buffer_size) { minimizers.sort_and_flush(buffer); buffer.clear(); From a66c3ab6f3aa14f6dbd1d08f971bfabc9244413c Mon Sep 17 00:00:00 2001 From: Giulio Ermanno Pibiri Date: Thu, 20 Aug 2026 13:55:29 +0000 Subject: [PATCH 08/14] scripts: drop the regular/canonical distinction There is one indexing modality now, so the drivers no longer run each configuration twice. The `canonical` parameter, the `--canonical` flag (gone from the build tool) and the `.canon` index suffix are removed, and the per-k result files lose their mode prefix: `build.json`, `bench.json`, `streaming-queries.json` and the matching `.log` files. The trade-off plots lose the Mode dimension: two series per k instead of four, drawn in the saturated colours the canonical series used. `print_csv.py` loses its `canonical` column and reads the un-prefixed files. It is single-mode only: result directories archived under `benchmarks/results-*` still hold `regular-`/`canon-` pairs and are left untouched, so re-reading them needs the version of the script from the corresponding commit. The benchmarks README says so, and also gets its path to `print_csv.py` corrected. --- benchmarks/README.md | 11 ++- benchmarks/print_csv.py | 47 ++++-------- script/bench.py | 18 ++--- script/build.py | 25 +++---- script/plot-trade-off-l.py | 144 ++++++++++++++++++------------------- script/plot-trade-off-m.py | 142 ++++++++++++++++++------------------ script/streaming-query.py | 18 ++--- script/sweep-m.py | 54 +++++--------- script/sweep-min-l.py | 56 ++++++--------- 9 files changed, 221 insertions(+), 294 deletions(-) diff --git a/benchmarks/README.md b/benchmarks/README.md index 5e65e01..68e531a 100644 --- a/benchmarks/README.md +++ b/benchmarks/README.md @@ -25,5 +25,12 @@ Queries were run using one thread, instead. The results can be exported to CSV format with - python3 ../script/print_csv.py ../benchmarks/results-21-01-26/k31 - python3 ../script/print_csv.py ../benchmarks/results-21-01-26/k63 + python3 ../benchmarks/print_csv.py ../benchmarks/results-21-01-26/k31 + python3 ../benchmarks/print_csv.py ../benchmarks/results-21-01-26/k63 + +Note that the scripts now produce a single set of result files per `k` +(`build.json`, `bench.json`, and `streaming-queries.json`), since the +regular/canonical distinction is gone: there is one indexing modality only. +Result directories archived before this change instead hold a `regular-` and a +`canon-` file for each of those; to re-read them, use the version of +`print_csv.py` from the corresponding commit. diff --git a/benchmarks/print_csv.py b/benchmarks/print_csv.py index 6e60816..adaae7d 100644 --- a/benchmarks/print_csv.py +++ b/benchmarks/print_csv.py @@ -12,7 +12,7 @@ def format_time(microseconds): seconds = int(seconds % 60) return f"{minutes}:{seconds:02d}" -def parse_build_file(path, canonical_flag): +def parse_build_file(path): """Parse build JSONL file.""" results = [] with open(path) as f: @@ -42,14 +42,13 @@ def parse_build_file(path, canonical_flag): "k": k, "Collection": collection, "m": d["m"], - "canonical": "yes" if canonical_flag else "no", "bits_per_kmer": f"{bits_per_kmer:.2f}", "total_GB": f"{gb:.2f}", "build_time": build_time_fmt }) return results -def parse_bench_file(path, canonical_flag): +def parse_bench_file(path): """Parse benchmark JSONL file and average per collection.""" lookup_data = {} with open(path) as f: @@ -67,9 +66,8 @@ def parse_bench_file(path, canonical_flag): collection = fname.split(".")[0].capitalize() m = d["m"] k = d["k"] - canonical = "yes" if canonical_flag else "no" - key = (collection, m, canonical) + key = (collection, m) entry = lookup_data.setdefault(key, { "k": k, "pos": [], "neg": [], "access": [], "iter": [] @@ -94,7 +92,7 @@ def parse_bench_file(path, canonical_flag): return lookup_data -def parse_streaming_file(path, canonical_flag): +def parse_streaming_file(path): """Parse streaming queries JSON file.""" stream_data = {} if not os.path.exists(path): @@ -113,9 +111,7 @@ def parse_streaming_file(path, canonical_flag): fname = os.path.basename(d["index_filename"]) collection = fname.split(".")[0].capitalize() - canonical = "yes" if canonical_flag else "no" - key = (collection, canonical) num_kmers = int(d["num_kmers"]) num_pos = int(d["num_positive_kmers"]) num_ext = int(d["num_extensions"]) @@ -125,7 +121,7 @@ def parse_streaming_file(path, canonical_flag): hit_rate = (num_pos / num_kmers) * 100 if num_kmers else 0 extension_rate = (num_ext / num_pos) * 100 if num_pos else 0 - stream_data[key] = { + stream_data[collection] = { "ns_per_kmer": f"{ns_per_kmer}", "hit_rate": f"{hit_rate:.2f}", "extension_rate": f"{extension_rate:.2f}" @@ -135,41 +131,26 @@ def parse_streaming_file(path, canonical_flag): def main(): if len(sys.argv) != 2: - print("Usage: print.py input_dir", file=sys.stderr) + print("Usage: print_csv.py input_dir", file=sys.stderr) sys.exit(1) input_dir = sys.argv[1] - reg_build_path = input_dir + "/regular-build.json" - canon_build_path = input_dir + "/canon-build.json" - reg_bench_path = input_dir + "/regular-bench.json" - canon_bench_path = input_dir + "/canon-bench.json" - reg_stream_path = input_dir + "/regular-streaming-queries-high-hit.json" - canon_stream_path = input_dir + "/canon-streaming-queries-high-hit.json" - - reg_build = parse_build_file(reg_build_path, False) - canon_build = parse_build_file(canon_build_path, True) - reg_bench = parse_bench_file(reg_bench_path, False) - canon_bench = parse_bench_file(canon_bench_path, True) - reg_stream = parse_streaming_file(reg_stream_path, False) - canon_stream = parse_streaming_file(canon_stream_path, True) - - # merge everything - all_builds = reg_build + canon_build - lookup_all = {**reg_bench, **canon_bench} - stream_all = {**reg_stream, **canon_stream} + builds = parse_build_file(input_dir + "/build.json") + lookup_all = parse_bench_file(input_dir + "/bench.json") + stream_all = parse_streaming_file(input_dir + "/streaming-queries.json") # CSV header - print("k,Collection,m,canonical,bits_per_kmer,total_GB,build_time,positive_lookup_ns,negative_lookup_ns,access_ns,iteration_ns,ns_per_kmer,hit_rate,extension_rate") + print("k,Collection,m,bits_per_kmer,total_GB,build_time,positive_lookup_ns,negative_lookup_ns,access_ns,iteration_ns,ns_per_kmer,hit_rate,extension_rate") - for r in sorted(all_builds, key=lambda x: (int(x["k"]), x["Collection"], x["canonical"])): + for r in sorted(builds, key=lambda x: (int(x["k"]), x["Collection"])): lookup = lookup_all.get( - (r["Collection"], r["m"], r["canonical"]), # key + (r["Collection"], r["m"]), # key {"pos": "NA", "neg": "NA", "access": "NA", "iter": "NA", "k": r["k"]}) stream = stream_all.get( - (r["Collection"], r["canonical"]), # key + r["Collection"], # key {"ns_per_kmer": "NA", "hit_rate": "NA", "extension_rate": "NA"}) - print(f"{r['k']},{r['Collection']},{r['m']},{r['canonical']},{r['bits_per_kmer']},{r['total_GB']},{r['build_time']},{lookup['pos']},{lookup['neg']},{lookup['access']},{lookup['iter']},{stream['ns_per_kmer']},{stream['hit_rate']},{stream['extension_rate']}") + print(f"{r['k']},{r['Collection']},{r['m']},{r['bits_per_kmer']},{r['total_GB']},{r['build_time']},{lookup['pos']},{lookup['neg']},{lookup['access']},{lookup['iter']},{stream['ns_per_kmer']},{stream['hit_rate']},{stream['extension_rate']}") if __name__ == "__main__": main() diff --git a/script/bench.py b/script/bench.py index 2b334e3..f082b41 100644 --- a/script/bench.py +++ b/script/bench.py @@ -47,19 +47,17 @@ def build_project(max_k63: bool): run_cmd(["make", "-j"]) -def run_bench(k, canonical, runs = 3): +def run_bench(k, runs = 3): """Run SSHASH benchmark for all datasets.""" - mode = "canon" if canonical else "regular" out_dir = results_dir / f"k{k}" out_dir.mkdir(parents=True, exist_ok=True) - log_file = out_dir / f"{mode}-bench.log" - json_file = out_dir / f"{mode}-bench.json" + log_file = out_dir / "bench.log" + json_file = out_dir / "bench.json" for dataset in datasets: - suffix = f".k{k}.canon.sshash" if canonical else f".k{k}.sshash" - index_path = index_dir / f"{dataset}{suffix}" + index_path = index_dir / f"{dataset}.k{k}.sshash" - print(f"\n>>> Benchmarking {dataset} (k={k}, mode={mode})\n") + print(f"\n>>> Benchmarking {dataset} (k={k})\n") for i in range(runs): print(f" ==> run {i+1}/{runs}") cmd = ["./sshash", "bench", "-i", str(index_path)] @@ -82,13 +80,11 @@ def run_bench(k, canonical, runs = 3): # --- Build for k=31 --- build_project(max_k63=False) -run_bench(31, False) -run_bench(31, True) +run_bench(31) # --- Build for k=63 --- build_project(max_k63=True) -run_bench(63, False) -run_bench(63, True) +run_bench(63) # --- Restore to default --- build_project(max_k63=False) diff --git a/script/build.py b/script/build.py index c80973b..8964a79 100644 --- a/script/build.py +++ b/script/build.py @@ -53,23 +53,20 @@ def build_project(max_k63: bool): run_cmd(["make", "-j"]) -def build_sshash(k, canonical, m_values): - mode_dir = results_dir / f"k{k}" - mode_dir.mkdir(parents=True, exist_ok=True) +def build_sshash(k, m_values): + out_dir = results_dir / f"k{k}" + out_dir.mkdir(parents=True, exist_ok=True) - mode = "canon" if canonical else "regular" - log_file = mode_dir / f"{mode}-build.log" - json_file = mode_dir / f"{mode}-build.json" - time_file = mode_dir / f"{mode}-build.time.log" + log_file = out_dir / "build.log" + json_file = out_dir / "build.json" + time_file = out_dir / "build.time.log" for dataset in datasets: m_val = m_values[dataset] input_file = datasets_dir / f"{dataset}.k{k}.eulertigs.fa.gz" output_file = index_dir / f"{dataset}.k{k}" - if canonical: - output_file = str(output_file) + ".canon" - print(f"\n>>> Building {dataset} (k={k}, m={m_val}, mode={mode})\n") + print(f"\n>>> Building {dataset} (k={k}, m={m_val})\n") # Clean tmp directory (should be empty after each build anyway) subprocess.run(f"rm -rf {tmp_dir}/*", shell=True, check=True) @@ -86,8 +83,6 @@ def build_sshash(k, canonical, m_values): "-d", str(tmp_dir), "-o", f"{output_file}.sshash" ] - if canonical: - cmd.append("--canonical") # Append stdout to .log, stderr to .json with open(log_file, "a") as log, open(json_file, "a") as js: @@ -100,13 +95,11 @@ def build_sshash(k, canonical, m_values): # k = 31 build_project(max_k63=False) -build_sshash(31, False, m_values_k31) -build_sshash(31, True, m_values_k31) +build_sshash(31, m_values_k31) # k = 63 build_project(max_k63=True) -build_sshash(63, False, m_values_k63) -build_sshash(63, True, m_values_k63) +build_sshash(63, m_values_k63) # rebuild back to default build_project(max_k63=False) diff --git a/script/plot-trade-off-l.py b/script/plot-trade-off-l.py index 8b56ea2..424f0c6 100644 --- a/script/plot-trade-off-l.py +++ b/script/plot-trade-off-l.py @@ -32,102 +32,96 @@ def parse_results(results_dir): if not k_match: continue k_val = int(k_match.group(1)) - for mode in ['regular', 'canon']: - build_json = k_dir / f"{mode}-build.json" - bench_json = k_dir / f"{mode}-bench.json" - - if not build_json.exists() or not bench_json.exists(): - continue - - # 1. Parse Build JSON for Space (bits/k-mer) - space_dict = {} - with open(build_json, 'r') as f: - for line in f: - try: - j = json.loads(line) - - # Extract dataset name from "/mnt/.../human.k31.eulertigs.fa.gz" - filename = os.path.basename(j.get("input_filename", "")) - ds = filename.split('.')[0] if filename else "unknown" - - if "index_size_in_bytes" in j and "num_kmers" in j: - bytes_size = float(j["index_size_in_bytes"]) - num_kmers = float(j["num_kmers"]) - # Calculate bits per k-mer - bits_per_kmer = (bytes_size * 8.0) / num_kmers - space_dict[ds] = bits_per_kmer - except json.JSONDecodeError: - continue - - # 2. Parse Bench JSON for Query Time (ns/kmer) - time_dict = {} - count_dict = {} - with open(bench_json, 'r') as f: - for line in f: - try: - j = json.loads(line) - - # Extract dataset name from ".../human.k31.l4.sshash" - filename = os.path.basename(j.get("index_filename", "")) - ds = filename.split('.')[0] if filename else "unknown" - - # Use positive lookup time - t_str = j.get("positive lookup (avg_nanosec_per_kmer)") - if t_str is not None: - t = float(t_str) - time_dict[ds] = time_dict.get(ds, 0.0) + t - count_dict[ds] = count_dict.get(ds, 0) + 1 - except json.JSONDecodeError: - continue - - # 3. Combine and store - for ds in space_dict.keys(): - if ds in time_dict and count_dict[ds] > 0: - # Average the benchmark runs - avg_time = time_dict[ds] / count_dict[ds] - data.append({ - 'Dataset': ds, - 'k': k_val, - 'l': l_val, - 'Mode': mode, - 'Space (bits/k-mer)': space_dict[ds], - 'Query Time (ns/k-mer)': avg_time - }) + build_json = k_dir / "build.json" + bench_json = k_dir / "bench.json" + + if not build_json.exists() or not bench_json.exists(): + continue + + # 1. Parse Build JSON for Space (bits/k-mer) + space_dict = {} + with open(build_json, 'r') as f: + for line in f: + try: + j = json.loads(line) + + # Extract dataset name from "/mnt/.../human.k31.eulertigs.fa.gz" + filename = os.path.basename(j.get("input_filename", "")) + ds = filename.split('.')[0] if filename else "unknown" + + if "index_size_in_bytes" in j and "num_kmers" in j: + bytes_size = float(j["index_size_in_bytes"]) + num_kmers = float(j["num_kmers"]) + # Calculate bits per k-mer + bits_per_kmer = (bytes_size * 8.0) / num_kmers + space_dict[ds] = bits_per_kmer + except json.JSONDecodeError: + continue + + # 2. Parse Bench JSON for Query Time (ns/kmer) + time_dict = {} + count_dict = {} + with open(bench_json, 'r') as f: + for line in f: + try: + j = json.loads(line) + + # Extract dataset name from ".../human.k31.l4.sshash" + filename = os.path.basename(j.get("index_filename", "")) + ds = filename.split('.')[0] if filename else "unknown" + + # Use positive lookup time + t_str = j.get("positive lookup (avg_nanosec_per_kmer)") + if t_str is not None: + t = float(t_str) + time_dict[ds] = time_dict.get(ds, 0.0) + t + count_dict[ds] = count_dict.get(ds, 0) + 1 + except json.JSONDecodeError: + continue + + # 3. Combine and store + for ds in space_dict.keys(): + if ds in time_dict and count_dict[ds] > 0: + # Average the benchmark runs + avg_time = time_dict[ds] / count_dict[ds] + data.append({ + 'Dataset': ds, + 'k': k_val, + 'l': l_val, + 'Space (bits/k-mer)': space_dict[ds], + 'Query Time (ns/k-mer)': avg_time + }) return pd.DataFrame(data) def plot_tradeoff(df, output_img="tradeoff_plot_l.pdf"): """ Generates a space-time trade-off plot. - Different lines for datasets/k/modes, points vary by 'l'. + Different lines for datasets/k, points vary by 'l'. """ if df.empty: print("No data parsed! Please check the JSON keys in the script.") return # Enforce categorical order so the legend is populated exactly how we want: - # Human before SE, and regular before canon + # Human before SE df['Dataset'] = pd.Categorical(df['Dataset'], categories=['human', 'se'], ordered=True) - df['Mode'] = pd.Categorical(df['Mode'], categories=['regular', 'canon'], ordered=True) - df = df.sort_values(by=['Dataset', 'k', 'Mode']) + df = df.sort_values(by=['Dataset', 'k']) plt.figure(figsize=(5, 10)) plt.style.use('seaborn-v0_8-whitegrid') - # Group by Dataset, k, and Mode with sort=False to preserve our categorical ordering - groups = df.groupby(['Dataset', 'k', 'Mode'], sort=False) + # Group by Dataset and k with sort=False to preserve our categorical ordering + groups = df.groupby(['Dataset', 'k'], sort=False) - for (dataset, k, mode), group in groups: + for (dataset, k), group in groups: # Sort by l to make the line connect logically group = group.sort_values(by='l') - label = f"{'Human' if dataset == 'human' else 'SE'} (k={k}, {mode})" - - # Color logic: Red for Human, Blue for SE. Darker if canonical. - if dataset == 'human': - color = 'firebrick' if mode == 'canon' else 'lightcoral' - else: # se - color = 'royalblue' if mode == 'canon' else 'lightskyblue' + label = f"{'Human' if dataset == 'human' else 'SE'} (k={k})" + + # Color logic: Red for Human, Blue for SE. + color = 'firebrick' if dataset == 'human' else 'royalblue' # Marker logic: Circle for k=31, Square for k=63 marker = 'o' if k == 31 else 's' diff --git a/script/plot-trade-off-m.py b/script/plot-trade-off-m.py index 0a1d079..cfffbee 100644 --- a/script/plot-trade-off-m.py +++ b/script/plot-trade-off-m.py @@ -32,102 +32,96 @@ def parse_results(results_dir): if not k_match: continue k_val = int(k_match.group(1)) - for mode in ['regular', 'canon']: - build_json = k_dir / f"{mode}-build.json" - bench_json = k_dir / f"{mode}-bench.json" - - if not build_json.exists() or not bench_json.exists(): - continue - - # 1. Parse Build JSON for Space (bits/k-mer) - space_dict = {} - with open(build_json, 'r') as f: - for line in f: - try: - j = json.loads(line) - - # Extract dataset name from "/mnt/.../human.k31.eulertigs.fa.gz" - filename = os.path.basename(j.get("input_filename", "")) - ds = filename.split('.')[0] if filename else "unknown" - - if "index_size_in_bytes" in j and "num_kmers" in j: - bytes_size = float(j["index_size_in_bytes"]) - num_kmers = float(j["num_kmers"]) - # Calculate bits per k-mer - bits_per_kmer = (bytes_size * 8.0) / num_kmers - space_dict[ds] = bits_per_kmer - except json.JSONDecodeError: - continue - - # 2. Parse Bench JSON for Query Time (ns/kmer) - time_dict = {} - count_dict = {} - with open(bench_json, 'r') as f: - for line in f: - try: - j = json.loads(line) - - # Extract dataset name from ".../human.k31.m17.sshash" - filename = os.path.basename(j.get("index_filename", "")) - ds = filename.split('.')[0] if filename else "unknown" - - # Use positive lookup time - t_str = j.get("positive lookup (avg_nanosec_per_kmer)") - if t_str is not None: - t = float(t_str) - time_dict[ds] = time_dict.get(ds, 0.0) + t - count_dict[ds] = count_dict.get(ds, 0) + 1 - except json.JSONDecodeError: - continue - - # 3. Combine and store - for ds in space_dict.keys(): - if ds in time_dict and count_dict[ds] > 0: - # Average the 3 benchmark runs - avg_time = time_dict[ds] / count_dict[ds] - data.append({ - 'Dataset': ds, - 'k': k_val, - 'm': m_val, - 'Mode': mode, - 'Space (bits/k-mer)': space_dict[ds], - 'Query Time (ns/k-mer)': avg_time - }) + build_json = k_dir / "build.json" + bench_json = k_dir / "bench.json" + + if not build_json.exists() or not bench_json.exists(): + continue + + # 1. Parse Build JSON for Space (bits/k-mer) + space_dict = {} + with open(build_json, 'r') as f: + for line in f: + try: + j = json.loads(line) + + # Extract dataset name from "/mnt/.../human.k31.eulertigs.fa.gz" + filename = os.path.basename(j.get("input_filename", "")) + ds = filename.split('.')[0] if filename else "unknown" + + if "index_size_in_bytes" in j and "num_kmers" in j: + bytes_size = float(j["index_size_in_bytes"]) + num_kmers = float(j["num_kmers"]) + # Calculate bits per k-mer + bits_per_kmer = (bytes_size * 8.0) / num_kmers + space_dict[ds] = bits_per_kmer + except json.JSONDecodeError: + continue + + # 2. Parse Bench JSON for Query Time (ns/kmer) + time_dict = {} + count_dict = {} + with open(bench_json, 'r') as f: + for line in f: + try: + j = json.loads(line) + + # Extract dataset name from ".../human.k31.m17.sshash" + filename = os.path.basename(j.get("index_filename", "")) + ds = filename.split('.')[0] if filename else "unknown" + + # Use positive lookup time + t_str = j.get("positive lookup (avg_nanosec_per_kmer)") + if t_str is not None: + t = float(t_str) + time_dict[ds] = time_dict.get(ds, 0.0) + t + count_dict[ds] = count_dict.get(ds, 0) + 1 + except json.JSONDecodeError: + continue + + # 3. Combine and store + for ds in space_dict.keys(): + if ds in time_dict and count_dict[ds] > 0: + # Average the 3 benchmark runs + avg_time = time_dict[ds] / count_dict[ds] + data.append({ + 'Dataset': ds, + 'k': k_val, + 'm': m_val, + 'Space (bits/k-mer)': space_dict[ds], + 'Query Time (ns/k-mer)': avg_time + }) return pd.DataFrame(data) def plot_tradeoff(df, output_img="tradeoff_plot.pdf"): """ Generates a space-time trade-off plot. - Different lines for datasets/k/modes, points vary by 'm'. + Different lines for datasets/k, points vary by 'm'. """ if df.empty: print("No data parsed! Please check the JSON keys in the script.") return # Enforce categorical order so the legend is populated exactly how we want: - # Human before SE, and regular before canon + # Human before SE df['Dataset'] = pd.Categorical(df['Dataset'], categories=['human', 'se'], ordered=True) - df['Mode'] = pd.Categorical(df['Mode'], categories=['regular', 'canon'], ordered=True) - df = df.sort_values(by=['Dataset', 'k', 'Mode']) + df = df.sort_values(by=['Dataset', 'k']) plt.figure(figsize=(5, 10)) plt.style.use('seaborn-v0_8-whitegrid') - # Group by Dataset, k, and Mode with sort=False to preserve our categorical ordering - groups = df.groupby(['Dataset', 'k', 'Mode'], sort=False) + # Group by Dataset and k with sort=False to preserve our categorical ordering + groups = df.groupby(['Dataset', 'k'], sort=False) - for (dataset, k, mode), group in groups: + for (dataset, k), group in groups: # Sort by m to make the line connect logically group = group.sort_values(by='m') - label = f"{'Human' if dataset == 'human' else 'SE'} (k={k}, {mode})" + label = f"{'Human' if dataset == 'human' else 'SE'} (k={k})" - # Color logic: Red for Human, Blue for SE. Darker if canonical. - if dataset == 'human': - color = 'firebrick' if mode == 'canon' else 'lightcoral' - else: # se - color = 'royalblue' if mode == 'canon' else 'lightskyblue' + # Color logic: Red for Human, Blue for SE. + color = 'firebrick' if dataset == 'human' else 'royalblue' # Marker logic: Circle for k=31, Square for k=63 marker = 'o' if k == 31 else 's' diff --git a/script/streaming-query.py b/script/streaming-query.py index e62dc08..5776e3b 100644 --- a/script/streaming-query.py +++ b/script/streaming-query.py @@ -57,19 +57,17 @@ def build_project(max_k63: bool): run_cmd(["make", "-j"]) -def run_bench(k, canonical, runs = 1): +def run_bench(k, runs = 1): """Run SSHASH benchmark for all datasets.""" - mode = "canon" if canonical else "regular" out_dir = results_dir / f"k{k}" out_dir.mkdir(parents=True, exist_ok=True) - log_file = out_dir / f"{mode}-streaming-queries.log" - json_file = out_dir / f"{mode}-streaming-queries.json" + log_file = out_dir / "streaming-queries.log" + json_file = out_dir / "streaming-queries.json" for dataset in datasets: - suffix = f".k{k}.canon.sshash" if canonical else f".k{k}.sshash" - index_path = index_dir / f"{dataset}{suffix}" + index_path = index_dir / f"{dataset}.k{k}.sshash" - print(f"\n>>> Benchmarking {dataset} (k={k}, mode={mode})\n") + print(f"\n>>> Benchmarking {dataset} (k={k})\n") for i in range(runs): print(f" ==> run {i+1}/{runs}") cmd = ["./sshash", "query", "-i", str(index_path), "-q", str(query_dir) + "/" + queries[dataset] + ".fastq.gz"] @@ -92,13 +90,11 @@ def run_bench(k, canonical, runs = 1): # --- Build for k=31 --- build_project(max_k63=False) -run_bench(31, False) -run_bench(31, True) +run_bench(31) # --- Build for k=63 --- build_project(max_k63=True) -run_bench(63, False) -run_bench(63, True) +run_bench(63) # --- Restore to default --- build_project(max_k63=False) diff --git a/script/sweep-m.py b/script/sweep-m.py index 2110f92..8c411be 100644 --- a/script/sweep-m.py +++ b/script/sweep-m.py @@ -53,24 +53,21 @@ def build_project(max_k63: bool): ]) run_cmd(["make", "-j"]) -def build_sshash(k, canonical, dataset, m_val): +def build_sshash(k, dataset, m_val): # Differentiate results dir by m_val - mode_dir = results_dir / f"m{m_val}" / f"k{k}" - mode_dir.mkdir(parents=True, exist_ok=True) + out_dir = results_dir / f"m{m_val}" / f"k{k}" + out_dir.mkdir(parents=True, exist_ok=True) - mode = "canon" if canonical else "regular" - log_file = mode_dir / f"{mode}-build.log" - json_file = mode_dir / f"{mode}-build.json" - time_file = mode_dir / f"{mode}-build.time.log" + log_file = out_dir / "build.log" + json_file = out_dir / "build.json" + time_file = out_dir / "build.time.log" input_file = datasets_dir / f"{dataset}.k{k}.eulertigs.fa.gz" - + # Append m_val to the output filename output_file = index_dir / f"{dataset}.k{k}.m{m_val}" - if canonical: - output_file = str(output_file) + ".canon" - print(f"\n>>> Building {dataset} (k={k}, m={m_val}, mode={mode})\n") + print(f"\n>>> Building {dataset} (k={k}, m={m_val})\n") # Clean tmp directory (should be empty after each build anyway) subprocess.run(f"rm -rf {tmp_dir}/*", shell=True, check=True) @@ -87,28 +84,23 @@ def build_sshash(k, canonical, dataset, m_val): "-d", str(tmp_dir), "-o", f"{output_file}.sshash" ] - if canonical: - cmd.append("--canonical") # Append stdout to .log, stderr to .json with open(log_file, "a") as log, open(json_file, "a") as js: subprocess.run(cmd, stdout=log, stderr=js, check=True) -def run_bench(k, canonical, dataset, m_val, runs=3): +def run_bench(k, dataset, m_val, runs=3): """Run SSHASH benchmark for a specific dataset and m.""" - mode = "canon" if canonical else "regular" - # Store results in the specific m_val / k folder out_dir = results_dir / f"m{m_val}" / f"k{k}" out_dir.mkdir(parents=True, exist_ok=True) - log_file = out_dir / f"{mode}-bench.log" - json_file = out_dir / f"{mode}-bench.json" + log_file = out_dir / "bench.log" + json_file = out_dir / "bench.json" - # Match the new index naming scheme that includes m_val - suffix = f".k{k}.m{m_val}.canon.sshash" if canonical else f".k{k}.m{m_val}.sshash" - index_path = index_dir / f"{dataset}{suffix}" + # Match the index naming scheme that includes m_val + index_path = index_dir / f"{dataset}.k{k}.m{m_val}.sshash" - print(f"\n>>> Benchmarking {dataset} (k={k}, m={m_val}, mode={mode})\n") + print(f"\n>>> Benchmarking {dataset} (k={k}, m={m_val})\n") for i in range(runs): print(f" ==> run {i+1}/{runs}") cmd = ["./sshash", "bench", "-i", str(index_path)] @@ -131,26 +123,16 @@ def run_bench(k, canonical, dataset, m_val, runs=3): for dataset in datasets: for current_m in m_sweeps_k31[dataset]: - # Regular - build_sshash(31, False, dataset, current_m) - run_bench(31, False, dataset, current_m) - - # Canonical - build_sshash(31, True, dataset, current_m) - run_bench(31, True, dataset, current_m) + build_sshash(31, dataset, current_m) + run_bench(31, dataset, current_m) # --- k = 63 Sweep --- build_project(max_k63=True) for dataset in datasets: for current_m in m_sweeps_k63[dataset]: - # Regular - build_sshash(63, False, dataset, current_m) - run_bench(63, False, dataset, current_m) - - # Canonical - build_sshash(63, True, dataset, current_m) - run_bench(63, True, dataset, current_m) + build_sshash(63, dataset, current_m) + run_bench(63, dataset, current_m) # Restore default compilation at the end print("\nRestoring default compilation (max_k63=False)...") diff --git a/script/sweep-min-l.py b/script/sweep-min-l.py index f16729a..c2bf862 100644 --- a/script/sweep-min-l.py +++ b/script/sweep-min-l.py @@ -87,26 +87,23 @@ def build_project(max_k63: bool): ]) run_cmd(["make", "-j"]) -def build_sshash(k, canonical, m_values, l_val): +def build_sshash(k, m_values, l_val): # Differentiate results dir by l_val - mode_dir = results_dir / f"l{l_val}" / f"k{k}" - mode_dir.mkdir(parents=True, exist_ok=True) + out_dir = results_dir / f"l{l_val}" / f"k{k}" + out_dir.mkdir(parents=True, exist_ok=True) - mode = "canon" if canonical else "regular" - log_file = mode_dir / f"{mode}-build.log" - json_file = mode_dir / f"{mode}-build.json" - time_file = mode_dir / f"{mode}-build.time.log" + log_file = out_dir / "build.log" + json_file = out_dir / "build.json" + time_file = out_dir / "build.time.log" for dataset in datasets: m_val = m_values[dataset] input_file = datasets_dir / f"{dataset}.k{k}.eulertigs.fa.gz" - + # 2. Append l_val to the output filename output_file = index_dir / f"{dataset}.k{k}.l{l_val}" - if canonical: - output_file = str(output_file) + ".canon" - print(f"\n>>> Building {dataset} (k={k}, m={m_val}, l={l_val}, mode={mode})\n") + print(f"\n>>> Building {dataset} (k={k}, m={m_val}, l={l_val})\n") # Clean tmp directory (should be empty after each build anyway) subprocess.run(f"rm -rf {tmp_dir}/*", shell=True, check=True) @@ -123,29 +120,24 @@ def build_sshash(k, canonical, m_values, l_val): "-d", str(tmp_dir), "-o", f"{output_file}.sshash" ] - if canonical: - cmd.append("--canonical") # Append stdout to .log, stderr to .json with open(log_file, "a") as log, open(json_file, "a") as js: subprocess.run(cmd, stdout=log, stderr=js, check=True) -def run_bench(k, canonical, l_val, runs=3): +def run_bench(k, l_val, runs=3): """Run SSHASH benchmark for all datasets.""" - mode = "canon" if canonical else "regular" - # Store results in the specific l_val / k folder out_dir = results_dir / f"l{l_val}" / f"k{k}" out_dir.mkdir(parents=True, exist_ok=True) - log_file = out_dir / f"{mode}-bench.log" - json_file = out_dir / f"{mode}-bench.json" + log_file = out_dir / "bench.log" + json_file = out_dir / "bench.json" for dataset in datasets: - # Match the new index naming scheme that includes l_val - suffix = f".k{k}.l{l_val}.canon.sshash" if canonical else f".k{k}.l{l_val}.sshash" - index_path = index_dir / f"{dataset}{suffix}" + # Match the index naming scheme that includes l_val + index_path = index_dir / f"{dataset}.k{k}.l{l_val}.sshash" - print(f"\n>>> Benchmarking {dataset} (k={k}, l={l_val}, mode={mode})\n") + print(f"\n>>> Benchmarking {dataset} (k={k}, l={l_val})\n") for i in range(runs): print(f" ==> run {i+1}/{runs}") cmd = ["./sshash", "bench", "-i", str(index_path)] @@ -166,23 +158,15 @@ def run_bench(k, canonical, l_val, runs=3): # Update the header file update_constants_hpp(current_l) - # Build and benchmark for k = 31 (Regular) + # Build and benchmark for k = 31 build_project(max_k63=False) - build_sshash(31, False, m_values_k31, current_l) - run_bench(31, False, current_l) - - # Build and benchmark for k = 31 (Canonical) - build_sshash(31, True, m_values_k31, current_l) - run_bench(31, True, current_l) + build_sshash(31, m_values_k31, current_l) + run_bench(31, current_l) - # Build and benchmark for k = 63 (Regular) + # Build and benchmark for k = 63 build_project(max_k63=True) - build_sshash(63, False, m_values_k63, current_l) - run_bench(63, False, current_l) - - # Build and benchmark for k = 63 (Canonical) - build_sshash(63, True, m_values_k63, current_l) - run_bench(63, True, current_l) + build_sshash(63, m_values_k63, current_l) + run_bench(63, current_l) # Restore default constants file at the end print("\nRestoring default constants.hpp (min_l = 6)...") From 87cab5ff0f22df2d747b0f47bc3aeaf09b07f16e Mon Sep 17 00:00:00 2001 From: Giulio Ermanno Pibiri Date: Fri, 28 Aug 2026 13:39:35 +0000 Subject: [PATCH 09/14] fix reverse buffer refill for kmer types wider than 64 bits `fill_buff_reverse` appended the words of [m_pos - uint_kmer_bits, m_pos) from the lowest address up, but `append64` appends at the low end of the buffer, so for the 128-bit kmer type the two words ended up swapped and `get_reverse` returned the kmer 32 characters below the requested position. The 64-bit type fills a single word and was unaffected. The bug could not produce wrong query results: it made every backward extension of a streaming query mis-compare under SSHASH_USE_MAX_KMER_LENGTH_63, so the query fell back to a full seed per kmer -- correct, but paying one lookup per kmer along backward matches. Caught by the ground-truth assertion in `streaming_query::lookup` once the same-minimizer memo started trusting the extension comparison. --- include/kmer_iterator.hpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/include/kmer_iterator.hpp b/include/kmer_iterator.hpp index 8672ab1..356db0e 100644 --- a/include/kmer_iterator.hpp +++ b/include/kmer_iterator.hpp @@ -72,7 +72,10 @@ struct kmer_iterator // inline void fill_buff_reverse() { static_assert(Kmer::uint_kmer_bits % 64 == 0); - for (int i = Kmer::uint_kmer_bits; i > 0; i -= 64) { + /* fill the buffer with the bits at [m_pos - uint_kmer_bits, m_pos), + in order: `append64` appends at the low end, so words must be + appended from the highest address down to the lowest */ + for (int i = 64; i <= Kmer::uint_kmer_bits; i += 64) { m_buff.append64(m_bv->get_word64(std::max(m_pos, Kmer::uint_kmer_bits) - i)); } m_avail = std::min(m_pos, Kmer::uint_kmer_bits); From c22c8979c582ad4c4f3b107bebea84887bf21147 Mon Sep 17 00:00:00 2001 From: Giulio Ermanno Pibiri Date: Fri, 28 Aug 2026 14:11:37 +0000 Subject: [PATCH 10/14] streaming queries: same-minimizer memo for negative-in-positive streams A sequencing error inside a kmer makes the k kmers covering it absent, while their minimizer is usually still present: each such kmer paid a full dictionary lookup (codeword access, offset decoding, text comparisons). Two memos in `streaming_query::seed` now answer most of these from the state of the last seed. `spectrum_preserving_string_set::lookup` is split so that its verification half, `lookup_from_positions`, can run on its own, and it learns to export the decoded locate set of the minimizer into a small `bucket_cache`, by decoding straight into it at no extra cost (buckets of size <= 8; HEAVYLOAD buckets, whose size is unknown, are never cached). `dictionary::lookup` forwards an optional cache pointer, and the streaming query keeps one such cache, filled at every real seed. When a seeded kmer carries the same minimizer as the last seed: - minimizer absent (existing memo): negative, unchanged; - singleton bucket, same minimizer occurrence (`pos_in_seq` equal: the scheme is forward, so equal sampled positions mean the occurrence was never abandoned), anchored at a positive match: negative for free. The only locus the kmer could occupy is the one the just-failed extension step already rejected, or it lies beyond the string boundary; the mirror locus would host the reverse complement of the minimizer occurrence, so it is excluded whenever the minimizer is not its own reverse complement (self-rc minimizers, possible for even m only, fall through to the next memo); - cached bucket (<= 8 positions): verify the cached locate set directly with `lookup_from_positions` -- same code path as a full lookup, with no codeword access and no offset decoding. Positive results re-anchor the singleton memo and set up extension as usual. Heavy buckets fall through to the full lookup. Error-free positive streams never enter the memos and are unaffected. The per-kmer ground-truth assertion in `lookup` (debug builds) is unchanged and covers every memo answer; `num_searches`/`num_extensions` and all reported counts are identical to the unmemoized code, with two new counters (`num_memo_singleton`, `num_memo_light`) in the query report. `seed` is outlined (it inlines into `lookup` at two call sites, and its body, grown by the memos, measurably slowed down streams that rarely seed when left inline). On simulated 150bp reads at 1% substitution error over the ecoli1 and penicillium test datasets (k = 31, m = 15), the memos answer 38% of the seed calls and streaming lookup improves by 9-11%; error-free and very low-hit streams pay a few percent from the residual code restructuring. --- .gitignore | 2 +- include/dictionary.hpp | 6 +- include/spectrum_preserving_string_set.hpp | 52 +++++++-- include/streaming_query.hpp | 117 ++++++++++++++++++--- include/util.hpp | 10 +- src/dictionary.cpp | 9 +- src/query.cpp | 6 ++ tools/query.cpp | 4 + 8 files changed, 180 insertions(+), 26 deletions(-) diff --git a/.gitignore b/.gitignore index 3094469..51855af 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,2 @@ .DS_Store -build +build* diff --git a/include/dictionary.hpp b/include/dictionary.hpp index 4eea0a1..b5fb618 100644 --- a/include/dictionary.hpp +++ b/include/dictionary.hpp @@ -11,6 +11,7 @@ template struct dictionary // { using kmer_type = Kmer; + using spss_type = spectrum_preserving_string_set; template friend struct dictionary_builder; @@ -161,12 +162,13 @@ struct dictionary // uint16_t m_m; hasher_type m_hasher; - spectrum_preserving_string_set m_spss; + spss_type m_spss; sparse_and_skew_index m_ssi; weights m_weights; - lookup_result lookup(Kmer uint_kmer, Kmer uint_kmer_rc, minimizer_info mini_info) const; + lookup_result lookup(Kmer uint_kmer, Kmer uint_kmer_rc, minimizer_info mini_info, + typename spss_type::bucket_cache* cache = nullptr) const; void forward_neighbours(Kmer suffix, neighbourhood& res, bool check_reverse_complement) const; diff --git a/include/spectrum_preserving_string_set.hpp b/include/spectrum_preserving_string_set.hpp index 5b88aae..b5fd5eb 100644 --- a/include/spectrum_preserving_string_set.hpp +++ b/include/spectrum_preserving_string_set.hpp @@ -26,17 +26,42 @@ struct spectrum_preserving_string_set // return util::read_kmer_at(strings, k - 1, Kmer::bits_per_char * (string_end - k + 1)); } + /* + The decoded locate set of one minimizer, filled by `lookup` (which + decodes straight into it, at no extra cost) and consumed by + `lookup_from_positions`: it lets streaming queries verify kmers + sharing the minimizer of their last seed without a codeword lookup + nor offset decoding. Only buckets of known size up to + `max_cached_size` are cached: HEAVYLOAD buckets have unknown size, + and larger buckets are too rare to be worth a bigger cache. + Buckets whose minimizer is not found are never cached (size = 0). + */ + struct bucket_cache { + static constexpr uint64_t max_cached_size = 8; + uint64_t size = 0; // 0 = nothing cached + std::array positions; + }; + template - lookup_result lookup(Iterator it, // - const Kmer kmer, const Kmer kmer_rc, // - const minimizer_info mini_info) const // + lookup_result lookup(Iterator it, // + const Kmer kmer, const Kmer kmer_rc, // + const minimizer_info mini_info, // + bucket_cache* cache = nullptr) const // { const uint64_t size = it.size(); assert(size > 0); static thread_local // std::array - v; + tl; + + auto* v = tl.data(); + if (cache != nullptr) { + cache->size = 0; + if (it.bucket_type() != bucket_t::HEAVYLOAD and size <= bucket_cache::max_cached_size) { + v = cache->positions.data(); + } + } for (uint64_t i = 0; i != size; ++i, ++it) { uint64_t minimizer_offset = *it; @@ -67,11 +92,26 @@ struct spectrum_preserving_string_set // return lookup_result(it.bucket_type() != bucket_t::HEAVYLOAD ? false : true); } + if (v != tl.data()) cache->size = size; + + return lookup_from_positions(v, size, kmer, kmer_rc, mini_info); + } + + /* + The verification half of `lookup`: try the candidate loci of each + position of the minimizer's locate set. The positions must be the + decoded locate set of `mini_info.minimizer`, whose presence at + `positions[0]` has already been checked. + */ + lookup_result lookup_from_positions(typename Offsets::decoded_offset const* positions, // + const uint64_t size, // + const Kmer kmer, const Kmer kmer_rc, // + const minimizer_info mini_info) const // + { lookup_result res; for (uint64_t i = 0; i != size; ++i) { - if (_lookup(res, v[i], kmer, kmer_rc, mini_info)) return res; + if (_lookup(res, positions[i], kmer, kmer_rc, mini_info)) return res; } - return lookup_result(); } diff --git a/include/streaming_query.hpp b/include/streaming_query.hpp index c039578..e82b10a 100644 --- a/include/streaming_query.hpp +++ b/include/streaming_query.hpp @@ -23,15 +23,19 @@ struct streaming_query // , m_minimizer_it(dict->m_k, dict->m_m, dict->m_hasher) , m_curr_mini_info() - , m_prev_mini_info() , m_it(dict->m_spss.strings, m_k) , m_remaining_string_bases(0) + , m_last_seed_mini_info() + , m_last_seed_positive(false) + , m_num_searches(0) , m_num_extensions(0) , m_num_invalid(0) , m_num_negative(0) + , m_num_memo_singleton(0) + , m_num_memo_light(0) {} @@ -40,6 +44,9 @@ struct streaming_query // m_remaining_string_bases = 0; m_res = lookup_result(); m_minimizer_it.reset(); + m_last_seed_mini_info = minimizer_info(); + m_last_seed_positive = false; + m_bucket_cache.size = 0; } lookup_result lookup(char const* kmer) // @@ -88,7 +95,6 @@ struct streaming_query // } /* 4. update state */ - m_prev_mini_info = m_curr_mini_info; m_start = false; assert(equal_lookup_result(m_dict->lookup(kmer), m_res)); @@ -100,6 +106,8 @@ struct streaming_query // uint64_t num_positive_lookups() const { return num_searches() + num_extensions(); } uint64_t num_negative_lookups() const { return m_num_negative; } uint64_t num_invalid_lookups() const { return m_num_invalid; } + uint64_t num_memo_singleton() const { return m_num_memo_singleton; } + uint64_t num_memo_light() const { return m_num_memo_light; } private: Dict const* m_dict; @@ -114,33 +122,114 @@ struct streaming_query // /* minimizer state */ minimizer_iterator m_minimizer_it; - minimizer_info m_curr_mini_info, m_prev_mini_info; + minimizer_info m_curr_mini_info; /* string state */ kmer_iterator m_it; uint64_t m_remaining_string_bases; + /* + Last real seed state, for the memos of `seed`: the minimizer of the + last `m_dict->lookup` call (with, for the singleton memo, the stream + position of its occurrence and whether a positive match anchors it), + and the cached decoded locate set of that minimizer. + */ + minimizer_info m_last_seed_mini_info; + bool m_last_seed_positive; + /* performance counts */ uint64_t m_num_searches; uint64_t m_num_extensions; uint64_t m_num_invalid; uint64_t m_num_negative; + uint64_t m_num_memo_singleton; + uint64_t m_num_memo_light; + + /* large and cold: keep it last, away from the per-kmer members above */ + typename Dict::spss_type::bucket_cache m_bucket_cache; - void seed() // + /* Whether the minimizer of the last seed is its own reverse complement + (possible for even m only). Rarely needed, so computed on demand. */ + bool self_rc_minimizer() const { + if constexpr (kmer_t::has_reverse_complement) { + return kmer_t::reverse_complement_mmer(m_last_seed_mini_info.minimizer, m_m) == + m_last_seed_mini_info.minimizer; + } + return false; + } + + /* Outlined: `seed` inlines into `lookup` at two call sites, and letting + its body (grown by the memos) bloat the extension fast path measurably + slows down streams that rarely seed. */ + __attribute__((noinline)) void seed() // { m_remaining_string_bases = 0; - /* if minimizer does not change and previous minimizer was not found, - surely any kmer having the same minimizer cannot be found as well */ - if (m_curr_mini_info.minimizer == m_prev_mini_info.minimizer and // - m_res.minimizer_found == false) // + if (m_curr_mini_info.minimizer == m_last_seed_mini_info.minimizer) // { - assert(m_res.kmer_id == constants::invalid_uint64); - m_num_negative += 1; - return; + /* The minimizer was absent at the last seed: any kmer having the + same minimizer is surely absent as well. */ + if (m_res.minimizer_found == false) { + assert(m_res.kmer_id == constants::invalid_uint64); + m_num_negative += 1; + return; + } + + /* + The minimizer is present but the kmer may not be -- the most + common negative in a positive stream, e.g., a sequencing error + inside the kmer. Two memos avoid the full lookup. + + 1. The last seed matched via a singleton bucket, at position L + say, and the current kmer still carries the very same + minimizer occurrence (same `pos_in_seq`: the scheme is + forward, so a sampled position, once abandoned, is never + re-selected). Then the only locus the kmer could occupy is + L - pos_in_kmer, which the just-failed extension step + already compared and rejected (or it lies beyond the string + boundary): negative, with no text access at all. The mirror + locus L - (k - m - pos_in_kmer) hosts the reverse + complement of the minimizer occurrence, not the occurrence + itself, so it is excluded too -- unless the minimizer is + its own reverse complement (possible for even m only), in + which case we fall through to the verification below. + */ + if (m_last_seed_positive and m_bucket_cache.size == 1 and + m_curr_mini_info.pos_in_seq == m_last_seed_mini_info.pos_in_seq and + !self_rc_minimizer()) // + { + m_res = lookup_result(); + m_num_negative += 1; + m_num_memo_singleton += 1; + return; + } + + /* 2. The decoded locate set of the bucket is cached: verify the + candidate loci directly against the text, with no codeword + lookup nor offset decoding. */ + if (m_bucket_cache.size != 0) { + m_res = m_dict->m_spss.lookup_from_positions( // + m_bucket_cache.positions.data(), m_bucket_cache.size, m_kmer, m_kmer_rc, + m_curr_mini_info); + m_num_memo_light += 1; + if (m_res.kmer_id == constants::invalid_uint64) { + m_num_negative += 1; + return; + } + /* keep the singleton memo anchored to the latest match */ + m_last_seed_mini_info = m_curr_mini_info; + m_last_seed_positive = true; + m_num_searches += 1; + begin_extension(); + return; + } + + /* heavy bucket: fall through to the full lookup */ } - m_res = m_dict->lookup(m_kmer, m_kmer_rc, m_curr_mini_info); + m_res = m_dict->lookup(m_kmer, m_kmer_rc, m_curr_mini_info, &m_bucket_cache); + m_last_seed_mini_info = m_curr_mini_info; + m_last_seed_positive = m_res.kmer_id != constants::invalid_uint64; if (m_res.kmer_id == constants::invalid_uint64) { m_num_negative += 1; @@ -149,6 +238,10 @@ struct streaming_query // assert(m_res.minimizer_found == true); m_num_searches += 1; + begin_extension(); + } + + void begin_extension() { uint64_t kmer_offset = 2 * (m_res.kmer_id + m_res.string_id * (m_k - 1)); m_remaining_string_bases = (m_res.string_end - m_res.string_begin - m_k) - m_res.kmer_id_in_string; diff --git a/include/util.hpp b/include/util.hpp index 9d4b18c..ff836fd 100644 --- a/include/util.hpp +++ b/include/util.hpp @@ -25,7 +25,9 @@ struct streaming_query_report { , num_negative_kmers(0) , num_invalid_kmers(0) , num_searches(0) - , num_extensions(0) {} + , num_extensions(0) + , num_memo_singleton(0) + , num_memo_light(0) {} uint64_t num_kmers; uint64_t num_positive_kmers; @@ -33,6 +35,12 @@ struct streaming_query_report { uint64_t num_invalid_kmers; uint64_t num_searches; uint64_t num_extensions; + + /* seeds answered from the state of the last seed, saving a full lookup: + singleton-bucket negatives answered for free, and light-bucket seeds + verified from the cached locate set (see `streaming_query::seed`) */ + uint64_t num_memo_singleton; + uint64_t num_memo_light; }; struct lookup_result { diff --git a/src/dictionary.cpp b/src/dictionary.cpp index ae5fa5d..711d2ca 100644 --- a/src/dictionary.cpp +++ b/src/dictionary.cpp @@ -26,14 +26,15 @@ lookup_result dictionary::lookup(Kmer uint_kmer, return res; } template -lookup_result dictionary::lookup(const Kmer uint_kmer, // - const Kmer uint_kmer_rc, // - const minimizer_info mini_info) const // +lookup_result dictionary::lookup(const Kmer uint_kmer, // + const Kmer uint_kmer_rc, // + const minimizer_info mini_info, // + typename spss_type::bucket_cache* cache) const // { assert(minimizer_info(mini_info.minimizer, mini_info.pos_in_kmer) == util::compute_minimizer(uint_kmer, uint_kmer_rc, m_k, m_m, m_hasher)); auto it = m_ssi.lookup(uint_kmer, uint_kmer_rc, mini_info); - return m_spss.lookup(it, uint_kmer, uint_kmer_rc, mini_info); + return m_spss.lookup(it, uint_kmer, uint_kmer_rc, mini_info, cache); } template diff --git a/src/query.cpp b/src/query.cpp index 29f6b1e..1024218 100644 --- a/src/query.cpp +++ b/src/query.cpp @@ -38,6 +38,8 @@ streaming_query_report streaming_query_from_fasta_file_multiline(Dict const* dic } report.num_searches = query.num_searches(); report.num_extensions = query.num_extensions(); + report.num_memo_singleton = query.num_memo_singleton(); + report.num_memo_light = query.num_memo_light(); report.num_positive_kmers = query.num_positive_lookups(); report.num_negative_kmers = query.num_negative_lookups(); report.num_invalid_kmers = query.num_invalid_lookups(); @@ -67,6 +69,8 @@ streaming_query_report streaming_query_from_fasta_file(Dict const* dict, std::is } report.num_searches = query.num_searches(); report.num_extensions = query.num_extensions(); + report.num_memo_singleton = query.num_memo_singleton(); + report.num_memo_light = query.num_memo_light(); report.num_positive_kmers = query.num_positive_lookups(); report.num_negative_kmers = query.num_negative_lookups(); report.num_invalid_kmers = query.num_invalid_lookups(); @@ -99,6 +103,8 @@ streaming_query_report streaming_query_from_fastq_file(Dict const* dict, std::is } report.num_searches = query.num_searches(); report.num_extensions = query.num_extensions(); + report.num_memo_singleton = query.num_memo_singleton(); + report.num_memo_light = query.num_memo_light(); report.num_positive_kmers = query.num_positive_lookups(); report.num_negative_kmers = query.num_negative_lookups(); report.num_invalid_kmers = query.num_invalid_lookups(); diff --git a/tools/query.cpp b/tools/query.cpp index 75bdcdd..ad6d98e 100644 --- a/tools/query.cpp +++ b/tools/query.cpp @@ -44,6 +44,8 @@ int query(int argc, char** argv) { query_stats.add("num_invalid_kmers", report.num_invalid_kmers); query_stats.add("num_searches", report.num_searches); query_stats.add("num_extensions", report.num_extensions); + query_stats.add("num_memo_singleton", report.num_memo_singleton); + query_stats.add("num_memo_light", report.num_memo_light); query_stats.add("elapsed_millisec", uint64_t(t.elapsed())); std::cout << "==== query report:\n"; @@ -60,6 +62,8 @@ int query(int argc, char** argv) { std::cout << "num_extensions = " << report.num_extensions << "/" << report.num_positive_kmers << " (" << (report.num_extensions * 100.0) / report.num_positive_kmers << "%)" << std::endl; + std::cout << "num_memo_singleton = " << report.num_memo_singleton << std::endl; + std::cout << "num_memo_light = " << report.num_memo_light << std::endl; std::cout << "elapsed = " << t.elapsed() / 1000 << " sec / "; std::cout << t.elapsed() / 1000 / 60 << " min / "; std::cout << (t.elapsed() * 1e6) / report.num_kmers << " ns/kmer" << std::endl; From 3fffd1afcda2eff3425770dcf7418f961c46e025 Mon Sep 17 00:00:00 2001 From: Giulio Ermanno Pibiri Date: Sat, 29 Aug 2026 07:45:33 +0000 Subject: [PATCH 11/14] rename the seed shortcut counters `num_memo_singleton` -> `num_skipped_singleton_lookups`: negatives implied by a singleton bucket, where the lookup is skipped entirely. Only a singleton bucket admits the free skip -- a larger bucket has candidate loci the failed extension never examined, which must be verified against the text. `num_memo_light` -> `num_bucket_cache_hits`: seeds resolved against the cached locate set of the last seed's minimizer, paying the text verification but not the codeword access nor the offset decoding. Same counters, clearer names; the "memo" wording is gone from comments too. No behavioural change. --- include/streaming_query.hpp | 29 ++++++++++++++++------------- include/util.hpp | 10 +++++----- src/query.cpp | 12 ++++++------ tools/query.cpp | 9 +++++---- 4 files changed, 32 insertions(+), 28 deletions(-) diff --git a/include/streaming_query.hpp b/include/streaming_query.hpp index e82b10a..d9c4a8d 100644 --- a/include/streaming_query.hpp +++ b/include/streaming_query.hpp @@ -34,8 +34,8 @@ struct streaming_query // , m_num_extensions(0) , m_num_invalid(0) , m_num_negative(0) - , m_num_memo_singleton(0) - , m_num_memo_light(0) + , m_num_skipped_singleton_lookups(0) + , m_num_bucket_cache_hits(0) {} @@ -106,8 +106,8 @@ struct streaming_query // uint64_t num_positive_lookups() const { return num_searches() + num_extensions(); } uint64_t num_negative_lookups() const { return m_num_negative; } uint64_t num_invalid_lookups() const { return m_num_invalid; } - uint64_t num_memo_singleton() const { return m_num_memo_singleton; } - uint64_t num_memo_light() const { return m_num_memo_light; } + uint64_t num_skipped_singleton_lookups() const { return m_num_skipped_singleton_lookups; } + uint64_t num_bucket_cache_hits() const { return m_num_bucket_cache_hits; } private: Dict const* m_dict; @@ -129,8 +129,8 @@ struct streaming_query // uint64_t m_remaining_string_bases; /* - Last real seed state, for the memos of `seed`: the minimizer of the - last `m_dict->lookup` call (with, for the singleton memo, the stream + Last real seed state, for the shortcuts of `seed`: the minimizer of + the last `m_dict->lookup` call (with, for the singleton shortcut, the stream position of its occurrence and whether a positive match anchors it), and the cached decoded locate set of that minimizer. */ @@ -142,8 +142,8 @@ struct streaming_query // uint64_t m_num_extensions; uint64_t m_num_invalid; uint64_t m_num_negative; - uint64_t m_num_memo_singleton; - uint64_t m_num_memo_light; + uint64_t m_num_skipped_singleton_lookups; + uint64_t m_num_bucket_cache_hits; /* large and cold: keep it last, away from the per-kmer members above */ typename Dict::spss_type::bucket_cache m_bucket_cache; @@ -159,7 +159,7 @@ struct streaming_query // } /* Outlined: `seed` inlines into `lookup` at two call sites, and letting - its body (grown by the memos) bloat the extension fast path measurably + its body (grown by the shortcuts) bloat the extension fast path measurably slows down streams that rarely seed. */ __attribute__((noinline)) void seed() // { @@ -178,7 +178,7 @@ struct streaming_query // /* The minimizer is present but the kmer may not be -- the most common negative in a positive stream, e.g., a sequencing error - inside the kmer. Two memos avoid the full lookup. + inside the kmer. Two shortcuts avoid the full lookup. 1. The last seed matched via a singleton bucket, at position L say, and the current kmer still carries the very same @@ -193,6 +193,9 @@ struct streaming_query // itself, so it is excluded too -- unless the minimizer is its own reverse complement (possible for even m only), in which case we fall through to the verification below. + Only a singleton bucket admits this free skip: a larger + bucket has candidate loci the failed extension never + examined, so those must be verified against the text. */ if (m_last_seed_positive and m_bucket_cache.size == 1 and m_curr_mini_info.pos_in_seq == m_last_seed_mini_info.pos_in_seq and @@ -200,7 +203,7 @@ struct streaming_query // { m_res = lookup_result(); m_num_negative += 1; - m_num_memo_singleton += 1; + m_num_skipped_singleton_lookups += 1; return; } @@ -211,12 +214,12 @@ struct streaming_query // m_res = m_dict->m_spss.lookup_from_positions( // m_bucket_cache.positions.data(), m_bucket_cache.size, m_kmer, m_kmer_rc, m_curr_mini_info); - m_num_memo_light += 1; + m_num_bucket_cache_hits += 1; if (m_res.kmer_id == constants::invalid_uint64) { m_num_negative += 1; return; } - /* keep the singleton memo anchored to the latest match */ + /* keep the singleton shortcut anchored to the latest match */ m_last_seed_mini_info = m_curr_mini_info; m_last_seed_positive = true; m_num_searches += 1; diff --git a/include/util.hpp b/include/util.hpp index ff836fd..3d7ee87 100644 --- a/include/util.hpp +++ b/include/util.hpp @@ -26,8 +26,8 @@ struct streaming_query_report { , num_invalid_kmers(0) , num_searches(0) , num_extensions(0) - , num_memo_singleton(0) - , num_memo_light(0) {} + , num_skipped_singleton_lookups(0) + , num_bucket_cache_hits(0) {} uint64_t num_kmers; uint64_t num_positive_kmers; @@ -37,10 +37,10 @@ struct streaming_query_report { uint64_t num_extensions; /* seeds answered from the state of the last seed, saving a full lookup: - singleton-bucket negatives answered for free, and light-bucket seeds + negatives implied by a singleton bucket (skipped entirely), and seeds verified from the cached locate set (see `streaming_query::seed`) */ - uint64_t num_memo_singleton; - uint64_t num_memo_light; + uint64_t num_skipped_singleton_lookups; + uint64_t num_bucket_cache_hits; }; struct lookup_result { diff --git a/src/query.cpp b/src/query.cpp index 1024218..3d96517 100644 --- a/src/query.cpp +++ b/src/query.cpp @@ -38,8 +38,8 @@ streaming_query_report streaming_query_from_fasta_file_multiline(Dict const* dic } report.num_searches = query.num_searches(); report.num_extensions = query.num_extensions(); - report.num_memo_singleton = query.num_memo_singleton(); - report.num_memo_light = query.num_memo_light(); + report.num_skipped_singleton_lookups = query.num_skipped_singleton_lookups(); + report.num_bucket_cache_hits = query.num_bucket_cache_hits(); report.num_positive_kmers = query.num_positive_lookups(); report.num_negative_kmers = query.num_negative_lookups(); report.num_invalid_kmers = query.num_invalid_lookups(); @@ -69,8 +69,8 @@ streaming_query_report streaming_query_from_fasta_file(Dict const* dict, std::is } report.num_searches = query.num_searches(); report.num_extensions = query.num_extensions(); - report.num_memo_singleton = query.num_memo_singleton(); - report.num_memo_light = query.num_memo_light(); + report.num_skipped_singleton_lookups = query.num_skipped_singleton_lookups(); + report.num_bucket_cache_hits = query.num_bucket_cache_hits(); report.num_positive_kmers = query.num_positive_lookups(); report.num_negative_kmers = query.num_negative_lookups(); report.num_invalid_kmers = query.num_invalid_lookups(); @@ -103,8 +103,8 @@ streaming_query_report streaming_query_from_fastq_file(Dict const* dict, std::is } report.num_searches = query.num_searches(); report.num_extensions = query.num_extensions(); - report.num_memo_singleton = query.num_memo_singleton(); - report.num_memo_light = query.num_memo_light(); + report.num_skipped_singleton_lookups = query.num_skipped_singleton_lookups(); + report.num_bucket_cache_hits = query.num_bucket_cache_hits(); report.num_positive_kmers = query.num_positive_lookups(); report.num_negative_kmers = query.num_negative_lookups(); report.num_invalid_kmers = query.num_invalid_lookups(); diff --git a/tools/query.cpp b/tools/query.cpp index ad6d98e..6793dcf 100644 --- a/tools/query.cpp +++ b/tools/query.cpp @@ -44,8 +44,8 @@ int query(int argc, char** argv) { query_stats.add("num_invalid_kmers", report.num_invalid_kmers); query_stats.add("num_searches", report.num_searches); query_stats.add("num_extensions", report.num_extensions); - query_stats.add("num_memo_singleton", report.num_memo_singleton); - query_stats.add("num_memo_light", report.num_memo_light); + query_stats.add("num_skipped_singleton_lookups", report.num_skipped_singleton_lookups); + query_stats.add("num_bucket_cache_hits", report.num_bucket_cache_hits); query_stats.add("elapsed_millisec", uint64_t(t.elapsed())); std::cout << "==== query report:\n"; @@ -62,8 +62,9 @@ int query(int argc, char** argv) { std::cout << "num_extensions = " << report.num_extensions << "/" << report.num_positive_kmers << " (" << (report.num_extensions * 100.0) / report.num_positive_kmers << "%)" << std::endl; - std::cout << "num_memo_singleton = " << report.num_memo_singleton << std::endl; - std::cout << "num_memo_light = " << report.num_memo_light << std::endl; + std::cout << "num_skipped_singleton_lookups = " << report.num_skipped_singleton_lookups + << std::endl; + std::cout << "num_bucket_cache_hits = " << report.num_bucket_cache_hits << std::endl; std::cout << "elapsed = " << t.elapsed() / 1000 << " sec / "; std::cout << t.elapsed() / 1000 / 60 << " min / "; std::cout << (t.elapsed() * 1e6) / report.num_kmers << " ns/kmer" << std::endl; From 4e2248c81094174c6be618c93de0fce434fe9514 Mon Sep 17 00:00:00 2001 From: Giulio Ermanno Pibiri Date: Sat, 29 Aug 2026 07:58:02 +0000 Subject: [PATCH 12/14] collapse the unreleased version bumps: current index version is 6.0.0 The branch bumped the index major version twice (5 -> 6 for the single canonical minimizer scheme, 6 -> 7 for the centre-closest tie-break), but neither was ever released: from the release lineage's viewpoint there is one breaking change since 5.1.1, so the next version is 6. Indexes stamped with the intermediate development versions need rebuilding, which `check_version_number` reports. --- include/constants.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/constants.hpp b/include/constants.hpp index b5f214d..215901a 100644 --- a/include/constants.hpp +++ b/include/constants.hpp @@ -20,7 +20,7 @@ constexpr int forward_orientation = 1; constexpr int backward_orientation = -1; namespace current_version_number { -constexpr uint8_t x = 7; +constexpr uint8_t x = 6; constexpr uint8_t y = 0; constexpr uint8_t z = 0; } // namespace current_version_number From 93947da02840f4105294ee740ba3e682bd91c9de Mon Sep 17 00:00:00 2001 From: Giulio Ermanno Pibiri Date: Sat, 29 Aug 2026 08:23:03 +0000 Subject: [PATCH 13/14] print_table.py: one table for both k, with a markdown option `print_csv.py` becomes `print_table.py`. It now processes a whole results directory at once, expecting the `k31` and `k63` subdirectories written by the benchmark scripts, and emits one table with the datasets in fixed order (Cod, Kestrel, Human, NCBI-v, SE, HPRC). The default output is CSV with every collected quantity; `--md` prints a markdown table with the main ones (space, building time, random lookup/access, streaming lookup), one block per k. The CSV header now labels the lookup and access columns `_us`: their values were always microseconds, despite the old `_ns` suffix. --- benchmarks/README.md | 15 ++- benchmarks/print_csv.py | 156 --------------------------- benchmarks/print_table.py | 214 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 224 insertions(+), 161 deletions(-) delete mode 100644 benchmarks/print_csv.py create mode 100644 benchmarks/print_table.py diff --git a/benchmarks/README.md b/benchmarks/README.md index 68e531a..9aa2b6f 100644 --- a/benchmarks/README.md +++ b/benchmarks/README.md @@ -23,14 +23,19 @@ Queries were run using one thread, instead. ![](results-21-01-26/results.png) -The results can be exported to CSV format with +The results of a run, for both k=31 and k=63 at once, can be exported to CSV +format (every collected quantity) or to a markdown table (the main quantities) +with - python3 ../benchmarks/print_csv.py ../benchmarks/results-21-01-26/k31 - python3 ../benchmarks/print_csv.py ../benchmarks/results-21-01-26/k63 + python3 ../benchmarks/print_table.py + python3 ../benchmarks/print_table.py --md + +where `` is the directory written by the scripts above (e.g., +`results-21-01-26`), containing the `k31` and `k63` subdirectories. Note that the scripts now produce a single set of result files per `k` (`build.json`, `bench.json`, and `streaming-queries.json`), since the regular/canonical distinction is gone: there is one indexing modality only. Result directories archived before this change instead hold a `regular-` and a -`canon-` file for each of those; to re-read them, use the version of -`print_csv.py` from the corresponding commit. +`canon-` file for each of those; to re-read them, use the `print_csv.py` +script from the corresponding commit. diff --git a/benchmarks/print_csv.py b/benchmarks/print_csv.py deleted file mode 100644 index adaae7d..0000000 --- a/benchmarks/print_csv.py +++ /dev/null @@ -1,156 +0,0 @@ -#!/usr/bin/env python3 - -import sys -import json -import os -from statistics import mean, StatisticsError -import math - -def format_time(microseconds): - seconds = microseconds / 1_000_000 - minutes = int(seconds // 60) - seconds = int(seconds % 60) - return f"{minutes}:{seconds:02d}" - -def parse_build_file(path): - """Parse build JSONL file.""" - results = [] - with open(path) as f: - for line in f: - line = line.strip() - if not line: - continue - try: - d = json.loads(line) - except json.JSONDecodeError: - print(f"Skipping invalid JSON line in {path}", file=sys.stderr) - continue - - num_kmers = int(d["num_kmers"]) - index_bytes = int(d["index_size_in_bytes"]) - build_time_us = int(d["total_build_time_in_microsec"]) - - bits_per_kmer = (index_bytes * 8) / num_kmers - gb = index_bytes / 1e9 - build_time_fmt = format_time(build_time_us) - - fname = os.path.basename(d["input_filename"]) - collection = fname.split(".")[0].capitalize() - k = d["k"] - - results.append({ - "k": k, - "Collection": collection, - "m": d["m"], - "bits_per_kmer": f"{bits_per_kmer:.2f}", - "total_GB": f"{gb:.2f}", - "build_time": build_time_fmt - }) - return results - -def parse_bench_file(path): - """Parse benchmark JSONL file and average per collection.""" - lookup_data = {} - with open(path) as f: - for line in f: - line = line.strip() - if not line: - continue - try: - d = json.loads(line) - except json.JSONDecodeError: - print(f"Skipping invalid JSON line in {path}", file=sys.stderr) - continue - - fname = os.path.basename(d["index_filename"]) - collection = fname.split(".")[0].capitalize() - m = d["m"] - k = d["k"] - - key = (collection, m) - entry = lookup_data.setdefault(key, { - "k": k, - "pos": [], "neg": [], "access": [], "iter": [] - }) - entry["pos"].append(float(d["positive lookup (avg_nanosec_per_kmer)"])) - entry["neg"].append(float(d["negative lookup (avg_nanosec_per_kmer)"])) - entry["access"].append(float(d["access (avg_nanosec_per_kmer)"])) - entry["iter"].append(float(d["iterator (avg_nanosec_per_kmer)"])) - - # average the results - for k, v in lookup_data.items(): - try: - lookup_data[k] = { - "k": v["k"], - "pos": f"{mean(v['pos'])/1000:.2f}", - "neg": f"{mean(v['neg'])/1000:.2f}", - "access": f"{mean(v['access'])/1000:.2f}", - "iter": f"{mean(v['iter']):.2f}", - } - except StatisticsError: - lookup_data[k] = {"k": v["k"], "pos": "NA", "neg": "NA", "access": "NA", "iter": "NA"} - return lookup_data - - -def parse_streaming_file(path): - """Parse streaming queries JSON file.""" - stream_data = {} - if not os.path.exists(path): - return stream_data - - with open(path) as f: - for line in f: - line = line.strip() - if not line: - continue - try: - d = json.loads(line) - except json.JSONDecodeError: - print(f"Skipping invalid JSON line in {path}", file=sys.stderr) - continue - - fname = os.path.basename(d["index_filename"]) - collection = fname.split(".")[0].capitalize() - - num_kmers = int(d["num_kmers"]) - num_pos = int(d["num_positive_kmers"]) - num_ext = int(d["num_extensions"]) - elapsed_ms = int(d["elapsed_millisec"]) - - ns_per_kmer = int(math.ceil(elapsed_ms * 1e6 / num_kmers)) - hit_rate = (num_pos / num_kmers) * 100 if num_kmers else 0 - extension_rate = (num_ext / num_pos) * 100 if num_pos else 0 - - stream_data[collection] = { - "ns_per_kmer": f"{ns_per_kmer}", - "hit_rate": f"{hit_rate:.2f}", - "extension_rate": f"{extension_rate:.2f}" - } - return stream_data - - -def main(): - if len(sys.argv) != 2: - print("Usage: print_csv.py input_dir", file=sys.stderr) - sys.exit(1) - - input_dir = sys.argv[1] - builds = parse_build_file(input_dir + "/build.json") - lookup_all = parse_bench_file(input_dir + "/bench.json") - stream_all = parse_streaming_file(input_dir + "/streaming-queries.json") - - # CSV header - print("k,Collection,m,bits_per_kmer,total_GB,build_time,positive_lookup_ns,negative_lookup_ns,access_ns,iteration_ns,ns_per_kmer,hit_rate,extension_rate") - - for r in sorted(builds, key=lambda x: (int(x["k"]), x["Collection"])): - lookup = lookup_all.get( - (r["Collection"], r["m"]), # key - {"pos": "NA", "neg": "NA", "access": "NA", "iter": "NA", "k": r["k"]}) - stream = stream_all.get( - r["Collection"], # key - {"ns_per_kmer": "NA", "hit_rate": "NA", "extension_rate": "NA"}) - - print(f"{r['k']},{r['Collection']},{r['m']},{r['bits_per_kmer']},{r['total_GB']},{r['build_time']},{lookup['pos']},{lookup['neg']},{lookup['access']},{lookup['iter']},{stream['ns_per_kmer']},{stream['hit_rate']},{stream['extension_rate']}") - -if __name__ == "__main__": - main() diff --git a/benchmarks/print_table.py b/benchmarks/print_table.py new file mode 100644 index 0000000..b79fda8 --- /dev/null +++ b/benchmarks/print_table.py @@ -0,0 +1,214 @@ +#!/usr/bin/env python3 + +""" +Print the benchmark results of both k=31 and k=63 as one table. + +The input directory must contain the two subdirectories `k31` and `k63`, +each with the files written by the benchmark scripts: +`build.json`, `bench.json`, and `streaming-queries.json`. + +The default output is CSV, with every collected quantity; with `--md`, +a markdown table with the main quantities is printed instead. +""" + +import argparse +import json +import os +import sys +import math +from statistics import mean, StatisticsError + +K_VALUES = [31, 63] + +# fixed dataset order and display names +DATASETS = ["cod", "kestrel", "human", "ncbi-virus", "se", "hprc"] +DISPLAY_NAME = { + "cod": "Cod", + "kestrel": "Kestrel", + "human": "Human", + "ncbi-virus": "NCBI-v", + "se": "SE", + "hprc": "HPRC", +} + +def format_time(microseconds): + seconds = microseconds / 1_000_000 + minutes = int(seconds // 60) + seconds = int(seconds % 60) + return f"{minutes}:{seconds:02d}" + +def dataset_of(filename): + return os.path.basename(filename).split(".")[0] + +def parse_build_file(path): + """Parse build JSONL file: one record per dataset.""" + results = {} + with open(path) as f: + for line in f: + line = line.strip() + if not line: + continue + try: + d = json.loads(line) + except json.JSONDecodeError: + print(f"Skipping invalid JSON line in {path}", file=sys.stderr) + continue + + num_kmers = int(d["num_kmers"]) + index_bytes = int(d["index_size_in_bytes"]) + build_time_us = int(d["total_build_time_in_microsec"]) + + results[dataset_of(d["input_filename"])] = { + "m": d["m"], + "bits_per_kmer": f"{(index_bytes * 8) / num_kmers:.2f}", + "total_GB": f"{index_bytes / 1e9:.2f}", + "build_time": format_time(build_time_us), + } + return results + +def parse_bench_file(path): + """Parse benchmark JSONL file and average the runs per dataset.""" + runs = {} + with open(path) as f: + for line in f: + line = line.strip() + if not line: + continue + try: + d = json.loads(line) + except json.JSONDecodeError: + print(f"Skipping invalid JSON line in {path}", file=sys.stderr) + continue + + entry = runs.setdefault(dataset_of(d["index_filename"]), + {"pos": [], "neg": [], "access": [], "iter": []}) + entry["pos"].append(float(d["positive lookup (avg_nanosec_per_kmer)"])) + entry["neg"].append(float(d["negative lookup (avg_nanosec_per_kmer)"])) + entry["access"].append(float(d["access (avg_nanosec_per_kmer)"])) + entry["iter"].append(float(d["iterator (avg_nanosec_per_kmer)"])) + + results = {} + for ds, v in runs.items(): + try: + results[ds] = { + "pos_us": f"{mean(v['pos']) / 1000:.2f}", + "neg_us": f"{mean(v['neg']) / 1000:.2f}", + "access_us": f"{mean(v['access']) / 1000:.2f}", + "iter_ns": f"{mean(v['iter']):.2f}", + } + except StatisticsError: + results[ds] = {"pos_us": "NA", "neg_us": "NA", "access_us": "NA", "iter_ns": "NA"} + return results + +def parse_streaming_file(path): + """Parse streaming queries JSONL file: one record per dataset.""" + results = {} + if not os.path.exists(path): + return results + + with open(path) as f: + for line in f: + line = line.strip() + if not line: + continue + try: + d = json.loads(line) + except json.JSONDecodeError: + print(f"Skipping invalid JSON line in {path}", file=sys.stderr) + continue + + num_kmers = int(d["num_kmers"]) + num_pos = int(d["num_positive_kmers"]) + num_ext = int(d["num_extensions"]) + elapsed_ms = int(d["elapsed_millisec"]) + + results[dataset_of(d["index_filename"])] = { + "ns_per_kmer": f"{int(math.ceil(elapsed_ms * 1e6 / num_kmers))}", + "hit_rate": f"{(num_pos / num_kmers) * 100 if num_kmers else 0:.2f}", + "extension_rate": f"{(num_ext / num_pos) * 100 if num_pos else 0:.2f}", + } + return results + +def collect(input_dir): + """One row per (k, dataset), in fixed order.""" + rows = [] + for k in K_VALUES: + k_dir = os.path.join(input_dir, f"k{k}") + if not os.path.isdir(k_dir): + print(f"Error: expected subdirectory '{k_dir}'", file=sys.stderr) + sys.exit(1) + build = parse_build_file(os.path.join(k_dir, "build.json")) + bench = parse_bench_file(os.path.join(k_dir, "bench.json")) + stream = parse_streaming_file(os.path.join(k_dir, "streaming-queries.json")) + first_row_of_k = len(rows) + for ds in DATASETS: + if ds not in build: + continue + row = {"k": str(k), "Collection": DISPLAY_NAME[ds], "first_of_k": False} + row.update(build[ds]) + row.update(bench.get(ds, {"pos_us": "NA", "neg_us": "NA", + "access_us": "NA", "iter_ns": "NA"})) + row.update(stream.get(ds, {"ns_per_kmer": "NA", "hit_rate": "NA", + "extension_rate": "NA"})) + rows.append(row) + if len(rows) > first_row_of_k: + rows[first_row_of_k]["first_of_k"] = True + return rows + +def print_csv(rows): + print("k,Collection,m,bits_per_kmer,total_GB,build_time," + "positive_lookup_us,negative_lookup_us,access_us,iteration_ns," + "streaming_ns_per_kmer,hit_rate,extension_rate") + for r in rows: + print(f"{r['k']},{r['Collection']},{r['m']},{r['bits_per_kmer']},{r['total_GB']}," + f"{r['build_time']},{r['pos_us']},{r['neg_us']},{r['access_us']},{r['iter_ns']}," + f"{r['ns_per_kmer']},{r['hit_rate']},{r['extension_rate']}") + +MD_COLUMNS = [ + # (header, row key, alignment: 'left' or 'center') + ("k", "k", "left"), + ("Collection", "Collection", "left"), + ("m", "m", "center"), + ("Space (bits/kmer)", "bits_per_kmer", "center"), + ("Space (total GB)", "total_GB", "center"), + ("Building time (m:ss)", "build_time", "center"), + ("Positive random lookup (µs/kmer)", "pos_us", "center"), + ("Negative random lookup (µs/kmer)", "neg_us", "center"), + ("Random Access (µs/kmer)", "access_us", "center"), + ("Streaming Lookup high-hit (ns/kmer)", "ns_per_kmer", "center"), +] + +def print_md(rows): + """The k value is shown on the first row of its block only, + and each block is preceded by an empty '||' row.""" + cells = [[str(r[key]) if key != "k" or r["first_of_k"] else "" + for _, key, _ in MD_COLUMNS] for r in rows] + widths = [max(len(header), *(len(row[i]) for row in cells)) + for i, (header, _, _) in enumerate(MD_COLUMNS)] + + def line(values): + return "| " + " | ".join(v.ljust(w) for v, w in zip(values, widths)) + " |" + + print(line([header for header, _, _ in MD_COLUMNS])) + print("|" + "|".join(f":{'-' * w}:" if align == "center" else "-" * (w + 2) + for w, (_, _, align) in zip(widths, MD_COLUMNS)) + "|") + for r, row_cells in zip(rows, cells): + if r["first_of_k"]: + print("||") + print(line(row_cells)) + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("input_dir", help="directory containing the k31 and k63 subdirectories") + parser.add_argument("--md", action="store_true", + help="print a markdown table instead of CSV") + args = parser.parse_args() + + rows = collect(args.input_dir) + if args.md: + print_md(rows) + else: + print_csv(rows) + +if __name__ == "__main__": + main() From 4d6d73b9cacb9a6aa4cf6fdf463cd12bf24a84b0 Mon Sep 17 00:00:00 2001 From: jermp Date: Sat, 29 Aug 2026 10:33:12 +0200 Subject: [PATCH 14/14] added new run logs and updated table with results --- benchmarks/README.md | 20 ++++- benchmarks/results-28-08-26/k31/bench.json | 18 ++++ benchmarks/results-28-08-26/k31/bench.log | 90 +++++++++++++++++++ .../k31/streaming-queries.json | 6 ++ .../k31/streaming-queries.log | 78 ++++++++++++++++ benchmarks/results-28-08-26/k63/bench.json | 18 ++++ benchmarks/results-28-08-26/k63/bench.log | 90 +++++++++++++++++++ .../k63/streaming-queries.json | 6 ++ .../k63/streaming-queries.log | 78 ++++++++++++++++ 9 files changed, 401 insertions(+), 3 deletions(-) create mode 100644 benchmarks/results-28-08-26/k31/bench.json create mode 100644 benchmarks/results-28-08-26/k31/bench.log create mode 100644 benchmarks/results-28-08-26/k31/streaming-queries.json create mode 100644 benchmarks/results-28-08-26/k31/streaming-queries.log create mode 100644 benchmarks/results-28-08-26/k63/bench.json create mode 100644 benchmarks/results-28-08-26/k63/bench.log create mode 100644 benchmarks/results-28-08-26/k63/streaming-queries.json create mode 100644 benchmarks/results-28-08-26/k63/streaming-queries.log diff --git a/benchmarks/README.md b/benchmarks/README.md index 9aa2b6f..c17b87e 100644 --- a/benchmarks/README.md +++ b/benchmarks/README.md @@ -14,14 +14,28 @@ To run the benchmarks, from within the `build` directory, run where `` should be replaced by a suitable basename, e.g., the current date. -These are the results obtained on 21/01/26 (see logs [here](results-21-01-26)) +These are the results obtained on 28/08/26 (see logs [here](results-28-08-26)) on a machine equipped with an AMD Ryzen Threadripper PRO 7985WX processor clocked at 5.40GHz. The code was compiled with `gcc` 13.3.0. - The indexes were build with a max RAM usage of 16 GB and 64 threads. Queries were run using one thread, instead. -![](results-21-01-26/results.png) +| k | Collection | m | Space (bits/kmer) | Space (total GB) | Building time (m:ss) | Positive random lookup (µs/kmer) | Negative random lookup (µs/kmer) | Random Access (µs/kmer) | Streaming Lookup high-hit (ns/kmer) | +|----|------------|:--:|:-----------------:|:----------------:|:--------------------:|:--------------------------------:|:--------------------------------:|:-----------------------:|:-----------------------------------:| +|| +| 31 | Cod | 20 | 7.94 | 0.50 | 0:27 | 0.43 | 0.36 | 0.27 | 24 | +| | Kestrel | 20 | 7.52 | 1.08 | 1:00 | 0.44 | 0.39 | 0.28 | 39 | +| | Human | 21 | 8.77 | 2.75 | 2:35 | 0.56 | 0.43 | 0.35 | 76 | +| | NCBI-v | 19 | 7.43 | 0.35 | 0:14 | 0.41 | 0.35 | 0.26 | 26 | +| | SE | 21 | 10.24 | 1.14 | 0:54 | 0.62 | 0.40 | 0.36 | 163 | +| | HPRC | 21 | 10.64 | 4.94 | 4:36 | 0.70 | 0.46 | 0.53 | 84 | +|| +| 63 | Cod | 24 | 4.62 | 0.32 | 0:24 | 0.57 | 0.45 | 0.29 | 37 | +| | Kestrel | 24 | 3.84 | 0.55 | 0:19 | 0.54 | 0.48 | 0.32 | 44 | +| | Human | 25 | 4.92 | 1.70 | 1:13 | 0.66 | 0.51 | 0.36 | 110 | +| | NCBI-v | 23 | 4.06 | 0.21 | 0:07 | 0.53 | 0.44 | 0.28 | 30 | +| | SE | 31 | 7.42 | 1.41 | 1:20 | 0.97 | 0.50 | 0.40 | 278 | +| | HPRC | 31 | 7.63 | 5.65 | 4:46 | 0.94 | 0.58 | 0.62 | 137 | The results of a run, for both k=31 and k=63 at once, can be exported to CSV format (every collected quantity) or to a markdown table (the main quantities) diff --git a/benchmarks/results-28-08-26/k31/bench.json b/benchmarks/results-28-08-26/k31/bench.json new file mode 100644 index 0000000..ffdc5ec --- /dev/null +++ b/benchmarks/results-28-08-26/k31/bench.json @@ -0,0 +1,18 @@ +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k31.sshash", "k": "31", "m": "20", "positive lookup (avg_nanosec_per_kmer)": "436.318757", "negative lookup (avg_nanosec_per_kmer)": "362.749596", "access (avg_nanosec_per_kmer)": "272.585695", "iterator (avg_nanosec_per_kmer)": "2.482751"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k31.sshash", "k": "31", "m": "20", "positive lookup (avg_nanosec_per_kmer)": "433.819361", "negative lookup (avg_nanosec_per_kmer)": "362.259503", "access (avg_nanosec_per_kmer)": "274.144289", "iterator (avg_nanosec_per_kmer)": "2.456011"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k31.sshash", "k": "31", "m": "20", "positive lookup (avg_nanosec_per_kmer)": "433.500928", "negative lookup (avg_nanosec_per_kmer)": "364.396456", "access (avg_nanosec_per_kmer)": "273.884328", "iterator (avg_nanosec_per_kmer)": "2.569083"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k31.sshash", "k": "31", "m": "20", "positive lookup (avg_nanosec_per_kmer)": "439.680768", "negative lookup (avg_nanosec_per_kmer)": "395.603409", "access (avg_nanosec_per_kmer)": "276.121400", "iterator (avg_nanosec_per_kmer)": "2.499332"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k31.sshash", "k": "31", "m": "20", "positive lookup (avg_nanosec_per_kmer)": "444.361067", "negative lookup (avg_nanosec_per_kmer)": "394.031758", "access (avg_nanosec_per_kmer)": "274.707766", "iterator (avg_nanosec_per_kmer)": "2.524888"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k31.sshash", "k": "31", "m": "20", "positive lookup (avg_nanosec_per_kmer)": "440.037087", "negative lookup (avg_nanosec_per_kmer)": "393.136266", "access (avg_nanosec_per_kmer)": "275.896044", "iterator (avg_nanosec_per_kmer)": "2.538501"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k31.sshash", "k": "31", "m": "21", "positive lookup (avg_nanosec_per_kmer)": "562.807039", "negative lookup (avg_nanosec_per_kmer)": "431.088550", "access (avg_nanosec_per_kmer)": "347.198096", "iterator (avg_nanosec_per_kmer)": "2.465738"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k31.sshash", "k": "31", "m": "21", "positive lookup (avg_nanosec_per_kmer)": "555.812279", "negative lookup (avg_nanosec_per_kmer)": "422.575776", "access (avg_nanosec_per_kmer)": "348.303889", "iterator (avg_nanosec_per_kmer)": "2.496474"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k31.sshash", "k": "31", "m": "21", "positive lookup (avg_nanosec_per_kmer)": "561.630325", "negative lookup (avg_nanosec_per_kmer)": "424.718263", "access (avg_nanosec_per_kmer)": "347.852823", "iterator (avg_nanosec_per_kmer)": "2.547386"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k31.sshash", "k": "31", "m": "19", "positive lookup (avg_nanosec_per_kmer)": "412.510594", "negative lookup (avg_nanosec_per_kmer)": "348.734417", "access (avg_nanosec_per_kmer)": "259.806846", "iterator (avg_nanosec_per_kmer)": "2.500924"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k31.sshash", "k": "31", "m": "19", "positive lookup (avg_nanosec_per_kmer)": "407.386907", "negative lookup (avg_nanosec_per_kmer)": "353.408523", "access (avg_nanosec_per_kmer)": "260.417569", "iterator (avg_nanosec_per_kmer)": "2.459740"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k31.sshash", "k": "31", "m": "19", "positive lookup (avg_nanosec_per_kmer)": "408.124222", "negative lookup (avg_nanosec_per_kmer)": "346.954378", "access (avg_nanosec_per_kmer)": "260.970479", "iterator (avg_nanosec_per_kmer)": "2.504823"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k31.sshash", "k": "31", "m": "21", "positive lookup (avg_nanosec_per_kmer)": "625.169634", "negative lookup (avg_nanosec_per_kmer)": "396.027242", "access (avg_nanosec_per_kmer)": "357.092978", "iterator (avg_nanosec_per_kmer)": "2.475710"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k31.sshash", "k": "31", "m": "21", "positive lookup (avg_nanosec_per_kmer)": "614.495341", "negative lookup (avg_nanosec_per_kmer)": "396.575080", "access (avg_nanosec_per_kmer)": "354.582941", "iterator (avg_nanosec_per_kmer)": "2.414229"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k31.sshash", "k": "31", "m": "21", "positive lookup (avg_nanosec_per_kmer)": "623.312167", "negative lookup (avg_nanosec_per_kmer)": "400.475688", "access (avg_nanosec_per_kmer)": "357.971520", "iterator (avg_nanosec_per_kmer)": "2.508970"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k31.sshash", "k": "31", "m": "21", "positive lookup (avg_nanosec_per_kmer)": "696.768429", "negative lookup (avg_nanosec_per_kmer)": "459.656809", "access (avg_nanosec_per_kmer)": "528.457965", "iterator (avg_nanosec_per_kmer)": "2.471353"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k31.sshash", "k": "31", "m": "21", "positive lookup (avg_nanosec_per_kmer)": "699.498620", "negative lookup (avg_nanosec_per_kmer)": "461.383063", "access (avg_nanosec_per_kmer)": "527.811735", "iterator (avg_nanosec_per_kmer)": "2.509292"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k31.sshash", "k": "31", "m": "21", "positive lookup (avg_nanosec_per_kmer)": "696.245840", "negative lookup (avg_nanosec_per_kmer)": "458.446995", "access (avg_nanosec_per_kmer)": "521.370161", "iterator (avg_nanosec_per_kmer)": "2.440193"} diff --git a/benchmarks/results-28-08-26/k31/bench.log b/benchmarks/results-28-08-26/k31/bench.log new file mode 100644 index 0000000..5259139 --- /dev/null +++ b/benchmarks/results-28-08-26/k31/bench.log @@ -0,0 +1,90 @@ +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 436.319 +negative lookup (avg_nanosec_per_kmer) 362.75 +access (avg_nanosec_per_kmer) = 272.586 +iterator (avg_nanosec_per_kmer) = 2.48275 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 433.819 +negative lookup (avg_nanosec_per_kmer) 362.26 +access (avg_nanosec_per_kmer) = 274.144 +iterator (avg_nanosec_per_kmer) = 2.45601 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 433.501 +negative lookup (avg_nanosec_per_kmer) 364.396 +access (avg_nanosec_per_kmer) = 273.884 +iterator (avg_nanosec_per_kmer) = 2.56908 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 439.681 +negative lookup (avg_nanosec_per_kmer) 395.603 +access (avg_nanosec_per_kmer) = 276.121 +iterator (avg_nanosec_per_kmer) = 2.49933 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 444.361 +negative lookup (avg_nanosec_per_kmer) 394.032 +access (avg_nanosec_per_kmer) = 274.708 +iterator (avg_nanosec_per_kmer) = 2.52489 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 440.037 +negative lookup (avg_nanosec_per_kmer) 393.136 +access (avg_nanosec_per_kmer) = 275.896 +iterator (avg_nanosec_per_kmer) = 2.5385 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 562.807 +negative lookup (avg_nanosec_per_kmer) 431.089 +access (avg_nanosec_per_kmer) = 347.198 +iterator (avg_nanosec_per_kmer) = 2.46574 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 555.812 +negative lookup (avg_nanosec_per_kmer) 422.576 +access (avg_nanosec_per_kmer) = 348.304 +iterator (avg_nanosec_per_kmer) = 2.49647 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 561.63 +negative lookup (avg_nanosec_per_kmer) 424.718 +access (avg_nanosec_per_kmer) = 347.853 +iterator (avg_nanosec_per_kmer) = 2.54739 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 412.511 +negative lookup (avg_nanosec_per_kmer) 348.734 +access (avg_nanosec_per_kmer) = 259.807 +iterator (avg_nanosec_per_kmer) = 2.50092 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 407.387 +negative lookup (avg_nanosec_per_kmer) 353.409 +access (avg_nanosec_per_kmer) = 260.418 +iterator (avg_nanosec_per_kmer) = 2.45974 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 408.124 +negative lookup (avg_nanosec_per_kmer) 346.954 +access (avg_nanosec_per_kmer) = 260.97 +iterator (avg_nanosec_per_kmer) = 2.50482 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 625.17 +negative lookup (avg_nanosec_per_kmer) 396.027 +access (avg_nanosec_per_kmer) = 357.093 +iterator (avg_nanosec_per_kmer) = 2.47571 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 614.495 +negative lookup (avg_nanosec_per_kmer) 396.575 +access (avg_nanosec_per_kmer) = 354.583 +iterator (avg_nanosec_per_kmer) = 2.41423 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 623.312 +negative lookup (avg_nanosec_per_kmer) 400.476 +access (avg_nanosec_per_kmer) = 357.972 +iterator (avg_nanosec_per_kmer) = 2.50897 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 696.768 +negative lookup (avg_nanosec_per_kmer) 459.657 +access (avg_nanosec_per_kmer) = 528.458 +iterator (avg_nanosec_per_kmer) = 2.47135 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 699.499 +negative lookup (avg_nanosec_per_kmer) 461.383 +access (avg_nanosec_per_kmer) = 527.812 +iterator (avg_nanosec_per_kmer) = 2.50929 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k31.sshash +positive lookup (avg_nanosec_per_kmer) = 696.246 +negative lookup (avg_nanosec_per_kmer) 458.447 +access (avg_nanosec_per_kmer) = 521.37 +iterator (avg_nanosec_per_kmer) = 2.44019 diff --git a/benchmarks/results-28-08-26/k31/streaming-queries.json b/benchmarks/results-28-08-26/k31/streaming-queries.json new file mode 100644 index 0000000..d7cbe12 --- /dev/null +++ b/benchmarks/results-28-08-26/k31/streaming-queries.json @@ -0,0 +1,6 @@ +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k31.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR12858649.fastq.gz", "num_kmers": "163287360", "num_positive_kmers": "132860997", "num_negative_kmers": "30426363", "num_invalid_kmers": "0", "num_searches": "6576340", "num_extensions": "126284657", "num_memo_singleton": "193441", "num_memo_light": "6812848", "elapsed_millisec": "3770"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k31.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR11449743_1.fastq.gz", "num_kmers": "695737535", "num_positive_kmers": "525542891", "num_negative_kmers": "170183654", "num_invalid_kmers": "10990", "num_searches": "12437476", "num_extensions": "513105415", "num_memo_singleton": "997817", "num_memo_light": "34475911", "elapsed_millisec": "26537"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k31.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz", "num_kmers": "1569974986", "num_positive_kmers": "1437949378", "num_negative_kmers": "130996597", "num_invalid_kmers": "1029011", "num_searches": "100222623", "num_extensions": "1337726755", "num_memo_singleton": "1072506", "num_memo_light": "27289697", "elapsed_millisec": "119021"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k31.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/ncbi-queries.fastq.gz", "num_kmers": "14092875", "num_positive_kmers": "13983775", "num_negative_kmers": "108161", "num_invalid_kmers": "939", "num_searches": "590894", "num_extensions": "13392881", "num_memo_singleton": "433", "num_memo_light": "107455", "elapsed_millisec": "356"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k31.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR27871075_1.fastq.gz", "num_kmers": "789838196", "num_positive_kmers": "764882549", "num_negative_kmers": "24935381", "num_invalid_kmers": "20266", "num_searches": "218875709", "num_extensions": "546006840", "num_memo_singleton": "85669", "num_memo_light": "7940520", "elapsed_millisec": "128190"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k31.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz", "num_kmers": "1569974986", "num_positive_kmers": "1485223278", "num_negative_kmers": "83722697", "num_invalid_kmers": "1029011", "num_searches": "135823240", "num_extensions": "1349400038", "num_memo_singleton": "565997", "num_memo_light": "17514165", "elapsed_millisec": "131342"} diff --git a/benchmarks/results-28-08-26/k31/streaming-queries.log b/benchmarks/results-28-08-26/k31/streaming-queries.log new file mode 100644 index 0000000..b2b5e62 --- /dev/null +++ b/benchmarks/results-28-08-26/k31/streaming-queries.log @@ -0,0 +1,78 @@ +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k31.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR12858649.fastq.gz +2026-08-28 20:34:21: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR12858649.fastq.gz'... +2026-08-28 20:34:24: DONE +==== query report: +num_kmers = 163287360 +num_positive_kmers = 132860997 (81.3664%) +num_negative_kmers = 30426363 (18.6336%) +num_invalid_kmers = 0 (0%) +num_searches = 6576340/132860997 (4.94979%) +num_extensions = 126284657/132860997 (95.0502%) +num_memo_singleton = 193441 +num_memo_light = 6812848 +elapsed = 3.77 sec / 0.0628333 min / 23.0881 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k31.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR11449743_1.fastq.gz +2026-08-28 20:34:25: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR11449743_1.fastq.gz'... +2026-08-28 20:34:52: DONE +==== query report: +num_kmers = 695737535 +num_positive_kmers = 525542891 (75.5375%) +num_negative_kmers = 170183654 (24.4609%) +num_invalid_kmers = 10990 (0.00157962%) +num_searches = 12437476/525542891 (2.3666%) +num_extensions = 513105415/525542891 (97.6334%) +num_memo_singleton = 997817 +num_memo_light = 34475911 +elapsed = 26.537 sec / 0.442283 min / 38.1423 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k31.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz +2026-08-28 20:34:53: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz'... +2026-08-28 20:36:52: DONE +==== query report: +num_kmers = 1569974986 +num_positive_kmers = 1437949378 (91.5906%) +num_negative_kmers = 130996597 (8.34387%) +num_invalid_kmers = 1029011 (0.0655431%) +num_searches = 100222623/1437949378 (6.96983%) +num_extensions = 1337726755/1437949378 (93.0302%) +num_memo_singleton = 1072506 +num_memo_light = 27289697 +elapsed = 119.021 sec / 1.98368 min / 75.8108 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k31.sshash -q /mnt/hd2/pibiri/DNA/queries/ncbi-queries.fastq.gz +2026-08-28 20:36:52: performing queries from file '/mnt/hd2/pibiri/DNA/queries/ncbi-queries.fastq.gz'... +2026-08-28 20:36:53: DONE +==== query report: +num_kmers = 14092875 +num_positive_kmers = 13983775 (99.2258%) +num_negative_kmers = 108161 (0.767487%) +num_invalid_kmers = 939 (0.00666294%) +num_searches = 590894/13983775 (4.22557%) +num_extensions = 13392881/13983775 (95.7744%) +num_memo_singleton = 433 +num_memo_light = 107455 +elapsed = 0.356 sec / 0.00593333 min / 25.261 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k31.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR27871075_1.fastq.gz +2026-08-28 20:36:53: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR27871075_1.fastq.gz'... +2026-08-28 20:39:01: DONE +==== query report: +num_kmers = 789838196 +num_positive_kmers = 764882549 (96.8404%) +num_negative_kmers = 24935381 (3.15702%) +num_invalid_kmers = 20266 (0.00256584%) +num_searches = 218875709/764882549 (28.6156%) +num_extensions = 546006840/764882549 (71.3844%) +num_memo_singleton = 85669 +num_memo_light = 7940520 +elapsed = 128.19 sec / 2.1365 min / 162.299 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k31.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz +2026-08-28 20:39:04: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz'... +2026-08-28 20:41:15: DONE +==== query report: +num_kmers = 1569974986 +num_positive_kmers = 1485223278 (94.6017%) +num_negative_kmers = 83722697 (5.33274%) +num_invalid_kmers = 1029011 (0.0655431%) +num_searches = 135823240/1485223278 (9.14497%) +num_extensions = 1349400038/1485223278 (90.855%) +num_memo_singleton = 565997 +num_memo_light = 17514165 +elapsed = 131.342 sec / 2.18903 min / 83.6587 ns/kmer diff --git a/benchmarks/results-28-08-26/k63/bench.json b/benchmarks/results-28-08-26/k63/bench.json new file mode 100644 index 0000000..49edb22 --- /dev/null +++ b/benchmarks/results-28-08-26/k63/bench.json @@ -0,0 +1,18 @@ +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k63.sshash", "k": "63", "m": "24", "positive lookup (avg_nanosec_per_kmer)": "567.271179", "negative lookup (avg_nanosec_per_kmer)": "447.896700", "access (avg_nanosec_per_kmer)": "292.035089", "iterator (avg_nanosec_per_kmer)": "2.723525"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k63.sshash", "k": "63", "m": "24", "positive lookup (avg_nanosec_per_kmer)": "567.158980", "negative lookup (avg_nanosec_per_kmer)": "442.781347", "access (avg_nanosec_per_kmer)": "290.059058", "iterator (avg_nanosec_per_kmer)": "2.716768"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k63.sshash", "k": "63", "m": "24", "positive lookup (avg_nanosec_per_kmer)": "577.632402", "negative lookup (avg_nanosec_per_kmer)": "445.099827", "access (avg_nanosec_per_kmer)": "289.093225", "iterator (avg_nanosec_per_kmer)": "2.713838"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k63.sshash", "k": "63", "m": "24", "positive lookup (avg_nanosec_per_kmer)": "536.517146", "negative lookup (avg_nanosec_per_kmer)": "479.223758", "access (avg_nanosec_per_kmer)": "324.438066", "iterator (avg_nanosec_per_kmer)": "2.757624"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k63.sshash", "k": "63", "m": "24", "positive lookup (avg_nanosec_per_kmer)": "537.308816", "negative lookup (avg_nanosec_per_kmer)": "482.519123", "access (avg_nanosec_per_kmer)": "327.085464", "iterator (avg_nanosec_per_kmer)": "2.705321"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k63.sshash", "k": "63", "m": "24", "positive lookup (avg_nanosec_per_kmer)": "543.875553", "negative lookup (avg_nanosec_per_kmer)": "484.284402", "access (avg_nanosec_per_kmer)": "321.921078", "iterator (avg_nanosec_per_kmer)": "2.708895"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k63.sshash", "k": "63", "m": "25", "positive lookup (avg_nanosec_per_kmer)": "655.823771", "negative lookup (avg_nanosec_per_kmer)": "512.844700", "access (avg_nanosec_per_kmer)": "355.444918", "iterator (avg_nanosec_per_kmer)": "2.719615"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k63.sshash", "k": "63", "m": "25", "positive lookup (avg_nanosec_per_kmer)": "657.564050", "negative lookup (avg_nanosec_per_kmer)": "514.922087", "access (avg_nanosec_per_kmer)": "356.297916", "iterator (avg_nanosec_per_kmer)": "2.773051"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k63.sshash", "k": "63", "m": "25", "positive lookup (avg_nanosec_per_kmer)": "659.735498", "negative lookup (avg_nanosec_per_kmer)": "510.912399", "access (avg_nanosec_per_kmer)": "355.614307", "iterator (avg_nanosec_per_kmer)": "2.715916"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k63.sshash", "k": "63", "m": "23", "positive lookup (avg_nanosec_per_kmer)": "531.164346", "negative lookup (avg_nanosec_per_kmer)": "434.532611", "access (avg_nanosec_per_kmer)": "281.460643", "iterator (avg_nanosec_per_kmer)": "2.733554"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k63.sshash", "k": "63", "m": "23", "positive lookup (avg_nanosec_per_kmer)": "525.985610", "negative lookup (avg_nanosec_per_kmer)": "437.954528", "access (avg_nanosec_per_kmer)": "277.253213", "iterator (avg_nanosec_per_kmer)": "2.727797"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k63.sshash", "k": "63", "m": "23", "positive lookup (avg_nanosec_per_kmer)": "531.091279", "negative lookup (avg_nanosec_per_kmer)": "434.230904", "access (avg_nanosec_per_kmer)": "283.331560", "iterator (avg_nanosec_per_kmer)": "2.727683"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k63.sshash", "k": "63", "m": "31", "positive lookup (avg_nanosec_per_kmer)": "965.438332", "negative lookup (avg_nanosec_per_kmer)": "495.346946", "access (avg_nanosec_per_kmer)": "402.682535", "iterator (avg_nanosec_per_kmer)": "2.755412"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k63.sshash", "k": "63", "m": "31", "positive lookup (avg_nanosec_per_kmer)": "969.421164", "negative lookup (avg_nanosec_per_kmer)": "494.920548", "access (avg_nanosec_per_kmer)": "404.439983", "iterator (avg_nanosec_per_kmer)": "2.749083"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k63.sshash", "k": "63", "m": "31", "positive lookup (avg_nanosec_per_kmer)": "961.154818", "negative lookup (avg_nanosec_per_kmer)": "498.388122", "access (avg_nanosec_per_kmer)": "400.903357", "iterator (avg_nanosec_per_kmer)": "2.820549"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k63.sshash", "k": "63", "m": "31", "positive lookup (avg_nanosec_per_kmer)": "941.367280", "negative lookup (avg_nanosec_per_kmer)": "581.844150", "access (avg_nanosec_per_kmer)": "622.161489", "iterator (avg_nanosec_per_kmer)": "2.827105"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k63.sshash", "k": "63", "m": "31", "positive lookup (avg_nanosec_per_kmer)": "944.656877", "negative lookup (avg_nanosec_per_kmer)": "588.830562", "access (avg_nanosec_per_kmer)": "621.531063", "iterator (avg_nanosec_per_kmer)": "2.838949"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k63.sshash", "k": "63", "m": "31", "positive lookup (avg_nanosec_per_kmer)": "942.756666", "negative lookup (avg_nanosec_per_kmer)": "582.045724", "access (avg_nanosec_per_kmer)": "624.022485", "iterator (avg_nanosec_per_kmer)": "2.915508"} diff --git a/benchmarks/results-28-08-26/k63/bench.log b/benchmarks/results-28-08-26/k63/bench.log new file mode 100644 index 0000000..dc29412 --- /dev/null +++ b/benchmarks/results-28-08-26/k63/bench.log @@ -0,0 +1,90 @@ +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 567.271 +negative lookup (avg_nanosec_per_kmer) 447.897 +access (avg_nanosec_per_kmer) = 292.035 +iterator (avg_nanosec_per_kmer) = 2.72352 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 567.159 +negative lookup (avg_nanosec_per_kmer) 442.781 +access (avg_nanosec_per_kmer) = 290.059 +iterator (avg_nanosec_per_kmer) = 2.71677 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 577.632 +negative lookup (avg_nanosec_per_kmer) 445.1 +access (avg_nanosec_per_kmer) = 289.093 +iterator (avg_nanosec_per_kmer) = 2.71384 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 536.517 +negative lookup (avg_nanosec_per_kmer) 479.224 +access (avg_nanosec_per_kmer) = 324.438 +iterator (avg_nanosec_per_kmer) = 2.75762 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 537.309 +negative lookup (avg_nanosec_per_kmer) 482.519 +access (avg_nanosec_per_kmer) = 327.085 +iterator (avg_nanosec_per_kmer) = 2.70532 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 543.876 +negative lookup (avg_nanosec_per_kmer) 484.284 +access (avg_nanosec_per_kmer) = 321.921 +iterator (avg_nanosec_per_kmer) = 2.70889 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 655.824 +negative lookup (avg_nanosec_per_kmer) 512.845 +access (avg_nanosec_per_kmer) = 355.445 +iterator (avg_nanosec_per_kmer) = 2.71961 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 657.564 +negative lookup (avg_nanosec_per_kmer) 514.922 +access (avg_nanosec_per_kmer) = 356.298 +iterator (avg_nanosec_per_kmer) = 2.77305 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 659.735 +negative lookup (avg_nanosec_per_kmer) 510.912 +access (avg_nanosec_per_kmer) = 355.614 +iterator (avg_nanosec_per_kmer) = 2.71592 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 531.164 +negative lookup (avg_nanosec_per_kmer) 434.533 +access (avg_nanosec_per_kmer) = 281.461 +iterator (avg_nanosec_per_kmer) = 2.73355 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 525.986 +negative lookup (avg_nanosec_per_kmer) 437.955 +access (avg_nanosec_per_kmer) = 277.253 +iterator (avg_nanosec_per_kmer) = 2.7278 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 531.091 +negative lookup (avg_nanosec_per_kmer) 434.231 +access (avg_nanosec_per_kmer) = 283.332 +iterator (avg_nanosec_per_kmer) = 2.72768 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 965.438 +negative lookup (avg_nanosec_per_kmer) 495.347 +access (avg_nanosec_per_kmer) = 402.683 +iterator (avg_nanosec_per_kmer) = 2.75541 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 969.421 +negative lookup (avg_nanosec_per_kmer) 494.921 +access (avg_nanosec_per_kmer) = 404.44 +iterator (avg_nanosec_per_kmer) = 2.74908 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 961.155 +negative lookup (avg_nanosec_per_kmer) 498.388 +access (avg_nanosec_per_kmer) = 400.903 +iterator (avg_nanosec_per_kmer) = 2.82055 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 941.367 +negative lookup (avg_nanosec_per_kmer) 581.844 +access (avg_nanosec_per_kmer) = 622.161 +iterator (avg_nanosec_per_kmer) = 2.8271 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 944.657 +negative lookup (avg_nanosec_per_kmer) 588.831 +access (avg_nanosec_per_kmer) = 621.531 +iterator (avg_nanosec_per_kmer) = 2.83895 +./sshash bench -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k63.sshash +positive lookup (avg_nanosec_per_kmer) = 942.757 +negative lookup (avg_nanosec_per_kmer) 582.046 +access (avg_nanosec_per_kmer) = 624.022 +iterator (avg_nanosec_per_kmer) = 2.91551 diff --git a/benchmarks/results-28-08-26/k63/streaming-queries.json b/benchmarks/results-28-08-26/k63/streaming-queries.json new file mode 100644 index 0000000..decc135 --- /dev/null +++ b/benchmarks/results-28-08-26/k63/streaming-queries.json @@ -0,0 +1,6 @@ +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k63.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR12858649.fastq.gz", "num_kmers": "97972416", "num_positive_kmers": "67275966", "num_negative_kmers": "30696450", "num_invalid_kmers": "0", "num_searches": "2280719", "num_extensions": "64995247", "num_memo_singleton": "1309608", "num_memo_light": "11103991", "elapsed_millisec": "3569"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k63.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR11449743_1.fastq.gz", "num_kmers": "461383839", "num_positive_kmers": "293470517", "num_negative_kmers": "167902332", "num_invalid_kmers": "10990", "num_searches": "6395458", "num_extensions": "287075059", "num_memo_singleton": "4872310", "num_memo_light": "56564763", "elapsed_millisec": "20048"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k63.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz", "num_kmers": "477818474", "num_positive_kmers": "406529529", "num_negative_kmers": "70615167", "num_invalid_kmers": "673778", "num_searches": "33432030", "num_extensions": "373097499", "num_memo_singleton": "3056220", "num_memo_light": "23887860", "elapsed_millisec": "52305"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k63.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/ncbi-queries.fastq.gz", "num_kmers": "10330949", "num_positive_kmers": "10230224", "num_negative_kmers": "99451", "num_invalid_kmers": "1274", "num_searches": "279697", "num_extensions": "9950527", "num_memo_singleton": "1727", "num_memo_light": "53202", "elapsed_millisec": "304"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k63.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR27871075_1.fastq.gz", "num_kmers": "541466405", "num_positive_kmers": "507202856", "num_negative_kmers": "34238416", "num_invalid_kmers": "25133", "num_searches": "131137192", "num_extensions": "376065664", "num_memo_singleton": "126573", "num_memo_light": "3425578", "elapsed_millisec": "150386"} +{"index_filename": "/mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k63.sshash", "query_filename": "/mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz", "num_kmers": "477818474", "num_positive_kmers": "434532302", "num_negative_kmers": "42612394", "num_invalid_kmers": "673778", "num_searches": "43895917", "num_extensions": "390636385", "num_memo_singleton": "877256", "num_memo_light": "13835358", "elapsed_millisec": "65020"} diff --git a/benchmarks/results-28-08-26/k63/streaming-queries.log b/benchmarks/results-28-08-26/k63/streaming-queries.log new file mode 100644 index 0000000..89312c6 --- /dev/null +++ b/benchmarks/results-28-08-26/k63/streaming-queries.log @@ -0,0 +1,78 @@ +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/cod.k63.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR12858649.fastq.gz +2026-08-28 20:41:29: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR12858649.fastq.gz'... +2026-08-28 20:41:32: DONE +==== query report: +num_kmers = 97972416 +num_positive_kmers = 67275966 (68.6683%) +num_negative_kmers = 30696450 (31.3317%) +num_invalid_kmers = 0 (0%) +num_searches = 2280719/67275966 (3.39009%) +num_extensions = 64995247/67275966 (96.6099%) +num_memo_singleton = 1309608 +num_memo_light = 11103991 +elapsed = 3.569 sec / 0.0594833 min / 36.4286 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/kestrel.k63.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR11449743_1.fastq.gz +2026-08-28 20:41:33: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR11449743_1.fastq.gz'... +2026-08-28 20:41:53: DONE +==== query report: +num_kmers = 461383839 +num_positive_kmers = 293470517 (63.6066%) +num_negative_kmers = 167902332 (36.391%) +num_invalid_kmers = 10990 (0.00238196%) +num_searches = 6395458/293470517 (2.17925%) +num_extensions = 287075059/293470517 (97.8207%) +num_memo_singleton = 4872310 +num_memo_light = 56564763 +elapsed = 20.048 sec / 0.334133 min / 43.4519 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/human.k63.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz +2026-08-28 20:41:53: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz'... +2026-08-28 20:42:46: DONE +==== query report: +num_kmers = 477818474 +num_positive_kmers = 406529529 (85.0803%) +num_negative_kmers = 70615167 (14.7787%) +num_invalid_kmers = 673778 (0.141011%) +num_searches = 33432030/406529529 (8.22376%) +num_extensions = 373097499/406529529 (91.7762%) +num_memo_singleton = 3056220 +num_memo_light = 23887860 +elapsed = 52.305 sec / 0.87175 min / 109.466 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/ncbi-virus.k63.sshash -q /mnt/hd2/pibiri/DNA/queries/ncbi-queries.fastq.gz +2026-08-28 20:42:46: performing queries from file '/mnt/hd2/pibiri/DNA/queries/ncbi-queries.fastq.gz'... +2026-08-28 20:42:46: DONE +==== query report: +num_kmers = 10330949 +num_positive_kmers = 10230224 (99.025%) +num_negative_kmers = 99451 (0.962651%) +num_invalid_kmers = 1274 (0.0123319%) +num_searches = 279697/10230224 (2.73403%) +num_extensions = 9950527/10230224 (97.266%) +num_memo_singleton = 1727 +num_memo_light = 53202 +elapsed = 0.304 sec / 0.00506667 min / 29.4261 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/se.k63.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR27871075_1.fastq.gz +2026-08-28 20:42:47: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR27871075_1.fastq.gz'... +2026-08-28 20:45:17: DONE +==== query report: +num_kmers = 541466405 +num_positive_kmers = 507202856 (93.6721%) +num_negative_kmers = 34238416 (6.32328%) +num_invalid_kmers = 25133 (0.00464165%) +num_searches = 131137192/507202856 (25.855%) +num_extensions = 376065664/507202856 (74.145%) +num_memo_singleton = 126573 +num_memo_light = 3425578 +elapsed = 150.386 sec / 2.50643 min / 277.738 ns/kmer +./sshash query -i /mnt/hd2/pibiri/DNA/sshash-indexes-new/hprc.k63.sshash -q /mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz +2026-08-28 20:45:20: performing queries from file '/mnt/hd2/pibiri/DNA/queries/SRR5833294.fastq.gz'... +2026-08-28 20:46:25: DONE +==== query report: +num_kmers = 477818474 +num_positive_kmers = 434532302 (90.9409%) +num_negative_kmers = 42612394 (8.91811%) +num_invalid_kmers = 673778 (0.141011%) +num_searches = 43895917/434532302 (10.1019%) +num_extensions = 390636385/434532302 (89.8981%) +num_memo_singleton = 877256 +num_memo_light = 13835358 +elapsed = 65.02 sec / 1.08367 min / 136.077 ns/kmer