103 simplify input of rq ranges

Cargo.toml

@@ -22,7 +22,7 @@ package.repository = "https://github.com/sebschmi/template-switch-aligner"
 
 [workspace.dependencies]
 serde = "1.0.219"
-compact-genome = "12.3.0"
+compact-genome = "12.5.0"
 traitsequence = "8.1.2"
 log = "0.4.27"
 num-traits = "0.2.19"

lib_tsalign/src/a_star_aligner.rs

@@ -154,7 +154,7 @@ pub fn template_switch_distance_a_star_align<
     query: &SubsequenceType,
     reference_name: &str,
     query_name: &str,
-    range: Option<AlignmentRange>,
+    range: AlignmentRange,
     config: &config::TemplateSwitchConfig<
         Strategies::Alphabet,
         <Strategies as AlignmentStrategySelector>::Cost,
@@ -167,6 +167,9 @@ pub fn template_switch_distance_a_star_align<
 where
     Strategies::Cost: From<u64>,
 {
+    debug!("Reference sequence: {}", reference.as_string());
+    debug!("Query sequence: {}", query.as_string());
+
     let memory = Memory {
         template_switch_min_length:
             <Strategies::TemplateSwitchMinLength as TemplateSwitchMinLengthStrategy<
@@ -180,6 +183,7 @@ where
         primary_match: (),
     };
 
+    info!("Calling aligner...");
     let mut result = a_star_align(template_switch_distance::Context::<
         SubsequenceType,
         Strategies,
@@ -195,12 +199,18 @@ where
         memory_limit,
         force_label_correcting,
     ));
-    debug!("CIGAR before extending: {}", result.cigar());
+    info!("Main alignment finished");
 
     info!("Extending template switches");
+    debug!("CIGAR before extending: {}", result.cigar());
     let mut range = range;
-    result.extend_beyond_range_with_equal_cost(reference, query, &mut range, config);
+    let extension_steps =
+        result.extend_beyond_range_without_increasing_cost(reference, query, &mut range, config);
     let range = range;
+    info!(
+        "Extended alignment {extension_steps} steps beyond the alignment range without increasing alignment costs"
+    );
+    info!("Alignment ranges after extension {range}");
 
     result.compute_ts_equal_cost_ranges(reference, query, &range, config);
     result

lib_tsalign/src/a_star_aligner/alignment_result.rs

@@ -237,32 +237,31 @@ impl<AlignmentType: IAlignmentType, Cost: AStarCost> AlignmentResult<AlignmentTy
 impl<Cost: AStarCost + From<u64>>
     AlignmentResult<super::template_switch_distance::AlignmentType, Cost>
 {
+    /// Heuristically extend the alignment beyond the alignment range as long as this does not increase the alignment cost.
+    ///
+    /// Returns the amount of steps the alignment was extended.
     #[allow(clippy::too_many_lines)]
-    pub fn extend_beyond_range_with_equal_cost<
+    pub fn extend_beyond_range_without_increasing_cost<
         AlphabetType: Alphabet,
         SubsequenceType: GenomeSequence<AlphabetType, SubsequenceType> + ?Sized,
     >(
         &mut self,
         reference: &SubsequenceType,
         query: &SubsequenceType,
-        range: &mut Option<AlignmentRange>,
+        range: &mut AlignmentRange,
         config: &TemplateSwitchConfig<AlphabetType, Cost>,
-    ) {
-        let Some(range) = range else {
-            trace!("No range given to extend beyond.");
-            return;
-        };
+    ) -> usize {
         let Self::WithTarget {
             alignment,
             statistics,
         } = self
         else {
             trace!("There is no alignment, therefore we cannot postprocess it.");
-            return;
+            return 0;
         };
         if config.left_flank_length > 0 || config.right_flank_length > 0 {
             warn!("Alignment extension does not support flanks");
-            return;
+            return 0;
         }
 
         // Compute cost before extending.
@@ -280,6 +279,7 @@ impl<Cost: AStarCost + From<u64>>
             "computed cost {current_cost} != alignment cost {alignment_cost}; {}",
             alignment.cigar()
         );
+        let mut extension_steps = 0;
 
         // Extend left with equal cost.
         while let Some(new_range) = range.move_offsets_left() {
@@ -312,8 +312,10 @@ impl<Cost: AStarCost + From<u64>>
                 config,
             );
 
+            trace!("Extending left to {new_range}. Cost {current_cost} -> {new_cost}");
+
             if new_cost > current_cost {
-                // If cost is not equal, revert alignment and break.
+                // If cost is increasing, revert alignment and break.
                 alignment.inner_mut()[0].0 -= 1;
                 if alignment.inner_mut()[0].0 == 0 {
                     alignment.inner_mut().remove(0);
@@ -322,6 +324,7 @@ impl<Cost: AStarCost + From<u64>>
             }
             current_cost = new_cost;
             *range = new_range;
+            extension_steps += 1;
         }
 
         // Extend right with equal cost.
@@ -364,23 +367,28 @@ impl<Cost: AStarCost + From<u64>>
                 config,
             );
 
+            trace!("Extending right to {new_range}. Cost {current_cost} -> {new_cost}");
+
             if new_cost > current_cost {
-                // If cost is not equal, revert alignment and break.
+                // If cost is increasing, revert alignment and break.
                 let Some((n, _)) = alignment.inner_mut().last_mut() else {
                     unreachable!("As we have just added at least one element")
                 };
                 *n -= 1;
                 if *n == 0 {
                     alignment.inner_mut().pop();
                 }
+                break;
             }
             current_cost = new_cost;
             *range = new_range;
+            extension_steps += 1;
         }
 
         // Update statistics from updated range
         statistics.reference_offset = range.reference_offset();
         statistics.query_offset = range.query_offset();
+        extension_steps
     }
 
     #[allow(clippy::too_many_lines)]
@@ -391,7 +399,7 @@ impl<Cost: AStarCost + From<u64>>
         &mut self,
         reference: &SubsequenceType,
         query: &SubsequenceType,
-        range: &Option<AlignmentRange>,
+        range: &AlignmentRange,
         config: &TemplateSwitchConfig<AlphabetType, Cost>,
     ) {
         let Self::WithTarget {
@@ -407,13 +415,8 @@ impl<Cost: AStarCost + From<u64>>
             return;
         }
 
-        let reference_offset = range.as_ref().map_or(
-            0,
-            super::alignment_geometry::AlignmentRange::reference_offset,
-        );
-        let query_offset = range
-            .as_ref()
-            .map_or(0, super::alignment_geometry::AlignmentRange::query_offset);
+        let reference_offset = range.reference_offset();
+        let query_offset = range.query_offset();
 
         for i in 0..alignment.inner_mut().len() {
             let alignment_type = alignment.inner_mut()[i].1;

lib_tsalign/src/a_star_aligner/configurable_a_star_align.rs

@@ -5,6 +5,7 @@ use compact_genome::{
     interface::sequence::{GenomeSequence, OwnedGenomeSequence},
 };
 use generic_a_star::cost::U64Cost;
+use traitsequence::interface::Sequence;
 
 use crate::{
     a_star_aligner::{
@@ -281,7 +282,8 @@ impl<AlphabetType: Alphabet> Aligner<AlphabetType> {
             query.as_genome_subsequence(),
             data.reference_name,
             data.query_name,
-            data.ranges,
+            data.ranges
+                .unwrap_or_else(|| AlignmentRange::new_complete(reference.len(), query.len())),
             &self.costs,
             cost_limit,
             data.memory_limit,

lib_tsalign/src/a_star_aligner/template_switch_distance/context.rs

@@ -80,7 +80,7 @@ impl<
         query: &'query SubsequenceType,
         reference_name: &str,
         query_name: &str,
-        range: Option<AlignmentRange>,
+        range: AlignmentRange,
         config: TemplateSwitchConfig<Strategies::Alphabet, Strategies::Cost>,
         memory: Memory<Strategies>,
         cost_limit: Option<Strategies::Cost>,
@@ -92,8 +92,7 @@ impl<
             query,
             reference_name: reference_name.to_owned(),
             query_name: query_name.to_owned(),
-            range: range
-                .unwrap_or_else(|| AlignmentRange::new_complete(reference.len(), query.len())),
+            range,
             config,
             a_star_buffers: Default::default(),
             memory,

lib_tsalign/src/a_star_aligner/template_switch_distance/lower_bounds/template_switch.rs

@@ -10,9 +10,11 @@ use compact_genome::{
 };
 use generic_a_star::{AStar, AStarNode, AStarResult, cost::AStarCost};
 use log::{debug, info, trace};
+use traitsequence::interface::Sequence;
 
 use crate::{
     a_star_aligner::{
+        alignment_geometry::AlignmentRange,
         alignment_result::IAlignmentType,
         template_switch_distance::{
             Context, Identifier, Node,
@@ -93,7 +95,7 @@ impl<Cost: AStarCost> TemplateSwitchLowerBoundMatrix<Cost> {
                     genome.as_genome_subsequence(),
                     "",
                     "",
-                    None,
+                    AlignmentRange::new_complete(genome.len(), genome.len()),
                     lower_bound_config.clone(),
                     Memory {
                         template_switch_min_length: (),

lib_tsalign/src/a_star_aligner/template_switch_distance/lower_bounds/template_switch_alignment.rs

@@ -11,9 +11,11 @@ use compact_genome::{
 use generic_a_star::{AStar, AStarNode, AStarResult, cost::AStarCost};
 use log::{debug, info, trace};
 use ndarray::Array2;
+use traitsequence::interface::Sequence;
 
 use crate::{
     a_star_aligner::{
+        alignment_geometry::AlignmentRange,
         alignment_result::IAlignmentType,
         template_switch_distance::{
             Context, Identifier,
@@ -95,7 +97,7 @@ impl<Cost: AStarCost> TemplateSwitchAlignmentLowerBoundMatrix<Cost> {
                 genome.as_genome_subsequence(),
                 "",
                 "",
-                None,
+                AlignmentRange::new_complete(genome.len(), genome.len()),
                 lower_bound_config.clone(),
                 Memory {
                     template_switch_min_length: (),

test_files/twin_embedded.fa

@@ -0,0 +1,4 @@
+>REF
+ACCA|GGAA|TTGG
+>QRY
+AGCA|GGAA|TAGG

tsalign-tests/tests/integration.rs

@@ -22,3 +22,8 @@ fn test_align_with_cost_limit() -> Result<()> {
 fn test_align_with_memory_limit() -> Result<()> {
     run_in_repo_root("align -p test_files/twin_100_0.01.fa --memory-limit 1000")
 }
+
+#[test]
+fn test_align_with_embedded_rq_ranges() -> Result<()> {
+    run_in_repo_root("align -p test_files/twin_embedded.fa --use-embedded-rq-ranges")
+}

tsalign/Cargo.toml

@@ -21,3 +21,4 @@ toml = "0.8.19"
 log.workspace = true
 simplelog = "0.12.2"
 anyhow = "1.0.97"
+utf8-chars = "3.0.5"