obikmer/docmd/architecture/siblings.md

Sibling annex — architecture (discussion)

Status: architecture decided (2026-08-14). Implementation not yet mandated.

Two index spaces, uncorrelated

Every kmer stored in a Layer lives in two independent index spaces:

• Iteration order: its position when enumerating unitigs.bin (the superkmer file), deterministic but arbitrary with respect to slot.
• MPHF slot: MphfLayer::index(kmer), the number the MPHF assigns.

The two are not correlated by any formula. Converting from one to the other requires either recomputing the MPHF (kmer → slot) or scanning the iteration stream (kmer → order). There is no slot → kmer operation: the MPHF is a one-way function, not an invertible bijection with a stored inverse. Any method that reconstructs a kmer from a bare slot number is wrong by construction, regardless of the mechanism used (MPHF re-hash, or evidence decode + direct unitig read). See MphfLayer::kmer_at (obilayeredmap/src/mphf_layer.rs) — flagged for removal, currently called from obikphylo/siblings/build.rs and family_scan.rs (since removed — see "Pending work" status below).

Two pipelines, never mixed

	origin of the kmer	membership known?	correct mapping
query pipeline	external (caller-supplied)	no	query/find/find_strict — MPHF + evidence check
iteration pipeline	enumerated from this layer's own unitigs.bin	yes, by construction	index/index_batch — MPHF only, no evidence

Evidence exists solely to answer "is this external kmer a member of the layer" for the query pipeline. Using it (or the MPHF) to go the other way — recover a kmer from a slot, or re-verify a kmer that was just produced by iterating the layer — is a conceptual error: evidence can be probabilistic (Approx mode), so any slot→kmer attempt is unsound in general, and pointless even in Exact/Hybrid mode since the kmer was already known.

Sibling annex: an iteration-pipeline artifact only

The sibling annex (FamilyMask/SiblingAnnex, .psib, obicompactvec/src/siblingannex.rs) records, per kmer, whether it is a family minorant and which family members are present in the index. Its only consumers (obikphylo/siblings/stats.rs, family_scan.rs) enumerate it exhaustively (0..annex.len()); no query-pipeline code path touches it.

Decision: the annex must be persisted in iteration order, not slot order. This lets readers zip-iterate Layer::iter_kmers() and the annex file directly — one linear, cache-friendly pass, no MPHF/slot indirection, no kmer_at. It also enables specialized iterators building on this zip: minorants-only iteration, batch-of-kmers → batch-of-family-members, etc.

Today the annex is built and stored in slot order (build_layer_sibling_annex, siblings/build.rs): slot_kmer is populated via (0..n_slots).map(|slot| mphf.kmer_at(slot)), and the origin slot is threaded through the whole cross-partition reconciliation pipeline (variant generation, query_partition_with, final mask[slot].fetch_or(...)). This must change to iterating iter_kmers()/enumerate_kmers() and threading the iteration index instead of the slot end to end — eliminating kmer_at from the build path entirely, not just the read path. No slot-indexed intermediate is needed even during construction; the iteration-order id is sufficient throughout.

The cross-partition side of the same pipeline is unaffected: checking whether a generated family-variant kmer exists in another partition is a genuine query-pipeline operation (the variant's membership in the target partition is unknown) and must keep going through KmerPartition::query_partition_with (MPHF + evidence), never a raw index().

Pending work — done

The plan above shipped: obikphylo (a new crate — phylo-domain extension traits over obikindex::KmerIndex/obilayeredmap::Layer<D>, replacing the old obikindex::siblings module) builds and reads the annex purely in iteration order (SiblingLayerExt::iter_siblings/iter_minorants, both with batch variants, mirroring Layer<D>'s own KmerIter/KmerBatchIter shape). MphfLayer::kmer_at has no remaining callers.

A separate, unrelated bug surfaced during this work and was fixed (2026-08-14): MphfLayer::enumerate_kmers_batch computed its batch_start_index via the stdlib .enumerate() adapter, which counts batches (0, 1, 2…), not the cumulative k-mer offset the annex is actually keyed on — every batch past the first wrote its mask/annex entries at the wrong iteration-order position. Fixed by tracking a running offset instead; regression tests added (sibling_annex_no_empty_masks_after_build, sibling_histogram_does_not_panic_on_partial_last_batch).

Performance: build_sibling_annex parallelism (2026-08-14)

Investigated on a real multi-genome run (phyloskims_sal_vac, k=31/m=11). Baseline: mostly one active core, with short multi-core bursts — average ~3 cores.

Fixes that helped, kept:

• CanonicalKmerOf::minimizer() (obikseq/src/kmer.rs) — a direct O(k) bit-arithmetic minimiser for a single isolated k-mer, replacing a RollingStat instance fed byte-by-byte through an ASCII round-trip (used by helpers::partition_of, called for every generated family variant). ~3x wall-clock improvement on its own, confirmed by sampling (obiskbuilder::rolling_stat/obikentropy frames disappeared from the hot path). CanonicalKmerOf::partition() added alongside it (wraps minimizer().seq_hash() & mask, the same routing rule KmerPartition/RoutableSuperKmer use).
• Cross-partition resolution (outgoing.par_iter() in build_layer_sibling_annex) parallelised at the partition level — one Rayon task per non-empty outgoing[dest] bucket. For k=31/m=11, a central-base substitution changes the winning minimiser (and thus the destination partition) only when that window overlaps the central base: ~11 of the 21 possible windows do, so ~10/21 (≈48%) of generated variants route right back to the partition already being built. That self bucket ends up far larger than any other, so the per-partition split pinned one thread to it alone while the rest of the pool finished instantly — confirmed by sampling: one thread solid in MphfLayer::find, everyone else idle. Fixed by splitting each non-empty bucket into total_queries / n_workers (capped 4096) chunks before par_iter(), preserving per-partition mmap locality (each chunk stays contiguous within one partition) while letting Rayon spread an oversized bucket across several threads. Net effect of both fixes together: ~3 cores average → ~10-13 cores average on the same run, and a projected total build time of ~1h15 down to ~30min on the real phyloskims_sal_vac run this was measured against.
• TracedBar's ETA (obisys/src/progress.rs) was silently starved: the custom progress message and the self-computed ETA text used to share one pb.set_message() slot, with the ETA holding off for 2s after any custom message — fine when custom messages are rare, broken once build_sibling_annex's per-partition callback fires more often than that. Fixed by keeping the two texts in separate fields, composed together on every render instead of one overwriting the other.

Tried and reverted — do not repeat blindly:

• Parallelising the outer partition loop in build_sibling_annex with obikindex::PartitionRunner (already used by merge/build_layers), splitting a fixed core budget between outer (partition) and inner (pipeline + resolution) concurrency so their product wouldn't exceed the budget. Measured worse: throughput dropped over time (26 partitions/5min → 38/11-12min) and peak resolution concurrency fell from ~11-12 cores to ~7-8. Cause: this capped the resolution burst — which scales very well on its own — to make room for outer concurrency, and running several partitions' resolution at once scatters access across multiple partitions' mmap regions at once, working against the locality outgoing's per-partition grouping exists for. PartitionRunner stayed exported from obikindex (new_capped too) since it's general-purpose, but nothing in obikphylo calls it.
• Splitting resolution chunks even finer (/(n_workers*8), cap 1024, instead of /n_workers, cap 4096) to smooth the residual sawtooth. Measured ~10% slower, wider dips, not narrower. Reverted to the original chunk sizing.

Cross-partition batch resolution — current state vs. the batched-accumulator design (discussion, 2026-08-14)

family_scan.rs::scan_layer_families (shared by snp_pseudo_alignment, sibling_annex_stats, cardinality_tally, scan_family_pairs) already implements most of a dispatch/accumulate/resolve pipeline: generation (cheap, CPU-only — builds outgoing[dest_partition] from FamilyMask and buckets cross-partition queries) runs on an obipipeline::throttle + make_pipe! stage, decoupled from resolution (I/O-bound, rayon::par_iter across partitions, one generated batch resolved at a time, never several concurrently — this ordering is deliberate, see the module's own docs on a reverted concurrent-batch-resolution attempt that scattered mmap access). The fast/slow mode gate (PartitionCache::fast_mode, cache.rs:162-163) already exists: n_layers <= 7 (checked once from the first non-empty partition's PartitionMeta::n_layers, documented as identical across every partition of an index — a structural, build-time property, never a per-partition state) decides whether FamilyMask's recorded layer_value can be trusted to skip straight to the right layer (find_presence_batch_fast) or must fall back to scanning every layer of the destination partition (find_presence_batch).

Real gap, confirmed not implemented: resolution is triggered by the source batch finishing (FAMILY_BATCH = 65536 minorants read from the scanned layer), not by an output accumulator filling up. Since most central-base variants of a family route back to the same partition being scanned (~48% per the k=31/m=11 measurement above), a FAMILY_BATCH's outgoing[dest] is large for the local/self partition and thin for the other ~255 (or however many) destination partitions — each of those gets resolved at low query density every batch instead of being accumulated across several source batches until resolving it is worthwhile. This is distinct from, and not fixed by, the fast/slow layer gate above.

Redesign sketched (not built): per-destination accumulators decoupled from FAMILY_BATCH, flushed on reaching a size threshold instead of on source-batch completion — a "hot" accumulator for the partition being scanned (sharded one-per-generation-worker, no lock, since all n_workers pipeline workers write to it concurrently — this differs from an earlier, simpler mental model of "one thread owns one layer's local collector," which doesn't hold here since n_workers threads cooperate on scanning one layer at a time, not one thread per layer) and "cold" mutex-per-partition accumulators for the rest, low contention expected since traffic to any single cold destination is a small fraction of total.

This breaks the current strict-iteration-order delivery of on_family (today: a reorder buffer keyed by batch, since a whole FAMILY_BATCH resolves atomically). With cross-batch accumulation, a family only becomes complete once every accumulator holding one of its outgoing queries has flushed, at unpredictable, independent times — no longer streamable strictly in order without a large, unbounded pending buffer. Resolution sketched: replace order-dependent consumers with coordinate-addressed writes instead of order-dependent appends (see PseudoAlignment idea below) wherever possible, since sibling_annex_stats's reduction (plain counts) is already order-independent and needs nothing here.

Pseudo-alignment at scale — pruning is unavoidable (discussion, 2026-08-14)

The reference run (phyloskims_sal_vac-scale bacterial test set, iqtree.fasta) produced a dense alignment for 13 genomes × 383,965 sites (4.8 MB) — trivially small. The in-progress plant index is expected to carry on the order of 9 billion minorant families; a dense byte-per-cell alignment at that column count is unbuildable regardless of genome count (hundreds of GB even at a handful of genomes). Long-term ambition is 6,000– 8,000 genomes on a large machine, which makes the per-cell cost dominant in the other dimension too. Pruning the retained family set before materializing anything is mandatory, not an optimization.

Already free: family_size() < 2 (no sibling variant registered at all) is a zero-cost structural filter, read directly off FamilyMask bits, already applied in snp_pseudo_alignment. Insufficient alone — per the ~80-85% mono-family estimate from earlier discussion, this only brings 9 billion down to roughly 1.3-1.8 billion, still unusable.

Criterion under discussion: fix a global cell budget (n_genomes × n_columns_retained ≤ threshold) and retain the highest-entropy families first until the budget is spent — column count self-adjusts to genome count automatically. Entropy of a family: Shannon entropy over the observed base distribution across genomes carrying that family, decided 2026-08-14: genomes where the family is absent are excluded from the calculation (denominator = genomes where present, not all genomes) — measures signal purity where it exists, independent of coverage rate.

Key consequence for storage: entropy needs per-variant genome counts, which cannot be derived from FamilyMask's bits alone — it requires the same cross-partition resolution work the alignment/stats pipeline already pays for. This favors computing it once at annex-build time and persisting it in an auxiliary vector alongside the annex (one value per retained minorant) over encoding a fixed keep/discard decision as one of FamilyMask's 3 remaining free bits (bits 13-15, siblingannex.rs:77): the bit approach bakes a single threshold in permanently (any different budget needs a full annex rebuild), the vector approach pays the expensive resolution once and lets budget/threshold be chosen freely per analysis afterward.

Open, explicitly deferred: scope of the top-K selection — global across the whole index (needs a first pass computing/storing every family's entropy, then a global threshold from the full distribution, e.g. a quantile, before building the final alignment: two full passes, but the size budget is honored precisely) vs. local per layer/partition (stays within the current single-pass streaming model, but the global budget is no longer exactly guaranteed — depends on how unevenly entropy is distributed across layers).

Known remaining limitation, not yet worth fixing: within one layer, the four stages (sequential unitigs.bin read → parallel generation → parallel resolution → sequential annex write) never overlap — confirmed by 1s-interval sampling: generation alone occupies ~17 threads evenly, but the next layer's read/generation never starts until the current layer's resolution and write are both done. This produces a real, periodic (~layer duration) alternation between "many cores" and "few cores" that neither of the fixes above touches, since both operate within one layer's resolution step. The only remaining lever is overlapping consecutive layers (e.g. a depth-2 pipeline: start layer N+1's read/generation while layer N's resolution/write is still running) — a real restructuring, not a parameter tweak, and explicitly not to be combined with the reverted budget-capping idea above (let each phase use however many cores it naturally wants; only the scheduling needs to overlap). Deferred, not started.