/*****************************************************************************
#   Copyright (C) 1994-2008 by David Gordon.
#   All rights reserved.                           
#                                                                           
#   This software is part of a beta-test version of the Consed/Autofinish
#   package.  It should not be redistributed or
#   used for any commercial purpose, including commercially funded
#   sequencing, without written permission from the author and the
#   University of Washington.
#   
#   This software is provided ``AS IS'' and any express or implied
#   warranties, including, but not limited to, the implied warranties of
#   merchantability and fitness for a particular purpose, are disclaimed.
#   In no event shall the authors or the University of Washington be
#   liable for any direct, indirect, incidental, special, exemplary, or
#   consequential damages (including, but not limited to, procurement of
#   substitute goods or services; loss of use, data, or profits; or
#   business interruption) however caused and on any theory of liability,
#   whether in contract, strict liability, or tort (including negligence
#   or otherwise) arising in any way out of the use of this software, even
#   if advised of the possibility of such damage.
#
#   Building Consed from source is error prone and not simple which is
#   why I provide executables.  Due to time limitations I cannot
#   provide any assistance in building Consed.  Even if you do not
#   modify the source, you may introduce errors due to using a
#   different version of the compiler, a different version of motif,
#   different versions of other libraries than I used, etc.  For this
#   reason, if you discover Consed bugs, I can only offer help with
#   those bugs if you first reproduce those bugs with an executable
#   provided by me--not an executable you have built.
# 
#   Modifying Consed is also difficult.  Although Consed is modular,
#   some modules are used by many other modules.  Thus making a change
#   in one place can have unforeseen effects on many other features.
#   It may takes months for you to notice these other side-effects
#   which may not seen connected at all.  It is not feasable for me to
#   provide help with modifying Consed sources because of the
#   potentially huge amount of time involved.
#
#*****************************************************************************/
#include    "rwcstring.h"
#include    <stdlib.h>
#include    "bool.h"
#include    "assert.h"
#include    <iostream.h>
#include    "mbtValOrderedVectorOfInt.h"
#include    "printAlignment.h"


static int nGap;
static int nMatch;
static int nMismatch;




#define MAX3( A, B, C ) ( (A) >= (B) && (A) >= (C) ) ? (A) : \
( (B) >= (C) ? (B) : (C) )
                        


typedef struct {
   int nX;
   int nY; } mbtPoint;





// used to convert nAcross from 0 based right to left coordinates
// to the coordinates used by pSeqAcross
inline int nParityX( const int nAcross, const mbtPoint ptEnd ) {
   return( ptEnd.nX - nAcross );
}

            
      
// used to convert nDown from 0 based bottom to top coordinates
// to coordinates used by pSeqDown
inline int nParityY( const int nDown, const mbtPoint ptEnd ) {
   return( ptEnd.nY - nDown );
}





inline RWCString soMakeStringFromInt( const int n ) {
   RWCString soTemp( (char) n );
   return( soTemp );
}






void getMidPointDivideByVerticalLine( const mbtPoint ptStart, 
                                      const mbtPoint ptEnd,
                                      const int nMidX,
                                      int& nBestY,
                                      int& nBestScore,
                                      char* pSeqAcross,
                                      char* pSeqDown ) {

   int nLengthDown = ptEnd.nY - ptStart.nY + 1;

   int* nScoreL = (int*) calloc( nLengthDown, sizeof( int ) );
   int* nScoreR = (int*) calloc( nLengthDown, sizeof( int ) );


   int nLengthAcrossL = nMidX - ptStart.nX + 1;
   int nLengthAcrossR = ptEnd.nX - nMidX + 1;

   // 1st column (nAcross = 0 )
   int nDown;
   for( nDown = 0; nDown < nLengthDown; ++nDown ) {
      nScoreL[ nDown ] = -nGap * nDown;
   }

   
   // other columsn (nAcross > 0 )
   int nAcross;
   for( nAcross = 1; nAcross < nLengthAcrossL; ++nAcross ) {
      
      int nToBeScoreDiagonallyAbove = nScoreL[0];
      
      // case of 1st row (nDown = 0)

      nScoreL[0] = nScoreL[0] - nGap;
      
      // now handle other rows (nDown > 0 )

      for( nDown = 1; nDown < nLengthDown; ++nDown ) {
         
         int nScoreDiagonallyAbove = nToBeScoreDiagonallyAbove;
         nToBeScoreDiagonallyAbove = nScoreL[ nDown ];

         int nScoreByDiagonal;

         // the -1 in each case is because the NW matrix
         // has its 0,0 before the start of the sequences
         // When moving from 0,0 to 1,1 you pass the
         // 0th char of each sequence

         if (pSeqAcross[ nAcross + ptStart.nX - 1 ] ==
             pSeqDown[ nDown + ptStart.nY - 1 ] ) 
            nScoreByDiagonal = nScoreDiagonallyAbove + nMatch;
         else
            nScoreByDiagonal = nScoreDiagonallyAbove - nMismatch;

         
         nScoreL[ nDown ] = MAX3(
                                 nScoreL[ nDown ] - nGap, // to left
                                 nScoreL[ nDown - 1 ] - nGap, // above
                                 nScoreByDiagonal
                                 );

      }


   }

   // so that is it for the left matrix
   // now compute the right matrix.  We are working from the bottom
   // right, but this time nDown goes up and nAcross goes to the left
   
   // 1st column (nAcross = 0 )
   for( nDown = 0; nDown < nLengthDown; ++nDown ) {
      nScoreR[ nDown ] = -nGap * nDown;
   }
    
   
   // subsequent columns (nAcross > 0 )

   for( nAcross = 1; nAcross < nLengthAcrossR; ++nAcross ) {
      int nToBeScoreDiagonallyBelow = nScoreR[ 0 ];


      // case of the bottom row (nDown = 0 )
      nScoreR[ 0 ] = nScoreR[ 0 ] - nGap;
      
      // other rows (nDown > 0 )

      for( nDown = 1; nDown < nLengthDown; ++nDown ) {
         int nScoreDiagonallyBelow = nToBeScoreDiagonallyBelow;
         nToBeScoreDiagonallyBelow = nScoreR[ nDown ];

         int nScoreByDiagonal;
         if ( pSeqAcross[ nParityX( nAcross, ptEnd ) ] ==
              pSeqDown[ nParityY( nDown, ptEnd ) ] ) 
            nScoreByDiagonal = nScoreDiagonallyBelow + nMatch;
         else
            nScoreByDiagonal = nScoreDiagonallyBelow - nMismatch;

         nScoreR[ nDown ] = MAX3(
                                 nScoreR[ nDown ] - nGap, // from right
                                 nScoreR[ nDown - 1 ] - nGap, // from below
                                 nScoreByDiagonal
                                 );

      }
   }

   // Now we have 2 columns of scores

   // Add them to find the point of maximum score, which will give us
   // a point that the optimal path goes through

   bool bBestScoreSet = false;
   int nBestScoreSoFar;
   int nBestDownSoFar;

   for( nDown = 0; nDown < nLengthDown; ++nDown ) {
      int nTotal = nScoreL[ nDown ] + nScoreR[ nLengthDown - 1 - nDown ];
      if ( !bBestScoreSet ) {
         bBestScoreSet = true;
         nBestScoreSoFar = nTotal;
         nBestDownSoFar = nDown;
      }
      else {
         if (nTotal > nBestScoreSoFar ) {
            nBestScoreSoFar = nTotal;
            nBestDownSoFar = nDown;
         }
      }
   }

   
   nBestScore = nBestScoreSoFar;
   nBestY = nBestDownSoFar + ptStart.nY;
}


         


void getMidPointDivideByHorizontalLine( const mbtPoint ptStart, 
                                      const mbtPoint ptEnd,
                                      const int nMidY,
                                      int& nBestX,
                                      int& nBestScore,
                                      char* pSeqAcross,
                                      char* pSeqDown ) {

   int nLengthWide = ptEnd.nX - ptStart.nX + 1;

   int nLengthHighT = nMidY - ptStart.nY + 1;
   int nLengthHighB = ptEnd.nY - nMidY + 1;

   int* nScoreT = (int*) calloc( nLengthWide, sizeof( int ) );
   int* nScoreB = (int*) calloc( nLengthWide, sizeof( int ) );

   // 1st row (nDown = 0 )
   int nAcross;
   for( nAcross = 0; nAcross < nLengthWide; ++nAcross ) {
      nScoreT[ nAcross ] = -nGap * nAcross;
   }

   
   // subsequent rows (nDown > 0 )
   int nDown;
   for( nDown = 1; nDown < nLengthHighT; ++nDown ) {
      
      int nToBeScoreDiagonallyAbove = nScoreT[0];
      
      // case of 1st column (nAcross = 0)

      nScoreT[0] = nScoreT[0] - nGap;
      
      // subsequent columns (nAcross > 0 )

      for( nAcross = 1; nAcross < nLengthWide; ++nAcross ) {
         
         int nScoreDiagonallyAbove = nToBeScoreDiagonallyAbove;
         nToBeScoreDiagonallyAbove = nScoreT[ nAcross ];

         int nScoreByDiagonal;
         // the -1 in each case is because the NW matrix
         // has its 0,0 before the start of the sequences
         // When moving from 0,0 to 1,1 you pass the
         // 0th char of each sequence
         if (pSeqAcross[ nAcross + ptStart.nX - 1] ==
             pSeqDown[ nDown + ptStart.nY - 1 ] ) 
            nScoreByDiagonal = nScoreDiagonallyAbove + nMatch;
         else
            nScoreByDiagonal = nScoreDiagonallyAbove - nMismatch;

         
         nScoreT[ nAcross ] = MAX3(
                                 nScoreT[ nAcross ] - nGap, // to left
                                 nScoreT[ nAcross - 1 ] - nGap, // above
                                 nScoreByDiagonal
                                 );

      }

   }

   // so that is it for the top matrix
   // now compute the bottom matrix.  We are working from the bottom
   // right, but this time nDown goes up and nAcross goes to the left


   // case of nDown == 0
   for( nAcross = 0; nAcross < nLengthWide; ++nAcross ) {
      nScoreB[ nAcross ] = -nGap * nAcross;
   }

   // subsequent rows (nDown > 0 

   for( nDown = 1; nDown < nLengthHighB; ++nDown ) {
      int nToBeScoreDiagonallyAbove = nScoreB[ 0 ];

      // case of 1st column (nAcross == 0 )

      nScoreB[0] = nScoreB[0] - nGap;

      // subsequent columns

      for( nAcross = 1; nAcross < nLengthWide; ++nAcross ) {
         int nScoreDiagonallyAbove = nToBeScoreDiagonallyAbove;
         nToBeScoreDiagonallyAbove = nScoreB[ nAcross ];

         int nScoreByDiagonal;
         if (pSeqAcross[ nParityX( nAcross, ptEnd ) ] ==
             pSeqDown[ nParityY( nDown, ptEnd ) ] )
            nScoreByDiagonal = nScoreDiagonallyAbove + nMatch;
         else
            nScoreByDiagonal = nScoreDiagonallyAbove - nMismatch;

         nScoreB[ nAcross ] = MAX3(
                                   nScoreB[ nAcross ] - nGap, // from above
                                   nScoreB[ nAcross - 1 ] - nGap, // from left
                                   nScoreByDiagonal
                                   );

      } // for (nAcross = 1 ...

   }    // for( nDown = 1 ...

   

   // now we have two rows of scores
   // Add them to find the point of maximum score, which will give us
   // a point that the optimal path goes through

   bool bBestScoreSet = false;
   int nBestScoreSoFar;
   int nBestAcrossSoFar;

   for( nAcross = 0; nAcross < nLengthWide; ++nAcross ) {
      int nTotal = nScoreT[ nAcross ] + nScoreB[ nLengthWide - 1 - nAcross ];

      if ( !bBestScoreSet ) {
         bBestScoreSet = true;
         nBestScoreSoFar = nTotal;
         nBestAcrossSoFar = nAcross;
      }
      else {
         if (nTotal > nBestScoreSoFar ) {
            nBestScoreSoFar = nTotal;
            nBestAcrossSoFar = nAcross;
         }
      }
   }

   
   nBestScore = nBestScoreSoFar;
   nBestX = nBestAcrossSoFar + ptStart.nX;
}


         


void getMidPoint( const mbtPoint ptStart, const mbtPoint ptEnd,
                  mbtPoint& ptMidPoint, int& nScore, char* pSeqAcross,
                  char* pSeqDown ) {
   

   if ( ptEnd.nY - ptStart.nY > 1 ) {
      int nMidY = ( ptEnd.nY + ptStart.nY ) / 2;
      int nBestX;

      getMidPointDivideByHorizontalLine( ptStart, ptEnd, nMidY, nBestX, nScore,
                                         pSeqAcross, pSeqDown );
      ptMidPoint.nX = nBestX;
      ptMidPoint.nY = nMidY;
   }
   else if ( ptEnd.nX - ptStart.nX > 1 ) {
      int nMidX = ( ptEnd.nX + ptStart.nX ) / 2;
      int nBestY;

      getMidPointDivideByVerticalLine( ptStart, ptEnd, nMidX, nBestY, nScore,
                                       pSeqAcross, pSeqDown );
      ptMidPoint.nX = nMidX;
      ptMidPoint.nY = nBestY;
   }
   else
      assert(false );


}




// this returns the alignment for the cases in which the figure is
// a line segment or a single cell (2 points by 2 points)

// Note that this does not return the upper left hand corner

void getPrimitiveAlignment( const mbtPoint ptStart,
                            const mbtPoint ptEnd,
                            char* pSeqAcross,
                            char* pSeqDown,
                            mbtValOrderedVectorOfInt& aAlignmentX,
                            mbtValOrderedVectorOfInt& aAlignmentY,
                            int& nScore) {

   aAlignmentX.clear();
   aAlignmentY.clear(); 
   nScore = 0;

   if (ptStart.nX == ptEnd.nX ) {
      for( int nDown = ptStart.nY + 1;
           nDown <= ptEnd.nY;
           ++nDown ) {
         aAlignmentX.append( ptStart.nX );
         aAlignmentY.append( nDown );
         nScore -= nGap;
      }
   }
   else if ( ptStart.nY == ptEnd.nY ) {
      for ( int nAcross = ptStart.nX + 1;
            nAcross <= ptEnd.nX;
            ++nAcross ) {
         
         aAlignmentX.append( nAcross );
         aAlignmentY.append( ptStart.nY );
         nScore -= nGap;
      }
   }
   else {
      // case in which the cell is 2 points by 2 points

      assert( ptStart.nX + 1 == ptEnd.nX );
      assert( ptStart.nY + 1 == ptEnd.nY );

         // the -1 in each case is because the NW matrix
         // has its 0,0 before the start of the sequences
         // When moving from 0,0 to 1,1 you pass the
         // 0th char of each sequence

      if (pSeqAcross[ ptEnd.nX - 1 ] == pSeqDown[ ptEnd.nY - 1 ] ) {
         aAlignmentX.append( ptEnd.nX );
         aAlignmentY.append( ptEnd.nY );
         nScore = nMatch;
      }
      else {
         if (nMismatch > (2*nGap) ) {
            // so prefer a gap
            aAlignmentX.append( ptStart.nX + 1 );
            aAlignmentY.append( ptStart.nY );

            aAlignmentX.append( ptEnd.nX );
            aAlignmentY.append( ptEnd.nY );
            nScore = -2*nGap;
         }
         else {
            // prefer a mismatch
            
            aAlignmentX.append( ptEnd.nX );
            aAlignmentY.append( ptEnd.nY );
            nScore = -nMismatch;
         }
      } //       if (pSeqAcross[ ptEnd.nX ] == pSeqDown[ ptEnd.nY ] ) {
   } //    if (ptStart.nX == ptEnd.nX ) {
} // getPrimitiveAlignment
   




void alignmentFromNWPath( const mbtValOrderedVectorOfInt& aAlignmentX,
                          const mbtValOrderedVectorOfInt& aAlignmentY,
                          char* pSeqAcross,
                          char* pSeqDown,
                          RWCString& soAlignmentTop,
                          RWCString& soAlignmentBottom )
{
   
   soAlignmentTop = "";
   soAlignmentBottom = "";

   assert( aAlignmentX.length() == aAlignmentY.length() );

   for( int n = 1; n < aAlignmentX.length(); ++n ) {
      if ( aAlignmentX[ n-1 ] == aAlignmentX[ n ] ) {
         // a vertical line, so a gap in the top sequence
         soAlignmentTop += "*";
         soAlignmentBottom += soMakeStringFromInt(
                         pSeqDown[ (int) aAlignmentY[ n-1 ] ]
                                                  );
      }
      else if (aAlignmentY[ n-1 ] == aAlignmentY[ n ] ) {
         // case of a horizontal line, so there is a gap in
         // the bottom sequence
         
         soAlignmentTop += soMakeStringFromInt(
                             pSeqAcross[ (int) aAlignmentX[ n-1 ] ]
                                               );
         soAlignmentBottom += "*";

      }
      else {
         // diagonal line so no gap.  The next base in the top
         // is aligned with the next base in the bottom.

         soAlignmentTop += soMakeStringFromInt(
                             pSeqAcross[ (int) aAlignmentX[ n-1 ] ]
                                               );

         soAlignmentBottom += soMakeStringFromInt(
                             pSeqDown[ (int) aAlignmentY[ n-1 ] ]
                                                  );

      }
   }
}



void checkAlignmentPath( const mbtValOrderedVectorOfInt& aAlignmentX,
                         const mbtValOrderedVectorOfInt& aAlignmentY,
                         char* pSeqAcross,
                         char* pSeqDown) {

   // This checks that the path goes diagonally from the upper
   // left hand corner to the lower right hand corner.

   assert( aAlignmentX.length() == aAlignmentY.length() );


   // case of n == 0
   assert( aAlignmentX[0] == 0 );
   assert( aAlignmentY[0] == 0 );

   for( int n = 1; n < aAlignmentX.length(); ++n ) {
      assert( ( aAlignmentX[n] == aAlignmentX[ n-1 ] + 1 )
              ||
              ( aAlignmentX[n] == aAlignmentX[ n-1 ] )
              );

      assert( ( aAlignmentY[n] == aAlignmentY[ n-1 ] + 1 )
              ||
              ( aAlignmentY[n] == aAlignmentY[ n-1 ] )
              );
   }
      

   // Now make sure that the path ends at the lower right-hand corner

   assert( (int) aAlignmentX[ aAlignmentX.length() - 1 ] ==
           strlen( pSeqAcross ) );

   assert( (int) aAlignmentY[ aAlignmentX.length() - 1 ] ==
           strlen( pSeqDown ) );

}



void checkAlignmentScore( const RWCString& soAlignmentTop,
                          const RWCString& soAlignmentBottom,
                          const int nSupposedScore ) {

   assert( soAlignmentTop.length() == soAlignmentBottom.length() );

   int nScore = 0;

   for( int n = 0; n < soAlignmentTop.length(); ++n ) {
      if (soAlignmentTop[n] == '*' )
         nScore -= nGap;
      else if ( soAlignmentBottom[n] == '*' )
         nScore -= nGap;
      else {
         if (soAlignmentTop[n] == soAlignmentBottom[n] )
            nScore += nMatch;
         else
            nScore -= nMismatch;
      }
   }


   assert( nScore == nSupposedScore );

}
      
            


void checkSameSequences2( const RWCString& soStringWithPads,
                          char *pSeq ) {

   int nSeq = -1;
   for( int n = 0; n < soStringWithPads.length(); ++n ) {
      if ( soStringWithPads[n] != '*' ) {
         ++nSeq;
         assert( pSeq[ nSeq ] == soStringWithPads[ n ] );
      }
   }


   // make sure all of pSeq is in soStringWithPads

   assert( nSeq == ( strlen( pSeq ) - 1 ) );
}


   
void checkSameSequences( const RWCString& soAlignmentTop,
                         const RWCString& soAlignmentBottom,
                         char* pSeqAcross,
                         char* pSeqDown ) {

   int nLength = strlen( pSeqAcross );

   checkSameSequences2( soAlignmentTop, pSeqAcross );
   checkSameSequences2( soAlignmentBottom, pSeqDown );
}




void getAlignment2( const mbtPoint ptStart, const mbtPoint ptEnd, 
                    char* pSeqAcross,
                    char* pSeqDown,
                    mbtValOrderedVectorOfInt& aPathX, mbtValOrderedVectorOfInt& aPathY,
                    int& nScore) {


   if (ptStart.nX == ptEnd.nX )
      getPrimitiveAlignment( ptStart, ptEnd, pSeqAcross, pSeqDown,
                             aPathX, aPathY, nScore );
   else if (ptStart.nY == ptEnd.nY )
      getPrimitiveAlignment( ptStart, ptEnd, pSeqAcross, pSeqDown,
                             aPathX, aPathY, nScore );
   else if ( (ptEnd.nX - ptStart.nX == 1 ) &&
             (ptEnd.nY - ptStart.nY == 1 ) )
      getPrimitiveAlignment( ptStart, ptEnd, pSeqAcross, pSeqDown,
                             aPathX, aPathY, nScore );
   else {
      mbtPoint ptMiddle;
      getMidPoint( ptStart, ptEnd, ptMiddle, nScore, pSeqAcross,
                   pSeqDown );


      mbtValOrderedVectorOfInt aPathFirstX;
      mbtValOrderedVectorOfInt aPathFirstY;
      int nPartialScore;
      getAlignment2( ptStart, ptMiddle, pSeqAcross, pSeqDown,
                     aPathFirstX, aPathFirstY, nPartialScore );

      mbtValOrderedVectorOfInt aPathLastX;
      mbtValOrderedVectorOfInt aPathLastY;
      getAlignment2( ptMiddle, ptEnd, pSeqAcross, pSeqDown,
                     aPathLastX, aPathLastY, nPartialScore );

      aPathX = aPathFirstX + aPathLastX;
      aPathY = aPathFirstY + aPathLastY;
   }

   assert( aPathX.length() == aPathY.length() );
}


      
   
void getAlignment( char* pSeqAcross, char* pSeqDown, 
                  const int nGapPenalty,
                  const int nMatchBoost,
                  const int nMismatchPenalty,
                  RWCString& soAlignmentTop,
                  RWCString& soAlignmentBottom,
                  int& nScore ) {

   
   nGap = nGapPenalty;
   nMatch = nMatchBoost;
   nMismatch = nMismatchPenalty;



   mbtPoint ptEnd;
   ptEnd.nX = strlen( pSeqAcross );
   ptEnd.nY = strlen( pSeqDown );

   mbtPoint ptStart;
   ptStart.nX = 0;
   ptStart.nY = 0;

   mbtValOrderedVectorOfInt aPathX;
   mbtValOrderedVectorOfInt aPathY;

   getAlignment2( ptStart, ptEnd, pSeqAcross, pSeqDown, aPathX, aPathY, 
                  nScore );

   // start the path at the upper left-hand corner
   // create a string with a single null in it

   mbtValOrderedVectorOfInt soOneNull;
   soOneNull.append( 0 );


   mbtValOrderedVectorOfInt aFullPathX = soOneNull + aPathX;
   mbtValOrderedVectorOfInt aFullPathY = soOneNull + aPathY;
   aPathX = aFullPathX;
   aPathY = aFullPathY;

   checkAlignmentPath( aPathX, aPathY, pSeqAcross, pSeqDown );

   alignmentFromNWPath( aPathX, aPathY, pSeqAcross, pSeqDown, 
                        soAlignmentTop, soAlignmentBottom );
   checkAlignmentScore( soAlignmentTop, soAlignmentBottom, nScore );

   checkSameSequences( soAlignmentTop, soAlignmentBottom, pSeqAcross, pSeqDown );

}







   


/*


int main() {
   RWCString soSequence1;
   RWCString soSequence2;

   while( true ) {

      RWCString soAlignmentTop;
      cout << "enter seq 1: ";
      cin >> soSequence1;

// "CTGATACATATTTAAAGCAGCACCCAAGTGTGTTCTAATAGAAATGCTGTGACCCTGAGGTCCTGGGGATTGAGAGAGGAAGTGATGTCACTGTGGGAACTGCCCTGTGGAGACAAGGACGGCCCTTATCCTCTGCTTCTGTTCACAGTGACACTGATCTGGTAAAGCCCCCATCCTGGCCTGACCCTGCCATGGGCACCAGGCTCCTCTGCTGGGTGGTCCTGGGTTTCCTAGGGACAGGTGAGTCCTCAGAACACCAAGTAGTTTCATTTTTTCTGTTTGTAGGTGTGTGTGTGTGAGAGAGTGGTTGTGTGTGTGTGTGTGTATGATGACTACAAATATTTTCCTTATTCTGTTGCCAAATTCTATTTCCACAGATCACACAGGTGCTGGAGTCTCCCAGTCCCCTAGGTACAAAGTCGCAAAGAGAGGACAGGATGTAGCTCTCAGGTGTGATCCAATTTCGGGTCATGTATCCCTTTTTTGGTACCAACAGGCCCTGGGGCAGGGGCCAGAGTTTCTGACTTATTTCCAGAATGAAGCTCAACTAGACAAATCGGGGCTGCCCAGTGATCGCTTCTTTGCAGAAAGGCCTGAGGGATCCGTCTCCACTCTGAAGATCCAGCGCACACAGCAGGAGGACTCCGCCGTGTATCTCTGTGCCAGCAGCTTAGCCACAGCATGGCACAGTTGCCTCCTTCCTGTTCACAAACCTCATCCTTCTCTCTCCTTGCAGCTCCTAGAGACCCTAAACAGAGGCCTCTCTTTGCTCCTCACTTTTCATGGGAAAGAATTACATCTGGACTTCAGCTGTTCTTTGGGTAGAAAGAGACCACAGATTCATTCCTGAAACACAGTGACTGAAAATGTAGGTGGTGAAAACAATCATGGGAGTCCTTGGA";

      RWCString soAlignmentBottom;
      cout << "enter seq 2: ";
      cin >> soSequence2;
// "CTGATAAATATTTAAAGCAGCACCCAACTGTGTTCTAATAGAAATGCTGTGATCCTGAGGTCCTGGGGATTGAGAGAGGAAGTGATGTCACTGTGGGAACTGCCCTGTGGAGACAAGGACATCCCTCATCCTCTGCTGCTGCTCACAGTGACACTGATCTGGTAAAGCCCTCATCCTGTCCTGACCCTGCCATGGGCACCAGTCTCCTATGCTGGGTGGTCCTGGGTTTCCTAGGGACAGGTGAGTCCTCAGAACACAAAGTAGTTTCATTTTTTTTTCTGTGTGTAGGCGTGTGGGCGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGCTGACGACAAATGTTTTCCTTATTCTGTTGCCAAATTCTATTTCCACAGATCACACAGGTGCTGGAGTCTCCCAGTCTCCCAGGTACAAAGTCACAAAGAGGGGACAGGATGTAGCTCTCAGGTGTGATCCAATTTCGGGTCATGTATCCCTTTATTGGTACCGACAGCCCTGGGGCAGGGCCCAGAGTTTCTGACTTACTTCAATTATGAAGCCCAACAAGACAAATCAGGGCTGCCCAATGATCGGTTCTCTGCAGAGAGGCCTGAGGGATCCATCTCCACTCTGACGATCCAGCGCACAGAGCAGCGGGACTCGGCCATGTATCGCTGTGCCAGCAGCTTAGCCACAGTGTGGCATAGTCGCCTCCTTCCTGTTCACAAACCTCATCCTTCTCTCTCCTTGCACCTCCTAGAGACCCTTAACAGAGGCCTCTCTTTGCTCCTCACTTTTGATGGGAAAGAAGTAGATTTGGACATCGGCTGTCCTTTGGGTAGAAAGAGACCACAGATTCATTCCTGAAACACAGTGACTGCAAATGTAGGTGGTGAAAACAATCACGTCCCACTGCCCTCTAGGAGGGCTCGGA";

      getAlignment( (char*) soSequence1.data(),
                    (char*) soSequence2.data(),
                    5,
                    2,
                    2,
                    soAlignmentTop,
                    soAlignmentBottom );


   printAlignment( soAlignmentTop, soAlignmentBottom );



   }
}
   
*/