nw_wrapper.cpp

#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/numpy.h>
#include <iostream>
#include <string>
#include <vector>
#include <map>
#define alphabets 26
#define SHIFTBITS 31
using namespace std;
#ifndef _LOGGER_HPP_
#define _LOGGER_HPP_

#include <iostream>
#include <sstream>

/* consider adding boost thread id since we'll want to know whose writting and
 * won't want to repeat it for every single call */

/* consider adding policy class to allow users to redirect logging to specific
 * files via the command line
 */

enum loglevel_e
    {logERROR, logWARNING, logINFO, logDEBUG, logDEBUG1, logDEBUG2, logDEBUG3, logDEBUG4};

class logIt
{
public:
    logIt(loglevel_e _loglevel = logERROR) {
        _buffer << _loglevel << " :" 
            << std::string(
                _loglevel > logDEBUG 
                ? (_loglevel - logDEBUG) * 4 
                : 1
                , ' ');
    }

    template <typename T>
    logIt & operator<<(T const & value)
    {
        _buffer << value;
        return *this;
    }

    ~logIt()
    {
        _buffer << std::endl;
        // This is atomic according to the POSIX standard
        // http://www.gnu.org/s/libc/manual/html_node/Streams-and-Threads.html
        std::cerr << _buffer.str();
    }

private:
    std::ostringstream _buffer;
};

extern loglevel_e loglevel;

#define log(level) \
if (level > loglevel) ; \
else logIt(level)

#endif
class NW{
private:
//const size_t alphabets = 26;
size_t score=0;
int align(const std::string &aa, const std::string &bb, int alpha_gap,int match,int mismatch, 
        int alpha[alphabets][alphabets], std::string &a_aligned, 
        std::string &b_aligned)
{   
    size_t n = aa.size();
    size_t mm = bb.size();
    cout << n << "  tamaño de la prot1\n";
    cout << mm << "tamaño de la prot2\n"; 
    long int nRows = mm + 1;
  long int nCols = n + 1;
  int a, b, c, m; // temps for max calculation
  
  // populate first row, column
  // s1 across top, s2 down side
  vector<int> t;
  t.reserve(nRows*nCols+2);
  t[0] = 0;
  

  // populate first row and column
  for (long int i = 1; i < nCols; ++i) { t[i] = t[i-1] + alpha_gap; }
  for (long int i = 1; i < nRows; ++i) { t[nCols*i] = t[nCols*(i-1)] + alpha_gap; }
  
  cout << "tamaño de t  " << (nRows*nCols)+2 << "\n"; 
  // populate remaining table, row-major order
  for (long int i = 1; i < nRows; ++i) {
    for (long int j = 1; j < nCols; ++j) {
      m = -(aa[j-1] == bb[i-1]);
      //a = t[((i-1) * nCols) + j-1] + ((m? match:alpha[aa[j-1]][bb[i-1]]));
      a = t[((i-1) * nCols) + j-1] + ((alpha[aa[j-1]-'A'][bb[i-1]-'A']));
      b = t[((i-1) * nCols) + j] + alpha_gap;
      c = t[(i * nCols) + j-1] + alpha_gap;

      // OPTIMIZATION: bitwise max 
      a = a - (((a - b) >> SHIFTBITS) & (a - b));
      a = a - (((a - c) >> SHIFTBITS) & (a - c));
      t[(i * nCols) + j] = a;
    }
  }
  cout << "aqui muere?" << "\n";
  // printTable(t, nRows, nCols);
    //for (int i = nRows*nCols+1; i >= 0; i--) 
    //cout << t[(nRows)*(nCols-1)] << "\n";
    //cout << t[(nRows*nCols)-1] << "\n";
    return t[(nRows*nCols)-1];
}
void print2DVector(const vector<vector<int> > &A)
{
    for (auto& i : A)
    {
        for (auto j : i)
            cout << j << "\n";
        cout << endl;
    }
}
/* 
 * Returns the Needleman-Wunsch score for the best alignment of a and b
 * and stores the aligned sequences in a_aligned and b_aligned
 */
  public:


      NW(std::string a1, std::string b1,std::map<char,std::map<char,float>> c, int match, int mismatch, int gap)
      {
    // The input strings that need to be aligned
      
        

    // Penalty for any alphabet letter matched with a gap
       int gap_penalty = gap;

    /* 
     * alpha[i][j] = penalty for matching the ith alphabet with the
     *               jth alphabet.
     * Here: Penalty for matching an alphabet with anoter one is 1 
     *       Penalty for matching an alphabet with itself is 0
     */
       //cout << match << " + "<< mismatch << "\n";
       int alpha[alphabets][alphabets];
       for (size_t i = 0; i < alphabets; ++i)
       {
         for (size_t j = 0; j < alphabets; ++j)
         {
            //is first set as mismatch and continues like that if a match is not found
            alpha[i][j]=mismatch;
            //if alphabet letter is equal the other alphabet letter match
            if (i == j) alpha[i][j] = match;
            //tries to set as match alphabet letters that belong to the same group
            else{if(c.find((char)'A'+i) != c.end() || c.find((char)'A'+j) != c.end() ){
                   float count=0;
                   if(c.find((char)'A'+i) != c.end()){
                     
                     if(c[(char)'A'+i].find((char)'A'+j) != c[(char)'A'+i].end() )
                      count= float(c[(char)'A'+i][(char)'A'+j]);
                      //cout << float(c[(char)'A'+i][(char)'A'+j]);
                   }
                   
                   else{
                      if(c[(char)'A'+j].find((char)'A'+i) != c[(char)'A'+j].end())
                         count= float(c[(char)'A'+j][(char)'A'+i]);}
                   if(count !=0){
                     alpha[i][j]= match;
                     alpha[j][i]= match;} 
            
            }

         }
       }
     }
    // Aligned sequences
       std::string a2, b2;
       cout << "aqui se llega" << "\n";
       int penalty = align(a1, b1, gap_penalty,match,mismatch, alpha, a2, b2);
    
    this->score = penalty;
    //delete alpha;
}

int get_score(){
    return this->score;
}



};
namespace py = pybind11;
PYBIND11_MODULE(nw_wrapper, m) {
    py::class_<NW>(m, "NW")
        .def(py::init<std::string, std::string, std::map<char, std::map<char, float>>, int, int, int>())
        .def("get_score", &NW::get_score);
}