refactor: extract partition iteration and unify group selection
Move partition iteration logic to obikdump, introducing a FilteredPartitionIter trait over IndexCache for batch-oriented scanning with configurable data retrieval and early termination. Consolidate ingroup and outgroup index storage in GroupQuorumFilter into a unified Selection struct driven by predicate matching. Update dependency manifests to include obikidxcache, rayon, and obikentropy, and remove the deprecated dump_layer module while adjusting public API re-exports.
This commit is contained in:
Generated
+4
@@ -1524,7 +1524,10 @@ name = "obikdump"
|
|||||||
version = "0.1.0"
|
version = "0.1.0"
|
||||||
dependencies = [
|
dependencies = [
|
||||||
"obikfilter",
|
"obikfilter",
|
||||||
|
"obikidxcache",
|
||||||
"obikindex",
|
"obikindex",
|
||||||
|
"obikseq",
|
||||||
|
"rayon",
|
||||||
]
|
]
|
||||||
|
|
||||||
[[package]]
|
[[package]]
|
||||||
@@ -1539,6 +1542,7 @@ name = "obikfilter"
|
|||||||
version = "0.1.0"
|
version = "0.1.0"
|
||||||
dependencies = [
|
dependencies = [
|
||||||
"obicompactvec",
|
"obicompactvec",
|
||||||
|
"obikentropy",
|
||||||
"obikindex",
|
"obikindex",
|
||||||
"obikseq",
|
"obikseq",
|
||||||
"obiskio",
|
"obiskio",
|
||||||
|
|||||||
@@ -4,5 +4,8 @@ version = "0.1.0"
|
|||||||
edition = "2024"
|
edition = "2024"
|
||||||
|
|
||||||
[dependencies]
|
[dependencies]
|
||||||
obikindex = { path = "../obikindex" }
|
obikindex = { path = "../obikindex" }
|
||||||
obikfilter = { path = "../obikfilter" }
|
obikfilter = { path = "../obikfilter" }
|
||||||
|
obikidxcache = { path = "../obikidxcache" }
|
||||||
|
obikseq = { path = "../obikseq" }
|
||||||
|
rayon = "1"
|
||||||
|
|||||||
+44
-31
@@ -5,9 +5,14 @@ use rayon::prelude::*;
|
|||||||
|
|
||||||
use obikindex::{OKIError, OKIResult};
|
use obikindex::{OKIError, OKIResult};
|
||||||
use obikindex::KmerIndex;
|
use obikindex::KmerIndex;
|
||||||
|
use obikidxcache::index_cache::IndexCache;
|
||||||
use obikfilter::KmerFilter;
|
use obikfilter::KmerFilter;
|
||||||
|
|
||||||
impl KmerIndex {
|
use crate::partition_iter::FilteredPartitionIter;
|
||||||
|
|
||||||
|
/// Raw content export of a `KmerIndex` — `KmerIndex` is a foreign type
|
||||||
|
/// (`obikindex`), so this is an extension trait rather than an inherent `impl`.
|
||||||
|
pub trait IndexDump {
|
||||||
/// Write a CSV table of all indexed kmers to `out`.
|
/// Write a CSV table of all indexed kmers to `out`.
|
||||||
///
|
///
|
||||||
/// Columns: `kmer`, then one column per genome (in index order).
|
/// Columns: `kmer`, then one column per genome (in index order).
|
||||||
@@ -17,11 +22,25 @@ impl KmerIndex {
|
|||||||
/// the output uses 0/1 presence columns.
|
/// the output uses 0/1 presence columns.
|
||||||
///
|
///
|
||||||
/// Partitions are scanned in parallel; each partition buffers its output locally
|
/// Partitions are scanned in parallel; each partition buffers its output locally
|
||||||
/// before the main thread writes the chunks in partition order.
|
/// before the main thread writes the chunks in partition order. Each partition's
|
||||||
|
/// layers are cached (`IndexCache`) only for the scan of that one partition —
|
||||||
|
/// `self` is a complete, read-only source index, never a destination.
|
||||||
///
|
///
|
||||||
/// The caller must have set the global kmer length (`obikseq::set_k`) before
|
/// The caller must have set the global kmer length (`obikseq::set_k`) before
|
||||||
/// calling this method.
|
/// calling this method.
|
||||||
pub fn dump<W: Write, F: Fn() + Send + Sync>(
|
fn dump<W: Write, F: Fn() + Send + Sync>(
|
||||||
|
&self,
|
||||||
|
out: &mut W,
|
||||||
|
force_presence: bool,
|
||||||
|
debug: bool,
|
||||||
|
head: Option<usize>,
|
||||||
|
filters: &[Box<dyn KmerFilter>],
|
||||||
|
on_partition: F,
|
||||||
|
) -> OKIResult<()>;
|
||||||
|
}
|
||||||
|
|
||||||
|
impl IndexDump for KmerIndex {
|
||||||
|
fn dump<W: Write, F: Fn() + Send + Sync>(
|
||||||
&self,
|
&self,
|
||||||
out: &mut W,
|
out: &mut W,
|
||||||
force_presence: bool,
|
force_presence: bool,
|
||||||
@@ -30,8 +49,8 @@ impl KmerIndex {
|
|||||||
filters: &[Box<dyn KmerFilter>],
|
filters: &[Box<dyn KmerFilter>],
|
||||||
on_partition: F,
|
on_partition: F,
|
||||||
) -> OKIResult<()> {
|
) -> OKIResult<()> {
|
||||||
let genomes = self.meta.genomes().map_err(OKIError::Io)?;
|
let genomes = self.meta().genomes().map_err(OKIError::Io)?;
|
||||||
let use_counts = self.meta.config.with_counts && !force_presence;
|
let use_counts = self.meta().config.with_counts && !force_presence;
|
||||||
let n_genomes = genomes.len().max(1);
|
let n_genomes = genomes.len().max(1);
|
||||||
let kmer_size = self.kmer_size();
|
let kmer_size = self.kmer_size();
|
||||||
|
|
||||||
@@ -68,20 +87,17 @@ impl KmerIndex {
|
|||||||
Ok(_) => { write_row(buf, row, prefix); true }
|
Ok(_) => { write_row(buf, row, prefix); true }
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
let cache = IndexCache::new(self, Some(vec![i]));
|
||||||
if debug {
|
if debug {
|
||||||
self
|
cache.iter_partition_kmers_located(i, use_counts, n_genomes, filters, |part, layer, kmer, row| {
|
||||||
.iter_partition_kmers_located(i, use_counts, n_genomes, filters, |part, layer, kmer, row| {
|
let seq = String::from_utf8(kmer.to_ascii()).unwrap_or_else(|_| "?".repeat(kmer_size));
|
||||||
let seq = String::from_utf8(kmer.to_ascii()).unwrap_or_else(|_| "?".repeat(kmer_size));
|
try_write(&mut buf, &row, &format!("{part},{layer},{seq}"))
|
||||||
try_write(&mut buf, &row, &format!("{part},{layer},{seq}"))
|
})?;
|
||||||
})
|
|
||||||
.map_err(OKIError::Partition)?;
|
|
||||||
} else {
|
} else {
|
||||||
self
|
cache.iter_partition_kmers(i, use_counts, n_genomes, filters, |kmer, row| {
|
||||||
.iter_partition_kmers(i, use_counts, n_genomes, filters, |kmer, row| {
|
let seq = String::from_utf8(kmer.to_ascii()).unwrap_or_else(|_| "?".repeat(kmer_size));
|
||||||
let seq = String::from_utf8(kmer.to_ascii()).unwrap_or_else(|_| "?".repeat(kmer_size));
|
try_write(&mut buf, &row, &seq)
|
||||||
try_write(&mut buf, &row, &seq)
|
})?;
|
||||||
})
|
|
||||||
.map_err(OKIError::Partition)?;
|
|
||||||
}
|
}
|
||||||
on_partition();
|
on_partition();
|
||||||
Ok(buf)
|
Ok(buf)
|
||||||
@@ -90,22 +106,19 @@ impl KmerIndex {
|
|||||||
// ── Unbounded: no atomic, no contention ───────────────────────────
|
// ── Unbounded: no atomic, no contention ───────────────────────────
|
||||||
(0..n).into_par_iter().map(|i| {
|
(0..n).into_par_iter().map(|i| {
|
||||||
let mut buf = Vec::<u8>::new();
|
let mut buf = Vec::<u8>::new();
|
||||||
|
let cache = IndexCache::new(self, Some(vec![i]));
|
||||||
if debug {
|
if debug {
|
||||||
self
|
cache.iter_partition_kmers_located(i, use_counts, n_genomes, filters, |part, layer, kmer, row| {
|
||||||
.iter_partition_kmers_located(i, use_counts, n_genomes, filters, |part, layer, kmer, row| {
|
let seq = String::from_utf8(kmer.to_ascii()).unwrap_or_else(|_| "?".repeat(kmer_size));
|
||||||
let seq = String::from_utf8(kmer.to_ascii()).unwrap_or_else(|_| "?".repeat(kmer_size));
|
write_row(&mut buf, &row, &format!("{part},{layer},{seq}"));
|
||||||
write_row(&mut buf, &row, &format!("{part},{layer},{seq}"));
|
true
|
||||||
true
|
})?;
|
||||||
})
|
|
||||||
.map_err(OKIError::Partition)?;
|
|
||||||
} else {
|
} else {
|
||||||
self
|
cache.iter_partition_kmers(i, use_counts, n_genomes, filters, |kmer, row| {
|
||||||
.iter_partition_kmers(i, use_counts, n_genomes, filters, |kmer, row| {
|
let seq = String::from_utf8(kmer.to_ascii()).unwrap_or_else(|_| "?".repeat(kmer_size));
|
||||||
let seq = String::from_utf8(kmer.to_ascii()).unwrap_or_else(|_| "?".repeat(kmer_size));
|
write_row(&mut buf, &row, &seq);
|
||||||
write_row(&mut buf, &row, &seq);
|
true
|
||||||
true
|
})?;
|
||||||
})
|
|
||||||
.map_err(OKIError::Partition)?;
|
|
||||||
}
|
}
|
||||||
on_partition();
|
on_partition();
|
||||||
Ok(buf)
|
Ok(buf)
|
||||||
|
|||||||
@@ -6,3 +6,6 @@
|
|||||||
//! reverse), same pattern as `obikindexer`/`obikquery`.
|
//! reverse), same pattern as `obikindexer`/`obikquery`.
|
||||||
|
|
||||||
mod dump;
|
mod dump;
|
||||||
|
mod partition_iter;
|
||||||
|
|
||||||
|
pub use partition_iter::FilteredPartitionIter;
|
||||||
|
|||||||
@@ -0,0 +1,126 @@
|
|||||||
|
//! Filtered, batch-oriented iteration over an already-cached index's
|
||||||
|
//! partitions/layers — the read side of `obikfilter`'s `KmerFilter`s.
|
||||||
|
//! `IndexCache` is a foreign type (`obikidxcache`), so this is an extension
|
||||||
|
//! trait rather than an inherent `impl`.
|
||||||
|
//!
|
||||||
|
//! Only meaningful on a *complete* source index: `IndexCache` panics if a
|
||||||
|
//! layer is missing, which a finished index never has. Never use this on a
|
||||||
|
//! destination index still being built (see `obikmerge::partition_merge`,
|
||||||
|
//! which follows the same source-only-cache rule).
|
||||||
|
|
||||||
|
use obikindex::OKIResult;
|
||||||
|
use obikindex::layer::{KmerLayer, LayerContent};
|
||||||
|
use obikidxcache::index_cache::IndexCache;
|
||||||
|
use obikseq::CanonicalKmer;
|
||||||
|
|
||||||
|
use obikfilter::{KmerFilter, passes_all};
|
||||||
|
|
||||||
|
/// Kmers pulled per batch from a layer before filtering — keeps matrix reads
|
||||||
|
/// grouped by (partition, layer) for locality instead of hopping row to row
|
||||||
|
/// across the index. Same convention as `obikphylo::siblings::build`.
|
||||||
|
const BATCH_SIZE: usize = 32768;
|
||||||
|
|
||||||
|
pub trait FilteredPartitionIter {
|
||||||
|
/// Iterate all indexed kmers in partition `part`, calling `cb(kmer, row)` for each
|
||||||
|
/// kmer that passes every filter in `filters`.
|
||||||
|
///
|
||||||
|
/// `use_counts = true` → reads count columns (u32 values per genome), only
|
||||||
|
/// meaningful for `Count` layers. `use_counts = false` → reads presence
|
||||||
|
/// columns, converted to 0/1 u32 (works for both `Count` and `Presence`
|
||||||
|
/// layers — counts collapse to presence).
|
||||||
|
///
|
||||||
|
/// Returns `Ok(true)` if all kmers were visited, `Ok(false)` if the callback halted.
|
||||||
|
fn iter_partition_kmers(
|
||||||
|
&self,
|
||||||
|
part: usize,
|
||||||
|
use_counts: bool,
|
||||||
|
n_genomes: usize,
|
||||||
|
filters: &[Box<dyn KmerFilter>],
|
||||||
|
cb: impl FnMut(CanonicalKmer, Box<[u32]>) -> bool,
|
||||||
|
) -> OKIResult<bool>;
|
||||||
|
|
||||||
|
/// Like [`iter_partition_kmers`](Self::iter_partition_kmers) but the callback
|
||||||
|
/// also receives `(partition, layer)` indices, enabling debug output that
|
||||||
|
/// identifies where each kmer was stored.
|
||||||
|
fn iter_partition_kmers_located(
|
||||||
|
&self,
|
||||||
|
part: usize,
|
||||||
|
use_counts: bool,
|
||||||
|
n_genomes: usize,
|
||||||
|
filters: &[Box<dyn KmerFilter>],
|
||||||
|
cb: impl FnMut(usize, usize, CanonicalKmer, Box<[u32]>) -> bool,
|
||||||
|
) -> OKIResult<bool>;
|
||||||
|
}
|
||||||
|
|
||||||
|
impl FilteredPartitionIter for IndexCache<'_> {
|
||||||
|
fn iter_partition_kmers(
|
||||||
|
&self,
|
||||||
|
part: usize,
|
||||||
|
use_counts: bool,
|
||||||
|
n_genomes: usize,
|
||||||
|
filters: &[Box<dyn KmerFilter>],
|
||||||
|
mut cb: impl FnMut(CanonicalKmer, Box<[u32]>) -> bool,
|
||||||
|
) -> OKIResult<bool> {
|
||||||
|
for l in 0..self.n_layer(part).unwrap_or(0) {
|
||||||
|
let layer = self.get_layer(part, l).expect("layer within n_layer(part)");
|
||||||
|
if !iter_layer_kmers(layer, use_counts, n_genomes, filters, &mut |kmer, row| cb(kmer, row))? {
|
||||||
|
return Ok(false);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
Ok(true)
|
||||||
|
}
|
||||||
|
|
||||||
|
fn iter_partition_kmers_located(
|
||||||
|
&self,
|
||||||
|
part: usize,
|
||||||
|
use_counts: bool,
|
||||||
|
n_genomes: usize,
|
||||||
|
filters: &[Box<dyn KmerFilter>],
|
||||||
|
mut cb: impl FnMut(usize, usize, CanonicalKmer, Box<[u32]>) -> bool,
|
||||||
|
) -> OKIResult<bool> {
|
||||||
|
for l in 0..self.n_layer(part).unwrap_or(0) {
|
||||||
|
let layer = self.get_layer(part, l).expect("layer within n_layer(part)");
|
||||||
|
if !iter_layer_kmers(layer, use_counts, n_genomes, filters, &mut |kmer, row| cb(part, l, kmer, row))? {
|
||||||
|
return Ok(false);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
Ok(true)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/// Batch-and-transpose one layer's kmers into per-kmer filtered rows.
|
||||||
|
/// Returns `Ok(false)` if `cb` asked to stop early.
|
||||||
|
fn iter_layer_kmers(
|
||||||
|
layer: &KmerLayer,
|
||||||
|
use_counts: bool,
|
||||||
|
n_genomes: usize,
|
||||||
|
filters: &[Box<dyn KmerFilter>],
|
||||||
|
cb: &mut dyn FnMut(CanonicalKmer, Box<[u32]>) -> bool,
|
||||||
|
) -> OKIResult<bool> {
|
||||||
|
let read_counts = use_counts && matches!(layer.content(), LayerContent::Count);
|
||||||
|
|
||||||
|
for kmers in layer.iter_kmers_batch(BATCH_SIZE) {
|
||||||
|
// Kmers come straight from this layer's own iterator, so every one
|
||||||
|
// is a guaranteed member — a raw hash is enough, no membership
|
||||||
|
// recheck needed (see `KmerLayer::hash_batch`'s own doc).
|
||||||
|
let slots = layer.hash_batch(&kmers);
|
||||||
|
|
||||||
|
let cols: Vec<Vec<u32>> = if read_counts {
|
||||||
|
let mut cols: Vec<Vec<u32>> = vec![Vec::new(); n_genomes];
|
||||||
|
layer.fill_sub_matrix(&slots, &mut cols);
|
||||||
|
cols
|
||||||
|
} else {
|
||||||
|
let mut bool_cols: Vec<Vec<bool>> = vec![Vec::new(); n_genomes];
|
||||||
|
layer.fill_sub_matrix_carries(&slots, &mut bool_cols);
|
||||||
|
bool_cols.iter().map(|c| c.iter().map(|&b| b as u32).collect()).collect()
|
||||||
|
};
|
||||||
|
|
||||||
|
for (i, kmer) in kmers.into_iter().enumerate() {
|
||||||
|
let row: Box<[u32]> = cols.iter().map(|c| c[i]).collect();
|
||||||
|
if passes_all(filters, kmer, &row, n_genomes) && !cb(kmer, row) {
|
||||||
|
return Ok(false);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
Ok(true)
|
||||||
|
}
|
||||||
@@ -9,3 +9,4 @@ obicompactvec = { path = "../obicompactvec" }
|
|||||||
obikseq = { path = "../obikseq" }
|
obikseq = { path = "../obikseq" }
|
||||||
obiskio = { path = "../obiskio" }
|
obiskio = { path = "../obiskio" }
|
||||||
obitaxonomy = { path = "../obitaxonomy" }
|
obitaxonomy = { path = "../obitaxonomy" }
|
||||||
|
obikentropy = { path = "../obikentropy" }
|
||||||
|
|||||||
@@ -1,192 +0,0 @@
|
|||||||
use obikindex::layer::MphfLayer;
|
|
||||||
use obicompactvec::{PersistentBitMatrix, PersistentCompactIntMatrix};
|
|
||||||
use obikseq::CanonicalKmer;
|
|
||||||
use obiskio::UnitigFileReader;
|
|
||||||
use obikindex::{OKIError, OKIResult};
|
|
||||||
|
|
||||||
use crate::filter::{KmerFilter, passes_all};
|
|
||||||
use obikindex::KmerIndex;
|
|
||||||
|
|
||||||
impl KmerIndex {
|
|
||||||
/// Iterate all indexed kmers in partition `part`, calling `cb(kmer, row)` for each
|
|
||||||
/// kmer that passes every filter in `filters`.
|
|
||||||
///
|
|
||||||
/// `use_counts = true` → reads count columns (u32 values per genome).
|
|
||||||
/// `use_counts = false` → reads presence columns, converted to 0/1 u32.
|
|
||||||
///
|
|
||||||
/// If no data matrix exists for a layer (pure set-membership, single genome),
|
|
||||||
/// a row of `n_genomes` ones is emitted for every kmer in that layer — unless
|
|
||||||
/// the filter rejects it, in which case the whole layer is skipped.
|
|
||||||
/// Like [`iter_partition_kmers`] but the callback returns `false` to stop early.
|
|
||||||
/// Returns `Ok(true)` if all kmers were visited, `Ok(false)` if the callback halted.
|
|
||||||
pub fn iter_partition_kmers(
|
|
||||||
&self,
|
|
||||||
part: usize,
|
|
||||||
use_counts: bool,
|
|
||||||
n_genomes: usize,
|
|
||||||
filters: &[Box<dyn KmerFilter>],
|
|
||||||
mut cb: impl FnMut(CanonicalKmer, Box<[u32]>) -> bool,
|
|
||||||
) -> OKIResult<bool> {
|
|
||||||
let index_dir = self.index_dir(part);
|
|
||||||
if !index_dir.exists() {
|
|
||||||
return Ok(true);
|
|
||||||
}
|
|
||||||
|
|
||||||
let mut l = 0;
|
|
||||||
loop {
|
|
||||||
let layer_dir = self.layer_dir(part, l)?;
|
|
||||||
if !layer_dir.exists() {
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
l += 1;
|
|
||||||
let mphf = MphfLayer::open(&layer_dir)?;
|
|
||||||
let reader = UnitigFileReader::open_sequential(&layer_dir.join("unitigs.bin"))?;
|
|
||||||
|
|
||||||
let counts_dir = layer_dir.join("counts");
|
|
||||||
let presence_dir = layer_dir.join("presence");
|
|
||||||
|
|
||||||
let cont = if use_counts && counts_dir.exists() {
|
|
||||||
let mat = PersistentCompactIntMatrix::open(&layer_dir).map_err(OKIError::Io)?;
|
|
||||||
let mut cont = true;
|
|
||||||
for (kmer, _, _) in reader.iter_indexed_canonical_kmers() {
|
|
||||||
if let Some(slot) = mphf.find(kmer) {
|
|
||||||
let row = mat.row(slot);
|
|
||||||
if passes_all(filters, kmer, &row, n_genomes) {
|
|
||||||
cont = cb(kmer, row);
|
|
||||||
if !cont {
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
cont
|
|
||||||
} else if !use_counts && presence_dir.exists() {
|
|
||||||
let mat = PersistentBitMatrix::open(&layer_dir).map_err(OKIError::Io)?;
|
|
||||||
let mut cont = true;
|
|
||||||
for (kmer, _, _) in reader.iter_indexed_canonical_kmers() {
|
|
||||||
if let Some(slot) = mphf.find(kmer) {
|
|
||||||
let row: Box<[u32]> = mat.row(slot).iter().map(|&b| b as u32).collect();
|
|
||||||
if passes_all(filters, kmer, &row, n_genomes) {
|
|
||||||
cont = cb(kmer, row);
|
|
||||||
if !cont {
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
cont
|
|
||||||
} else {
|
|
||||||
// No data matrix: implicit presence — all values are 1. `row`
|
|
||||||
// is identical for every kmer, but a filter can still depend
|
|
||||||
// on the kmer's own sequence (e.g. MinComplexity), so this
|
|
||||||
// cannot be evaluated once for the whole layer — filters must
|
|
||||||
// still be tested per kmer.
|
|
||||||
let all_present: Box<[u32]> = vec![1u32; n_genomes].into();
|
|
||||||
let mut cont = true;
|
|
||||||
for (kmer, _, _) in reader.iter_indexed_canonical_kmers() {
|
|
||||||
if mphf.find(kmer).is_some()
|
|
||||||
&& passes_all(filters, kmer, &all_present, n_genomes)
|
|
||||||
{
|
|
||||||
cont = cb(kmer, all_present.clone());
|
|
||||||
if !cont {
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
cont
|
|
||||||
};
|
|
||||||
|
|
||||||
if !cont {
|
|
||||||
return Ok(false);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
Ok(true)
|
|
||||||
}
|
|
||||||
|
|
||||||
/// Like [`iter_partition_kmers`] but the callback also receives `(partition, layer)`
|
|
||||||
/// indices, enabling debug output that identifies where each kmer was stored.
|
|
||||||
/// Returns `Ok(true)` if all kmers were visited, `Ok(false)` if the callback halted.
|
|
||||||
pub fn iter_partition_kmers_located(
|
|
||||||
&self,
|
|
||||||
part: usize,
|
|
||||||
use_counts: bool,
|
|
||||||
n_genomes: usize,
|
|
||||||
filters: &[Box<dyn KmerFilter>],
|
|
||||||
mut cb: impl FnMut(usize, usize, CanonicalKmer, Box<[u32]>) -> bool,
|
|
||||||
) -> OKIResult<bool> {
|
|
||||||
let index_dir = self.index_dir(part);
|
|
||||||
if !index_dir.exists() {
|
|
||||||
return Ok(true);
|
|
||||||
}
|
|
||||||
|
|
||||||
let mut layer = 0;
|
|
||||||
loop {
|
|
||||||
let layer_dir = self.layer_dir(part, layer);
|
|
||||||
if !layer_dir.exists() {
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
let mphf = MphfLayer::open(&layer_dir)?;
|
|
||||||
let reader = UnitigFileReader::open_sequential(&layer_dir.join("unitigs.bin"))?;
|
|
||||||
|
|
||||||
let counts_dir = layer_dir.join("counts");
|
|
||||||
let presence_dir = layer_dir.join("presence");
|
|
||||||
|
|
||||||
let cont = if use_counts && counts_dir.exists() {
|
|
||||||
let mat = PersistentCompactIntMatrix::open(&layer_dir).map_err(OKIError::Io)?;
|
|
||||||
let mut cont = true;
|
|
||||||
for (kmer, _, _) in reader.iter_indexed_canonical_kmers() {
|
|
||||||
if let Some(slot) = mphf.find(kmer) {
|
|
||||||
let row = mat.row(slot);
|
|
||||||
if passes_all(filters, kmer, &row, n_genomes) {
|
|
||||||
cont = cb(part, layer, kmer, row);
|
|
||||||
if !cont {
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
cont
|
|
||||||
} else if !use_counts && presence_dir.exists() {
|
|
||||||
let mat = PersistentBitMatrix::open(&layer_dir).map_err(OKIError::Io)?;
|
|
||||||
let mut cont = true;
|
|
||||||
for (kmer, _, _) in reader.iter_indexed_canonical_kmers() {
|
|
||||||
if let Some(slot) = mphf.find(kmer) {
|
|
||||||
let row: Box<[u32]> = mat.row(slot).iter().map(|&b| b as u32).collect();
|
|
||||||
if passes_all(filters, kmer, &row, n_genomes) {
|
|
||||||
cont = cb(part, layer, kmer, row);
|
|
||||||
if !cont {
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
cont
|
|
||||||
} else {
|
|
||||||
// Same as iter_partition_kmers: row is constant but a filter
|
|
||||||
// may still depend on the kmer's own sequence, so this must
|
|
||||||
// be tested per kmer, not once for the whole layer.
|
|
||||||
let all_present: Box<[u32]> = vec![1u32; n_genomes].into();
|
|
||||||
let mut cont = true;
|
|
||||||
for (kmer, _, _) in reader.iter_indexed_canonical_kmers() {
|
|
||||||
if mphf.find(kmer).is_some()
|
|
||||||
&& passes_all(filters, kmer, &all_present, n_genomes)
|
|
||||||
{
|
|
||||||
cont = cb(part, layer, kmer, all_present.clone());
|
|
||||||
if !cont {
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
cont
|
|
||||||
};
|
|
||||||
|
|
||||||
if !cont {
|
|
||||||
return Ok(false);
|
|
||||||
}
|
|
||||||
layer += 1;
|
|
||||||
}
|
|
||||||
|
|
||||||
Ok(true)
|
|
||||||
}
|
|
||||||
}
|
|
||||||
@@ -1,6 +1,8 @@
|
|||||||
use obicompactvec::FilterMask;
|
use obicompactvec::FilterMask;
|
||||||
use obikseq::CanonicalKmer;
|
use obikseq::CanonicalKmer;
|
||||||
|
|
||||||
|
use crate::predicate::Selection;
|
||||||
|
|
||||||
/// Trait for kmer filters.
|
/// Trait for kmer filters.
|
||||||
///
|
///
|
||||||
/// `kmer` is the k-mer's own canonical sequence, reconstructed from the
|
/// `kmer` is the k-mer's own canonical sequence, reconstructed from the
|
||||||
@@ -173,14 +175,12 @@ impl KmerFilter for MaxTotalCount {
|
|||||||
|
|
||||||
// ── Group-based quorum filter ─────────────────────────────────────────────────
|
// ── Group-based quorum filter ─────────────────────────────────────────────────
|
||||||
|
|
||||||
/// Quorum filter operating on pre-classified genome groups.
|
/// Quorum filter operating on a pre-classified genome [`Selection`].
|
||||||
///
|
///
|
||||||
/// `ingroup_idx` / `outgroup_idx` are column indices into the per-genome row.
|
/// `selection.ingroup_idx` / `selection.outgroup_idx` are column indices into
|
||||||
/// When `ingroup_idx` is empty, no ingroup quorum is checked.
|
/// the per-genome row. When empty, the corresponding quorum is not checked.
|
||||||
/// When `outgroup_idx` is empty, no outgroup quorum is checked.
|
|
||||||
pub struct GroupQuorumFilter {
|
pub struct GroupQuorumFilter {
|
||||||
pub ingroup_idx: Vec<usize>,
|
pub selection: Selection,
|
||||||
pub outgroup_idx: Vec<usize>,
|
|
||||||
pub threshold: u32,
|
pub threshold: u32,
|
||||||
pub min_count: usize,
|
pub min_count: usize,
|
||||||
pub max_count: usize,
|
pub max_count: usize,
|
||||||
@@ -225,22 +225,22 @@ impl GroupQuorumFilter {
|
|||||||
|
|
||||||
impl KmerFilter for GroupQuorumFilter {
|
impl KmerFilter for GroupQuorumFilter {
|
||||||
fn passes(&self, _kmer: CanonicalKmer, row: &[u32], _n_genomes: usize) -> bool {
|
fn passes(&self, _kmer: CanonicalKmer, row: &[u32], _n_genomes: usize) -> bool {
|
||||||
if !self.ingroup_idx.is_empty() {
|
if !self.selection.ingroup_idx.is_empty() {
|
||||||
let n = self.ingroup_idx.iter()
|
let n = self.selection.ingroup_idx.iter()
|
||||||
.filter(|&&i| row.get(i).copied().unwrap_or(0) > self.threshold)
|
.filter(|&&i| row.get(i).copied().unwrap_or(0) > self.threshold)
|
||||||
.count();
|
.count();
|
||||||
let denom = self.ingroup_idx.len();
|
let denom = self.selection.ingroup_idx.len();
|
||||||
if n < self.min_count { return false; }
|
if n < self.min_count { return false; }
|
||||||
if n > self.max_count { return false; }
|
if n > self.max_count { return false; }
|
||||||
let frac = n as f64 / denom as f64;
|
let frac = n as f64 / denom as f64;
|
||||||
if frac < self.min_frac { return false; }
|
if frac < self.min_frac { return false; }
|
||||||
if frac > self.max_frac { return false; }
|
if frac > self.max_frac { return false; }
|
||||||
}
|
}
|
||||||
if !self.outgroup_idx.is_empty() {
|
if !self.selection.outgroup_idx.is_empty() {
|
||||||
let n = self.outgroup_idx.iter()
|
let n = self.selection.outgroup_idx.iter()
|
||||||
.filter(|&&i| row.get(i).copied().unwrap_or(0) > self.threshold)
|
.filter(|&&i| row.get(i).copied().unwrap_or(0) > self.threshold)
|
||||||
.count();
|
.count();
|
||||||
let denom = self.outgroup_idx.len();
|
let denom = self.selection.outgroup_idx.len();
|
||||||
if n < self.min_outgroup_count { return false; }
|
if n < self.min_outgroup_count { return false; }
|
||||||
if n > self.max_outgroup_count { return false; }
|
if n > self.max_outgroup_count { return false; }
|
||||||
let frac = n as f64 / denom as f64;
|
let frac = n as f64 / denom as f64;
|
||||||
@@ -253,17 +253,17 @@ impl KmerFilter for GroupQuorumFilter {
|
|||||||
fn column_mask_expr(&self, _n_genomes: usize) -> Option<FilterMask> {
|
fn column_mask_expr(&self, _n_genomes: usize) -> Option<FilterMask> {
|
||||||
let t = self.threshold.checked_add(1)?;
|
let t = self.threshold.checked_add(1)?;
|
||||||
let mut parts: Vec<FilterMask> = Vec::new();
|
let mut parts: Vec<FilterMask> = Vec::new();
|
||||||
if !self.ingroup_idx.is_empty() {
|
if !self.selection.ingroup_idx.is_empty() {
|
||||||
Self::group_mask_parts(
|
Self::group_mask_parts(
|
||||||
&self.ingroup_idx, t,
|
&self.selection.ingroup_idx, t,
|
||||||
self.min_count, self.max_count,
|
self.min_count, self.max_count,
|
||||||
self.min_frac, self.max_frac,
|
self.min_frac, self.max_frac,
|
||||||
&mut parts,
|
&mut parts,
|
||||||
);
|
);
|
||||||
}
|
}
|
||||||
if !self.outgroup_idx.is_empty() {
|
if !self.selection.outgroup_idx.is_empty() {
|
||||||
Self::group_mask_parts(
|
Self::group_mask_parts(
|
||||||
&self.outgroup_idx, t,
|
&self.selection.outgroup_idx, t,
|
||||||
self.min_outgroup_count, self.max_outgroup_count,
|
self.min_outgroup_count, self.max_outgroup_count,
|
||||||
self.min_outgroup_frac, self.max_outgroup_frac,
|
self.min_outgroup_frac, self.max_outgroup_frac,
|
||||||
&mut parts,
|
&mut parts,
|
||||||
|
|||||||
@@ -4,19 +4,18 @@
|
|||||||
//! operates on already-retained k-mers.
|
//! operates on already-retained k-mers.
|
||||||
//!
|
//!
|
||||||
//! [`filter`] (the `KmerFilter` trait + its implementations) depends only
|
//! [`filter`] (the `KmerFilter` trait + its implementations) depends only
|
||||||
//! on `obicompactvec`/`obikseq`, not on `obikindex`. [`dump_layer`] is the
|
//! on `obicompactvec`/`obikseq`, not on `obikindex`. The partition/layer
|
||||||
//! extension over `obikindex::KmerIndex` that actually iterates a
|
//! iteration that actually runs these filters over an index
|
||||||
//! partition's k-mers through those filters (`iter_partition_kmers`,
|
//! (`iter_partition_kmers`, `iter_partition_kmers_located`) lives in
|
||||||
//! `iter_partition_kmers_located`) — every k-mer, filtered or not, goes
|
//! `obikdump` instead (`FilteredPartitionIter`) — it needs `obikidxcache`'s
|
||||||
//! through this same path (`passes_all` on an empty filter list is always
|
//! `IndexCache` to read a *complete* source index, which this crate has no
|
||||||
//! `true`), so this crate depends one-way on `obikindex`, not the reverse.
|
//! reason to depend on.
|
||||||
|
|
||||||
mod filter;
|
mod filter;
|
||||||
mod dump_layer;
|
|
||||||
mod predicate;
|
mod predicate;
|
||||||
|
|
||||||
pub use filter::{
|
pub use filter::{
|
||||||
GroupQuorumFilter, KmerFilter, MaxGenomeCount, MaxGenomeFraction, MaxTotalCount,
|
GroupQuorumFilter, KmerFilter, MaxGenomeCount, MaxGenomeFraction, MaxTotalCount,
|
||||||
MinComplexity, MinGenomeCount, MinGenomeFraction, MinTotalCount, passes_all,
|
MinComplexity, MinGenomeCount, MinGenomeFraction, MinTotalCount, passes_all,
|
||||||
};
|
};
|
||||||
pub use predicate::{GroupFilterParams, MetaPred};
|
pub use predicate::{GenomeSelector, GroupFilterParams, MetaPred, Selection};
|
||||||
|
|||||||
@@ -68,14 +68,6 @@ impl MetaPred {
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
impl GenomeInfo {
|
|
||||||
/// Evaluate a single metadata predicate against this genome.
|
|
||||||
/// Returns `None` when the predicate's key is absent (NA propagation).
|
|
||||||
pub fn matches(&self, pred: &MetaPred) -> Option<bool> {
|
|
||||||
pred.eval(&self.meta)
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// ── Path matching ─────────────────────────────────────────────────────────────
|
// ── Path matching ─────────────────────────────────────────────────────────────
|
||||||
|
|
||||||
/// True if the stored taxonomy `value` matches `pattern`.
|
/// True if the stored taxonomy `value` matches `pattern`.
|
||||||
@@ -153,19 +145,49 @@ pub struct GroupFilterParams {
|
|||||||
pub max_outgroup_frac: Option<f64>,
|
pub max_outgroup_frac: Option<f64>,
|
||||||
}
|
}
|
||||||
|
|
||||||
impl IndexMeta {
|
// ── Genome selector ──────────────────────────────────────────────────────────
|
||||||
/// Returns indices of genomes matching `pred_str` (single predicate).
|
|
||||||
pub fn matching_genome_indices(&self, pred_str: &str) -> Result<Vec<usize>, String> {
|
/// Result of running a [`GenomeSelector`] against an index's metadata.
|
||||||
let pred = MetaPred::parse(pred_str)?;
|
pub struct Selection {
|
||||||
let genomes = self.genomes().map_err(|e| e.to_string())?;
|
pub ingroup_idx: Vec<usize>,
|
||||||
Ok(genomes.iter().enumerate()
|
pub outgroup_idx: Vec<usize>,
|
||||||
.filter_map(|(i, g)| {
|
}
|
||||||
if g.matches(&pred) == Some(true) { Some(i) } else { std::option::Option::None }
|
|
||||||
})
|
pub struct GenomeSelector {
|
||||||
.collect())
|
pub(crate) ingroup: Vec<MetaPred>,
|
||||||
|
pub(crate) outgroup: Vec<MetaPred>,
|
||||||
|
}
|
||||||
|
|
||||||
|
impl GenomeSelector {
|
||||||
|
/// Parse ingroup (AND'd) and outgroup (OR'd) predicate strings.
|
||||||
|
pub fn parse(ingroup: &[String], outgroup: &[String]) -> Result<Self, String> {
|
||||||
|
let ingroup = ingroup.iter().map(|s| MetaPred::parse(s)).collect::<Result<Vec<_>, _>>()?;
|
||||||
|
let outgroup = outgroup.iter().map(|s| MetaPred::parse(s)).collect::<Result<Vec<_>, _>>()?;
|
||||||
|
Ok(Self { ingroup, outgroup })
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Build a `GroupQuorumFilter` from parsed predicates, evaluated against `self.genomes`.
|
/// Classify `meta`'s genomes into ingroup/outgroup indices.
|
||||||
|
///
|
||||||
|
/// - No predicates at all: every genome is (implicitly) ingroup.
|
||||||
|
/// - Otherwise: ingroup wins on overlap; uncategorized genomes are dropped.
|
||||||
|
pub fn run(&self, meta: &IndexMeta) -> Result<Selection, String> {
|
||||||
|
let genomes = meta.genomes().map_err(|e| e.to_string())?;
|
||||||
|
|
||||||
|
if self.ingroup.is_empty() && self.outgroup.is_empty() {
|
||||||
|
return Ok(Selection { ingroup_idx: (0..genomes.len()).collect(), outgroup_idx: vec![] });
|
||||||
|
}
|
||||||
|
|
||||||
|
let members = classify(&genomes, &self.ingroup, &self.outgroup);
|
||||||
|
let ingroup_idx: Vec<usize> = members.iter().enumerate()
|
||||||
|
.filter(|(_, m)| matches!(m, Membership::Ingroup))
|
||||||
|
.map(|(i, _)| i).collect();
|
||||||
|
let outgroup_idx: Vec<usize> = members.iter().enumerate()
|
||||||
|
.filter(|(_, m)| matches!(m, Membership::Outgroup))
|
||||||
|
.map(|(i, _)| i).collect();
|
||||||
|
Ok(Selection { ingroup_idx, outgroup_idx })
|
||||||
|
}
|
||||||
|
|
||||||
|
/// Build a `GroupQuorumFilter` from this selector's classification of `meta`.
|
||||||
///
|
///
|
||||||
/// - No groups defined: `ingroup_idx` = all genomes (implicit ingroup).
|
/// - No groups defined: `ingroup_idx` = all genomes (implicit ingroup).
|
||||||
/// - `ingroup` predicates only: outgroup indices are empty.
|
/// - `ingroup` predicates only: outgroup indices are empty.
|
||||||
@@ -173,34 +195,20 @@ impl IndexMeta {
|
|||||||
/// - Both defined: ingroup wins on overlap; uncategorized genomes are ignored.
|
/// - Both defined: ingroup wins on overlap; uncategorized genomes are ignored.
|
||||||
pub fn build_group_filter(
|
pub fn build_group_filter(
|
||||||
&self,
|
&self,
|
||||||
ingroup_preds: &[MetaPred],
|
meta: &IndexMeta,
|
||||||
outgroup_preds: &[MetaPred],
|
p: GroupFilterParams,
|
||||||
p: GroupFilterParams,
|
|
||||||
) -> Result<GroupQuorumFilter, String> {
|
) -> Result<GroupQuorumFilter, String> {
|
||||||
let genomes = self.genomes().map_err(|e| e.to_string())?;
|
let selection = self.run(meta)?;
|
||||||
let (ingroup_idx, outgroup_idx) = if ingroup_preds.is_empty() && outgroup_preds.is_empty() {
|
let in_size = selection.ingroup_idx.len();
|
||||||
((0..genomes.len()).collect(), vec![])
|
let out_size = selection.outgroup_idx.len();
|
||||||
} else {
|
|
||||||
let members = classify(&genomes, ingroup_preds, outgroup_preds);
|
|
||||||
let in_idx: Vec<usize> = members.iter().enumerate()
|
|
||||||
.filter(|(_, m)| matches!(m, Membership::Ingroup))
|
|
||||||
.map(|(i, _)| i).collect();
|
|
||||||
let out_idx: Vec<usize> = members.iter().enumerate()
|
|
||||||
.filter(|(_, m)| matches!(m, Membership::Outgroup))
|
|
||||||
.map(|(i, _)| i).collect();
|
|
||||||
(in_idx, out_idx)
|
|
||||||
};
|
|
||||||
|
|
||||||
let in_size = ingroup_idx.len();
|
|
||||||
let out_size = outgroup_idx.len();
|
|
||||||
|
|
||||||
let ingroup_quorum_explicit = p.min_count.is_some() || p.max_count.is_some()
|
let ingroup_quorum_explicit = p.min_count.is_some() || p.max_count.is_some()
|
||||||
|| p.min_frac.is_some() || p.max_frac.is_some();
|
|| p.min_frac.is_some() || p.max_frac.is_some();
|
||||||
let outgroup_quorum_explicit = p.min_outgroup_count.is_some() || p.max_outgroup_count.is_some()
|
let outgroup_quorum_explicit = p.min_outgroup_count.is_some() || p.max_outgroup_count.is_some()
|
||||||
|| p.min_outgroup_frac.is_some() || p.max_outgroup_frac.is_some();
|
|| p.min_outgroup_frac.is_some() || p.max_outgroup_frac.is_some();
|
||||||
|
|
||||||
let default_min_frac = if !ingroup_preds.is_empty() && !ingroup_quorum_explicit { 1.0 } else { 0.0 };
|
let default_min_frac = if !self.ingroup.is_empty() && !ingroup_quorum_explicit { 1.0 } else { 0.0 };
|
||||||
let default_max_outgroup_count = if !outgroup_preds.is_empty() && !outgroup_quorum_explicit { 0 } else { out_size };
|
let default_max_outgroup_count = if !self.outgroup.is_empty() && !outgroup_quorum_explicit { 0 } else { out_size };
|
||||||
|
|
||||||
// Resolve a signed count: negative means an offset from the group size
|
// Resolve a signed count: negative means an offset from the group size
|
||||||
// (e.g. -1 = all but one), floored at 1 so the negative form always keeps
|
// (e.g. -1 = all but one), floored at 1 so the negative form always keeps
|
||||||
@@ -238,8 +246,7 @@ impl IndexMeta {
|
|||||||
}
|
}
|
||||||
|
|
||||||
Ok(GroupQuorumFilter {
|
Ok(GroupQuorumFilter {
|
||||||
ingroup_idx,
|
selection,
|
||||||
outgroup_idx,
|
|
||||||
threshold: p.threshold,
|
threshold: p.threshold,
|
||||||
min_count,
|
min_count,
|
||||||
max_count,
|
max_count,
|
||||||
@@ -252,3 +259,4 @@ impl IndexMeta {
|
|||||||
})
|
})
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|||||||
Reference in New Issue
Block a user