#include "ancestral.h"

void compare(char *outFile, char *predFile);

NODE *node;
FILE *resultFile=stdout;
int seqNum, seqLength;
long idum=(-64);
char aa[21]={'-','A','R','N','D','C','Q','E','G','H','I','L','K','M','F','P','S','T','W','Y','V'};

void usage()
{
	printf("Usage: Compare_out alignment outFile cutoff\n");
	exit(1);
}
	

int main(int argc, char **argv)
{
	int i;
	int site;
	FILE *input;
	char seqFile[50],outFile[50],predFile[50];

	if(argc<4) usage();
	strcpy(seqFile,argv[1]);
	strcpy(outFile,argv[2]);
	strcpy(predFile,argv[3]);
	input=fopen(seqFile,"r");
	if(!input) {
		printf("can not open file:%s\n",seqFile);
		exit(1);
	}
	fclose(input);
	input=fopen(outFile,"r");
	if(!input) {
		printf("can not open file:%s\n",outFile);
		exit(1);
	}
	fclose(input);
	input=fopen(predFile,"r");
	if(!input) {
		printf("can not open file:%s\n",predFile);
		exit(1);
	}
	fclose(input);

	//read S_ij matrix
	readali(seqFile);
	compare(outFile,predFile);
}

void compare(char *outFile, char *predFile)
{
	int i,j,k;
	int *cutoff5,*cutoff8,index5=1,index8=1;
//	cutoff5=ivector(1,seqLength);
//	cutoff8=ivector(1,seqLength);
	cutoff5 =new int [seqLength+1];
	cutoff8 =new int [seqLength+1];
	FILE *cutoffFile;
	char line[1000];
	char seps[]=" \n";
	char *token;
	if((cutoffFile=fopen(predFile,"r"))!=NULL){
		if(fgets(line,1000,cutoffFile)!=NULL){
			token=strtok(line,seps);
			while(token!=NULL){
				if(strcmp(token,"cutoff5:")!=0){
					cutoff5[index5++]=atoi(token);
				}
				token=strtok(NULL,seps);
			}
		}
		if(fgets(line,1000,cutoffFile)!=NULL){
			token=strtok(line,seps);
			while(token!=NULL){
				if(strcmp(token,"cutoff8:")!=0){
					cutoff8[index8++]=atoi(token);
				}
				token=strtok(NULL,seps);
			}
		}
	}
	fclose(cutoffFile);
	printf("Finish reading cutoff file\n");

	bool *siteFlag,*invFlag;
	siteFlag= new bool [seqLength+1];
	invFlag= new bool [seqLength+1];
	for(i=1;i<=seqLength;i++) {
		siteFlag[i]=false; int count=0;
		for(j=1;j<=seqNum;j++){
			if(node[j].seq_aa[i]=='-') count++;
		}
		if(2*count>=seqNum) siteFlag[i]=true;
	}
	for(i=1;i<=seqLength;i++){
		invFlag[i]=true;
		for(j=2;j<=seqNum;j++) if(node[j].seq_aa[i]!=node[1].seq_aa[i]) {invFlag[i]=false;break;}
	}
	printf("Gap site:\n");
	for(i=1;i<=seqLength;i++) if(siteFlag[i]) printf("%3d ",i);
	printf("\nConserve site:\n");
	for(i=1;i<=seqLength;i++) if(invFlag[i]) printf("%3d ",i);
	printf("\n");

	int *index1,*index2,*index3,*index4;
	double *value1,*value2,*value3,*value4;
	index1=ivector(1,seqLength); index2=ivector(1,seqLength);
	index3=ivector(1,seqLength); index4=ivector(1,seqLength);
	value1=dvector(1,seqLength); value2=dvector(1,seqLength);
	value3=dvector(1,seqLength); value4=dvector(1,seqLength);
	//read out file
	FILE *out;
	out=fopen(outFile,"r");
	if(!out) {printf("can not open file: %s\n",outFile); exit(0); }
	else {
		int tmpCount=0;
		while(fgets(line,1000,out)!=NULL){
			if(strlen(line)<=1) continue;
//			printf("%d %s",strlen(line),line);
//			getchar();
			if(strlen(line)>0) {
				token=strtok(line,seps);
//				printf("%s\n",token);
//				getchar();
				if(strcmp(token,"index1")==0){
					printf("find same string\n");
					tmpCount++;
				}
				if(tmpCount==3){
					for(i=1;i<=seqLength;i++){
						fgets(line,1000,out);
						token=strtok(line,seps); index1[i]=atoi(token);
						token=strtok(NULL,seps); value1[i]=atof(token);
						token=strtok(NULL,seps); index2[i]=atoi(token);
						token=strtok(NULL,seps); value2[i]=atof(token);
						token=strtok(NULL,seps); index3[i]=atoi(token);
						token=strtok(NULL,seps); value3[i]=atof(token);
						token=strtok(NULL,seps); index4[i]=atoi(token);
						token=strtok(NULL,seps); value4[i]=atof(token);
//						printf("%s\n",line);
					}
					break;
				}
			}

		}
	}
	fclose(out);

//	for(i=1;i<=seqLength;i++) printf("%3d %10.5f %3d %10.5f %3d %10.5f %3d %10.5f\n",index1[i],value1[i],index2[i],value2[i],index3[i],value3[i],index4[i],value4[i]);

	//calculate the hit
	char hitFile[50];
	strcpy(hitFile,predFile);
	strcat(hitFile,".hit1");
	out=fopen(hitFile,"a+");
//	out=stdout;

	bool calFlag[5];
	if(value1[1]<1.0e-10) calFlag[1]=false; else calFlag[1]=true;
	if(value2[1]<1.0e-10) calFlag[2]=false; else calFlag[2]=true;
	if(value3[1]<1.0e-10) calFlag[3]=false; else calFlag[3]=true;
	if(value4[1]<1.0e-10) calFlag[4]=false; else calFlag[4]=true;

	int A5_5=0, A5_8=0,A10_5=0,A10_8=0,A15_5=0,A15_8=0,A20_5=0,A20_8=0;
	int count1=0,count2=0,count3=0,count4=0;
	for(i=seqLength;i>=0;i--){
		if(!siteFlag[index1[i]] && ((value1[i]>1.0e-10) ||(value1[i]<=1.0e-10 && invFlag[index1[i]]))){
			count1++; count2++; count3++; count4++;
			if(count1<=5) {
				for(j=1;j<index5;j++) {
					if(index1[i]==cutoff5[j]) { A5_5++; break; }
				}
				for(j=1;j<index8;j++) {
					if(index1[i]==cutoff8[j]) {A5_8++; break; }
				}
			}
			if(count2<=10) {
				for(j=1;j<index5;j++) {
					if(index1[i]==cutoff5[j]) { A10_5++; break; }
				}
				for(j=1;j<index8;j++) {
					if(index1[i]==cutoff8[j]) {A10_8++; break; }
				}
			}
			if(count2<=15) {
				for(j=1;j<index5;j++) {
					if(index1[i]==cutoff5[j]) { A15_5++; break; }
				}
				for(j=1;j<index8;j++) {
					if(index1[i]==cutoff8[j]) {A15_8++; break; }
				}
			}
			if(count2<=20) {
				for(j=1;j<index5;j++) {
					if(index1[i]==cutoff5[j]) { A20_5++; break; }
				}
				for(j=1;j<index8;j++) {
					if(index1[i]==cutoff8[j]) {A20_8++; break; }
				}
			}
		}
//		printf("%3d %3d %3d %3d %3d\n", index1[i],count1,count2,count3,count4);
//		getchar();
		if(count1>5 && count2>10 && count3>15 && count4>20) break;
	}
	if(!calFlag[1]) A5_5=A5_8=A10_5=A10_8=A15_5=A15_8=A20_5=A20_8=0;
	fprintf(out," %c  %c  %c  %c  %d  %d ",outFile[strlen(outFile)-8],outFile[strlen(outFile)-7],outFile[strlen(outFile)-6],outFile[strlen(outFile)-5],3,1);
	fprintf(out,"%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n",A5_5,A5_8,A5_8-A5_5,A10_5,A10_8,A10_8-A10_5,A15_5,A15_8,A15_8-A15_5,A20_5,A20_8,A20_8-A20_5);

	A5_5=A5_8=A10_5=A10_8=A15_5=A15_8=A20_5=A20_8=0;
	count1=count2=count3=count4=0;
	for(i=seqLength;i>=0;i--){
		if(!siteFlag[index1[i]] && ((value1[i]>1.0e-10) ||(value1[i]<=1.0e-10 && invFlag[index1[i]]))){
			count1++; count2++; count3++; count4++;
			if(count1<=5) {
				for(j=1;j<index5;j++) {
					if(index2[i]==cutoff5[j]) { A5_5++; break; }
				}
				for(j=1;j<index8;j++) {
					if(index2[i]==cutoff8[j]) {A5_8++; break; }
				}
			}
			if(count2<=10) {
				for(j=1;j<index5;j++) {
					if(index2[i]==cutoff5[j]) { A10_5++; break; }
				}
				for(j=1;j<index8;j++) {
					if(index2[i]==cutoff8[j]) {A10_8++; break; }
				}
			}
			if(count2<=15) {
				for(j=1;j<index5;j++) {
					if(index2[i]==cutoff5[j]) { A15_5++; break; }
				}
				for(j=1;j<index8;j++) {
					if(index2[i]==cutoff8[j]) {A15_8++; break; }
				}
			}
			if(count2<=20) {
				for(j=1;j<index5;j++) {
					if(index2[i]==cutoff5[j]) { A20_5++; break; }
				}
				for(j=1;j<index8;j++) {
					if(index2[i]==cutoff8[j]) {A20_8++; break; }
				}
			}
		}
		if(count1>5 && count2>10 && count3>15 && count4>20) break;
	}
	if(!calFlag[2]) A5_5=A5_8=A10_5=A10_8=A15_5=A15_8=A20_5=A20_8=0;
	fprintf(out," %c  %c  %c  %c  %d  %d ",outFile[strlen(outFile)-8],outFile[strlen(outFile)-7],outFile[strlen(outFile)-6],outFile[strlen(outFile)-5],3,2);
	fprintf(out,"%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n",A5_5,A5_8,A5_8-A5_5,A10_5,A10_8,A10_8-A10_5,A15_5,A15_8,A15_8-A15_5,A20_5,A20_8,A20_8-A20_5);
	
	A5_5=A5_8=A10_5=A10_8=A15_5=A15_8=A20_5=A20_8=0;
	count1=count2=count3=count4=0;
	for(i=seqLength;i>=0;i--){
		if(!siteFlag[index1[i]] && ((value1[i]>1.0e-10) ||(value1[i]<=1.0e-10 && invFlag[index1[i]]))){
			count1++; count2++; count3++; count4++;
			if(count1<=5) {
				for(j=1;j<index5;j++) {
					if(index3[i]==cutoff5[j]) { A5_5++; break; }
				}
				for(j=1;j<index8;j++) {
					if(index3[i]==cutoff8[j]) {A5_8++; break; }
				}
			}
			if(count2<=10) {
				for(j=1;j<index5;j++) {
					if(index3[i]==cutoff5[j]) { A10_5++; break; }
				}
				for(j=1;j<index8;j++) {
					if(index3[i]==cutoff8[j]) {A10_8++; break; }
				}
			}
			if(count2<=15) {
				for(j=1;j<index5;j++) {
					if(index3[i]==cutoff5[j]) { A15_5++; break; }
				}
				for(j=1;j<index8;j++) {
					if(index3[i]==cutoff8[j]) {A15_8++; break; }
				}
			}
			if(count2<=20) {
				for(j=1;j<index5;j++) {
					if(index3[i]==cutoff5[j]) { A20_5++; break; }
				}
				for(j=1;j<index8;j++) {
					if(index3[i]==cutoff8[j]) {A20_8++; break; }
				}
			}
		}
		if(count1>5 && count2>10 && count3>15 && count4>20) break;
	}
	if(!calFlag[3]) A5_5=A5_8=A10_5=A10_8=A15_5=A15_8=A20_5=A20_8=0;
	fprintf(out," %c  %c  %c  %c  %d  %d ",outFile[strlen(outFile)-8],outFile[strlen(outFile)-7],outFile[strlen(outFile)-6],outFile[strlen(outFile)-5],3,3);
	fprintf(out,"%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n",A5_5,A5_8,A5_8-A5_5,A10_5,A10_8,A10_8-A10_5,A15_5,A15_8,A15_8-A15_5,A20_5,A20_8,A20_8-A20_5);

	A5_5=A5_8=A10_5=A10_8=A15_5=A15_8=A20_5=A20_8=0;
	count1=count2=count3=count4=0;
	for(i=seqLength;i>=0;i--){
		if(!siteFlag[index1[i]] && ((value1[i]>1.0e-10) ||(value1[i]<=1.0e-10 && invFlag[index1[i]]))){
			count1++; count2++; count3++; count4++;
			if(count1<=5) {
				for(j=1;j<index5;j++) {
					if(index4[i]==cutoff5[j]) { A5_5++; break; }
				}
				for(j=1;j<index8;j++) {
					if(index4[i]==cutoff8[j]) {A5_8++; break; }
				}
			}
			if(count2<=10) {
				for(j=1;j<index5;j++) {
					if(index4[i]==cutoff5[j]) { A10_5++; break; }
				}
				for(j=1;j<index8;j++) {
					if(index4[i]==cutoff8[j]) {A10_8++; break; }
				}
			}
			if(count2<=15) {
				for(j=1;j<index5;j++) {
					if(index4[i]==cutoff5[j]) { A15_5++; break; }
				}
				for(j=1;j<index8;j++) {
					if(index4[i]==cutoff8[j]) {A15_8++; break; }
				}
			}
			if(count2<=20) {
				for(j=1;j<index5;j++) {
					if(index4[i]==cutoff5[j]) { A20_5++; break; }
				}
				for(j=1;j<index8;j++) {
					if(index4[i]==cutoff8[j]) {A20_8++; break; }
				}
			}
		}
		if(count1>5 && count2>10 && count3>15 && count4>20) break;
	}
	if(!calFlag[4]) A5_5=A5_8=A10_5=A10_8=A15_5=A15_8=A20_5=A20_8=0;
	fprintf(out," %c  %c  %c  %c  %d  %d ",outFile[strlen(outFile)-8],outFile[strlen(outFile)-7],outFile[strlen(outFile)-6],outFile[strlen(outFile)-5],3,4);
	fprintf(out,"%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n",A5_5,A5_8,A5_8-A5_5,A10_5,A10_8,A10_8-A10_5,A15_5,A15_8,A15_8-A15_5,A20_5,A20_8,A20_8-A20_5);
	fclose(out);
}



static char str[MAXSTR+1];

char **aname, **aseq;
int nal, alilen, *astart, *alen;


static char *mymalloc(int size)
//int size;
{
	char *buf;

	if ((buf = (char *)malloc((size_t) (size*sizeof(char)))) == NULL) {
		fprintf(stderr, "Not enough memory: %d\n", size);
		exit(1);
	}
	return buf;
}

char *ccvector(long nl, long nh)
/* allocate an char vector with subscript range v[nl..nh] */
{
	char *v;

	v=(char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(char)));
	if (!v) nrerror("allocation failure in cvector()");
	return v-nl+NR_END;
}

char *strsave(char *str)
//char *str;
{
	char *buf;
	int l;

	l = strlen(str);
	buf = mymalloc(l + 1);
//	buf=ccvector(0,l+1);
	strcpy(buf, str);
	return buf;
}

char *strnsave(char *str, int l)
//char *str;
//int l;
{
	char *buf;

	buf = mymalloc(l + 1);
//	buf=ccvector(0,l+1);
	memcpy(buf, str, l);
	buf[l] = '\0';
	return buf;
}

static char **incbuf(int n, char **was)
{
	char **buf;
	buf = (char **) malloc((size_t) (n+1) * sizeof(buf[0]));
	if(!buf) {fprintf(stderr, "Not enough memory: %d\n",n+1); exit(1); }
//	buf=cmatrix(0,(n+1),0,(n+1));
	if (n > 0) {
		memcpy(buf, was, n * sizeof(was[0]));
		free(was);
	}
	buf[n] = NULL;
	return buf;
}

static int *incibuf(int n, int *was)
//int n, *was;
{
	int *ibuf;

	ibuf = (int *) malloc((size_t) (n+1) * sizeof(int));
	if(!ibuf) {fprintf(stderr, "Not enough memory: %d\n",n+1); exit(1); }
//	ibuf=ivector(0,(n+1)*sizeof(ibuf[0]));
	if (n > 0) {
		memcpy(ibuf, was, n * sizeof(was[0]));
		free(was);
	}
	ibuf[n] = 0;
	return ibuf;
}


void readali(char *filename)
/* fills alignment from the file *filename to
nal - number of sequences in the alignment;
alilen - alignment length;
aname[0..nal-1][0..until '\0'] - names of sequences;
aseq[0..nal-1][0..alilen-1] - sequences;
astart[0..nal-1] - starting residue number;
alen[0..nal-1] - sequence length;
*/

{
	FILE *fp;
	char *s, *ss, *seqbuf;
	int i,j,k, n, l, len, len0;
	int ii,mark=1;

	printf("Reading sequences\n");
	if ((fp = fopen(filename, "r")) == NULL) {
		fprintf(stderr, "No such file: \"%s\"\n", filename);
//		err_readali(1);
		exit(1);
	}
	alilen = 0;
	nal = 0;
	n = 0;
	if(fgets(str, MAXSTR, fp) != NULL) {
		if(strncmp(str,"CLUSTAL W",9)!=0){rewind(fp);}
					   }
	while (fgets(str, MAXSTR, fp) != NULL) {
		for (ss = str; isspace(*ss); ss++) ;
//		printf("SS1:%s\n",ss);
		if ((ii<=ss-str)&&(mark==0)) {continue;}
		if (*ss == '\0') {
			if (n == 0) {
				continue;
			}
			if (nal == 0) {
				if (n == 0) {
					fprintf(stderr, "No alignments read\n");
//					err_readali(2);
					exit(1);
				}
				nal = n;
			} else if (n != nal) {
				fprintf(stderr, "Wrong nal, was: %d, now: %d\n", nal, n);
//				err_readali(3); 
				exit(1);
			}
			n = 0;
			continue;
		}
		for (s = ss; *s != '\0' && !isspace(*s); s++) ;
//		printf("S1:%s\n",s);
		*s++ = '\0';
		if (nal == 0) {
			astart = incibuf(n, astart);
			alen = incibuf(n, alen);
			aseq = incbuf(n, aseq);
			aname = incbuf(n, aname);
			aname[n] = strsave(ss);
		} else {
			if (n < 0 || n >= nal) {
				fprintf(stderr, "Bad sequence number: %d of %d\n", n, nal);
//				err_readali(4);  
				exit(1);
			}
			if (strcmp(ss, aname[n]) != 0) {
				fprintf(stderr, "Names do not match");
				fprintf(stderr, ", was: %s, now: %s\n", aname[n], ss);
//				err_readali(5);  
				exit(1);
			}
		}
		for (ss = s; isspace(*ss); ss++);
//		printf("SS2:%s\n",ss);
		if(mark==1){
		ii = ss-str;
		mark=0;}
		
		for (s = ss; isdigit(*s); s++) ;
//		printf("S2:%s\n",s);
		if (isspace(*s)) {
			if (nal == 0) {
				astart[n] = atoi(ss);
			}
			for (ss = s; isspace(*ss); ss++);
		}
		for (s = ss, l = 0; *s != '\0' && !isspace(*s); s++) {
			if (isalpha(*s)) {
				l++;
			}
		}
		len = s - ss;
//		printf("n:%d  nal:%d\n",n,nal);
		if (n == 0) {
			len0 = len;
			alilen += len;
		} else if (len != len0) {
			fprintf(stderr, "wrong len for %s", aname[n]);
			fprintf(stderr, ", was: %d, now: %d\n", len0, len);
//			err_readali(6); 
			exit(1);
		}
		alen[n] += l;
		if (aseq[n] == NULL) {
			aseq[n] = strnsave(ss, len);
		} else {
//			printf("save sequence\n");
//			seqbuf=ccvector(0,alilen);
			seqbuf=mymalloc(alilen+1);
			memcpy(seqbuf, aseq[n], alilen-len);
			free(aseq[n]);
			aseq[n] = seqbuf;
			memcpy(seqbuf+alilen-len, ss, len);
			seqbuf[alilen] = '\0';
//			printf("after save sequence\n");
		}
		n++;
	}
	if (nal == 0) {
		if (n == 0) {
			fprintf(stderr, "No alignments read\n");
//			err_readali(7);  
			exit(1);
		}
		nal = n;
	} else if (n != 0 && n != nal) {
		fprintf(stderr, "Wrong nal, was: %d, now: %d\n", nal, n);
//		err_readali(8);  
		exit(1);
	}
	fclose(fp);
//	nal=n;

//	printf("n=%d alilen=%d nal=%d\n",n,alilen,nal);
	//define new structures
	seqNum=nal; seqLength=alilen;
	int DIM= seqNum*2+1;
	node=new NODE [DIM];
	for(i=0;i<DIM;i++) { 
		node[i].seq_aa=new char [seqLength+2];
		node[i].seq_num=new int [seqLength+1];
		node[i].prob= dmatrix(1,seqLength+1,1,NP+1);
	}
	//transfer sequence to NODE structure
	for(i=1;i<=(seqNum*2-2);i++) node[i].nodeNum=i;
	for(i=1;i<=seqNum;i++) strcpy(node[i].nodeName,aname[i-1]);
	for(i=0;i<DIM;i++) { node[i].seq_aa[0]=' ';node[i].seq_aa[seqLength+1]='\0';}
//	for(i=0;i<seqNum;i++) printf("%s\n",aseq[i]);
	for(i=1;i<=seqNum;i++) for(j=1;j<=seqLength;j++) node[i].seq_aa[j]=aseq[i-1][j-1];	
	//print the sequence info
	fprintf(resultFile,"The input sequences:\n");
	//assign seq_num
	for(i=1;i<=seqNum;i++){
		for(j=1;j<=seqLength;j++){
			bool unknownChar=true;
			for(k=0;k<=NP;k++){
				if(node[i].seq_aa[j]==aa[k]) { unknownChar=false; node[i].seq_num[j]=k; break; }
			}
			if(unknownChar) { node[i].seq_num[j]=0; node[i].seq_aa[j]='-'; }
		}
	}
	//compare sequence name to find duplicates
	printf("try to find the duplicate name of different sequences\n");
	bool dupFlag=false;
	for(i=1;i<seqNum;i++){
		for(j=i+1;j<=seqNum;j++){
			if(strcmp(node[i].nodeName, node[j].nodeName)==0){ 
				dupFlag=true;
				printf("duplication: seqNo:%d seqName:%s  seqNo:%d seqName:%s\n",i,node[i].nodeName,j,node[j].nodeName);
			}
		}
	}
	if(dupFlag) exit(1);
	

}


char *itoa(int b)
{
	char *convert;
	int count=0,copy,c,index=0,d;
	copy= b;
	while(copy>0){
		count++;
		copy/=10;
	}
	convert= new char [count+1];
	for(int i=count;i>=0;i--)
	{	d=1;
		for(int j=0;j<i;j++) d=d*10;
		c=b/d;
		b=b-c*d;
		if(c>0)
		{
			if(c==1) convert[index]='1';
			else if(c==2) convert[index]='2';
			else if(c==3) convert[index]='3';
			else if(c==4) convert[index]='4';
			else if(c==5) convert[index]='5';
			else if(c==6) convert[index]='6';
			else if(c==7) convert[index]='7';
			else if(c==8) convert[index]='8';
			else if(c==9) convert[index]='9';
			index++;
		}
		else if(c==0)
		{
			if(index>0) 
			{
				convert[index]='0';
				index++;
			}
		}
	}
	convert[index]='\0';
	return convert;
}	
