Changeset 875


Ignore:
Timestamp:
Sep 19, 2007, 2:17:05 AM (14 years ago)
Author:
Peter
Message:

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

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/test/alignment_test.cc

    r870 r875  
    114114    }
    115115
     116  std::string a("AGGUUGUCCGUGGUGAGUUCGCA");
     117  std::string b("GAGGUUGUCCGUGGUGAGUUCG");
     118  utility::matrix m(a.size(), b.size());
     119  for (size_t j=0; j<a.size(); ++j)
     120    for (size_t k=0; k<b.size(); ++k)
     121      m(j,k) = (a[j]==b[k] ? 1 : -1000);
     122  double score=utility::SmithWaterman(m, 100, 100);
     123  if (score!=21)
     124    ok=false;
     125
    116126  // testing ssearch
    117127  if (utility::ssearch("Hello", "Hll", 0.0, 1.0)!=2){
     
    124134    ok=false;
    125135
     136  if (ok)
     137    *error << "Test is ok." << std::endl;
     138  else
     139    *error << "Test fails." << std::endl;
     140   
     141  if (error!=&std::cerr)
     142    delete error;
     143
    126144  if (ok) {
    127145    return 0;
  • trunk/yat/utility/Alignment.cc

    r869 r875  
    101101                       double gap, double open_gap)
    102102  {
     103    enum direction { none, right, down, diagonal };
     104
    103105    // Calculating S-W matrix
    104     matrix m(s.rows()+2,s.columns()+2);
    105     for (size_t i=2; i<m.rows(); ++i)
    106       for (size_t j=2; j<m.columns(); ++j)
    107         m(i,j) = max(0.0, m(i-1,j-1)+s(i-2,j-2),
    108                      m(i-1,j)-gap-open_gap, m(i-2,j)-2*gap,
    109                      m(i,j-1)-gap-open_gap, m(i,j-2)-2*gap);
    110 
     106    matrix m(s.rows()+1,s.columns()+1);
     107    matrix array(m);
     108    for (size_t i=1; i<m.rows(); ++i)
     109      for (size_t j=1; j<m.columns(); ++j){
     110        m(i,j)=0;
     111        array(i,j) = none;
     112        if (m(i-1,j-1)+s(i-1,j-1)>m(i,j)){
     113          m(i,j) = m(i-1,j-1)+s(i-1,j-1);
     114          array(i,j) = diagonal;
     115        }         
     116        if (array(i-1,j)==right && m(i-1,j)-gap > m(i,j)){
     117          m(i,j) = m(i-1,j)-gap;
     118          array(i,j) = right;
     119        }         
     120        if (array(i-1,j)!=right && m(i-1,j)-gap-open_gap > m(i,j)){
     121          m(i,j) = m(i-1,j)-gap-open_gap;
     122          array(i,j) = right;
     123        }         
     124        if (array(i,j-1)==down && m(i,j-1)-gap > m(i,j)){
     125          m(i,j) = m(i,j-1)-gap;
     126          array(i,j) = down;
     127        }         
     128        if (array(i,j-1)!=down && m(i,j-1)-gap-open_gap > m(i,j)){
     129          m(i,j) = m(i,j-1)-gap-open_gap;
     130          array(i,j) = down;
     131        }         
     132      }
    111133    return max(m);
    112134  }
Note: See TracChangeset for help on using the changeset viewer.