#include "ancestral.h"

void shaffle(char *seqFile);
void random(char *seqFile);

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: shaffle alignment\n");
	exit(1);
}
	

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

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

	//read S_ij matrix
	readali(seqFile);
	shaffle(seqFile);
	random(seqFile);
}

void shaffle(char *seqFile)
{
	printf("In shaffle function\n");
	double ran1(long *idum);
	int i,j,k,m;
	char shaffleFile[50];
	strcpy(shaffleFile,seqFile);
	FILE *output;
	for(i=1;i<=seqNum;i++) strcpy(node[i+seqNum].seq_aa,node[i].seq_aa);
	for(i=seqNum+1;i<=seqNum*2;i++) sprintf(node[i].nodeName,"%d",i);
	for(i=1;i<=10000;i++) for(j=1;j<=8000;j++) double y=0.1*i*j;	
	srand(time(NULL));
	idum=-(rand()%1000);
	printf("idum=%d\n",idum);
	for(i=1;i<=seqLength;i++){
		for(j=seqNum+1;j<=seqNum*2-1;j++){
			int shaffleNum=(int) (ran1(&idum)*(seqNum*2-j+1)+j);
			char shaffleChar=node[shaffleNum].seq_aa[i];
			if(shaffleChar=='-' || node[j].seq_aa[i]=='-') continue;
			int tmpK=0;
			for(k=shaffleNum;k>=j+1;k--){ //node[k].seq_aa[i]=node[k-1].seq_aa[i];
				if(node[k].seq_aa[i]=='-') continue;
				tmpK=0;
				for(m=k-1; m>=j; m--) if(node[m].seq_aa[i]!='-') {tmpK=m; break;}
				if(tmpK!=0) node[k].seq_aa[i]=node[tmpK].seq_aa[i];
			}
			node[j].seq_aa[i]=shaffleChar;
		}
	}
	for(i=0;i<strlen(shaffleFile);i++) if(shaffleFile[i]=='.' ||shaffleFile[i]=='_') {shaffleFile[i]='\0'; break;}
	strcat(shaffleFile,"_3.aln");
	output=fopen(shaffleFile,"w");
	for(i=seqNum+1;i<=seqNum*2;i++) fprintf(output,"%-20s %s\n",node[i].nodeName,node[i].seq_aa);
	fclose(output);
}

void random(char *seqFile)
{
	int i,j,k;
	char randomFile[50];
	FILE *output;
	double ran1(long *idum);
	strcpy(randomFile,seqFile);
	double PI[21];
	int aaCount[21];
	for(i=0;i<=20;i++) aaCount[i]=0;
	//calculate PI
	for(i=1;i<=seqLength;i++){
		int count=0;
		for(j=1;j<=seqNum;j++)	if(node[j].seq_num[i]==0) count++;
		if(count<=seqNum-3) for(j=1;j<=seqNum;j++) aaCount[node[j].seq_num[i]]++;
	}
	int total=0;
	for(i=0;i<=20;i++) total+=aaCount[i];
	for(i=0;i<=20;i++) PI[i]=1.0*aaCount[i]/total;
	double accumulate[21];
	accumulate[0]=PI[0];
	for(i=1;i<=20;i++) accumulate[i]=PI[i]+accumulate[i-1];
//	for(i=0;i<=20;i++) printf("%8.4f   %8.4f\n",PI[i],accumulate[i]);
	srand(time(NULL));
	idum=-(rand()%1000);
	printf("in random function, idum = %d\n",idum);
	for(i=seqNum+1;i<=seqNum*2;i++){
		for(j=1;j<=seqLength;j++){
			double randNum=ran1(&idum);
			if(randNum<=accumulate[0]) {node[i].seq_aa[j]='-';node[i].seq_num[j]=0;}
			else {
				for(k=1;k<=20;k++) {
					if(randNum>accumulate[k-1] && randNum<=accumulate[k]) {
						node[i].seq_num[j]=k;
						node[i].seq_aa[j]=aa[k];
						break;
					}
				}
			}
//			printf("%8.4f  %d  %c ",randNum,node[i].seq_num[j],node[i].seq_aa[j]);
//			getchar();
		}
	}
	for(i=0;i<strlen(randomFile);i++) if(randomFile[i]=='-' || randomFile[i]=='.') {randomFile[i]='\0'; break;}
	strcat(randomFile,"_4.aln");
	output=fopen(randomFile,"w");
	for(i=seqNum+1;i<=seqNum*2;i++) fprintf(output,"%-20s %s\n",node[i].nodeName,node[i].seq_aa);
	fclose(output);

}




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);
	

}


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;
}	
