/*------------------------------------------------------------
  Copyright (C) 1999 by Rolf Backofen, Sebastian Will.
  All Rights Reserved.
  
  Permission to use, copy, modify, and distribute this
  software and its documentation for NON-COMMERCIAL purposes
  and without fee is hereby granted provided that this
  copyright notice appears in all copies.
  
  THE AUTHORS AND PUBLISHER MAKE NO REPRESENTATIONS OR
  WARRANTIES ABOUT THE SUITABILITY OF THE SOFTWARE, EITHER
  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
  IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
  PARTICULAR PURPOSE, OR NON-INFRINGEMENT. THE AUTHORS
  AND PUBLISHER SHALL NOT BE LIABLE FOR ANY DAMAGES SUFFERED
  BY LICENSEE AS A RESULT OF USING, MODIFYING OR DISTRIBUTING
  THIS SOFTWARE OR ITS DERIVATIVES.

  ------------------------------------------------------------*/


#include <stdlib.h>
#include <algo.h>
#include <math.h>
#include <stdio.h>
#include <list.h>

#include "matrix.hh"
#include "matrices.hh"
#include "trace.hh"
#include "tracearrow.hh"
#include "printarray.hh"


// ---------------------------------------------------------
// Implementation of Needleman/Wunsch local
//
int main (int argc, char *argv[]) {
  if (argc != 6) 
    {
      cout << "Usage: " << argv[0] << " A B match_cost mismatch_cost indel_cost" << endl;
      exit(1);
    }

  // ---------------------------------------------------------  
  // a, b
  //
  Sequence a = argv[1];
  Sequence b = argv[2];
  cost_t match = atof(argv[3]);
  cost_t mismatch = atof(argv[4]);
  cost_t indel = atof(argv[5]);
  
  cout << "indel: " << indel << " match: " << match << " mismatch: " << mismatch << endl;

  // ---------------------------------------------------------
  // auxiliary variables
  //
  int n = a.length();
  int m = b.length();
  cost_t v1, v2, v3;

  // ---------------------------------------------------------
  // define the distance matrix and trace matrix
  //
  DistanceMatrix D(a,b);
  TraceMatrix TraceD(n,m,D);

  // ---------------------------------------------------------
  // the different arrows
  //
  // define arrows
  //                              
  SimilarityArrow diag_arrow(TraceD,1,match,mismatch);
  AffineIndelArrow left_arrow(TraceD,0,1,indel);
  AffineIndelArrow up_arrow(TraceD,1,0,indel);

  // ---------------------------------------------------------
  // Intialization
  //
  D[0][0]=0;
  for (int i=1;i<=n;++i)
    {
      D[i][0] = 0;
    }
  for (int j=1;j<=m;++j)
    {
      D[0][j] = 0;
    }
  
  
  cost_t maxD=MINUSINFTY;
  int    max_i=0, max_j=0;
  
  // ---------------------------------------------------------
  // recursive definition of D-Matrix
  //
  for (int i=1;i<=n;++i)
    for (int j=1;j<=m;++j)
      {
        // calculate D_ij
        //
        //  for debuging
        // cout << "hier i: " << i << "  j:" << j << endl;
        v1= D.get(diag_arrow.end_position(i,j))
          + diag_arrow.cost(i,j,a,b);
        v2= D.get(up_arrow.end_position(i,j))
          + up_arrow.cost(i,j);
        v3= D.get(left_arrow.end_position(i,j))
          + left_arrow.cost(i,j);
        D[i][j] = max(0.0,max(max(v1,v2),v3));

        // set trace entries
        //
        if ( approx_equal( v1, D[i][j] ) )
          TraceD[i][j].push_back(  TMatrixEntryElem(diag_arrow) );
        if ( approx_equal( v2, D[i][j] ) )
          TraceD[i][j].push_back(  TMatrixEntryElem(up_arrow) );
        if ( approx_equal( v3, D[i][j] ) )
          TraceD[i][j].push_back(  TMatrixEntryElem(left_arrow) );
	
	// calculate maximal element.
	//
        if (D[i][j] > maxD)
          {
            max_i = i;
            max_j = j;
            maxD = D[i][j];
          }
      }
  
  // ------------------------------------------------------------
  // do backtrace
  BackTracer bt;

  // ---------------------------------------------------------
  // calc traceback
  //
  bt.calc_traceback(TraceD,max_i,max_j,1);   // 1 for end at 0 value
  
  // ---------------------------------------------------------
  // print D and traceback
  // we have to calculate how many characters are printed maximally 
  // in lines and columns
  // 
  PrintArray pr(1+2*(n+1),2+(m+1)*(ROWFILL+PRECISION));

  // first, the trace is drawn on D (connecting all cells in
  // backtrace with
  //
  bt.draw_traceback(pr);
  // then draw the matrix overwriteing part of the trace
  //
  D.draw(pr);

  // now print the distance matrix and trace stored in pr
  //
  pr.print();

  // finally, print out alignment
  //
  bt.print_alignment(a,b);
  
  // ---------------------------------------------------------
  // done
  exit(0);
}