Changeset 256

Timestamp:

Mar 3, 2005, 7:33:01 PM (17 years ago)

Author:

Peter

Message:

added functionality to get information about best path

Location:

trunk/src

Files:

• Alignment.cc (modified) (2 diffs)
• Alignment.h (modified) (2 diffs)

trunk/src/Alignment.cc

@@ -3 +3 @@
 #include "Alignment.h"
 #include "matrix.h"
+
+#include <utility>
+#include <vector>
 
 namespace theplu {
@@ -9 +12 @@
 
   double NeedlemanWunsch(const gslapi::matrix& dot_matrix,
+                         std::vector<std::pair<size_t, size_t> >& path,
                          const double mismatch)
   {
     gslapi::matrix s(dot_matrix.rows()+1,dot_matrix.columns()+1);
-    s(0,0)=0;
+    // Init upper and left border of matrix 
     for (size_t i=1; i<s.rows(); i++) 
       s(i,0)=-i*mismatch;
     for (size_t i=1; i<s.columns(); i++) 
       s(0,i)=-i*mismatch;
+    // choice(i,j) tells us how we came to s(i,j). 1 is diagonal, 2
+    // vertical, and 3 horizontal,
+    gslapi::matrix choice(s.rows()-1,s.columns()-1);
 
+    // Calculating NeedlemanWunsch matrix
     for (size_t i=1; i<s.rows(); i++) 
       for (size_t j=1; j<s.columns(); j++){ 
-        if (s(i-1,j-1) + dot_matrix(i-1,j-1) > s(i-1,j) && 
-            s(i-1,j-1) + dot_matrix(i-1,j-1) > s(i,j-1))
+        if (s(i-1,j-1) + dot_matrix(i-1,j-1) > s(i-1,j)-mismatch && 
+            s(i-1,j-1) + dot_matrix(i-1,j-1) > s(i,j-1)-mismatch){
           s(i,j)=s(i-1,j-1) + dot_matrix(i-1,j-1);
-        else if (s(i-1,j) > s(i,j-1))
-          s(i,j)=s(i-1,j);
-        else
-          s(i,j)=s(i,j-1);
+          choice(i,j)=1;
+        }
+        else if (s(i-1,j) > s(i,j-1)){
+          s(i,j)=s(i-1,j)-mismatch;
+          choice(i,j)=2;
+        }
+        else{
+          s(i,j)=s(i,j-1)-mismatch;
+          choice(i,j)=3;
+        }
       }
+    // Going backwards to find best path
+    size_t i = s.rows()-1;
+    size_t j= s.columns()-1;
+    path.resize(0);
+    while (i && j){
+      if (choice(i,j)==1){
+        i--;
+        j--;
+        path.push_back(std::make_pair(i,j));
+      }
+      else if (choice(i,j)==2)
+        i--;
+      else
+        j--;
+    }
+
     return s(s.rows()-1,s.columns()-1);
   }

trunk/src/Alignment.h

@@ -3 +3 @@
 #ifndef _theplu_cpptools_alignment_
 #define _theplu_cpptools_alignment_
+
+#include <utility>
+#include <vector>
 
 namespace theplu {
@@ -27 +30 @@
   /// sequences. The latter is penalized with \a gap, which means if
   /// \a gap is large there will be no gaps in the best aligning (if
-  /// the sequences have the same length).  @return score from the
-  /// best aligning
+  /// the sequences have the same length). \a path contains what
+  /// elements in \A that is used in the best path. In case of
+  /// degeneracy only one path is returned.
+  ///
+  /// @return score from the best aligning
   ///
   double NeedlemanWunsch(const gslapi::matrix& A,
+                         std::vector<std::pair<size_t, size_t> >& path,
                          const double gap);