aligner.h

Go to the documentation of this file.

/*
 * String aligner for the Phonetic Distance Package
 * + 
 * + April 2006 pcomas@lsi.upc.edu
 *
 */
#ifndef _aligner_h
#define _aligner_h

#include <iostream>
#include <fstream>
#include <sstream>
#include <map>
#include <set>
#include <string>
#include <math.h>

#include "freeling/morfo/phd.h"

#define GLOBAL    1  // The aligner produces local alignments
#define SEMILOCAL 2  // The aligner produces semilocal alignments
#define LOCAL     3  // The aligner produces global alignments

namespace freeling {

  template<typename T=int> 
    class aligner{

  private:
  phd<T>* sc;
  T score;
  int debug;

  public:

  struct alin {
    T score;   // score of this alignment
    double scoren; // score normalized
    double context;
    int Psubstitutions;
    int Pinserts;
    int Pdeletions;
    int Wsubstitutions;
    int Winserts;
    int Wdeletions;
    unsigned int kword; //The number of keyword involved in this alignment
    int begin; // index in b[] where the alignment starts
    int end;   // index in b[] where the alignment ends
    int beginW; //number of word in b[] where the alignment starts (each ' ' changes the word)
    int endW;   //number of word in b[] where the alignment ends
    wchar_t* seg; // Segment of b[] with the letters
    wchar_t* a;   // Alignment of a
    wchar_t* b;   // Alignment of b
    bool good; // Marks if this alignment is selected as relevant

    ~alin(){
      delete[](a);
      delete[](b);
      delete[](seg);
    }

  alin(T _score, int _begin, int _end, int _beginW, int _endW, int _substitutions, int _inserts, int _deletions, wchar_t* _a, wchar_t* _b) :
    score(_score), Psubstitutions(_substitutions), Pinserts(_inserts), Pdeletions(_deletions), begin(_begin), end(_end), beginW(_beginW), endW(_endW), seg(NULL), a(_a), b(_b) {}
  };


  aligner(const std::wstring &fname, int const _debug = 0){
    sc = new phd<T>(fname);
    debug=_debug;
    // if(debug>3){ sc->show(cerr);}
  } //constructor

  ~aligner(){
    delete(sc);
  }


  /*
   * Alinia la cadena A contra la cadena B. Conceptualment es considera que A és la query
   */
  alin* align(const std::wstring &wa, const std::wstring &wb, const int mode = SEMILOCAL){
    // a is the short string with length tj
    // b is the long string with length ti
    // The algorithm searchs for the best match of a against b in a semi-local or global point of view.
    // Usage of global matching is discouraged for the sake of coherence :-) because it doesn't try
    // to keep the letters together so 'a' will be meaningless scattered through all 'b'

    const wchar_t *a=wa.c_str(); const int tj=wa.size();
    const wchar_t *b=wb.c_str(); const int ti=wb.size();

    /*
    //  Build the align matrix
    */
    int const W = ti+1;

    if( ti==0 || tj == 0 ){ return new alin(0,0,0,0,0,0,0,0,0,0); }
    int i,j;
    int* m = new int[W*(tj+1)+ti+1]; //int m[tj+1][ti+1];

    // Variables for alignment reconstruction
    wchar_t* answerA = new wchar_t[tj+ti];
    wchar_t* answerB = new wchar_t[tj+ti];
    int pA = 0;
    int pB = 0;
    int insertions=0, deletions=0, substitutions=0;
    int spacesA=0;


    // INITIALIZATION STEP

    // Decide the number of word that each character belongs to
    int nwords=0;
    for(int i=0; i<ti; i++){ 
      if( b[i] == L' ' || b[i] == L'_' ){
        // word change at _i_
        nwords++;
      }
    }

    int* words = new int[nwords];
    j=0;
    for(int i=0; i<ti; i++){ 
      if( b[i] == L' ' || b[i] == L'_' ){
        words[j++] = i;
      }
    }

    switch(mode){
    case GLOBAL:  //Start values are skips
    m[0]=0;
    for(int j=1; j<=tj; j++){ m[j*W] = m[(j-1)*W] + sc->dSkip(a[j-1]); }
    for(int i=1; i<=ti; i++){ m[i]   = m[i-1]     + sc->dSkip(b[i-1]); }
    break;
    default: //Start values are 0
    for(int j=0;j<=tj;j++){ m[j*W] = 0; }
    for(int i=1;i<=ti;i++){ m[i] = 0; }
    break;
    }


    /*
    // Throw the scoring process
    */
    int i1, i2, i3, i4, i5, indexInit, indexEnd, bestJ, bestI;
    indexEnd  = 0;
    indexInit = 0;
    int initWord  = 0;
    int endWord   = 0;
    score = -100000000;
    bestJ = tj;
    bestI = ti;

    // FILL THE FIRST COLUMN & ROW
    for(j=1;j<=tj;j++){
      m[j*W+1] = std::max(  m[W*(j-1)]+sc->dSub(a[j-1],b[0]) , m[W*(j-1)+1]+sc->dSkip(a[j-1]) );
      if(a[j-1]==L' '|| a[j-1]==L'_' ){spacesA++;}
      if(score < m[W*j+1]){
        score = m[W*j+1];
        bestJ = j;
        bestI = 0;
      }
    }
    for(i=1;i<=ti;i++){
      m[ W+i ] = std::max(  m[ i-1 ]+sc->dSub(a[0],b[i-1])  ,  m[ W+i-1 ]+sc->dSkip(b[i-1])  );
      if(score < m[W+i]){
        score = m[W+i];
        bestJ = 0;
        bestI = i;
      }
    }

    //FILL THE REST OF THE MATRIX
    //int step = ti/4;
    int lowLimit = 0;
    int upperLimit = 0;
    int threshold = -10000000;
    if( mode == LOCAL ){ threshold = 0; }

    for(j=2;j<=tj;j++){
      //if( mode==GLOBAL && tj/ti> 0.75){
      // // If the scoring is GLOBAL we do not fill the entire matrix
      // lowLimit   = max(  2 , j-step );
      // upperLimit = min( ti , j+step );
      // m[W*j+lowLimit+1] = 0;
      //      } else {
      lowLimit   = 2;
      upperLimit = ti;
      //      }

      for(i=lowLimit;i<=upperLimit;i++){
        i1 = m[W*(j-1)+i-1] + sc->dSub(a[j-1],b[i-1]);
        i2 = m[W*j+i-1] + sc->dSkip(b[i-1]);
        i3 = m[W*(j-1)+i-2] + sc->dExp(a[j-1],b[i-2],b[i-1]);
        i4 = m[W*(j-1)+i] + sc->dSkip(a[j-1]);
        i5 = m[W*(j-2)+i-1] + sc->dExp(b[i-1],a[j-2],a[j-1]);
        m[W*j+i] = std::max(threshold,std::max(i1,std::max(i1,std::max(i2,std::max(i3,std::max(i4,i5))))));

        if(score < m[W*j+i]){
          score = m[W*j+i];
          bestJ = j;
          bestI = i;
        }

      }
    }

    //Prints the full score matrix if needed
    /*if(debug>5){
      cerr << endl << "\t";
      for(i=0;i<ti;i++){ cerr << "\t" << b[i]; }
      cerr << endl << "\t";
      for(j=0;j<=tj;j++){
      if(j>0) cerr << a[j-1] << "\t";
      for(i=0;i<=ti;i++){
      cerr << m[W*j+i] << "\t";
      }
      cerr << endl;
      }
      }*/


    /*
      || Reconstruction Step
    */
    bool final = true;    //Exit condition
    int lastScore = score;
    int steps = 0; //Number of steps uset to align 'a' with 'b'
    int lastChar = 0;


    int ni=0;
    int mymax=0;

    switch(mode){
    case SEMILOCAL:
      /* In semi-local alignment, the best-path reconstructions doesn't start from the 
       * botom-right corner of the matrix but from the best-scoring cell among last
       * column and last row scores.
       * In this implementation only the last row is taken in account since it is 
       * the "long" string
       */
      for(i=0;i<=ti;i++){
        if( mymax < m[W*tj+i] ){
          mymax = m[W*tj+i];
          ni=i;
        }
      }
      score = mymax;
      for(i=ti-1;i>=ni;i--){ // Carry on half of the string for debugging purpouses
        answerA[pA++]= L'-';
        answerB[pB++]= b[i];
        insertions++;
      }
      i=ni; j=tj;
      break;
    case GLOBAL: // now this is global alignment
      i = ti; j = tj;
      break;
    case LOCAL: 
      /* In LOCAL alignment, the best-path reconstruction starts with at the best scoring 
       * position in the whole matrix.
       */
      j = bestJ;
      i = bestI;
      break;
    }


    while( final ){

      /*if(debug>5){
        cerr << "Looking for " << m[ W*j +i] << " in ("<<j<<","<<i<<")"<<endl;
        cerr << "m[j-1][i-1]=" << m[ W*(j-1) +i-1] << ": dSub=" << sc->dSub(a[j-1],b[i-1]) << endl;
        cerr << "m[j][i-1]=" <<   m[ W*j +i-1] << ": dSkipB=" << sc->dSkip(b[i-1]) << endl;
        cerr << "m[j-1][i]=" <<   m[ W*(j-1) +i] << ": dSkipA=" << sc->dSkip(a[j-1]) << endl;
        if(j>1){ cerr << "m[j-1][i-2]=" << m[ W*(j-1) +i-2] << ": dExp=" << sc->dExp(a[j-1],b[i-2],b[i-1]) << endl; }
        if(i>1){ cerr << "m[j-2][i-1]=" << m[ W*(j-2) +i-1] << ": dExp=" << sc->dExp(b[i-1],a[j-2],a[j-1]) << endl; }
        }*/


      if( j>0 && i>0 && m[W*j+i] ==  m[W*(j-1)+i-1] + sc->dSub(a[j-1],b[i-1]) ){
        //Aliniar b[i] amb a[j]
        lastScore = sc->dSub(a[j-1],b[i-1]);
        i--; j--;
        //cerr << "Sub a[" << j << "]=" << a[j] << " b[" << i << "]=" << b[i] << " @ " << pA << endl;
        answerA[pA++] = a[j];
        answerB[pB++] = b[i];
        indexInit= i;
        indexEnd = std::max(i,indexEnd);
        if( a[j]!=b[i] ){ substitutions++;  }
        steps++;
        lastChar = steps;

      } else if ( i>0 && j>0 && m[W*j+i] == m[W*j+i-1] + sc->dSkip(b[i-1]) ) {
        lastScore = sc->dSkip(b[i-1]);
        i--;
        answerA[pA++] = L'-';
        answerB[pB++] = b[i];
        insertions++;
        if(steps==0 && b[i]==L' ' && b[i]==L'_' ){
          score = m[W*j+i];
        } else {
          steps++;
        }

      } else if( j>1 && i>0 && m[W*j+i] ==  m[W*(j-2)+i-1] + sc->dExp(b[i-1],a[j-2],a[j-1]) ){
        lastScore =  sc->dExp(b[i-1],a[j-2],a[j-1]);
        i--; j-=2;
        answerA[pA++] = a[j];
        answerA[pA++] = a[j+1];
        answerB[pB++] = b[i];
        answerB[pB++] = L'+';
        indexInit = i;
        indexEnd = std::max(i,indexEnd);
        deletions++;
        if( a[j]!=b[i] ){ substitutions++; }
        steps+=2;
        lastChar = steps;

      } else if( j>0 && i>0 && m[W*j+i] ==  m[W*(j-1)+i] + sc->dSkip(a[j-1]) ){
        lastScore = sc->dSkip(a[j-1]);
        j--;
        answerA[pA++] = a[j];
        answerB[pB++] = L'-';
        indexInit = i;
        indexEnd = std::max(i,indexEnd);
        deletions++;
        steps++;
        lastChar = steps;

      } else if ( i>1 && j>0 && m[W*j+i] == m[ W*(j-1) +i-2 ] + sc->dExp(a[j-1],b[i-2],b[i-1]) ) {
        lastScore = sc->dExp(a[j-1],b[i-2],b[i-1]);
        j--; i--;
        answerA[pA++] = a[j];
        answerA[pA++] = L'+';
        answerB[pB++] = b[i];
        answerB[pB++] = b[i-1];
        indexInit = i;
        indexEnd = std::max(i,indexEnd);
        insertions++;
        if( a[j]!=b[i] ){ substitutions++; }
        i--;
        steps+=2;
        lastChar = steps;

      } else if ( j==0 ){
        i--;
        answerA[pA++] = L'-';
        answerB[pB++] = b[i];
        insertions++;
        if(steps!=0 || b[i]!=L' ' || b[i]!=L'_' ) steps++;

      } else if ( i==0 ){
        j--;
        answerA[pA++] = a[j];
        answerB[pB++] = L'-';
        deletions++;

        if(steps!=0 || b[i]!=L' ' || b[i]!=L'_' ) steps++;
        lastChar = steps;
        
      } else {
        std::wcerr << L"BOINK! Error at "<< j << L"," << i << std::endl;
        break;
      }
      
      // This is the termination condition
      switch(mode){
      case SEMILOCAL:
      final = j!=0;
      if(!final){
        i--;
        while(i>=0){
          answerA[pA++] = L'-';
          answerB[pB++] = b[i--];
          insertions++;
        }
      }
      break;

      case GLOBAL:
      final = i!=0 || j!=0;
      break;

      case LOCAL:
      final = m[W*j+i]!=0;
      break;
      }
    }


    // Computing of the score
    // Score= avg.points/phoneme
    steps = lastChar;

    switch(mode){
    case GLOBAL:
      // average of similarity per phonem (without spaces?)
      score /= (pA-spacesA);
      break;
    case SEMILOCAL:
      //if(debug>2) cerr << "     Score = " << score << " Last Score = " << lastScore << " steps = " << steps <<  endl;
      //score = (score-lastScore) / steps;
      score = score / steps;
      break;
    case LOCAL:
      score = score / steps;
      break;
    }

    // Print the alignment
    wchar_t* newA = new wchar_t[pA+1];
    wchar_t* newB = new wchar_t[pB+1];
    newA[pA]=0;
    newB[pB]=0;
    --pA;
    --pB;

    for(int i=0; pA>=0; ++i, --pA){
      newA[i] = answerA[pA];
    }
    for(int i=0; pB>=0; ++i, --pB){
      newB[i] = answerB[pB];
    }

    delete[](answerA);
    delete[](answerB);
    delete[] m;

    // Finds in wich word does indexInit and indexEnd fall
    for(int i=0; i<nwords; i++){ 
      if( words[i] == indexInit ){
        initWord = i+1;
        break;
      } else if( words[i] < indexInit ) {
        initWord = i+1;
      } else if( words[i] > indexInit ) {
        initWord = i;
        break;
      }

    }

    for(int i=0; i<nwords; i++){ 
      if( words[i] == indexEnd ){
        endWord = i;
        break;
      } else if( words[i] < indexEnd ) {
        endWord = i+1;
      } else if( words[i] > indexEnd ) {
        endWord = i;
        break;
      }

    }

    delete[] words;

    return new alin(score,indexInit,indexEnd,initWord,endWord,substitutions,insertions,deletions,newA,newB);

  }

  };

} // namespace

#endif