Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
97 changes: 71 additions & 26 deletions crates/fgumi-umi/src/assigner.rs
Original file line number Diff line number Diff line change
Expand Up @@ -224,10 +224,7 @@ use std::sync::atomic::{AtomicU64, Ordering};
use crate::MoleculeId;
use fgumi_dna::BitEnc;

/// Default minimum number of unique UMIs to use N-gram/BK-tree index for candidate search.
/// Below this threshold, linear scan is used. Based on benchmarks: index provides ~10%
/// speedup even at 100 UMIs/position with negligible construction overhead.
pub const DEFAULT_INDEX_THRESHOLD: usize = 100;
use crate::index_threshold::{DEFAULT_INDEX_THRESHOLD, IndexThreshold};

/// A node in the BK-tree.
///
Expand Down Expand Up @@ -545,6 +542,26 @@ pub enum Strategy {
}

impl Strategy {
/// Whether this strategy consults the UMI index at all, at `edits`.
///
/// `Identity` and `Edit` have no index code path. `Adjacency` indexes only at
/// `edits == 1` — `build_adjacency_graph_bitenc` gates on it, and the BK-tree branch
/// that would serve larger distances is reachable only from the generic
/// `build_adjacency_graph` the paired strategy uses. `Paired` therefore indexes at
/// every edit distance.
///
/// Callers use this to reject an [`IndexThreshold`] that
/// [demands indexing](IndexThreshold::demands_indexing) under a configuration that
/// cannot deliver it.
#[must_use]
pub fn can_use_index(&self, edits: u32) -> bool {
match self {
Strategy::Identity | Strategy::Edit => false,
Strategy::Adjacency => edits == 1,
Strategy::Paired => true,
}
}

/// Create a new UMI assigner for this strategy
///
/// # Arguments
Expand All @@ -560,7 +577,7 @@ impl Strategy {
/// Panics if `edits` is non-zero when using the Identity strategy.
#[must_use]
pub fn new_assigner(&self, edits: u32) -> Box<dyn UmiAssigner> {
self.new_assigner_full(edits, 1, DEFAULT_INDEX_THRESHOLD)
self.new_assigner_full(edits, 1, IndexThreshold::default())
}

/// Create a new UMI assigner with thread count specified
Expand All @@ -575,7 +592,7 @@ impl Strategy {
/// A boxed trait object implementing the `UmiAssigner` trait
#[must_use]
pub fn new_assigner_with_threads(&self, edits: u32, threads: usize) -> Box<dyn UmiAssigner> {
self.new_assigner_full(edits, threads, DEFAULT_INDEX_THRESHOLD)
self.new_assigner_full(edits, threads, IndexThreshold::default())
}

/// Create a new UMI assigner with all parameters specified
Expand All @@ -596,12 +613,13 @@ impl Strategy {
///
/// Panics if `edits` is non-zero with the `Identity` strategy.
#[must_use]
pub fn new_assigner_full(
pub fn new_assigner_full<T: Into<IndexThreshold>>(
&self,
edits: u32,
threads: usize,
index_threshold: usize,
index_threshold: T,
) -> Box<dyn UmiAssigner> {
let index_threshold = index_threshold.into();
match self {
Strategy::Identity => {
assert_eq!(edits, 0, "Edits should be zero when using the identity UMI assigner");
Expand Down Expand Up @@ -1176,7 +1194,7 @@ struct Node {
pub struct AdjacencyUmiAssigner {
max_mismatches: u32,
threads: usize,
index_threshold: usize,
index_threshold: IndexThreshold,
counter: AtomicU64,
}

Expand All @@ -1194,8 +1212,17 @@ impl AdjacencyUmiAssigner {
///
/// A new `AdjacencyUmiAssigner` instance
#[must_use]
pub fn new(max_mismatches: u32, threads: usize, index_threshold: usize) -> Self {
Self { max_mismatches, threads, index_threshold, counter: AtomicU64::new(0) }
pub fn new<T: Into<IndexThreshold>>(
max_mismatches: u32,
threads: usize,
index_threshold: T,
) -> Self {
Self {
max_mismatches,
threads,
index_threshold: index_threshold.into(),
counter: AtomicU64::new(0),
}
}

/// Whether a position group of `num_umis` distinct UMIs is *eligible* to be
Expand All @@ -1216,7 +1243,7 @@ impl AdjacencyUmiAssigner {
/// and edit-distance conditions permit indexing", not "an index was built".
#[must_use]
pub fn uses_index(&self, num_umis: usize) -> bool {
num_umis >= self.index_threshold && self.max_mismatches == 1
self.index_threshold.admits(num_umis) && self.max_mismatches == 1
}

/// Generate the next unique molecule ID
Expand Down Expand Up @@ -1461,7 +1488,7 @@ impl AdjacencyUmiAssigner {
// Unlike `build_adjacency_graph_bitenc`, this path indexes at every edit
// distance: N-gram for k=1 (most efficient), BK-tree for k>1. Only the
// `index_threshold` minimum applies -- see `PairedUmiAssigner::uses_index`.
let (ngram_index, bktree) = if nodes.len() >= self.index_threshold {
let (ngram_index, bktree) = if self.index_threshold.admits(nodes.len()) {
let umis: Vec<(usize, &str)> =
umi_counts.iter().enumerate().map(|(i, (umi, _, _))| (i, umi.as_str())).collect();

Expand Down Expand Up @@ -1793,7 +1820,11 @@ impl PairedUmiAssigner {
///
/// A new `PairedUmiAssigner` instance
#[must_use]
pub fn new_with_threads(max_mismatches: u32, threads: usize, index_threshold: usize) -> Self {
pub fn new_with_threads<T: Into<IndexThreshold>>(
max_mismatches: u32,
threads: usize,
index_threshold: T,
) -> Self {
let prefix_len = max_mismatches as usize + 1;
Self {
adjacency: AdjacencyUmiAssigner::new(max_mismatches, threads, index_threshold),
Expand All @@ -1816,7 +1847,7 @@ impl PairedUmiAssigner {
/// then falls back to the linear scan.
#[must_use]
pub fn uses_index(&self, num_umis: usize) -> bool {
num_umis >= self.adjacency.index_threshold
self.adjacency.index_threshold.admits(num_umis)
}

/// Get the prefix for lower-ordered reads in a pair
Expand Down Expand Up @@ -3859,25 +3890,38 @@ mod tests {
// ========================================================================

/// Pins what `--index-threshold` actually does, so the CLI help can be checked
/// against it. The gate is `num_umis >= index_threshold`, so `0` means the index
/// is ALWAYS built -- it does not disable the index. Disabling it requires a
/// threshold larger than any position group, e.g. `usize::MAX`.
/// against it. A bare integer is a MINIMUM group size, not a switch, so `0` turns
/// the index fully on. `never` is the only way to turn it off.
#[rstest]
#[case::zero_always_indexes(0, 1, true)]
#[case::zero_indexes_empty_group(0, 0, true)]
#[case::below_threshold_scans(100, 99, false)]
#[case::at_threshold_indexes(100, 100, true)]
#[case::above_threshold_indexes(100, 101, true)]
#[case::usize_max_never_indexes(usize::MAX, 1_000_000, false)]
#[case::zero_always_indexes(IndexThreshold::MinUmis(0), 1, true)]
#[case::zero_indexes_empty_group(IndexThreshold::MinUmis(0), 0, true)]
#[case::below_threshold_scans(IndexThreshold::MinUmis(100), 99, false)]
#[case::at_threshold_indexes(IndexThreshold::MinUmis(100), 100, true)]
#[case::above_threshold_indexes(IndexThreshold::MinUmis(100), 101, true)]
#[case::always_indexes_singleton(IndexThreshold::Always, 1, true)]
#[case::never_skips_huge_group(IndexThreshold::Never, 1_000_000, false)]
fn test_adjacency_index_threshold_is_a_minimum_not_a_switch(
#[case] index_threshold: usize,
#[case] index_threshold: IndexThreshold,
#[case] num_umis: usize,
#[case] expected: bool,
) {
let assigner = AdjacencyUmiAssigner::new(1, 1, index_threshold);
assert_eq!(assigner.uses_index(num_umis), expected);
}

/// `never` wins over every other consideration, including an edit distance that
/// would otherwise index (paired) -- there is no combination that resurrects it.
#[rstest]
#[case::adjacency_edits_one(1)]
#[case::adjacency_edits_two(2)]
fn test_never_suppresses_the_index_at_every_edit_distance(#[case] max_mismatches: u32) {
let adjacency = AdjacencyUmiAssigner::new(max_mismatches, 1, IndexThreshold::Never);
assert!(!adjacency.uses_index(1_000_000));

let paired = PairedUmiAssigner::new_with_threads(max_mismatches, 1, IndexThreshold::Never);
assert!(!paired.uses_index(1_000_000));
}

/// The adjacency assigner only indexes at `--edits 1`; every other value falls
/// back to a linear scan over all UMI pairs regardless of `--index-threshold`.
/// `build_adjacency_graph_bitenc` gates on `self.max_mismatches == 1`, and the
Expand Down Expand Up @@ -3908,7 +3952,8 @@ mod tests {
#[case] num_umis: usize,
#[case] expected: bool,
) {
let assigner = PairedUmiAssigner::new_with_threads(max_mismatches, 1, 100);
let assigner =
PairedUmiAssigner::new_with_threads(max_mismatches, 1, IndexThreshold::MinUmis(100));
assert_eq!(assigner.uses_index(num_umis), expected);
}
}
178 changes: 178 additions & 0 deletions crates/fgumi-umi/src/index_threshold.rs
Original file line number Diff line number Diff line change
@@ -0,0 +1,178 @@
//! When the UMI-matching index is used instead of a linear scan.

use std::fmt;
use std::str::FromStr;

/// Default minimum number of unique UMIs to use N-gram/BK-tree index for candidate search.
/// Below this threshold, linear scan is used. Based on benchmarks: index provides ~10%
/// speedup even at 100 UMIs/position with negligible construction overhead.
pub const DEFAULT_INDEX_THRESHOLD: usize = 100;

/// Value of `--index-threshold`: when to match UMIs through the N-gram/BK-tree index
/// rather than by comparing every pair.
///
/// The index is a pure optimization — it narrows the candidate set, it does not change
/// which UMIs group together — so every variant yields the same grouping and differs only
/// in how much work reaching it takes.
///
/// [`MinUmis`](Self::MinUmis) is a *minimum group size*, not a switch, and that has caught
/// people out: `MinUmis(0)` admits every group, so a plain `0` turns the index fully **on**
/// — the opposite of what the CLI help claimed for a long time. [`Always`](Self::Always) and
/// [`Never`](Self::Never) let that intent be stated in a word rather than encoded in a
/// number no one can read back. `Never` is not otherwise reachable: suppressing the index
/// through `MinUmis` alone would take a value larger than any position group.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum IndexThreshold {
/// Index every position group, whatever its size.
Always,
/// Never index; always compare all UMI pairs.
Never,
/// Index a position group once it holds at least this many distinct UMIs.
MinUmis(usize),
}

impl IndexThreshold {
/// Whether a position group of `num_umis` distinct UMIs is matched through the index.
#[must_use]
pub fn admits(self, num_umis: usize) -> bool {
match self {
Self::Always => true,
Self::Never => false,
Self::MinUmis(min) => num_umis >= min,
}
}

/// Whether this value *asserts* that indexing must happen, rather than merely tuning
/// when it starts.
///
/// Only the `always` keyword asserts. A number — including `0`, which admits every
/// group exactly as `Always` does — is a tuning knob allowed to end up inert, which is
/// what lets the default `100` coexist with strategies that never index. Callers use
/// this to reject a request they cannot honour instead of silently ignoring it.
#[must_use]
pub fn demands_indexing(self) -> bool {
matches!(self, Self::Always)
}
}

impl Default for IndexThreshold {
fn default() -> Self {
Self::MinUmis(DEFAULT_INDEX_THRESHOLD)
}
}

impl From<usize> for IndexThreshold {
fn from(min_umis: usize) -> Self {
Self::MinUmis(min_umis)
}
}

impl FromStr for IndexThreshold {
type Err = String;

fn from_str(s: &str) -> std::result::Result<Self, Self::Err> {
if s.eq_ignore_ascii_case("always") {
Ok(Self::Always)
} else if s.eq_ignore_ascii_case("never") {
Ok(Self::Never)
} else {
s.parse::<usize>().map(Self::MinUmis).map_err(|e| {
format!("expected 'always', 'never', or a non-negative integer, got '{s}': {e}")
})
}
}
}

impl fmt::Display for IndexThreshold {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
match self {
Self::Always => f.write_str("always"),
Self::Never => f.write_str("never"),
Self::MinUmis(min) => write!(f, "{min}"),
}
}
}

#[cfg(test)]
mod tests {
use super::*;
use rstest::rstest;

/// `--index-threshold` accepts the two keywords and any non-negative integer.
#[rstest]
#[case::always("always", IndexThreshold::Always)]
#[case::always_mixed_case("Always", IndexThreshold::Always)]
#[case::never("never", IndexThreshold::Never)]
#[case::never_upper("NEVER", IndexThreshold::Never)]
#[case::zero("0", IndexThreshold::MinUmis(0))]
#[case::default_value("100", IndexThreshold::MinUmis(100))]
fn test_index_threshold_parses(#[case] input: &str, #[case] expected: IndexThreshold) {
assert_eq!(input.parse::<IndexThreshold>(), Ok(expected));
}

/// Rejected input names the accepted forms, so the error is actionable.
#[rstest]
#[case::negative("-1")]
#[case::word("off")]
#[case::empty("")]
#[case::float("1.5")]
fn test_index_threshold_rejects_invalid(#[case] input: &str) {
let err = input.parse::<IndexThreshold>().expect_err("must reject");
assert!(
err.contains("always") && err.contains("never"),
"error should name the accepted forms, got: {err}"
);
}

/// `admits` is the gate the assigners consult. `MinUmis` is a MINIMUM, so `0`
/// admits every group — the reading the CLI help had backwards for a long time.
#[rstest]
#[case::always_admits_empty(IndexThreshold::Always, 0, true)]
#[case::always_admits_huge(IndexThreshold::Always, 1_000_000, true)]
#[case::never_rejects_huge(IndexThreshold::Never, 1_000_000, false)]
#[case::zero_admits_singleton(IndexThreshold::MinUmis(0), 1, true)]
#[case::below_minimum(IndexThreshold::MinUmis(100), 99, false)]
#[case::at_minimum(IndexThreshold::MinUmis(100), 100, true)]
#[case::above_minimum(IndexThreshold::MinUmis(100), 101, true)]
fn test_index_threshold_admits(
#[case] threshold: IndexThreshold,
#[case] num_umis: usize,
#[case] expected: bool,
) {
assert_eq!(threshold.admits(num_umis), expected);
}

/// Only the `always` keyword asserts that indexing must happen. A number is a
/// tuning knob allowed to end up inert, which is what lets the default `100`
/// coexist with strategies that never index.
#[rstest]
#[case::always_demands(IndexThreshold::Always, true)]
#[case::never_does_not(IndexThreshold::Never, false)]
#[case::zero_does_not(IndexThreshold::MinUmis(0), false)]
#[case::default_does_not(IndexThreshold::MinUmis(100), false)]
fn test_only_always_demands_indexing(
#[case] threshold: IndexThreshold,
#[case] expected: bool,
) {
assert_eq!(threshold.demands_indexing(), expected);
}

/// Values round-trip through `Display`, so logs and help text stay parseable.
#[rstest]
#[case::always(IndexThreshold::Always, "always")]
#[case::never(IndexThreshold::Never, "never")]
#[case::min_umis(IndexThreshold::MinUmis(100), "100")]
fn test_index_threshold_displays_and_round_trips(
#[case] threshold: IndexThreshold,
#[case] expected: &str,
) {
assert_eq!(threshold.to_string(), expected);
assert_eq!(expected.parse::<IndexThreshold>(), Ok(threshold));
}

/// The default matches the documented `--index-threshold` default.
#[test]
fn test_default_is_the_documented_threshold() {
assert_eq!(IndexThreshold::default(), IndexThreshold::MinUmis(DEFAULT_INDEX_THRESHOLD));
}
}
2 changes: 2 additions & 0 deletions crates/fgumi-umi/src/lib.rs
Original file line number Diff line number Diff line change
Expand Up @@ -8,12 +8,14 @@
//! - UMI grouping and family analysis

pub mod assigner;
pub mod index_threshold;

// Re-export commonly used items
pub use assigner::{
AdjacencyUmiAssigner, IdentityUmiAssigner, PairedUmiAssigner, SimpleErrorUmiAssigner, Strategy,
Umi, UmiAssigner,
};
pub use index_threshold::{DEFAULT_INDEX_THRESHOLD, IndexThreshold};

use std::collections::HashSet;

Expand Down
Loading
Loading