#include "ancestral.h"

void reproduce(char *seqFile,char *headFile,int last, int first,int at);

NODE *node;
FILE *resultFile=stdout;
int seqNum, seqLength;
char aa[NP+1]={'-','A','R','N','D','C','Q','E','G','H','I','L','K','M','F','P','S','T','W','Y','V'};
long idum=(-64);
void usage()
{
//	printf("Usage: chooseSeq -i sequencefile -l listfile -n sequenceNumber\n");
	printf("Usage: chooseSeq -i sequencefile -L (1~7) -C (1~4) -o output alignment\n");
	exit(1);
}


int main(int argc, char **argv)
{
	int i;
	FILE *input;
	int last=0, first=0,at=0;
	char seqFile[50], headFile[50];

	if(argc<=1) usage();
//	for(i=1;i<argc;i++) printf("%s\n",argv[i]);
	for(i=1;i<argc;++i){
//		printf("i=%d\n",i);
		char option,optidx;
		if(argv[i][0]=='-') option=argv[i][1];
		else usage();
//		printf("option:%c\n",option);
		switch(option){
			case 'L':
				last=atoi(argv[++i]);
				break;
			case 'F':
				first=atoi(argv[++i]);
				break;
			case 'A':
				at=atoi(argv[++i]);
				break;

			case 'i':
				++i;
				optidx=0;
				if(i<argc) {
					strcpy(seqFile,argv[i]);
//					printf("%s\n",seqFile);
					input=fopen(seqFile,"r");
					if(!input) {
						printf("can not open file:%s\n",seqFile);
						exit(1);
					}
					fclose(input);
				}
				else usage();
				break;
			case 'H':
				++i;
				optidx=0;
				if(i<argc) {
					strcpy(headFile,argv[i]);
//					printf("%s\n",seqFile);
					input=fopen(headFile,"r");
					if(!input) {
						printf("can not open file:%s\n",headFile);
						exit(1);
					}
					fclose(input);
				}
				else usage();
				break;
			default:
				usage();
		}
	}

	//read S_ij matrix
	readali(seqFile);
	reproduce(seqFile,headFile,last,first,at);
}


void reproduce(char *seqFile,char *headFile,int last, int first,int at)
{
	int i,j,k;
	FILE *input;
	input=fopen(headFile,"r");
	char name[50],line[3000];
	fgets(name,50,input);
	fgets(line,3000,input);
	fclose(input);
	//deal with name
//	printf("deal with name\n");
	strcpy(name,"pdb_");
	strcat(name,headFile);
	
	for(i=0;i<=strlen(name);i++) if(name[i]=='.') {name[i]='\0'; break;}

	//prepare output
	for(i=0;i<=strlen(seqFile);i++) if(seqFile[i]=='.') {seqFile[i]='\0'; break;}
	strcat(seqFile,"_pdb.aln");
	input=fopen(seqFile,"w");
	fprintf(input,"%-20s %s",name,line);
//	printf("output first sequence\n");
	//deal with $
	for(i=1;i<=seqNum;i++) {
		for(j=0;j<=strlen(node[i].nodeName);j++) if(node[i].nodeName[j]=='$') node[i].nodeName[j]='_';
	}
//	printf("deal with $ \n");
	for(i=1;i<=seqNum;i++){
		fprintf(input,"%-20s ",node[i].nodeName);
		for(j=1;j<=seqLength;j++){
			if(j>first && j<=(seqLength-last) && j!=at) fprintf(input,"%c",node[i].seq_aa[j]);
		}
		fprintf(input,"\n");
	}
	fclose(input);
}
	

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-1)*2+1;
//	node=new NODE [DIM];
	int DIM=seqNum+1;
	node=new NODE [seqNum+1];
	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<=seqNum/*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];
	for(i=1;i<=seqNum;i++) for(j=1;j<=seqLength;j++) node[i].seq_aa[j]=toupper(node[i].seq_aa[j]);
	//deal with non-amino acid character
	for(i=1;i<=seqNum;i++){
		for(j=1;j<=seqLength;j++){
			bool flag=false;
			for(int k=0;k<=NP;k++) if(node[i].seq_aa[j]==aa[k]) flag=true;
			if(!flag) 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);

	//put template sequence into node[0]
	for(i=0;i<=seqNum-1;i++) node[i]=node[i+1];
	seqNum--;
//	printf("seqNum:%d\n",seqNum);
	//find the template sequence non gap length
	int stdLength=0;
	for(i=1;i<=seqLength;i++) if(node[0].seq_aa[i]!='-') stdLength++;
	//find those sequences that can not cover 80% of template sequence non-gap position
	int *lengthDist; lengthDist=ivector(1,seqNum);
	for(i=1;i<=seqNum;i++){
		lengthDist[i]=0;
		for(j=1;j<=seqLength;j++){
			if(node[0].seq_aa[j]!='-'){
				if(node[i].seq_aa[j]!='-') lengthDist[i]++;
			}
		}
		lengthDist[i]=lengthDist[i]*100/stdLength;
	}
	//histogram
	int *binCount;
	binCount=ivector(1,20);
	for(i=1;i<=20;i++) binCount[i]=0;
	for(i=1;i<=seqNum;i++){
			for(k=0;k<=19;k++){
				if((lengthDist[i]>=k*5)&&(lengthDist[i]<(k+1)*5)) {
					binCount[k+1]++;
					break;
				}
				else if(lengthDist[i]>95) {
					binCount[20]++;
					break;
				}
			}

	}
//	for(i=1;i<=20;i++) printf("%5d\n",binCount[i]);
	//delete sequence
	for(i=1;i<=seqNum;i++){
		if(lengthDist[i]<75){
/*			for(j=i;j<seqNum;j++) { node[j]=node[j+1]; lengthDist[j]=lengthDist[j+1]; }
			i--;
			seqNum--;
*/			while(lengthDist[seqNum]<75) { seqNum--;}
			node[i]=node[seqNum];
			seqNum--;

		}
	}
//	printf("seqNum:%d\n",seqNum);
	//print the sequence info
//	fprintf(resultFile,"The input sequences:\n");
	for(i=1;i<=seqNum;i++){
//		fprintf(resultFile,"%5d %20s %s",i,node[i].nodeName, node[i].seq_aa);
//		fprintf(resultFile,"%20s %s",node[i].nodeName, node[i].seq_aa);
//		fprintf(resultFile,"\n");
	}

}


