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

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

file structure modifications. NOTE, this revision is not working, please wait for the next...

  • 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 295 2005-04-29 09:15:58Z 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 cpptools {
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 cpptools and namespace theplu
Note: See TracBrowser for help on using the repository browser.