PartitionedIndex
-├── LayeredPartition (one per minimiser bucket)
-│ ├── MphfLayer 0 kmer → slot (immutable bijection)
-│ │ ├── DataStore A slot → T (e.g. counts)
-│ │ └── DataStore B slot → T (e.g. presence/absence, derived)
-│ ├── MphfLayer 1
-│ │ └── DataStore A
-│ └── ...
-├── LayeredPartition
-│ └── ...
+KmerIndex (index.meta + KmerPartition)
+├── partition_0/index/ one directory per minimiser bucket
+│ ├── meta.json PartitionMeta { n_layers }
+│ ├── layer_0/
+│ │ ├── layer_meta.json LayerMeta { evidence: EvidenceKind }
+│ │ ├── mphf.bin PtrHash MPHF
+│ │ ├── unitigs.bin unitig spine (never overwritten)
+│ │ ├── evidence.bin exact evidence (Exact only)
+│ │ ├── unitigs.bin.idx block index (Exact only)
+│ │ ├── fingerprint.bin fingerprints (Approx only)
+│ │ ├── counts/ PersistentCompactIntMatrix (with_counts = true)
+│ │ └── presence/ PersistentBitMatrix
+│ └── layer_1/
+│ └── ...
+└── partition_1/index/
+ └── ...
PartitionedIndex : routes queries to partitions via canonical minimiser hash. Owns the partition count and routing scheme (fixed at creation). Dispatches aggregations across partitions in parallel.
-LayeredPartition : one directory per minimiser bucket. Holds a Vec<MphfLayer>. Each layer covers a disjoint kmer set — layer 0 is built from dataset A; layer 1 covers kmers in B absent from layer 0; and so on. Layers within a partition are always disjoint.
-MphfLayer : the MPHF + evidence + unitig spine. Maps kmer → slot for its disjoint kmer set. Immutable once built. Independent of any data attached to it.
-DataStore : a slot-indexed data array (e.g. PersistentCompactIntMatrix, PersistentBitMatrix). Attached to a MphfLayer externally. Multiple stores of different types can coexist on the same MphfLayer.
+KmerIndex : root entry point. Owns IndexMeta (written to index.meta) and a KmerPartition that routes canonical kmers to partition directories. All partition-level operations are dispatched in parallel via rayon.
+Partition directory : one directory per minimiser bucket. PartitionMeta (stored as meta.json) records n_layers. Layers within a partition cover disjoint kmer sets.
+Layer directory : one MphfLayer plus optional data stores. LayerMeta (stored as layer_meta.json) records which EvidenceKind was used. The MPHF and unitigs.bin are immutable once built; evidence files are the only part replaced by reindex.
MphfLayer — autonomous mapping layer
-MphfLayer :: find ( kmer : CanonicalKmer ) -> Option < usize > // slot, or None if absent
-MphfLayer :: n () -> usize // number of slots
-MphfLayer :: build ( dir : & Path ) -> OLMResult < ( Self , usize ) > // from unitigs.bin
-MphfLayer :: open ( dir : & Path ) -> OLMResult < Self >
-find returns Some(slot) only if the kmer is actually in this layer (evidence check included). Returns None for kmers present in other layers or absent from the index.
-The MPHF (mphf.bin, evidence.bin, unitigs.bin) is built once and never rebuilt. All data derivation operations (count → presence, thresholding, merging) reuse the same MphfLayer.
-DataStore — slot-indexed data
-trait DataStore {
- type Item ;
- fn get ( & self , slot : usize ) -> Self :: Item ;
- fn n ( & self ) -> usize ;
+
+pub struct IndexConfig {
+ pub kmer_size : usize ,
+ pub minimizer_size : usize ,
+ pub n_bits : usize , // log2(n_partitions)
+ pub with_counts : bool ,
+ pub evidence : EvidenceKind ,
+ pub block_bits : u8 , // .idx granularity: 2^block_bits unitigs/block; 0 = one entry per unitig
+}
+
+pub struct IndexMeta {
+ pub version : u32 ,
+ pub config : IndexConfig ,
+ pub genomes : Vec < GenomeInfo > , // ordered; index = genome column number
+}
+
+pub struct GenomeInfo {
+ pub label : String ,
+ pub meta : HashMap < String , String > , // arbitrary categorical metadata
}
Concrete types from obicompactvec:
+IndexMeta is serialised as index.meta (JSON). It is the authority for the ordered list of genomes and for the parameters that govern all subsequent operations on the index.
+EvidenceKind
+pub enum EvidenceKind {
+ Exact ,
+ Approx { b : u8 , z : u8 },
+}
+Controls which files are written per layer and which query path is taken:
+
+
+
+Variant
+Files written
+False-positive rate
+
+
+
+
+Exactevidence.bin, unitigs.bin.idx0
+
+
+Approx { b, z }fingerprint.bin≈ W / 2^(b·z) per read (Findere)
+
+
+
+EvidenceKind is stored both in IndexConfig (index-wide default, updated by reindex) and in each LayerMeta (per-layer record of what was actually built).
+MphfLayer — autonomous kmer → slot mapping
+pub struct MphfLayer {
+ mphf : PtrHash < … > ,
+ ev : LayerEvidence , // Exact { evidence, unitigs } | Approx { fingerprint }
+ n : usize ,
+}
+MphfLayer::find(kmer) dispatches transparently to find_exact or find_approx based on the evidence loaded at open time (read from layer_meta.json). Returns Some(slot) only if the kmer is confirmed present; None for absent or out-of-range.
+find_exact: slot = mphf(kmer); decode evidence → (chunk_id, rank); verify kmer in unitigs
+find_approx: slot = mphf(kmer); check fingerprint[slot] == seq_hash(kmer)
+block_bits controls the .idx file written alongside evidence.bin. At block_bits = 0, every unitig chunk has an index entry, giving O(1) random access; larger values trade access time for a smaller .idx.
+The MPHF and unitigs.bin are never rebuilt by any post-build operation.
+Layer\<D> — MPHF + data payload
+pub struct Layer < D : LayerData = () > {
+ mphf : MphfLayer ,
+ data : D ,
+}
+D selects the attached data payload:
+
+
+
+DData directory
+Item returned by query
+
+
+
+()—
+() (set membership only)
+
+PersistentCompactIntMatrixcounts/Box<[u32]> (counts per genome)
+
+PersistentBitMatrixpresence/Box<[bool]> (presence per genome)
+
+
+Layer::query(kmer) delegates to MphfLayer::find, then calls data.read(slot) if a slot is returned. Both exact and approximate evidence are handled transparently; the caller sees only Option<Hit<D::Item>>.
+Build-time entry points:
+Layer < () > :: build ( out_dir , block_bits ) // set membership
+Layer < PersistentCompactIntMatrix > :: build ( out_dir , block_bits , count_of )
+Layer < PersistentBitMatrix > :: build_presence ( out_dir , block_bits , n_genomes , present_in )
+
+Layer :: < () > :: build_evidence ( layer_dir , kind , block_bits ) // evidence only (reindex path)
+DataStore — slot-indexed data
+PersistentCompactIntMatrix and PersistentBitMatrix are slot-indexed stores. They know nothing about kmers or MPHFs.
Type
ItemColumn stats
+Aggregation method
Use
@@ -1471,120 +1552,46 @@
PersistentCompactIntMatrixBox<[u32]>sum() -> Array1<u64>count per sample per slot
+sum() → Array1<u64>counts per genome per slot
PersistentBitMatrixBox<[bool]>count_ones() -> Array1<u64>presence per sample per slot
+count_ones() → Array1<u64>presence per genome per slot
-sum() and count_ones() are the bridge between the per-matrix level and cross-layer aggregation: they give the total weight of each column within one (partition, layer) pair, which can be summed to get global column weights.
-A DataStore knows nothing about kmers or MPHFs. It is indexed by usize slot only.
Distance matrix API on DataStore types
-Both PersistentCompactIntMatrix and PersistentBitMatrix expose two families of distance matrix methods.
-Full distance matrices
-Compute the final n_cols × n_cols distance matrix from data within a single matrix. Internally parallelised over the upper triangle via rayon.
-// PersistentCompactIntMatrix
-fn bray_dist_matrix ( & self ) -> Array2 < f64 >
-fn relfreq_bray_dist_matrix ( & self ) -> Array2 < f64 >
-fn euclidean_dist_matrix ( & self ) -> Array2 < f64 >
-fn relfreq_euclidean_dist_matrix ( & self ) -> Array2 < f64 >
-fn hellinger_dist_matrix ( & self ) -> Array2 < f64 >
-fn jaccard_dist_matrix ( & self ) -> Array2 < f64 >
-fn threshold_jaccard_dist_matrix ( & self , threshold : u32 ) -> Array2 < f64 >
-
-// PersistentBitMatrix
-fn jaccard_dist_matrix ( & self ) -> Array2 < f64 >
-fn hamming_dist_matrix ( & self ) -> Array2 < u64 >
-These are convenience methods. For a LayeredDataStore or PartitionedDataStore they cannot be used directly — the partial API is required.
-Partial distance matrices
-Return additive components that can be summed element-wise across (partition, layer) pairs before computing the final distance. This is what makes cross-layer and cross-partition aggregation possible.
-Category 1 — self-contained partials : additive without any external parameter.
-// PersistentCompactIntMatrix
-fn partial_bray_dist_matrix ( & self )
- -> ( Array2 < u64 > , // sum_min[i,j]
- Array1 < u64 > ) // col_sums[k]
-
-fn partial_euclidean_dist_matrix ( & self ) -> Array2 < f64 > // sum of squared diffs
-fn partial_threshold_jaccard_dist_matrix ( & self , threshold : u32 )
- -> ( Array2 < u64 > , // inter[i,j]
- Array2 < u64 > ) // union[i,j]
-
-// PersistentBitMatrix
-fn partial_jaccard_dist_matrix ( & self )
- -> ( Array2 < u64 > , // inter[i,j]
- Array2 < u64 > ) // union[i,j]
-fn partial_hamming_dist_matrix ( & self ) -> Array2 < u64 > // differing bits
-Category 2 — normalised partials : require global column sums as input, computed beforehand across all (partition, layer) pairs.
-// PersistentCompactIntMatrix only
-fn partial_relfreq_bray_dist_matrix ( & self , col_sums : & Array1 < u64 > )
- -> Array2 < f64 > // Σ_slot min(a_slot/sum_i, b_slot/sum_j)
-
-fn partial_relfreq_euclidean_dist_matrix ( & self , col_sums : & Array1 < u64 > )
- -> Array2 < f64 > // Σ_slot (a_slot/sum_i - b_slot/sum_j)²
-
-fn partial_hellinger_euclidean_dist_matrix ( & self , col_sums : & Array1 < u64 > )
- -> Array2 < f64 > // Σ_slot (√(a/sum_i) - √(b/sum_j))²
-The col_sums parameter must reflect the GLOBAL count across all layers and all partitions — passing a per-layer sum would give a wrong result. This constraint drives the two-pass algorithm described below.
-Progressive aggregation principle
-Aggregation is hierarchical : each level computes its contribution by aggregating from the level immediately below it. No level skips a level or collects raw data from two levels down.
-PersistentCompactIntMatrix::col_weights() — column sums for one (partition, layer) matrix
- ↓ Σ across layers
-LayeredStore<PersistentCompactIntMatrix>::col_weights() — column sums for one partition
- ↓ Σ across partitions
-LayeredStore<LayeredStore<…>>::col_weights() — global column sums
-The same cascade applies to every partial:
-PersistentCompactIntMatrix::partial_bray() — one (partition, layer)
- ↓ element-wise Σ across layers
-LayeredStore<PersistentCompactIntMatrix>::partial_bray() — one partition
- ↓ element-wise Σ across partitions
-LayeredStore<LayeredStore<…>>::partial_bray() — global partial → final dist
-Each level presents a stable trait surface to the level above; no level reaches two levels down.
-Traits — obicompactvec::traits
-Three traits unify the aggregation API across all levels of the hierarchy.
+Aggregation traits — obicompactvec::traits
+Three traits unify the aggregation API across all hierarchy levels.
trait ColumnWeights : Send + Sync {
fn col_weights ( & self ) -> Array1 < u64 > ;
}
trait CountPartials : ColumnWeights {
- // self-contained partials (additive, no parameter)
- fn partial_bray ( & self ) -> Array2 < u64 > ;
- fn partial_euclidean ( & self ) -> Array2 < f64 > ;
- fn partial_threshold_jaccard ( & self , threshold : u32 ) -> ( Array2 < u64 > , Array2 < u64 > );
- // normalised partials (global col_weights passed in cascade)
- fn partial_relfreq_bray ( & self , global : & Array1 < u64 > ) -> Array2 < f64 > ;
- fn partial_relfreq_euclidean ( & self , global : & Array1 < u64 > ) -> Array2 < f64 > ;
- fn partial_hellinger ( & self , global : & Array1 < u64 > ) -> Array2 < f64 > ;
- // provided finalisation methods (default implementations)
- fn bray_dist_matrix ( & self ) -> Array2 < f64 > { … }
- fn euclidean_dist_matrix ( & self ) -> Array2 < f64 > { … }
- fn threshold_jaccard_dist_matrix ( & self , threshold : u32 ) -> Array2 < f64 > { … }
- fn relfreq_bray_dist_matrix ( & self ) -> Array2 < f64 > { … }
- fn relfreq_euclidean_dist_matrix ( & self ) -> Array2 < f64 > { … }
- fn hellinger_dist_matrix ( & self ) -> Array2 < f64 > { … }
+ fn partial_bray ( & self ) -> Array2 < u64 > ;
+ fn partial_euclidean ( & self ) -> Array2 < f64 > ;
+ fn partial_threshold_jaccard ( & self , threshold : u32 ) -> ( Array2 < u64 > , Array2 < u64 > );
+ fn partial_relfreq_bray ( & self , global : & Array1 < u64 > ) -> Array2 < f64 > ;
+ fn partial_relfreq_euclidean ( & self , global : & Array1 < u64 > ) -> Array2 < f64 > ;
+ fn partial_hellinger ( & self , global : & Array1 < u64 > ) -> Array2 < f64 > ;
+ // provided finalisation methods with default impls
+ fn bray_dist_matrix ( & self ) -> Array2 < f64 > { … }
+ fn relfreq_bray_dist_matrix ( & self ) -> Array2 < f64 > { … }
+ // …
}
trait BitPartials : ColumnWeights {
- fn partial_jaccard ( & self ) -> ( Array2 < u64 > , Array2 < u64 > );
- fn partial_hamming ( & self ) -> Array2 < u64 > ;
+ fn partial_jaccard ( & self ) -> ( Array2 < u64 > , Array2 < u64 > );
+ fn partial_hamming ( & self ) -> Array2 < u64 > ;
// provided
fn jaccard_dist_matrix ( & self ) -> Array2 < f64 > { … }
fn hamming_dist_matrix ( & self ) -> Array2 < u64 > { … }
}
Leaf implementors (in obicompactvec):
+Leaf implementors:
@@ -1595,40 +1602,73 @@ LayeredStore<LayeredStore<…>>::partial_bray() — glo
PersistentCompactIntMatrixColumnWeights (via sum()), CountPartialsColumnWeights, CountPartials
PersistentBitMatrixColumnWeights (via count_ones()), BitPartialsColumnWeights, BitPartials
-PersistentCompactIntVec and PersistentBitVec do not implement these traits — they are single-column primitives, not matrix-level aggregators.
LayeredStore<S> — obilayeredmapA single generic wrapper replaces the need for named LayeredDataStore and PartitionedDataStore types:
+LayeredStore\<S> — recursive aggregation wrapper
pub struct LayeredStore < S > ( Vec < S > );
Three blanket impls propagate the traits up the hierarchy:
-impl < S : ColumnWeights > ColumnWeights for LayeredStore < S > { … } // Σ across inner stores
-impl < S : CountPartials > CountPartials for LayeredStore < S > { … } // same pattern
-impl < S : BitPartials > BitPartials for LayeredStore < S > { … } // same pattern
+Three blanket impls propagate all traits up the hierarchy:
+impl < S : ColumnWeights > ColumnWeights for LayeredStore < S > { … }
+impl < S : CountPartials > CountPartials for LayeredStore < S > { … }
+impl < S : BitPartials > BitPartials for LayeredStore < S > { … }
Because the blanket impl is recursive, LayeredStore<LayeredStore<S>>S does — no separate PartitionedStore type is needed:
-PersistentCompactIntMatrix implements CountPartials
-LayeredStore<PersistentCompactIntMatrix> via blanket impl (= one partition)
-LayeredStore<LayeredStore<…>> via blanket impl (= partitioned index)
+This makes LayeredStore<LayeredStore<PersistentCompactIntMatrix>> automatically implement CountPartials — no separate PartitionedStore type is needed:
+PersistentCompactIntMatrix leaf (one layer)
+LayeredStore<PersistentCompactIntMatrix> one partition (layers are disjoint)
+LayeredStore<LayeredStore<…>> whole index (partitions are independent)
Normalised metrics — two-pass cascade
-The normalised finalisation methods call col_weights() first (pass 1), then the normalised partial (pass 2). Both calls go through the same blanket impl, so the cascade is automatic:
-// called on LayeredStore<LayeredStore<PersistentCompactIntMatrix>>
+Normalised metrics require global column sums — computed in a two-pass cascade:
+// on LayeredStore<LayeredStore<PersistentCompactIntMatrix>>
fn relfreq_bray_dist_matrix ( & self ) -> Array2 < f64 > {
- let global = self . col_weights (); // pass 1 — progressive sum at every level
- let p = self . partial_relfreq_bray ( & global ); // pass 2 — global passed in cascade
- p . mapv ( | v | 1.0 - v ) // finalise (diagonal zeroed separately)
+ let global = self . col_weights (); // pass 1 — sums up hierarchy
+ let p = self . partial_relfreq_bray ( & global ); // pass 2 — global broadcast read-only
+ p . mapv ( | v | 1.0 - v )
}
global is exact: each kmer belongs to exactly one (partition, layer) pair, so there is no double-counting across the hierarchy.
+Because each kmer belongs to exactly one (partition, layer) pair, col_weights() has no double-counting across the hierarchy.
+Progressive aggregation principle
+No level reaches two levels down. Each level sums contributions from the level immediately below:
+PersistentCompactIntMatrix::col_weights() — one (partition, layer)
+ ↓ Σ across layers
+LayeredStore<PersistentCompactIntMatrix>::col_weights() — one partition
+ ↓ Σ across partitions
+LayeredStore<LayeredStore<…>>::col_weights() — global
+The same cascade applies to every partial method.
+Multi-genome column invariant
+After any merge, every layer in every partition has exactly n_genomes columns, where n_genomes is the current total in index.meta. This holds for both PersistentCompactIntMatrix and PersistentBitMatrix.
+Maintained by three coordinated operations:
+Existing layers — column append. Layer::append_genome_column appends one column to each existing layer. Slots matching the incoming genome receive its count or true; all other slots receive 0 or false.
+New layers — absent columns prepended. When a new layer is created for kmers unique to the incoming genome, n_existing_genomes absent columns are prepended before the incoming genome's column, so the new layer immediately has the same column count as all other layers.
+First merge, Presence mode — init_presence_matrix. The initial single-genome index has no presence/ directory (presence is implicit). On the first merge, Layer<()>::init_presence_matrix materialises genome 0's presence column (all true) retroactively, raising the column count from 0 to 1 before appending column 1.
+This invariant is the precondition for correct progressive aggregation: every level can blindly sum matrices from below because all matrices have the same shape.
+Query model
+Point query
+minimiser(kmer) → partition p
+for each layer l in p:
+ if let Some(slot) = MphfLayer_l.find(kmer):
+ return data_l.read(slot)
+return None
+O(n_layers) MPHF probes worst case; O(1) expected. The result comes from exactly one (partition, layer).
+Aggregation
+result = reduce(
+ for p in partitions: // parallel
+ for l in layers(p): // parallel
+ partial(data_p_l)
+)
+For normalised metrics, replace with the two-pass cascade.
Parallelism model
@@ -1642,12 +1682,12 @@ LayeredStore<LayeredStore<…>> via blanket impl (= par
Across partitions
-LayeredStore<LayeredStore<S>> inner storesnone — fully independent
+inner stores of LayeredStore<LayeredStore<S>>
+none
Across layers within a partition
-LayeredStore<S> inner storesinner stores of LayeredStore<S>
none — disjoint kmer sets
@@ -1667,55 +1707,17 @@ LayeredStore<LayeredStore<…>> via blanket impl (= par
-All levels use rayon par_iter internally; reduce_with performs a parallel tree reduction.
Query model
-Point query — kmer → Option<Item>
-minimiser(kmer) → partition p
-for each layer l in p:
- slot = MphfLayer_l.find(kmer)
- if slot is Some:
- return DataStore_l.get(slot)
-return None
-O(n_layers) MPHF probes worst case; O(1) expected. No cross-layer fusion — the result comes from exactly one (partition, layer).
-Aggregation — → Result
-result = reduce(
- for p in partitions: // parallel
- for l in layers(p): // parallel
- partial(DataStore_p_l)
-)
-For normalised metrics replace with the two-pass scheme above.
-DataStore derivation
-Because the MphfLayer is independent of its data stores, new stores can be derived from existing ones without rebuilding the MPHF:
-// count → presence/absence, parallel across (partition, layer)
-for (p, l) in all_partition_layer_pairs().par_iter():
- count_store = open PersistentCompactIntMatrix at (p, l)
- presence_store = PersistentBitMatrix::from_count_matrix(count_store, threshold, dir)
-Other derivations: threshold a count matrix → binary presence matrix; union two presence matrices; merge two count matrices (saturating add, column-wise). All are local to one (partition, layer) pair.
-Relationship to current implementation
-What is implemented
+reindex — evidence conversion in place
+KmerIndex::reindex(target, block_bits) converts every layer's evidence bundle to target without touching the MPHF or unitigs.bin:
-obicompactvec::traitsColumnWeights, CountPartials, BitPartials are defined and implemented on PersistentCompactIntMatrix and PersistentBitMatrix.obilayeredmap::LayeredStore<S>LayeredStore<LayeredStore<S>> is the partitioned level — no separate type needed. Tests confirm that splitting data across layers and across partitions gives the same distance matrices as computing on flat combined data.→ Exact: builds evidence.bin + unitigs.bin.idx; removes fingerprint.bin→ Approx { b, z }: builds fingerprint.bin; removes evidence.bin + unitigs.bin.idx
-What is not yet implemented
-
-Layer<D: LayerData> still fuses MphfLayer and one DataStore. Multiple data stores on the same MPHF are not supported.LayeredMap is a single-partition structure without distance matrix API; it does not yet use LayeredStore.No PartitionedIndex type for point queries with parallel partition dispatch.
-
-Planned refactoring
-
-Extract MphfLayer from Layer<D> as an autonomous type.
-Replace LayerData trait with the DataStore / ColumnWeights / CountPartials / BitPartials system.
-Rewire LayeredMap to hold LayeredStore<PersistentCompactIntMatrix> (or bit variant) alongside the MPHF layers.
-Implement PartitionedIndex using LayeredStore<LayeredStore<S>> for data and parallel dispatch for queries.
-
+On success, IndexConfig::evidence and IndexConfig::block_bits are updated in index.meta. Each layer's layer_meta.json is also rewritten with the new EvidenceKind.
+estimate — parameter dry-run
+estimate resolves approximate-evidence parameters (z, b, target FP rate) and prints the resulting effective kmer size and per-kmer / per-z-window false-positive rates without touching any index. Used to calibrate Approx { b, z } before building or reindexing.
diff --git a/doc/architecture/query.refs/index.html b/doc/architecture/query.refs/index.html
new file mode 100644
index 0000000..ab56d5a
--- /dev/null
+++ b/doc/architecture/query.refs/index.html
@@ -0,0 +1,1062 @@
+
+
+
+
+
+ Query.refs - obikmer
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Coverage: architecture/query.md
+Code couvert
+
+obikmer/src/cmd/query.rs — commande query, format de sortieobikpartitionner/src/query_layer.rs — routage de la requête à travers les partitionsobiread/src/lib.rs — lecture des séquences d'entrée pour la requête
+Notes
+RISQUE DE DÉRIVE. Vérifier :
+- La commande unitig a été modifiée pour utiliser open_sequential() — vérifier si query est concerné
+- find_exact / find_approx / find générique ont été ajoutés dans MphfLayer — le chemin de requête a changé
+- Si l'index est approximatif (Approx), la requête peut produire des faux positifs : la doc le mentionne-t-elle ?
+- Format de sortie CSV (obikindex/src/csv.rs ou équivalent) à vérifier
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Query system - obikmer
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Query system
+Goal
+Given a set of query sequences, determine for each sequence how many of its k-mers are found in the index and, for each indexed genome, how many k-mers match. The query system is the foundation for read classification and sequence-to-genome mapping.
+
+Query sequences in FASTA or FASTQ format (gzip supported, streaming stdin supported). GenBank flat files are not supported at query time (only at index time).
+Sequences shorter than k bases are silently skipped.
+Non-ACGT characters are handled by the superkmer decomposition layer: they act as hard breaks, producing shorter superkmers (identical to the behaviour at indexing time).
+
+Algorithm
+The query follows the same superkmer-based partitioning strategy used at indexing time.
+for each chunk of sequences (parallel workers via obipipeline):
+ build QueryBatch: decompose all sequences into s-mers via superkmers, deduplicate
+ allocate seq_results[seq_idx][smer_pos] = None ← per-sequence s-mer result vectors
+ split superkmers by partition via minimiser hash
+ for each partition p:
+ query_partition(p, superkmers_routed_to_p)
+ → load QueryLayer(s) for p
+ → for each s-mer in each superkmer: MphfLayer::find(smer)
+ fill seq_results[seq_idx][kmer_offset + j] from partition results
+ for each sequence:
+ apply_findere(seq_results[seq_idx], effective_z) ← per full sequence
+ accumulate confirmed k-mer results into acc and cov
+ emit annotated sequences
+Superkmers that appear more than once in the batch (same sequence or across sequences) are deduplicated: each unique RoutableSuperKmer is queried once per partition, and the result is broadcast to every SKDesc entry that references it.
+Findere requires full-sequence aggregation. apply_findere is applied once per sequence on the complete s-mer result vector, after all partitions have contributed. Applying it per superkmer would produce false negatives at superkmer boundaries, where the z-window spans two superkmers.
+Batches are processed in parallel via obipipeline workers; the --threads flag controls the number of worker threads.
+Findere z-window filter
+For approximate index modes, the index physically stores s-mers of size s = k_user − z + 1. At query time, set_k(s) is in effect, so queries naturally produce s-mer results. apply_findere then aggregates z consecutive s-mer results into one k_user-mer answer:
+fn apply_findere (
+ results : & [ Option < Box < [ u32 ] >> ], // N s-mer results
+ z : usize ,
+ n_genomes : usize ,
+) -> Vec < Option < Box < [ u32 ] >>> // N − z + 1 k_user-mer results
+Input length N (s-mers), output length N − z + 1 (k_user-mers).
+For each genome g independently, a sliding window of size z scans the input. Output position i is confirmed for genome g iff all z values results[i..i+z][g] are nonzero (None counts as zero for all genomes). The scan is O(n) per genome.
+Output values come from results[i] (leftmost s-mer of each window); genomes not confirmed are zeroed. If all genomes are zero, the position is returned as None.
+Short sequences : when the s-mer count is less than z, no complete window can form — apply_findere returns an empty vector. K-mers from sequences shorter than k_user are not emitted.
+Exact indexes : z = 1, apply_findere is a passthrough (output length = input length).
+Effective z at query time
+effective_z is resolved at the start of run():
+let effective_z = args . findere_z . unwrap_or_else ( || match idx . meta (). config . evidence {
+ IndexMode :: Approx { z , .. } | IndexMode :: Hybrid { z , .. } => z as usize ,
+ IndexMode :: Exact => 1 ,
+});
+The -z CLI option overrides the index metadata value. A higher z increases stringency (lower FP, some true positives may be discarded at sequence ends); a lower z increases sensitivity.
+Layer lookup: MphfLayer::find
+MphfLayer::open(dir, mode: &IndexMode) receives the mode from PartitionMeta — no per-layer file is read. The caller (QueryLayer) never chooses the dispatch path: it is fixed at open time by LayerEvidence. See obilayeredmap for the full find / find_strict API.
+QueryLayer variant selectionQueryLayer::open in query_layer.rs selects the data matrix to pair with MphfLayer:
+
+
+
+Condition
+Variant
+Data returned per k-mer
+
+
+
+
+with_counts=true and counts/ existsCountraw count per genome
+
+
+presence/ existsPresence0/1 per genome (bit matrix)
+
+
+only counts/ exists
+Countcounts used as-is
+
+
+neither exists
+SetOnly1 for every genome
+
+
+
+Presence / count mode at query time
+The --force-presence flag and --presence-threshold control how per-genome values are accumulated, independently of what the index stores:
+genome_totals[g] += if presence { u32::from(v >= threshold) } else { v }
+presence is true when --force-presence is set or when the index has no counts (!with_counts). The default presence_threshold is 1, so any nonzero count counts as a match.
+Coverage vectors (--detail)
+When --detail is requested, a 3-D accumulator cov[seq_idx][genome][kmer_pos] is allocated after all partitions are queried, with dimensions derived from n_kmers_out = n_smers − z + 1 (k_user-mer positions, not s-mer positions):
+cov[seq_idx][g][pos] += contribution
+where pos is the k_user-mer index in the filtered (post-Findere) vector
+Coverage reflects confirmed k_user-mers only. The vectors are emitted in the JSON annotation under the key "coverage".
+kmer_missing semanticskmer_missing counts k_user-mer positions where the first s-mer (seq_results[seq_idx][pos]) is None — i.e. absent from the index entirely. K-mers where the z-window fails because a later s-mer is absent or zero are not counted as missing (the first s-mer being present is used as proxy for index membership).
+Output sequences are written in OBITools4 format : the original sequence with a JSON annotation map in the title line.
+>read_id {"kmer_count":59,"kmer_strict_matches":{"genome_a":42,"genome_b":7}}
+ATCGATCG...
+With --detail:
+>read_id {"kmer_count":59,"kmer_strict_matches":{...},"coverage":{"genome_a":[0,1,2,...],...}}
+ATCGATCG...
+Genome keys follow the iteration order of meta.genomes.
+Annotation schema
+
+
+
+Key
+Type
+Condition
+Semantics
+
+
+
+
+kmer_countint
+always
+k-mers confirmed (post-Findere) with at least one genome match
+
+
+kmer_missingint
+--count-missingk-mers absent from the index entirely (pre-Findere None)
+
+
+kmer_strict_matchesobject
+always
+per-genome accumulated value (label → count or 0/1)
+
+
+coverageobject
+--detailper-genome array of per-position contributions (label → [u32])
+
+
+
+kmer_count + kmer_missing ≤ total k_user-mers in the sequence. The gap corresponds to k_user-mers whose z-window was not fully confirmed (at least one s-mer absent or zero for all genomes) but whose first s-mer was present in the index.
+CLI
+obikmer query <index> [--detail] [--mismatch] [--count-missing]
+ [--force-presence] [--presence-threshold <n>]
+ [-z <z>] [-T <threads>]
+ <query.fa> [<query2.fa> ...]
+
+
+
+Option
+Default
+Semantics
+
+
+
+
+-z / --findere-zfrom index metadata
+Override Findere z parameter
+
+
+--detailoff
+Emit per-position coverage vectors in JSON
+
+
+--count-missingoff
+Add kmer_missing field to JSON
+
+
+--force-presenceoff
+Report 0/1 per genome regardless of index counts
+
+
+--presence-threshold1
+Minimum count to declare genome present
+
+
+-T / --threadsall CPUs
+Worker threads
+
+
+
+--mismatch is accepted but currently ignored with a warning on stderr.
+Future work
+
+--mismatch3·k single-substitution variants per k-mer, look each up independently.Read classification (--classify): assign each read to the genome with the highest match score.Whitelist / blacklist filtering : threshold-based accept/reject on per-genome match scores.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Invariant.refs - obikmer
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Coverage: architecture/sequences/invariant.md
+Code couvert
+
+obikseq/src/sequence.rs — invariants de représentation des séquences (ACGT, longueur max)obikseq/src/unitig.rs — type Unitig, contrainte MAX_KMERS_PER_CHUNK (255 kmers par chunk)
+Notes
+Document court et stable. Vérifier que la limite de 256 nucléotides (ou 255 kmers) par chunk
+est toujours la même dans obiskio::MAX_KMERS_PER_CHUNK.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Evidence elimination (discussion)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
@@ -725,6 +753,62 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Merge command
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Kmer filtering (rebuild/dump/unitig)
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/doc/implementation/chunkreader.refs/index.html b/doc/implementation/chunkreader.refs/index.html
new file mode 100644
index 0000000..ec8a30e
--- /dev/null
+++ b/doc/implementation/chunkreader.refs/index.html
@@ -0,0 +1,1058 @@
+
+
+
+
+
+ Chunkreader.refs - obikmer
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Coverage: implementation/chunkreader.md
+Code couvert
+
+obiread/src/chunk.rs — SeqChunkIter, détection de frontières FASTA/FASTQ, state machinesobikrope/src/lib.rs — type Rope (Vec), opérations zero-copy
+Notes
+Document stable (la stratégie de chunking rope ne devrait pas avoir changé).
+Vérifier que le split FASTA/FASTQ reste correct si de nouveaux formats ont été ajoutés.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Evidence elimination (discussion)
+
+
+
@@ -383,6 +404,30 @@
+
+
+
+
+
+
+
+
+ Merge command
+
+
+
+
+
+
+
+
+
+
+
+ Kmer filtering (rebuild/dump/unitig)
+
+
+
@@ -454,19 +499,28 @@
-
+
- Output type: rope
+ Two reading paths
-
+
- Allocation policy
+ Record path: chunk reader
+
+
+
+
+
+
+
+
+ Output type: Rope
@@ -506,68 +560,77 @@
Chunk reader — implementation
-The obiread crate provides a streaming iterator that reads FASTA or FASTQ files in fixed-size blocks and yields self-contained chunks, each ending on a complete sequence record boundary. Chunks are consumed in parallel by downstream workers.
-Output type: rope
-Each chunk is a Vec<Bytes> — a rope : a list of reference-counted byte slices that are not necessarily contiguous in memory. The consumer iterates over the slices in order.
-Using bytes::Bytes means the split at the record boundary is O(1): Bytes::split_to(n) adjusts a reference counter, not memory. No memcpy in the common case.
-Allocation policy
+obiread exposes two distinct sequence reading paths, each optimised for a different use case.
+Two reading paths
-Case
-Cost
+Path
+API
+Output unit
+Per-record identity
+Use case
-Boundary found in the current block (common)
-zero extra allocation — split_to only
+Record path read_sequence_chunks → parse_chunkSeqRecord (id + raw sequence + normalised rope)yes
+query — must read complete records
-Boundary straddles multiple blocks (sequence > block size, rare)
-one allocation to pack the rope into a flat buffer
-
-
-EOF flush
-zero extra allocation
+Stream path open_nuc_streamNucPage (flat normalised byte buffer)no
+index, superkmer — bulk throughput
+The record path uses Rope-backed chunks and is described in detail below.
+The stream path (NucStream / NucPage) is described in the scatter section of pipeline .
+Record path: chunk reader
+The chunk reader reads FASTA or FASTQ files in fixed-size blocks and yields self-contained chunks, each ending on a complete sequence record boundary. parse_chunk then converts each chunk into a Vec<SeqRecord>, where each record carries its identifier, raw sequence bytes, and a normalised rope ready for superkmer building.
+This path is mandatory for query, where superkmers must be tracked back to their originating sequence (id, kmer offset) for output annotation.
+Output type: Rope
+Each chunk is a Rope — a segmented byte sequence: a Vec of blocks, where each block is a Vec<Cell<u8>>. The consumer iterates over the blocks via a forward or backward cursor.
+Rope::split_off(pos) splits at an absolute byte offset in O(log n) (binary search over block-start index). If pos falls inside a block, that block is split in two via Vec::split_off — no memcpy in the common case.
SeqChunkIter
pub struct SeqChunkIter < R : Read > { /* private */ }
impl < R : Read > Iterator for SeqChunkIter < R > {
- type Item = io :: Result < Vec < Bytes >> ;
+ type Item = io :: Result < Rope > ;
}
pub fn fasta_chunks < R : Read > ( source : R ) -> SeqChunkIter < R >
pub fn fastq_chunks < R : Read > ( source : R ) -> SeqChunkIter < R >
next() loop:
-1. read one block of block_size bytes → push onto rope
-2. probe check: if the boundary marker ("\n>" or "\n@") is absent from the
- last block, skip the splitter (avoids a full backward scan for nothing)
-3. call splitter on last block
- if found at offset n:
- remainder = last_block.split_to(n) ← O(1), zero copy
- return std::mem::take(&mut self.rope) ← the chunk
-4. if rope.len() > 1 (multi-block accumulation):
- pack rope into one flat buffer ← one alloc
- retry splitter on flat buffer
-5. if EOF: flush remaining rope as final chunk
+1. read one block of block_size bytes → push onto Rope
+2. call splitter(rope) → Option<abs_offset>
+ if Some(pos):
+ tail = rope.split_off(pos) ← O(log n), may split one block
+ chunk = mem::replace(&mut rope, tail)
+ return Some(Ok(chunk))
+3. if EOF and rope non-empty: return Some(Ok(rope)) as final chunk
+4. if EOF and rope empty: return None
The Splitter function signature is fn(&Rope) -> Option<usize>. It returns the absolute byte offset of the start of the last complete record, or None if no boundary was found in the accumulated rope (need more data).
Boundary detection — FASTA
-Backward scan with a 2-state machine. Searches for > immediately preceded by \n or \r:
+Backward scan with a 2-state machine. Searches (right to left) for > followed by \n or \r (i.e., a > that is preceded by a newline in forward order):
stateDiagram-v2
direction LR
[*] --> Scanning
Scanning --> FoundGt : '>'
FoundGt --> Scanning : other
FoundGt --> [*] : '\\n' / '\\r' ✓Returns the byte offset of the > that starts the last complete record.
+Returns the byte offset of the > that starts the last complete record. Returns None if only one > is found (cannot confirm there is a prior complete record).
Boundary detection — FASTQ
FASTQ records have a rigid 4-line structure (@header, sequence, +, quality). The @ character (ASCII 64, Phred score 31) can appear legitimately in quality lines, making any forward heuristic unreliable. The backward scanner verifies the full structural context before accepting a candidate @.
-7-state machine (port of Go's EndOfLastFastqEntry), scanning from right to left . Each time a + is found, its position is saved as restart; any state mismatch resets the scan to that position.
+7-state machine (states 0–6), scanning from right to left . Each time a + is found, its position is saved as restart; any state mismatch resets the scan to that position.
stateDiagram-v2
direction LR
diff --git a/doc/implementation/evidence_elimination.refs/index.html b/doc/implementation/evidence_elimination.refs/index.html
new file mode 100644
index 0000000..565edff
--- /dev/null
+++ b/doc/implementation/evidence_elimination.refs/index.html
@@ -0,0 +1,1068 @@
+
+
+
+
+
+ Evidence elimination.refs - obikmer
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Coverage: implementation/evidence_elimination.md
+Code couvert
+
+obilayeredmap/src/fingerprint.rs — FingerprintVec, FingerprintVecWriter, stockage b bits/slot, matches()obilayeredmap/src/mphf_layer.rs — build_approx_evidence(dir, b, z), find_approx()obilayeredmap/src/meta.rs — EvidenceKind::Approx { b, z }, LayerMetaobikindex/src/reindex.rs — KmerIndex::reindex(), conversion exact↔approx en placeobikmer/src/cmd/reindex.rs — CLI reindex, options --approx, -z, --evidence-bits, --fp, --block-sizeobikmer/src/cmd/index.rs — resolve_approx_params(), options --approx, -z, --evidence-bits, --fpobikmer/src/cmd/estimate.rs — commande estimate (dry-run des paramètres)
+Notes
+Ce document était à l'origine une discussion de design (4 approches). L'implémentation
+a maintenant convergé vers l'approche fingerprint (Findere-style).
+FORT RISQUE DE DÉRIVE — le contenu est probablement un mélange de design et d'implémentation :
+- Le modèle FP = 1/2^(b·z) et les règles de résolution (2-of-3 parmi b, z, fp) sont implémentés
+- La commande reindex permet la conversion a posteriori exact↔approx
+- La commande estimate fait le dry-run des paramètres
+Cette page doit être réécrite pour documenter l'implémentation Findere réelle plutôt que les alternatives abandonnées.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Evidence elimination (discussion) - obikmer
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Approximate evidence: fingerprint-based index
+Motivation
+evidence.bin maps each MPHF slot to the position of the k-mer that owns it,
+enabling zero-FP verification. On the bacterial BCT dataset (2048 partitions,
+k=31, ~33 M k-mers/partition) it accounts for 66 % of the lookup-layer footprint:
+
+
+
+file
+size/partition
+fraction
+
+
+
+
+evidence.bin
+132 MB
+66 %
+
+
+unitigs.bin
+58 MB
+29 %
+
+
+mphf.bin
+10 MB
+5 %
+
+
+
+evidence.bin is a bijection from MPHF-space to unitig-position-space and
+costs at minimum ⌈log₂ N⌉ bits per slot — an information-theoretic floor with
+only ~22 % packing headroom. Compression is not a path to elimination.
+The approximate index replaces evidence.bin + unitigs.bin.idx with a
+fingerprint.bin file. The MPHF and unitigs.bin are kept unchanged. Set
+operations still require an exact index; the approximate index targets query
+workloads that can tolerate a bounded false-positive rate.
+The Findere model
+A B-bit fingerprint stored per MPHF slot provides the discrimination that
+evidence.bin would otherwise provide through full k-mer reconstruction.
+For a foreign k-mer query, the MPHF maps it to some slot s. The fingerprint
+stored at s belongs to the legitimate k-mer at that slot. The FP event is:
+P(FP per k-mer) = 1 / 2^b
+The Findere trick reduces the indexed k-mer size. When the user specifies k_user
+and z, the index physically stores k-mers of size s = k_user − z + 1. At query
+time, the same s-mer size is used. After collecting per-position s-mer results
+over the full query sequence, a sliding window of size z aggregates z consecutive
+s-mer hits into one confirmed k_user-mer hit, reducing the per-window FP rate:
+P(FP per k_user-mer) = 1 / 2^(b·z)
+IndexConfig::kmer_size stores s = k_user − z + 1, not k_user. Both indexing
+and querying use this stored size via set_k(idx.kmer_size()).
+Parameters b and z are stored in layer_meta.json (EvidenceKind::Approx { b, z }).
+FingerprintVec on diskfingerprint.bin layout:
+magic: b"FPVF" (4 bytes)
+b: u8 (bits per slot, 1..=64)
+padding: [0u8; 3]
+n: u64 LE (number of slots)
+data: packed bits, ceil(n·b/8) bytes, Lsb0 order
+FingerprintVec is memory-mapped. The match check against a query k-mer:
+fn matches ( & self , slot : usize , fingerprint : u64 ) -> bool {
+ self . get ( slot ) == ( fingerprint & self . mask )
+}
+build_approx_evidence iterates unitigs.bin sequentially, writes
+kmer.seq_hash() into the slot assigned by the MPHF, then saves fingerprint.bin
+and layer_meta.json. No .idx file is produced; random access into
+unitigs.bin is not needed.
+At build time, find_approx in MphfLayer:
+let slot = self . mphf . index ( & kmer . raw ());
+if fingerprint . matches ( slot , kmer . seq_hash ()) { Some ( slot ) } else { None }
+layer_meta.json records which evidence bundle is present:
+pub enum EvidenceKind {
+ Exact ,
+ Approx { b : u8 , z : u8 },
+}
+MphfLayer::open reads this tag and dispatches find to find_exact or
+find_approx transparently. find_exact panics on an approximate layer;
+find_approx panics on an exact layer — mode mixing is a programming error.
+Parameter resolution (resolve_approx_params)
+The identity b·z = ⌈−log₂(fp)⌉ lets any two of (b, z, fp) derive the third.
+resolve_approx_params implements a 2-of-3 rule with conservative ceiling
+rounding:
+
+
+
+given
+derived
+
+
+
+
+b, z
+fp = 1/2^(b·z)
+
+
+z, fp
+b = ⌈−log₂(fp) / z⌉
+
+
+b, fp
+z = ⌈−log₂(fp) / b⌉
+
+
+z only
+b = 8 (default), fp derived
+
+
+b only
+z = 1 (default), fp derived
+
+
+fp only
+b = 8 (default), z derived
+
+
+none
+b = 8, z = 1, fp = 1/256
+
+
+
+When all three are given, b and z are authoritative and fp is recomputed.
+CLI flags
+Both index and reindex accept the same flags:
+
+
+
+flag
+type
+meaning
+
+
+
+
+--approxbool
+enable fingerprint evidence
+
+
+--evidence-bits (b)u8
+fingerprint bits per slot
+
+
+-zu8
+Findere z parameter
+
+
+--fpf64
+target FP rate per z-window
+
+
+--block-sizeusize
+unitig block size for exact .idx; ignored in approx mode
+
+
+
+--approx must be set explicitly; the other three flags are optional and
+resolved by the 2-of-3 rule. Omitting all three produces b=8, z=1.
+reindex commandreindex converts an existing index between exact and approximate evidence
+in-place across all partitions and layers, running partitions in parallel via
+Rayon.
+Conversion to approximate (--approx):
+
+Builds fingerprint.bin from unitigs.bin + mphf.bin.
+Removes evidence.bin and unitigs.bin.idx.
+Updates layer_meta.json with EvidenceKind::Approx { b, z }.
+
+Conversion to exact (default, no --approx):
+
+Builds evidence.bin + unitigs.bin.idx from unitigs.bin + mphf.bin.
+Removes fingerprint.bin.
+Updates layer_meta.json with EvidenceKind::Exact.
+
+The root index.meta is updated with the new evidence kind on success.
+mphf.bin and unitigs.bin are never modified.
+estimate commandestimate is a dry-run that resolves and prints (b, z, fp) without touching
+any index. It accepts the same --evidence-bits, -z, and --fp flags and
+additionally accepts -k to display the effective indexed k-mer length:
+k (user): 31
+k (indexed, s=k-z+1): 27
+z: 5
+evidence bits (b): 8
+FP per s-mer: 3.906e-3 (1/2^8)
+FP per k-mer window: 9.537e-7 (1/2^(8·5))
+Useful for choosing parameters before committing to an index build.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Kmer.refs - obikmer
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Coverage: implementation/kmer.md
+Code couvert
+
+obikseq/src/kmer.rs — layout mémoire (repr(transparent) u64), encodage/décodage, revcomp, forme canoniqueobikseq/src/params.rs — k global (set_k / k())
+Notes
+Document d'implémentation stable. L'algorithme de revcomp bit-à-bit est décrit —
+vérifier qu'il correspond à revcomp_raw dans obiskio/src/unitig_index.rs (copie locale)
+et à l'implémentation dans obikseq/src/kmer.rs.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
+
- Memory layout
+ Types and layout
+
+
+
+
+
+
+
+
+
+
+ Global parameters
@@ -1017,10 +1145,32 @@
-
+
- Canonical form
+ Canonical form and CanonicalKmerOf
+
+
+
+
+
+
+
+
+
+
+ Sliding window helpers
+
+
+
+
+
+
+
+
+
+
+ Hashing
@@ -1045,12 +1195,43 @@
Kmer — implementation
-Memory layout
-Kmer is a #[repr(transparent)] newtype over u64:
+Types and layout
+KmerOf<L> is a #[repr(transparent)] newtype over u64 parameterized by a KmerLength marker:
#[repr(transparent)]
-pub struct Kmer ( u64 );
+pub struct KmerOf < L : KmerLength > ( u64 , PhantomData < L > );
Nucleotides are packed 2 bits each, left-aligned , MSB-first. Nucleotide 0 occupies bits 63–62; nucleotide i occupies bits 63−2i and 62−2i. The low 64−2k bits are always zero. k is not stored — it is a parameter of every operation that needs it, and will be owned by the future collection-level indexer.
+Three marker types implement KmerLength:
+
+
+
+Marker
+len() sourceUsed for
+
+
+
+
+KLenparams::k()k-mers
+
+
+MLenparams::m()minimizers
+
+
+ConstLen<N>const generic N
+tests
+
+
+
+Public aliases:
+pub type Kmer = KmerOf < KLen > ; // k-mer, global k
+pub type Minimizer = CanonicalKmerOf < MLen > ; // canonical m-mer, global m
+Nucleotides are packed 2 bits each, left-aligned , MSB-first. Nucleotide 0 occupies bits 63–62; nucleotide i occupies bits 63−2i and 62−2i. The low 64−2·len bits are always zero. The length is not stored — every operation reads it from L::len().
+Global parameters
+params::set_k(k) / params::k() and params::set_m(m) / params::m() are backed by OnceLock<usize> in production (write-once, panic on conflict) and by thread_local! { Cell<usize> } in test builds (per-thread, freely writable). params::init(k, m) sets both in one call.
Encoding
-Kmer::from_ascii(ascii, k) encodes the first k bytes of an ASCII slice using the shared ENC table (see SuperKmer — ASCII encoding ):
+KmerOf::<L>::from_ascii(ascii) encodes the first L::len() bytes using the shared ENC table (see SuperKmer — ASCII encoding ):
for i in 0 .. k {
val = ( val << 2 ) | encode_base ( ascii [ i ]) as u64 ;
}
-Kmer ( val << ( 64 - 2 * k ))
+KmerOf ( val << ( 64 - 2 * k ), PhantomData )
Zero allocation — result lives on the stack.
Decoding
-write_ascii(k, buf) appends k ASCII characters to a caller-supplied Vec<u8> using the shared DEC4 table: one lookup per 4 nucleotides, two partial-byte lookups for the remainder. No allocation in the hot path.
-to_ascii(k) is a convenience wrapper that allocates and returns a Vec<u8>; intended for tests and display only.
+write_ascii(writer) writes k ASCII characters to any W: Write using the shared DEC4 table: one lookup per 4 nucleotides, one partial lookup for the remainder. No allocation in the hot path.
+to_ascii() is a convenience wrapper that allocates and returns a Vec<u8>; intended for tests and display only.
Reverse complement
Computed as pure arithmetic — no lookup table, no memory access:
let x = ! self . 0 ; // complement
let x = x . swap_bytes (); // reverse bytes
let x = (( x >> 4 ) & 0x0F0F0F0F0F0F0F0F ) | (( x & 0x0F0F0F0F0F0F0F0F ) << 4 ); // swap nibbles
let x = (( x >> 2 ) & 0x3333333333333333 ) | (( x & 0x3333333333333333 ) << 2 ); // swap 2-bit groups
-Kmer ( x << ( 64 - 2 * k ))
+KmerOf ( x << ( 64 - 2 * k ), PhantomData )
After complementing, bytes are reversed (swap_bytes), then nibbles, then 2-bit groups — restoring 2-bit nucleotides to their correct positions in reverse order. A final left-shift realigns to MSB. Zero allocation — result lives on the stack.
-
-pub fn canonical ( & self , k : usize ) -> Self {
- let rc = self . revcomp ( k );
- if self . 0 <= rc . 0 { * self } else { rc }
+
+canonical() returns a CanonicalKmerOf<L> — a distinct newtype that carries the same u64 layout but enforces the invariant that the stored value equals min(kmer, revcomp):
+pub fn canonical ( & self ) -> CanonicalKmerOf < L > {
+ let rc = self . revcomp ();
+ CanonicalKmerOf ( if self . 0 <= rc . 0 { self . 0 } else { rc . 0 }, PhantomData )
}
Lexicographic minimum of forward and reverse-complement, comparing the raw u64 values directly (left-aligned encoding makes this equivalent to nucleotide-wise comparison). Zero allocation — result lives on the stack.
+CanonicalKmerOf::from_raw_unchecked(raw) is the only other public constructor, for trusted paths such as deserialisation.
+Sliding window helpers
+push_right(nuc) / push_left(nuc) shift the window by one base in O(1). is_overlapping(other) checks whether the last k−1 nucleotides of self equal the first k−1 of other.
+Hashing
+hash_kmer(raw: u64) -> u64 computes mix64(raw ^ 0x9e3779b97f4a7c15), the seeded splitmix64 finalizer. CanonicalKmerOf::seq_hash() delegates to hash_kmer.
diff --git a/doc/implementation/merge.refs/index.html b/doc/implementation/merge.refs/index.html
new file mode 100644
index 0000000..ff3be53
--- /dev/null
+++ b/doc/implementation/merge.refs/index.html
@@ -0,0 +1,1065 @@
+
+
+
+
+
+ Merge.refs - obikmer
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Coverage: implementation/merge.md
+Code couvert
+
+obikindex/src/merge.rs — KmerIndex::merge(), validation de compatibilité d'évidence, validate_evidence_compat()obikpartitionner/src/merge_layer.rs — merge_partition(), construction de la nouvelle layer, paramètre block_bitsobikpartitionner/src/rebuild_layer.rs — rebuild_partition(), paramètre block_bitsobilayeredmap/src/layer.rs — Layer::append_genome_column() (PersistentCompactIntMatrix et PersistentBitMatrix)obicompactvec/src/intmatrix.rs — append_column pour PersistentCompactIntMatrixobicompactvec/src/bitmatrix.rs — append_column pour PersistentBitMatrix
+Notes
+FORT RISQUE DE DÉRIVE. Changements récents :
+- Ajout de la validation de compatibilité d'évidence : merge exact+approx → erreur (OKIError::IncompatibleEvidence)
+- merge_partition reçoit maintenant block_bits: u8
+- La commande reindex a été ajoutée comme outil de conversion exact↔approx avant merge
+Vérifier que la doc décrit la politique de merge mixed-evidence et le recours à reindex.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Merge command - obikmer
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Merge command
+Purpose
+obikmer merge combines multiple existing kmer indexes into a single index. The result contains all kmers from all sources, with per-genome presence/absence or count data for every genome across every layer.
+Modes
+pub enum MergeMode { Presence , Count }
+Default mode is Presence. Count mode requires all source indexes to have with_counts=true; mixing count and non-count sources is rejected at validation.
+
+
+
+Mode
+Column type
+Constraint
+
+
+
+
+PresencePersistentBitMatrix (one bit per genome per slot)none
+
+
+CountPersistentCompactIntMatrix (one u32 per genome per slot)all sources with_counts=true
+
+
+
+All source indexes must satisfy:
+
+IndexState::Indexed (fully built — index.done sentinel present)Same kmer_size, minimizer_size, n_partitions
+Same evidence kind: all Exact, or all Approx with identical (b, z) parameters
+If Count mode: all sources must have with_counts=true
+
+--force: if the output directory already exists, it is deleted before the merge begins.
+Evidence compatibility
+validate_evidence_compat(sources) is called before any I/O. It compares each source's EvidenceKind against sources[0]:
+
+All Exact → accepted, output uses Exact
+All Approx { b, z } with same (b, z) → accepted, output uses those parameters
+Any other combination → OKIError::IncompatibleEvidence, with a message directing the user to run reindex first
+
+Mixed exact/approx is a hard error, not a silent conversion.
+fn validate_evidence_compat ( sources : & [ & KmerIndex ]) -> OKIResult < EvidenceKind >
+Genome label deduplication
+compute_labels(sources, rename_duplicates) assigns final genome labels across all sources before any file is written. The first occurrence of a label keeps the original name. Subsequent occurrences receive .1, .2, … suffixes when rename_duplicates is true, or trigger OKIError::DuplicateGenomeLabel otherwise.
+Algorithm
+1. Validation
+Check all sources against the constraints above. Abort on any mismatch.
+2. Bootstrap output from first source
+Recursive file copy of sources[0] → output. Immediately after the copy:
+
+index.meta is rewritten with the final genome list (all sources, possibly renamed) and the effective evidence kind.In Presence mode, any counts/ directories inherited from source_0 are removed.
+spectrums/ from source_0 is removed and rebuilt from scratch across all sources, applying the (possibly renamed) labels.
+This establishes the partition layout, all existing MPHFs, unitigs, and evidence files. The first source's genomes occupy columns 0 … n_dst_genomes - 1 in the destination.
+3. For each subsequent source (parallel across partitions)
+KmerPartition::merge_partition(i, sources, mode, n_dst_genomes, block_bits) is called for each partition index i. block_bits is taken from dst.meta.config.block_bits.
+Each entry in sources is (&KmerPartition, n_genomes) where n_genomes is the column count that source contributes (> 1 when the source is itself a merged index).
+First merge, Presence mode : when n_dst_genomes == 1, Layer::<()>::init_presence_matrix is called on every existing destination layer before any source column is appended. This creates presence/col_000000.pbiv set all-true (genome 0 is present in every slot).
+Pass 1 — classify kmers
+Iterate all kmers from all source partitions (via UnitigFileReader + canonical kmer iteration). For each kmer, probe the destination LayeredMap<()>:
+
+Hit : kmer already in the destination; record for Pass 2.Miss : push kmer into a GraphDeBruijn accumulator.
+New layer construction
+If the accumulator is non-empty, compute de Bruijn unitigs and call Layer::<()>::build(&new_layer_dir, block_bits). All kmers absent from the destination — across all sources — accumulate into a single graph, producing one new layer per merge operation (not one per source).
+Pass 2 — fill column builders
+For each source and each of its layers, re-iterate unitigs and look up stored values via SrcLayerData::lookup(kmer, src_n):
+
+SrcLayerData::SetMembership — no data directory exists; every kmer returns vec![1; n_genomes]SrcLayerData::Presence — reads PersistentBitMatrix from presence/SrcLayerData::Count — reads PersistentCompactIntMatrix from counts/
+Hits are routed to exist_builders[dst_layer][src_col]; misses are routed to new_src_builders[src_col].
+Column prepending for new layers
+Before source columns are written to the new layer, n_dst_genomes absent columns (all-zero / all-false) are prepended — one per genome already in the index — so the column count invariant holds immediately after layer creation.
+Close and update metadata
+Close all builders; update presence/meta.json or counts/meta.json with {"n": N, "n_cols": n_dst_genomes + n_src_total}; increment PartitionMeta::n_layers if a new layer was added.
+
+index.meta was already written during bootstrap with the complete genome list and evidence kind. No further update is needed after the partition loop.
+append_genome_columnDefined on two concrete specialisations of Layer<D>:
+impl Layer < PersistentCompactIntMatrix > {
+ pub fn append_genome_column ( layer_dir : & Path , value_of : impl Fn ( usize ) -> u32 ) -> OLMResult < () >
+}
+
+impl Layer < PersistentBitMatrix > {
+ pub fn append_genome_column ( layer_dir : & Path , value_of : impl Fn ( usize ) -> bool ) -> OLMResult < () >
+}
+Each appends one column file to the matrix subdirectory (counts/ or presence/). In merge_partition, columns are written directly via PersistentBitVecBuilder / PersistentCompactIntVecBuilder rather than through these helpers, but the invariant they enforce is the same.
+Column count invariant
+After any merge, every layer in every partition has exactly n_genomes columns , where n_genomes is the total genome count in the index at that point.
+Maintained by three mechanisms:
+
+Existing layers : n_src_total columns appended (one per source genome).New layers created during merge : n_dst_genomes absent columns prepended before source columns.First merge, Presence mode : init_presence_matrix retroactively creates presence/col_0 all-true for genome 0.
+The invariant is a precondition of LayeredStore aggregation traits: col_weights() and all partial distance methods assume every inner store has the same column count.
+Error variants relevant to merge
+
+
+
+Variant
+Condition
+
+
+
+
+OKIError::NotIndexed(path)Source not in Indexed state
+
+
+OKIError::IncompatibleConfigMismatched kmer_size, minimizer_size, or n_partitions
+
+
+OKIError::MismatchedModeCount mode but a source has with_counts=false
+
+
+OKIError::IncompatibleEvidence(msg)Mixed exact/approx or different approx (b, z)
+
+
+OKIError::DuplicateGenomeLabel(label)Duplicate label and rename_duplicates=false
+
+
+
+On-disk impact
+After merging G genomes (sources_0 contributes G0, subsequent sources the rest):
+partitions/
+ part_00000/
+ index/
+ meta.json ← n_layers updated if new layer added
+ layer_0/
+ mphf.bin ← unchanged
+ unitigs.bin ← unchanged
+ evidence.bin ← unchanged
+ presence/ ← created on first merge (Presence mode)
+ meta.json {"n": N, "n_cols": G}
+ col_000000.pbiv ← all-true (genome 0 … G0-1)
+ col_000001.pbiv ← next source
+ ...
+ counts/ ← extended (Count mode)
+ meta.json {"n": N, "n_cols": G}
+ col_000000.pciv ← genome 0 counts (from original build)
+ col_000001.pciv ← next source
+ ...
+ layer_N/ ← new layer (if new kmers found)
+ mphf.bin
+ unitigs.bin
+ evidence.bin
+ presence/ or counts/
+ meta.json {"n": N1, "n_cols": G}
+ col_000000.pbiv ← all-false (absent for existing genomes)
+ ...
+spectrums/
+ <label>.json ← one file per genome, rebuilt from all sources
+index.meta ← complete genome list + evidence kind written at bootstrap
+
+
+
+
+
+
+
+
+
+
Mphf.refs - obikmer
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Coverage: implementation/mphf.md
+Code couvert
+
+obilayeredmap/src/mphf_layer.rs — type Mphf (PtrHash + CubicEps + CachelineEfVec + Xx64), construction en 2 passes, build(), build_exact_evidence(), build_approx_evidence(), build_evidence()obikpartitionner/src/index_layer.rs — build_index_layer() avec passage de block_bits
+Notes
+FORT RISQUE DE DÉRIVE. Changements récents :
+- build_exact_evidence(dir, block_bits) — block_bits maintenant paramétrisé (défaut 0)
+- build_approx_evidence(dir, b, z) — nouvelle fonction pour l'évidence fingerprint
+- build_evidence(dir, kind, block_bits) — dispatch selon EvidenceKind
+- Construction en 2 phases : pass 1 (Rayon parallèle) + pass 2 (callback fill_slot)
+Vérifier que la doc décrit correctement les deux nouvelles routes d'évidence et le paramètre block_bits.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Evidence modes
+
+
+
+
+
+
+
+
+
+
+ Build functions
+
+
+
+
@@ -840,6 +862,34 @@
+
+
+
+
+
+
+
+
+ Evidence elimination (discussion)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
@@ -918,6 +968,62 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Merge command
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Kmer filtering (rebuild/dump/unitig)
+
+
+
+
+
+
+
+
+
+
+
+
+
@@ -1165,6 +1271,28 @@
+
+
+
+
+
+
+ Evidence modes
+
+
+
+
+
+
+
+
+
+
+ Build functions
+
+
+
+
@@ -1226,26 +1354,26 @@
Why two phases are needed
Kmer indexing per partition proceeds in two phases. The separation is necessary because the exact number of surviving unique kmers is not known until after counting and filtering low-abundance kmers.
Phase 1 — provisional MPHF + kmer spectrum
-Implemented in obikpartitionner::KmerPartition::count_kmer().
+Implemented in obikpartitionner::KmerPartition::count_kmer() → count_partition().
-Pass 1 : read the dereplicated superkmer file; enumerate all unique canonical kmers into a HashSet. Exact count known after this pass.Build a provisional MPHF (GOFunction from the ph crate) over the exact kmer set. Produces mphf1.bin.Create counts1.bin : one zero-initialised u32 per MPHF slot (mmap'd).Pass 2 : re-read the dereplicated file; for each kmer, query mphf1.get(kmer) and atomically accumulate the superkmer count into counts1[slot].External sort : read the dereplicated superkmer file; extract the raw u64 canonical kmer value for every kmer of every superkmer. Sort in RAM-bounded chunks (adaptive budget: 40% of available RAM ÷ n_threads, minimum 1 M kmers per chunk), then k-way merge with inline dedup. Result: sorted_unique.bin — a flat array of f0 distinct sorted u64 values. Exact kmer count f0 is known at this point.Build provisional MPHF (ptr_hash, same configuration as phase 2) over sorted_unique.bin using new_from_par_iter. Delete sorted_unique.bin immediately after. Persist to mphf1.bin.Create counts1.bin : PersistentCompactIntVec with f0 slots, zero-initialised.Accumulation pass : re-read the dereplicated superkmer file; for each kmer in each superkmer, compute slot = mphf.index(kmer.raw()) and increment counts1[slot] by the superkmer's COUNT.Build kmer frequency spectrum from counts1: histogram {count → n_kmers}, totals f0 (distinct kmers) and f1 (total abundance). Written to kmer_spectrum_raw.json per partition, then merged globally.
Files produced per partition:
part_XXXXX/
- mphf1.bin — GOFunction (provisional MPHF, discarded after phase 2)
- counts1.bin — [u32; n_kmers] kmer counts, mmap'd
+ mphf1.bin — ptr_hash provisional MPHF (discarded after phase 2)
+ counts1.bin — PersistentCompactIntVec, f0 × u32 kmer counts
kmer_spectrum_raw.json — local frequency spectrum
Phase 2 — definitive MPHF
After filtering (applying a min-count threshold derived from the spectrum) and building the local De Bruijn graph + unitigs (see Construction pipeline ), the exact filtered kmer set is available via unitigs.bin.
-MphfLayer::build is called on the unitig file:
+MphfLayer::build(dir, block_bits, mode: &IndexMode, fill_slot) is called on the unitig directory:
-Pass 1 : iterate all canonical kmers from unitigs.bin in parallel, build and store mphf.bin (ptr_hash).Pass 2 : iterate sequentially, fill evidence.bin, call the mode-specific fill_slot callback.Pass 1 (parallel): a CanonicalKmerIter — clonable via Arc<Mmap>, no file reopening — is passed directly to new_from_par_iter via par_bridge(). No .idx is read or created at this stage; parallelism is at partition/layer level, not within a single MPHF. Produces mphf.bin.Pass 2 (sequential): iterate with iter_indexed_canonical_kmers; fill evidence files; call fill_slot(slot, kmer) callback per kmer. For Exact/Hybrid, .idx is written at the end of this pass — never earlier.
mphf1.bin and counts1.bin are no longer needed after phase 2 and can be deleted.
FMPH/FMPHGO (ph crate, Beling, ACM JEA 2023):
~2.1 bits/key — most compact; good query speed; deterministic construction
-Works well from an exact or slightly overestimated count
-GOFunction (group-oriented variant) is the specific type usedGOFunction (group-oriented variant) was the original phase-1 choice; eliminated when the external sort made the exact count available at phase 1 as well
MPHF choice per phase
-Phase 1 (provisional, discarded after spectrum computation): ph::fmph::GOFunction. Compact, fast to build from the exact post-dedup kmer set. Query speed is secondary — the structure is only used during pass 2 of count_kmer.
-Phase 2 (persistent, queried repeatedly): ptr_hash . Exact key count is available from the unitig index; ptr_hash query speed (≥2.1×) and construction speed (≥3.1× over FMPH) are the decisive factors. The 2.4 bits/key overhead is acceptable.
-boomphf is eliminated: largest space overhead, streaming advantage does not apply.
+Both phases : ptr_hash , same type alias and construction parameters. The external sort (phase 1) and the unitig index (phase 2) both provide the exact key count before MPHF construction, so ptr_hash's requirement is satisfied in both cases. Using a single MPHF implementation removes the ph crate dependency.
+boomphf: eliminated — largest space overhead, streaming advantage no longer needed. FMPH/GOFunction: eliminated — exact count available, ptr_hash is faster at equivalent compactness.
Space at scale
For 1 024 partitions × 100 M kmers/partition (phase 2 index, after filtering):
@@ -1320,9 +1446,12 @@
Layer structure
Each layer is a self-contained unit. See obilayeredmap for the full on-disk layout. The MPHF-relevant files are:
layer_i/
- unitigs.bin — packed 2-bit nucleotide sequences (kmer evidence)
+ unitigs.bin — packed 2-bit nucleotide sequences (kmer evidence source)
+ unitigs.bin.idx — random-access block index (block_bits controls granularity)
mphf.bin — ptr_hash phase-2 MPHF
- evidence.bin — n × u32: (chunk_id: 25 bits | rank: 7 bits) per slot
+ evidence.bin — n × (chunk_id: 25 bits | rank: 7 bits) per slot [exact mode]
+ fingerprint.bin — n × b-bit fingerprints per slot [approx mode]
+ [no layer_meta.json — mode stored once in partition-level meta.json]
Layers are disjoint : a canonical kmer belongs to exactly one layer. Layer 0 is built from dataset A. Adding dataset B:
@@ -1330,17 +1459,43 @@
Collect kmers of B not present in any layer → set B \ A.
Build layer 1 from B \ A (dereplicate → count → De Bruijn → unitigs → MphfLayer::build).
+Evidence modes
+Three evidence modes are supported via IndexMode, stored once in PartitionMeta at partition root. There is no layer_meta.json.
+Exact (IndexMode::Exact): evidence.bin stores one (chunk_id, rank) pair per MPHF slot. Verification reconstructs the kmer and compares to the query. Zero false positives. .idx required at query time.
+Approx (IndexMode::Approx { b, z }): fingerprint.bin stores a b-bit hash per slot. False-positive rate 1/2^b per query; Findere z-parameter reduces window FP to ≈ 1/2^(b·z). No .idx written or needed.
+Hybrid (IndexMode::Hybrid { b, z }): both fingerprint.bin and evidence.bin + .idx. find() uses the fingerprint (O(1)); find_strict() uses exact evidence (O(1)).
+Build functions
+MphfLayer::build(dir, block_bits, mode: &IndexMode, fill_slot)
+ Pass 1: CanonicalKmerIter + par_bridge() → build mphf.bin (no .idx used)
+ Pass 2: sequential iter → fill evidence files + call fill_slot
+ .idx written last for Exact/Hybrid (query-time only)
+
+MphfLayer::build_exact_evidence(dir, block_bits)
+ Post-hoc: builds evidence.bin + .idx from existing mphf.bin + unitigs.bin
+ Uses open_sequential(); no .idx required on entry
+
+MphfLayer::build_approx_evidence(dir, b, z)
+ Post-hoc: builds fingerprint.bin from existing mphf.bin + unitigs.bin
+ Uses open_sequential(); never writes .idx
+There is no build_evidence dispatch wrapper. Callers choose the appropriate post-hoc build directly.
+In obikpartitionner, build_index_layer receives block_bits: u8 from IndexConfig::block_bits and forwards it directly to Layer::build and Layer::build_approx_evidence.
Membership verification
-ptr_hash maps any input to a valid slot — it does not natively detect absent keys. Membership is verified using the evidence entry: decode the kmer from (chunk_id, rank) and compare to the query. A mismatch means the kmer is absent from this layer; probe the next layer.
+ptr_hash maps any input to a valid slot — it does not natively detect absent keys. Membership is verified using the evidence entry:
+
+Exact : decode (chunk_id, rank) from evidence.bin; reconstruct the kmer via unitigs.verify_canonical_kmer; compare to query.Approx : compare kmer.seq_hash() to the b-bit fingerprint stored at the slot.
+A mismatch in either mode means the kmer is absent from this layer; probe the next layer.
Query algorithm
fn query(kmer) → Option<(layer_index, slot)>:
for (i, layer) in layers.iter().enumerate():
slot = layer.mphf.index(kmer)
- if layer.evidence.decode(slot) matches kmer:
+ if layer.evidence.matches(slot, kmer): // exact or approx dispatch
return Some((i, slot))
return None
Expected probe depth: 1 for kmers in layer 0. Each probe is a ptr_hash lookup (~10 ns) plus one evidence decode.
+MphfLayer::find dispatches on LayerEvidence at O(1) — no panicking find_exact/find_approx methods. find_strict always performs an exact check: O(1) for Exact/Hybrid, O(n) sequential scan for Approx. Expected probe depth: 1 for kmers in layer 0. Each probe is a ptr_hash lookup (~10 ns) plus one evidence check.
Merging layers
Two layer chains can be merged by re-indexing their union through the full pipeline. This is expensive (full rebuild) but produces an optimal single-layer index. Merge is a maintenance operation, not a query-path requirement.
diff --git a/doc/implementation/obilayeredmap.refs/index.html b/doc/implementation/obilayeredmap.refs/index.html
new file mode 100644
index 0000000..2a11dce
--- /dev/null
+++ b/doc/implementation/obilayeredmap.refs/index.html
@@ -0,0 +1,1068 @@
+
+
+
+
+
+ Obilayeredmap.refs - obikmer
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Coverage: implementation/obilayeredmap.md
+Code couvert
+
+obilayeredmap/src/mphf_layer.rs — MphfLayer, LayerEvidence enum (Exact/Approx), find(), find_exact(), find_approx()obilayeredmap/src/layer.rs — Layer, trait LayerData, modes () / PersistentCompactIntMatrix / PersistentBitMatrix, build(), build_evidence(), append_genome_column() obilayeredmap/src/map.rs — LayeredMap, push_layer(), query() obilayeredmap/src/evidence.rs — Evidence, EvidenceWriter, encodage chunk_id:rankobilayeredmap/src/fingerprint.rs — FingerprintVec, FingerprintVecWriter, matches()obilayeredmap/src/meta.rs — LayerMeta, EvidenceKind (Exact / Approx { b, z })
+Notes
+FORT RISQUE DE DÉRIVE. C'est le fichier le plus affecté par les changements récents :
+- EvidenceKind (Exact / Approx) est désormais un concept de premier plan — toute la sémantique de query en dépend
+- LayerEvidence enum interne à MphfLayer : dispatch transparent find() → find_exact() ou find_approx()
+- fingerprint.rs : module entièrement nouveau (FingerprintVec + FingerprintVecWriter)
+- build_evidence() / build_exact_evidence() / build_approx_evidence() sont nouveaux
+- block_bits dans les fonctions build : O(1) garanti avec le chemin chaud explicit pour block_bits=0
+- Séparation open() (accès aléatoire, requiert .idx) vs open_sequential() (itération seule)
+Pratiquement toute cette page est à réécrire.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Evidence elimination (discussion)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
@@ -729,6 +757,17 @@
+
+
+
+
+
+
+ Index mode (homogeneity invariant)
+
+
+
+
@@ -740,6 +779,34 @@
+
+
+
+
@@ -751,6 +818,73 @@
+
+
+
+
+
+
+
+
+
+
+ FingerprintVec and FingerprintVecWriter
+
+
+
+
+
+
+
+
+
+
+ LayeredMap\<D> — collection of layers
+
+
+
+
+
+
+
+
@@ -776,10 +910,10 @@
-
+
- Evidence encoding
+ Evidence encoding (exact)
@@ -798,14 +932,53 @@
-
+
- Query path
+ Column append and merge support
+
+
+
+
@@ -895,6 +1068,62 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Merge command
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Kmer filtering (rebuild/dump/unitig)
+
+
+
+
+
+
+
+
+
+
+
+
+
@@ -1058,6 +1287,17 @@
+
+
+
+
+
+
+ Index mode (homogeneity invariant)
+
+
+
+
@@ -1069,6 +1309,34 @@
+
+
+
+
@@ -1080,6 +1348,73 @@
+
+
+
+
+
+
+
+
+
+
+ FingerprintVec and FingerprintVecWriter
+
+
+
+
+
+
+
+
+
+
+ LayeredMap\<D> — collection of layers
+
+
+
+
+
+
+
+
@@ -1105,10 +1440,10 @@
-
+
- Evidence encoding
+ Evidence encoding (exact)
@@ -1127,14 +1462,53 @@
-
+
- Query path
+ Column append and merge support
+
+
+
+
@@ -1178,7 +1552,7 @@
obilayeredmap — layered kmer index crate
Purpose
-obilayeredmap implements a persistent, incrementally extensible kmer index. The index is organised in three levels: index root → partition → layer . Each layer covers a disjoint kmer set and wraps a ptr_hash MPHF with associated per-slot data. Adding a new dataset never rebuilds existing layers.
+obilayeredmap implements a persistent, incrementally extensible kmer index. Each layer covers a disjoint kmer set and wraps a ptr_hash MPHF with associated per-slot data. Adding a new dataset never rebuilds existing layers.
Three usage modes
The MPHF + evidence infrastructure is the same for all modes. The payload varies.
@@ -1214,34 +1588,65 @@
Both PersistentCompactIntMatrix and PersistentBitMatrix come from the obicompactvec crate.
Index mode (homogeneity invariant)
+A partitioned index is homogeneous: every layer within a partition shares the same mode. The mode is determined once at LayeredMap::open() from PartitionMeta.mode and passed to each Layer::open() — no per-layer file is read.
+#[derive(Serialize, Deserialize, Default)]
+#[serde(tag = "type" , rename_all = "snake_case" )]
+pub enum IndexMode {
+ #[default]
+ Exact ,
+ Approx { b : u8 , z : u8 },
+ Hybrid { b : u8 , z : u8 },
+}
+IndexMode is stored once in PartitionMeta (meta.json at partition root). There is no layer_meta.json.
+
+Exact : writes evidence.bin + unitigs.bin.idx. Zero false positives.Approx : writes fingerprint.bin only. FP rate per kmer = 1/2^b; with Findere z-parameter, z consecutive kmers must all match → effective window FP ≈ 1/2^(b·z). No .idx written or required.Hybrid : writes both fingerprint.bin and evidence.bin + .idx. find() uses the fingerprint (fast, O(1)); find_strict() uses exact evidence.
+MphfLayer — autonomous kmer → slot mapping
-MphfLayer encapsulates the MPHF + evidence + unitig spine for one layer. It is independent of any payload data.
+MphfLayer encapsulates the MPHF and evidence store for one layer. It is independent of any payload.
pub struct MphfLayer {
- mphf : Mphf ,
- evidence : Evidence ,
- unitigs : UnitigFileReader ,
- n : usize , // number of indexed kmers = number of MPHF slots
+ mphf : Mphf ,
+ ev : LayerEvidence , // loaded at open() time
+ n : usize ,
}
Public API:
-impl MphfLayer {
- pub fn open ( dir : & Path ) -> OLMResult < Self >
- pub fn find ( & self , kmer : CanonicalKmer ) -> Option < usize > // Some(slot) or None
- pub fn n ( & self ) -> usize
- pub fn unitig_writer ( dir : & Path ) -> OLMResult < UnitigFileWriter >
- pub ( crate ) fn build (
- dir : & Path ,
- fill_slot : & mut impl FnMut ( usize , CanonicalKmer ) -> OLMResult < () > ,
- ) -> OLMResult < usize >
+LayerEvidence is an internal enum, not public:
+enum LayerEvidence {
+ Exact { evidence : Evidence , unitigs : UnitigFileReader },
+ Approx { fingerprint : FingerprintVec , unitigs_path : PathBuf },
+ Hybrid { evidence : Evidence , unitigs : UnitigFileReader , fingerprint : FingerprintVec },
}
find returns Some(slot) only after verifying via evidence that the kmer is actually indexed. It returns None for absent keys (ptr_hash maps any input to a valid slot; evidence verification is the only correct-membership test).
-build runs two sequential passes over unitigs.bin:
+MphfLayer::open(dir, mode: &IndexMode) receives the mode from PartitionMeta — no per-layer file is read.
+Query API
+Two public query methods, both returning Option<usize> (slot index):
+pub fn find ( & self , kmer : CanonicalKmer ) -> Option < usize >
+pub fn find_strict ( & self , kmer : CanonicalKmer ) -> Option < usize >
+
+find: O(1) auto-dispatch. Exact/Hybrid → exact evidence check. Approx/Hybrid → fingerprint comparison.find_strict: always exact. Exact/Hybrid → O(1) evidence check. Approx → O(n) sequential scan (no .idx).
+There are no find_exact/find_approx methods; panicking dispatch is eliminated.
+Build surface
+// Full MPHF + evidence build (two-pass)
+pub ( crate ) fn build ( dir , block_bits , mode : & IndexMode , fill_slot ) -> OLMResult < usize >
+
+// Evidence-only post-hoc builds (MPHF already present)
+pub fn build_exact_evidence ( dir , block_bits ) -> OLMResult < usize >
+pub fn build_approx_evidence ( dir , b , z ) -> OLMResult < usize >
+MphfLayer::build runs two passes over unitigs.bin:
-Pass 1 : iterate all canonical kmers in parallel via rayon, construct and store mphf.bin. new_from_par_iter avoids materialising a full key Vec.Pass 2 : iterate again sequentially, fill evidence.bin, call fill_slot(slot, kmer) once per kmer for payload population. A compact n/8-byte seen-bitset verifies MPHF injectivity inline.Pass 1 (parallel via rayon): a CanonicalKmerIter (clonable, Arc<Mmap>, no file reopening) is passed to new_from_par_iter via par_bridge(). Produces mphf.bin. No .idx is read or created at this stage.Pass 2 (sequential): fill evidence files; call fill_slot(slot, kmer) per kmer. .idx is written last for Exact/Hybrid modes (query-time only).
-For empty layers (n = 0), build returns Ok(0) immediately after creating empty mphf.bin and evidence.bin.
+There is no build_evidence dispatch wrapper — callers invoke build_exact_evidence or build_approx_evidence directly.
+For empty layers (n = 0), all build variants return Ok(0) immediately after creating empty output files.
Layer\<D: LayerData> — MPHF + payload
Layer<D> pairs an MphfLayer with one payload store.
@@ -1261,7 +1666,7 @@
pub data : T ,
}
LayerData covers the read path only (open + read). Build signatures differ between modes and are not in the trait.
+LayerData covers the read path only (open + read). Build signatures differ between modes and are not part of the trait.
-Build signatures:
+Build signatures
// mode 1
impl Layer < () > {
- pub fn build ( out_dir : & Path ) -> OLMResult < usize >
+ pub fn build ( out_dir : & Path , block_bits : u8 , mode : & IndexMode ) -> OLMResult < usize >
}
// mode 2
impl Layer < PersistentCompactIntMatrix > {
- pub fn build ( out_dir : & Path , count_of : impl Fn ( CanonicalKmer ) -> u32 ) -> OLMResult < usize >
- pub fn build_from_map ( out_dir : & Path , counts : & HashMap < CanonicalKmer , u32 > ) -> OLMResult < usize >
+ pub fn build ( out_dir : & Path , block_bits : u8 , mode : & IndexMode ,
+ count_of : impl Fn ( CanonicalKmer ) -> u32 ) -> OLMResult < usize >
+ pub fn build_from_map ( out_dir : & Path , block_bits : u8 , mode : & IndexMode ,
+ counts : & HashMap < CanonicalKmer , u32 > ) -> OLMResult < usize >
}
// mode 3
impl Layer < PersistentBitMatrix > {
- pub fn build_presence (
- out_dir : & Path ,
- n_genomes : usize ,
- present_in : impl Fn ( CanonicalKmer , usize ) -> bool ,
- ) -> OLMResult < usize >
+ pub fn build_presence ( out_dir : & Path , block_bits : u8 , mode : & IndexMode ,
+ n_genomes : usize ,
+ present_in : impl Fn ( CanonicalKmer , usize ) -> bool ) -> OLMResult < usize >
}
All build impls delegate MPHF + evidence construction to MphfLayer::build via a mode-specific fill_slot callback. Mode 2 pre-reads n_kmers from unitigs.bin to size the PersistentCompactIntMatrixBuilder before calling MphfLayer::build. Mode 3 does the same for PersistentBitMatrixBuilder.
+All build impls delegate to MphfLayer::build via a mode-specific fill_slot callback. The mode parameter is forwarded directly — no LayerMeta is written.
+Evidence-only post-hoc builds are accessible directly on Layer<D>:
+impl < D : LayerData > Layer < D > {
+ pub fn build_exact_evidence ( layer_dir : & Path , block_bits : u8 ) -> OLMResult < usize >
+ pub fn build_approx_evidence ( layer_dir : & Path , b : u8 , z : u8 ) -> OLMResult < usize >
+}
+There is no build_evidence dispatch wrapper.
+FingerprintVec and FingerprintVecWriter
+Approximate evidence is stored as a packed b-bit array, one fingerprint per MPHF slot.
+fingerprint.bin format:
+ magic: b"FPVF" (4 bytes)
+ b: u8 (bits per fingerprint, 1..=64)
+ padding: [0u8; 3]
+ n: u64 LE (number of slots)
+ data: packed bits, ceil(n*b/8) bytes, Lsb0 order
+impl FingerprintVec {
+ pub fn open ( path : & Path ) -> OLMResult < Self >
+ pub fn get ( & self , slot : usize ) -> u64
+ pub fn matches ( & self , slot : usize , fingerprint : u64 ) -> bool
+ pub fn n ( & self ) -> usize
+ pub fn b ( & self ) -> u8
+}
+matches(slot, hash) extracts the b-bit fingerprint stored at slot and compares it to the low b bits of hash. It is the core operation of find_approx.
+LayeredMap\<D> — collection of layers
+LayeredMap<D> wraps Vec<Layer<D>> for a single partition directory.
+pub struct LayeredMap < D : LayerData = () > {
+ root : PathBuf ,
+ meta : PartitionMeta ,
+ layers : Vec < Layer < D >> ,
+}
+PartitionMeta (meta.json at the partition root) stores n_layers.
+Common methods
+pub fn open ( root : & Path ) -> OLMResult < Self >
+pub fn create ( root : & Path , mode : IndexMode ) -> OLMResult < Self >
+pub fn n_layers ( & self ) -> usize
+pub fn layer ( & self , i : usize ) -> & Layer < D >
+pub fn mode ( & self ) -> & IndexMode
+pub fn query ( & self , kmer : CanonicalKmer ) -> Option < ( usize , Hit < D :: Item > ) >
+pub fn next_layer_writer ( & self ) -> OLMResult < UnitigFileWriter >
+open reads PartitionMeta once, extracts mode, and passes it to every Layer::open — no per-layer file is read. create stores the given mode in PartitionMeta.
+query probes layers in order and returns (layer_index, Hit) on the first match. Expected probe depth: 1 for kmers in layer 0.
+push_layer
+push_layer builds the next layer from a unitigs.bin already written via next_layer_writer, using DEFAULT_BLOCK_BITS:
+// mode 1
+impl LayeredMap < () > {
+ pub fn push_layer ( & mut self ) -> OLMResult < usize >
+}
+
+// mode 2
+impl LayeredMap < PersistentCompactIntMatrix > {
+ pub fn push_layer ( & mut self , count_of : impl Fn ( CanonicalKmer ) -> u32 ) -> OLMResult < usize >
+ pub fn push_layer_from_map ( & mut self , counts : & HashMap < CanonicalKmer , u32 > ) -> OLMResult < usize >
+}
+Mode 3 (PersistentBitMatrix) has no push_layer on LayeredMap; callers build directly via Layer<PersistentBitMatrix>::build_presence.
LayeredStore\<S> and aggregation traits
LayeredStore<S> is a generic aggregation wrapper over Vec<S>. It propagates three traits from obicompactvec::traits up the hierarchy via blanket impls:
@@ -1320,11 +1786,6 @@
impl < S : BitPartials > BitPartials for LayeredStore < S > { … } // element-wise Σ partials
Because blanket impls compose, LayeredStore<LayeredStore<S>> automatically inherits all three traits when S does — providing the partitioned level without a separate type.
-Aggregation hierarchy:
-PersistentCompactIntMatrix implements CountPartials
-LayeredStore<PersistentCompactIntMatrix> via blanket impl (one partition)
-LayeredStore<LayeredStore<…>> via blanket impl (partitioned index)
-Leaf implementors (in obicompactvec):
@@ -1344,69 +1805,77 @@ LayeredStore<LayeredStore<…>> via blanket impl
-PersistentCompactIntVec and PersistentBitVec do not implement these traits — they are single-column primitives, not matrix-level aggregators.
See Kmer index architecture for the full trait API and the two-pass normalised-metric pattern.
On-disk structure
-