From ea1446dfae21dd22ce9bb13df341075cd3295629 Mon Sep 17 00:00:00 2001 From: Kristoffer Sahlin Date: Wed, 30 Sep 2026 23:04:46 +0200 Subject: [PATCH 1/2] Break out of the chaining lookback loop earlier collinear_chaining looks back at up to max_lookback previous anchors. It stopped at an anchor on a different contig or, for collinear anchors only, at an anchor at least max_ref_gap away. Because anchors are sorted by reference start, all remaining anchors are then even further away and cannot be chained. Both checks are now combined into a single comparison of the reference start that is done first, for every anchor, so the loop also stops early when it encounters non-collinear anchors. The same bound is used when checking the best anchor so far. --- src/chainer.rs | 33 ++++++++++++++++++++------------- 1 file changed, 20 insertions(+), 13 deletions(-) diff --git a/src/chainer.rs b/src/chainer.rs index 6e02a78b..cb721275 100644 --- a/src/chainer.rs +++ b/src/chainer.rs @@ -145,12 +145,17 @@ impl Chainer { contig_index += 1; } let ref_contig_start = starts[contig_index]; + // Smallest ref_start of an anchor that is on the same contig as a[i] + // and less than max_ref_gap away from it + let min_ref_start = + ref_contig_start.max((ai.ref_start() + 1).saturating_sub(max_ref_gap)); for j in (lookup_end..i).rev() { let aj = anchors[j]; - // Do not attempt to chain to an anchor on a different contig. - // This check works because we go through anchors in reverse. - if aj.ref_start() < ref_contig_start { + // Do not attempt to chain to an anchor on a different contig or + // too far away. This check works because we go through anchors + // in reverse, so all remaining anchors are even further away. + if aj.ref_start() < min_ref_start { break; } @@ -166,10 +171,6 @@ impl Chainer { let dr = ai.ref_start() - aj.ref_start(); - if dr >= max_ref_gap { - break; - } - let diagonal_ratio = dq.max(dr) as f32 / dq.min(dr) as f32; if diagonal_ratio > self.parameters.max_diagonal_ratio { continue; @@ -194,15 +195,12 @@ impl Chainer { if best_index != usize::MAX { let aj = anchors[best_index]; - if aj.ref_start() >= ref_contig_start && ai.query_start() > aj.query_start() { + if aj.ref_start() >= min_ref_start && ai.query_start() > aj.query_start() { let dq = ai.query_start() - aj.query_start(); let dr = ai.ref_start() - aj.ref_start(); let diagonal_ratio = dq.max(dr) as f32 / dq.min(dr) as f32; - if dr < max_ref_gap - && dr > 0 - && diagonal_ratio <= self.parameters.max_diagonal_ratio - { + if dr > 0 && diagonal_ratio <= self.parameters.max_diagonal_ratio { let score = self.compute_score_cached(dq, dr); let new_score = dp[best_index] + score; if new_score > dp[i] { @@ -555,13 +553,22 @@ pub mod test { (95, 35), ]; let starts = ContigStarts::new(vec![0], 200); - let chaining_result = chainer.collinear_chaining(anchors, &starts, 2000); + let chaining_result = chainer.collinear_chaining(anchors.clone(), &starts, 2000); // The best chain has score 42.342842 and uses anchors 0, 1, 3. // When using the heuristic that breaks early if the predecessor is on the // same diagonal, a suboptimal chain is found that has score 38.2 and // consists of anchors 2 and 3. assert!(chaining_result.best_score > 42.0); + + // Anchors 1 and 3 are 65 apart on the reference, just within max_ref_gap + let parameters = ChainingParameters { + max_ref_gap: Some(66), + ..ChainingParameters::default() + }; + let chainer = Chainer::new(20, parameters); + let chaining_result = chainer.collinear_chaining(anchors, &starts, 2000); + assert!(chaining_result.best_score > 42.0); } #[test] From 927038706aefb7f00494c0821d11009a55dbdb4b Mon Sep 17 00:00:00 2001 From: Kristoffer Sahlin Date: Wed, 30 Sep 2026 23:04:46 +0200 Subject: [PATCH 2/2] Let the diagonal ratio test skip a division sometimes The division can be skipped when both distances are equal, in which case the ratio is 1. --- CHANGES.md | 4 ++++ src/chainer.rs | 14 +++++++++++--- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/CHANGES.md b/CHANGES.md index 0cc1db73..e9b3b85f 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,5 +1,9 @@ # Strobealign Changelog +## development version + +* #646: Speed up chaining slightly by avoiding some computations. + ## v0.19.0 (2026-10-08) This release improves mapping speed, in particular when mapping paired-end diff --git a/src/chainer.rs b/src/chainer.rs index cb721275..f197c1e7 100644 --- a/src/chainer.rs +++ b/src/chainer.rs @@ -171,9 +171,17 @@ impl Chainer { let dr = ai.ref_start() - aj.ref_start(); - let diagonal_ratio = dq.max(dr) as f32 / dq.min(dr) as f32; - if diagonal_ratio > self.parameters.max_diagonal_ratio { - continue; + if dq == dr { + // Same diagonal (the most common case): Avoid the division since the + // ratio is 1 (or NaN if dq == dr == 0, which never exceeds the maximum) + if dq > 0 && self.parameters.max_diagonal_ratio < 1.0 { + continue; + } + } else { + let diagonal_ratio = dq.max(dr) as f32 / dq.min(dr) as f32; + if diagonal_ratio > self.parameters.max_diagonal_ratio { + continue; + } } let score = self.compute_score_cached(dq, dr);