SmithWaterman.h

#ifndef _theplu_yat_utility_smith_waterman
#define _theplu_yat_utility_smith_waterman

// $Id: SmithWaterman.h 3338 2014-11-05 12:10:16Z peter $

/*
  Copyright (C) 2014 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 "Aligner.h"
#include "Matrix.h"

#include <cstddef> // size_t

namespace theplu {
namespace yat {
namespace utility {

  class SmithWaterman
  {
  public:
    explicit SmithWaterman(double gap=0, double open_gap=0);

    SmithWaterman(double vertical_gap, double open_vertical_gap,
                  double horizon_gap, double open_horizon_gap);

    const Aligner::Cigar& cigar(void) const;

    size_t position(void) const;

    const Matrix& score(void) const;

    double operator()(const Matrix& dot);

    template<typename RandomAccessIterator1, typename RandomAccessIterator2>
    double operator()(RandomAccessIterator1 reference_begin,
                      RandomAccessIterator1 reference_end,
                      RandomAccessIterator2 query_begin,
                      RandomAccessIterator2 query_end,
                      double mismatch=0);

  private:
    Aligner aligner_;
    Aligner::Cigar cigar_;
    size_t position_;
    Matrix score_;
  };

  // template implementaion

  template<typename RandomAccessIterator1, typename RandomAccessIterator2>
  double SmithWaterman::operator()(RandomAccessIterator1 reference_begin,
                                   RandomAccessIterator1 reference_end,
                                   RandomAccessIterator2 query_begin,
                                   RandomAccessIterator2 query_end,
                                   double mismatch)
  {
    Matrix dot(reference_end-reference_begin, query_end-query_begin, -mismatch);
    for (size_t i=0; i<dot.rows(); ++i)
      for (size_t j=0; j<dot.columns(); ++j)
        if (reference_begin[i] == query_begin[j])
          dot(i, j) = 1;
    return this->operator()(dot);
  }

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

#endif