#include "ancestral.h"

//Global Variables
NODE *node;
double S_ij[NP+1][NP+1];
double Q_ij[NP+1][NP+1];
double PI[NP+1],**sitePI;
int seqLength,seqNum,DIM;
int *treeRoot;
bool *dashFlag,*calFlag;
char aa[NP+1]={'-','A','R','N','D','C','Q','E','G','H','I','L','K','M','F','P','S','T','W','Y','V'};
//used for calculate the eignvalue and eignvector of matrix Q
double **a, **v, **inv, *wr, *wi, *w, **svdv, *colu;
//used for print the tree
int **levelFlag;
//for output file
FILE *resultFile=stdout;
topology **gapTree;


void usage()
{
	printf("Usage: ancestral [-O/-C] [-D] -i sequencefile [-t treefile] [-o outputfile]\n");
	printf("       DEFAULT: use constant alfa (1.0) and Goldstein method to find the optimal solution (lower precision)\n");
	printf("       -O: optimize alfa (more time and higher precision)\n");
	printf("       -C: calculate correlated alfa (less time and high precision)\n");
	printf("       -D:  User Dan method to find the optimal solution\n");
	printf("       -i name: read sequences from file\n");
	printf("       -t name: read tree structure from file or use weighbor method to produce tree file\n");
	printf("       -o name: write the result into file or standard output device\n");
	exit(1);
}
	

int main(int argc, char **argv)
{
	int i;
	int refineFlag=1, conditionalProb=0;
	bool flag=true,treefileFlag=false;
	FILE *input; int number;
	char seqFile[20], treeFile[20];
	unsigned int calTime1, calTime2;
	strcpy(treeFile,"DISTMATRIX.tre");
	//initiate the matrix and vectors
	a=dmatrix(1,NP,1,NP);//copy matrix Q will be destroyed
	v=dmatrix(1,NP,1,NP);//store eignvector matrix
	inv=dmatrix(1,NP,1,NP);//store inverse of eignvector matrix
	wr=dvector(1,NP);//store eignvalue
	wi=dvector(1,NP);//store eignvalue
	w=dvector(1,NP);//used for SVD algo
	svdv=dmatrix(1,NP,1,NP);//used for SVD algo
	colu=dvector(1,NP);//used for inverse of eignvector matrix

	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 'O':
				refineFlag=0;
				break;
			case 'C':
				refineFlag=2;
				break;
			case 'D':
				conditionalProb=1;
				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 't':
				++i;
				optidx=0;
				if(i<argc) {
					treefileFlag=true;
					strcpy(treeFile,argv[i]);
				}
				else usage();
				break;
			case 'o':
				++i;
				optidx=0;
				if(i<argc) {
					resultFile=fopen(argv[i],"w");
					if(!resultFile) {
						printf("can not open file:%s\n",seqFile);
						exit(1);
					}
				}
				else usage();
				break;
			default:
				usage();
		}
	}

	//read S_ij matrix
	calTime1=time(NULL);
	init_S();
	readali(seqFile);
	searchGap();
	cal_PI();
	cal_Q();//calculate PI vector and Q(ij) matrix
	cal_Eign();//calculate eignvalue vector and eignvector matrix of Q(ij)
	if(!treefileFlag) cal_dist();
	else readTree(treeFile);
	//assign tree structure to gapTree
	assignTree();
	//for all that have at least 3 nongap leave nodes, rebuild the tree
/*	printf("Begin rebuilding tree\n");
	for(i=1;i<=seqLength;i++){
		if(dashFlag[i]&&calFlag[i]) rebuildTree(i);
	}
	printf("Finish rebuilding tree\n");
*/
	calTime2=time(NULL);
	printf("Prepare time: %ds\n", calTime2-calTime1);
	calTime1=calTime2;
	if((refineFlag==0)||(refineFlag==2)) optimalAlfa(refineFlag,conditionalProb);
	else noOptimalAlfa(conditionalProb);
				//also calculate the inverse matrix of eignvector matrix
	calTime2=time(NULL);
	printf("Time to search optimal result: %ds\n", calTime2-calTime1);
	printf("Print tree structure...\n");
	fprintf(resultFile,"\n\n\nThe tree structure:\n\n");
	print_tree(seqNum+1);
	printf("Finish print tree structure\n");
	fclose(resultFile);

}



void assignTree()
{
	int i,j;
	for(i=1;i<=(seqNum*2-2);i++){
		for(j=1;j<=seqLength;j++){
			gapTree[j][i].nodeLevel=node[i].nodeLevel;
			gapTree[j][i].childL=node[i].childL;
			gapTree[j][i].childR=node[i].childR;
			gapTree[j][i].childE=node[i].childE;
			gapTree[j][i].parent=node[i].parent;
			gapTree[j][i].rootFlag=node[i].rootFlag;
			gapTree[j][i].dist=node[i].dist;
			if(i<=seqNum){
				if(node[i].seq_aa[j]=='-') gapTree[j][i].gapFlag=true;
				else gapTree[j][i].gapFlag=false;
			}
			else gapTree[j][i].gapFlag=false;
		}
	}
	for(i=1;i<=seqLength;i++){
		if(!dashFlag[i]) treeRoot[i]=seqNum+1;
		else if(calFlag[i]) treeRoot[i]=rebuildTree(i);
		else treeRoot[i]=0;
	}
}

void readTree(char treeFile[])
{
	int i,j,k,index=seqNum+1;
	char c;
	FILE *input;
	//for tree
	int *treeL,*treeR,treeE,max,left,right;
	char tmpName[20];
	int tmpIndex;
	float tmpValue;
	int levelIndex=-1,levelP;

	printf("Reading treefile: %s\n",treeFile);
	//read *.tre file
	treeL=ivector(1,(seqNum-1)*2+1);
	treeR=ivector(1,(seqNum-1)*2+1);
	for(i=1;i<=(seqNum-1)*2;i++) { treeL[i]=0; treeR[i]=0; node[i].rootFlag=false; }
//	if((input=fopen("pii2.tre","r"))==NULL)
	if((input=fopen(treeFile,"r"))==NULL)
		nrerror("Data file pii2.tre not found\n");
	i=index; max=i;left=1;right=0;
	while(fgetc(input)!='(') {}
	node[i].nodeNum=i;
	while((c=fgetc(input))!=EOF) {
//		printf("read the letter: %c",c);
//		getchar();
		if(c!=' ' && c!=',' && c!='\n') {
			if(c=='(') {
				left=left+1;
				if(treeL[i]==0) { treeL[i]=i+1; i=i+1; max=max+1; }
				else {
					if(treeR[i]==0) { treeR[i]=max+1;max=max+1;i=max; }
					else { treeE=max+1;max=max+1;i=max; }
				}
			}
			else {//1
				if(c==')') {
					right=right+1;
					node[i].nodeNum=i;
					if((c=fgetc(input))=='\n') { }
					else ungetc(c,input);
					if((c=fgetc(input))==':'){ 
//						printf("read dist for node: %d\n",i);
						fscanf(input,"%f",&tmpValue); node[i].dist=tmpValue;
//						printf("read dist for node: %d and dist: %f\n",i,node[i].dist);
						for(j=index;j<=(index-1)*2-2;j++) {
							if(treeL[j]==i || treeR[j]==i) {i=j; break;}
						}
						
					}
				}
				else {//2
//					printf("read node name\n");
					ungetc(c,input);
					if(treeL[i]==0) {
						tmpIndex=0;
						while((c=fgetc(input))!=':') { tmpName[tmpIndex]=c; tmpIndex++;}
						for(j=1;j<=index;j++) {
							if(memcmp(node[j].nodeName,tmpName,tmpIndex)==0) {
//								printf("Yes\n");
								treeL[i]=node[j].nodeNum;
								fscanf(input,"%f",&tmpValue); 
								node[j].dist=tmpValue;
//								printf("node: %d dist: %f\n",j,node[j].dist);
//								printf("in tree Left: i=%d, j=%d, nodeNum=%d, dist=%f\n",i,j,node[j].nodeNum,node[j].dist);
								break;							}
//							else printf("No\n");
						}
					}
					else {//3
						if(treeR[i]==0) {
							tmpIndex=0;
							while((c=fgetc(input))!=':') { tmpName[tmpIndex]=c; tmpIndex++;}
//							printf("i=%d\n",i);
							for(j=1;j<=index;j++) {
								if(memcmp(node[j].nodeName,tmpName,tmpIndex)==0) {
									treeR[i]=node[j].nodeNum;
									fscanf(input,"%f",&tmpValue); node[j].dist=tmpValue;
//									printf("in tree Right: i=%d, j=%d, nodeNum=%d, dist=%f\n",i,j,node[j].nodeNum,node[j].dist);
									break;
								}
							}
						}
						else {//4
							tmpIndex=0;
							while((c=fgetc(input))!=':') { tmpName[tmpIndex]=c; tmpIndex++;}
							for(j=1;j<=index;j++) {
								if(memcmp(node[j].nodeName,tmpName,tmpIndex)==0) {
									treeE=node[j].nodeNum;
									fscanf(input,"%f",&tmpValue); node[j].dist=tmpValue;
//									printf("in tree E: i=%d, j=%d, nodeNum=%d, dist=%f\n",i,j,node[j].nodeNum,node[j].dist);
									break;
								}
							}
						}//end else 4
					}//end else 3
				}//end else 2
			}//end else 1
		}//end if
	}//end while
	fclose(input);
//	printf("Finish reading tree file\n");
/*	for(i=index;i<=(index-1)*2-2;i++){
		printf("i: %d  treeLeft: %d  treeRight: %d\n",i,treeL[i],treeR[i]);
	}
	printf("treeE: %d\n",treeE);
	for(i=1;i<=(index-1)*2-2;i++){
		printf("node: %d, nodeName: %10s distance: %f\n", node[i].nodeNum, node[i].nodeName, node[i].dist);
	}
*/
	//input child, parent, level information into node[]
	//root node
	node[index].rootFlag=true;
	node[index].childL=treeL[index];node[index].childR=treeR[index];node[index].childE=treeE;
	node[index].parent=node[index].nodeNum;
	node[treeL[index]].parent=index;node[treeR[index]].parent=index;node[treeE].parent=index;
	//other nodes
	for(i=index+1;i<=(index-1)*2-2;i++){
		node[i].childL=treeL[i]; node[i].childR=treeR[i];
		node[treeL[i]].parent=i; node[treeR[i]].parent=i;
	}
//	printf("Finish assigning node, before assign level\n");
	setLevel(index);
/*
	for(i=1;i<=(index-1)*2-2;i++){
		levelIndex=0;
		levelP=i;
		while(node[levelP].parent!=index){
			levelIndex++;
			levelP=node[levelP].parent;
		}
		if(i!=index) levelIndex++;
		node[i].nodeLevel=levelIndex;
	}

	for(i=1;i<=(index-1)*2-2;i++){
		printf("nodeNum:%2d nodeName:%10s level:%d childL:%2d childR:%2d",node[i].nodeNum,node[i].nodeName,node[i].nodeLevel,node[i].childL,node[i].childR);
		if(i!=index)
			printf(" parent:%2d dist:%f\n",node[i].parent,node[i].dist);
		else 
			printf(" childE:%2d parent:%2d dist:%f\n",node[i].childE,node[i].parent,node[i].dist);
	}
*/
	printf("Finish reading tree\n");
}

void setLevel(int nodeNum)
{
	if(node[nodeNum].rootFlag){
		node[nodeNum].nodeLevel=0;
		setLevel(node[nodeNum].childL);
		setLevel(node[nodeNum].childR);
		setLevel(node[nodeNum].childE);
	}
	else if((node[nodeNum].childL!=0)&&(node[nodeNum].childR!=0)){
		node[nodeNum].nodeLevel=node[node[nodeNum].parent].nodeLevel+1;
		setLevel(node[nodeNum].childL);
		setLevel(node[nodeNum].childR);
	}
	else node[nodeNum].nodeLevel=node[node[nodeNum].parent].nodeLevel+1;
}		

void cal_dist()
{
	FILE *fp;
	char cmd[80];
	int i,j,k;
	double **dist_matrix;
	int **compare;
	dist_matrix=dmatrix(1,seqNum,1,seqNum);
	compare=imatrix(1,seqNum,1,seqNum);
	for(i=2;i<=seqNum;i++){
		for(j=1;j<i;j++){
			compare[i][j]=1;
			for(k=1;k<=seqLength;k++) {
				if(node[i].seq_aa[k]!=node[j].seq_aa[k]) {
					compare[i][j]=0;break;
				}
			}
		}
	}
    	fp = fopen("DISTMATRIX.txt", "w");
	printf("Calculating distance matrix and produce tree file\n");
	for(i=1;i<=seqNum;i++) for(j=1;j<=seqNum;j++) dist_matrix[i][j]=0.0;
	for(i=2;i<=seqNum;i++){
		for(j=1;j<i;j++){
			if(compare[i][j]==1) dist_matrix[i][j]=0.0;
			else dist_matrix[i][j]=dist(i,j);
			if(dist_matrix[i][j]>10.0)
				printf("dist_matrix[%d][%d]=%8.5f\n",i,j,dist_matrix[i][j]);
		}
	}

	for(i=1;i<=seqNum;i++) 	dist_matrix[i][i]=0.0;
	for(i=1;i<seqNum;i++) for(j=i+1;j<=seqNum;j++) dist_matrix[i][j]=dist_matrix[j][i];
	//print the distance matrix
	fprintf(resultFile,"\n\nThe distance matrix:\n");
	for(i=1;i<=seqNum;i++){
		for(j=1;j<=seqNum;j++) fprintf(resultFile,"%10.5f ",dist_matrix[i][j]);
		fprintf(resultFile,"\n");
	}

	fprintf(fp,"%d\n",seqNum);
	for(i=1;i<=seqNum;i++){
		fprintf(fp,"%20s\t",node[i].nodeName);
		for(j=1;j<=seqNum;j++) fprintf(fp,"%10.5f ",dist_matrix[i][j]);
		fprintf(fp,"\n");
	}
	fclose(fp);

	//use weighbor method to build a tree
	printf("Use weighbor software to build tree\n");
/*	strcpy(cmd,"./weighbor_1.2.1_linux ");
//	strcpy(cmd,"./weighbor ");
	strcat(cmd,"-L ");
	strcat(cmd, itoa(seqLength));
	strcat(cmd, " -b 20 ");
	strcat(cmd, "-i DISTMATRIX.txt -o DISTMATRIX.tre");
	system(cmd);
*/
	weighbor(seqNum,seqLength,node,dist_matrix);
	printf("Finish building tree\n");
	free_imatrix(compare,1,seqNum,1,seqNum);
	free_dmatrix(dist_matrix,1,seqNum,1,seqNum);
	setLevel(seqNum+1);

/*
	for(i=1;i<=seqNum*2-2;i++){
		printf("nodeNum:%2d nodeName:%10s level:%d childL:%2d childR:%2d",node[i].nodeNum,node[i].nodeName,node[i].nodeLevel,node[i].childL,node[i].childR);
		if(i!=(seqNum+1))
			printf(" parent:%2d dist:%f\n",node[i].parent,node[i].dist);
		else 
			printf(" childE:%2d parent:%2d dist:%f\n",node[i].childE,node[i].parent,node[i].dist);
	}
*/
}

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;
}	
	
	
	
double dist(int node1, int node2)
{
	int i,j,JMAX=40;
	int **Nab;
	double xLeft,xRight,dx,xMid,xacc=0.000001;
	double leftBracket,rightBracket,midBracket;
	Nab=imatrix(0,NP,0,NP);
	for(i=0;i<=NP;i++) for(j=0;j<=NP;j++) Nab[i][j]=0;
	xLeft=0.000001;
	for(i=1;i<=seqLength;i++) Nab[node[node1].seq_num[i]][node[node2].seq_num[i]]++;
	for(i=2;i<=NP;i++) for(j=1;j<i;j++) Nab[i][j]=Nab[i][j]+Nab[j][i];
//	for(i=1;i<NP;i++) for(j=i+1;j<=NP;j++) Nab[i][j]=Nab[j][i];
/*
	printf("Nab matrix\n");
	for(i=1;i<=NP;i++){
		for(j=1;j<=NP;j++) printf("%4d",Nab[i][j]);
		printf("\n");
	}
	for(i=0;i<=300;i++) {
		printf("%15.8f",functionValue1(Nab,i*0.1+1.0e-8));
		if((i+1)%8==0) printf("\n");
	}
	getchar();
*/
	//find the two brackets
//	printf("node1:%d, node2:%d\n",node1,node2);
	leftBracket=functionValue(Nab,xLeft);
	if(fabs(leftBracket)<0.000001) return xLeft;
	xRight=1.0;
	rightBracket=functionValue(Nab,xRight);
//	printf("rightBracket:%8.5f\n",rightBracket);
	while((leftBracket*rightBracket)>0.0){
		xLeft=xRight; leftBracket=rightBracket;
		xRight+=1.0;
		rightBracket=functionValue(Nab,xRight);
		if(fabs(rightBracket)<1.0e-6) return xRight;
//		printf("in while, xLeft:%f, xRight:%f; leftBracket:%12.8f; rightBracket:%12.8f\n",xLeft,xRight,leftBracket,rightBracket);
//		getchar();
	}
//	printf("leftBracket:%8.5f; rightBracket:%8.5f\n",leftBracket,rightBracket);
	//use bisection to solve this function
	for(j=1;j<=JMAX;j++){
		dx=xRight-xLeft;
		xMid=xLeft+dx/2.0;
		midBracket=functionValue(Nab,xMid);
//		printf("xLeft:%12.8f xRight:%12.8f xMid:%12.8f\n",xLeft,xRight,xMid);
//		printf("left:%12.8f Right:%12.8f Mid:%12.8f\n",leftBracket,rightBracket,midBracket);
//		getchar();
		if(fabs(dx)<xacc || midBracket==0.0) return xMid;
		if((leftBracket*midBracket)>0.0) { leftBracket=midBracket; xLeft=xMid; }
		else if((rightBracket*midBracket)>0.0) { rightBracket=midBracket; xRight=xMid; }
	}
	printf("Too many iterations to find the distance between node:%d and node:%d\n",node1,node2);
	return 1000000.0;
}

double functionValue(int **Nab, double d)
{
	int i,j;
	double value=0.0;
	for(i=1;i<=NP;i++){
		for(j=1;j<=i;j++){
			if(Nab[i][j]!=0)
				value+=Nab[i][j]*dpij(i,j,d,1.0)/pij(i,j,d,1.0);
		}
	}
	return value;
}

double functionValue1(int **Nab, double d)
{
	int i,j;
	double value=0.0;
	for(i=1;i<=NP;i++){
		for(j=1;j<=i;j++){
			if(Nab[i][j]!=0)
				value+=Nab[i][j]*log(pij(i,j,d,1.0));
		}
	}
	return value;
}

//calculate dPij(t)
//double dpij(double **v,double **inv,double *wr,int s1,int s2,double d, double alfa)
double dpij(int s1,int s2,double d, double alfa)
{
	int i;
	double pij_t=0.0;
	if((s2<1)||(s2>20)) { printf("Error in dpij function\n"); exit(1); }
	for(i=1;i<=NP;i++){
//		pij_t+=v[s1][i]*inv[i][s2]*pow(E, d*alfa*wr[i])*wr[i];
		pij_t+=v[s1][i]*inv[i][s2]*exp(d*alfa*wr[i])*wr[i];
	}
	return pij_t;
}

void searchGap()
{
	int i,j,k;
	int dashCount; //used for gap
	//find all gap sites
	for(i=1;i<=seqLength;i++) dashFlag[i]=false;
	for(i=1;i<=seqNum;i++){
		for(j=1;j<=seqLength;j++){
			if(node[i].seq_aa[j]=='-') dashFlag[j]=true;
		}
	}
	//release the site if it has just one gap
	for(i=1;i<=seqLength;i++){
		if(dashFlag[i]){
			dashCount=0;
			for(j=1;j<=seqNum;j++){
				if(node[j].seq_aa[i]=='-') dashCount++;
			}
			if(dashCount<=(seqNum-3)) calFlag[i]=true;
			else calFlag[i]=false;
		}
	}


	//find all representative amino acid number
	for(i=1;i<=seqNum;i++){
		for(j=1;j<=seqLength;j++){
			for(k=0;k<=NP;k++){
				if(node[i].seq_aa[j]==aa[k]) { node[i].seq_num[j]=k; break; }
			}
		}
	}
        //assign all site that dashFlag=true or calFlag=false
        for(i=1;i<=seqLength;i++) for(j=seqNum+1;j<=seqNum*2-2;j++) node[j].seq_num[i]=0;

}

void cal_PI()
{
	int i,j,k;
	int totalAA=0, aaNum[NP+1];//used for calculating PI vector
	sitePI=dmatrix(1,seqLength,1,NP);
	//calculate PI vector
	for(i=0;i<=NP;i++) aaNum[i]=0;
	for(i=1;i<=seqNum;i++){
		for(j=1;j<=seqLength;j++){
			if(!dashFlag[j]||calFlag[j]) { aaNum[node[i].seq_num[j]]++; if(node[i].seq_num[j]>0) totalAA++; }
		}
	}
	for(i=1;i<=NP;i++){
		PI[i]=1.0*aaNum[i]/totalAA;
	}
	//assign PI for every site
	for(i=1;i<=seqLength;i++){
		totalAA=0;
		for(j=0;j<=NP;j++) aaNum[j]=0;
		if(!dashFlag[i]||calFlag[i]){
			for(j=1;j<=seqNum;j++){
			 	aaNum[node[j].seq_num[i]]++; if(node[j].seq_num[i]>0) totalAA++; 
			}
			for(j=1;j<=NP;j++) sitePI[i][j]=1.0*aaNum[j]/totalAA;
		}
		else {
			for(j=1;j<=NP;j++) sitePI[i][j]=0.0;
		}
	}
	fprintf(resultFile,"PI vector for each site (*10^(-5))\n");
	fprintf(resultFile,"site  ");
	for(i=1;i<=NP;i++) fprintf(resultFile, "     %c",aa[i]);
	fprintf(resultFile,"\n");
	fprintf(resultFile,"All:  ");
	for(i=1;i<=NP;i++) fprintf(resultFile,"%6.0f ",PI[i]*100000.0);
	fprintf(resultFile,"\n");
	for(i=1;i<=seqLength;i++){
		fprintf(resultFile,"%5d ",i);
		for(j=1;j<=NP;j++) fprintf(resultFile,"%6.0f ",sitePI[i][j]*100000.0);
		fprintf(resultFile,"\n");
	}
}
void cal_Q()
{
	int i,j,k;
	double normalization;//used for Q matrix normalization
	printf("sequence number: %d  sequence length: %d\n",seqNum,seqLength);
	//calculate Q_ij
	for(i=1;i<=NP;i++){
		for(j=1;j<=NP;j++) Q_ij[i][j]=S_ij[i][j]*PI[j];
	}
	//calculate Q_ii
	for(i=1;i<=NP;i++){
		for(j=1;j<=NP;j++){
			if(j!=i) Q_ij[i][i]+=(-1.0)*Q_ij[i][j];
		}
	}
	//normalize Q_ij
	normalization=0.0;
	for(i=1;i<=NP;i++){
		normalization+=(-1.0)*PI[i]*Q_ij[i][i];
	}
	for(i=1;i<=NP;i++){
		for(j=1;j<=NP;j++) Q_ij[i][j]=Q_ij[i][j]/normalization;
	}
	//used for test, print some result
/*	printf("S_ij Matrix:\n");
	for(i=1;i<=NP;i++){
		for(j=1;j<=NP;j++) printf("%12.8f",S_ij[i][j]);
		printf("\n");
	}
	printf("PI Vector:\n");
	for(i=1;i<=NP;i++) printf("%12.8f",PI[i]);
	printf("\n");
	printf("Normalization: %10.8f\n",normalization);
	printf("Q_ij Matrix (after normalization) :\n");
	for(i=1;i<=NP;i++){
		for(j=1;j<=NP;j++) printf("%12.8f",Q_ij[i][j]);
		printf("\n");
	}
*/
}

void cal_Eign()
{
	int i,j,k;
	double wmax,wmin;//used for SVD algo
	double d=0.0; //used for finding inverse matrix
	int *indx;
	double check[NP+1][NP+1];
	double norSVD; //used for normalize eignvector matrix

	//calculate eignvalue
	for(i=1;i<=NP;i++) for(j=1;j<=NP;j++) a[i][j]=Q_ij[i][j];
	balanc(a,NP);
	elmhes(a,NP);
	hqr(a,NP,wr,wi);
/*	cout<<"The eignvalue\n";
	for(i=1;i<=NP;i++) printf("%12.8f",wr[i]);
	printf("\n");
*/
	//calculate eignvectors use SVD algo
//	printf("The eignvector matrix:\n");
	for(k=1;k<=NP;k++){
		for(i=1;i<=NP;i++) for(j=1;j<=NP;j++) a[i][j]=Q_ij[i][j];
		for(i=1;i<=NP;i++) a[i][i]-=wr[k];
		svdcmp(a,NP,NP,w,svdv);
		wmax=0.0;
		for(i=1;i<=NP;i++) if(w[i]>wmax) wmax=w[i];
		wmin=wmax*(1.0e-6);
		for(i=1;i<=NP;i++) if((w[i]<wmin)&&(w[i]>(-1.0*wmin))){
			for(j=1;j<=NP;j++) v[j][k]=svdv[j][i];
		}
		norSVD=0.0;
		for(i=1;i<=NP;i++) norSVD+=v[i][k];
		for(i=1;i<=NP;i++) v[i][k]=v[i][k]/norSVD;
	}
/*	for(i=1;i<=NP;i++) {
		for(j=1;j<=NP;j++) printf("%12.8f",v[i][j]);
		printf("\n");
	}
*/
	//calcuate inverse of eignvector matrix
	indx=ivector(1,NP);
	for(i=1;i<=NP;i++) for(j=1;j<=NP;j++) a[i][j]=v[i][j];
	ludcmp(a,NP,indx,&d);
	for(j=1;j<=NP;j++){
		for(i=1;i<=NP;i++) colu[i]=0.0;
		colu[j]=1.0;
		lubksb(a,NP,indx,colu);
		for(i=1;i<=NP;i++) inv[i][j]=colu[i];
	}
/*	printf("Inverse matrix of eignvector matrix\n");
	for(i=1;i<=NP;i++) {
		for(j=1;j<=NP;j++) printf("%12.8f",inv[i][j]);
		printf("\n");
	}
*/
	//check if the eignvector matrix is correct
/*	for(i=1; i<=NP; i++) for(j=1; j<=NP; j++) check[i][j]=0.0;
	for(k=1;k<=NP;k++) {
		for(i=1;i<=NP;i++) for(j=1;j<=NP;j++) a[i][j]=Q_ij[i][j];
		for(i=1;i<=NP;i++) a[i][i]-=wr[k];
		for(j=1;j<=NP;j++){
			for(i=1;i<=NP;i++){
				check[j][k]+=a[j][i]*v[i][k];
			}
		}
	}
	printf("the check matrix for eignvector\n");
	for(i=1;i<=NP;i++) {
		for(j=1;j<=NP;j++) printf("%12.8f",check[i][j]);
		printf("\n");
	}
*/
	//check if the inverse matrix is correct
/*	for(i=1; i<=NP; i++) for(j=1; j<=NP; j++) check[i][j]=0.0;
	for(i=1; i<=NP; i++) 
		for(j=1; j<=NP; j++) 
			for(k=1; k<=NP; k++){
				check[i][j]+=v[i][k]*inv[k][j];
			}
	printf("the check matrix for inverse\n");
	for(i=1;i<=NP;i++) {
		for(j=1;j<=NP;j++) printf("%12.8f",check[i][j]);
		printf("\n");
	}
*/
}

//calculate Pij(t)
//double pij(double **v,double **inv,double *wr,int s1,int s2,double d, double alfa)
double pij(int s1,int s2,double d, double alfa)
{
	int i;
	double pij_t=0.0;
	if((s2<1)||(s2>20)) return 1.0;
	for(i=1;i<=NP;i++){
//		pij_t+=v[s1][i]*inv[i][s2]*pow(E, d*alfa*wr[i]);
		pij_t+=v[s1][i]*inv[i][s2]*exp(d*alfa*wr[i]);
	}
	return pij_t;
}


struct left_right{
	bool exist;
	double data;
};

struct LR_index{
	bool flag;
	double high;
	int highIndex;
	double low;
	int lowIndex;
};


void optimalAlfa(int refineFlag, int condProb)
{
	int i,j,k;
	int gapNum=0,pairNum;
	double *correlation,norCor=0.0, *alfa;
	double totalLikely;
	double **corP;//correlation matrix of P(ij)
	bool invariantFlag;
	int nongapIndex;
	//used for optimization
	left_right *left, *right;
	LR_index leftIndex, rightIndex;
	bool flag=true,Lflag=true;//Lflag control the calculation of Li value
	int high_index,low_index,recalIndex1=-1,recalIndex2=-1;
	double alfa_low, alfa_high,high,low,highHigh,highLow,lowHigh;
	double delta=0.01;
	double *result;//store Log likelyhood for each site

	correlation= new double [seqLength+1];
	alfa= new double [seqLength+1];
	result= new double [seqLength+1];
	left =new left_right [seqLength+1];
	right= new left_right [seqLength+1];
	corP=dmatrix(1,NP,1,NP);

	printf("Optimize alfa value\n");
	//find the number of gap sites
	for(i=1;i<=seqLength;i++) if(dashFlag[i]&&!calFlag[i]) gapNum++;

	//find the correlation matrix (d=1.0 && Alfa=1.0) and make this matrix symmetric by Pij=(pij+pji)/2
	for(i=1;i<=NP;i++) for(j=1;j<=NP;j++) corP[i][j]=pij(i,j,1.0,1.0);
	for(i=1;i<=NP;i++) for(j=1;j<=i;j++) { corP[i][j]=(corP[i][j]+corP[j][i])/2.0; corP[j][i]=corP[i][j]; }

	//find the correlation vector
	for(i=1;i<=seqLength;i++){
		correlation[i]=0.0;
		if(!dashFlag[i]||calFlag[i]) {
			pairNum=0;
			for(j=1;j<seqNum;j++){
				for(k=j+1;k<=seqNum;k++){
					//deal with gaps
					if((node[j].seq_num[i]!=0)&&(node[k].seq_num[i]!=0))
					{ correlation[i]+=corP[node[j].seq_num[i]][node[k].seq_num[i]]; pairNum++;}
				}
			}
			correlation[i]=1.0*pairNum/correlation[i];
			norCor+=correlation[i];
		}
	}
	//set set invariant site alfa=0.0
	for(i=1;i<=seqLength;i++){
		invariantFlag=true;
		if(!dashFlag[i]){
			for(j=2;j<=seqNum;j++){
				if(node[j].seq_num[i]!=node[1].seq_num[i]) { invariantFlag=false; break;}
			}
			if(invariantFlag){
				norCor-=correlation[i];
				correlation[i]=0.0;
			}
			invariantFlag=true;
		}
		else if(calFlag[i]){
			//find the first non-gap sequence
			for(j=1;j<=seqNum;j++){
				if(node[j].seq_num[i]!=0) {nongapIndex=j;break;}
			}
			for(j=nongapIndex+1;j<=seqNum;j++){
				if((node[j].seq_num[i]!=node[nongapIndex].seq_num[i])&&(node[j].seq_num[i]!=0)){
					invariantFlag=false; break;
				}
			}
			if(invariantFlag){
				norCor-=correlation[i];
				correlation[i]=0.0;
			}
			invariantFlag=true;
		}
	}

	//initial alfa vector
	for(i=1;i<=seqLength;i++) {
		if(!dashFlag[i]||calFlag[i]) { 
			correlation[i]=correlation[i]*(1.0*(seqLength-gapNum))/norCor;
			alfa[i]=correlation[i];
		}
		else alfa[i]=0.0;
	}

	//calculate the LOG likelyhood
	totalLikely=0.0;
	for(i=1;i<=seqLength;i++){
		result[i]=0.0;
		if(!dashFlag[i]||calFlag[i]) {result[i]=calTree(i,treeRoot[i],alfa[i]); totalLikely+=result[i];}
//		printf("siteNo: %d Likelihood: %12.5f\n", i, result[i]);
	}
	fprintf(resultFile,"\n\nThe original LOG Likelyhood: %21.15f\n",totalLikely);

	//if refineFlag==0, refine alfa
	if(refineFlag==0){

		//optimize the solution
		while(flag){
			//calculate the Li(alfa-delta) and Li(alfa+delta) value
			for(i=1;i<=seqLength;i++){
				if(Lflag||(i==recalIndex1)||(i==recalIndex2)){
					left[i].exist=false; right[i].exist=false;
					if(!dashFlag[i]||calFlag[i]){
						if(alfa[i]>=delta) {
							left[i].exist=true;
							right[i].exist=true;
							alfa_low=alfa[i]-delta;
							alfa_high=alfa[i]+delta;
							left[i].data=calTree(i,treeRoot[i],alfa_low);
							right[i].data=calTree(i,treeRoot[i],alfa_high);
						}
						else{
							right[i].exist=true;
							alfa_high=alfa[i]+delta;
							right[i].data=calTree(i,treeRoot[i],alfa_high);
						}
					}//end if(dashFlag...
				}//end if(Lflag...
			}//end for
			Lflag=false;

			//search maximum uphill and minimum downhill for both sides
			leftIndex.flag=false; leftIndex.high=0.0; leftIndex.low=-1000.0;
			rightIndex.flag=false; rightIndex.high=0.0; rightIndex.low=-1000.0;
			for(i=1;i<=seqLength;i++){
				if(!dashFlag[i]||calFlag[i]){
					if(left[i].exist){
						leftIndex.flag=true;
						if((result[i]-left[i].data)<0) { 
							if((left[i].data-result[i])>leftIndex.high){
								leftIndex.highIndex=i;
								leftIndex.high=left[i].data-result[i];
							}
						}
						else if((result[i]-left[i].data)>0) {
							if((left[i].data-result[i])>leftIndex.low){
								leftIndex.lowIndex=i;
								leftIndex.low=left[i].data-result[i];
							}
						}
					}
					if(right[i].exist){
						rightIndex.flag=true;
						if((result[i]-right[i].data)<0) { 
							if((right[i].data-result[i])>rightIndex.high){
								rightIndex.highIndex=i;
								rightIndex.high=right[i].data-result[i];
							}
						}
						else if((result[i]-right[i].data)>0) {
							if((right[i].data-result[i])>leftIndex.low){
								rightIndex.lowIndex=i;
								rightIndex.low=right[i].data-result[i];
							}
						}
					}
				}
			}
			//pick up two sequence, one alfa minuses delta and one alfa adds delta,  
			//keep total alfa equals to K, but increase ln(L) value
			if(leftIndex.flag&&rightIndex.flag){
				highHigh=leftIndex.high+rightIndex.high; highLow=leftIndex.high+rightIndex.low;
				lowHigh=leftIndex.low+rightIndex.high;
				if((highHigh>0.0)||(highLow>0.0)||(lowHigh>0.0)){
					if((highHigh>lowHigh)&&(highHigh>highLow)&&(leftIndex.high>0.00000000001)&&(rightIndex.high>0.0000000001)) {
						alfa[rightIndex.highIndex]+=delta;
						result[rightIndex.highIndex]=right[rightIndex.highIndex].data;
						alfa[leftIndex.highIndex]-=delta;
						result[leftIndex.highIndex]=left[leftIndex.highIndex].data;
						recalIndex1=rightIndex.highIndex; recalIndex2=leftIndex.highIndex;
					}
					else{
						if((highLow>lowHigh)&&(highLow>0.0)){
							alfa[leftIndex.highIndex]-=delta;
							result[leftIndex.highIndex]=left[leftIndex.highIndex].data;
							alfa[rightIndex.lowIndex]+=delta;
							result[rightIndex.lowIndex]=right[rightIndex.lowIndex].data;
							recalIndex1=leftIndex.highIndex; recalIndex2=rightIndex.lowIndex;
						}
						else if((lowHigh>0.0)&&(lowHigh>highLow)){
							alfa[rightIndex.highIndex]+=delta;
							result[rightIndex.highIndex]=right[rightIndex.highIndex].data;
							alfa[leftIndex.lowIndex]-=delta;
							result[leftIndex.lowIndex]=left[leftIndex.lowIndex].data;
							recalIndex1=rightIndex.highIndex; recalIndex2=leftIndex.lowIndex;
						}
						else flag=false;
					}
				}
				else flag=false;
			}
			else flag=false;

			totalLikely=0.0;
			for(i=1;i<=seqLength;i++){
				totalLikely+=result[i];
			}
		}
	}

	//print correlation and solution (Alfa) vector
/*	for(i=1;i<=seqLength;i++){
		printf("%12.5f %12.5f\n",correlation[i],alfa[i]);
	}
*/
	//print LOG likelyhood for every site
	fprintf(resultFile,"\n\nThe final LOG likelyhood: %18.12f\n",totalLikely);
	fprintf(resultFile," site");
	for(j=1;j<=seqNum;j++) fprintf(resultFile,"%c ",' ');
	fprintf(resultFile,"    alfa    LOG(likelihood)\n");
	for(i=1;i<=seqLength;i++){
		fprintf(resultFile,"%5d ",i);
		for(j=1;j<=seqNum;j++) fprintf(resultFile,"%c ",node[j].seq_aa[i]);
		fprintf(resultFile,"%8.5f   %12.8f\n",alfa[i],result[i]);
	}

	//find the optimal letter for each internal node
	for(i=1;i<=seqLength;i++){
		if((!dashFlag[i]||calFlag[i])&&(condProb==1)) optTree(i,alfa[i]);
		if((!dashFlag[i]||calFlag[i])&&(condProb==0)) optTree_Goldstein(i,alfa[i]);
	}
	printf("finish optTree\n");
	for(i=seqNum+1;i<=(seqNum-1)*2;i++){
		node[i].seq_aa[0]=' ';
		for(j=1;j<=seqLength;j++){
//			if(!dashFlag[j]) node[i].seq_aa[j]=aa[node[i].seq_num[j]];
//			else node[i].seq_aa[j]='-';
			node[i].seq_aa[j]=aa[node[i].seq_num[j]];
		}
	}
	printf("finish assign letter\n");
	//print the transition probability of root node
	if(condProb==0){
		fprintf(resultFile,"\n\n\nThe transition probability:\n");
		for(k=seqNum+1;k<=(seqNum-1)*2;k++){
			fprintf(resultFile,"Sequence Number: %4d  (probability*10^(-4))\n",k);
			fprintf(resultFile,"  SN  AA");
			for(i=1;i<=NP;i++) fprintf(resultFile,"     %c",aa[i]);
			fprintf(resultFile,"\n");
			for(i=1;i<=seqLength;i++){
				fprintf(resultFile,"%4d  %c ",i,node[k].seq_aa[i]);
				for(j=1;j<=NP;j++){
					fprintf(resultFile,"%6.0f",node[k].prob[i][j]*10000);
				}
				fprintf(resultFile,"\n");
			}
			fprintf(resultFile,"\n\n");
		}
	}
	//print output for correlation test
/*	if(condProb==0){
		for(k=seqNum+1;k<=(seqNum-1)*2;k++){
			for(i=1;i<=seqLength;i++){
				for(j=1;j<=NP;j++) printf("%8.5f\n",node[k].prob[i][j]);
			}
		}
	}
*/


	//print the optimal solution
	fprintf(resultFile,"\n\nThe prediction for all ancestral nodes:\n");
	for(i=1;i<=(seqNum-1)*2;i++){
		fprintf(resultFile,"%5d %20s %s\n",i,node[i].nodeName,node[i].seq_aa);
	}

}

//find the solution if alfa equals 1
void noOptimalAlfa(int condProb)
{
	int i,j,k;
	double *alfa;
	double totalLikely;
	int gapNum=0;
	alfa= new double [seqLength];

	printf("Use alfa=1.0 to find result\n");

	//find the number of gap sites
	for(i=1;i<=seqLength;i++) if(dashFlag[i]) gapNum++;

	//initial alfa equals to 1
	for(i=1;i<=seqLength;i++) alfa[i]=1.0;

	//calculate the total likelyhood
	totalLikely=0.0;
	for(i=1;i<=seqLength;i++){
		if(!dashFlag[i]||calFlag[i]) totalLikely+=calTree(i,treeRoot[i],alfa[i]);
	}
	fprintf(resultFile,"\n\nThe LOG likelihood: %12.5f\n\n",totalLikely);
	//find the optimal letter for each internal node
	for(i=1;i<=seqLength;i++){
//		printf("find the result\n");
		if((!dashFlag[i]||calFlag[i])&&(condProb==1)) optTree(i,alfa[i]);
		if((!dashFlag[i]||calFlag[i])&&(condProb==0)) optTree_Goldstein(i,alfa[i]);
	}
	for(i=seqNum+1;i<=(seqNum-1)*2;i++){
		node[i].seq_aa[0]=' ';
		for(j=1;j<=seqLength;j++){
//			if(!dashFlag[j]) node[i].seq_aa[j]=aa[node[i].seq_num[j]];
//			else node[i].seq_aa[j]='-';
			node[i].seq_aa[j]=aa[node[i].seq_num[j]];
		}
	}

	//print the transition probability of internal node
	if(condProb==0){
		fprintf(resultFile,"\n\n\nThe transition probability:\n");
		for(k=seqNum+1;k<=(seqNum-1)*2;k++){
			fprintf(resultFile,"Sequence Number: %4d  (probability*10^(-4))\n",k);
			fprintf(resultFile,"  SN  AA");
			for(i=1;i<=NP;i++) fprintf(resultFile,"     %c",aa[i]);
			fprintf(resultFile,"\n");
			for(i=1;i<=seqLength;i++){
				fprintf(resultFile,"%4d  %c ",i,node[k].seq_aa[i]);
				for(j=1;j<=NP;j++){
					fprintf(resultFile,"%6.0f",node[k].prob[i][j]*10000);
				}
				fprintf(resultFile,"\n");
			}
			fprintf(resultFile,"\n\n\n");
		}
	}

	//print output for correlation test
/*	if(condProb==0){
		for(k=seqNum+1;k<=(seqNum-1)*2;k++){
			for(i=1;i<=seqLength;i++){
				for(j=1;j<=NP;j++) printf("%8.5f\n",node[k].prob[i][j]);
			}
		}
	}
*/

	//print the optimal solution
	fprintf(resultFile,"\n\nThe prediction for all ancestral nodes:\n");
	for(i=1;i<=(seqNum-1)*2;i++){
		fprintf(resultFile,"%5d %20s %s\n",i,node[i].nodeName,node[i].seq_aa);
	}


}


//calculate the LOG likelyhood for the tree at one site based on alfa value
double calTree(int siteNo, int root, double alfa)
{
	int i,j,k,l;
	int maxLevel=0;
	double siteLikely=0.0;
//	int root;
//	node[seqNum+1].amp=0.0;
//	root=treeRoot[siteNo];
	calRL(siteNo,root,alfa);
        for(i=1;i<=NP;i++){
		gapTree[siteNo][root].L[i]=PI[i]*gapTree[siteNo][gapTree[siteNo][root].childL].L[i]*gapTree[siteNo][gapTree[siteNo][root].childR].L[i]*gapTree[siteNo][gapTree[siteNo][root].childE].L[i];
            	siteLikely+=gapTree[siteNo][root].L[i];
	}
	gapTree[siteNo][root].amp=gapTree[siteNo][gapTree[siteNo][root].childL].amp+gapTree[siteNo][gapTree[siteNo][root].childR].amp;
	gapTree[siteNo][root].amp+=gapTree[siteNo][gapTree[siteNo][root].childE].amp;
	siteLikely=log(siteLikely)-gapTree[siteNo][root].amp*log(10.0);
	return siteLikely;
}


void calRL(int siteNo,int nodeNo,double alfa)
{
	if(gapTree[siteNo][nodeNo].rootFlag){
		calRL(siteNo,gapTree[siteNo][nodeNo].childL,alfa);
		calRL(siteNo,gapTree[siteNo][nodeNo].childR,alfa);
		calRL(siteNo,gapTree[siteNo][nodeNo].childE,alfa);
	}
	else if((gapTree[siteNo][nodeNo].childL!=0)&&(gapTree[siteNo][nodeNo].childR!=0)){
		gapTree[siteNo][nodeNo].amp=0.0;
		double big;
		calRL(siteNo,gapTree[siteNo][nodeNo].childL,alfa);
		calRL(siteNo,gapTree[siteNo][nodeNo].childR,alfa);
		for(int k=1;k<=NP;k++){
			gapTree[siteNo][nodeNo].L[k]=0.0;
			for(int l=1;l<=NP;l++){
				gapTree[siteNo][nodeNo].L[k]+=pij(k,l,gapTree[siteNo][nodeNo].dist,alfa)*gapTree[siteNo][gapTree[siteNo][nodeNo].childL].L[l]*gapTree[siteNo][gapTree[siteNo][nodeNo].childR].L[l];
			}
			if(k==1) big=gapTree[siteNo][nodeNo].L[k]; 
			else {
				if(gapTree[siteNo][nodeNo].L[k]>big) big=gapTree[siteNo][nodeNo].L[k]; 
			}
		}
		//find the amplifier
		while(fabs(big)<0.1){
			gapTree[siteNo][nodeNo].amp+=1.0;
			big*=10.0;
		}
		for(int i=1;i<=NP;i++) gapTree[siteNo][nodeNo].L[i]*=pow(10.0,gapTree[siteNo][nodeNo].amp);
		gapTree[siteNo][nodeNo].amp+=gapTree[siteNo][gapTree[siteNo][nodeNo].childL].amp+gapTree[siteNo][gapTree[siteNo][nodeNo].childR].amp;
	}
	else {
		gapTree[siteNo][nodeNo].amp=0.0;
		double big;
		for(int k=1;k<=NP;k++){
			gapTree[siteNo][nodeNo].L[k]=pij(k,node[nodeNo].seq_num[siteNo],gapTree[siteNo][nodeNo].dist,alfa);
			if(k==1) big=gapTree[siteNo][nodeNo].L[k]; 
			else {
				if(gapTree[siteNo][nodeNo].L[k]>big) big=gapTree[siteNo][nodeNo].L[k]; 
			}
		}
		//find the amplifier
		while(fabs(big)<0.1){
			gapTree[siteNo][nodeNo].amp+=1.0;
			big*=10.0;
		}
		for(int i=1;i<=NP;i++) gapTree[siteNo][nodeNo].L[i]*=pow(10.0,gapTree[siteNo][nodeNo].amp);
	}
}
				
//find the optimal solution for the tree at one site, based on alfa value
//implement the algorithm from paper "A Fast Algorithm for joint Reconstruction of Ancestral Amino Acid Sequences"
void optTree(int siteNo, double alfa)
{
	int i,j,k,l;
	int maxLevel=0,indexMax;
	int root,amp;
//	double **L, **C;
	double tmpL[NP+1],tmpMax;
//	L=dmatrix(1,(seqNum-1)*2,1,NP);
//	C=imatrix(1,(seqNum-1)*2,1,NP);
//	printf("In optTree function\n");
	if(!dashFlag[siteNo]){
		//calculate the L and C value for leaf node
		for(i=1;i<=seqNum;i++){
			for(k=1;k<=NP;k++)	node[i].L[k]=pij(k,node[i].seq_num[siteNo],node[i].dist,alfa);
		}
		//find the maximum level
		for(i=1;i<=(seqNum-1)*2;i++){
			if(maxLevel<node[i].nodeLevel) maxLevel=node[i].nodeLevel;
		}
		//find the max prob for each node on the tree
		if(maxLevel>1){
			for(i=maxLevel-1;i>=1;i--){
				for(j=seqNum+1;j<=(seqNum-1)*2;j++){
					if(node[j].nodeLevel==i){
						for(k=1;k<=NP;k++){
							for(l=1;l<=NP;l++){
								tmpL[l]=pij(k,l,node[j].dist,alfa)*node[node[j].childL].L[l]*node[node[j].childR].L[l];
							}
							//find the max prob
							tmpMax=tmpL[1]; indexMax=1;
							for(l=2;l<=NP;l++){
								if(tmpL[l]>tmpMax) {tmpMax=tmpL[l]; indexMax=l;}
							}
							node[j].L[k]=tmpMax;
							node[j].C[k]=indexMax;
						}
						tmpMax=node[j].L[1];
						for(k=2;k<=NP;k++) if(tmpMax<node[j].L[k]) tmpMax=node[j].L[k];
						amp=0;
						while(tmpMax<0.1){ amp++; tmpMax*=10.0;}
						for(k=1;k<=NP;k++) node[j].L[k]*=pow(10.0,amp);
					}
				}
			}
		}
		//find the max prob for root node
		for(i=1;i<=NP;i++){
			tmpL[i]=PI[i]*node[node[seqNum+1].childL].L[i]*node[node[seqNum+1].childR].L[i]*node[node[seqNum+1].childE].L[i];
		}
		tmpMax=tmpL[1]; indexMax=1;
		for(i=2;i<=NP;i++){
			if(tmpMax<tmpL[i]) { tmpMax=tmpL[i]; indexMax=i; }
		}
		
		//assign the optimal Letter to each node
		node[seqNum+1].seq_num[siteNo]=indexMax;
		for(i=1;i<maxLevel;i++){
			for(j=seqNum+1;j<=(seqNum-1)*2;j++){
				if(node[j].nodeLevel==i){
					node[j].seq_num[siteNo]=node[j].C[node[node[j].parent].seq_num[siteNo]];
				}
			}
		}
	}
	else if(calFlag[siteNo]){
		//find the root
		for(i=1;i<=(seqNum-1)*2;i++){
			if(gapTree[siteNo][i].rootFlag) { root=i; break; }
		}

		for(i=1;i<=seqNum;i++){
			if(!gapTree[siteNo][i].gapFlag){
				for(k=1;k<=NP;k++) node[i].L[k]=pij(k,node[i].seq_num[siteNo],gapTree[siteNo][i].dist,alfa);
			}
		}
		//find the maximum level
		maxLevel=0;
		for(i=1;i<=(seqNum-1)*2;i++){
			if(!gapTree[siteNo][i].gapFlag){
				if(maxLevel<gapTree[siteNo][i].nodeLevel) maxLevel=gapTree[siteNo][i].nodeLevel;
			}
		}

		if(maxLevel>1){
			for(i=maxLevel-1;i>=1;i--){
				for(j=seqNum+1;j<=(seqNum-1)*2;j++){
					if(!gapTree[siteNo][j].gapFlag){
						if(gapTree[siteNo][j].nodeLevel==i){
							for(k=1;k<=NP;k++){
//								node[j].L[k]=0.0;
								for(l=1;l<=NP;l++){
									tmpL[l]=pij(k,l,gapTree[siteNo][j].dist,alfa)*node[gapTree[siteNo][j].childL].L[l]*node[gapTree[siteNo][j].childR].L[l];
								}
								//find the max prob
								tmpMax=tmpL[1];indexMax=1;
								for(l=2;l<=NP;l++){
									if(tmpL[l]>tmpMax) {tmpMax=tmpL[l]; indexMax=l;}
								}
								node[j].L[k]=tmpMax;
								node[j].C[k]=indexMax;
							}
							tmpMax=node[j].L[1];
							for(k=2;k<=NP;k++) if(tmpMax<node[j].L[k]) tmpMax=node[j].L[k];
							amp=0;
							while(tmpMax<0.1){ amp++; tmpMax*=10.0;}
							for(k=1;k<=NP;k++) node[j].L[k]*=pow(10.0,amp);
						}
					}
				}
			}
		}
		//find the max prob for root
		for(i=1;i<=NP;i++){
			tmpL[i]=PI[i]*node[gapTree[siteNo][root].childL].L[i]*node[gapTree[siteNo][root].childR].L[i]*node[gapTree[siteNo][root].childE].L[i];
		}
		tmpMax=tmpL[1]; indexMax=1;
		for(i=2;i<=NP;i++){
			if(tmpMax<tmpL[i]) { tmpMax=tmpL[i]; indexMax=i; }
		}
		
		//assign the optimal Letter to each node
		node[root].seq_num[siteNo]=indexMax;
		for(i=1;i<maxLevel;i++){
			for(j=seqNum+1;j<=(seqNum-1)*2;j++){
				if(gapTree[siteNo][j].gapFlag) node[j].seq_num[siteNo]=0;
				else {
					if(gapTree[siteNo][j].nodeLevel==i){
						node[j].seq_num[siteNo]=node[j].C[node[gapTree[siteNo][j].parent].seq_num[siteNo]];
					}
				}
			}
		}
	}


}


//find the optimal solution for the tree at one site, based on alfa value
//likelyhood calculation from " Evolutionary Trees from DNA Sequences: A Maximum Likelihood Approach"
//conditional probability from "Probabilistic Reconstruction of Ancestral Protein Sequences"
void optTree_Goldstein(int siteNo, double alfa)
{
	int i,j;
//	int root;
	double tmpSum,tmpMax;
	int indexMax;
	topology *treeCopy;
	treeCopy = new topology [seqNum*2-1];
//	printf("siteNo: %d ",siteNo);
	if(dashFlag[siteNo]){
		for(i=1;i<=seqNum*2-2;i++){
			treeCopy[i].nodeLevel=gapTree[siteNo][i].nodeLevel;
			treeCopy[i].childL=gapTree[siteNo][i].childL;
			treeCopy[i].childR=gapTree[siteNo][i].childR;
			treeCopy[i].childE=gapTree[siteNo][i].childE;
			treeCopy[i].parent=gapTree[siteNo][i].parent;
			treeCopy[i].rootFlag=gapTree[siteNo][i].rootFlag;
			treeCopy[i].dist=gapTree[siteNo][i].dist;
			treeCopy[i].gapFlag=gapTree[siteNo][i].gapFlag;
		}
	}
/*	if(dashFlag[siteNo]) {
		for(j=seqNum+1;j<=seqNum*2-2;j++)
			printf("seqNo: %3d childL: %3d childR: %3d childE: %3d dist: %8.5f\n",j,gapTree[siteNo][j].childL,gapTree[siteNo][j].childR,gapTree[siteNo][j].childE,gapTree[siteNo][j].dist);
	} */
	for(i=seqNum+1;i<=seqNum*2-2;i++){
		if(!gapTree[siteNo][i].gapFlag){
			if(!gapTree[siteNo][i].rootFlag) changeRoot(siteNo,i,i,i);
			calTree(siteNo,i,alfa);
			//normalize the probability and put into prob[] vector
			tmpSum=0.0;
			for(j=1;j<=NP;j++){
				tmpSum+=gapTree[siteNo][i].L[j];
			}
			for(j=1;j<=NP;j++) node[i].prob[siteNo][j]=gapTree[siteNo][i].L[j]/tmpSum;
			//find the max prob for root node and assign the letter to root node
			tmpMax=node[i].prob[siteNo][1]; indexMax=1;
			for(j=2;j<=NP;j++){
				if(node[i].prob[siteNo][j]>tmpMax) { tmpMax=node[i].prob[siteNo][j]; indexMax=j;}
			}
			for(j=1;j<=NP;j++) if(node[i].prob[siteNo][j]<0.00000001) node[i].prob[siteNo][j]=0.0;
			node[i].seq_num[siteNo]=indexMax;
		}
	}
	if(dashFlag[siteNo]){
		for(i=1;i<=seqNum*2-2;i++){
			gapTree[siteNo][i].nodeLevel=treeCopy[i].nodeLevel;
			gapTree[siteNo][i].childL=treeCopy[i].childL;
			gapTree[siteNo][i].childR=treeCopy[i].childR;
			gapTree[siteNo][i].childE=treeCopy[i].childE;
			gapTree[siteNo][i].parent=treeCopy[i].parent;
			gapTree[siteNo][i].rootFlag=treeCopy[i].rootFlag;
			gapTree[siteNo][i].dist=treeCopy[i].dist;
			gapTree[siteNo][i].gapFlag=treeCopy[i].gapFlag;
		}
	}
}

void changeRoot(int siteNo, int nodeNo, int root, int child)
{
	int parent;
	if(nodeNo==root){
		parent=gapTree[siteNo][nodeNo].parent;
		changeRoot(siteNo,parent,root,nodeNo);
		gapTree[siteNo][nodeNo].childE=parent;
		gapTree[siteNo][nodeNo].dist=0.0;
		gapTree[siteNo][nodeNo].rootFlag=true;
		gapTree[siteNo][nodeNo].parent=nodeNo;
	}
	else if(!gapTree[siteNo][nodeNo].rootFlag){
		changeRoot(siteNo, gapTree[siteNo][nodeNo].parent,root,nodeNo);
		if(child==gapTree[siteNo][nodeNo].childL)
			gapTree[siteNo][nodeNo].childL=gapTree[siteNo][nodeNo].parent;
		else gapTree[siteNo][nodeNo].childR=gapTree[siteNo][nodeNo].parent;
		gapTree[siteNo][nodeNo].parent=child;
		gapTree[siteNo][nodeNo].dist=gapTree[siteNo][child].dist;
	}
	else {
		if(child==gapTree[siteNo][nodeNo].childL)
			gapTree[siteNo][nodeNo].childL=gapTree[siteNo][nodeNo].childE;
		else if(child==gapTree[siteNo][nodeNo].childR)
			gapTree[siteNo][nodeNo].childR=gapTree[siteNo][nodeNo].childE;
		gapTree[siteNo][nodeNo].childE=0;
		gapTree[siteNo][nodeNo].parent=child;
		gapTree[siteNo][nodeNo].dist=gapTree[siteNo][child].dist;
		gapTree[siteNo][nodeNo].rootFlag=false;
	}

}

int rebuildTree(int siteNo)
{
	int i,j;
	int maxLevel=0;
	int ChildL,ChildR,ChildE,CCR,CCL;
	bool childLFlag,childRFlag,childEFlag;
	int root,levelIndex,levelP;

//	printf("Rebuild tree structure for site: %d\n",siteNo);
	//initiate the tree
/*	for(i=1;i<=(seqNum-1)*2;i++){
		gapTree[siteNo][i].nodeLevel=node[i].nodeLevel;
		gapTree[siteNo][i].childL=node[i].childL;
		gapTree[siteNo][i].childR=node[i].childR;
		gapTree[siteNo][i].childE=node[i].childE;
		gapTree[siteNo][i].parent=node[i].parent;
		gapTree[siteNo][i].rootFlag=node[i].rootFlag;
		gapTree[siteNo][i].dist=node[i].dist;
	}
	for(i=1;i<=seqNum;i++){
		if(node[i].seq_aa[siteNo]=='-') gapTree[siteNo][i].gapFlag=true;
		else gapTree[siteNo][i].gapFlag=false;
	}
	for(i=seqNum+1;i<=(seqNum-1)*2;i++) gapTree[siteNo][i].gapFlag=false;
*/	//find the maximum level
	for(i=1;i<=(seqNum-1)*2;i++){
		if(maxLevel<node[i].nodeLevel) maxLevel=node[i].nodeLevel;
	}

	//rebuild tree
	if(maxLevel>1){
		for(i=maxLevel-1; i>=1; i--){
			for(j=seqNum+1;j<=(seqNum-1)*2;j++){
				if(gapTree[siteNo][j].nodeLevel==i){
					if(gapTree[siteNo][gapTree[siteNo][j].childL].gapFlag&&gapTree[siteNo][gapTree[siteNo][j].childR].gapFlag){
						gapTree[siteNo][j].childL=0; gapTree[siteNo][j].childR=0;
						gapTree[siteNo][j].gapFlag=true;
					}
					else if(!gapTree[siteNo][gapTree[siteNo][j].childL].gapFlag&&gapTree[siteNo][gapTree[siteNo][j].childR].gapFlag){
						gapTree[siteNo][gapTree[siteNo][j].childL].parent=gapTree[siteNo][j].parent;
						gapTree[siteNo][gapTree[siteNo][j].childL].dist+=gapTree[siteNo][j].dist;
						gapTree[siteNo][gapTree[siteNo][j].childL].nodeLevel--;
						if(gapTree[siteNo][gapTree[siteNo][j].parent].childL==j)
							gapTree[siteNo][gapTree[siteNo][j].parent].childL=gapTree[siteNo][j].childL;
						else if(gapTree[siteNo][gapTree[siteNo][j].parent].childR==j)
							gapTree[siteNo][gapTree[siteNo][j].parent].childR=gapTree[siteNo][j].childL;
						else gapTree[siteNo][gapTree[siteNo][j].parent].childE=gapTree[siteNo][j].childL;
						gapTree[siteNo][j].childL=0; gapTree[siteNo][j].childR=0;
						gapTree[siteNo][j].gapFlag=true;
					}
					else if(gapTree[siteNo][gapTree[siteNo][j].childL].gapFlag&&!gapTree[siteNo][gapTree[siteNo][j].childR].gapFlag){
						gapTree[siteNo][gapTree[siteNo][j].childR].parent=gapTree[siteNo][j].parent;
						gapTree[siteNo][gapTree[siteNo][j].childR].dist+=gapTree[siteNo][j].dist;
						gapTree[siteNo][gapTree[siteNo][j].childR].nodeLevel--;
						if(gapTree[siteNo][gapTree[siteNo][j].parent].childL==j)
							gapTree[siteNo][gapTree[siteNo][j].parent].childL=gapTree[siteNo][j].childR;
						else if(gapTree[siteNo][gapTree[siteNo][j].parent].childR==j)
							gapTree[siteNo][gapTree[siteNo][j].parent].childR=gapTree[siteNo][j].childR;
						else gapTree[siteNo][gapTree[siteNo][j].parent].childE=gapTree[siteNo][j].childR;
						gapTree[siteNo][j].childL=0; gapTree[siteNo][j].childR=0;
						gapTree[siteNo][j].gapFlag=true;
					}
					else gapTree[siteNo][j].gapFlag=false;
				}
			}
		}
	}

	//deal with root node
	ChildL=gapTree[siteNo][seqNum+1].childL; childLFlag=gapTree[siteNo][ChildL].gapFlag;
	ChildR=gapTree[siteNo][seqNum+1].childR; childRFlag=gapTree[siteNo][ChildR].gapFlag;
	ChildE=gapTree[siteNo][seqNum+1].childE; childEFlag=gapTree[siteNo][ChildE].gapFlag;
	if(childLFlag || childRFlag || childEFlag) {
		gapTree[siteNo][seqNum+1].gapFlag=true;
		gapTree[siteNo][seqNum+1].rootFlag=false;
		if(childLFlag){
			if(!childRFlag&&!childEFlag){
				if(ChildR<=seqNum){
					gapTree[siteNo][ChildR].parent=ChildE;
					gapTree[siteNo][ChildE].childE=ChildR;
					gapTree[siteNo][ChildE].rootFlag=true;
					gapTree[siteNo][ChildR].dist+=gapTree[siteNo][ChildE].dist;
					gapTree[siteNo][ChildE].dist=0.0;
				}
				else {
					gapTree[siteNo][ChildE].parent=ChildR;
					gapTree[siteNo][ChildR].childE=ChildE;
					gapTree[siteNo][ChildR].rootFlag=true;
					gapTree[siteNo][ChildE].dist+=gapTree[siteNo][ChildR].dist;
					gapTree[siteNo][ChildR].dist=0.0;
				}
			}
			else if(!childRFlag){
				gapTree[siteNo][ChildR].gapFlag=true;
				CCL=gapTree[siteNo][ChildR].childL;
				CCR=gapTree[siteNo][ChildR].childR;
				if(CCL<=seqNum){
					gapTree[siteNo][CCL].parent=CCR;
					gapTree[siteNo][CCR].childE=CCL;
					gapTree[siteNo][CCR].rootFlag=true;
					gapTree[siteNo][CCL].dist+=gapTree[siteNo][CCR].dist;
					gapTree[siteNo][CCR].dist=0.0;
				}
				else {
					gapTree[siteNo][CCR].parent=CCL;
					gapTree[siteNo][CCL].childE=CCR;
					gapTree[siteNo][CCL].rootFlag=true;
					gapTree[siteNo][CCR].dist+=gapTree[siteNo][CCL].dist;
					gapTree[siteNo][CCL].dist=0.0;
				}
			}
			else if(!childEFlag){
				gapTree[siteNo][ChildE].gapFlag=true;
				CCL=gapTree[siteNo][ChildE].childL;
				CCR=gapTree[siteNo][ChildE].childR;
				if(CCL<=seqNum){
					gapTree[siteNo][CCL].parent=CCR;
					gapTree[siteNo][CCR].childE=CCL;
					gapTree[siteNo][CCR].rootFlag=true;
					gapTree[siteNo][CCL].dist+=gapTree[siteNo][CCR].dist;
					gapTree[siteNo][CCR].dist=0.0;
				}
				else {
					gapTree[siteNo][CCR].parent=CCL;
					gapTree[siteNo][CCL].childE=CCR;
					gapTree[siteNo][CCL].rootFlag=true;
					gapTree[siteNo][CCR].dist+=gapTree[siteNo][CCL].dist;
					gapTree[siteNo][CCL].dist=0.0;
				}
			}
		}
		else if(!childLFlag){
			if(childRFlag&&childEFlag){
				gapTree[siteNo][ChildL].gapFlag=true;
				CCL=gapTree[siteNo][ChildL].childL;
				CCR=gapTree[siteNo][ChildL].childR;
				if(CCL<=seqNum){
					gapTree[siteNo][CCL].parent=CCR;
					gapTree[siteNo][CCR].childE=CCL;
					gapTree[siteNo][CCR].rootFlag=true;
					gapTree[siteNo][CCL].dist+=gapTree[siteNo][CCR].dist;
					gapTree[siteNo][CCR].dist=0.0;
				}
				else {
					gapTree[siteNo][CCR].parent=CCL;
					gapTree[siteNo][CCL].childE=CCR;
					gapTree[siteNo][CCL].rootFlag=true;
					gapTree[siteNo][CCR].dist+=gapTree[siteNo][CCL].dist;
					gapTree[siteNo][CCL].dist=0.0;
				}
			}
			else if(!childRFlag){
				if(ChildL<=seqNum){
					gapTree[siteNo][ChildL].parent=ChildR;
					gapTree[siteNo][ChildR].childE=ChildL;
					gapTree[siteNo][ChildR].rootFlag=true;
					gapTree[siteNo][ChildL].dist+=gapTree[siteNo][ChildR].dist;
					gapTree[siteNo][ChildE].dist=0.0;
				}
				else {
					gapTree[siteNo][ChildR].parent=ChildL;
					gapTree[siteNo][ChildL].childE=ChildR;
					gapTree[siteNo][ChildL].rootFlag=true;
					gapTree[siteNo][ChildR].dist+=gapTree[siteNo][ChildL].dist;
					gapTree[siteNo][ChildL].dist=0.0;
				}
			}
			else if(!childEFlag){
				if(ChildL<=seqNum){
					gapTree[siteNo][ChildL].parent=ChildE;
					gapTree[siteNo][ChildE].childE=ChildL;
					gapTree[siteNo][ChildE].rootFlag=true;
					gapTree[siteNo][ChildL].dist+=gapTree[siteNo][ChildE].dist;
					gapTree[siteNo][ChildE].dist=0.0;
				}
				else {
					gapTree[siteNo][ChildE].parent=ChildL;
					gapTree[siteNo][ChildL].childE=ChildE;
					gapTree[siteNo][ChildL].rootFlag=true;
					gapTree[siteNo][ChildE].dist+=gapTree[siteNo][ChildL].dist;
					gapTree[siteNo][ChildL].dist=0.0;
				}
			}
		}
	}
	//redefine every one's node level
	//find the root
	for(i=1;i<=(seqNum-1)*2;i++){
		if(gapTree[siteNo][i].rootFlag) {root=i; break;}
	}
//	printf("root is %3d\n",root);
	//set up root's parent is itself
	gapTree[siteNo][root].parent=root;

//	printf("Finish rebuild tree, before reassign level\n");
	setRLevel(siteNo,root);
	return root;

/*	if(siteNo==125){
	for(i=1;i<=(seqNum-1)*2;i++){
		if(gapTree[siteNo][i].gapFlag)	fprintf(resultFile,"node: %3d has a gap\n",i);
		else {
			if(i<=seqNum) fprintf(resultFile,"node: %3d parent %3d distance %8.5f",i,gapTree[siteNo][i].parent,gapTree[siteNo][i].dist);
			else if(i==root)
				fprintf(resultFile,"node: %3d childL: %3d childR: %3d childE: %3d",i,gapTree[siteNo][i].childL,gapTree[siteNo][i].childR,gapTree[siteNo][i].childE);
			else 
				fprintf(resultFile,"node: %3d childL: %3d childR: %3d parent: %3d",i,gapTree[siteNo][i].childL,gapTree[siteNo][i].childR,gapTree[siteNo][i].parent);
			fprintf(resultFile," Level: %3d\n",gapTree[siteNo][i].nodeLevel);
		}
	}
	}
*/
}

void setRLevel(int siteNo,int nodeNum)
{
	if(gapTree[siteNo][nodeNum].rootFlag){
		gapTree[siteNo][nodeNum].nodeLevel=0;
		setRLevel(siteNo,gapTree[siteNo][nodeNum].childL);
		setRLevel(siteNo,gapTree[siteNo][nodeNum].childR);
		setRLevel(siteNo,gapTree[siteNo][nodeNum].childE);
	}
	else if((gapTree[siteNo][nodeNum].childL!=0)&&(gapTree[siteNo][nodeNum].childR!=0)) {
		gapTree[siteNo][nodeNum].nodeLevel=gapTree[siteNo][gapTree[siteNo][nodeNum].parent].nodeLevel+1;
		setRLevel(siteNo,gapTree[siteNo][nodeNum].childL);
		setRLevel(siteNo,gapTree[siteNo][nodeNum].childR);
	}
	else 
		gapTree[siteNo][nodeNum].nodeLevel=gapTree[siteNo][gapTree[siteNo][nodeNum].parent].nodeLevel+1;
}
		

//*levelFlag=new bool [DIM];
//for(int i=0;i<DIM;i++) levelFlag[i]=new bool [DIM];
	
void print_tree(int nodeNo)
{
	int i,j;
//	printf("in printTree function, nodeNo:%d childL:%d childR:%d\n",nodeNo,node[nodeNo].childL,node[nodeNo].childR);
	if((nodeNo<1)||(nodeNo>(seqNum*2-2))) { printf("There is something wrong in printing tree\n"); exit(1); }
	if((nodeNo>=1)&&(nodeNo<=seqNum)) {
		for(j=0;j<=node[nodeNo].nodeLevel-2;j++){
			if(levelFlag[nodeNo][j]==1) fprintf(resultFile,"|");
			else fprintf(resultFile," ");
			for(i=1;i<=9;i++) fprintf(resultFile," ");
		}
//		for(i=1;i<=node[node[nodeNo].parent].nodeLevel*10;i++) fprintf(resultFile," ");
		fprintf(resultFile,"+");
		for(i=1;i<=9;i++) fprintf(resultFile,"-");
		fprintf(resultFile,"%d:(%s)\n",nodeNo,node[nodeNo].nodeName);
	}
	else {
//		printf("go to subtree, childL:%d, childR:%d, childE:%d\n",node[nodeNo].childL,node[nodeNo].childR,node[nodeNo].childE);
		if(node[nodeNo].rootFlag){
			print_tree(node[nodeNo].childL);
			fprintf(resultFile,"%d\n",nodeNo);
			print_tree(node[nodeNo].childR);
			print_tree(node[nodeNo].childE);
		}
		else {
			if(nodeNo==node[node[nodeNo].parent].childL){
				for(i=0;i<node[node[nodeNo].parent].nodeLevel;i++)
				{ 
					levelFlag[node[nodeNo].childL][i]=levelFlag[nodeNo][i];
					levelFlag[node[nodeNo].childR][i]=levelFlag[nodeNo][i];
				}
				levelFlag[node[nodeNo].childL][node[nodeNo].nodeLevel-1]=0;//false;
				print_tree(node[nodeNo].childL);
				for(j=0;j<=node[nodeNo].nodeLevel-2;j++){
					if(levelFlag[nodeNo][j]) fprintf(resultFile,"|");
					else fprintf(resultFile," ");
					for(i=1;i<=9;i++) fprintf(resultFile," ");
				}
				fprintf(resultFile,"+");
				for(i=1;i<=9;i++) fprintf(resultFile,"-");
				fprintf(resultFile,"%d\n",nodeNo);
				levelFlag[node[nodeNo].childR][node[nodeNo].nodeLevel-1]=1;//true;
				print_tree(node[nodeNo].childR);
			}
			else if(nodeNo==node[node[nodeNo].parent].childR){
				for(i=0;i<node[node[nodeNo].parent].nodeLevel;i++)
				{ 
					levelFlag[node[nodeNo].childL][i]=levelFlag[nodeNo][i];
					levelFlag[node[nodeNo].childR][i]=levelFlag[nodeNo][i];
				}
				levelFlag[node[nodeNo].childL][node[nodeNo].nodeLevel-1]=1;//true;
				print_tree(node[nodeNo].childL);
				for(j=0;j<=node[nodeNo].nodeLevel-2;j++){
					if(levelFlag[nodeNo][j]==1) fprintf(resultFile,"|");
					else fprintf(resultFile," ");
					for(i=1;i<=9;i++) fprintf(resultFile," ");
				}
				fprintf(resultFile,"+");
				for(i=1;i<=9;i++) fprintf(resultFile,"-");
				fprintf(resultFile,"%d\n",nodeNo);
				if(node[nodeNo].parent==(seqNum+1))
					levelFlag[node[nodeNo].childR][node[nodeNo].nodeLevel-1]=1;//true;
				else levelFlag[node[nodeNo].childR][node[nodeNo].nodeLevel-1]=0;//false;
				print_tree(node[nodeNo].childR);
			}
			else if(nodeNo==node[node[nodeNo].parent].childE){
				for(i=0;i<node[node[nodeNo].parent].nodeLevel;i++)
				{ 
					levelFlag[node[nodeNo].childL][i]=levelFlag[nodeNo][i];
					levelFlag[node[nodeNo].childR][i]=levelFlag[nodeNo][i];
				}
				levelFlag[node[nodeNo].childL][node[nodeNo].nodeLevel-1]=1;//true;
				print_tree(node[nodeNo].childL);
				for(j=0;j<=node[nodeNo].nodeLevel-2;j++){
					if(levelFlag[nodeNo][j]) fprintf(resultFile,"|");
					else fprintf(resultFile," ");
					for(i=1;i<=9;i++) fprintf(resultFile," ");
				}
				fprintf(resultFile,"+");
				for(i=1;i<=9;i++) fprintf(resultFile,"-");
				fprintf(resultFile,"%d\n",nodeNo);
				levelFlag[node[nodeNo].childR][node[nodeNo].nodeLevel-1]=0;//false;
				print_tree(node[nodeNo].childR);
			}
		}
//		for(i=1;i<=node[node[nodeNo].parent].nodeLevel*10;i++) printf(" ");
/*		printf("+");
		for(i=1;i<=9;i++) printf("-");

		print_tree(node[nodeNo].childL);
		if(nodeNo!=seqNum+1){
			for(i=1;i<=node[node[nodeNo].parent].nodeLevel*10;i++) printf(" ");
			printf("+");
			for(i=1;i<=9;i++) printf("-");
		}
		printf("%d\n",nodeNo);
		print_tree(node[nodeNo].childR);
		if(node[nodeNo].childE!=0) printTree(node[nodeNo].childE);
*/
	}
}

void init_S()
{
	int i,j;
	S_ij[2][1]=  0.551571;S_ij[3][1]=  0.509848;S_ij[3][2]=  0.635346;S_ij[4][1]=  0.738998;S_ij[4][2]=  0.147304;
	S_ij[4][3]=  5.429420;S_ij[5][1]=  1.027040;S_ij[5][2]=  0.528191;S_ij[5][3]=  0.265256;S_ij[5][4]=  0.030295;
	S_ij[6][1]=  0.908598;S_ij[6][2]=  3.035500;S_ij[6][3]=  1.543640;S_ij[6][4]=  0.616783;S_ij[6][5]=  0.098818;
	S_ij[7][1]=  1.582850;S_ij[7][2]=  0.439157;S_ij[7][3]=  0.947198;S_ij[7][4]=  6.174160;S_ij[7][5]=  0.021352;
	S_ij[7][6]=  5.469470;S_ij[8][1]=  1.416720;S_ij[8][2]=  0.584665;S_ij[8][3]=  1.125560;S_ij[8][4]=  0.865584;
	S_ij[8][5]=  0.306674;S_ij[8][6]=  0.330052;S_ij[8][7]=  0.567717;S_ij[9][1]=  0.316954;S_ij[9][2]=  2.137150;
	S_ij[9][3]=  3.956290;S_ij[9][4]=  0.930676;S_ij[9][5]=  0.248972;S_ij[9][6]=  4.294110;S_ij[9][7]=  0.570025;
	S_ij[9][8]=  0.249410;S_ij[10][1]=  0.193335;S_ij[10][2]=  0.186979;S_ij[10][3]=  0.554236;S_ij[10][4]=  0.039437;
	S_ij[10][5]=  0.170135;S_ij[10][6]=  0.113917;S_ij[10][7]=  0.127395;S_ij[10][8]=  0.030450;S_ij[10][9]=  0.138190;
	S_ij[11][1]=  0.397915;S_ij[11][2]=  0.497671;S_ij[11][3]=  0.131528;S_ij[11][4]=  0.084805;S_ij[11][5]=  0.384287;
	S_ij[11][6]=  0.869489;S_ij[11][7]=  0.154263;S_ij[11][8]=  0.061304;S_ij[11][9]=  0.499462;S_ij[11][10]=  3.170970;
	S_ij[12][1]=  0.906265;S_ij[12][2]=  5.351420;S_ij[12][3]=  3.012010;S_ij[12][4]=  0.479855;S_ij[12][5]=  0.074034;
	S_ij[12][6]=  3.894900;S_ij[12][7]=  2.584430;S_ij[12][8]=  0.373558;S_ij[12][9]=  0.890432;S_ij[12][10]=  0.323832;
	S_ij[12][11]=  0.257555;S_ij[13][1]=  0.893496;S_ij[13][2]=  0.683162;S_ij[13][3]=  0.198221;S_ij[13][4]=  0.103754;
	S_ij[13][5]=  0.390482;S_ij[13][6]=  1.545260;S_ij[13][7]=  0.315124;S_ij[13][8]=  0.174100;S_ij[13][9]=  0.404141;
	S_ij[13][10]=  4.257460;S_ij[13][11]=  4.854020;S_ij[13][12]=  0.934276;S_ij[14][1]=  0.210494;S_ij[14][2]=  0.102711;
	S_ij[14][3]=  0.096162;S_ij[14][4]=  0.046730;S_ij[14][5]=  0.398020;S_ij[14][6]=  0.099921;S_ij[14][7]=  0.081134;
	S_ij[14][8]=  0.049931;S_ij[14][9]=  0.679371;S_ij[14][10]=  1.059470;S_ij[14][11]=  2.115170;S_ij[14][12]=  0.088836;
	S_ij[14][13]=  1.190630;S_ij[15][1]=  1.438550;S_ij[15][2]=  0.679489;S_ij[15][3]=  0.195081;S_ij[15][4]=  0.423984;
	S_ij[15][5]=  0.109404;S_ij[15][6]=  0.933372;S_ij[15][7]=  0.682355;S_ij[15][8]=  0.243570;S_ij[15][9]=  0.696198;
	S_ij[15][10]=  0.099929;S_ij[15][11]=  0.415844;S_ij[15][12]=  0.556896;S_ij[15][13]=  0.171329;S_ij[15][14]=  0.161444;
	S_ij[16][1]=  3.370790;S_ij[16][2]=  1.224190;S_ij[16][3]=  3.974230;S_ij[16][4]=  1.071760;S_ij[16][5]=  1.407660;
	S_ij[16][6]=  1.028870;S_ij[16][7]=  0.704939;S_ij[16][8]=  1.341820;S_ij[16][9]=  0.740169;S_ij[16][10]=  0.319440;
	S_ij[16][11]=  0.344739;S_ij[16][12]=  0.967130;S_ij[16][13]=  0.493905;S_ij[16][14]=  0.545931;S_ij[16][15]=  1.613280;
	S_ij[17][1]=  2.121110;S_ij[17][2]=  0.554413;S_ij[17][3]=  2.030060;S_ij[17][4]=  0.374866;S_ij[17][5]=  0.512984;
	S_ij[17][6]=  0.857928;S_ij[17][7]=  0.822765;S_ij[17][8]=  0.225833;S_ij[17][9]=  0.473307;S_ij[17][10]=  1.458160;
	S_ij[17][11]=  0.326622;S_ij[17][12]=  1.386980;S_ij[17][13]=  1.516120;S_ij[17][14]=  0.171903;S_ij[17][15]=  0.795384;
	S_ij[17][16]=  4.378020;S_ij[18][1]=  0.113133;S_ij[18][2]=  1.163920;S_ij[18][3]=  0.071917;S_ij[18][4]=  0.129767;
	S_ij[18][5]=  0.717070;S_ij[18][6]=  0.215737;S_ij[18][7]=  0.156557;S_ij[18][8]=  0.336983;S_ij[18][9]=  0.262569;
	S_ij[18][10]=  0.212483;S_ij[18][11]=  0.665309;S_ij[18][12]=  0.137505;S_ij[18][13]=  0.515706;S_ij[18][14]=  1.529640;
	S_ij[18][15]=  0.139405;S_ij[18][16]=  0.523742;S_ij[18][17]=  0.110864;S_ij[19][1]=  0.240735;S_ij[19][2]=  0.381533;
	S_ij[19][3]=  1.086000;S_ij[19][4]=  0.325711;S_ij[19][5]=  0.543833;S_ij[19][6]=  0.227710;S_ij[19][7]=  0.196303;
	S_ij[19][8]=  0.103604;S_ij[19][9]=  3.873440;S_ij[19][10]=  0.420170;S_ij[19][11]=  0.398618;S_ij[19][12]=  0.133264;
	S_ij[19][13]=  0.428437;S_ij[19][14]=  6.454280;S_ij[19][15]=  0.216046;S_ij[19][16]=  0.786993;S_ij[19][17]=  0.291148;
	S_ij[19][18]=  2.485390;S_ij[20][1]=  2.006010;S_ij[20][2]=  0.251849;S_ij[20][3]=  0.196246;S_ij[20][4]=  0.152335;
	S_ij[20][5]=  1.002140;S_ij[20][6]=  0.301281;S_ij[20][7]=  0.588731;S_ij[20][8]=  0.187247;S_ij[20][9]=  0.118358;
	S_ij[20][10]=  7.821300;S_ij[20][11]=  1.800340;S_ij[20][12]=  0.305434;S_ij[20][13]=  2.058450;S_ij[20][14]=  0.649892;
	S_ij[20][15]=  0.314887;S_ij[20][16]=  0.232739;S_ij[20][17]=  1.388230;S_ij[20][18]=  0.365369;S_ij[20][19]=  0.314730;

	//symmetrize S_ij
	for(i=1;i<=NP;i++){
		S_ij[i][i]=0.0;
		for(j=1;j<i;j++) S_ij[j][i]=S_ij[i][j];
	}
}


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;
	DIM=(seqNum-1)*2+1;
	node=new NODE [DIM];
	for(i=1;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);
	}
	gapTree=topmatrix(0,seqLength+1,0,(seqNum*2-2));//seqLength);
	levelFlag=imatrix(1,DIM,1,DIM);
        treeRoot=ivector(1,seqLength);
	dashFlag=new bool [seqLength+1];
	calFlag= new bool [seqLength+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=1;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");
	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");
	}
//	printf("index=%d\n",index);
	//assign seq_num
	for(i=1;i<=seqNum;i++){
		for(j=1;j<=seqLength;j++){
			for(k=0;k<=NP;k++){
				if(node[i].seq_aa[j]==aa[k]) { node[i].seq_num[j]=k; break; }
			}
		}
	}
	//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);
	

}


topology **topmatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
{
	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
	topology **m;

	/* allocate pointers to rows */
	m=(topology **) malloc((size_t)((nrow+NR_END)*sizeof(topology*)));
	if (!m) nrerror("allocation failure 1 in matrix()");
	m += NR_END;
	m -= nrl;


	/* allocate rows and set pointers to them */
	m[nrl]=(topology *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(topology)));
	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
	m[nrl] += NR_END;
	m[nrl] -= ncl;

	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

	/* return pointer to array of pointers to rows */
	return m;
}

