#include "ancestral.h"
#define TIMES 10

void calConsensus(bool prototype);
void calGapPercentage();

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: calConsensus alignment [Y]\n");
	printf("Y: the prototype sequence is the first sequence in alignment\n");
	exit(1);
}
	

int main(int argc, char **argv)
{
	int i;
	int site;
	FILE *input;
	char seqFile[50];
	bool prototype=false;

	if(argc<2) usage();
	strcpy(seqFile,argv[1]);
	if(argc==3) {
		if(argv[2][0]=='Y') prototype=true;
	}

	input=fopen(seqFile,"r");
	if(!input) {
		printf("can not open file:%s\n",seqFile);
		exit(1);
	}
	fclose(input);

	//read S_ij matrix
	readali(seqFile);
	calGapPercentage();
//	calConsensus(prototype);
}

void calGapPercentage()
{
	int i,j,k;
	double *gapPercentage;
	gapPercentage=dvector(1,seqLength);
	for(i=1;i<=seqLength;i++) gapPercentage[i]=0.0;
	for(i=1;i<=seqLength;i++) {
		for(j=1;j<=seqNum;j++) {
			bool gapFlag=true;
			for(k=1;k<=20;k++) {
				if(node[j].seq_aa[i]==aa[k]) { gapFlag=false; break; }
			}
			if(!gapFlag) gapPercentage[i]+=1.0;
		}
		gapPercentage[i]/=(1.0*seqNum);
		gapPercentage[i]=1.0-gapPercentage[i];
	}
	for(i=1;i<=seqLength;i++) printf("%-4d %6.4f\n",i,gapPercentage[i]);
}

void calConsensus(bool prototype)
{
	int i,j,k;
	double **prob;
	prob=dmatrix(1,seqLength,1,20);
	int *equal;
	equal=ivector(1,seqLength); for(i=1;i<=seqLength;i++) equal[i]=0;
	for(i=1;i<=seqLength;i++) for(j=1;j<=20;j++) prob[i][j]=0.0;
	double *max, *max_prototype;
	max=dvector(1,seqLength); max_prototype=dvector(1,seqLength);
	if(prototype) {
		//find the consensus sequence and store in node[0]
		for(i=1;i<=seqLength;i++) for(j=2;j<=seqNum;j++) prob[i][node[j].seq_num[i]]+=1.0/(seqNum-1);
		for(i=1;i<=seqLength;i++) {
			double tmpProb=prob[i][1]; int tmpIndex=1;
			for(j=2;j<=20;j++) {
				if(prob[i][j]>tmpProb) {tmpProb=prob[i][j]; tmpIndex=j; }
			}
			node[0].seq_num[i]=tmpIndex;
			node[0].seq_aa[i]=aa[tmpIndex];
			max[i]=tmpProb;
			max_prototype[i]=prob[i][node[1].seq_num[i]];
			for(j=1;j<=20;j++) {
				if(fabs(prob[i][j]-tmpProb)<1.0e-8 && tmpIndex!=j) equal[i]=1;
			}
		}
		double *identity;
		identity=dvector(2,seqNum);
		for(i=2;i<=seqNum;i++) {
			for(j=1;j<=seqLength;j++) {
				if(node[i].seq_num[j]==node[1].seq_num[j]) identity[i]+=1.0/seqLength;
			}
		}
		double identity_consensus=0.0;
		for(j=1;j<=seqLength;j++) {
			if(node[0].seq_num[j]==node[1].seq_num[j]) identity_consensus+=1.0/seqLength;
		}
		double identity_consensus1=0.0;
		for(j=1;j<=seqLength;j++) {
			if(fabs(max[j]-max_prototype[j])<1.0e-8){
				if(equal[j]!=1) identity_consensus1+=1.0;
				else {
					double equal_aa=0.0;
					for(i=1;i<=NP;i++) {
						if(fabs(prob[j][i]-max[j])<1.0e-8) equal_aa+=1.0;
					}
					identity_consensus1+=1.0/equal_aa;
				}
			}
		}
		identity_consensus1/=seqLength;
		//output result
		fprintf(resultFile,"Consensus sequence:\n");
		fprintf(resultFile,"%-20s %s\n","consensus",node[0].seq_aa);
		fprintf(resultFile,"%-20s %s\n","prototype",node[1].seq_aa);
		fprintf(resultFile,"probability for each position\n");
		for(i=1;i<=NP;i++) fprintf(resultFile, "   %c ",aa[i]);
		fprintf(resultFile,"  max   max_prototype equal?\n");
		for(i=1;i<=seqLength;i++){
			for(j=1;j<=NP;j++) fprintf(resultFile,"%3.2f ",prob[i][j]);
			fprintf(resultFile," %3.2f  %3.2f ",max[i],max_prototype[i]);
			if(equal[i]==1) fprintf(resultFile,"   *\n");
			else fprintf(resultFile,"\n");
		}
		fprintf(resultFile,"Identity between consensus and prototype sequence: %6.2f\n",identity_consensus);
		fprintf(resultFile,"Identity between consensus and prototype sequence <after consider tie prob>:%6.2f\n",identity_consensus1);
		fprintf(resultFile,"Identity between reconstructed sequences and prototype sequence:\n");
		for(i=2;i<=seqNum;i++) fprintf(resultFile,"%3.2f\n",identity[i]);
	}
	else {
                //find the consensus sequence and store in node[0]
                for(i=1;i<=seqLength;i++) for(j=1;j<=seqNum;j++) prob[i][node[j].seq_num[i]]+=1.0/seqNum;
                for(i=1;i<=seqLength;i++) {
                        double tmpProb=prob[i][1]; int tmpIndex=1;
                        for(j=2;j<=20;j++) {
                                if(prob[i][j]>tmpProb) {tmpProb=prob[i][j]; tmpIndex=j; }
                        }
                        node[0].seq_num[i]=tmpIndex;
                        node[0].seq_aa[i]=aa[tmpIndex];
                        max[i]=tmpProb;
                        for(j=1;j<=20;j++) {
                                if(fabs(prob[i][j]-tmpProb)<1.0e-8 && tmpIndex!=j) equal[i]=1;
                        }

                }
                //output result
                fprintf(resultFile,"Consensus sequence:\n");
                fprintf(resultFile,"%-20s %s\n","consensus",node[0].seq_aa);
                for(i=1;i<=NP;i++) fprintf(resultFile, "   %c ",aa[i]);
                fprintf(resultFile,"  max equal?\n");
		for(i=1;i<=seqLength;i++){
                        for(j=1;j<=NP;j++) fprintf(resultFile,"%3.2f ",prob[i][j]);
			fprintf(resultFile," %3.2f ",max[i]);
                        if(equal[i]==1) fprintf(resultFile,"   *\n");
			else fprintf(resultFile,"\n");
                }
	}
}



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;
		}
	}
	for(i=1;i<=seqNum;i++) {
		if(strlen(node[i].nodeName)>16) {
			node[i].nodeName[0]='s'; node[i].nodeName[1]='e';node[i].nodeName[2]='q';
		}
		else {
			char tmpName[50];
			strcpy(tmpName,"seq");
			strcat(tmpName,node[i].nodeName);
			if(strlen(tmpName)>20) tmpName[19]='\0';
			strcpy(node[i].nodeName,tmpName);
		}
	}
	//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);
*/
	int composition[21];	
	for(i=0;i<=20;i++) composition[i]=0;
	for(i=1;i<=seqNum;i++) for(j=1;j<=seqLength;j++) composition[node[i].seq_num[j]]++;
	for(i=1;i<=20;i++) printf("%6d ",composition[i]);
	printf("\n");
}

