3
\$\begingroup\$

This is a follow up of this question

Optimal substructure of ED:

Here is the reasoning behind my solution:

let \$x = (\alpha _{1},\alpha _{2},\alpha _{3},...,\alpha _{m})\$ and \$y = (\beta_{1},\beta_{2},\beta_{3},...,\beta_{n})\$.

We define \$c(i,j)\$ for \$ 0\leq i\leq m\$ and \$0\leq j\leq n\$ as the minimum edit distance to transform \$(\alpha _{1},...,\alpha _{i})\$ to \$(\beta_{1},...,\beta _{j})\$ without kill:

For \$1\leq i\leq m\$ and \$1\leq j\leq n\$:

\$c(i,j) = min\left\{\begin{matrix} c(i-1,j)+cost(delete) \phantom{-----------} \\ c(i,j-1)+cost(insert) \phantom{-----------} \\ \left\{\begin{matrix}c(i-1,j-1) + cost(copy) \hspace{5mm} if \hspace{2mm}\alpha _{i} = \beta _{j} \\ \infty \hspace{5cm} otherwise \end{matrix}\right. \phantom{----} \\ c(i-1,j-1)+cost(replace) \phantom{---------} \\\left\{\begin{matrix}c(i-2,j-2) + cost(twiddle)\phantom{-------}\newline if \hspace{2mm} i > 1 \hspace{2mm} \wedge \hspace{2mm}j > 1 \hspace{2mm} \wedge \hspace{2mm} \alpha _{i} = \beta _{j-1} \hspace{2mm} \wedge \hspace{2mm}\alpha _{i-1} = \beta _{j} \\ \infty \hspace{6cm} otherwise \end{matrix} \right.\phantom{--} \end{matrix} \right.\$

For \$1\leq i\leq m\$ and \$j = 0\$:

\$c(i,j)=c(i-1,j)\;+\cos t(delete)\$

For \$1\leq j\leq n\$ and \$i = 0\$:

\$c(i,j)=c(i,j-1)\;+\cos t(insert)\$

For \$i=0\$ and \$j = 0\$:

\$c(i,j)=0\$

Then let \$k(i,j)\$ be the minimum edit distance with kill:

\$k(i,j)= min\left\{\begin{matrix}c(i,j) \phantom{------------------} \\ min(\{c(1,j), c(2,j), ...,c(i-1,j)\})+cost(kill) \end{matrix}\right.\$

Question part 2:

This is taken from the book Introduction to Algorithms, Third Edition By Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest and Clifford Stein and is the second part of the same problem:

The edit-distance problem generalizes the problem of aligning two DNA sequences. There are several methods for measuring the similarity of two DNA sequences by aligning them. One such method to align two sequences \$x\$ and \$y\$ consists of inserting spaces at arbitrary locations in the two sequences (including at either end) so that the resulting sequences \$x'\$ and \$y'\$ have the same length but do not have a space in the same position (i.e., for no position \$j\$ are both \$x'[j]\$ and \$y'[j]\$ a space). Then we assign a “score” to each position. Position \$j\$ receives a score as follows:

  • \$+1\$ if \$x'[j] = y'[j]\$ and neither is a space,
  • \$-1\$ if \$x'[j] ≠ y'[j]\$ and neither is a space,
  • \$-2\$ if either \$x'[j]\$ or \$y'[j]\$ is a space.

The score for the alignment is the sum of the scores of the individual positions. For example, given the sequences x = GATCGGCAT and y = CAATGTGAATC, one alignment is

G ATCG GCAT
CAAT GTGAATC
- *++*+*+ -++*

A \$+\$ under a position indicates a score of \$+1\$ for that position, a \$-\$ indicates a score of \$-1\$, and a \$*\$ indicates a score of \$-2\$, so that this alignment has a total score of
\$6·1- 2·1 - 4·2 = 4\$.

b. Explain how to cast the problem of finding an optimal alignment as an edit distance problem using a subset of the transformation operations copy, replace,delete, insert, twiddle, and kill.

Code:

use std::ops::{Index, IndexMut};

#[derive(Debug)]
struct Matrix<T> {
    vec: Vec<T>,
    rows: usize,
    columns: usize,
}

impl<T: Clone> Matrix<T> {
    fn new(default: T, r: usize, c: usize) -> Matrix<T> {
        Matrix {
            vec: vec![default; r * c],
            rows: r,
            columns: c,
        }
    }
}

impl<T> Matrix<T> {
    fn check(&self, row: usize, column: usize) {
        assert!(row < self.rows, "row out of range");
        assert!(column < self.columns, "column out of range");
    }
}

impl<T> IndexMut<(usize, usize)> for Matrix<T> {
    fn index_mut<'a>(&'a mut self, index: (usize, usize)) -> &'a mut T {
        let (r, c) = index;
        self.check(r, c);
        &mut self.vec[r * self.columns + c]
    }
}

impl<T> Index<(usize, usize)> for Matrix<T> {
    type Output = T;
    fn index<'a>(&'a self, index: (usize, usize)) -> &'a T {
        let (r, c) = index;
        self.check(r, c);
        &self.vec[r * self.columns + c]
    }
}

#[derive(Debug, Copy, Clone)]
struct EditPosition {
    idx: usize,
    prev: char,
    current: char,
}

impl EditPosition {
    fn new(i: usize, p: char, c: char) -> EditPosition {
        EditPosition {
            idx: i,
            prev: p,
            current: c,
        }
    }
    fn valid_twiddle(&self, other: &EditPosition) -> bool {
        self.idx > 1 && other.idx > 1 && self.current == other.prev && self.prev == other.current
    }

    fn valid_copy(&self, other: &EditPosition) -> bool {
        self.current == other.current
    }
}

use std::collections::HashMap;

#[derive(Debug, Copy, Clone, PartialEq, Eq, Hash)]
enum EditOp {
    Copy,
    Replace,
    Delete,
    Insert,
    Twiddle,
    Kill,
}

impl EditOp {
    fn valid(&self, i: EditPosition, j: EditPosition) -> bool {
        use EditOp::*;
        match *self {
            Copy    => i.valid_copy(&j),
            Twiddle => i.valid_twiddle(&j),
            _ => true,
        }
    }

    fn displacement(&self) -> (usize, usize) {
        use EditOp::*;
        match *self {
            Copy | Replace => (1, 1),
            Delete         => (1, 0),
            Insert         => (0, 1),
            Twiddle        => (2, 2),
            Kill           => (0, 0),
        }
    }

    fn dec(&self, d: (usize, usize)) -> (usize, usize) {
        let (i, j) = self.displacement();
        (d.0 - i, d.1 - j)
    }

    fn align<F, T>(&self,
                   mut from_itr: F,
                   mut to_itr: T,
                   new_from: &mut String,
                   new_to: &mut String)
        where F: Iterator<Item = char>,
              T: Iterator<Item = char>
    {
        fn add_white_space<I>(a: &mut String, b: &mut String, mut a_itr: I)
            where I: Iterator<Item = char>
        {
            a.push(a_itr.next().unwrap());
            b.push(' ');
        }

        use EditOp::*;
        match *self {
            Copy | Replace => {
                new_from.push(from_itr.next().unwrap());
                new_to.push(to_itr.next().unwrap());
            }
            Delete => add_white_space(new_from, new_to, from_itr),
            Insert => add_white_space(new_to, new_from, to_itr),
            _ => panic!("Operation is not suported"),
        }
    }
}

type CostMap = HashMap<EditOp, isize>;
type EditTable = Matrix<(EditOp, isize)>;

fn create_table(from_len: usize, to_len: usize, costs: &CostMap) -> EditTable {
    use EditOp::*;
    let i_cost = *(costs.get(&Insert)
        .expect("Insert most be in costs map"));
    let d_cost = *(costs.get(&Delete)
        .expect("Delete most be in costs map"));


    let mut table = EditTable::new((Kill, 0), from_len, to_len);
    for i in 1..from_len {
        table[(i, 0)] = (Delete, table[(i - 1, 0)].1 + d_cost);
    }
    for j in 1..to_len {
        table[(0, j)] = (Insert, table[(0, j - 1)].1 + i_cost);
    }
    table
}


fn edit_distance(from: &str, to: &str, costs: &CostMap) -> (isize, Operations) {
    let from_len  = from.chars().count();
    let to_len    = to.chars().count();
    let mut table = create_table(from_len + 1, to_len + 1, costs);

    let (mut ilast, mut jlast) = (' ', ' ');
    for (i, ichar) in from.chars().enumerate() {
        for (j, jchar) in to.chars().enumerate() {
            let ij   = (i + 1, j + 1);
            let ipos = EditPosition::new(ij.0, ilast, ichar);
            let jpos = EditPosition::new(ij.1, jlast, jchar);

            table[ij] = costs.iter()
                .filter(|&(o, _)| o != &EditOp::Kill)
                .filter(|&(o, _)| o.valid(ipos, jpos))
                .map(|(o, c)| (*o, table[o.dec(ij)].1 + c))
                .min_by_key(|&(_, c)| c)
                .expect("no operation was valid");
            jlast = jchar;
        }
        ilast = ichar;
    }
    add_kill(&table, costs, from_len, to_len)
}

fn add_kill(table: &EditTable,
            costs: &CostMap,
            from_len: usize,
            to_len: usize)
            -> (isize, Operations) {
    let no_kill_cost    = table[(from_len, to_len)].1;
    let (i, cost, kill) = costs.get(&EditOp::Kill)
        .map_or((from_len, no_kill_cost, None),
                |kc| min_with_kill(&table, from_len, to_len, *kc));

    (cost, build_operations(kill, &table, i, to_len))
}

type Operations = Vec<EditOp>;

fn min_with_kill(table: &EditTable,
                 from_size: usize,
                 to_size: usize,
                 kill_cost: isize)
                 -> (usize, isize, Option<EditOp>) {
    let no_kill_cost = table[(from_size, to_size)].1;
    (1..from_size)
        .map(|i| (i, table[(i, to_size)].1 + kill_cost))
        .map(|(i, c)| (i, c, Some(EditOp::Kill)))
        .chain(Some((from_size, no_kill_cost, None)).into_iter())
        .min_by_key(|&(_, cost, _)| cost)
        .unwrap()
}

fn build_operations(kill: Option<EditOp>, table: &EditTable, i: usize, j: usize) -> Operations {
    let itr = std::iter::repeat(())
        .scan((i, j), |s, _| {
            let op = table[*s].0;
            *s     = op.dec(*s);
            Some(op)
        })
        .take_while(|op| op != &EditOp::Kill);

    let mut stack: Vec<_> = kill.into_iter()
        .chain(itr)
        .collect();
    stack.reverse();
    stack
}


fn optimal_alignment(from: &str, to: &str, cost_map: &CostMap) -> (String, String) {
    let (_, ops)     = edit_distance(from, to, cost_map);
    let mut new_from = String::new();
    let mut new_to   = String::new();
    let mut from_itr = from.chars();
    let mut to_itr   = to.chars();

    for op in ops {;
        op.align(&mut from_itr, 
                 &mut to_itr, 
                 &mut new_from, 
                 &mut new_to);
    }
    (new_from, new_to)
}

fn main() {
    let mut cost_map = HashMap::new();
    cost_map.insert(EditOp::Delete, 2);
    cost_map.insert(EditOp::Insert, 2);
    cost_map.insert(EditOp::Copy, -1);
    cost_map.insert(EditOp::Replace, 1);

    let from   = "CBADDKDDBBTPPEP";
    let to     = "CGAEEDDKDDBBTEEFGKKK";
    let (f, t) = optimal_alignment(from, to, &cost_map);
    println!("{}\n{}", f, t);
}
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.