Repository navigation
Conversation
The per-mate encoding loop fused the nt4 base translation with building the reverse-complement strand. The translation is a 256-entry table gather the compiler cannot vectorize, and fusing it kept the reverse-complement scalar as well. Splitting into two loops leaves the gather unchanged but lets the compiler auto-vectorize the reverse-complement (rev64+ext reverse, compare, and blend), about 1.4x faster on the encode at 150 bp reads, with no intrinsics. Output is unchanged: minibwa map produces byte-identical SAM on a 6000-pair sampled set, and the loop computes exactly the same values in the same order.
|
Retired to the graveyard; not shipping. Measured on The premise holds: the revcomp really does auto-vectorize, and identically on both architectures. The problem is that splitting the loop makes loop 2 reload with 16-byte vector loads the buffer that scalar loop 1 just wrote byte-by-byte, which defeats x86 store-to-load forwarding. c6a16 is SMT, so siblings contend for the store buffer; c8g16 has no SMT and tolerates the forwarding better, which is why it reads flat there. Moot either way: a profile puts the whole nt4 encode path under 0.1% of runtime -- no matching symbol appears at all -- so the region was never worth optimizing. The manifest now carries this as The branch stays so the disproof remains reproducible. |
Add two CLI flags that gate Pass-2 sub-SMEM reseeding (the inner loop
in mb_seed_intv / mb_seed_intv_batch) on the SA-interval size of each
Pass-1 SMEM:
--max-sub-occ N (default 10) upper bound: skip Pass-2 for SMEMs with
occurrence > N (existing internal default,
now user-configurable; N=0 disables Pass-2
entirely)
--min-sub-occ N (default 1) lower bound: skip Pass-2 for SMEMs with
occurrence < N; intended as an ablation
for the "skip Pass-2 on uniquely-anchored
SMEMs" hypothesis
Default values preserve current behavior bit-for-bit (Pass-1 only ever
appends SMEMs with size > 0, so size < 1 never fires). Plumbed through
the mb_seed_intv / mb_seed_intv_batch signatures (mbpriv.h, seed.c) and
the three callers in map-algo.c (mb_map_seq + mb_map_batch) and
map-main.c (streaming worker_for_se_batch).
CLI option codes are 312 / 313 (310 and 311 are taken by --meth and
--hic on this branch). Documented in usage() output. Both flags gate on
the same SMEM SA-interval size field; symmetric naming
(max-/min-sub-occ) reflects that.
Smoke-tested on test/chrM-human.fa.gz: default behavior is bit-identical
whether the flags are omitted or set to their defaults explicitly.
Empirical measurement on WGS-1M R1-only against GRCh38 (1,000,000 SE
Illumina reads, hs38 ref, 8 threads), reproducing the numbers from the
prior implementation on the methylation branch:
Pass-2 ON Pass-2 OFF predA (size>=2)
(default) (--max=0) (--min=2)
Wall 14.6 s 6.4 s 12.4 s
CPU 54.5 s 33.0 s 39.6 s
shift60 vs ON - 34,675 reads 316 of 34,675 (0.91%)
shift60 = same-placement reads where Pass-2 ON deflates MAPQ from 60 to
<60. The lower-bound predicate preserves only 0.91% of those deflations
-- the dominant case is "uniquely-seeded SMEM with paralog signal hidden
inside it" -- so --min-sub-occ is a measurement knob, not a recommended
production setting.
Follow-up to issue #10 (Pass-2 sub-SMEM reseeding reduces MAPQ60
recall). Implementation ported from the methylation branch.
What
The per-mate query encoding in
mb_map_pe1(pe.c) fused two operations into one loop: thekom_nt4_tablebase translation and building the reverse-complement strand. The translation is a 256-entry table gather the compiler cannot vectorize, and fusing it kept the reverse-complement scalar too. This splits the loop in two — the gather is unchanged, but the reverse-complement now auto-vectorizes (rev64+extreverse, compare-for-N, blend) under-O3. No intrinsics.Why
This started as an attempt to hand-vectorize the 2-bit reverse-complement with
vqtbl1q/pshufb+XOR 3. That turned out to be a dead end: Apple clang (and gcc) already auto-vectorize the scalar reverse-complement, choosingrev64+extover a table lookup, unrolling 4×, and vectorizing the tail — so hand-written intrinsics lost at 150 bp. Inspecting the generated asm showed the only place the compiler couldn't vectorize was here, where thent4gather poisoned the whole fused loop. Splitting it is the entire fix.Of the four per-read reverse-complement sites in the codebase, this is the only one that benefits: the
align.candseed.cloops already auto-vectorize as written (and un-fusingalign.cactually regresses it ~0.8×, since its fused copy+revcomp vectorizes in a single pass).Benefit
Validation
-O3confirmed to emit vectorized NEON for the reverse-complement loop after the split (and scalar before).