Skip to content

Commit

Permalink
slight tweaks, reformatting
Browse files Browse the repository at this point in the history
  • Loading branch information
@anesthetice
anesthetice committed Jul 16, 2024
1 parent 7b4152b commit 25cfdac
Show file tree
Hide file tree
Showing 5 changed files with 216 additions and 148 deletions.
2 changes: 1 addition & 1 deletion Cargo.toml
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
[package]
name = "syracuse"
version = "2.3.0"
version = "2.3.1"
edition = "2021"

[dependencies]
Expand Down
147 changes: 90 additions & 57 deletions src/algorithms.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,53 +9,67 @@ pub fn smith_waterman(seq_1: &str, seq_2: &str) -> f64 {
// initialization, seq_1 on the left and seq_2 on the top
let seq_1: Vec<char> = seq_1.chars().collect();
let seq_2: Vec<char> = seq_2.chars().collect();
let mut matrix : Vec<Vec<i16>> = vec![vec![0; seq_2.len()+1]; seq_1.len()+1];
let mut matrix: Vec<Vec<i16>> = vec![vec![0; seq_2.len() + 1]; seq_1.len() + 1];

// matrix filling
let mut total_score: i16 = 0;
let mut largest_score_index: (usize, usize) = (0,0);
for i in 1..seq_1.len()+1 { for j in 1..seq_2.len()+1 {
matrix[i][j] = {
let score = [
0,
matrix[i][j-1] + gap_penalty,
matrix[i-1][j] + gap_penalty,
if seq_1[i-1] == seq_2[j-1] {
matrix[i-1][j-1] + match_score
} else {
matrix[i-1][j-1] + mismatch_penalty
},
].into_iter().max().unwrap(); // safe unwrap
if score > total_score {
total_score = score;
largest_score_index = (i, j);
let mut largest_score_index: (usize, usize) = (0, 0);
for i in 1..seq_1.len() + 1 {
for j in 1..seq_2.len() + 1 {
matrix[i][j] = {
let score = [
0,
matrix[i][j - 1] + gap_penalty,
matrix[i - 1][j] + gap_penalty,
if seq_1[i - 1] == seq_2[j - 1] {
matrix[i - 1][j - 1] + match_score
} else {
matrix[i - 1][j - 1] + mismatch_penalty
},
]
.into_iter()
.max()
.unwrap(); // safe unwrap
if score > total_score {
total_score = score;
largest_score_index = (i, j);
}
score
}
score
}
}}
}

// traceback
let (mut i, mut j) = largest_score_index;
let mut score = total_score;
while score != 0 {
score = [
(i, j-1, matrix[i][j-1]),
(i-1, j, matrix[i-1][j]),
(i-1, j-1, matrix[i-1][j-1])
].into_iter().fold(0, |acc, (_i, _j, val)| {
(i, j - 1, matrix[i][j - 1]),
(i - 1, j, matrix[i - 1][j]),
(i - 1, j - 1, matrix[i - 1][j - 1]),
]
.into_iter()
.fold(0, |acc, (_i, _j, val)| {
if val > acc {
i = _i; j = _j;
i = _i;
j = _j;
val
} else {acc}
} else {
acc
}
});
total_score += score;
}

// normalisation
let max_score = {
let val = [match_score, mismatch_penalty, gap_penalty].into_iter().max().unwrap().abs();
let val = [match_score, mismatch_penalty, gap_penalty]
.into_iter()
.max()
.unwrap()
.abs();
let n = std::cmp::min(seq_1.len(), seq_2.len()) as i16 * val;
(n*(n+1))/2
(n * (n + 1)) / 2
};

total_score as f64 / max_score as f64
Expand All @@ -70,65 +84,84 @@ pub fn needleman_wunsch(seq_1: &str, seq_2: &str) -> f64 {
// initialization, seq_1 on the left and seq_2 on the top
let seq_1: Vec<char> = seq_1.chars().collect();
let seq_2: Vec<char> = seq_2.chars().collect();
let mut matrix : Vec<Vec<i16>> = vec![vec![0; seq_2.len()+1]; seq_1.len()+1];
let mut matrix: Vec<Vec<i16>> = vec![vec![0; seq_2.len() + 1]; seq_1.len() + 1];

// matrix filling
for j in 1..seq_2.len()+1 {matrix[0][j] = j as i16 * gap_penalty}
for i in 1..seq_1.len()+1 {matrix[i][0] = i as i16 * gap_penalty}
for i in 1..seq_1.len()+1 { for j in 1..seq_2.len()+1 {
matrix[i][j] = {
[
matrix[i][j-1] + gap_penalty,
matrix[i-1][j] + gap_penalty,
if seq_1[i-1] == seq_2[j-1] {
matrix[i-1][j-1] + match_score
} else {
matrix[i-1][j-1] + mismatch_penalty
}
].into_iter().max().unwrap() // safe unwrap
}
}}
for j in 1..seq_2.len() + 1 {
matrix[0][j] = j as i16 * gap_penalty
}
for i in 1..seq_1.len() + 1 {
matrix[i][0] = i as i16 * gap_penalty
}
for i in 1..seq_1.len() + 1 {
for j in 1..seq_2.len() + 1 {
matrix[i][j] = {
[
matrix[i][j - 1] + gap_penalty,
matrix[i - 1][j] + gap_penalty,
if seq_1[i - 1] == seq_2[j - 1] {
matrix[i - 1][j - 1] + match_score
} else {
matrix[i - 1][j - 1] + mismatch_penalty
},
]
.into_iter()
.max()
.unwrap() // safe unwrap
}
}
}

// traceback
let (mut i, mut j) = (seq_1.len(), seq_2.len());
let mut total_score: i16 = matrix[i][j];
while j != 0 && i != 0 {
if seq_1[i-1] == seq_2[j-1] {
i -= 1; j -= 1;
}
else {
if seq_1[i - 1] == seq_2[j - 1] {
i -= 1;
j -= 1;
} else {
// this only considers a single path, gives priority to going 'upwards'
if matrix[i][j-1] > matrix[i-1][j] {
if matrix[i][j - 1] > matrix[i - 1][j] {
j -= 1;
} else {i-=1}
} else {
i -= 1
}
}
total_score += matrix[i][j]
}
// goes to (0, 0) when reaching an edge
let i_gp = i as i16 * gap_penalty;
let j_gp = j as i16 * gap_penalty;
total_score += i_gp*(i_gp+1)/2 + j_gp*(j_gp+1)/2;
total_score += i_gp * (i_gp + 1) / 2 + j_gp * (j_gp + 1) / 2;

// normalisation
// paranoid, but this should mathematically ensure the result is between -1 and 1 even if the user does something absurd like set the match_score to a negative integer
// or if they have absurdly high penalites, etc...
match total_score {
1.. => {
let max_pos_score = {
let val = [match_score, mismatch_penalty, gap_penalty].into_iter().max().unwrap().abs();
let val = [match_score, mismatch_penalty, gap_penalty]
.into_iter()
.max()
.unwrap()
.abs();
let n = std::cmp::max(seq_1.len(), seq_2.len()) as i16 * val;
(n*(n+1))/2
(n * (n + 1)) / 2
};
total_score as f64 / max_pos_score as f64
},
}
0 => 0.0,
..=-1 => {
let max_neg_score = {
let val = [match_score, mismatch_penalty, gap_penalty].into_iter().min().unwrap().abs();
let val = [match_score, mismatch_penalty, gap_penalty]
.into_iter()
.min()
.unwrap()
.abs();
let n = std::cmp::max(seq_1.len(), seq_2.len()) as i16 * val;
(n*(n+1))/2
(n * (n + 1)) / 2
};
total_score as f64 / max_neg_score as f64
}
}
}
}
}
51 changes: 35 additions & 16 deletions src/animation.rs
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
use std::io::{self, Write};
use crossterm::style::Stylize;
use std::io::{self, Write};

use crate::warn;

Expand All @@ -12,28 +12,47 @@ pub struct Animation {
}

impl Animation {
pub fn construct(builder: AnimationBuilder, max_focus_len: usize, min_focus_len: usize) -> Self {
let max_frame_len = builder.iter().map(|(l, r)| l.len() + r.len()).max().unwrap_or(0_usize) + max_focus_len;
pub fn construct(
builder: AnimationBuilder,
max_focus_len: usize,
min_focus_len: usize,
) -> Self {
let max_frame_len = builder
.iter()
.map(|(l, r)| l.len() + r.len())
.max()
.unwrap_or(0_usize)
+ max_focus_len;

Self {
index: 0,
frames: builder.into_iter()
.map(|(l, r)| {
let num_repeats = max_frame_len - min_focus_len - l.len() - r.len();
("\r".to_string() + l.as_str(), r + " ".repeat(num_repeats).as_str())
})
.collect()
frames: builder
.into_iter()
.map(|(l, r)| {
let num_repeats = max_frame_len - min_focus_len - l.len() - r.len();
(
"\r".to_string() + l.as_str(),
r + " ".repeat(num_repeats).as_str(),
)
})
.collect(),
}
}
pub fn step(&mut self, stdout: &mut io::Stdout, focus: &str) {
if let Some((l, r)) = self.frames.get(self.index) {
let _ = stdout.write_all(
&[l.as_bytes(), focus.as_bytes(), r.as_bytes()].concat()
).map_err(|err| {warn!("failed to write animation to stdout, {}", err);});
let _ = stdout.flush().map_err(|err| {warn!("failed to flush stdout, {}", err);});
if self.index < self.frames.len()-1 {
let _ = stdout
.write_all(&[l.as_bytes(), focus.as_bytes(), r.as_bytes()].concat())
.map_err(|err| {
warn!("failed to write animation to stdout, {}", err);
});
let _ = stdout.flush().map_err(|err| {
warn!("failed to flush stdout, {}", err);
});
if self.index < self.frames.len() - 1 {
self.index += 1
} else {self.index = 0}
} else {
self.index = 0
}
}
}
}
}
Loading

0 comments on commit 25cfdac

Please sign in to comment.