/*****************************************************************************
#   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    <stdlib.h>
#include    <stdio.h>
#include "extendAlignment.h"

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

#include    "numutil.h"
#include    <iostream.h>
#include    "assert.h"
#include    <string.h>



void extendAlignment( char* pSeqDown, char* pSeqAcross, int nSeqDownLength, 
	int nSeqAcrossLength,  int* pnBestDown, int* pnBestAcross, 
	int* pnBestScore,
    const int nGap,
    const int nMatch,
    const int nMismatch ) {

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

  assert( strlen( pSeqDown ) == nSeqDownLength );
  assert( strlen( pSeqAcross ) == nSeqAcrossLength );



  // 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

  for ( nAcross = 0; nAcross <= nSeqAcrossLength; ++nAcross ) {
    nScore[ nAcross ] = -nGap * nAcross;
  }


  for( nDown = 1; nDown <= nSeqDownLength; ++nDown ) {
    nToBeScoreDiagonallyAbove = nScore[0];


    /* case of first column */
    nScore[0] = nScore[0] - nGap;

    /* now handle case of other columns */
    for( nAcross = 1; nAcross <= nSeqAcrossLength; ++nAcross ) {
      nScoreDiagonallyAbove = nToBeScoreDiagonallyAbove;
      nToBeScoreDiagonallyAbove = nScore[nAcross];
      
         // 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 - 1] == pSeqDown[nDown-1] ) 
         nScoreByDiagonal = nScoreDiagonallyAbove + nMatch;
      else
         nScoreByDiagonal = nScoreDiagonallyAbove - nMismatch;

      nScore[nAcross] = MAX( nScore[nAcross] - nGap,
			     MAX( nScore[nAcross - 1] - nGap,
				  nScoreByDiagonal
				)
                           );


      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 );
}