nicksenger / hirschberg

Generic implementation of Hirschberg's algorithm in Rust

Hirschberg

Hirschberg provides a generic implementation of Hirschberg's algorithm for finding the optimal global alignment of two sequences.

Example usage:

let a = "ACCCGGTCGTCAATTA".chars().collect::<Vec<_>>();
let b = "ACCACCGGTTGGTCCAATAA".chars().collect::<Vec<_>>();

let output = hirschberg::Config::default().compute(&a, &b);

let (aligned_a, aligned_b): (String, String) = output
    .alignment()
    .iter()
    .map(|(a, b)| (a.copied().unwrap_or('_'), b.copied().unwrap_or('_')))
    .unzip();

assert_eq!(output.score(), 8);
assert_eq!(
    [aligned_a, aligned_b],
    [
        "A_C_CCGG_TCGT_CAATTA".to_string(),
        "ACCACCGGTTGGTCCAATAA".to_string()
    ]
);

The match, mismatch and gap contribution to the score may be configured:

let a = "ACCCGGTCGTCAATTA".chars().collect::<Vec<_>>();
let b = "ACCACCGGTTGGTCCAATAA".chars().collect::<Vec<_>>();

let output = hirschberg::Config::default().mismatch_score(-4).compute(&a, &b);

let (aligned_a, aligned_b): (String, String) = output
    .alignment()
    .iter()
    .map(|(a, b)| (a.copied().unwrap_or('_'), b.copied().unwrap_or('_')))
    .unzip();

assert_eq!(output.score(), 6);
assert_eq!(
    [aligned_a, aligned_b],
    [
        "A_C_CCGG_T_CGT_CAAT_TA".to_string(),
        "ACCACCGGTTG_GTCCAATA_A".to_string()
    ]
);

hirschberg's Issues

Fraud! This is a Needleman wish and not a Hirschberg algorithm!

A proper hirschberg would look like this:

fn hirschberg(s1: &str, s1_len: usize, s2: &str, s2_len: usize) -> String {
    if s1_len == 0 || s2_len == 0 {
        return String::new();
    }

    if s1_len == 1 || s2_len == 1 {
        let mut lcs = String::new();
        for c in s1.chars() {
            if s2.contains(c) {
                lcs.push(c);
                break;
            }
        }
        return lcs;
    }

    let mid1 = s1_len / 2;
    let mid2 = s2_len / 2;

    let lcs_rev = hirschberg(&s1[mid1..].to_string().chars().rev().collect::<String>(), s1_len - mid1,
                             &s2[mid2..].to_string().chars().rev().collect::<String>(), s2_len - mid2);

    let lcs_fwd = hirschberg(&s1[..mid1].to_string(), mid1,
                             &s2[..mid2].to_string(), mid2);

    if lcs_fwd.len() > lcs_rev.chars().rev().collect::<String>().len() {
        lcs_fwd
    } else {
        lcs_rev.chars().rev().collect::<String>()
    }
}

fn main() {
    let s1 = "AGGTAB";
    let s2 = "GXTXAYB";
    let lcs = hirschberg(s1, s1.len(), s2, s2.len());
    println!("LCS: {}", lcs); // Output: LCS: GTAB
}

As can be seen, the Hirschberg algorithm is based on the principle of divide and conquer in order to enable a linear runtime on its data structure. Recursive calls are necessary for this. Your algorithm has a quadratic runtime on its data structure. This is completely understandable, as you have implemented a needleman algorithm. This is therefore cheating, as you are using the wrong algorithm.