/*****************************************************************************
#   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    <cstdlib>
using namespace std;
#include    <cstdio>
#include    "extendAlignmentBanded.h"


#ifndef WALRUS_BUILD
#include    <malloc.h>
#endif

#include    "numutil.h"
#include    <iostream>
#include    "assert.h"
#include    <cstring>
#include    "setMatchMatrix.h"
#include    "max.h"
#include    "min.h"
#include    "consedParameters.h"


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


void extendAlignmentBanded( 
    char* pSeqDown, 
    char* pSeqAcross, 
    int nSeqDownLength, 
	int nSeqAcrossLength,  
    int* pnBestDown, 
    int* pnBestAcross, 
	int* pnBestScore,
    const int nGapPenalty,
    const int nMatchBoost,
    const int nMismatchPenalty,
    const int nBandSize ) {

  int nBestScore;
  int nBestDown;
  int nBestAcross;
  int nScoreDiagonallyAbove;
  int nToBeScoreDiagonallyAbove;
  int nScoreByDiagonal;
  int nAcross;
  int nDown;
  int* nScore;
  int nMinAcross;
  int nMaxAcross;



  setMatchMatrix( nMatchBoost,
                  nMismatchPenalty );

  // the +1 is because the size of the NW matrix is 1 greater
  // than that of the sequence (in either dimension)
  nScore = (int*) calloc( nSeqAcrossLength + 1, sizeof(int) );

  nBestScore = 0;
  nBestDown = 0;
  nBestAcross = 0;

  // scores for for nDown = 0.
  // This could be removed to allow free slippage of the pSeqAcross
  // sequence to the right (but not to the left) with respect
  // to pSeqDown sequence

  
  nDown = 0;
  nMinAcross = 0; // would be nDown - nBandSize, but bumped into left edge
  nMaxAcross = MIN( nBandSize, nSeqAcrossLength); // nDown + nBandSize
  // keep in mind that the bandsize could be greater than nSeqAcrossLength

  for ( nAcross = 0; nAcross <= nMaxAcross; ++nAcross ) {
    nScore[ nAcross ] = -nGapPenalty * nAcross;

  }

  for( nDown = 1; nDown <= nSeqDownLength; ++nDown ) {
    nMinAcross = MAX( nDown - nBandSize, 0 );
    nMaxAcross = MIN( nDown + nBandSize, nSeqAcrossLength );

    if ( nMinAcross > nMaxAcross )
       break;


    nToBeScoreDiagonallyAbove = nScore[nMinAcross];

    /* case of first column */

    if ( nMinAcross == 0 ) {
       // if bumped into left edge
       nScore[0] = nScore[0] - nGapPenalty;
    }
    else {
       nScore[nMinAcross] = 
          MAX( nScore[nMinAcross] - nGapPenalty, // move down
                nScore[nMinAcross - 1 ] + 
                nMatchMatrix[ pSeqAcross[nMinAcross - 1]][ pSeqDown[ nDown - 1 ]] );

    }

    if (nScore[nMinAcross] > nBestScore ) {
       nBestScore = nScore[nMinAcross];
       nBestAcross = nMinAcross;
       nBestDown = nDown;
    }


    bool bMaximizeUsingPointAbove = 
       ( nDown > ( nSeqAcrossLength - nBandSize ) ? true : false );

    /* now handle case of other columns */
    for( nAcross = nMinAcross + 1; nAcross <= nMaxAcross; ++nAcross ) {
      nScoreDiagonallyAbove = nToBeScoreDiagonallyAbove;
      nToBeScoreDiagonallyAbove = nScore[nAcross];

      if ( nAcross != nMaxAcross ||
           bMaximizeUsingPointAbove ) {
         nScore[ nAcross ] = 
         MAX3( nScore[nAcross] - nGapPenalty, // move down
               nScore[nAcross-1] - nGapPenalty, // move right
               nScoreDiagonallyAbove + 
               nMatchMatrix[ pSeqAcross[ nAcross - 1]][ pSeqDown[ nDown - 1 ]] );
      }
      else {
         nScore[ nAcross ] =
            MAX( nScore[ nAcross-1] - nGapPenalty, // move right
                 nScoreDiagonallyAbove +
                 nMatchMatrix[ pSeqAcross[ nAcross - 1]][ pSeqDown[ nDown - 1 ]] );
      }

      if (nScore[nAcross] > nBestScore ) {
        nBestScore = nScore[nAcross];
        nBestAcross = nAcross;
        nBestDown = nDown;
      }
    } /* for( nAcross = 1 ... */


    /*    for( nAcross = 0; nAcross < nSeqAcrossLength; ++nAcross ) */
  } /* for( nDown = 1 ... */


  *pnBestAcross = nBestAcross;
  *pnBestDown = nBestDown;
  *pnBestScore = nBestScore;
  free( (char*) nScore );
}