Merge pull request #194 from algbio/193-failing-debug-assertion-in-compute_ts_equal_cost_ranges

193 failing debug assertion in compute ts equal cost ranges

Author: sebschmi

lib_tsalign/src/a_star_aligner/alignment_result/alignment/template_switch_specifics.rs

@@ -759,12 +759,11 @@ impl Alignment<AlignmentType> {
                     descendant = ts_descendant;
                     ancestor = ts_ancestor;
                     direction = ts_direction;
-                    let cost_increment = config.base_cost.get(descendant, ancestor, direction);
-                    let Some(cost_increment) = cost_increment.checked_add(
-                        &config
-                            .offset_costs(descendant, ancestor)
-                            .evaluate(&first_offset),
-                    ) else {
+                    let base_cost = config.base_cost.get(descendant, ancestor, direction);
+                    let offset_cost = config
+                        .offset_costs(descendant, ancestor)
+                        .evaluate(&first_offset);
+                    let Some(cost_increment) = base_cost.checked_add(&offset_cost) else {
                         return Cost::max_value();
                     };
                     descendant_index = match descendant {
@@ -1298,7 +1297,6 @@ mod tests {
         for (expected_alignment, expected_cost) in
             START_ALIGNMENTS[1..].iter().zip(&START_COSTS[1..])
         {
-            println!("{}", alignment.cigar());
             assert!(alignment.move_template_switch_start_backwards(
                 reference.as_genome_subsequence(),
                 query.as_genome_subsequence(),

lib_tsalign/src/a_star_aligner/configurable_a_star_align.rs

@@ -62,6 +62,9 @@ use super::{
 pub use compact_genome::implementation::alphabets;
 pub use compact_genome::interface::alphabet::Alphabet;
 
+#[cfg(test)]
+mod tests;
+
 #[derive(Debug, Default)]
 #[cfg_attr(feature = "serde", derive(serde::Deserialize, serde::Serialize))]
 #[cfg_attr(feature = "serde", serde(rename_all = "snake_case"))]
@@ -434,34 +437,3 @@ impl<AlphabetType: Alphabet> Aligner<AlphabetType> {
         )
     }
 }
-
-#[cfg(test)]
-mod tests {
-    use num_traits::real::Real;
-
-    use crate::a_star_aligner::alignment_geometry::AlignmentCoordinates;
-
-    use super::*;
-
-    /// This test case uses real data and used to panic.
-    #[test]
-    fn test_panic() {
-        let mut aligner = Aligner::new();
-        aligner.set_min_length_strategy(MinLengthStrategySelector::PreprocessFilter);
-        let res = aligner.align(
-            "reference",
-            b"TACCGAGACTGCAGAAAGTGAAAGCTATACTAA",
-            "query",
-            b"TAACTTTTAATGCCAAATATTTTATCCAAATAGGAAATTGTTTTCCGGTAAAATTTAACAAAAGAACCAGTTTACCCCCTTCAATGATTTATTTTTCTTCTTAGATTGAACTCTCGGGTTAGATCTCATTTTAACTGAAATTTGGTAAAAAATCCATATTACGGTTCAAGCCTAACCGAGACTGCAGAAAGTGAAAGCTAAAAGCTAATTTTTTTTTTTTTTTTGTATTTCACACCTATCGCAATACATCCTGGACAACACTGTATATTGAAACATTTTTTGCCTACAGCAATGGGCCTATAATTTTTTCTCGGCATTAGCTCTACAATCCAATTCTATCCTGCTTCTTCTTGTAAACAGGGATAACTTTAACTAACATTCAGTTTGCTTGGGAAAGAACCGATTGATAATGTA",
-            AlignmentRange::new_offset_limit(
-                AlignmentCoordinates::new(27, 200),
-                AlignmentCoordinates::new(33, 214)
-            ).into(),
-            &[],
-            None,
-            None,
-        true);
-        println!("{res:#?}");
-        assert!(res.statistics().cost.is_sign_positive());
-    }
-}

lib_tsalign/src/a_star_aligner/configurable_a_star_align/tests.rs

@@ -0,0 +1,41 @@
+use num_traits::real::Real;
+
+use crate::a_star_aligner::alignment_geometry::AlignmentCoordinates;
+
+use super::*;
+
+/// This test case uses real data and used to panic.
+#[test]
+fn test_panic() {
+    let mut aligner = Aligner::new();
+    aligner.set_min_length_strategy(MinLengthStrategySelector::PreprocessFilter);
+    let res = aligner.align(
+            "reference",
+            b"TACCGAGACTGCAGAAAGTGAAAGCTATACTAA",
+            "query",
+            b"TAACTTTTAATGCCAAATATTTTATCCAAATAGGAAATTGTTTTCCGGTAAAATTTAACAAAAGAACCAGTTTACCCCCTTCAATGATTTATTTTTCTTCTTAGATTGAACTCTCGGGTTAGATCTCATTTTAACTGAAATTTGGTAAAAAATCCATATTACGGTTCAAGCCTAACCGAGACTGCAGAAAGTGAAAGCTAAAAGCTAATTTTTTTTTTTTTTTTGTATTTCACACCTATCGCAATACATCCTGGACAACACTGTATATTGAAACATTTTTTGCCTACAGCAATGGGCCTATAATTTTTTCTCGGCATTAGCTCTACAATCCAATTCTATCCTGCTTCTTCTTGTAAACAGGGATAACTTTAACTAACATTCAGTTTGCTTGGGAAAGAACCGATTGATAATGTA",
+            AlignmentRange::new_offset_limit(
+                AlignmentCoordinates::new(27, 200),
+                AlignmentCoordinates::new(33, 214)
+            ).into(),
+            &[],
+            None,
+            None,
+        true);
+    println!("{res:#?}");
+    assert!(res.statistics().cost.is_sign_positive());
+}
+
+#[test]
+fn test_193() {
+    //        [-- 10 --][-- 11 ---]1
+    let r = b"AAAAAAAAAACCCCCCCCCCCC";
+    let q = b"AAAAAAAAAATTTTTTTTTTTC";
+
+    let mut tsa = Aligner::<DnaAlphabetOrN>::default();
+
+    let mut costs = TemplateSwitchConfig::default();
+    costs.base_cost.qrr = 4u64.into();
+    tsa.set_costs(costs.clone());
+    let _tsa_res = tsa.align("ref", r, "qry", q, None, &[], None, None, false);
+}

lib_tsalign/src/a_star_aligner/template_switch_distance.rs

@@ -2,7 +2,7 @@ use std::fmt::Display;
 
 use compact_genome::interface::sequence::GenomeSequence;
 use generic_a_star::{AStarNode, cost::AStarCost};
-use num_traits::{Bounded, Zero};
+use num_traits::{Bounded, CheckedAdd, Zero};
 use strategies::{
     AlignmentStrategiesNodeMemory, AlignmentStrategySelector, node_ord::NodeOrdStrategy,
     primary_match::PrimaryMatchStrategy,
@@ -25,7 +25,7 @@ pub use identifier::{
 
 use crate::{
     a_star_aligner::template_switch_distance::strategies::descendant::TemplateSwitchDescendantStrategy,
-    config::BaseCost,
+    config::BaseCost, costs::cost_function::CostFunction,
 };
 
 #[derive(Debug, Clone, PartialEq, Eq)]
@@ -223,8 +223,8 @@ impl<Strategies: AlignmentStrategySelector> Node<Strategies> {
         SubsequenceType: GenomeSequence<Strategies::Alphabet, SubsequenceType> + ?Sized,
     >(
         &'result self,
-        rq_qr_cost_increment: Strategies::Cost,
-        rr_qq_cost_increment: Strategies::Cost,
+        rq_qr_offset_costs: &'result CostFunction<isize, Strategies::Cost>,
+        rr_qq_offset_costs: &'result CostFunction<isize, Strategies::Cost>,
         base_cost: &'result BaseCost<Strategies::Cost>,
         context: &'result Context<SubsequenceType, Strategies>,
     ) -> impl 'result + Iterator<Item = Self> {
@@ -234,10 +234,6 @@ impl<Strategies: AlignmentStrategySelector> Node<Strategies> {
         ) {
             unreachable!("This method is only called on primary nodes.")
         }
-        debug_assert!(
-            rq_qr_cost_increment < Strategies::Cost::max_value()
-                || rr_qq_cost_increment < Strategies::Cost::max_value(),
-        );
 
         self.node_data
             .identifier
@@ -268,34 +264,42 @@ impl<Strategies: AlignmentStrategySelector> Node<Strategies> {
                     *template_switch_ancestor,
                     *template_switch_direction,
                 );
-                let cost_increment = match (*template_switch_descendant, *template_switch_ancestor)
-                {
+                let cost_function = match (*template_switch_descendant, *template_switch_ancestor) {
                     (TemplateSwitchDescendant::Reference, TemplateSwitchAncestor::Query)
                     | (TemplateSwitchDescendant::Query, TemplateSwitchAncestor::Reference) => {
-                        rq_qr_cost_increment
+                        rq_qr_offset_costs
                     }
                     (TemplateSwitchDescendant::Reference, TemplateSwitchAncestor::Reference)
                     | (TemplateSwitchDescendant::Query, TemplateSwitchAncestor::Query) => {
-                        rr_qq_cost_increment
+                        rr_qq_offset_costs
                     }
                 };
+                let cost_increment = cost_function.evaluate(template_switch_first_offset);
 
-                (base_cost != Strategies::Cost::max_value()
-                    && cost_increment != Strategies::Cost::max_value())
-                .then(|| {
-                    self.generate_successor(
-                        identifier,
-                        base_cost + cost_increment,
-                        AlignmentType::TemplateSwitchEntrance {
-                            descendant: *template_switch_descendant,
-                            ancestor: *template_switch_ancestor,
-                            direction: *template_switch_direction,
-                            equal_cost_range: EqualCostRange::new_invalid(),
-                            first_offset: *template_switch_first_offset,
-                        },
-                        context,
-                    )
-                })
+                if base_cost != Strategies::Cost::max_value()
+                    && cost_increment != Strategies::Cost::max_value()
+                {
+                    if let Some(cost_increment) = base_cost.checked_add(&cost_increment)
+                        && cost_increment != Strategies::Cost::max_value()
+                    {
+                        Some(self.generate_successor(
+                            identifier,
+                            cost_increment,
+                            AlignmentType::TemplateSwitchEntrance {
+                                descendant: *template_switch_descendant,
+                                ancestor: *template_switch_ancestor,
+                                direction: *template_switch_direction,
+                                equal_cost_range: EqualCostRange::new_invalid(),
+                                first_offset: *template_switch_first_offset,
+                            },
+                            context,
+                        ))
+                    } else {
+                        None
+                    }
+                } else {
+                    None
+                }
             })
     }
 

lib_tsalign/src/a_star_aligner/template_switch_distance/context.rs

@@ -208,8 +208,8 @@ impl<
                 opened_nodes_output.extend(
                     primary_node
                         .generate_initial_template_switch_entrance_successors(
-                            config.rq_qr_offset_costs.evaluate(&0),
-                            config.rr_qq_offset_costs.evaluate(&0),
+                            &config.rq_qr_offset_costs,
+                            &config.rr_qq_offset_costs,
                             &config.base_cost,
                             self,
                         )
@@ -444,22 +444,15 @@ impl<
 
                 // Template switches are always allowed, as long as we have a left flank.
                 if flank_index == config.left_flank_length && can_start_another_template_switch {
-                    let rq_qr_offset_costs = config.rq_qr_offset_costs.evaluate(&0);
-                    let rr_qq_offset_costs = config.rr_qq_offset_costs.evaluate(&0);
-
-                    if rq_qr_offset_costs != Strategies::Cost::max_value()
-                        || rr_qq_offset_costs != Strategies::Cost::max_value()
-                    {
-                        opened_nodes_output.extend(
-                            node.generate_initial_template_switch_entrance_successors(
-                                rq_qr_offset_costs,
-                                rr_qq_offset_costs,
-                                &config.base_cost,
-                                self,
-                            )
-                            .map(Into::into),
-                        );
-                    }
+                    opened_nodes_output.extend(
+                        node.generate_initial_template_switch_entrance_successors(
+                            &config.rq_qr_offset_costs,
+                            &config.rr_qq_offset_costs,
+                            &config.base_cost,
+                            self,
+                        )
+                        .map(Into::into),
+                    );
                 }
             }
 

lib_tsalign/src/config.rs

@@ -223,10 +223,10 @@ impl<AlphabetType: Alphabet, Cost: AStarCost> Default for TemplateSwitchConfig<A
             right_flank_length: 0,
             template_switch_min_length: 5,
             base_cost: BaseCost {
-                rrf: 4.into(),
-                rqf: 4.into(),
-                qrf: 4.into(),
-                qqf: 4.into(),
+                rrf: 40.into(),
+                rqf: 40.into(),
+                qrf: 40.into(),
+                qqf: 40.into(),
                 rrr: 3.into(),
                 rqr: 2.into(),
                 qrr: 2.into(),