#include "variance/Aligner.hh"

#include <algorithm>
#include <iostream>
#include <vector>

/**
 * Project     : CASAVA
 * Module      : $RCSfile: Aligner.cpp,v $
 * @author     : Tony Cox
 * Copyright   : Copyright (c) Illumina 2008, 2009. All rights reserved.
 *
 ** This software is covered by the "Illumina Genome Analyzer Software
 ** License Agreement" and the "Illumina Source Code License Agreement",
 ** and certain third party copyright/licenses, and any user of this
 ** source file is bound by the terms therein (see accompanying files
 ** Illumina_Genome_Analyzer_Software_License_Agreement.pdf and
 ** Illumina_Source_Code_License_Agreement.pdf and third party
 ** copyright/license notices).
 *
 */

namespace ca {
namespace variance_detection {

typedef unsigned char NScoreType;

typedef std::vector<std::vector<ScoreType> > DPMatrix;
typedef std::vector<std::vector<NScoreType> > DPMatrixN;



static
void max3(ScoreType& max,
          NScoreType& which,
          const ScoreType v0,
          const ScoreType v1,
          const ScoreType v2) {
    max=v0;
    which=0;
    if (v1>v0) {
        max=v1;
        which=1;
    }
    if (v2>max) {
        max=v2;
        which=2;
    }
    //    cout << v0 << " " << v1 << " " << v2 << " - " << max << " " << which << endl;
} // ~max3



ScoreType Aligner::operator()
( const char* x, const char* y, const unsigned x_size, const unsigned y_size )
{
    //  const ScoreType w_match(3),w_mismatch(1),w_open(2),w_extend(2);
    xt_="";
    yt_="";
    at_.clear();

    //std::cout << x_size << std::endl;
    //std::cout << y_size << std::endl;
    //std::cout << x << std::endl;
    //std::cout << y << std::endl;

    const unsigned ys(y_size+1);
    const unsigned xs(x_size+1);

    DPMatrix E_(xs,std::vector<ScoreType>(ys,0));
    DPMatrix F_(xs,std::vector<ScoreType>(ys,0));
    DPMatrix G_(xs,std::vector<ScoreType>(ys,0));

    DPMatrixN TE_(xs,std::vector<NScoreType>(ys,0));
    DPMatrixN TF_(xs,std::vector<NScoreType>(ys,0));
    DPMatrixN TG_(xs,std::vector<NScoreType>(ys,0));

	for (unsigned i(0);i<=x_size;i++)
	{
		//E_[i][0]=-w_open_-(i*w_extend_);
	        //E_[i][0]=0;
	        E_[i][0]=-10000;
		F_[i][0]=-10000;
		G_[i][0]=-10000;
	}

	for (unsigned j(0);j<=y_size;j++)
	{
	  // Chris's suggestion
	  //E_[0][j]=-w_open_;
	  //F_[0][j]=-10000;

	  //F_[0][j]=-w_open_;
	  F_[0][j]=0;
	  E_[0][j]=-10000;

	  //F_[0][j]=-10000;
	  //F_[0][j]=0;
	  //F_[0][j]=-10000;
	  //E_[0][j]=0;

	  //E_[0][j]=-10000;
	  //F_[0][j]=-w_open_-(j*w_extend_);
	  //F_[0][j]=0;
	  G_[0][j]=0;
	}
	for (unsigned i(1);i<=x_size;i++)
	{
#ifdef DEBUG
		cout << i << ": ";
#endif
		for (unsigned j(1);j<=y_size;j++)
		{
			// update G_
		  //max3(G_[i][j],TG_[i][j],
		  //G_[i-1][j-1],E_[i-1][j-1],F_[i-1][j-1]);
			max3(G_[i][j],TG_[i][j],
					G_[i-1][j-1],E_[i-1][j-1] - w_open_, F_[i-1][j-1] - w_open_);

			//      G_[i][j]=V_[i-1][j-1];
			if (x[i-1]==y[j-1])
			{
				G_[i][j]+=w_match_;
			}
			else
			{
				G_[i][j]-=w_mismatch_;
			}
			//      max =G_[i][j];
			//      dir=0;

			// update E_
			max3(E_[i][j],TE_[i][j],
					G_[i][j-1]-w_open_,
					E_[i][j-1],
					F_[i][j-1]-w_open_);

			E_[i][j]-=w_extend_;
			//      E_[i][j]=V_[i][j-1]-w_open_;
			//      if (E_[i][j-1]>E_[i][j])
			//	E_[i][j]=E_[i][j-1];
			//   E_[i][j]-=w_extend_;

			//      if (E_[i][j]>max)
			//      {
			//	max=E_[i][j];
			//	dir=1;
			//      }

			// update F_
			max3(F_[i][j],TF_[i][j],
					G_[i-1][j]-w_open_,
					E_[i-1][j]-w_open_,
					F_[i-1][j]);

			F_[i][j]-=w_extend_;
			//      F_[i][j]=V_[i-1][j]-w_open_;
			//    if (F_[i-1][j]>F_[i][j])
			//	F_[i][j]=F_[i-1][j];
			//   F_[i][j]-=w_extend_;

			//      if (F_[i][j]>max)
			//  {
			//	max=F_[i][j];
			//	dir=2;
			//     }
			//      V_[i][j]=max;
			//    T_[i][j]=dir;
#ifdef DEBUG
			//      cout << i << " " << j;
			printf(" %3d:%3d:%3d/%1d%1d%1d",G_[i][j],E_[i][j],F_[i][j],TG_[i][j],TE_[i][j],TF_[i][j]);
			//      cout << '\t' << G_[i][j]<< ":"<< E_[i][j] << ":" << F_[i][j];
			//   cout << "/";
			//  cout << TG_[i][j]<< ":"<< TE_[i][j] << ":" << TF_[i][j];
			//      cout << endl;
			//      cout << " " << max << "/" << dir;
#endif
		} // ~for j
#ifdef DEBUG
		cout << endl;
#endif
	} // ~for i

	//  cout << TG_[6][7] << " " << TE_[6][7] << " " << TF_[6][7] << endl;
        ScoreType max(0),thisMax;
        NScoreType whichMatrix(0),thisMatrix;
        int ii(0),jj(0);

#if 0
	for (unsigned i(0);i<=x_size;i++)
	{
		max3( thisMax, thisMatrix,
				G_[i][y_size], E_[i][y_size], F_[i][y_size] );
		if ((i==0) or (thisMax>max))
		{
			max=thisMax;
			ii=i;
			jj=y_size;
			whichMatrix=thisMatrix;
		} // ~if

		//    if (V_[i][y_size]>max)
		//   {
		//ii=i;
		// jj=y_size;
		//  max=V_[i][y_size];
		//  }
	} // ~for i
#endif
	for (unsigned j(0);j<=y_size;j++)
	{
	  //max3( thisMax, thisMatrix,
	  //		G_[x_size][j], E_[x_size][j], F_[x_size][j] );

	  if (G_[x_size][j] > F_[x_size][j]) {
	    thisMax=G_[x_size][j];
	    thisMatrix=0;
	  } else {
	    thisMax=F_[x_size][j];
	    thisMatrix=2;
	  }

	  //max3( thisMax, thisMatrix,
	  //G_[x_size][j], -10000, F_[x_size][j] );
	  if ((j==0) or (thisMax>max))
	    {
	      max=thisMax;
	      ii=x_size;
	      jj=j;
	      whichMatrix=thisMatrix;
	    } // ~if

		//    if (V_[x_size][j]>max)
		//  {
		//    ii=x_size;
		//   jj=j;
		//   max=V_[x_size][j];
		//  }
	}


	score_=max;
#ifdef DEBUG
	cout << "start cell = " << ii << " " << jj << " " << whichMatrix << endl;
#endif
	x_start_=ii;
	y_start_=jj;

	DPMatrixN* pT_ =&TE_;

	while ((ii>0)&&(jj>0))
	{
		if (whichMatrix==0)
		{
			pT_=&TG_;
		}
		else if (whichMatrix==1)
		{
			pT_=&TE_;
		}
		else
		{
			//assert(whichMatrix==2);
			pT_=&TF_;
		}
		//  cout << TG_[6][7] << " " << TE_[6][7] << " " << TF_[6][7] << endl;

#ifdef DEBUG
		cout << ii << " " << jj << " " << whichMatrix << " ";
		cout << " "<< G_[ii][jj]<< ":"<< E_[ii][jj] << ":" << F_[ii][jj];
		cout << "/";
		cout << " "<< TG_[ii][jj]<< ":"<< TE_[ii][jj] << ":" << TF_[ii][jj];
		cout << " - " << (*pT_)[ii][jj] << endl;
#endif
		const NScoreType nextMatrix=(*pT_)[ii][jj];
		if (whichMatrix==0)
		{
			xt_+=x[ii-1];
			yt_+=y[jj-1];
			at_+=((x[ii-1]==y[jj-1])?'|':'.');
			ii--; jj--;
			// whichMatrix=0;
		} // ~if
		else if (whichMatrix==1)
		{
			xt_+='-';
			yt_+=y[jj-1];
			at_+=' ';
			jj--;
			//    whichMatrix=1;
		} // ~else if
		else
		{
			//assert(whichMatrix==2);
			xt_+=x[ii-1];
			yt_+='-';
			at_+=' ';
			ii--;
			//      whichMatrix=2;
		} // ~else
		whichMatrix=nextMatrix;
	} // ~while
#ifdef DEBUG
	cout << "end cell = " << ii << " " << jj << endl;
#endif
	x_end_=ii;
	y_end_=jj;

	reverse(xt_.begin(),xt_.end());
	reverse(at_.begin(),at_.end());
	reverse(yt_.begin(),yt_.end());
	//cout << xt_ <<endl;
	// cout << at_ <<endl;
	//  cout << yt_ <<endl;
	return 0;
} // ~ScoreType Aligner::operator()( const string& x, const string& y)


}
} // end namespace casava{ namespace { applications