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

Consolidate compare_sparse as example and clean up project artifacts

Browse Source 
Restructure the project by moving the standalone compare_sparse utility into an example directory, removing Sankoff parameter configurations and benchmark scripts, updating version control ignores, and expanding the test suite with diagnostic checks and performance benchmarks.
This commit is contained in:
Eric Coissac
2026-08-20 12:48:42 +02:00
parent dc3d82f8db
commit dbd8af376c
13 changed files with 12 additions and 1184 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/obikphylo/src/siblings/tests.rs
+1 -632
View File
@@ -15,7 +15,7 @@ use super::cardinality::CardinalityExt;
use super::distance::DistanceExt;
use super::entropy::ShannonEntropyExt;
use super::entropy_annex::{EntropyAnnex, ENTROPY_ANNEX_FILE_NAME};
use super::helpers::{central_base, is_minorant};
use super::helpers::is_minorant;
use super::sankoff_bundle::SankoffBundleExt;
use super::stats::SiblingStatsExt;
use super::subsample::EntropyBias;
@@ -550,634 +550,3 @@ fn entropy_annex_builds_on_demand_and_biases_selection() {
assert!(find_layer(&selections_far).is_none(), "mu far from the true entropy must never select it, in any layer");
}
#[test]
#[ignore]
fn diag_real_index_layer_distribution() {
let idx = KmerIndex::open("/Users/coissac/Sync/travail/__MOI__/obikmer/benchmark/global_index_presence")
.expect("open real index");
let layer_dirs = super::family_scan::sibling_layer_dirs(&idx).expect("layer dirs");
let counts = super::subsample::non_monomorphic_counts(&layer_dirs).expect("counts");
let total: u64 = counts.iter().sum();
let n_zero = counts.iter().filter(|&&c| c == 0).count();
let max = counts.iter().max().unwrap();
let min_nonzero = counts.iter().filter(|&&c| c > 0).min().unwrap_or(&0);
println!("n_layers={} total={} n_zero_layers={} max={} min_nonzero={}",
counts.len(), total, n_zero, max, min_nonzero);
let n = 100_000usize;
let selections = super::subsample::compute_selections(&idx, &layer_dirs, Some(n), None).expect("selections");
let selected_total: usize = selections.iter().map(|s| match s {
None => 0, // shouldn't happen when total_count > n
Some(set) => set.len(),
}).sum();
println!("requested={n} selected_total={selected_total}");
let mut sorted_counts = counts.clone();
sorted_counts.sort_unstable_by(|a, b| b.cmp(a));
println!("top10 counts: {:?}", &sorted_counts[..10.min(sorted_counts.len())]);
}
#[test]
#[ignore]
fn diag_real_index_base_a_representation() {
// Diagnostic for a user-reported observation, reproduced on two
// unrelated real datasets (a plant and this bacterial index):
// `composition_transitions`'s row/column for base 'A' is entirely
// zero, even though 'A' appears at real (~6-22%) frequency in the
// pseudo-alignment itself. A minimal synthetic reproduction
// (`base_pair_tally_accumulates_base_a_diagnostic`) showed the
// pairwise tally mechanism itself works correctly for base A in
// isolation — so this checks one level lower, purely structural
// (annex bits only, no cross-partition/genome-level resolution): does
// base A even show up as *present* in minorant families' own
// `FamilyMask`s at all, at a rate proportional to the other 3 bases?
// If yes, the anomaly is specific to the pairwise/`single_form`
// resolution stage; if `has_base[0]` is itself near-zero here, the
// anomaly originates earlier, in `build_sibling_annex`/`central_base`.
let idx = KmerIndex::open("/Users/coissac/Sync/travail/__MOI__/obikmer/benchmark/global_index_presence")
.expect("open real index");
let layer_dirs = super::family_scan::sibling_layer_dirs(&idx).expect("layer dirs");
let mut has_base = [0u64; 4];
let mut minorant_total = 0u64;
let mut minorant_central_a = 0u64;
for layer_dir in &layer_dirs {
let annex = SiblingAnnex::open(&layer_dir.join(ANNEX_FILE_NAME)).unwrap();
for slot in 0..annex.len() {
let Some(mask) = annex.get(slot) else { continue };
if !mask.is_minorant() {
continue;
}
minorant_total += 1;
for b in 0..4u8 {
if mask.has(b) {
has_base[b as usize] += 1;
}
}
}
// Cross-check: for a sample of this layer's minorants, is the
// slot's *own* central base ever actually 'A' (bit 0)? — reads
// the real k-mer via the layer's MPHF, not just the mask.
let meta = PartitionMeta::load(layer_dir.parent().unwrap()).unwrap();
let mphf = MphfLayer::open(layer_dir, &meta.mode).unwrap();
for (order, kmer) in mphf.enumerate_kmers().take(2_000_000) {
let Some(mask) = annex.get(order) else { continue };
if mask.is_minorant() && central_base(kmer, idx.kmer_size()) == 0 {
minorant_central_a += 1;
}
}
}
println!(
"minorant_total={minorant_total} has_base(A,C,G,T)={has_base:?} minorant_central_a_sampled={minorant_central_a}"
);
}
#[test]
#[ignore]
fn diag_real_index_sankoff_bundle_base_a() {
// Structural check (`diag_real_index_base_a_representation`) shows
// base A fully, proportionally represented at the family level (~23M
// minorants, the *highest* of the four bases) — so the anomaly must
// be downstream, in `sankoff_bundle`'s actual pairwise resolution.
// Reproduces through the real cross-partition code path (not the
// minimal synthetic fixture, which showed the mechanism working in
// isolation), bounded by `--subsample` to stay fast.
let idx = KmerIndex::open("/Users/coissac/Sync/travail/__MOI__/obikmer/benchmark/global_index_presence")
.expect("open real index");
let n_genomes = idx.meta().genomes.len();
let exclude_mask = vec![false; n_genomes];
let bundle = idx.sankoff_bundle(Some(1_000_000), None, 0.5, &exclude_mask).expect("sankoff_bundle");
println!("base_pair_tally.same={:?}", bundle.base_pair_tally.same);
println!("base_pair_tally.counts={:?}", bundle.base_pair_tally.counts);
let alignment_a_count: usize = bundle.alignment.sequences.iter()
.map(|seq| seq.iter().filter(|&&b| b == b'A').count())
.sum();
println!("alignment raw 'A' byte count={alignment_a_count}");
}
#[test]
#[ignore]
fn diag_real_index_genome_mask_for_base_a_families() {
// Directly inspects the raw `genome_mask` array `sankoff_bundle`'s
// closures see, for the first few variable families where base A is
// present — to check whether A ever co-occurs as `single_form` in two
// *different* genomes at the same family at all (the precondition for
// `bp_same`/`bp_counts` to ever increment for base A), rather than
// reasoning about it further.
use std::sync::Arc;
let idx = KmerIndex::open("/Users/coissac/Sync/travail/__MOI__/obikmer/benchmark/global_index_presence")
.expect("open real index");
let n_parts = idx.n_partitions();
let n_genomes = idx.meta().genomes.len();
let with_counts = idx.meta().config.with_counts;
let k = idx.kmer_size();
let n_bits = n_parts.trailing_zeros() as usize;
let partition = obikpartitionner::KmerPartition::open_with_config(
idx.root_path(), idx.kmer_size(), idx.minimizer_size(), n_bits,
).unwrap();
let cache = Arc::new(super::cache::PartitionCache::build(&partition, n_parts, with_counts).unwrap());
let layer_dirs = super::family_scan::sibling_layer_dirs(&idx).unwrap();
let mut printed = 0usize;
for layer_dir in &layer_dirs {
super::family_scan::scan_layer_families(
layer_dir, n_parts, n_genomes, with_counts, k, &cache, &super::family_scan::Selection::All,
|_family_idx, mask, genome_mask| {
if printed >= 15 || !mask.has(0) || mask.family_size() < 2 {
return;
}
let single_a: Vec<usize> = (0..n_genomes).filter(|&g| genome_mask[g] == 1).collect();
let others: Vec<(usize, u8)> = (0..n_genomes)
.filter(|&g| genome_mask[g] != 0 && genome_mask[g] != 1)
.map(|g| (g, genome_mask[g]))
.collect();
println!(
"family bits={:#06b} size={} single_A_genomes={:?} other_nonzero_genomes={:?}",
mask.bits(), mask.family_size(), single_a, others,
);
printed += 1;
},
).unwrap();
if printed >= 15 {
break;
}
}
}
#[test]
#[ignore]
fn diag_plant_index_presence_matrix_sparsity() {
let idx = KmerIndex::open("/Users/coissac/travail/obiskim/data/phyloalps/phyloskims_sal_vac")
.expect("open real plant index");
let layer_dirs = super::family_scan::sibling_layer_dirs(&idx).expect("layer dirs");
println!("n_layers={}", layer_dirs.len());
// Sample a spread of layers rather than just the first few (partition
// order, not necessarily representative) — every 37th layer, capped at 20.
let sample: Vec<&std::path::PathBuf> = layer_dirs.iter().step_by(37).take(20).collect();
let mut total_ones: u64 = 0;
let mut total_cells: u64 = 0;
for layer_dir in &sample {
let mat = obicompactvec::PersistentBitMatrix::open(layer_dir).expect("open presence matrix");
let n = mat.n() as u64;
let n_cols = mat.n_cols() as u64;
let ones: u64 = mat.count_ones().iter().sum();
let cells = n * n_cols;
let density = ones as f64 / cells as f64;
println!(
"{}: n_slots={n} n_cols={n_cols} ones={ones} cells={cells} density={:.6} sparsity={:.6}",
layer_dir.display(), density, 1.0 - density,
);
total_ones += ones;
total_cells += cells;
}
let overall_density = total_ones as f64 / total_cells as f64;
println!(
"OVERALL sampled_layers={} total_ones={total_ones} total_cells={total_cells} density={:.6} sparsity={:.6}",
sample.len(), overall_density, 1.0 - overall_density,
);
}
#[test]
#[ignore]
fn diag_plant_index_multigenome_set_duplication() {
use std::collections::HashMap;
let idx = KmerIndex::open("/Users/coissac/travail/obiskim/data/phyloalps/phyloskims_sal_vac")
.expect("open real plant index");
let layer_dirs = super::family_scan::sibling_layer_dirs(&idx).expect("layer dirs");
// Same spread-sampling strategy as the sparsity diagnostic — every
// 37th layer, capped at 8 (this one is heavier per layer: builds a
// hash map of every distinct multi-genome set).
let sample: Vec<&std::path::PathBuf> = layer_dirs.iter().step_by(37).take(8).collect();
let mut total_rows: u64 = 0;
let mut total_multi_rows: u64 = 0;
let mut total_distinct_multi_sets: u64 = 0;
let mut total_singleton_rows: u64 = 0;
for layer_dir in &sample {
let mat = obicompactvec::PersistentBitMatrix::open(layer_dir).expect("open presence matrix");
let n = mat.n();
let n_cols = mat.n_cols();
let mut seen: HashMap<Vec<u32>, u32> = HashMap::new();
let mut singleton = 0u64;
let mut multi = 0u64;
let mut buf = vec![0u32; n_cols];
for slot in 0..n {
mat.fill_row(slot, &mut buf);
let set: Vec<u32> = (0..n_cols).filter(|&c| buf[c] != 0).map(|c| c as u32).collect();
match set.len() {
0 => {}
1 => singleton += 1,
_ => {
multi += 1;
*seen.entry(set).or_insert(0) += 1;
}
}
}
let distinct = seen.len() as u64;
println!(
"{}: n_slots={n} n_cols={n_cols} singleton_rows={singleton} multi_rows={multi} distinct_multi_sets={distinct} dedup_ratio={:.4}",
layer_dir.display(),
if multi > 0 { distinct as f64 / multi as f64 } else { 0.0 },
);
total_rows += n as u64;
total_singleton_rows += singleton;
total_multi_rows += multi;
total_distinct_multi_sets += distinct;
}
println!(
"OVERALL sampled_layers={} total_rows={total_rows} singleton_rows={total_singleton_rows} multi_rows={total_multi_rows} distinct_multi_sets={total_distinct_multi_sets} dedup_ratio={:.4}",
sample.len(),
if total_multi_rows > 0 { total_distinct_multi_sets as f64 / total_multi_rows as f64 } else { 0.0 },
);
}
#[test]
#[ignore]
fn diag_plant_index_cardinality_distribution() {
let idx = KmerIndex::open("/Users/coissac/travail/obiskim/data/phyloalps/phyloskims_sal_vac")
.expect("open real plant index");
let layer_dirs = super::family_scan::sibling_layer_dirs(&idx).expect("layer dirs");
let sample: Vec<&std::path::PathBuf> = layer_dirs.iter().step_by(37).take(8).collect();
let mut hist: Vec<u64> = Vec::new(); // hist[k] = number of rows with cardinality k
let mut total_rows: u64 = 0;
for layer_dir in &sample {
let mat = obicompactvec::PersistentBitMatrix::open(layer_dir).expect("open presence matrix");
let n = mat.n();
let n_cols = mat.n_cols();
if hist.len() < n_cols + 1 {
hist.resize(n_cols + 1, 0);
}
let mut buf = vec![0u32; n_cols];
for slot in 0..n {
mat.fill_row(slot, &mut buf);
let card = buf.iter().filter(|&&v| v != 0).count();
hist[card] += 1;
}
total_rows += n as u64;
}
println!("total_rows={total_rows}");
let mut cum: u64 = 0;
for (card, &count) in hist.iter().enumerate() {
if count == 0 { continue; }
cum += count;
println!(
"cardinality={card:3} count={count:10} frac={:.6} cumulative_frac={:.6}",
count as f64 / total_rows as f64,
cum as f64 / total_rows as f64,
);
}
}
#[test]
#[ignore]
fn bench_sparse_matrix_inmemory_construction_ram() {
use std::collections::HashMap;
use std::time::Instant;
let idx = KmerIndex::open("/Users/coissac/travail/obiskim/data/phyloalps/phyloskims_sal_vac")
.expect("open real plant index");
let layer_dirs = super::family_scan::sibling_layer_dirs(&idx).expect("layer dirs");
// A `layer_1` (the larger, merged layer) — pick the one already
// profiled: part_00018/index/layer_1, ~30.2M rows.
let layer_dir = layer_dirs.iter()
.find(|p| p.to_string_lossy().contains("part_00018") && p.to_string_lossy().contains("layer_1"))
.expect("expected layer not found — layer_dirs order may have changed");
let mat = obicompactvec::PersistentBitMatrix::open(layer_dir).expect("open presence matrix");
let n = mat.n();
let n_cols = mat.n_cols();
println!("layer={} n_slots={n} n_cols={n_cols}", layer_dir.display());
let rss_before = obisys::peak_rss_bytes();
let t0 = Instant::now();
// Plan design: is_multi flag (1 bit/row) + two SEPARATE arrays
// (singleton-only, multi-only), each sized to its own value range —
// not one array covering the whole dict_id space. In-memory prototype
// (u32 per entry, not yet bit-packed — the fixed-width primitive
// isn't built yet); this benchmark measures RAM + the split's real
// final size, not the construction-time representation's size.
let mut is_multi: Vec<bool> = Vec::with_capacity(n);
let mut singleton_array: Vec<u32> = Vec::new();
let mut multi_array: Vec<u32> = Vec::new();
let mut dedup: HashMap<Vec<u32>, u32> = HashMap::new();
let mut next_dict_id: u32 = 0;
let mut dict_values_bytes: u64 = 0; // running varint-encoded size estimate
let mut buf = vec![0u32; n_cols];
for slot in 0..n {
mat.fill_row(slot, &mut buf);
let set: Vec<u32> = (0..n_cols).filter(|&c| buf[c] != 0).map(|c| c as u32).collect();
match set.len() {
0 => { is_multi.push(false); singleton_array.push(0); } // shouldn't happen on a real built index; placeholder
1 => {
is_multi.push(false);
singleton_array.push(set[0]);
}
_ => {
is_multi.push(true);
let id = if let Some(&id) = dedup.get(&set) {
id
} else {
let id = next_dict_id;
next_dict_id += 1;
// varint size estimate: 1 byte per index < 128 (always
// true here, n_cols=91), matching the plan's encoding.
dict_values_bytes += set.iter().map(|&v| if v < 128 { 1 } else { 2 }).sum::<u64>();
dedup.insert(set, id);
id
};
multi_array.push(id);
}
};
}
let elapsed = t0.elapsed();
let rss_after = obisys::peak_rss_bytes();
let n_distinct_multi = dedup.len() as u64;
let n_singleton = singleton_array.len() as u64;
let n_multi = multi_array.len() as u64;
let is_multi_bytes = (n as u64).div_ceil(8);
let singleton_width_bits = (32 - (n_cols as u32).max(1).leading_zeros()).max(1);
let multi_width_bits = (32 - (n_distinct_multi as u32).max(1).leading_zeros()).max(1);
let singleton_array_bytes = (n_singleton * singleton_width_bits as u64).div_ceil(8);
let multi_array_bytes = (n_multi * multi_width_bits as u64).div_ceil(8);
let split_total = is_multi_bytes + singleton_array_bytes + multi_array_bytes + dict_values_bytes;
let dense_bytes = (n as u64 * n_cols as u64).div_ceil(8);
println!(
"elapsed={:?} rss_before={} rss_after={} rss_delta={}",
elapsed, fmt_mb(rss_before), fmt_mb(rss_after), fmt_mb(rss_after.saturating_sub(rss_before)),
);
println!(
"n_singleton={n_singleton} n_multi={n_multi} n_distinct_multi={n_distinct_multi} \
singleton_width_bits={singleton_width_bits} multi_width_bits={multi_width_bits}",
);
println!(
"is_multi_bytes={} singleton_array_bytes={} multi_array_bytes={} dict_values_bytes={} \
split_total_bytes={} ({:.1}x vs dense) dense_bytes={}",
fmt_mb(is_multi_bytes), fmt_mb(singleton_array_bytes), fmt_mb(multi_array_bytes),
fmt_mb(dict_values_bytes), fmt_mb(split_total),
dense_bytes as f64 / split_total as f64,
fmt_mb(dense_bytes),
);
}
fn fmt_mb(bytes: u64) -> String {
format!("{:.1}MB", bytes as f64 / (1024.0 * 1024.0))
}
#[test]
#[ignore]
fn bench_real_persistent_sparse_bit_matrix_on_disk_size() {
use std::time::Instant;
let idx = KmerIndex::open("/Users/coissac/travail/obiskim/data/phyloalps/phyloskims_sal_vac")
.expect("open real plant index");
let layer_dirs = super::family_scan::sibling_layer_dirs(&idx).expect("layer dirs");
let layer_dir = layer_dirs.iter()
.find(|p| p.to_string_lossy().contains("part_00018") && p.to_string_lossy().contains("layer_1"))
.expect("expected layer not found — layer_dirs order may have changed");
let dense = obicompactvec::PersistentBitMatrix::open(layer_dir).expect("open presence matrix");
let n = dense.n();
let n_cols = dense.n_cols();
println!("layer={} n_slots={n} n_cols={n_cols}", layer_dir.display());
let out_dir = std::env::temp_dir().join(format!("sparse_bench_{}", std::process::id()));
let _ = std::fs::remove_dir_all(&out_dir);
let rss_before = obisys::peak_rss_bytes();
let t0 = Instant::now();
let sparse = obicompactvec::PersistentSparseBitMatrixBuilder::build_from_dense(&dense, &out_dir)
.expect("build_from_dense")
.finish()
.expect("finish");
let elapsed = t0.elapsed();
let rss_after = obisys::peak_rss_bytes();
// Real on-disk size: sum of every file actually written under out_dir.
let mut sparse_bytes: u64 = 0;
for entry in std::fs::read_dir(&out_dir).unwrap() {
let entry = entry.unwrap();
let size = entry.metadata().unwrap().len();
println!(" {}: {}", entry.file_name().to_string_lossy(), fmt_mb(size));
sparse_bytes += size;
}
let dense_bytes = (n as u64 * n_cols as u64).div_ceil(8);
println!(
"elapsed={:?} rss_before={} rss_after={} rss_delta={}",
elapsed, fmt_mb(rss_before), fmt_mb(rss_after), fmt_mb(rss_after.saturating_sub(rss_before)),
);
println!(
"REAL on-disk: sparse_total={} dense={} ({:.1}x)",
fmt_mb(sparse_bytes), fmt_mb(dense_bytes), dense_bytes as f64 / sparse_bytes as f64,
);
// Sanity: reopen from disk (fresh mmap, not the just-built handle) and
// spot-check a handful of rows against the dense original.
drop(sparse);
let reopened = obicompactvec::PersistentSparseBitMatrix::open(&out_dir).expect("reopen");
let mut dense_buf = vec![0u32; n_cols];
for &slot in &[0usize, 1, 1000, n / 2, n - 1] {
dense.fill_row(slot, &mut dense_buf);
let sparse_row = reopened.row(slot);
for c in 0..n_cols {
assert_eq!(sparse_row[c], dense_buf[c] != 0, "slot {slot}, col {c}");
}
}
println!("spot-check rows: OK");
// ── Access-time comparison ──────────────────────────────────────────
// Both patterns matter in practice: sequential (a full-layer scan, the
// existing `scan_layer_families` shape) and random (a single-family
// cross-partition lookup, the entropy/`--shannon` shape). Same slot
// sequence used against both structures for a fair comparison.
const N_ACCESS: usize = 2_000_000;
// Deterministic xorshift, no extra dependency — fine for a benchmark's
// access pattern, not for anything security- or correctness-sensitive.
let mut rng_state: u64 = 0x9E3779B97F4A7C15;
let mut next_rand = move || {
rng_state ^= rng_state << 13;
rng_state ^= rng_state >> 7;
rng_state ^= rng_state << 17;
rng_state
};
let random_slots: Vec<usize> = (0..N_ACCESS).map(|_| (next_rand() as usize) % n).collect();
let sequential_slots: Vec<usize> = (0..N_ACCESS).map(|i| i % n).collect();
let mut buf = vec![0u32; n_cols];
let t = Instant::now();
for &slot in &sequential_slots {
dense.fill_row(slot, &mut buf);
}
let dense_seq = t.elapsed();
let t = Instant::now();
for &slot in &sequential_slots {
reopened.fill_row(slot, &mut buf);
}
let sparse_seq = t.elapsed();
let t = Instant::now();
for &slot in &random_slots {
dense.fill_row(slot, &mut buf);
}
let dense_rand = t.elapsed();
let t = Instant::now();
for &slot in &random_slots {
reopened.fill_row(slot, &mut buf);
}
let sparse_rand = t.elapsed();
println!(
"ACCESS ({N_ACCESS} rows) sequential: dense={:?} ({:.0}ns/row) sparse={:?} ({:.0}ns/row) — {:.2}x",
dense_seq, dense_seq.as_nanos() as f64 / N_ACCESS as f64,
sparse_seq, sparse_seq.as_nanos() as f64 / N_ACCESS as f64,
sparse_seq.as_secs_f64() / dense_seq.as_secs_f64(),
);
println!(
"ACCESS ({N_ACCESS} rows) random: dense={:?} ({:.0}ns/row) sparse={:?} ({:.0}ns/row) — {:.2}x",
dense_rand, dense_rand.as_nanos() as f64 / N_ACCESS as f64,
sparse_rand, sparse_rand.as_nanos() as f64 / N_ACCESS as f64,
sparse_rand.as_secs_f64() / dense_rand.as_secs_f64(),
);
// ── Column-major access, "just for fun" ─────────────────────────────
// The whole point of `DevDocMD/architecture/siblings.md`'s "Explicitly
// deferred" section: dense is genome-major (native, contiguous column
// access), sparse is k-mer-major (no column method at all — reading
// one column means decoding every row and keeping one bit each time).
// Extract one full column (all n rows) both ways.
let col = n_cols / 2;
let t = Instant::now();
let dense_col_view = dense.col_view(col);
let dense_col_ones: u64 = (0..n).filter(|&s| dense_col_view.get(s)).count() as u64;
let dense_col_time = t.elapsed();
let t = Instant::now();
let mut buf2 = vec![0u32; n_cols];
let sparse_col_ones: u64 = (0..n)
.filter(|&s| { reopened.fill_row(s, &mut buf2); buf2[col] != 0 })
.count() as u64;
let sparse_col_time = t.elapsed();
assert_eq!(dense_col_ones, sparse_col_ones, "column {col} popcount disagrees");
println!(
"COLUMN-MAJOR (col {col}, {n} rows) dense={:?} ({:.0}ns/row) sparse={:?} ({:.0}ns/row) — {:.1}x SLOWER on sparse",
dense_col_time, dense_col_time.as_nanos() as f64 / n as f64,
sparse_col_time, sparse_col_time.as_nanos() as f64 / n as f64,
sparse_col_time.as_secs_f64() / dense_col_time.as_secs_f64(),
);
let _ = std::fs::remove_dir_all(&out_dir);
}
#[test]
#[ignore]
fn diag_real_index_verify_kmer_resolution() {
use std::sync::Arc;
let idx = KmerIndex::open("/Users/coissac/Sync/travail/__MOI__/obikmer/benchmark/global_index_presence")
.expect("open real index");
let n_parts = idx.n_partitions();
let n_genomes = idx.meta().genomes.len();
let with_counts = idx.meta().config.with_counts;
let k = idx.kmer_size();
let n_bits = n_parts.trailing_zeros() as usize;
let partition = obikpartitionner::KmerPartition::open_with_config(
idx.root_path(), idx.kmer_size(), idx.minimizer_size(), n_bits,
).unwrap();
let cache = Arc::new(super::cache::PartitionCache::build(&partition, n_parts, with_counts).unwrap());
// The k-mer with A+other (mask 0b0011 = A+C)
let kmer_str = "AGCTAGCTATTGAGCCTGGTCCGTATGAAAC";
let kmer = canonical(kmer_str.as_bytes());
let rc_str = kmer_revcomp_to_string(&kmer, k);
println!("kmer_forward={} kmer_revcomp={}", kmer_str, rc_str);
println!("kmer partition={}", kmer.partition(n_parts));
// Check if this k-mer is found in its own partition
let dest_part = kmer.partition(n_parts);
if let Some(layer) = cache.find(dest_part, kmer) {
println!("FOUND in partition {} layer {}", dest_part, layer);
} else {
println!("NOT FOUND in partition {}", dest_part);
}
// Now check the variant with A at central position
// The mask is 0b0011 = A+C, so there should be a variant with A
// Let's generate the canonical neighbors and check which one has A
let mut found_variants = Vec::new();
for variant in kmer.central_canonical_neighbors() {
if variant == kmer {
continue;
}
let base = central_base(variant, k);
if base == 0 { // A
let v_part = variant.partition(n_parts);
println!("variant with A: partition={} found={}", v_part, cache.find(v_part, variant).is_some());
found_variants.push(variant);
}
}
// Also check the C variant
for variant in kmer.central_canonical_neighbors() {
if variant == kmer {
continue;
}
let base = central_base(variant, k);
if base == 1 { // C
let v_part = variant.partition(n_parts);
println!("variant with C: partition={} found={}", v_part, cache.find(v_part, variant).is_some());
}
}
}
fn kmer_to_string(kmer: &CanonicalKmer, k: usize) -> String {
let mut s = String::with_capacity(k);
for i in 0..k {
let n = kmer.nucleotide(i);
s.push(match n {
0 => 'A',
1 => 'C',
2 => 'G',
3 => 'T',
_ => unreachable!(),
});
}
s
}
fn kmer_revcomp_to_string(kmer: &CanonicalKmer, k: usize) -> String {
let rc = kmer.revcomp();
let mut s = String::with_capacity(k);
for i in 0..k {
let n = rc.nucleotide(i);
s.push(match n {
0 => 'A',
1 => 'C',
2 => 'G',
3 => 'T',
_ => unreachable!(),
});
}
s
}