diff --git a/src/obikindex/src/siblings/alignment.rs b/src/obikindex/src/siblings/alignment.rs index f321d671..a0409617 100644 --- a/src/obikindex/src/siblings/alignment.rs +++ b/src/obikindex/src/siblings/alignment.rs @@ -1,5 +1,3 @@ -use rayon::prelude::*; - use obikpartitionner::KmerPartition; use obisys::progress_bar; @@ -74,24 +72,26 @@ impl KmerIndex { let layer_dirs = self.sibling_layer_dirs()?; let pb = progress_bar("snp_pseudo_alignment", layer_dirs.len() as u64, "layers"); - // `Vec>` per layer, one entry (column) per variable family; - // `par_iter().map(...).collect()` on this indexed source preserves - // input order, so concatenating the results below in order gives a - // single deterministic column order across the whole index. - let partials: Vec>> = layer_dirs - .par_iter() - .map(|layer_dir| -> OKIResult>> { - let families = scan_layer_families(layer_dir, n_parts, n_genomes, with_counts, k, &cache)?; - let columns = families - .into_iter() - .filter(|f| f.mask.family_size() >= 2) // monomorphic family — no signal, skip - .map(|f| f.genome_mask.iter().map(|&m| iupac_code(m)).collect()) - .collect(); - - pb.inc(1); - Ok(columns) - }) - .collect::>>()?; + // One layer at a time, not `par_iter()` over layers — same + // rationale as `build_sibling_annex`: running many layers' + // `scan_layer_families` concurrently would each group their own + // lookups by partition internally, but interleave those sweeps + // across layers at the OS level, scattering page-cache access over + // every partition at once again and defeating the whole point of + // the grouping. `Vec>` per layer, one entry (column) per + // variable family, appended in layer order for a single + // deterministic column order across the whole index. + let mut partials: Vec>> = Vec::with_capacity(layer_dirs.len()); + for layer_dir in &layer_dirs { + let families = scan_layer_families(layer_dir, n_parts, n_genomes, with_counts, k, &cache)?; + let columns = families + .into_iter() + .filter(|f| f.mask.family_size() >= 2) // monomorphic family — no signal, skip + .map(|f| f.genome_mask.iter().map(|&m| iupac_code(m)).collect()) + .collect(); + partials.push(columns); + pb.inc(1); + } pb.finish_and_clear(); let mut sequences: Vec> = vec![Vec::new(); n_genomes]; diff --git a/src/obikindex/src/siblings/cardinality.rs b/src/obikindex/src/siblings/cardinality.rs index 84bd5b38..652cc09b 100644 --- a/src/obikindex/src/siblings/cardinality.rs +++ b/src/obikindex/src/siblings/cardinality.rs @@ -1,5 +1,4 @@ use ndarray::Array2; -use rayon::prelude::*; use obikpartitionner::KmerPartition; use obisys::progress_bar; @@ -77,45 +76,46 @@ impl KmerIndex { let layer_dirs = self.sibling_layer_dirs()?; let pb = progress_bar("cardinality_tally", layer_dirs.len() as u64, "layers"); - let partials: Vec<[[u64; 5]; 5]> = layer_dirs - .par_iter() - .map(|layer_dir| -> OKIResult<[[u64; 5]; 5]> { - let families = scan_layer_families(layer_dir, n_parts, n_genomes, with_counts, k, &cache)?; - let mut counts = [[0u64; 5]; 5]; + // One layer at a time — see `snp_pseudo_alignment`'s comment for why + // `par_iter()` over layers would defeat `scan_layer_families`'s + // partition-grouped locality. + let mut partials: Vec<[[u64; 5]; 5]> = Vec::with_capacity(layer_dirs.len()); + 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 { - if family.mask.family_size() < 2 { - // Fully invariant family (never varies anywhere in - // the index) — genome-wide background, not - // SNP-adjacent signal; would otherwise swamp the - // diagonal (`c=1/c=1` etc.), which needs to reflect - // the same variable-families-only population the - // `+ASC`-corrected alignment/likelihood actually - // models. See `base_pair_tally`'s `variable` gate - // on its own `same` diagonal for the matching fix. - continue; - } + for family in &families { + if family.mask.family_size() < 2 { + // Fully invariant family (never varies anywhere in + // the index) — genome-wide background, not + // SNP-adjacent signal; would otherwise swamp the + // diagonal (`c=1/c=1` etc.), which needs to reflect + // the same variable-families-only population the + // `+ASC`-corrected alignment/likelihood actually + // models. See `base_pair_tally`'s `variable` gate + // on its own `same` diagonal for the matching fix. + continue; + } - let genome_mask = &family.genome_mask; - for i in 0..n_genomes { - let card_i = genome_mask[i].count_ones() as usize; - for j in (i + 1)..n_genomes { - if !included[[i, j]] { - continue; - } - let card_j = genome_mask[j].count_ones() as usize; - counts[card_i][card_j] += 1; - if card_i != card_j { - counts[card_j][card_i] += 1; - } + let genome_mask = &family.genome_mask; + for i in 0..n_genomes { + let card_i = genome_mask[i].count_ones() as usize; + for j in (i + 1)..n_genomes { + if !included[[i, j]] { + continue; + } + let card_j = genome_mask[j].count_ones() as usize; + counts[card_i][card_j] += 1; + if card_i != card_j { + counts[card_j][card_i] += 1; } } } + } - pb.inc(1); - Ok(counts) - }) - .collect::>>()?; + partials.push(counts); + pb.inc(1); + } pb.finish_and_clear(); let mut total = [[0u64; 5]; 5]; diff --git a/src/obikindex/src/siblings/distance.rs b/src/obikindex/src/siblings/distance.rs index 6b5d37e0..746bb062 100644 --- a/src/obikindex/src/siblings/distance.rs +++ b/src/obikindex/src/siblings/distance.rs @@ -1,5 +1,4 @@ use ndarray::Array2; -use rayon::prelude::*; use obikpartitionner::KmerPartition; use obisys::progress_bar; @@ -50,8 +49,11 @@ impl KmerIndex { /// should only reflect genuine SNP-adjacent agreement, not the /// genome-wide invariant background, need it (see /// [`base_pair_tally`](Self::base_pair_tally)'s `same` field). Layers - /// are processed in parallel (rayon); each gets its own accumulator - /// from `zero()`, combined pairwise via `combine`. + /// are processed one at a time, not in parallel — see + /// `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( &self, label: &str, @@ -81,44 +83,41 @@ impl KmerIndex { let layer_dirs = self.sibling_layer_dirs()?; let pb = progress_bar(label, layer_dirs.len() as u64, "layers"); - let partials: Vec = layer_dirs - .par_iter() - .map(|layer_dir| -> OKIResult { - let mut acc = zero(); - let families = scan_layer_families(layer_dir, n_parts, n_genomes, with_counts, k, &cache)?; + // One layer at a time — see `snp_pseudo_alignment`'s comment for why + // `par_iter()` over layers would defeat `scan_layer_families`'s + // partition-grouped locality. + let mut total = zero(); + 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 { - let variable = family.mask.family_size() >= 2; - let genome_mask = &family.genome_mask; + for family in &families { + let variable = family.mask.family_size() >= 2; + let genome_mask = &family.genome_mask; - // Per genome: which single form (if exactly one) it - // carries — `None` (a `popcount != 1` mask) once a - // second form is seen, ambiguous/not single-copy, - // ineligible for either side of a pair. - let single_form = |g: usize| -> Option { - let m = genome_mask[g]; - (m.count_ones() == 1).then(|| m.trailing_zeros() as u8) - }; + // Per genome: which single form (if exactly one) it + // carries — `None` (a `popcount != 1` mask) once a + // second form is seen, ambiguous/not single-copy, + // ineligible for either side of a pair. + let single_form = |g: usize| -> Option { + let m = genome_mask[g]; + (m.count_ones() == 1).then(|| m.trailing_zeros() as u8) + }; - for i in 0..n_genomes { - let Some(bi) = single_form(i) else { continue }; - for j in (i + 1)..n_genomes { - let Some(bj) = single_form(j) else { continue }; - on_pair(&mut acc, i, j, bi, bj, variable); - } + for i in 0..n_genomes { + let Some(bi) = single_form(i) else { continue }; + for j in (i + 1)..n_genomes { + let Some(bj) = single_form(j) else { continue }; + on_pair(&mut acc, i, j, bi, bj, variable); } } + } - pb.inc(1); - Ok(acc) - }) - .collect::>>()?; + total = combine(total, acc); + pb.inc(1); + } pb.finish_and_clear(); - let mut total = zero(); - for partial in partials { - total = combine(total, partial); - } Ok(total) } diff --git a/src/obikindex/src/siblings/stats.rs b/src/obikindex/src/siblings/stats.rs index ac659773..9d8c7b71 100644 --- a/src/obikindex/src/siblings/stats.rs +++ b/src/obikindex/src/siblings/stats.rs @@ -1,5 +1,3 @@ -use rayon::prelude::*; - use obikpartitionner::KmerPartition; use obisys::progress_bar; @@ -58,52 +56,32 @@ impl KmerIndex { let cache = PartitionCache::build(&partition, n_parts, with_counts)?; let layer_dirs = self.sibling_layer_dirs()?; - // One layer's worth of work, parallelised across layers with Rayon - // — independent, read-only, each producing its own partial tally - // merged at the end. + // One layer at a time, not parallelised across layers — see + // `snp_pseudo_alignment`'s comment for why `par_iter()` over layers + // would defeat `scan_layer_families`'s partition-grouped locality. let pb = progress_bar("sibling_annex_stats", layer_dirs.len() as u64, "layers"); - let partials: Vec = layer_dirs - .par_iter() - .map(|layer_dir| -> OKIResult { - 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::>>()?; - pb.finish_and_clear(); - let mut stats = SiblingAnnexStats { per_genome: vec![[0u64; 4]; n_genomes], ..Default::default() }; - for part in partials { - for s in 0..4 { - stats.counts[s] += part.counts[s]; - } - for g in 0..n_genomes { - for s in 0..4 { - stats.per_genome[g][s] += part.per_genome[g][s]; + for layer_dir in &layer_dirs { + 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); } + pb.finish_and_clear(); + Ok(stats) } }