#include #include #include #include #include #include #include #define alphabets 26 #define SHIFTBITS 31 using namespace std; #ifndef _LOGGER_HPP_ #define _LOGGER_HPP_ #include #include /* 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 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 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 > &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> 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_(m, "NW") .def(py::init>, int, int, int>()) .def("get_score", &NW::get_score); }