Fix batch enumeration offsets and refactor sibling annex construction

Shifts sibling annex construction from slot-indexed enumeration to iteration-order traversal by correcting cumulative k-mer offset tracking in batch enumeration. Replaces coarse per-partition parallelism with chunked work distribution to prevent thread starvation on skewed partitions. Decouples custom progress messages from ETA updates to eliminate display clobbering during high-frequency callbacks. Adds regression tests validating batch offset correctness, partial batch handling, and iterator-order consistency across layer builds.
This commit is contained in:
Eric Coissac
2026-08-16 20:51:19 +02:00
parent d2548e8c33
commit 32d6720f50
11 changed files with 413 additions and 88 deletions
+53 -8
View File
@@ -174,8 +174,16 @@ fn build_layer_sibling_annex(
// — a plain 1-to-1 (`|`, not `||`) transform, batch in, batch of
// variants out, one message either way — keeps the per-message
// synchronisation cost amortised over thousands of lookups instead
// of one to three. ──────────────────────────────────────────────
const BATCH_SIZE: usize = 4096;
// of one to three. `BATCH_SIZE` was originally tuned back when the
// per-k-mer cost inside this stage (minimiser via `RollingStat`)
// was ~3x higher than it is now (see `CanonicalKmerOf::minimizer`,
// a direct O(k) bit-arithmetic replacement) — at the old per-k-mer
// cost, dispatch overhead (the scheduler's single-threaded `Select`
// loop: one channel round-trip per batch) was negligible next to
// the work each batch represented; now that the work itself is
// cheaper, that fixed per-batch overhead is proportionally larger,
// so a bigger batch amortises it over more k-mers again. ─────────
const BATCH_SIZE: usize = 32768;
let n_workers = obisys::effective_parallelism();
let capacity = 256;
@@ -203,6 +211,16 @@ fn build_layer_sibling_annex(
})
.collect::<Vec<_>>()
});
// Diagnostic: count still-empty mask slots after seeding + cross-partition
// resolution. A non-zero count here means some k-mer's own base was never
// written (should be impossible after the batch_offset fix above).
let check_empty = |label: &str| {
let empty = mask.iter().filter(|v| v.load(Ordering::Relaxed) == 0).count();
if empty > 0 {
tracing::warn!("{label}: {empty} mask slots still empty after resolution (layer {layer_dir:?})");
}
};
let throttled = obipipeline::throttle(batches, n_workers).map(|t| SourceBatch {
items: t.item,
_permit: t.guard,
@@ -246,18 +264,45 @@ fn build_layer_sibling_annex(
}
}
// ── Resolve each partition's batch against the cache in one
// contiguous pass; parallelised across partitions (independent,
// read-only) so this keeps using multiple cores without giving up
// the per-partition locality above. ─────────────────────────────
outgoing.par_iter().enumerate().filter(|(_, q)| !q.is_empty()).for_each(|(dest, queries)| {
for &(variant, source_order, base) in queries {
check_empty("after_seeding");
// ── Resolve each partition's batch against the cache, parallelised
// over roughly equal-sized *chunks*, not over partitions — a plain
// `outgoing.par_iter()` over `n_parts` buckets gives one thread the
// whole of one partition's bucket, however large, and a lot of
// buckets are far from equal: a central-base substitution changes
// the minimiser (and thus routes to a different partition) only
// when the winning minimiser window overlaps the central position.
// For k=31/m=11 that is ~11 of the 21 possible windows — so *most*
// of the remaining ~10/21 send the variant right back to the
// partition already being built. That self-partition bucket ends up
// far larger than any other, so the naive per-partition split
// leaves one thread grinding through it alone long after every
// other partition's (tiny) bucket is done — confirmed by sampling a
// real run: one thread solely in `MphfLayer::find` while the rest
// of the pool sits idle. Splitting each non-empty bucket into
// `chunk_size`-sized pieces first keeps this pass's per-partition
// locality (each chunk is still contiguous within one partition,
// still resolved via `cache.find` grouped by that partition) while
// letting Rayon spread a single oversized bucket across several
// threads instead of pinning it to one. ──────────────────────────
let chunk_size = ((outgoing.iter().map(Vec::len).sum::<usize>() / n_workers).max(1)).min(4096);
let work: Vec<(usize, &[(CanonicalKmer, usize, u8)])> = outgoing
.iter()
.enumerate()
.filter(|(_, q)| !q.is_empty())
.flat_map(|(dest, q)| q.chunks(chunk_size).map(move |c| (dest, c)))
.collect();
work.par_iter().for_each(|&(dest, chunk)| {
for &(variant, source_order, base) in chunk {
if cache.find(dest, variant) {
mask[source_order].fetch_or(1 << base, Ordering::Relaxed);
}
}
});
check_empty("after_resolution");
// ── Write the layer's annex file, indexed by iteration order — a
// second streamed pass over `unitigs.bin` (via `enumerate_kmers`),
// now that every entry's mask is final; never a `Vec` hold of the
+40
View File
@@ -258,3 +258,43 @@ fn family_scan_consumers_agree_on_one_sibling_each() {
assert_eq!(total, 2, "no other base pair should be tallied");
assert_eq!(base_pairs.same, [0, 0, 0, 0], "the two genomes never agree at this locus");
}
#[test]
fn sibling_annex_no_empty_masks_after_build() {
let dir = tempdir().unwrap();
let seq = b"ACGTACGTACGT".repeat(500);
let g1 = build_single_genome_index(dir.path(), "g1", &seq);
g1.build_sibling_annex().expect("build_sibling_annex");
let index_dir = g1.partition().part_dir(0).join(INDEX_SUBDIR);
let meta = PartitionMeta::load(&index_dir).expect("partition meta");
for l in 0..meta.n_layers {
let layer_dir = index_dir.join(format!("layer_{l}"));
let annex = SiblingAnnex::open(&layer_dir.join(ANNEX_FILE_NAME)).expect("annex open");
for slot in 0..annex.len() {
let mask = annex.get(slot).expect("slot must have an entry");
assert!(
mask.family_size() >= 1,
"layer {l} slot {slot}: empty family mask (bits={:04b}) after build — \
indicates batch-offset bug or unseeded mask slot",
mask.bits()
);
}
}
}
#[test]
fn sibling_histogram_does_not_panic_on_partial_last_batch() {
// Regression test for the bug where enumerate_kmers_batch().enumerate()
// used the batch index instead of the sequence offset, leaving the
// trailing slots of a non-multiple-of-BATCH_SIZE layer unseeded and
// producing the minorant-only sentinel (0x10) that caused
// sibling_family_size_histogram to panic with index u32::MAX.
let dir = tempdir().unwrap();
let seq = b"ACGTACGTACGT".repeat(500);
let g1 = build_single_genome_index(dir.path(), "g1", &seq);
g1.build_sibling_annex().expect("build_sibling_annex");
let hist = g1.sibling_family_size_histogram().expect("sibling_family_size_histogram");
let total: u64 = hist.iter().sum();
assert!(total > 0, "histogram should contain at least one family");
}