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

  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 Gotoh
//
int main (int argc, char *argv[]) {
  if (argc != 6) 
    {
      cout << "Usage: " << argv[0] << " A B mismatch_cost indel_var indel_const" << endl; 
      exit(1);
    }


  // ---------------------------------------------------------  
  // a, b
  //
  Sequence a = argv[1];
  Sequence b = argv[2];
  cost_t subst = atof(argv[3]);
  cost_t indel_var = atof(argv[4]);
  cost_t indel_const = atof(argv[5]);
  cout << " mismatch: " << subst 
       << " indel-const: " << indel_const 
       << " indel-var: " << indel_var
       << endl;
  // ---------------------------------------------------------
  // auxiliary variables
  //
  int n = a.length();
  int m = b.length();
  cost_t d1, p1, p2, q1, q2;

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

  // ---------------------------------------------------------
  // the different arrows
  //
  // define arrows
  //                              
  SubstArrow diag_arrow(TraceD,1,subst);
  AffineIndelArrow leftD_arrow(TraceD,0,1,indel_const+1*indel_var);
  AffineIndelArrow leftQ_arrow(TraceQ,0,1,indel_var);
  AffineIndelArrow upD_arrow(TraceD,1,0,indel_const+1*indel_var);
  AffineIndelArrow upP_arrow(TraceP,1,0,indel_var);
  TraceArrow P_arrow(TraceP,0,0);
  TraceArrow Q_arrow(TraceQ,0,0);

  // ---------------------------------------------------------
  // Intialization
  //
  D[0][0]=0;
  P[0][0]=INFTY;
  Q[0][0]=INFTY;
  for (int i=1;i<=n;++i)
    {
      Q[i][0]=INFTY;
      D[i][0] = indel_const + indel_var*i;
      TraceD[i][0].push_back( TMatrixEntryElem(upD_arrow) );
    }
  for (int j=1;j<=m;++j)
    {
      P[0][j]=INFTY;
      D[0][j]= indel_const + indel_var*j;
      TraceD[0][j].push_back( TMatrixEntryElem(leftD_arrow) );
    }
  
  // ---------------------------------------------------------
  // recursive definition of D-Matrix
  //
  for (int i=1;i<=n;++i)
    for (int j=1;j<=m;++j)
      {
        // calculate P_ij
        //
        //  for debuging
        // cout << "hier i: " << i << "  j:" << j << endl;
        p1= P.get(upP_arrow.end_position(i,j))
          + upP_arrow.cost(i,j);
        p2= D.get(upD_arrow.end_position(i,j))
          + upD_arrow.cost(i,j);
        P[i][j]=min(p1,p2);
        // set trace entries
        //
        if (p1 == P[i][j])
          TraceP[i][j].push_back(  TMatrixEntryElem(upP_arrow) );
        if (p2 == P[i][j])
          TraceP[i][j].push_back(  TMatrixEntryElem(upD_arrow) );

        // calculate Q_ij
        //
        q1= Q.get(leftQ_arrow.end_position(i,j))
          + leftQ_arrow.cost(i,j);
        q2= D.get(leftD_arrow.end_position(i,j))
          + leftD_arrow.cost(i,j);
        Q[i][j]=min(q1,q2);
        // set trace entries
        //
        if (q1 == Q[i][j])
          TraceQ[i][j].push_back(  TMatrixEntryElem(leftQ_arrow) );
        if (q2 == Q[i][j])
          TraceQ[i][j].push_back(  TMatrixEntryElem(leftD_arrow) );


        // calculate D_ij
        //
        d1= D.get(diag_arrow.end_position(i,j))
          + diag_arrow.cost(i,j,a,b);
        D[i][j] = min(d1,min(P[i][j],Q[i][j]));
        

        // set trace entries
        //
        if (d1 == D[i][j])
          TraceD[i][j].push_back(  TMatrixEntryElem(diag_arrow) );
        if (P[i][j] == D[i][j])
          TraceD[i][j].push_back(  TMatrixEntryElem(P_arrow) );
        if (Q[i][j] == D[i][j])
          TraceD[i][j].push_back(  TMatrixEntryElem(Q_arrow) );
      }
    
  // ------------------------------------------------------------
  // do backtrace
  BackTracer bt;

  // ---------------------------------------------------------
  // calc traceback
  //
  bt.calc_traceback(TraceD,n,m,0);   // 0 mean stop at cell (0,0);
  
  // ---------------------------------------------------------
  // 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();

  // print out alignment
  //
  bt.print_alignment(a,b);
  
  // print out the other matrices
  //
  cout <<  endl;
  cout << "P-matrix:" << endl << P << endl;
  cout << "Q-matrix:" << endl << Q << endl;

  // ---------------------------------------------------------
  // done
  exit(0);
}