source: trunk/lib/utility/Alignment.cc @ 301

Last change on this file since 301 was 301, checked in by Peter, 17 years ago

modified includes in tests

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 1.5 KB
Line 
1// $Id: Alignment.cc 301 2005-04-30 13:39:27Z peter $
2
3#include <c++_tools/utility/Alignment.h>
4#include <c++_tools/gslapi/matrix.h>
5
6#include <utility>
7#include <vector>
8
9namespace theplu {
10namespace utility {
11namespace alignment {
12
13  double NeedlemanWunsch(const gslapi::matrix& s,
14                         std::vector<std::pair<size_t, size_t> >& path,
15                         const double gap)
16  {
17    gslapi::matrix m(s.rows()+1,s.columns()+1);
18    // Init upper and left border of matrix
19    for (size_t i=1; i<m.rows(); i++) 
20      m(i,0)=-i*gap;
21    for (size_t i=1; i<m.columns(); i++) 
22      m(0,i)=-i*gap;
23    // choice(i,j) tells us how we came to s(i,j). 1 is diagonal, 2
24    // vertical, and 3 horizontal,
25    gslapi::matrix choice(m.rows(),m.columns());
26
27    // Calculating NeedlemanWunsch matrix
28    for (size_t i=1; i<m.rows(); i++) 
29      for (size_t j=1; j<m.columns(); j++){ 
30        if (m(i-1,j-1) + s(i-1,j-1) > m(i-1,j)-gap && 
31            m(i-1,j-1) + s(i-1,j-1) > m(i,j-1)-gap){
32          m(i,j)=m(i-1,j-1) + s(i-1,j-1);
33          choice(i,j)=1;
34        }
35        else if (m(i-1,j) > m(i,j-1)){
36          m(i,j)=m(i-1,j)-gap;
37          choice(i,j)=2;
38        }
39        else{
40          m(i,j)=m(i,j-1)-gap;
41          choice(i,j)=3;
42        }
43      }
44    // Going backwards to find best path
45    size_t i = m.rows()-1;
46    size_t j= m.columns()-1;
47    path.resize(0);
48    while (i && j){
49      if (choice(i,j)==1){
50        i--;
51        j--;
52        path.push_back(std::make_pair(i,j));
53      }
54      else if (choice(i,j)==2)
55        i--;
56      else
57        j--;
58    }
59
60    return m(m.rows()-1,m.columns()-1);
61  }
62
63}}} // of namespace alignment namespace utility and namespace theplu
Note: See TracBrowser for help on using the repository browser.