#include <string.h>
#include <time.h>
#include <unistd.h>
#include <math.h>
#include "FASTA.h"
#include "needle.h"
#include "kmer.h"
#include "util.h"
#include "ProbModel.h"
#include <fstream>
#include <iostream>
using namespace std;

int SeqCount=10000;
int ComputeAll=0;
int trim_space=1;

char **SeqStrs;
int SeqNum;
int *SeqLens;

REAL gap_o=10;
REAL gap_e=0.5;

REAL calcDistance(char *seq1, char *seq2) {


  unsigned int i, j, k;
  REAL residuecount, distance;
	
  residuecount = distance = 0.0;

	int gap_alert1 = 0;
	int gap_alert2 = 0;

	int start_seq1 = 0;
  int start_seq2 = 0;
  int end_seq1 = strlen(seq1);
  int end_seq2 = strlen(seq2);

	k = 0;
  while(seq1[k] == '-') k++;
  start_seq1 =  k;

  k = 0;
  while(seq2[k] == '-') k++;
  start_seq2 = k;

	k = end_seq1-1;
  while(seq1[k] == '-') k--;
  end_seq1 = k+1; 
  	
	k = end_seq2-1;
  while(seq2[k] == '-') k--;
  end_seq2 = k+1; 
	
  int start_seq = start_seq1;
  if(start_seq < start_seq2)
     {
        start_seq = start_seq2;
     }

  int end_seq = end_seq1;
  if(end_seq > end_seq2)
     {
         end_seq = end_seq2;
     }
	gap_alert1 = 0;
	gap_alert2 = 0;

	for( k=start_seq; k <end_seq; k++) {

	   if (seq1[k] == '-' && seq2[k] == '-')
		{
			continue;
		}

  	if(seq1[k] == 'N' && seq2[k] == 'N')
		{
			gap_alert1 = 0;
			gap_alert2 = 0;
			continue;
		}


  	if( seq1[k] == '.' &&seq2[k] == '.')
		{
			gap_alert1 = 0;
			gap_alert2 = 0;
			continue;
		}


/* multiple consecutive gaps  are counted as 1 mismatch */
/*  but residue is shortened. */
/* you have to check both sequences... */


    if ( seq1[k] == '-' )
		{
			if(gap_alert1 == 0)
			{
        	residuecount+=1.0;
         	distance+=1.0;
			}
			gap_alert1 = 1;
      gap_alert2 = 0;
			continue;
		}

    if ( seq2[k] == '-' )
		{
			if(gap_alert2 == 0)
			{
        	 	residuecount+=1.0;
           	distance+=1.0;
			}
			gap_alert1 = 0;
      gap_alert2 = 1;
			continue;
		}

     /* ok no gaps , check for mismatch */

    if ( seq1[k] !=seq2[k])
		{
    	 distance += 1.0;
       residuecount+=1.0;
			 gap_alert1 = 0;
			 gap_alert2 = 0;
	 		continue;
	 	}

   if ( seq1[k] ==seq2[k])
		{
      residuecount+= 1.0;
			gap_alert1 = 0;
			gap_alert2 = 0;
			continue;	
		}
  } /* end of loop on k */ 

  if (residuecount > 0) {
		distance = distance / residuecount;
  }
  else {
		distance = 1.0;
  }
  return distance;  
}


void LoadSeqs(char *seqfile)
{
	FASTA fasta;
	char *seq;
	char *label;
	int idx=0;
	int seqlen;
	
	SeqStrs=(char **)Malloc(SeqCount*sizeof(char *));
	SeqLens=(int *)Malloc(SeqCount*sizeof(int));
	fasta.OpenRead(seqfile);
	while ((seqlen=fasta.GetFastaSeq(&label,&seq,0))>0)
	{
			SeqLens[idx]=seqlen;
			SeqStrs[idx++]=seq;
			free(label);		
	}
	fasta.Close();
	SeqNum=idx;
}

void FreeSeqs()
{
	int i;
	for (i=0;i<SeqNum;i++)
		free(SeqStrs[i]);
	free(SeqStrs);
	free(SeqLens);
}

void CalcDist(char *ofile,char *distfile)
{
	Needle needle(gap_o,gap_e,1);

	Kmer kmer(6,4);
	clock_t time_1,time_2;
	
	ProbString al1,al2;
	char *als1,*als2;
	int id1,id2;
	REAL dist;
	
	als1=(char *)Malloc(1000*sizeof(char));
	als2=(char *)Malloc(1000*sizeof(char));

	time_1=clock();

	int **KmerTabs=(int**)Malloc(SeqNum*sizeof(int *));
	int **KmerSeqs=(int**)Malloc(SeqNum*sizeof(int *));
	for (int i=0;i<SeqNum;i++)
	{
		KmerSeqs[i]=(int *)Malloc((SeqLens[i]-5+1)*sizeof(int));
		KmerTabs[i]=kmer.AllocCodeTable();
		kmer.KmerCount(SeqStrs[i],KmerTabs[i],KmerSeqs[i]);
	}
	
	for (id1=0;id1 < SeqNum-1;id1++)
	{
		for (id2=id1+1;id2<SeqNum;id2++)
		{
			 if (SeqLens[id1]> SeqLens[id2])
				 dist=1 - kmer.KmerComp(KmerTabs[id2],KmerTabs[id1],KmerSeqs[id2],SeqLens[id2]);
			else	 
				 dist=1 - kmer.KmerComp(KmerTabs[id1],KmerTabs[id2],KmerSeqs[id1],SeqLens[id1]);
		}
	}
	time_2=clock();
	printf("Time Int Kmer %10.3g\n",(REAL) (time_2-time_1)/CLOCKS_PER_SEC);
	for (int i=0;i<SeqNum;i++)
	{
		free(KmerSeqs[i]);
		free(KmerTabs[i]);
	}	
	free(KmerSeqs);
	free(KmerTabs);

	ProbString *Strs=new ProbString[SeqNum];
	for (id1=0;id1<SeqNum;id1++)
		Strs[id1].FromString(SeqStrs[id1]);

	time_1=clock();

	REAL **KcodeTable=(REAL **)Malloc(SeqNum*sizeof(REAL *));
	for (id1=0;id1 < SeqNum;id1++)
	{
		KcodeTable[id1]=kmer.KmerCount(Strs[id1]);
	}
	
	for (id1=0;id1 < SeqNum-1;id1++)
		for (id2=id1+1;id2<SeqNum;id2++)
		{
			dist=1-kmer.KmerComp( KcodeTable[id1], KcodeTable[id2]);
		}
	time_2=clock();	
	printf("Time REAL Kmer %10.3g\n",(REAL) (time_2-time_1)/CLOCKS_PER_SEC);
	for (id1=0;id1<SeqNum;id1++) kmer.FreeCodeTable(KcodeTable[id1]);
	free(KcodeTable);	
	
	time_1=clock();
		
	for (id1=0;id1 < SeqNum-1;id1++)
	{
		for (id2=id1+1;id2<SeqNum;id2++)
		{
			needle.Align(SeqStrs[id1],SeqStrs[id2],als1,als2,0.25);
			dist=calcDistance(als1,als2);
		}
	}
	time_2=clock();
	printf("Time Str Needle %10.3g\n",(REAL) (time_2-time_1)/CLOCKS_PER_SEC);
	time_1=clock();
		
	for (id1=0;id1 < SeqNum-1;id1++)
	{
		for (id2=id1+1;id2<SeqNum;id2++)
		{
			needle.Align(SeqStrs[id1],Strs[id2],als1,al2,0.25);
			dist=al2.AveDist(als1);
		}
	}
	time_2=clock();
	printf("Time Str-Pr Needle %10.3g\n",(REAL) (time_2-time_1)/CLOCKS_PER_SEC);
	time_1=clock();

	for (id1=0;id1 < SeqNum-1;id1++)
	{
		for (id2=id1+1;id2<SeqNum;id2++)
		{
			needle.Align(Strs[id1],Strs[id2],al1,al2,0.25);
			dist=al2.AveDist(al1);
		}
	}
	time_2=clock();
	printf("Time Pr Needle %10.3g\n",(REAL) (time_2-time_1)/CLOCKS_PER_SEC);
	
	delete [] Strs;
		
  free(als1);
  free(als2);
}


int main(int argc, char **argv)
{
	char *seqfile, *ofile, *distfile;
	
	seqfile=argv[1];
	ofile=argv[2];
	distfile=argv[3];
	
	LoadSeqs(seqfile);
	CalcDist(ofile,distfile);
	FreeSeqs();	
	return 0;		
}
