use super::{Alignment, AlignmentType, Scoring};
use crate::graph::{EdgeId, Graph};
#[inline]
fn backtrack_match_step(
graph: &Graph,
node_id_to_rank: &[u32],
sequence_profile: &[i32],
h: &[i32],
matrix_width: usize,
i: usize,
j: usize,
) -> Option<(usize, usize)> {
if i == 0 || j == 0 {
return None;
}
let h_ij = h[i * matrix_width + j];
let node = &graph.nodes[graph.rank_to_node[i - 1].0 as usize];
let match_cost = sequence_profile[node.code as usize * matrix_width + j];
let pred_row = |edge_id: EdgeId| -> usize {
let tail = graph.edges[edge_id.0 as usize].tail;
node_id_to_rank[tail.0 as usize] as usize + 1
};
let pred_first = if node.inedges.is_empty() {
0
} else {
pred_row(node.inedges[0])
};
if h_ij == h[pred_first * matrix_width + (j - 1)] + match_cost {
return Some((pred_first, j - 1));
}
for p in 1..node.inedges.len() {
let pred_i = pred_row(node.inedges[p]);
if h_ij == h[pred_i * matrix_width + (j - 1)] + match_cost {
return Some((pred_i, j - 1));
}
}
None
}
#[allow(clippy::too_many_arguments)]
#[inline]
pub(crate) fn backtrack_linear(
graph: &Graph,
node_id_to_rank: &[u32],
sequence_profile: &[i32],
h: &[i32],
matrix_width: usize,
alignment_type: AlignmentType,
scoring: &Scoring,
max_i: usize,
max_j: usize,
max_score: i32,
) -> Alignment {
if max_i == 0 && max_j == 0 {
return Vec::new();
}
debug_assert_eq!(
h[max_i * matrix_width + max_j],
max_score,
"fill's max_score must match H at (max_i, max_j)"
);
let g = i32::from(scoring.g);
let pred_row = |edge_id: EdgeId| -> usize {
let tail = graph.edges[edge_id.0 as usize].tail;
node_id_to_rank[tail.0 as usize] as usize + 1
};
let mut alignment: Alignment = Vec::new();
let mut i = max_i;
let mut j = max_j;
loop {
let keep_going = match alignment_type {
AlignmentType::Local => h[i * matrix_width + j] != 0,
AlignmentType::Global => !(i == 0 && j == 0),
AlignmentType::Overlap => !(i == 0 || j == 0),
};
if !keep_going {
break;
}
let h_ij = h[i * matrix_width + j];
let mut prev_i = 0usize;
let mut prev_j = 0usize;
let mut predecessor_found = false;
if let Some((pi, pj)) = backtrack_match_step(
graph,
node_id_to_rank,
sequence_profile,
h,
matrix_width,
i,
j,
) {
prev_i = pi;
prev_j = pj;
predecessor_found = true;
}
if !predecessor_found && i != 0 {
let node = &graph.nodes[graph.rank_to_node[i - 1].0 as usize];
let pred_first = if node.inedges.is_empty() {
0
} else {
pred_row(node.inedges[0])
};
if h_ij == h[pred_first * matrix_width + j] + g {
prev_i = pred_first;
prev_j = j;
predecessor_found = true;
} else {
for p in 1..node.inedges.len() {
let pred_i = pred_row(node.inedges[p]);
if h_ij == h[pred_i * matrix_width + j] + g {
prev_i = pred_i;
prev_j = j;
predecessor_found = true;
break;
}
}
}
}
if !predecessor_found && j != 0 && h_ij == h[i * matrix_width + (j - 1)] + g {
prev_i = i;
prev_j = j - 1;
}
let node_slot = if i == prev_i {
-1
} else {
graph.rank_to_node[i - 1].0 as i32
};
let seq_slot = if j == prev_j { -1 } else { (j - 1) as i32 };
alignment.push((node_slot, seq_slot));
i = prev_i;
j = prev_j;
}
alignment.reverse();
alignment
}
#[allow(clippy::too_many_arguments)]
#[inline]
pub(crate) fn backtrack_affine(
graph: &Graph,
node_id_to_rank: &[u32],
sequence_profile: &[i32],
h: &[i32],
e: &[i32],
f: &[i32],
matrix_width: usize,
alignment_type: AlignmentType,
scoring: &Scoring,
max_i: usize,
max_j: usize,
max_score: i32,
) -> Alignment {
if max_i == 0 && max_j == 0 {
return Vec::new();
}
debug_assert_eq!(
h[max_i * matrix_width + max_j],
max_score,
"fill's max_score must match H at (max_i, max_j)"
);
let g = i32::from(scoring.g);
let e_penalty = i32::from(scoring.e);
let pred_row = |edge_id: EdgeId| -> usize {
let tail = graph.edges[edge_id.0 as usize].tail;
node_id_to_rank[tail.0 as usize] as usize + 1
};
let mut alignment: Alignment = Vec::new();
let mut i = max_i;
let mut j = max_j;
loop {
let keep_going = match alignment_type {
AlignmentType::Local => h[i * matrix_width + j] != 0,
AlignmentType::Global => !(i == 0 && j == 0),
AlignmentType::Overlap => !(i == 0 || j == 0),
};
if !keep_going {
break;
}
let h_ij = h[i * matrix_width + j];
let mut prev_i = 0usize;
let mut prev_j = 0usize;
let mut predecessor_found = false;
let mut extend_left = false;
let mut extend_up = false;
if let Some((pi, pj)) = backtrack_match_step(
graph,
node_id_to_rank,
sequence_profile,
h,
matrix_width,
i,
j,
) {
prev_i = pi;
prev_j = pj;
predecessor_found = true;
}
if !predecessor_found && i != 0 {
let node = &graph.nodes[graph.rank_to_node[i - 1].0 as usize];
let pred_first = if node.inedges.is_empty() {
0
} else {
pred_row(node.inedges[0])
};
let a = h_ij == f[pred_first * matrix_width + j] + e_penalty;
extend_up = a;
if a || h_ij == h[pred_first * matrix_width + j] + g {
prev_i = pred_first;
prev_j = j;
predecessor_found = true;
} else {
for p in 1..node.inedges.len() {
let pred_i = pred_row(node.inedges[p]);
let a = h_ij == f[pred_i * matrix_width + j] + e_penalty;
extend_up = a;
if a || h_ij == h[pred_i * matrix_width + j] + g {
prev_i = pred_i;
prev_j = j;
predecessor_found = true;
break;
}
}
}
}
if !predecessor_found && j != 0 {
let a = h_ij == e[i * matrix_width + (j - 1)] + e_penalty;
extend_left = a;
if a || h_ij == h[i * matrix_width + (j - 1)] + g {
prev_i = i;
prev_j = j - 1;
}
}
let node_slot = if i == prev_i {
-1
} else {
graph.rank_to_node[i - 1].0 as i32
};
let seq_slot = if j == prev_j { -1 } else { (j - 1) as i32 };
alignment.push((node_slot, seq_slot));
i = prev_i;
j = prev_j;
if extend_left {
loop {
alignment.push((-1, (j - 1) as i32));
j -= 1;
if e[i * matrix_width + j] + e_penalty != e[i * matrix_width + (j + 1)] {
break;
}
}
} else if extend_up {
loop {
let mut stop = false;
prev_i = 0;
for &edge_id in &graph.nodes[graph.rank_to_node[i - 1].0 as usize].inedges {
let pred_i = pred_row(edge_id);
let s = f[i * matrix_width + j] == h[pred_i * matrix_width + j] + g;
stop = s;
if s || f[i * matrix_width + j] == f[pred_i * matrix_width + j] + e_penalty {
prev_i = pred_i;
break;
}
}
alignment.push((graph.rank_to_node[i - 1].0 as i32, -1));
i = prev_i;
if stop || i == 0 {
break;
}
}
}
}
alignment.reverse();
alignment
}
#[allow(clippy::too_many_arguments)]
#[inline]
pub(crate) fn backtrack_convex(
graph: &Graph,
node_id_to_rank: &[u32],
sequence_profile: &[i32],
h: &[i32],
e: &[i32],
f: &[i32],
o: &[i32],
q: &[i32],
matrix_width: usize,
alignment_type: AlignmentType,
scoring: &Scoring,
max_i: usize,
max_j: usize,
max_score: i32,
) -> Alignment {
if max_i == 0 && max_j == 0 {
return Vec::new();
}
debug_assert_eq!(
h[max_i * matrix_width + max_j],
max_score,
"fill's max_score must match H at (max_i, max_j)"
);
let g = i32::from(scoring.g);
let e_penalty = i32::from(scoring.e);
let q_penalty = i32::from(scoring.q);
let c_penalty = i32::from(scoring.c);
let pred_row = |edge_id: EdgeId| -> usize {
let tail = graph.edges[edge_id.0 as usize].tail;
node_id_to_rank[tail.0 as usize] as usize + 1
};
let mut alignment: Alignment = Vec::new();
let mut i = max_i;
let mut j = max_j;
loop {
let keep_going = match alignment_type {
AlignmentType::Local => h[i * matrix_width + j] != 0,
AlignmentType::Global => !(i == 0 && j == 0),
AlignmentType::Overlap => !(i == 0 || j == 0),
};
if !keep_going {
break;
}
let h_ij = h[i * matrix_width + j];
let mut prev_i = 0usize;
let mut prev_j = 0usize;
let mut predecessor_found = false;
let mut extend_left = false;
let mut extend_up = false;
if let Some((pi, pj)) = backtrack_match_step(
graph,
node_id_to_rank,
sequence_profile,
h,
matrix_width,
i,
j,
) {
prev_i = pi;
prev_j = pj;
predecessor_found = true;
}
if !predecessor_found && i != 0 {
let node = &graph.nodes[graph.rank_to_node[i - 1].0 as usize];
let pred_first = if node.inedges.is_empty() {
0
} else {
pred_row(node.inedges[0])
};
let a = h_ij == f[pred_first * matrix_width + j] + e_penalty;
extend_up |= a;
let cond = a
|| h_ij == h[pred_first * matrix_width + j] + g
|| {
let b = h_ij == o[pred_first * matrix_width + j] + c_penalty;
extend_up |= b;
b
}
|| h_ij == h[pred_first * matrix_width + j] + q_penalty;
if cond {
prev_i = pred_first;
prev_j = j;
predecessor_found = true;
} else {
for p in 1..node.inedges.len() {
let pred_i = pred_row(node.inedges[p]);
let a = h_ij == f[pred_i * matrix_width + j] + e_penalty;
extend_up |= a;
let cond = a
|| h_ij == h[pred_i * matrix_width + j] + g
|| {
let b = h_ij == o[pred_i * matrix_width + j] + c_penalty;
extend_up |= b;
b
}
|| h_ij == h[pred_i * matrix_width + j] + q_penalty;
if cond {
prev_i = pred_i;
prev_j = j;
predecessor_found = true;
break;
}
}
}
}
if !predecessor_found && j != 0 {
let a = h_ij == e[i * matrix_width + (j - 1)] + e_penalty;
extend_left |= a;
let cond = a
|| h_ij == h[i * matrix_width + (j - 1)] + g
|| {
let b = h_ij == q[i * matrix_width + (j - 1)] + c_penalty;
extend_left |= b;
b
}
|| h_ij == h[i * matrix_width + (j - 1)] + q_penalty;
if cond {
prev_i = i;
prev_j = j - 1;
}
}
let node_slot = if i == prev_i {
-1
} else {
graph.rank_to_node[i - 1].0 as i32
};
let seq_slot = if j == prev_j { -1 } else { (j - 1) as i32 };
alignment.push((node_slot, seq_slot));
i = prev_i;
j = prev_j;
if extend_left {
loop {
alignment.push((-1, (j - 1) as i32));
j -= 1;
let e_stops = e[i * matrix_width + j] + e_penalty != e[i * matrix_width + (j + 1)];
let q_stops = q[i * matrix_width + j] + c_penalty != q[i * matrix_width + (j + 1)];
if e_stops && q_stops {
break;
}
}
} else if extend_up {
loop {
let mut stop = true;
prev_i = 0;
for &edge_id in &graph.nodes[graph.rank_to_node[i - 1].0 as usize].inedges {
let pred_i = pred_row(edge_id);
if f[i * matrix_width + j] == f[pred_i * matrix_width + j] + e_penalty
|| o[i * matrix_width + j] == o[pred_i * matrix_width + j] + c_penalty
{
prev_i = pred_i;
stop = false;
break;
}
}
if stop {
for &edge_id in &graph.nodes[graph.rank_to_node[i - 1].0 as usize].inedges {
let pred_i = pred_row(edge_id);
if f[i * matrix_width + j] == h[pred_i * matrix_width + j] + g
|| o[i * matrix_width + j] == h[pred_i * matrix_width + j] + q_penalty
{
prev_i = pred_i;
break;
}
}
}
alignment.push((graph.rank_to_node[i - 1].0 as i32, -1));
i = prev_i;
if stop || i == 0 {
break;
}
}
}
}
alignment.reverse();
alignment
}
#[cfg(test)]
mod tests {
use super::*;
use crate::graph::Graph;
fn linear_scoring() -> Scoring {
Scoring::new(5, -4, -8, -8, -8, -8).unwrap()
}
fn affine_scoring() -> Scoring {
Scoring::new(5, -4, -8, -6, -8, -6).unwrap()
}
fn convex_scoring() -> Scoring {
Scoring::new(5, -4, -8, -6, -10, -4).unwrap()
}
#[test]
fn backtrack_linear_matches_single_node_on_hand_built_matrix() {
let mut g = Graph::new();
g.add_alignment_weight(&[], b"A", 1).unwrap();
let scoring = linear_scoring();
let node_id_to_rank = vec![0u32];
let sequence_profile = vec![0, 5];
let matrix_width = 2;
let h = vec![0, -8, -8, 5];
let alignment = backtrack_linear(
&g,
&node_id_to_rank,
&sequence_profile,
&h,
matrix_width,
AlignmentType::Global,
&scoring,
1,
1,
5,
);
assert_eq!(alignment, vec![(0, 0)]);
}
#[test]
fn backtrack_affine_unwinds_e_run_then_deletes() {
use crate::align::sisd::NEG_INF;
let mut g = Graph::new();
g.add_alignment_weight(&[], b"A", 1).unwrap();
let node_id = g.rank_to_node[0].0 as i32;
let scoring = affine_scoring(); let node_id_to_rank = vec![0u32];
let sequence_profile = vec![0, 5, 5];
let matrix_width = 3;
let h = vec![0, -8, -16, -8, NEG_INF, -26];
let e = vec![0, -8, -16, 0, -20, NEG_INF];
let f = vec![0, NEG_INF, NEG_INF, NEG_INF, NEG_INF, NEG_INF];
let alignment = backtrack_affine(
&g,
&node_id_to_rank,
&sequence_profile,
&h,
&e,
&f,
matrix_width,
AlignmentType::Global,
&scoring,
1,
2,
-26,
);
assert_eq!(alignment, vec![(node_id, -1), (-1, 0), (-1, 1)]);
}
#[test]
fn backtrack_convex_matches_single_node_on_hand_built_matrix() {
use crate::align::sisd::NEG_INF;
let mut g = Graph::new();
g.add_alignment_weight(&[], b"A", 1).unwrap();
let scoring = convex_scoring();
let node_id_to_rank = vec![0u32];
let sequence_profile = vec![0, 5];
let matrix_width = 2;
let h = vec![0, -8, -8, 5];
let e = vec![0, -8, NEG_INF, NEG_INF];
let f = vec![0, NEG_INF, NEG_INF, NEG_INF];
let o = vec![0, NEG_INF, NEG_INF, NEG_INF];
let q = vec![0, -10, NEG_INF, NEG_INF];
let alignment = backtrack_convex(
&g,
&node_id_to_rank,
&sequence_profile,
&h,
&e,
&f,
&o,
&q,
matrix_width,
AlignmentType::Global,
&scoring,
1,
1,
5,
);
assert_eq!(alignment, vec![(0, 0)]);
}
}