source: trunk/yat/utility/Aligner.cc @ 3194

// $Id: Aligner.cc 3194 2014-04-22 09:24:56Z peter $

/*
  Copyright (C) 2012 Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 3 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include <config.h>

#include "Aligner.h"

#include "Matrix.h"

#include <cassert>
#include <limits>

namespace theplu {
namespace yat {
namespace utility {

  Aligner::Aligner(double gap, double open_gap)
    : vertical_gap_(gap), open_vertical_gap_(open_gap), horizon_gap_(gap),
      open_horizon_gap_(open_gap)
  {
  }


  Aligner::Aligner(double vertical_gap, double open_vertical_gap,
                   double horizon_gap, double open_horizon_gap)
    : vertical_gap_(vertical_gap), open_vertical_gap_(open_vertical_gap),
      horizon_gap_(horizon_gap), open_horizon_gap_(open_horizon_gap)
  {
  }


  double Aligner::needleman_wunsch(const Matrix& s)
  {
    Matrix m(s.rows()+1,s.columns()+1, -std::numeric_limits<double>::max());
    m(0, 0) = 0.0;
    // Init upper and left border of matrix
    for (size_t i=1; i<m.rows(); i++)
      m(i,0) = -i * vertical_gap_ - open_vertical_gap_;
    for (size_t i=1; i<m.columns(); i++)
      m(0,i) = -i * horizon_gap_ - open_horizon_gap_;
    (*this)(s, m);
    // return lower right corner element
    return m(s.rows(), s.columns());
  }


  double Aligner::smith_waterman(const Matrix& s)
  {
    Matrix m(s.rows()+1,s.columns()+1);
    (*this)(s, m);
    return max(m);
  }


  void Aligner::operator()(const Matrix& score, Matrix& x)
  {
    assert(x.columns() == score.columns()+1);
    assert(x.rows() == score.rows()+1);

    columns_ = x.columns();
    alignment_.resize(columns_*x.rows(), none);

    for (size_t i=1; i<x.rows(); ++i)
      for (size_t j=1; j<x.columns(); ++j) {
        // match (or mismatch)
        if (x(i-1, j-1) + score(i-1, j-1) > x(i,j)) {
          x(i,j) = x(i-1, j-1) + score(i-1, j-1);
          directions(i, j) = diagonal;
        }

        if (directions(i, j-1) == right) {
          if (x(i, j-1)-horizon_gap_ > x(i,j)) {
            x(i,j) = x(i,j-1) - horizon_gap_;
            directions(i,j) = right;
          }
        }
        else if(x(i,j-1)-horizon_gap_-open_horizon_gap_ > x(i,j)) {
          x(i,j) = x(i,j-1) - horizon_gap_ - open_horizon_gap_;
          directions(i,j) = right;
        }

        if (directions(i-1,j) == down) {
          if (x(i-1,j) - vertical_gap_ > x(i,j)) {
            x(i,j) = x(i-1,j) - vertical_gap_;
            directions(i,j) = down;
          }
        }
        else if (x(i-1,j) - vertical_gap_ - open_vertical_gap_ > x(i,j)) {
          x(i,j) = x(i-1,j) - vertical_gap_ - open_vertical_gap_;
          directions(i,j) = down;
        }
      }
  }


  const Aligner::direction& Aligner::alignment(size_t i, size_t j) const
  {
    return alignment_[i*columns_ + j];
  }


  Aligner::direction& Aligner::directions(size_t i, size_t j)
  {
    return alignment_[i*columns_ + j];
  }

}}} // of namespace utility, yat, and theplu