Changeset 256


Ignore:
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:
2 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  }
  • trunk/src/Alignment.h

    r185 r256  
    33#ifndef _theplu_cpptools_alignment_
    44#define _theplu_cpptools_alignment_
     5
     6#include <utility>
     7#include <vector>
    58
    69namespace theplu {
     
    2730  /// sequences. The latter is penalized with \a gap, which means if
    2831  /// \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
    3137  ///
    3238  double NeedlemanWunsch(const gslapi::matrix& A,
     39                         std::vector<std::pair<size_t, size_t> >& path,
    3340                         const double gap);
    3441
Note: See TracChangeset for help on using the changeset viewer.