OBIKmers/obikmer
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 27 Wiki Activity

Push zpwxxpnpktps #67

Merged
coissac merged 46 commits from push-zpwxxpnpktps into main 2026-08-17 09:41:42 +00:00
Conversation 0 Commits 46 Files Changed 405
+104 -127
4 changed files with 104 additions and 127 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Showing only changes of commit 0e2e3b5bae - Show all commits
src/obikindex/src/siblings/alignment.rs
+20 -20
View File
@@ -1,5 +1,3 @@
use rayon::prelude::*;
use obikpartitionner::KmerPartition; use obikpartitionner::KmerPartition;
use obisys::progress_bar; use obisys::progress_bar;
@@ -74,24 +72,26 @@ impl KmerIndex {
let layer_dirs = self.sibling_layer_dirs()?; let layer_dirs = self.sibling_layer_dirs()?;
let pb = progress_bar("snp_pseudo_alignment", layer_dirs.len() as u64, "layers"); let pb = progress_bar("snp_pseudo_alignment", layer_dirs.len() as u64, "layers");
// `Vec<Vec<u8>>` per layer, one entry (column) per variable family; // One layer at a time, not `par_iter()` over layers — same
// `par_iter().map(...).collect()` on this indexed source preserves // rationale as `build_sibling_annex`: running many layers'
// input order, so concatenating the results below in order gives a // `scan_layer_families` concurrently would each group their own
// single deterministic column order across the whole index. // lookups by partition internally, but interleave those sweeps
let partials: Vec<Vec<Vec<u8>>> = layer_dirs // across layers at the OS level, scattering page-cache access over
.par_iter() // every partition at once again and defeating the whole point of
.map(|layer_dir| -> OKIResult<Vec<Vec<u8>>> { // the grouping. `Vec<Vec<u8>>` per layer, one entry (column) per
let families = scan_layer_families(layer_dir, n_parts, n_genomes, with_counts, k, &cache)?; // variable family, appended in layer order for a single
let columns = families // deterministic column order across the whole index.
.into_iter() let mut partials: Vec<Vec<Vec<u8>>> = Vec::with_capacity(layer_dirs.len());
.filter(|f| f.mask.family_size() >= 2) // monomorphic family — no signal, skip for layer_dir in &layer_dirs {
.map(|f| f.genome_mask.iter().map(|&m| iupac_code(m)).collect()) let families = scan_layer_families(layer_dir, n_parts, n_genomes, with_counts, k, &cache)?;
.collect(); let columns = families
.into_iter()
pb.inc(1); .filter(|f| f.mask.family_size() >= 2) // monomorphic family — no signal, skip
Ok(columns) .map(|f| f.genome_mask.iter().map(|&m| iupac_code(m)).collect())
}) .collect();
.collect::<OKIResult<Vec<_>>>()?; partials.push(columns);
pb.inc(1);
}
pb.finish_and_clear(); pb.finish_and_clear();
let mut sequences: Vec<Vec<u8>> = vec![Vec::new(); n_genomes]; let mut sequences: Vec<Vec<u8>> = vec![Vec::new(); n_genomes];
src/obikindex/src/siblings/cardinality.rs
+34 -34
View File
@@ -1,5 +1,4 @@
use ndarray::Array2; use ndarray::Array2;
use rayon::prelude::*;
use obikpartitionner::KmerPartition; use obikpartitionner::KmerPartition;
use obisys::progress_bar; use obisys::progress_bar;
@@ -77,45 +76,46 @@ impl KmerIndex {
let layer_dirs = self.sibling_layer_dirs()?; let layer_dirs = self.sibling_layer_dirs()?;
let pb = progress_bar("cardinality_tally", layer_dirs.len() as u64, "layers"); let pb = progress_bar("cardinality_tally", layer_dirs.len() as u64, "layers");
let partials: Vec<[[u64; 5]; 5]> = layer_dirs // One layer at a time — see `snp_pseudo_alignment`'s comment for why
.par_iter() // `par_iter()` over layers would defeat `scan_layer_families`'s
.map(|layer_dir| -> OKIResult<[[u64; 5]; 5]> { // partition-grouped locality.
let families = scan_layer_families(layer_dir, n_parts, n_genomes, with_counts, k, &cache)?; let mut partials: Vec<[[u64; 5]; 5]> = Vec::with_capacity(layer_dirs.len());
let mut counts = [[0u64; 5]; 5]; for layer_dir in &layer_dirs {
let families = scan_layer_families(layer_dir, n_parts, n_genomes, with_counts, k, &cache)?;
let mut counts = [[0u64; 5]; 5];
for family in &families { for family in &families {
if family.mask.family_size() < 2 { if family.mask.family_size() < 2 {
// Fully invariant family (never varies anywhere in // Fully invariant family (never varies anywhere in
// the index) — genome-wide background, not // the index) — genome-wide background, not
// SNP-adjacent signal; would otherwise swamp the // SNP-adjacent signal; would otherwise swamp the
// diagonal (`c=1/c=1` etc.), which needs to reflect // diagonal (`c=1/c=1` etc.), which needs to reflect
// the same variable-families-only population the // the same variable-families-only population the
// `+ASC`-corrected alignment/likelihood actually // `+ASC`-corrected alignment/likelihood actually
// models. See `base_pair_tally`'s `variable` gate // models. See `base_pair_tally`'s `variable` gate
// on its own `same` diagonal for the matching fix. // on its own `same` diagonal for the matching fix.
continue; continue;
} }
let genome_mask = &family.genome_mask; let genome_mask = &family.genome_mask;
for i in 0..n_genomes { for i in 0..n_genomes {
let card_i = genome_mask[i].count_ones() as usize; let card_i = genome_mask[i].count_ones() as usize;
for j in (i + 1)..n_genomes { for j in (i + 1)..n_genomes {
if !included[[i, j]] { if !included[[i, j]] {
continue; continue;
} }
let card_j = genome_mask[j].count_ones() as usize; let card_j = genome_mask[j].count_ones() as usize;
counts[card_i][card_j] += 1; counts[card_i][card_j] += 1;
if card_i != card_j { if card_i != card_j {
counts[card_j][card_i] += 1; counts[card_j][card_i] += 1;
}
} }
} }
} }
}
pb.inc(1); partials.push(counts);
Ok(counts) pb.inc(1);
}) }
.collect::<OKIResult<Vec<_>>>()?;
pb.finish_and_clear(); pb.finish_and_clear();
let mut total = [[0u64; 5]; 5]; let mut total = [[0u64; 5]; 5];
src/obikindex/src/siblings/distance.rs
+32 -33
View File
@@ -1,5 +1,4 @@
use ndarray::Array2; use ndarray::Array2;
use rayon::prelude::*;
use obikpartitionner::KmerPartition; use obikpartitionner::KmerPartition;
use obisys::progress_bar; use obisys::progress_bar;
@@ -50,8 +49,11 @@ impl KmerIndex {
/// should only reflect genuine SNP-adjacent agreement, not the /// should only reflect genuine SNP-adjacent agreement, not the
/// genome-wide invariant background, need it (see /// genome-wide invariant background, need it (see
/// [`base_pair_tally`](Self::base_pair_tally)'s `same` field). Layers /// [`base_pair_tally`](Self::base_pair_tally)'s `same` field). Layers
/// are processed in parallel (rayon); each gets its own accumulator /// are processed one at a time, not in parallel — see
/// from `zero()`, combined pairwise via `combine`. /// `snp_pseudo_alignment`'s comment for why `par_iter()` over layers
/// would defeat `scan_layer_families`'s partition-grouped locality;
/// each layer gets its own accumulator from `zero()`, combined
/// pairwise via `combine`.
fn scan_family_pairs<Acc, F, C>( fn scan_family_pairs<Acc, F, C>(
&self, &self,
label: &str, label: &str,
@@ -81,44 +83,41 @@ impl KmerIndex {
let layer_dirs = self.sibling_layer_dirs()?; let layer_dirs = self.sibling_layer_dirs()?;
let pb = progress_bar(label, layer_dirs.len() as u64, "layers"); let pb = progress_bar(label, layer_dirs.len() as u64, "layers");
let partials: Vec<Acc> = layer_dirs // One layer at a time — see `snp_pseudo_alignment`'s comment for why
.par_iter() // `par_iter()` over layers would defeat `scan_layer_families`'s
.map(|layer_dir| -> OKIResult<Acc> { // partition-grouped locality.
let mut acc = zero(); let mut total = zero();
let families = scan_layer_families(layer_dir, n_parts, n_genomes, with_counts, k, &cache)?; for layer_dir in &layer_dirs {
let families = scan_layer_families(layer_dir, n_parts, n_genomes, with_counts, k, &cache)?;
let mut acc = zero();
for family in &families { for family in &families {
let variable = family.mask.family_size() >= 2; let variable = family.mask.family_size() >= 2;
let genome_mask = &family.genome_mask; let genome_mask = &family.genome_mask;
// Per genome: which single form (if exactly one) it // Per genome: which single form (if exactly one) it
// carries — `None` (a `popcount != 1` mask) once a // carries — `None` (a `popcount != 1` mask) once a
// second form is seen, ambiguous/not single-copy, // second form is seen, ambiguous/not single-copy,
// ineligible for either side of a pair. // ineligible for either side of a pair.
let single_form = |g: usize| -> Option<u8> { let single_form = |g: usize| -> Option<u8> {
let m = genome_mask[g]; let m = genome_mask[g];
(m.count_ones() == 1).then(|| m.trailing_zeros() as u8) (m.count_ones() == 1).then(|| m.trailing_zeros() as u8)
}; };
for i in 0..n_genomes { for i in 0..n_genomes {
let Some(bi) = single_form(i) else { continue }; let Some(bi) = single_form(i) else { continue };
for j in (i + 1)..n_genomes { for j in (i + 1)..n_genomes {
let Some(bj) = single_form(j) else { continue }; let Some(bj) = single_form(j) else { continue };
on_pair(&mut acc, i, j, bi, bj, variable); on_pair(&mut acc, i, j, bi, bj, variable);
}
} }
} }
}
pb.inc(1); total = combine(total, acc);
Ok(acc) pb.inc(1);
}) }
.collect::<OKIResult<Vec<_>>>()?;
pb.finish_and_clear(); pb.finish_and_clear();
let mut total = zero();
for partial in partials {
total = combine(total, partial);
}
Ok(total) Ok(total)
} }
src/obikindex/src/siblings/stats.rs
+18 -40
View File
@@ -1,5 +1,3 @@
use rayon::prelude::*;
use obikpartitionner::KmerPartition; use obikpartitionner::KmerPartition;
use obisys::progress_bar; use obisys::progress_bar;
@@ -58,52 +56,32 @@ impl KmerIndex {
let cache = PartitionCache::build(&partition, n_parts, with_counts)?; let cache = PartitionCache::build(&partition, n_parts, with_counts)?;
let layer_dirs = self.sibling_layer_dirs()?; let layer_dirs = self.sibling_layer_dirs()?;
// One layer's worth of work, parallelised across layers with Rayon // One layer at a time, not parallelised across layers — see
// — independent, read-only, each producing its own partial tally // `snp_pseudo_alignment`'s comment for why `par_iter()` over layers
// merged at the end. // would defeat `scan_layer_families`'s partition-grouped locality.
let pb = progress_bar("sibling_annex_stats", layer_dirs.len() as u64, "layers"); let pb = progress_bar("sibling_annex_stats", layer_dirs.len() as u64, "layers");
let partials: Vec<SiblingAnnexStats> = layer_dirs
.par_iter()
.map(|layer_dir| -> OKIResult<SiblingAnnexStats> {
let mut stats = SiblingAnnexStats {
per_genome: vec![[0u64; 4]; n_genomes],
..Default::default()
};
let families = scan_layer_families(layer_dir, n_parts, n_genomes, with_counts, k, &cache)?;
for family in &families {
// "Genome g represents this family" means g carries
// *any* of its members, not just the minorant's own —
// `genome_mask[g] != 0` is exactly that.
let s = family.mask.siblings() as usize;
stats.counts[s] += 1;
for (g, &m) in family.genome_mask.iter().enumerate() {
if m != 0 {
stats.per_genome[g][s] += 1;
}
}
}
pb.inc(1);
Ok(stats)
})
.collect::<OKIResult<Vec<_>>>()?;
pb.finish_and_clear();
let mut stats = SiblingAnnexStats { let mut stats = SiblingAnnexStats {
per_genome: vec![[0u64; 4]; n_genomes], per_genome: vec![[0u64; 4]; n_genomes],
..Default::default() ..Default::default()
}; };
for part in partials { for layer_dir in &layer_dirs {
for s in 0..4 { let families = scan_layer_families(layer_dir, n_parts, n_genomes, with_counts, k, &cache)?;
stats.counts[s] += part.counts[s]; for family in &families {
} // "Genome g represents this family" means g carries
for g in 0..n_genomes { // *any* of its members, not just the minorant's own —
for s in 0..4 { // `genome_mask[g] != 0` is exactly that.
stats.per_genome[g][s] += part.per_genome[g][s]; let s = family.mask.siblings() as usize;
stats.counts[s] += 1;
for (g, &m) in family.genome_mask.iter().enumerate() {
if m != 0 {
stats.per_genome[g][s] += 1;
}
} }
} }
pb.inc(1);
} }
pb.finish_and_clear();
Ok(stats) Ok(stats)
} }
} }