Changeset 256 for trunk/src/Alignment.cc


Ignore:
Timestamp:
Mar 3, 2005, 7:33:01 PM (17 years ago)
Author:
Peter
Message:

added functionality to get information about best path

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/Alignment.cc

    r183 r256  
    33#include "Alignment.h"
    44#include "matrix.h"
     5
     6#include <utility>
     7#include <vector>
    58
    69namespace theplu {
     
    912
    1013  double NeedlemanWunsch(const gslapi::matrix& dot_matrix,
     14                         std::vector<std::pair<size_t, size_t> >& path,
    1115                         const double mismatch)
    1216  {
    1317    gslapi::matrix s(dot_matrix.rows()+1,dot_matrix.columns()+1);
    14     s(0,0)=0;
     18    // Init upper and left border of matrix
    1519    for (size_t i=1; i<s.rows(); i++)
    1620      s(i,0)=-i*mismatch;
    1721    for (size_t i=1; i<s.columns(); i++)
    1822      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);
    1926
     27    // Calculating NeedlemanWunsch matrix
    2028    for (size_t i=1; i<s.rows(); i++)
    2129      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){
    2432          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        }
    2943      }
     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
    3060    return s(s.rows()-1,s.columns()-1);
    3161  }
Note: See TracChangeset for help on using the changeset viewer.