Changeset 256
- Timestamp:
- Mar 3, 2005, 7:33:01 PM (18 years ago)
- Location:
- trunk/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/Alignment.cc
r183 r256 3 3 #include "Alignment.h" 4 4 #include "matrix.h" 5 6 #include <utility> 7 #include <vector> 5 8 6 9 namespace theplu { … … 9 12 10 13 double NeedlemanWunsch(const gslapi::matrix& dot_matrix, 14 std::vector<std::pair<size_t, size_t> >& path, 11 15 const double mismatch) 12 16 { 13 17 gslapi::matrix s(dot_matrix.rows()+1,dot_matrix.columns()+1); 14 s(0,0)=0;18 // Init upper and left border of matrix 15 19 for (size_t i=1; i<s.rows(); i++) 16 20 s(i,0)=-i*mismatch; 17 21 for (size_t i=1; i<s.columns(); i++) 18 22 s(0,i)=-i*mismatch; 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(s.rows()-1,s.columns()-1); 19 26 27 // Calculating NeedlemanWunsch matrix 20 28 for (size_t i=1; i<s.rows(); i++) 21 29 for (size_t j=1; j<s.columns(); j++){ 22 if (s(i-1,j-1) + dot_matrix(i-1,j-1) > s(i-1,j) &&23 s(i-1,j-1) + dot_matrix(i-1,j-1) > s(i,j-1) )30 if (s(i-1,j-1) + dot_matrix(i-1,j-1) > s(i-1,j)-mismatch && 31 s(i-1,j-1) + dot_matrix(i-1,j-1) > s(i,j-1)-mismatch){ 24 32 s(i,j)=s(i-1,j-1) + dot_matrix(i-1,j-1); 25 else if (s(i-1,j) > s(i,j-1)) 26 s(i,j)=s(i-1,j); 27 else 28 s(i,j)=s(i,j-1); 33 choice(i,j)=1; 34 } 35 else if (s(i-1,j) > s(i,j-1)){ 36 s(i,j)=s(i-1,j)-mismatch; 37 choice(i,j)=2; 38 } 39 else{ 40 s(i,j)=s(i,j-1)-mismatch; 41 choice(i,j)=3; 42 } 29 43 } 44 // Going backwards to find best path 45 size_t i = s.rows()-1; 46 size_t j= s.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 30 60 return s(s.rows()-1,s.columns()-1); 31 61 } -
trunk/src/Alignment.h
r185 r256 3 3 #ifndef _theplu_cpptools_alignment_ 4 4 #define _theplu_cpptools_alignment_ 5 6 #include <utility> 7 #include <vector> 5 8 6 9 namespace theplu { … … 27 30 /// sequences. The latter is penalized with \a gap, which means if 28 31 /// \a gap is large there will be no gaps in the best aligning (if 29 /// the sequences have the same length). @return score from the 30 /// best aligning 32 /// the sequences have the same length). \a path contains what 33 /// elements in \A that is used in the best path. In case of 34 /// degeneracy only one path is returned. 35 /// 36 /// @return score from the best aligning 31 37 /// 32 38 double NeedlemanWunsch(const gslapi::matrix& A, 39 std::vector<std::pair<size_t, size_t> >& path, 33 40 const double gap); 34 41
Note: See TracChangeset
for help on using the changeset viewer.