Refactor SSHash indexing to unify parsing modalities and optimize performance - #93
Merged
Conversation
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.
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.
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.
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.
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.
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.
`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.
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.
`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.
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.
`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.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This pull request removes the canonical/regular indexing distinction, simplifying both the codebase and user experience. Now, the index always uses a forward minimizer scheme, so each minimizer position corresponds to a single tuple, and there is no
--canonicalflag or related space/speed trade-off. The build and benchmark scripts, documentation, and internal data structures are all updated to reflect this unified approach. Additionally, the code now enforces the forwardness property during tuple processing.Indexing Modality Simplification:
Removed all code, build options, and documentation related to the canonical/regular distinction; there is now only one indexing modality using a forward minimizer scheme (
--canonicalflag and related logic removed, documentation updated). [1] [2] [3] [4] [5] [6] [7] [8] [9]Updated the version number to 6.0.0 to reflect this breaking change.
Code and Data Structure Updates:
Refactored
bucket_typeandminimizers_tuplesto remove all logic and fields related to super-kmers and canonical minimizer positions; now, the number of tuples equals the number of minimizer positions, and this is checked for all input data. [1] [2] [3] [4] [5] [6]Added runtime checks to ensure the minimizer scheme is forward: no (minimizer, position) pair can occur more than once.
Benchmarking and Script Adjustments:
benchmarks/print_csv.pyto handle only the new single-modality result files, removing all code and CSV columns related to canonical/regular distinction. [1] [2] [3] [4] [5] [6] [7] [8]Build System:
test_minimizerto the build system.