Changeset 875

Timestamp:

Sep 19, 2007, 2:17:05 AM (16 years ago)

Author:

Peter

Message:

rewriting SW algo, making code clearer, and perhaps removing a few bugs...?

Location:

trunk

Files:

• test/alignment_test.cc (modified) (2 diffs)
• yat/utility/Alignment.cc (modified) (1 diff)

trunk/test/alignment_test.cc

@@ -114 +114 @@
     }
 
+  std::string a("AGGUUGUCCGUGGUGAGUUCGCA");
+  std::string b("GAGGUUGUCCGUGGUGAGUUCG");
+  utility::matrix m(a.size(), b.size());
+  for (size_t j=0; j<a.size(); ++j)
+    for (size_t k=0; k<b.size(); ++k)
+      m(j,k) = (a[j]==b[k] ? 1 : -1000);
+  double score=utility::SmithWaterman(m, 100, 100);
+  if (score!=21)
+    ok=false;
+
   // testing ssearch
   if (utility::ssearch("Hello", "Hll", 0.0, 1.0)!=2){
@@ -124 +134 @@
     ok=false;
 
+  if (ok)
+    *error << "Test is ok." << std::endl;
+  else
+    *error << "Test fails." << std::endl;
+    
+  if (error!=&std::cerr)
+    delete error;
+
   if (ok) {
     return 0;

trunk/yat/utility/Alignment.cc

@@ -101 +101 @@
                        double gap, double open_gap)
   {
+    enum direction { none, right, down, diagonal };
+
     // Calculating S-W matrix
-    matrix m(s.rows()+2,s.columns()+2);
-    for (size_t i=2; i<m.rows(); ++i) 
-      for (size_t j=2; j<m.columns(); ++j) 
-        m(i,j) = max(0.0, m(i-1,j-1)+s(i-2,j-2), 
-                     m(i-1,j)-gap-open_gap, m(i-2,j)-2*gap, 
-                     m(i,j-1)-gap-open_gap, m(i,j-2)-2*gap);
-
+    matrix m(s.rows()+1,s.columns()+1);
+    matrix array(m);
+    for (size_t i=1; i<m.rows(); ++i) 
+      for (size_t j=1; j<m.columns(); ++j){ 
+        m(i,j)=0;
+        array(i,j) = none;
+        if (m(i-1,j-1)+s(i-1,j-1)>m(i,j)){
+          m(i,j) = m(i-1,j-1)+s(i-1,j-1);
+          array(i,j) = diagonal;
+        }         
+        if (array(i-1,j)==right && m(i-1,j)-gap > m(i,j)){
+          m(i,j) = m(i-1,j)-gap;
+          array(i,j) = right;
+        }         
+        if (array(i-1,j)!=right && m(i-1,j)-gap-open_gap > m(i,j)){
+          m(i,j) = m(i-1,j)-gap-open_gap;
+          array(i,j) = right;
+        }         
+        if (array(i,j-1)==down && m(i,j-1)-gap > m(i,j)){
+          m(i,j) = m(i,j-1)-gap;
+          array(i,j) = down;
+        }         
+        if (array(i,j-1)!=down && m(i,j-1)-gap-open_gap > m(i,j)){
+          m(i,j) = m(i,j-1)-gap-open_gap;
+          array(i,j) = down;
+        }         
+      }
     return max(m);
   }