#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <unistd.h>
#include <math.h>

#include "FASTA.h"
#include "kmer.h"
#include "needle.h"
#include "util.h"
#include "global.h"

#define MAXPATH 2000
int KmerLen=6;

double gap_o=10;
double gap_e=0.5;
int SeqCount=1000;
double step=0.03;
double pval=0.05;

char **SeqStrs;
int *SeqLens;
int seqmax=1000;

typedef struct{
	double kdist;
	double ndist;
	double krate;
}KrateRec;

KrateRec *kdrec;

void LoadSeqs(char *seqfile)
{
	FASTA fasta;
	char *seq;
	char *label;
	int idx=0;
	int seqlen;
	int i;
	
	
	SeqLens=(int *)Malloc(SeqCount*sizeof(int));
	SeqStrs=(char **)Malloc(SeqCount*sizeof(char *));

	fasta.OpenRead(seqfile);
	while ((seqlen=fasta.GetFastaSeq(&label,&seq,0))>0)
	{
			if (2*seqlen >seqmax)
				seqmax=2*seqlen;
			if (seq[0]!=0)
			{
				SeqLens[idx]=seqlen;
				SeqStrs[idx++]=seq;
			}
			else
			{
				seq[0]=1;
				free(seq);
			}
			free(label);
			if (idx>=SeqCount) break;		
	}
	fasta.Close();
	if (SeqCount > idx) 
		{
			for (i=idx;i<SeqCount;i++)
				free(SeqStrs[i]);
			SeqCount = idx;
		}
	printf("%d Sequences Read\n",SeqCount);
}

void KmerCorr()
{
	int **SeqTabs;
	int **SeqCodes;
	Needle *needle;
	Kmer *kmer;
	int i,j;
	int *kmertab;
	int *kmerseq;
	double kdist,ndist;
	char *str1,*str2;
	kmer=new Kmer(KmerLen);
	needle= new Needle(gap_o,gap_e);
		
	kdrec=(KrateRec *)Malloc(SeqCount*(SeqCount-1)/2*sizeof(KrateRec));
		
	SeqTabs=(int**)Malloc(SeqCount*sizeof(int *));
	SeqCodes=(int**)Malloc(SeqCount*sizeof(int *));
	for (i=0;i<SeqCount;i++)
	{
		kmerseq=(int *)Malloc((SeqLens[i]-KmerLen+1)*sizeof(int));
		kmertab=kmer->AllocCodeTable();
		kmer->KmerCount(SeqStrs[i], kmertab,kmerseq);
		SeqTabs[i]=kmertab;		
		SeqCodes[i]=kmerseq;
	}
	
	str1=(char *)Malloc(seqmax*sizeof(char));
	str2=(char *)Malloc(seqmax*sizeof(char));
	
	int ptr=0;
	for (i=0;i<SeqCount-1;i++)
	{
		printf("\rComputing %d of %d",i+1,SeqCount);
		fflush(stdout);
		for (j=i+1;j<SeqCount;j++)
		{
			if (SeqLens[i]>SeqLens[j])
				kdist =1.0 - kmer->KmerComp(SeqTabs[j],SeqTabs[i],SeqCodes[j],SeqLens[j]);
			else
				kdist =1.0 - kmer->KmerComp(SeqTabs[i],SeqTabs[j],SeqCodes[i],SeqLens[i]);
			needle->Align(SeqStrs[i],SeqStrs[j],str1, str2, 1.0);
			ndist =	calcDistance(str1,str2);
			kdrec[ptr].kdist=kdist;
			kdrec[ptr].ndist=ndist;
			kdrec[ptr++].krate=kdist/ndist;
		}
	}
	printf("\n");		
	free(str1);
	free(str2);
	for (i=0;i<SeqCount;i++)
	{
		free(SeqTabs[i]);
		free(SeqCodes[i]);
		free(SeqStrs[i]);
	}
	
	free(SeqTabs);
	free(SeqCodes);
	free(SeqLens);
	free(SeqStrs);
	delete kmer;
	delete needle;
}

int CompNDist(const void * a, const void * b)
{
	KrateRec *rec1=(KrateRec *)a;
	KrateRec *rec2=(KrateRec *)b;
	if (rec1->ndist < rec2->ndist) return -1;	
	if (rec1->ndist > rec2->ndist) return 1;
	return 0;	
}

int CompRate(const void * a, const void * b)
{
	KrateRec *rec1=(KrateRec *)a;
	KrateRec *rec2=(KrateRec *)b;
	if (rec1->krate < rec2->krate) return -1;
	if (rec1->krate > rec2->krate) return 1;
	return 0;
}

void WriteKrate(char * ofile)
{
	FILE *fp=fopen(ofile,"w");
	if (fp==NULL)
	{
		fprintf(stderr,"Cannot Create Output\n");
		exit(1);
	}
	int numrec= SeqCount*(SeqCount-1)/2;
	
	qsort(kdrec,numrec,sizeof(KrateRec),CompNDist);
	
	float maxd=kdrec[numrec-1].ndist;
	int numit=(int) ceil(maxd/step);
	fprintf(fp,"%d\n",numit);
	fprintf(fp,"%.2f\n",step);
	
	KrateRec *headptr;
	int ptr=0;
	for (int i=0;i<numit;i++)
	{
		double cutoff=step*(i+1);
		headptr=(kdrec+ptr);
		int numpt=0;
		while (ptr < numrec && kdrec[ptr].ndist <cutoff) 
		{
			ptr++;
			numpt++;
		}
		qsort(headptr,numpt,sizeof(KrateRec),CompRate);
		int numcut=(int) ceil(numpt* (1-pval));
		if (numcut > numpt) numcut=numpt;
		if (numcut <1) numcut=1;
		fprintf(fp,"%.6f\n",headptr[numcut-1].krate);
	}
	fclose(fp);
	free(kdrec);
}

void Usage()
{
	printf("FindKrate : Generate Kmer Configure File for given data and parameter.\n");
	printf("FindKrate [-c seqnum] [-k klen] [-g gap_open] [-e gap_extend] [-s step] <input> <output>\n");
	printf("\t seqnum:\t Number of seqs to calculate. Default 1000.\n");
	printf("\t -k\t:\t Kmer length.\n");
	printf("\t -g,e\t:\t gap_open/gap_extend parameter for NW align\n");
	printf("\t -s\t:\t stepsize of distance levels. [0.01-0.1] Default 0.03.\n");
}

int main(int argc, char **argv)
{
	char *ifile, *ofile;
	char c;
	clock_t time_1,time_2;
	
	ifile=(char *)Malloc(sizeof(char)*MAXPATH);
	ofile=(char *)Malloc(sizeof(char)*MAXPATH);
	
  while ((c = getopt(argc, argv, "hc:k:g:e:s:p:")) != -1)
    {
        switch (c)
        {
        case 'h':
        	Usage();
        	exit(0);
        	break;
        case 'c':
        	SeqCount = atoi(optarg);
        	break;
		case 'k':
        	KmerLen = atoi(optarg);
			break;
		case 'g':
			gap_o = atof(optarg);
			break;
		case 'e':
			gap_e = atof(optarg);
			break;		
		case 's':
			step = atof(optarg);
			break;
		case 'p':
			pval = atof(optarg);
			break;
        default:
        	fprintf(stderr,"Unknown option -%c\n",c);
					break;
        }
    }
    argc -= optind;
    argv += optind;
	if (argc<2) 
		{
			Usage();
			exit(0);
		}
		
	strcpy(ifile,argv[0]);
	strcpy(ofile,argv[1]);
	
	time_1=clock();
	LoadSeqs(ifile);
	KmerCorr();
	WriteKrate(ofile);
	time_2=clock();
	#ifdef CONSOLE
	printf("\rDist Correlation Computed in %10.3g secs\n",(double)(time_2-time_1)/CLOCKS_PER_SEC);
	#endif
	free(ifile);
	free(ofile);
	return 0;		
}
