// 08/19/02 add -G parameter to correct estimated alpha calcualation
#include "mpi.h"
#include "ancestral.h"
#include <sys/ddi.h>

struct BRANCH{
//	int No;
	double length;
	int LPoint;
	int RPoint;
	int LPointNumber;
	int RPointNumber;
};
struct CUT{
	int branchNo;
	double leftDist;
	double **prob;
	bool *gapFlag;
};

void weighbor(int N, int sequenceLength, NODE *node, double **distMatrix, char *treeFile);
void cal_dist(char *seqFile);
void assignBranch(int nodeNo);
void findMidpoint();
int reconstructTree(int siteNo, int root, bool rootedFlag);
void reconstruct(int siteNo, int nodeNo);
void calRL(int siteNo,int nodeNo, int leafSite, double alfa);
double dpij_alfa(int s1,int s2,double d, double alfa);
double calTree_alfa(int siteNo, int root, double alfa, int *sign);
void calRL_alfa(int siteNo,int nodeNo,double alfa,int deri_nodeNo);
double func_alfa(int i, int root, double alfa);
void optimalAlfa(int refineFlag, int condProb);
double optimizeAlpha(int site);
void alpha();
double dbrent1(double ax, double bx, double cx, double (*f)(int,int,double), double tol,
        double *xmin,int site,int root);
void worker();
void manager(int numprocs);
int main_function(int argc, char **argv);
char *produceTree(int nodeNo);

//test
void find_dist();

//Global Variables
NODE *node;
BRANCH *branch;
double **S_ij,**Q_ij,*PI,**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;
topology **copyTree;
int branchNo=0; //used for assignBranch
double *alfa;
//double **idMatrix;
int *gapNumber;
double gapForEst;
int D_A_P_Flag,Anneal_site_func1;
int cutIndex;
long idum=(-64);
//long idum=(-3);
bool rootOnly=false,treefileFlag=false,bioRoot=false;
int refineFlag=1, conditionalProb=0;
char seqFile[50], treeFile[50], alignFile[50];
bool treeOnlyFlag=false;


void usage()
{
	printf("Usage: ancestral [-O/-C] [-D] [-R/-RO [-Z]] [-G] -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:  Use Dan method to find the optimal solution\n");
	printf("       -R:  Reconstruct sequence for biological root and all internal nodes\n");
	printf("       -RO: Only reconstruct sequence for the biological root\n");
	printf("       -G:  Non-gap number or percentage (>=3 or =0~1), no estimated alpha used at that site if non-gap number is less\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 numprocs, myrank;
	bool treefileFlag;
	MPI_Init(&argc, &argv);
	MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
	MPI_Comm_rank(MPI_COMM_WORLD, &myrank);

	main_function(argc, argv);

	if(myrank!=0) worker();
	else manager(numprocs);

	MPI_Finalize();
	return 0;
}

void manager(int numprocs)
{
	int i,j,k;
	int sender, numsent=0;
	int send;
	MPI_Status status;
	int namelen; char processor_name[MPI_MAX_PROCESSOR_NAME];
	int myrank;
	MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
	MPI_Get_processor_name(processor_name, &namelen);

	if(!treefileFlag){
		FILE *fp;
		double **dist_matrix;
		dist_matrix=dmatrix(1,seqNum,1,seqNum);
		double *receive;
		receive=dvector(0,seqNum);
		fp=fopen("DISTMATRIX_mpi.txt","w");
		for(i=1;i<=seqNum;i++) for(j=1;j<=seqNum;j++) dist_matrix[i][j]=0.0;
		for(i=1;i<min(seqNum,numprocs);i++){
			send=i;
			MPI_Send(&send,1,MPI_INT,i,i,MPI_COMM_WORLD);
//			printf("server: %s send message for seqNum %d\n",processor_name,send);
		}
		numsent=min(seqNum,numprocs);
		
		for(i=1;i<=seqNum-1;i++) {
			MPI_Recv(receive,seqNum+1,MPI_DOUBLE, MPI_ANY_SOURCE, MPI_ANY_TAG,MPI_COMM_WORLD, &status);
			sender=status.MPI_SOURCE;
			int site=status.MPI_TAG;
			for(j=1;j<=seqNum;j++) dist_matrix[site][j]=receive[j];
			if(numsent<=seqNum-1){
				send=numsent;
				MPI_Send(&send,1,MPI_INT,sender,numsent,MPI_COMM_WORLD);
				numsent++;
			}
			else
				MPI_Send(MPI_BOTTOM,0,MPI_INT,sender,0,MPI_COMM_WORLD);
		}

		//symmetrize dist_matrix and output DISTMATRIX.txt
		for(i=1;i<seqNum;i++) for(j=i+1;j<=seqNum;j++) dist_matrix[j][i]=dist_matrix[i][j];
		for(i=1;i<=seqNum;i++) dist_matrix[i][i]=0.0;
		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);
		//build tree
		char outputTreeFile[50];
		strcpy(outputTreeFile,seqFile);
		for(i=strlen(outputTreeFile);i>=0;i--) if(outputTreeFile[i]=='.') { outputTreeFile[i]='\0'; break; }
		strcat(outputTreeFile,".tre");
		weighbor(seqNum,seqLength,node,dist_matrix,outputTreeFile);
		setLevel(seqNum+1);

		//define the root flag
		for(i=1;i<=seqNum*2-2;i++){
			if(i==seqNum+1) node[i].rootFlag=true;
			else node[i].rootFlag=false;
		}
		//define new date type
		int *treeData;
		treeData=new int [6*(seqNum*2-2)+1];
		for(i=1;i<=seqNum*2-2;i++){
			treeData[(i-1)*6+1]=node[i].childL;
			treeData[(i-1)*6+2]=node[i].childR;
			if(node[i].rootFlag) treeData[(i-1)*6+3]=1;
			else treeData[(i-1)*6+3]=0;
			treeData[(i-1)*6+4]=node[i].childE;
			treeData[(i-1)*6+5]=node[i].parent;
			treeData[(i-1)*6+6]=(int) (node[i].dist*pow(10.0,8));
		}
		//send tree information to all workers
		for(i=1;i<numprocs;i++) MPI_Send(treeData,(seqNum*2-2)*6+1,MPI_INT,i,i,MPI_COMM_WORLD);
	}
	else {
		readTree(treeFile);
		setLevel(seqNum+1);

		//define the root flag
		for(i=1;i<=seqNum*2-2;i++){
			if(i==seqNum+1) node[i].rootFlag=true;
			else node[i].rootFlag=false;
		}
		//define new date type
		int *treeData;
		treeData=new int [6*(seqNum*2-2)+1];
		for(i=1;i<=seqNum*2-2;i++){
			treeData[(i-1)*6+1]=node[i].childL;
			treeData[(i-1)*6+2]=node[i].childR;
			if(node[i].rootFlag) treeData[(i-1)*6+3]=1;
			else treeData[(i-1)*6+3]=0;
			treeData[(i-1)*6+4]=node[i].childE;
			treeData[(i-1)*6+5]=node[i].parent;
			treeData[(i-1)*6+6]=(int) (node[i].dist*pow(10.0,8));
		}
		//send tree information to all workers
		for(i=1;i<numprocs;i++) MPI_Send(treeData,(seqNum*2-2)*6+1,MPI_INT,i,i,MPI_COMM_WORLD);

	}

	assignTree();
	branchNo=0;
	assignBranch(seqNum+1);
/*
	for(i=1;i<=seqLength;i++) {
		printf("site:%4d root:%4d ",i,treeRoot[i]);
		if(gapTree[i][treeRoot[i]].gapFlag) printf("gapFlag: 1\n");
		else printf("gapFlag: 0\n");
	}
	fflush(stdout);
*/
//	optimalAlfa(refineFlag,conditionalProb);
	if(refineFlag==0){
		double receiveAlpha;
		double *optimizedAlpha; optimizedAlpha=new double [seqLength+1];
		for(i=0;i<=seqLength;i++) optimizedAlpha[i]=0.0;
		for(i=1;i<min(seqLength, numprocs);i++){
			send=i;
			MPI_Send(&send,1,MPI_INT,i,i,MPI_COMM_WORLD);
			printf("server: %s send message for site %d\n",processor_name,send);
		}
		numsent=min(seqLength,numprocs);
		
		for(i=1;i<=seqLength;i++){ 
			MPI_Recv(&receiveAlpha,1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
			sender=status.MPI_SOURCE;
			int site=status.MPI_TAG;
			printf("site: %3d original alpha: %10.6f optimized alpha: %10.6f logL: %10.6f\n",site,alfa[site],receiveAlpha,calTree(site,treeRoot[site],receiveAlpha));
//			printf("site: %3d original alpha: %10.6f optimized alpha: %10.6f\n",site,alfa[site],receiveAlpha);
			fflush(stdout);
			optimizedAlpha[site]=receiveAlpha;
			if(numsent<=seqLength) {
				send=numsent;
				MPI_Send(&send,1,MPI_INT,sender,numsent,MPI_COMM_WORLD);
				numsent++;
			}
			else 
				MPI_Send(MPI_BOTTOM,0,MPI_INT,sender,0,MPI_COMM_WORLD);
		}

		//remove outliers, not done yet
		bool *estFlag; estFlag=new bool [seqLength+1]; 
		for(i=1;i<=seqLength;i++){
			estFlag[i]=true;
			if(fabs(optimizedAlpha[i]-alfa[i])<1.0e-8 || alfa[i]<1.0e-10) estFlag[i]=false;
		}
		double *amplier, mean, sigma, *zScore;
		amplier=dvector(1,seqLength);
		zScore=dvector(1,seqLength);
		int length=0;
		for(i=1;i<=seqLength;i++) {
			amplier[i]=1.0; zScore[i]=0.0;
			if(estFlag[i]) { amplier[i]=optimizedAlpha[i]/alfa[i]; length++; }
		}
		bool outlierFlag=true;
		int loop=0;
		while(outlierFlag) {
			loop++;
			outlierFlag=false;
			mean=0.0; sigma=0.0;
			for(i=1;i<=seqLength;i++) if(estFlag[i]) mean+=amplier[i]/(1.0*length);
			for(i=1;i<=seqLength;i++) if(estFlag[i]) sigma+=(amplier[i]-mean)*(amplier[i]-mean);
			sigma=sqrt(sigma/(length-1));
			for(i=1;i<=seqLength;i++) if(estFlag[i]) {
				zScore[i]=(amplier[i]-mean)/sigma;
				if(zScore[i]>3.0) {
					amplier[i]=3.0*sigma+mean;
					outlierFlag=true;
				}
			}
			if(loop>20) outlierFlag=false;
		}
		fprintf(resultFile, "result from optimizing alpha\n");
		for(i=1;i<=seqLength;i++) 
			fprintf(resultFile,"site: %3d original alpha: %10.6f ML alpha: %10.6f final alpha: %10.6f\n",i,alfa[i],optimizedAlpha[i],amplier[i]*alfa[i]);

		for(i=1;i<=seqLength;i++) if(estFlag[i]) optimizedAlpha[i]=alfa[i]*amplier[i];
		//finish remove outliers

		for(i=1;i<numprocs;i++) {
			MPI_Send(optimizedAlpha,seqLength+1,MPI_DOUBLE,i,i,MPI_COMM_WORLD);
		}

		fflush(resultFile);
		for(i=1;i<=seqLength;i++) alfa[i]=optimizedAlpha[i];
	}


		for(i=1;i<=seqLength;i++) 
			fprintf(resultFile,"site: %3d original alpha: %10.6f logL: %10.6f\n",i,alfa[i],calTree(i,treeRoot[i],alfa[i]));
		for(i=1;i<=seqLength;i++) {
			if(fabs(alfa[i]-1.0)<1.0e-8) fprintf(resultFile,"%12.4f ",0.0);
			else fprintf(resultFile,"%12.4f ",calTree(i,treeRoot[i],alfa[i]));
		}
		fprintf(resultFile,"\n");
		for(i=1;i<=seqLength;i++) fprintf(resultFile,"%10.4f ",alfa[i]);
		fprintf(resultFile,"\n");
		fflush(resultFile);

	if(refineFlag==1) { for(i=1;i<=seqLength;i++) alfa[i]=1.0; }

	if(!rootOnly){
		for(i=1;i<min(seqLength,numprocs);i++){
			send=i;
			MPI_Send(&send,1,MPI_INT,i,i,MPI_COMM_WORLD);
//			printf("server: %s send message for seqNum %d\n",processor_name,send);
		}
		numsent=min(seqLength,numprocs);
		
//		char *tmpChar;
//		tmpChar = new char [seqNum-1];
		double *dataStream;
		dataStream=new double [0,seqNum-1+(seqNum-2)*NP];
		for(i=1;i<=seqLength;i++) {
//			MPI_Recv(tmpChar,seqNum-1,MPI_DOUBLE, MPI_ANY_SOURCE, MPI_ANY_TAG,MPI_COMM_WORLD, &status);
			MPI_Recv(dataStream,seqNum-1+(seqNum-2)*NP,MPI_DOUBLE, MPI_ANY_SOURCE, MPI_ANY_TAG,MPI_COMM_WORLD, &status);
			sender=status.MPI_SOURCE;
			int site=status.MPI_TAG;
//			for(j=1;j<=seqNum-2;j++) node[j+seqNum].seq_aa[site]=tmpChar[j];
			for(j=1;j<=seqNum-2;j++) node[j+seqNum].seq_aa[site]=aa[(int) dataStream[j]];
			for(j=1;j<=seqNum-2;j++) {
				for(k=1;k<=NP;k++) {
					node[j+seqNum].prob[site][k]=dataStream[seqNum-2+(j-1)*NP+k];
				}
			}
			if(numsent<=seqLength){
				send=numsent;
				MPI_Send(&send,1,MPI_INT,sender,numsent,MPI_COMM_WORLD);
				numsent++;
			}
			else
				MPI_Send(MPI_BOTTOM,0,MPI_INT,sender,0,MPI_COMM_WORLD);
		}
	}
	if(bioRoot) findMidpoint();

	for(k=seqNum+1;k<=seqNum*2-2;k++) {
        	fprintf(resultFile,"\n\n\nThe transition probability:\n");
        	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 reconstructed ancestral sequences
	fprintf(resultFile,"\n\nThe prediction for all ancestral nodes:\n");
	if(!rootOnly){
		for(i=1;i<=(seqNum-1)*2;i++){
			fprintf(resultFile,"%5d %20s %s\n",i,node[i].nodeName,node[i].seq_aa);
		}
	}
	if(bioRoot) fprintf(resultFile,"%5d %20s %s\n",i,node[0].nodeName,node[0].seq_aa);

	print_tree(seqNum+1);
	fclose(resultFile);

	if(strlen(alignFile)>3) {
		FILE *align;
		align=fopen(alignFile,"w");
		if(!align) { printf("can not open file: %s\n",alignFile); exit(1); }
		if(!rootOnly){
			for(i=seqNum+1;i<=(seqNum-1)*2;i++){
				fprintf(align,"%-20s %s\n",node[i].nodeName,node[i].seq_aa);
			}
		}
		if(bioRoot) fprintf(align,"%-20s %s\n",node[0].nodeName,node[0].seq_aa);
	}
}

void worker()
{
	int i,j,k;
	int myrank;
	MPI_Status status;
	int send;
	int namelen; char processor_name[MPI_MAX_PROCESSOR_NAME];
	MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
	MPI_Get_processor_name(processor_name, &namelen);
	printf("processor: %s rank: %d is running\n",processor_name,myrank);
	if(!treefileFlag){
		double *receive;
		receive=dvector(0,seqNum);
		int **compare;
		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;
				}
			}
		}
		for(i=2;i<=seqNum;i++) for(j=1;j<i;j++) compare[j][i]=compare[i][j];
		for(i=1;i<=seqNum;i++) compare[i][i]=0;
		if(myrank<=seqNum-1){
			MPI_Recv(&send,1,MPI_INT,0,MPI_ANY_TAG, MPI_COMM_WORLD, &status);
			while(status.MPI_TAG>0) {
				for(i=1;i<=seqNum;i++){
					if(i<=send) receive[i]=0.0;
					else if(compare[send][i]==1) receive[i]=0.0;
					else receive[i]=dist(send,i);
				}
				MPI_Send(receive,seqNum+1,MPI_DOUBLE,0,send,MPI_COMM_WORLD);
				MPI_Recv(&send,1,MPI_INT,0,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
			}
		}

		//receive tree information
		int *treeData;
		treeData=new int [6*(seqNum*2-2)+1];
		MPI_Recv(treeData,(seqNum*2-2)*6+1,MPI_INT,0,MPI_ANY_TAG,MPI_COMM_WORLD,&status);

		//assign tree structure
		for(i=1;i<=seqNum*2-2;i++){
			node[i].childL=treeData[(i-1)*6+1];
			node[i].childR=treeData[(i-1)*6+2];
			if(treeData[(i-1)*6+3]==1) node[i].rootFlag=true;
			else node[i].rootFlag=false;
			node[i].childE=treeData[(i-1)*6+4];
			node[i].parent=treeData[(i-1)*6+5];
			node[i].dist=treeData[(i-1)*6+6]/pow(10.0,8);
		}

		setLevel(seqNum+1);

	}
//	else readTree(treeFile);
	else {
		int *treeData;
		treeData=new int [6*(seqNum*2-2)+1];
		MPI_Recv(treeData,(seqNum*2-2)*6+1,MPI_INT,0,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
		//assign tree structure
		for(i=1;i<=seqNum*2-2;i++){
			node[i].childL=treeData[(i-1)*6+1];
			node[i].childR=treeData[(i-1)*6+2];
			if(treeData[(i-1)*6+3]==1) node[i].rootFlag=true;
			else node[i].rootFlag=false;
			node[i].childE=treeData[(i-1)*6+4];
			node[i].parent=treeData[(i-1)*6+5];
			node[i].dist=treeData[(i-1)*6+6]/pow(10.0,8);
		}
		setLevel(seqNum+1);
	}


	assignTree();
	branchNo=0;
	assignBranch(seqNum+1);

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

	if(refineFlag==0){
		double sendAlpha;
		double *changedAlpha; changedAlpha=new double [seqLength+1];
		if(myrank<=seqLength) {
			MPI_Recv(&send,1,MPI_INT,0,MPI_ANY_TAG,MPI_COMM_WORLD, &status);
			while(status.MPI_TAG>0){
				int row=send;
				sendAlpha=optimizeAlpha(send);
				MPI_Send(&sendAlpha,1,MPI_DOUBLE,0,row,MPI_COMM_WORLD);
				MPI_Recv(&send,1,MPI_INT,0,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
			}
		}
		MPI_Recv(changedAlpha,seqLength+1,MPI_DOUBLE,0,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
		for(j=1;j<=seqLength;j++) alfa[j]=changedAlpha[j];
//		if(myrank==1) for(j=1;j<=seqLength;j++) fprintf(resultFile,"rank=1 alfa: %10.6f\n",alfa[j]);
//		fflush(resultFile);
	}

	if(refineFlag==1) { for(i=1;i<=seqLength;i++) alfa[i]=1.0; }

	if(!rootOnly){
		if(myrank<=seqLength){
//			char *tmpChar;
//			tmpChar = new char [seqNum-1];
                	double *dataStream;
                	dataStream=new double [0,seqNum-1+(seqNum-2)*NP];
			MPI_Recv(&send,1,MPI_INT,0,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
			while(status.MPI_TAG>0){
				if((!dashFlag[i] || calFlag[i]) && (conditionalProb==1)) optTree(send,alfa[send]);
				if((!dashFlag[i] || calFlag[i]) && (conditionalProb==0)) optTree_Goldstein(send,alfa[send]);
/*				for(i=seqNum+1;i<=seqNum*2-2;i++) {
					tmpChar[i-seqNum]=aa[node[i].seq_num[send]];
				}
*/
				for(i=seqNum+1;i<=seqNum*2-2;i++) dataStream[i-seqNum]=node[i].seq_num[send];
				for(i=seqNum+1;i<=seqNum*2-2;i++) for(j=1;j<=NP;j++) dataStream[seqNum-2+(i-seqNum-1)*NP+j]=node[i].prob[send][j];
//				MPI_Send(tmpChar,seqNum-1,MPI_CHAR,0,send,MPI_COMM_WORLD);
				MPI_Send(dataStream,seqNum-1+(seqNum-2)*NP,MPI_DOUBLE,0,send,MPI_COMM_WORLD);
				MPI_Recv(&send,1,MPI_INT,0,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
			}
		}
	}
}


int main_function(int argc, char **argv)
{
//	bool rootOnly=false,treefileFlag=false,bioRoot=false;
	int i;
//	int refineFlag=1, conditionalProb=0;
//	bool flag=true;
	FILE *input; int number;
//	char seqFile[50], treeFile[50];
	strcpy(alignFile,"");
	unsigned int calTime1, calTime2;
	strcpy(treeFile,"DISTMATRIX.tre");
	gapForEst=0.5;
	//initiate the matrix and vectors
	S_ij=dmatrix(1,NP,1,NP);
	Q_ij=dmatrix(1,NP,1,NP);
	PI=dvector(1,NP);
	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 'G':
				++i;
				optidx=0;
				gapForEst=atof(argv[i]);
//				printf("%20.16f\n",gapForEst);
				if(gapForEst<0.0 || (gapForEst>1.0000000001 && gapForEst<2.999999999)) usage();
				break;
			case 'C':
				refineFlag=2;
				break;
			case 'D':
				conditionalProb=1;
				break;
			case 'R':
				bioRoot=true;
				if(argv[i][2]=='O') rootOnly=true;
				break;
			case 'T':
				if(argv[i][2]=='O') treeOnlyFlag=true;
				else usage();
				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 'a':
				++i;
				optidx=0;
				if(i<argc) {
					strcpy(alignFile,argv[i]);
					input=fopen(alignFile,"w");
					if(!input) {
                                                printf("can not open file:%s\n",alignFile);
                                                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(PI);//calculate PI vector and Q(ij) matrix
	cal_Eign();//calculate eignvalue vector and eignvector matrix of Q(ij)
	alpha(); //calculate Alignment Based Alpha (AB alpha)
//	if(!treefileFlag) cal_dist();
/*	if(!treefileFlag) {
		cal_dist(seqFile); 
	}
	else readTree(treeFile);
    assignTree();//assign tree structure to gapTree
//	printf("Finish assign tree\n");

	branchNo=0;
	assignBranch(seqNum+1);
*/
//	printf("Finish assign branch\n");
	calTime2=time(NULL);
//	printf("Prepare time: %ds\n", calTime2-calTime1);
	calTime1=calTime2;
//	if(!rootOnly || refineFlag==2) optimalAlfa(refineFlag,conditionalProb,rootOnly);
//	if(!rootOnly) optimalAlfa(refineFlag,conditionalProb,rootOnly);
/*	optimalAlfa(refineFlag,conditionalProb,rootOnly);
	if(conditionalProb==0){
		if(bioRoot){
			findMidpoint();
			if(ZFlag){
//				Z_score1();
				Z_score2(refineSeqItem,cutDistItem,predFile,seqFile);
	//			searchFunctionalSite();
			}
		}
	}
	
	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);

	//free memory

	free_dmatrix(S_ij,1,NP,1,NP);
//	for(i=1;i<=NP;i++) { printf("%3d ",i); for(int j=1;j<=NP;j++) printf("%8.3f ",Q_ij[i][j]); printf("\n"); }
	free_dmatrix(Q_ij,1,NP,1,NP);
//	for(int j=1;j<=NP;j++) printf("%8.3f ",PI[j]); printf("\n\n");
	free_dvector(PI,1,NP);
//	for(i=1;i<=NP;i++) { printf("%3d ",i); for(int j=1;j<=NP;j++) printf("%8.3f ",a[i][j]); printf("\n"); }
	free_dmatrix(a,1,NP,1,NP);
	free_dmatrix(v,1,NP,1,NP);
	free_dmatrix(inv,1,NP,1,NP);
	free_dvector(wr,1,NP);
	free_dvector(wi,1,NP);
	free_dvector(w,1,NP);
	free_dmatrix(svdv,1,NP,1,NP);
	free_dvector(colu,1,NP);
*/
}

/*
void find_dist()
{
	int i,j;
	double **dist;
	dist=dmatrix(1,seqNum,1,seqNum);
	cout<<"in find_dist\n";
	for(i=1;i<seqNum;i++){
		changeRoot(0,i,gapTree[0][i].parent,i);
//		cout<<"ok for root change"<<endl;
		for(j=i+1;j<=seqNum;j++){
			double dist1=0.0;
			int parent=j;
			while(!gapTree[0][parent].rootFlag){
				dist1+=gapTree[0][parent].dist;
				parent=gapTree[0][parent].parent;
			}
			dist[i][j]=dist1+gapTree[0][i].dist;
//			cout<<dist[i][j]<<endl;
		}
	}
	for(i=1;i<seqNum;i++) for(j=i+1;j<=seqNum;j++) {
//		printf("original dist: %10.5f\n",dist[i][j]);
		dist[i][j]=(int) (dist[i][j]*seqLength);
//		printf("multiple seqLength dist: %10.5f\n",dist[i][j]);
		dist[i][j]/=seqLength;
//		printf("final dist: %10.5f\n",dist[i][j]);
//		getchar();
	}
		// dist[i][j]= 1.0*((int) dist[i][j]*seqLength)/seqLength;
	for(i=1;i<=seqNum;i++) dist[i][i]=0.0;
	for(i=1;i<seqNum;i++) for(j=i+1;j<=seqNum;j++) dist[j][i]=dist[i][j];
//	weighbor(seqNum,seqLength,node,dist,"test.tre");
	for(i=1;i<seqNum;i++) for(j=i+1;j<=seqNum;j++) fprintf(resultFile,"%10.5f\n",dist[i][j]);
}
*/

void assignBranch(int nodeNo)
{
	if(node[nodeNo].rootFlag){
		branchNo++;
		branch[branchNo].length=node[node[nodeNo].childL].dist;
		branch[branchNo].LPoint=nodeNo;
		branch[branchNo].RPoint=node[nodeNo].childL;
		assignBranch(node[nodeNo].childL);

		branchNo++;
		branch[branchNo].length=node[node[nodeNo].childR].dist;
		branch[branchNo].LPoint=nodeNo;
		branch[branchNo].RPoint=node[nodeNo].childR;
		assignBranch(node[nodeNo].childR);

		branchNo++;
		branch[branchNo].length=node[node[nodeNo].childE].dist;
		branch[branchNo].LPoint=nodeNo;
		branch[branchNo].RPoint=node[nodeNo].childE;
		assignBranch(node[nodeNo].childE);
	}
	else if((node[nodeNo].childL!=0)&&(node[nodeNo].childR!=0)){
		branchNo++;
		branch[branchNo].length=node[node[nodeNo].childL].dist;
		branch[branchNo].LPoint=nodeNo;
		branch[branchNo].RPoint=node[nodeNo].childL;
		assignBranch(node[nodeNo].childL);

		branchNo++;
		branch[branchNo].length=node[node[nodeNo].childR].dist;
		branch[branchNo].LPoint=nodeNo;
		branch[branchNo].RPoint=node[nodeNo].childR;
		assignBranch(node[nodeNo].childR);
	}
}

void findMidpoint()
{
	double ran1(long *idum);
	int i,j;
	double **length1,**length2,**length3;
	double *mean,*x,*y;
	double a,b,lengthx;
	int parent,indexY,valueY;
//	topology **copyTree;
//	copyTree = topmatrix(0,1,0,seqNum*2-2);
//	copyTree = new topology * [2];
//	for(i=0;i<2;i++) copyTree[i] =new topology [seqNum*2-1];

	length1=dmatrix(1,2,1,seqNum);
	length2=dmatrix(1,2,1,seqNum);
	length3=dmatrix(1,2,1,seqNum);
	mean=dvector(1,2);
	x=dvector(1,seqNum*2-3);
	y=dvector(1,seqNum*2-3);
	//count each internal node has how many leaf node as his descedent
	int *countNode;
	countNode=ivector(1,seqNum*2-2);
	for(i=1;i<=seqNum*2-2;i++) countNode[i]=0;

	//find out all LPointNumber and RPointNumber for each branch
	for(i=1;i<=seqNum;i++){
		parent=i;
		while(!gapTree[0][parent].rootFlag){
			countNode[parent]++;
			parent=gapTree[0][parent].parent;
		}
	}

//	for(i=1;i<=seqNum*2-2;i++) printf("%4d %4d\n",i,countNode[i]);

	for(i=1;i<=seqNum*2-3;i++){
		branch[i].RPointNumber=countNode[branch[i].RPoint];
		branch[i].LPointNumber=seqNum-branch[i].RPointNumber;
	}
	//test 
	
	for(i=1;i<=seqNum*2-3;i++){
		x[i]=0.0;y[i]=0.0;
		changeRoot(0,branch[i].LPoint,branch[i].LPoint,branch[i].LPoint);
		for(j=1;j<=seqNum;j++){
			parent=j; lengthx=0.0;
			length1[2][j]=1.0;
			while(!gapTree[0][parent].rootFlag){
				lengthx+=gapTree[0][parent].dist;
				if(parent==branch[i].RPoint) length1[2][j]=-1.0;
				parent=gapTree[0][parent].parent;
			}
			length1[1][j]=lengthx;
//			printf("leaf node:%3d a:%8.5f b:%4.1f\n",j,length1[1][j],length1[2][j]);
		}
		mean[1]=mean[2]=0.0;
		for(j=1;j<=seqNum;j++) {
			mean[1]+=length1[1][j]/(1.0*seqNum);
			mean[2]+=length1[2][j]/(1.0*seqNum);
		}
//		printf("mean1:%8.5f mean2:%8.5f\n",mean[1],mean[2]);
		a=0.0; b=0.0;
		for(j=1;j<=seqNum;j++){
			length2[1][j]=length1[1][j]-mean[1];
			length2[2][j]=length1[2][j]-mean[2];
			length3[1][j]=2.0*length2[1][j]*length2[2][j];
			length3[2][j]=2.0*length2[2][j]*length2[2][j];
			a+=length3[1][j]; b+=length3[2][j];
//			printf("leaf node:%3d a:%8.5f b:%8.5f  ",j,length2[1][j],length2[2][j]);
//			printf("leaf node:%3d a:%8.5f b:%8.5f\n",j,length3[1][j],length3[2][j]);
		}
//		printf("a:%8.5f b:%8.5f\n",a,b);
		x[i]=-1.0*a/b;
		if((x[i]<0.0)||(x[i]>branch[i].length)) y[i]=1.0e6;
		else{
			for(j=1;j<=seqNum;j++) 
				y[i]+=(length2[1][j]+length2[2][j]*x[i])*(length2[1][j]+length2[2][j]*x[i]);
		}
		if(y[i]<1.0e6){
			printf("branch:%3d x value:%10.6f variance:%10.6f\t LeftLeaf number:%4d RightLeaf number:%4d ",i,x[i],y[i],branch[i].LPointNumber,branch[i].RPointNumber);
			printf("difference: %4d\n",abs(branch[i].LPointNumber-branch[i].RPointNumber));
		}
	}

//	valueY=y[1];indexY=1;
/*	valueY=10000.0;indexY=0;
	for(i=1;i<=seqNum*2-3;i++){
		if(y[i]<valueY) { valueY=y[i]; indexY=i; }
	}
	if(indexY==0) { printf("there is something wrong with the root find\n"); exit(0);}
*/
	//find the root with balance leaf node at each side
	valueY=seqNum*2; indexY=0;
	for(i=1;i<=seqNum*2-3;i++){
		if(y[i]<1.0e5){
			if(abs(branch[i].LPointNumber-branch[i].RPointNumber)<valueY){
				valueY=abs(branch[i].LPointNumber-branch[i].RPointNumber);
				indexY=i;
			}
			else if(abs(branch[i].LPointNumber-branch[i].RPointNumber)==valueY){
				if(y[i]<y[indexY]) indexY=i;
			}
		}
	}

	fprintf(resultFile,"best point: branch No:%3d distance from left point:%10.6f\n",indexY,x[indexY]);
	fprintf(resultFile,"left point: %d: %s ",branch[indexY].LPoint,node[branch[indexY].LPoint].nodeName);
	fprintf(resultFile,"right point: %d: %s\n",branch[indexY].RPoint,node[branch[indexY].RPoint].nodeName); 
	free_dmatrix(length1,1,2,1,seqNum);
	free_dmatrix(length2,1,2,1,seqNum);
	free_dmatrix(length3,1,2,1,seqNum);
	free_dvector(mean,1,2);
	//add root node to the tree
	changeRoot(0,branch[indexY].LPoint,branch[indexY].LPoint,0);
	if(branch[indexY].RPoint==gapTree[0][branch[indexY].LPoint].childL){
		gapTree[0][branch[indexY].LPoint].childL=gapTree[0][branch[indexY].LPoint].childE;
		gapTree[0][branch[indexY].LPoint].childE=branch[indexY].RPoint;
	}
	else if(branch[indexY].RPoint==gapTree[0][branch[indexY].LPoint].childR){
		gapTree[0][branch[indexY].LPoint].childR=gapTree[0][branch[indexY].LPoint].childE;
		gapTree[0][branch[indexY].LPoint].childE=branch[indexY].RPoint;
	}

//	double inputV;
//	cout<<"input:";
//	cin>>inputV;
	gapTree[0][branch[indexY].LPoint].parent=0;
	gapTree[0][branch[indexY].LPoint].childE=0;
	gapTree[0][branch[indexY].LPoint].rootFlag=false;
	gapTree[0][branch[indexY].LPoint].dist=x[indexY];
//	gapTree[0][branch[indexY].LPoint].dist=inputV;

	gapTree[0][branch[indexY].RPoint].parent=0;
	gapTree[0][branch[indexY].RPoint].dist=branch[indexY].length-x[indexY];
//	gapTree[0][branch[indexY].RPoint].dist=branch[indexY].length-inputV;

	gapTree[0][0].parent=0;
	gapTree[0][0].childL=branch[indexY].LPoint;
	gapTree[0][0].childR=branch[indexY].RPoint;
	gapTree[0][0].childE=0;
	gapTree[0][0].rootFlag=true;

	free_dvector(x,1,seqNum*2-3);
	free_dvector(y,1,seqNum*2-3);

	//output the rooted tree to output file
	fprintf(resultFile,"Rooted Tree: %s\n",produceTree(0));

	//calculate the mean distance from root to all leaf nodes
	double mean_root_leaf=0.0, *dist_root_leaf;
	dist_root_leaf=dvector(1,seqNum);
	for(i=1;i<=seqNum;i++) dist_root_leaf[i]=0.0;
	for(i=1;i<=seqNum;i++) {
		int parent = i;
		while(!gapTree[0][parent].rootFlag) {
			dist_root_leaf[i]+=gapTree[0][parent].dist;
			parent=gapTree[0][parent].parent;
		}
		mean_root_leaf+=dist_root_leaf[i]/seqNum;
//		printf("%10.5f\n",dist_root_leaf[i]);
	}
	printf("The mean distance from root to leaf nodes: %10.5f\n",mean_root_leaf);
	fprintf(resultFile,"The mean distance from root to leaf nodes: %10.5f\n",mean_root_leaf);

	if(treeOnlyFlag) return;
	

	for(i=1;i<=seqLength;i++){
		//copy gapTree[0] to copyTree
		for(j=0;j<=seqNum*2-2;j++) copyTree[0][j]=gapTree[0][j];
		//set gapFlag for each leaf node
		for(j=1;j<=seqNum;j++) if(node[j].seq_aa[i]=='-') gapTree[0][j].gapFlag=true;
		//rebuild tree
		reconstructTree(0,0,true);
		if(!gapTree[0][0].gapFlag){
			if(D_A_P_Flag==4){
				for(j=1;j<=NP;j++) PI[j]=sitePI[i][j];
//				printf("%3d ",i);
//				for(j=1;j<=NP;j++) printf("%6.0f ",PI[j]*100000.0);
//				printf("\n");
				cal_Q(PI);
				cal_Eign();
			}

			//calculate L vector for root node
			calRL(0,gapTree[0][0].childL,i,alfa[i]);
			calRL(0,gapTree[0][0].childR,i,alfa[i]);
			for(int k=1;k<=NP;k++){
				gapTree[0][0].L[k]=PI[k]*gapTree[0][gapTree[0][0].childL].L[k]*gapTree[0][gapTree[0][0].childR].L[k];
			}
			double tmpsum=0.0;
			for(int k=1;k<=NP;k++) tmpsum+=gapTree[0][0].L[k];
			for(int k=1;k<=NP;k++) node[0].prob[i][k]=gapTree[0][0].L[k]/tmpsum;
			double tmpMax; int indexMax=1;
			tmpMax=node[0].prob[i][1];
//			tmpMax=-1.0; indexMax=0;
			for(int k=2;k<=NP;k++){
				if(node[0].prob[i][k]>tmpMax) {
					tmpMax=node[0].prob[i][k];
					indexMax=k;
				}
			}
//			if(indexMax==0) printf("Warning: something wrong in find optimal root\n");
			for(int k=1;k<=NP;k++)
				if(node[0].prob[i][k]<0.00000001) node[0].prob[i][k]=0.0;
			node[0].seq_num[i]=indexMax;
			node[0].seq_aa[i]=aa[indexMax];
			//used for second best
			double tmpSec=-1.0; int indexSec=0;
			for(int k=1;k<=NP;k++) {
				if(node[0].prob[i][k]>tmpSec && k!=indexMax) {tmpSec=node[0].prob[i][k]; indexSec=k; }
			}
			if(node[0].prob[i][indexSec]<1.0e-10) indexSec=(int) (ran1(&idum)*20+1);
//			secondBest[0].seq_num[i]=indexSec;
//			secondBest[0].seq_aa[i]=aa[indexSec];
		}
		else {
			node[0].seq_num[i]=0;
			node[0].seq_aa[i]=aa[0];
//			secondBest[0].seq_num[i]=0; secondBest[0].seq_aa[i]=aa[0];
		}
		//copy gapTree[0] to copyTree
		for(j=0;j<=seqNum*2-2;j++) gapTree[0][j]=copyTree[0][j];
	}
	//print out result
	//print the optimal solution
	fprintf(resultFile,"\n\nThe prediction for midpoint:\n");
	fprintf(resultFile,"%5d %20s %s\n",0,"Root-MidPoint",node[0].seq_aa);

	fprintf(resultFile,"\n\n\nThe transition probability:\n");
		fprintf(resultFile,"Sequence Number: %4d  (probability*10^(-4))\n",0);
		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[0].seq_aa[i]);
			for(j=1;j<=NP;j++){
				fprintf(resultFile,"%6.0f",node[0].prob[i][j]*10000);
			}
			fprintf(resultFile,"\n");
		}
		fprintf(resultFile,"\n\n");

//	for(i=0;i<2;i++) delete copyTree[i];
//	delete copyTree;

	//test for root locate at different place on the branch near one internal node
/*		double **test;
		double tmpTotal=0.0;
		test=dmatrix(1,seqLength,1,NP);
		for(i=1;i<=seqLength;i++){
			for(j=1;j<=NP;j++){
				test[i][j]=node[0].prob[i][j]-node[branch[indexY].LPoint].prob[i][j];
			}
			//find the maximum for one site
			double tmpTest=test[i][1];
			for(j=2;j<=NP;j++) if(test[i][j]>tmpTest) tmpTest=test[i][j];
			tmpTotal+=tmpTest;
		}
		printf("%f %20.10f\n",inputV,tmpTotal);
*/	

}

int reconstructTree(int siteNo, int root, bool rootedFlag)
{
	int CCL,CCR,CCE=-1;
	CCL=gapTree[siteNo][root].childL;
	CCR=gapTree[siteNo][root].childR;
	if(!rootedFlag) CCE=gapTree[siteNo][root].childE;
	reconstruct(siteNo,CCL);
	reconstruct(siteNo,CCR);
	if(!rootedFlag) reconstruct(siteNo,CCE);
	if(rootedFlag){
		CCL=gapTree[siteNo][root].childL;
		CCR=gapTree[siteNo][root].childR;
		if((gapTree[siteNo][CCL].gapFlag)||(gapTree[siteNo][CCR].gapFlag))
			gapTree[siteNo][root].gapFlag=true;
	}
}

void reconstruct(int siteNo, int nodeNo)
{
	int CCL,CCR;
	int parent;
	if((gapTree[siteNo][nodeNo].childL!=0)&&(gapTree[siteNo][nodeNo].childR!=0)){
		if(gapTree[siteNo][nodeNo].childL>seqNum)
			reconstruct(siteNo,gapTree[siteNo][nodeNo].childL);
		if(gapTree[siteNo][nodeNo].childR>seqNum)
			reconstruct(siteNo,gapTree[siteNo][nodeNo].childR);
		CCL=gapTree[siteNo][nodeNo].childL;
		CCR=gapTree[siteNo][nodeNo].childR;

		if((gapTree[siteNo][CCL].gapFlag)&&(gapTree[siteNo][CCR].gapFlag))
		{ gapTree[siteNo][nodeNo].gapFlag=true; gapTree[siteNo][nodeNo].childL=gapTree[siteNo][nodeNo].childR=-1;}
		else if((!gapTree[siteNo][CCL].gapFlag)&&(gapTree[siteNo][CCR].gapFlag)){
			parent=gapTree[siteNo][nodeNo].parent;
			gapTree[siteNo][CCL].parent=parent;
			if(nodeNo==gapTree[siteNo][parent].childL) gapTree[siteNo][parent].childL=CCL;
			else if(nodeNo==gapTree[siteNo][parent].childR) gapTree[siteNo][parent].childR=CCL;
			else if(nodeNo==gapTree[siteNo][parent].childE) gapTree[siteNo][parent].childE=CCL;
			gapTree[siteNo][CCL].dist+=gapTree[siteNo][nodeNo].dist;
			gapTree[siteNo][nodeNo].gapFlag=true;
			gapTree[siteNo][nodeNo].childL=gapTree[siteNo][nodeNo].childR=-1;
		}
		else if((gapTree[siteNo][CCL].gapFlag)&&(!gapTree[siteNo][CCR].gapFlag)){
			parent=gapTree[siteNo][nodeNo].parent;
			gapTree[siteNo][CCR].parent=parent;
			if(nodeNo==gapTree[siteNo][parent].childL) gapTree[siteNo][parent].childL=CCR;
			else if(nodeNo==gapTree[siteNo][parent].childR) gapTree[siteNo][parent].childR=CCR;
			else if(nodeNo==gapTree[siteNo][parent].childE) gapTree[siteNo][parent].childE=CCR;
			gapTree[siteNo][CCR].dist+=gapTree[siteNo][nodeNo].dist;
			gapTree[siteNo][nodeNo].gapFlag=true;
			gapTree[siteNo][nodeNo].childL=gapTree[siteNo][nodeNo].childR=-1;
		}
	}
}

void calRL(int siteNo,int nodeNo, int leafSite, 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,leafSite,alfa);
		calRL(siteNo,gapTree[siteNo][nodeNo].childR,leafSite,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[leafSite],gapTree[siteNo][nodeNo].dist,alfa);
//			if((gapTree[siteNo][nodeNo].L[k]<0.0)&&(fabs(gapTree[siteNo][nodeNo].L[k])<1.0e-10))
//				gapTree[siteNo][nodeNo].L[k]=0.0;
			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);
		
	}
}


void assignTree()
{
	int i,j;
//	alfa= dvector(0,seqLength);
	for(i=1;i<=(seqNum*2-2);i++){
		for(j=0;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;
			for(j=0;j<=seqNum*2-2;j++) gapTree[i][j].gapFlag=true;
		}
	}
	//assign alfa vector
//	for(i=1;i<=seqLength;i++) alfa[i]=1.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++;}
						bool nameFlag=false;
						for(j=1;j<=index;j++) {
							if(tmpIndex==strlen(node[j].nodeName) && memcmp(node[j].nodeName,tmpName,tmpIndex)==0) {
//								printf("Yes\n");
								treeL[i]=node[j].nodeNum;
								fscanf(input,"%f",&tmpValue); 
								node[j].dist=tmpValue;
								nameFlag=true;
//								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");
						}
						if(!nameFlag) { printf("can not find the match name for: %s\n",tmpName);}
					}
					else {//3
						if(treeR[i]==0) {
							tmpIndex=0;
							while((c=fgetc(input))!=':') { tmpName[tmpIndex]=c; tmpIndex++;}
//							printf("i=%d\n",i);
							bool nameFlag=false;
							for(j=1;j<=index;j++) {
								if(tmpIndex==strlen(node[j].nodeName) && memcmp(node[j].nodeName,tmpName,tmpIndex)==0) {
									treeR[i]=node[j].nodeNum;
									fscanf(input,"%f",&tmpValue); node[j].dist=tmpValue;
									nameFlag=true;
//									printf("in tree Right: i=%d, j=%d, nodeNum=%d, dist=%f\n",i,j,node[j].nodeNum,node[j].dist);
									break;
								}
							}
							if(!nameFlag) { printf("can not find the match name for: %s\n",tmpName);}
						}
						else {//4
							tmpIndex=0;
							while((c=fgetc(input))!=':') { tmpName[tmpIndex]=c; tmpIndex++;}
							bool nameFlag=false;
							for(j=1;j<=index;j++) {
								if(tmpIndex==strlen(node[j].nodeName) && memcmp(node[j].nodeName,tmpName,tmpIndex)==0) {
									treeE=node[j].nodeNum;
									fscanf(input,"%f",&tmpValue); node[j].dist=tmpValue;
									nameFlag=true;
//									printf("in tree E: i=%d, j=%d, nodeNum=%d, dist=%f\n",i,j,node[j].nodeNum,node[j].dist);
									break;
								}
							}
							if(!nameFlag) { printf("can not find the match name for: %s\n",tmpName);}
						}//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:%4d nodeName:%10s level:%4d childL:%4d childR:%4d",node[i].nodeNum,node[i].nodeName,node[i].nodeLevel,node[i].childL,node[i].childR);
		if(i!=index)
			printf(" parent:%4d dist:%10.6f\n",node[i].parent,node[i].dist);
		else 
			printf(" childE:%4d parent:%4d dist:%10.6f\n",node[i].childE,node[i].parent,node[i].dist);
	}

	double multiply=0.0;
	printf("Please input: ");
	cin>>multiply;
	cout<<multiply<<endl;
	for(i=1;i<=seqNum*2-2;i++) node[i].dist*=multiply;
*/
	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(char *seqFile)
{
	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);

//	for(i=1;i<seqNum;i++) for(j=i+1;j<=seqNum;j++) fprintf(resultFile,"%10.5f\n",dist_matrix[i][j]);

/*	fprintf(stdout,"%d\n",seqNum);
	for(i=1;i<=seqNum;i++){
		fprintf(stdout,"%20s\t",node[i].nodeName);
		for(j=1;j<=seqNum;j++) fprintf(stdout,"%10.5f ",dist_matrix[i][j]);
		fprintf(stdout,"\n");
	}
*/
	//use weighbor method to build a tree
	printf("Use weighbor software to build tree\n");
	char outputTreeFile[50];
	strcpy(outputTreeFile,seqFile);
	for(i=strlen(outputTreeFile);i>=0;i--) if(outputTreeFile[i]=='.') outputTreeFile[i]='\0';
	strcat(outputTreeFile,".tre");
	weighbor(seqNum,seqLength,node,dist_matrix,outputTreeFile);
	//define the root flag
	for(i=1;i<=seqNum*2-2;i++) {
		if(i==seqNum+1) node[i].rootFlag=true;
		else node[i].rootFlag=false;
	}
	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)
double dist(int node1, int node2)
{
	int i,j,JMAX=60;
	int **Nab;
	double xLeft,xRight,dx,xMid,xacc=0.000001;
	double leftBracket,rightBracket,midBracket;
//	Nab=imatrix(0,NP,0,NP);
	Nab=imatrix(1,2,1,seqLength);
//	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<=seqLength;i++) Nab[1][i]=node[node1].seq_num[i];
	for(i=1;i<=seqLength;i++) Nab[2][i]=node[node2].seq_num[i];
/*	if(node1==seqNum && node2==1){
		for(i=1;i<=seqLength;i++) printf("%3d",node[node1].seq_num[i]);
		printf("\n");
		printf("Nab matrix for %d and %d\n",node1,node2);
		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) {/*free_imatrix(Nab,0,NP,0,NP);*/ 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)||(xLeft>=10.0)){
			printf("Warning:distance between %s and %s is inf, set distance 10.0\n",node[node1].nodeName,node[node2].nodeName);
//			free_imatrix(Nab,0,NP,0,NP);
			return 10.0;
		}
//		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) { /*free_imatrix(Nab,0,NP,0,NP);*/ 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);
//	free_imatrix(Nab,0,NP,0,NP);
	return 11.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);
		}
	}
*/
	for(i=1;i<=seqLength;i++){
		if(Nab[1][i]!=0 && Nab[2][i]!=0 && alfa[i]>1.0e-8)
			value+=dpij(Nab[1][i],Nab[2][i],d,alfa[i])/pij(Nab[1][i],Nab[2][i],d,alfa[i]);
	}
	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]*alfa;
	}
	return pij_t;
}

double dpij_alfa(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]*exp(d*alfa*wr[i])*wr[i]*d;
	}
	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 at least 3 non 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);
	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]) { 
				if(node[i].seq_num[j]<0 || node[i].seq_num[j]>20) printf("%d %d %d\n",i,j,node[i].seq_num[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]=PI[j];
		}
	}
	fprintf(resultFile,"PI vector for each site (*10^(-4))\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]*10000.0);
	fprintf(resultFile,"\n");
	for(i=1;i<=seqLength;i++){
		//fprintf(resultFile,"%5d ",i);
		//for(j=1;j<=NP;j++) fprintf(resultFile,"%6d ",aliMatrix[i][j]);
		//fprintf(resultFile,"\n");
		fprintf(resultFile,"%5d ",i);
		for(j=1;j<=NP;j++) fprintf(resultFile,"%6.0f ",sitePI[i][j]*10000.0);
		fprintf(resultFile,"\n");
	}

/*
	char *am="-WFYMLIVACGPTSNQDEHRKBZX*.wfymlivacgptsnqdehrkbzx";
	int **aliMatrix; double **aliPseudo;
	aliMatrix = imatrix(1,seqNum,1,seqLength);
	aliPseudo = dmatrix(1,seqLength,0,NP);
	for(i=1;i<=seqNum;i++){
		for(j=1;j<=seqLength;j++) {
			for(k=0;k<=48;k++) {
				if(node[i].seq_aa[j]==am[k]) { aliMatrix[i][j]=k; break; }
				if(k==48) printf("Something wrong in getting aliMatrix\n");
			}
		}
	}
	prfprfmatrix(aliMatrix,seqLength,seqNum,aliPseudo);
	for(i=1;i<=seqLength;i++){
		for(j=1;j<=NP;j++){
			if(j==1) sitePI[i][18]=aliPseudo[i][j];
			else if(j==2) sitePI[i][14]=aliPseudo[i][j];
			else if(j==3) sitePI[i][19]=aliPseudo[i][j];
			else if(j==4) sitePI[i][13]=aliPseudo[i][j];
			else if(j==5) sitePI[i][11]=aliPseudo[i][j];
			else if(j==6) sitePI[i][10]=aliPseudo[i][j];
			else if(j==7) sitePI[i][20]=aliPseudo[i][j];
			else if(j==8) sitePI[i][1]=aliPseudo[i][j];
			else if(j==9) sitePI[i][5]=aliPseudo[i][j];
			else if(j==10) sitePI[i][8]=aliPseudo[i][j];
			else if(j==11) sitePI[i][15]=aliPseudo[i][j];
			else if(j==12) sitePI[i][17]=aliPseudo[i][j];
			else if(j==13) sitePI[i][16]=aliPseudo[i][j];
			else if(j==14) sitePI[i][3]=aliPseudo[i][j];
			else if(j==15) sitePI[i][6]=aliPseudo[i][j];
			else if(j==16) sitePI[i][4]=aliPseudo[i][j];
			else if(j==17) sitePI[i][7]=aliPseudo[i][j];
			else if(j==18) sitePI[i][9]=aliPseudo[i][j];
			else if(j==19) sitePI[i][2]=aliPseudo[i][j];
			else if(j==20) sitePI[i][12]=aliPseudo[i][j];
		}
	}
	//normalize sitePI for each site
	for(i=1;i<=seqLength;i++){
		double normal=0.0;
		for(j=1;j<=NP;j++) normal+=sitePI[i][j];
		if(normal>0.0) for(j=1;j<=NP;j++) sitePI[i][j]/=normal;
	}

	fprintf(resultFile,"Psudocount 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,"%6d ",aliMatrix[i][j]);
		//fprintf(resultFile,"\n");
		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(double *pi)
{
	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++){
		Q_ij[i][i]=0.0;
		for(j=1;j<=NP;j++){
			if(j!=i) Q_ij[i][i]+=(-1.0)*Q_ij[i][j];
		}
	}
/*
	printf("Q_ij Matrix (before normalization) :\n");
	for(i=1;i<=NP;i++){
		for(j=1;j<=NP;j++) printf("%15.8f ",Q_ij[i][j]);
		printf("\n");
	}
*/
	//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("%15.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("%15.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;
	}
/*
	printf("Matrix of eignvector matrix\n");
	for(i=1;i<=NP;i++) {
		for(j=1;j<=NP;j++) printf("%15.6e ",v[i][j]);
		printf("\n\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("%15.6e ",inv[i][j]);
		printf("\n\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("%15.6e ",check[i][j]);
		printf("\n\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("%15.6e ",check[i][j]);
		printf("\n");
	}

    //check if the Q=U*lenta*invU
    double **check1; check1=dmatrix(1,NP,1,NP);
    for(i=1;i<=NP;i++) for(j=1;j<=NP;j++) { check[i][j]=0.0; check1[i][j]=0.0; }
    for(i=1;i<=NP;i++) for(j=1;j<=NP;j++) check1[i][j]=v[i][j]*wr[j];
    for(i=1; i<=NP; i++)
            for(j=1; j<=NP; j++)
                    for(k=1; k<=NP; k++){
                            check[i][j]+=check1[i][k]*inv[k][j];
                    }
    printf("the check matrix for Q\n");
    for(i=1;i<=NP;i++) {
            for(j=1;j<=NP;j++) printf("%15.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)) {printf("Warning: can not calculate pij\n node1: %d and node2: %d\n",s1,s2); exit(0);}//return 1.0;
	if(PI[s2]<1.0e-10 && s1!=s2) return 0.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]);
	}
	if((pij_t<0.0)&&(fabs(pij_t)<1.0e-12)) pij_t=0.0;
	return pij_t;
}

double func_alfa(int i, int root, double alfa)
{
	double func_value=-1.0*calTree(i,root,alfa);
	return func_value;
}

struct left_right{
	bool exist;
	double data;
};

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

void alpha()
{
	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;
	int gapForEst1=0;
	bool *estFlag;
	estFlag=new bool [seqLength+1];
//	double idAVG=0.0,idSite;

	correlation= new double [seqLength+1];
	corP=dmatrix(1,NP,1,NP);

	//calculate the average identity
//	for(i=1;i<seqNum;i++) for(j=i+1;j<=seqNum;j++) idAVG+=idMatrix[i][j];
//	idAVG=idAVG/(1.0*seqNum*(seqNum-1)/2.0);
//	fprintf(resultFile, "average identity: %4.0f\n", idAVG*100);

//	if(gapForEst==0) gapForEst=(int) (seqNum*0.5);
	if(gapForEst<1.000000001) { gapForEst1=(int) (gapForEst*seqNum); if(gapForEst1<3) gapForEst1=3;}
	else if(gapForEst>2.99999999) gapForEst1=(int) (gapForEst+0.000001); //if(gapForEst1>seqNum) gapForEst=seqNum;
//	printf("gapForEst1: %d\n",gapForEst1);
	for(i=1;i<=seqLength;i++) { if(seqNum-gapNumber[i]>=gapForEst1) estFlag[i]=true; else estFlag[i]=false; }
	for(i=1;i<=seqLength;i++) if(!estFlag[i]) gapNum++;
//	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]) {
			if(estFlag[i]) {
				pairNum=0; //idSite=0.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]]; /*idSite+=idMatrix[j][k];*/ pairNum++;}
					}
				}
				correlation[i]=pow(1.0*pairNum/correlation[i],1.3);
//				correlation[i]=correlation[i]/pairNum;
//				idSite=idSite/pairNum;
//				correlation[i]=idAVG/(idSite*correlation[i]);
//				correlation[i]=pow(correlation[i],1.3);
//				correlation[i]=pow(1.0*idSite/(correlation[i]*idAVG), 1.3);
				norCor+=correlation[i];
//				if(i==67) fprintf(resultFile,"site 67: idsite: %4.0f\n", idSite*100);
			}
			else correlation[i]=0.0;
		}
		//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;
//					printf("invariant site: %d\n",i);
				}
				invariantFlag=true;
			}
			else if(estFlag[i] && 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;
//					printf("invariant site: %d\n",i);
				}
				invariantFlag=true;
			}
			else if(!estFlag[i] && calFlag[i]) alfa[i]=1.0;
			else alfa[i]=0.0;

		}

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

double optimizeAlpha(int site)
{
	int i,j,k;
	int gapForEst1;

	if(gapForEst<1.000000001) { gapForEst1=(int) (gapForEst*seqNum); if(gapForEst1<3) gapForEst1=3;}
	else if(gapForEst>2.99999999) gapForEst1=(int) (gapForEst+0.000001);

	//check if gapNum > 0.5*seqNum
//	printf("site:%3d gapForEst1:%3d gapNumber:%3d seqNum:%3d\n",site,gapForEst1,gapNumber[site],seqNum);
	if((seqNum-gapNumber[site])<gapForEst1) return 1.0;
//	printf("This is site: %3d\n",site); fflush(stdout);
	//check if this site is invariant site
	char aminoacid; bool aaFlag=true;
	for(i=1;i<=seqNum;i++) if(node[i].seq_aa[site]!='-') { aminoacid=node[i].seq_aa[site]; break; }
	for(i=1;i<=seqNum;i++) if(node[i].seq_aa[site]!='-' && node[i].seq_aa[site]!=aminoacid) { aaFlag=false; break; }
	if(aaFlag) return 1.0e-8;
	//optimize alpha
		double maxL,ax,bx,cx,ay,by,cy,maxX,xmin;
		int dimension=10;
		double *xvalue,*yvalue;
		xvalue=new double [dimension];
		yvalue=new double [dimension];
//		maxL=dvector(1,seqLength);
//		maxX=dvector(1,seqLength);
		for(i=site;i<=site;i++) {
				xvalue[0]=1.0e-8; yvalue[0]=calTree(i,treeRoot[i],xvalue[0]);
/*
			for(j=0;j<dimension;j++) {
					xvalue[j]=exp((log(10)-log(1.0e-8))/(1.0*(dimension-1))*j+log(1.0e-8));
					yvalue[j]=calTree(i,treeRoot[i],xvalue[j]);
				}
*/
				xvalue[1]=1.0e-3; yvalue[1]=calTree(i,treeRoot[i],xvalue[1]);
	//			for(j=2;j<dimension;j++) {xvalue[j]=1.0*(j-1)*seqLength/(dimension-2); yvalue[j]=calTree(i,treeRoot[i],xvalue[j]); }
				for(j=2;j<=5;j++) {xvalue[j]=1.0+3.0*(j-2);  yvalue[j]=calTree(i,treeRoot[i],xvalue[j]); }
				for(j=6;j<=9;j++) { xvalue[j]=10.0+1.0*(j-5)*(seqLength-10)/4; yvalue[j]=calTree(i,treeRoot[i],xvalue[j]); }
				double maxY; int indexY;
				maxY=yvalue[0]; indexY=0;
				for(j=1;j<dimension;j++) if(yvalue[j]>maxY) { maxY=yvalue[j]; indexY=j; }
	//			for(j=0;j<=9;j++) printf("%10.5e  %10.5e\n",xvalue[j],yvalue[j]);
	//			printf("max xvalue: %10.5e  yvalue: %10.5e\n",xvalue[indexY],maxY);
				if(indexY==dimension-1) {//printf("site %d can not converge\n",i); maxL[i]=maxY; maxX[i]=xvalue[indexY]; continue; }
					//calculate the deritive of L
					int sign=1; double deriL=0;
					deriL=exp(calTree_alfa(i,treeRoot[i],xvalue[dimension-1],&sign)-calTree(i,treeRoot[i],xvalue[dimension-1]));
					deriL*=sign;
					if(deriL>=0.0) { maxX=xvalue[dimension-1]; maxL=calTree(i,treeRoot[i],maxX); }
					else { 
						printf("Warning: difficult to find the max likelihood at site %d\n",i);
						ax=xvalue[0];bx=xvalue[dimension/2];cx=xvalue[dimension-1];
						maxL=-1.0*dbrent1(ax,bx,cx,func_alfa,FTOL,&xmin,i,treeRoot[i]);
						maxX=xmin;
					}
				}
					
				else {
					ax=xvalue[indexY-1];bx=xvalue[indexY];cx=xvalue[indexY+1];
	//				printf("ax: %10.5e  bx: %10.5e  cx: %10.5e\n",ax,bx,cx);
					maxL=-1.0*dbrent1(ax,bx,cx,func_alfa,FTOL,&xmin,i,treeRoot[i]);
					maxX=xmin;
				}
		}

		return maxX;
}



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
	double CValue;//used for normalize optimal alpha

	correlation= 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");
	if((refineFlag==0)||(refineFlag==2)){
		//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
		if(D_A_P_Flag!=4){
			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;
				if(D_A_P_Flag==4){
					for(j=1;j<=NP;j++) PI[j]=sitePI[i][j];
					cal_Q(PI);
					cal_Eign();
					for(j=1;j<=NP;j++) for(k=1;k<=NP;k++) corP[j][k]=pij(j,k,1.0,1.0);
					for(j=1;j<=NP;j++) for(k=1;k<=j;k++){
						corP[j][k]=(corP[j][k]+corP[k][j])/2.0; corP[k][j]=corP[j][k];
					}
				}
				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]=pow(1.0*pairNum/correlation[i],1.3);
				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;
		}
*/
	}
	else if(refineFlag==1)	for(i=1;i<=seqLength;i++) alfa[i]=1.0;

//	for(i=1;i<=seqLength;i++) printf("%12.8f\n",alfa[i]);

	//test for shape of likelihood function
/*	double **logL; int dim=2001;
	int sign=1;
	logL=dmatrix(1,seqLength,1,dim);
	for(i=8;i<=8;i++){
		for(j=1;j<=dim;j++){
			logL[i][j]=exp(calTree_alfa(i,treeRoot[i],(0.000000001+0.1*(j-1)),&sign)-calTree(i,treeRoot[i],(0.000000001+0.1*(j-1))));
			logL[i][j]*=sign;
			printf("%4d %4d %12.5e %12.5e %12.5e\n",i,j,logL[i][j],calTree_alfa(i,treeRoot[i],(0.000000001+0.005*(j-1)),&sign),calTree(i,treeRoot[i],(0.000000001+0.005*(j-1))));

			logL[i][j]=calTree(i,treeRoot[i],(1.0e-8+0.1*(j-1)));
//			printf("%4d %4d %12.5e\n",i,j,logL[i][j]);
		}
	}
	getchar();
	fprintf(resultFile,"print out tree root for each site\n");
	for(i=1;i<=seqLength;i++) fprintf(resultFile,"%3d ",treeRoot[i]);
	fprintf(resultFile,"\n");
	fprintf(resultFile,"print the test result for shape of likelihood function\n");
	for(j=1;j<=dimension;j++){
		for(i=1;i<=seqLength;i++) fprintf(resultFile,"%10.5e ",logL[i][j]);
		fprintf(resultFile,"\n");
	}
*/
	//calculate the LOG likelyhood
	totalLikely=0.0;
	for(i=1;i<=seqLength;i++){
		result[i]=0.0;
		Anneal_site_func1=i;
		if(!dashFlag[i]||calFlag[i]) {
			if(D_A_P_Flag==4){
				for(j=1;j<=20;j++) PI[j]=sitePI[i][j];
				cal_Q(PI);
				cal_Eign();
			}
			result[i]=calTree(i,treeRoot[i],alfa[i]); totalLikely+=result[i];
//			fprintf(resultFile,"site:%4d likelihood: %10.5f  %10.5f\n",i,result[i],func1(PI));
		}
//		if(i==83) printf("siteNo: %d Likelihood: %12.5f\n", i, result[i]);
	}
	fprintf(resultFile,"\n\nThe original LOG Likelyhood: %21.15f\n",totalLikely);

	if(refineFlag==0){
		printf("optimizing alpha vector...\n");
		//search for maximum for each likelihood function
		double *maxL,ax,bx,cx,ay,by,cy,*maxX,xmin;
		bool *invantFlag;
		int dimension=10;
		double *xvalue,*yvalue;
		xvalue=new double [dimension];
		yvalue=new double [dimension];
		maxL=dvector(1,seqLength);
		maxX=dvector(1,seqLength);
		invantFlag=new bool [seqLength+1];
		for(i=1;i<=seqLength;i++){
			int invantBegin=-1;
			invantFlag[i]=true;
			for(j=1;j<=seqNum;j++) if(node[j].seq_aa[i]!='-') {invantBegin=j; break;}
			if(invantBegin==-1) invantFlag[i]=true;
			else if(invantBegin==seqLength) invantFlag[i]=false;
			else 
				for(j=invantBegin+1;j<=seqNum;j++) if(node[j].seq_aa[i]!=node[invantBegin].seq_aa[i] && node[j].seq_aa[i]!='-') { invantFlag[i]=false; break; }
		}
		for(i=1;i<=seqLength;i++) {
			if(!invantFlag[i] && (!dashFlag[i] || calFlag[i])){
				xvalue[0]=1.0e-8; yvalue[0]=calTree(i,treeRoot[i],xvalue[0]);
				xvalue[1]=1.0e-3; yvalue[1]=calTree(i,treeRoot[i],xvalue[1]);
	//			for(j=2;j<dimension;j++) {xvalue[j]=1.0*(j-1)*seqLength/(dimension-2); yvalue[j]=calTree(i,treeRoot[i],xvalue[j]); }
				for(j=2;j<=5;j++) {xvalue[j]=1.0+3.0*(j-2);  yvalue[j]=calTree(i,treeRoot[i],xvalue[j]); }
				for(j=6;j<=9;j++) { xvalue[j]=10.0+1.0*(j-5)*(seqLength-10)/4; yvalue[j]=calTree(i,treeRoot[i],xvalue[j]); }
				double maxY; int indexY;
				maxY=yvalue[0]; indexY=0;
				for(j=1;j<dimension;j++) if(yvalue[j]>maxY) { maxY=yvalue[j]; indexY=j; }
	//			for(j=0;j<=9;j++) printf("%10.5e  %10.5e\n",xvalue[j],yvalue[j]);
	//			printf("max xvalue: %10.5e  yvalue: %10.5e\n",xvalue[indexY],maxY);
				if(indexY==dimension-1) {//printf("site %d can not converge\n",i); maxL[i]=maxY; maxX[i]=xvalue[indexY]; continue; }
					//calculate the deritive of L
					int sign=1; double deriL=0;
					deriL=exp(calTree_alfa(i,treeRoot[i],xvalue[dimension-1],&sign)-calTree(i,treeRoot[i],xvalue[dimension-1]));
					deriL*=sign;
					if(deriL>=0.0) { maxX[i]=xvalue[dimension-1]; maxL[i]=calTree(i,treeRoot[i],maxX[i]); }
					else { 
						printf("Warning: difficult to find the max likelihood at site %d\n",i);
						ax=xvalue[0];bx=xvalue[dimension/2];cx=xvalue[dimension-1];
						maxL[i]=-1.0*dbrent1(ax,bx,cx,func_alfa,FTOL,&xmin,i,treeRoot[i]);
						maxX[i]=xmin;
					}
				}
					
				else {
					ax=xvalue[indexY-1];bx=xvalue[indexY];cx=xvalue[indexY+1];
	//				printf("ax: %10.5e  bx: %10.5e  cx: %10.5e\n",ax,bx,cx);
					maxL[i]=-1.0*dbrent1(ax,bx,cx,func_alfa,FTOL,&xmin,i,treeRoot[i]);
					maxX[i]=xmin;
				}
			}
			else if(invantFlag[i] && (!dashFlag[i] || calFlag[i]))
				{ maxX[i]=1.0e-8; maxL[i]=calTree(i,treeRoot[i],0.00000001); }
			else {maxX[i]=1.0e-8; maxL[i]=0.0;}
			printf("site: %4d  E alpha: %10.5f alfa value: %10.5f  L value: %10.5e\n",i,alfa[i],maxX[i],maxL[i]);
		}
//		getchar();
		printf("Finish finding each alpha value for max L\n");
		double total_alfa=0.0;
		for(i=1;i<=seqLength;i++) { total_alfa+=maxX[i]; totalLikely+=maxL[i]; }
//		printf("total_alfa: %15.5f    totalLikely: %12.5f\n",total_alfa, totalLikely);
		double *amplier, mean, var,*Zscore;
		amplier=dvector(1,seqLength); Zscore=dvector(1,seqLength);
		for(i=1;i<=seqLength;i++) amplier[i]=Zscore[i]=0.0;
		for(i=1;i<=seqLength;i++) if(!dashFlag[i] || calFlag[i]) {
			if(invantFlag[i]) amplier[i]=0.0;
			else amplier[i]=maxX[i]/alfa[i];
		}
		bool outlierFlag=true;
		int loop=0;
		while(outlierFlag){
			loop++;
	//		printf("loop: %d\n",loop++);
			outlierFlag=false;
			mean=0.0;
			int length=seqLength;
			for(i=1;i<=seqLength;i++){
				if(!invantFlag[i] && (!dashFlag[i] || calFlag[i])) mean+=amplier[i];
				else length--;
			}
	//		for(i=1;i<=seqLength;i++) printf("%12.5f  %12.5f  %12.5f\n",maxX[i],alfa[i],amplier[i]);
			mean/=length;
//			printf("mean: %12.5f\n",mean);
			var=0.0;
			for(i=1;i<=seqLength;i++) if(!invantFlag[i] && (!dashFlag[i] || calFlag[i])) var+=(amplier[i]-mean)*(amplier[i]-mean);
			var=sqrt(var/(length-1));
//			printf("variance: %12.5f\n",var);
			for(i=1;i<=seqLength;i++) if(!invantFlag[i] && (!dashFlag[i] || calFlag[i])) Zscore[i]=(amplier[i]-mean)/var;
			for(i=1;i<=seqLength;i++) if(!invantFlag[i] && (!dashFlag[i] || calFlag[i])) if(Zscore[i]>3.0) { amplier[i]=3.0*var+mean; outlierFlag=true;}
			if(loop>20) outlierFlag=false;
//			printf("Loop: %d\n",loop++);
//			getchar();
		}
		for(i=1;i<=seqLength;i++) if(!invantFlag[i] && (!dashFlag[i] || calFlag[i])) maxX[i]=amplier[i]*alfa[i];
//		for(i=1;i<=seqLength;i++) printf("%12.5f   %12.5f\n",alfa[i],maxX[i]);
		total_alfa=0.0; totalLikely=0.0;
		for(i=1;i<=seqLength;i++) { 
			total_alfa+=maxX[i]; 
			if(!dashFlag[i] || calFlag[i]) {
				result[i]=calTree(i,treeRoot[i],maxX[i]);
				totalLikely+=result[i];
			}
			else result[i]=0.0;
		}
//		printf("total_alfa: %15.5f    totalLikely: %12.5f\n",total_alfa, totalLikely);
		//assign new PI
		CValue=total_alfa/(seqLength-gapNum);
		printf("CValue: %12.8f\n",CValue);
		for(i=1;i<=seqLength;i++) fprintf(resultFile,"%12.6f   %12.6f  %12.6f\n",alfa[i],maxX[i]/CValue,maxX[i]);
		for(i=1;i<=seqLength;i++) alfa[i]=maxX[i];
		printf("Finish optimizing alpha vector...\n");
	}
//	exit(0);
//	for(i=1;i<=seqLength;i++) fprintf(resultFile,"site: %4d  alfa value: %10.5f  L value: %10.5e\n",i,maxX[i],maxL[i]);



	if((refineFlag==0)||(refineFlag==2)){
		//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]);
			if(refineFlag==2) fprintf(resultFile,"%8.5f   %12.8f\n",alfa[i],result[i]);
			else fprintf(resultFile,"%8.5f   %12.8f\n",(alfa[i]/CValue),result[i]);
		}
	}

	if(rootOnly) return;
/*
	//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 finding optimal solution for each internal node\n");
	for(i=seqNum+1;i<=(seqNum-1)*2;i++){
		node[i].seq_aa[0]=' ';
//		secondBest[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]='-';
//			if(node[i].seq_num[j]<0 || node[i].seq_num[j]>20) printf("%d %d %d\n",i,j,node[i].seq_num[j]);
//			if(secondBest[i].seq_num[j]<0 || secondBest[i].seq_num[j]>20) printf("%d %d %d\n",i,j,secondBest[i].seq_num[j]);
			node[i].seq_aa[j]=aa[node[i].seq_num[j]];
//			secondBest[i].seq_aa[j]=aa[secondBest[i].seq_num[j]];
		}
	}
	printf("finish assigning letter for each internal node\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 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_alfa(int siteNo, int root, double alfa, int *sign)
{
	int i,j,k,l;
	int maxLevel=0;
	double *siteLikely,*amplify;
	siteLikely=dvector(1,seqNum*2-2);
	amplify=dvector(1,seqNum*2-2);
	for(i=1;i<=seqNum*2-2;i++) siteLikely[i]=amplify[i]=0.0;
	for(j=1;j<=seqNum*2-2;j++) if(j!=root && gapTree[siteNo][j].dist>1.0e-10) {
//		printf("derive at node: %d\n",j);
		calRL_alfa(siteNo,gapTree[siteNo][root].childL,alfa,j);
		calRL_alfa(siteNo,gapTree[siteNo][root].childR,alfa,j);
		calRL_alfa(siteNo,gapTree[siteNo][root].childE,alfa,j);
		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[j]+=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;
		amplify[j]=gapTree[siteNo][root].amp;
	}
	//find the largest amplify
	double maxAmp=0;
	for(i=1;i<=seqNum*2-2;i++) if(amplify[i]>maxAmp) maxAmp=amplify[i];
	for(i=1;i<=seqNum*2-2;i++) siteLikely[i]=siteLikely[i]*pow(10.0,maxAmp-amplify[i]);
	double total=0;
	for(i=1;i<=seqNum*2-2;i++) total+=siteLikely[i];
	if(total<0.0) {*sign=-1; total=-1.0*total; }
	else *sign=1;
	total=log(total)-maxAmp*log(10.0);

	return total;
}

void calRL_alfa(int siteNo,int nodeNo,double alfa,int deri_nodeNo)
{
/*	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_alfa(siteNo,gapTree[siteNo][nodeNo].childL,alfa,deri_nodeNo);
		calRL_alfa(siteNo,gapTree[siteNo][nodeNo].childR,alfa,deri_nodeNo);
		for(int k=1;k<=NP;k++){
			gapTree[siteNo][nodeNo].L[k]=0.0;
			for(int l=1;l<=NP;l++){
				if(deri_nodeNo==nodeNo)
					gapTree[siteNo][nodeNo].L[k]+=dpij_alfa(k,l,gapTree[siteNo][nodeNo].dist,alfa)*gapTree[siteNo][gapTree[siteNo][nodeNo].childL].L[l]*gapTree[siteNo][gapTree[siteNo][nodeNo].childR].L[l];
				else 
					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++){
			if(deri_nodeNo==nodeNo){ //printf("nodeNo: %d  site amino acid: %d ",nodeNo,node[nodeNo].seq_num[siteNo]);
				gapTree[siteNo][nodeNo].L[k]=dpij_alfa(k,node[nodeNo].seq_num[siteNo],gapTree[siteNo][nodeNo].dist,alfa);
//				printf("%10.5e\n",gapTree[siteNo][nodeNo].L[k]);
			}

			else
				gapTree[siteNo][nodeNo].L[k]=pij(k,node[nodeNo].seq_num[siteNo],gapTree[siteNo][nodeNo].dist,alfa);
//			if((gapTree[siteNo][nodeNo].L[k]<0.0)&&(fabs(gapTree[siteNo][nodeNo].L[k])<1.0e-10))
//				gapTree[siteNo][nodeNo].L[k]=0.0;
			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);
	}
}

/*
void findLeafProb(int siteNo, int nodeNo, double alfa)
{
	if(gapTree[siteNo][nodeNo].rootFlag){
		for(int i=1;i<=NP;i++) {
			gapTree[siteNo][gapTree[siteNo][nodeNo].childL].L[i]=0.0; gapTree[siteNo][gapTree[siteNo][nodeNo].childL].amp=0.0;
			gapTree[siteNo][gapTree[siteNo][nodeNo].childR].L[i]=0.0; gapTree[siteNo][gapTree[siteNo][nodeNo].childR].amp=0.0;
			gapTree[siteNo][gapTree[siteNo][nodeNo].childE].L[i]=0.0; gapTree[siteNo][gapTree[siteNo][nodeNo].childE].amp=0.0;
			for(int j=1;j<=NP;j++){
				gapTree[siteNo][gapTree[siteNo][nodeNo].childL].L[i]+=PI[j]*pij(j,i,gapTree[siteNo][gapTree[siteNo][nodeNo].childL].dist,alfa);
				gapTree[siteNo][gapTree[siteNo][nodeNo].childR].L[i]+=PI[j]*pij(j,i,gapTree[siteNo][gapTree[siteNo][nodeNo].childR].dist,alfa);
				gapTree[siteNo][gapTree[siteNo][nodeNo].childE].L[i]+=PI[j]*pij(j,i,gapTree[siteNo][gapTree[siteNo][nodeNo].childE].dist,alfa);
			}
		}
		findLeafProb(siteNo,gapTree[siteNo][nodeNo].childL,alfa);
		findLeafProb(siteNo,gapTree[siteNo][nodeNo].childR,alfa);
		findLeafProb(siteNo,gapTree[siteNo][nodeNo].childE,alfa);
	}
	else if(nodeNo>seqNum){
*/

//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);
	if(gapTree[siteNo][root].gapFlag) return 0.0;
	calRL(siteNo,gapTree[siteNo][root].childL,alfa);
	calRL(siteNo,gapTree[siteNo][root].childR,alfa);
	calRL(siteNo,gapTree[siteNo][root].childE,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);
/*
	printf("Root:%3d ",root);
	for(int k=1;k<=NP;k++) printf("%12.6e ",gapTree[siteNo][root].L[k]);
	printf("AMP:%3.0f\n",gapTree[siteNo][root].amp);
*/
	return siteLikely;
}


void calRL(int siteNo,int nodeNo,double alfa)
{
	
	if((gapTree[siteNo][nodeNo].childL!=0)&&(gapTree[siteNo][nodeNo].childR!=0)){
//		printf("nodeNo:%3d, childL:%3d childR:%3d\n",nodeNo,gapTree[siteNo][nodeNo].childL,gapTree[siteNo][nodeNo].childR);
		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
		if(big>1.0e-200){
			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 {
//		printf("leaf node:%3d, parent:%3d\n",nodeNo,gapTree[siteNo][nodeNo].parent);
		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((gapTree[siteNo][nodeNo].L[k]<0.0)&&(fabs(gapTree[siteNo][nodeNo].L[k])<1.0e-10))
//				gapTree[siteNo][nodeNo].L[k]=0.0;
			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
		if(big>1.0e-200){
			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);
		
	}
/*
	for(int k=1;k<=NP;k++) {
		for(int m=1;m<=NP;m++) printf("%12.6e ",pij(k,m,gapTree[siteNo][nodeNo].dist,1.0));
		printf("\n");
	}
	printf("Node:%3d ",nodeNo);
	for(int k=1;k<=NP;k++) printf("%12.6e ",gapTree[siteNo][nodeNo].L[k]);
	printf("AMP:%3.0f\n\n",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)
{
	double ran1(long *idum);
	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);
	} */
/*
	if(D_A_P_Flag==4){
		for(j=1;j<=20;j++) PI[j]=sitePI[siteNo][j];
		cal_Q(PI);
		cal_Eign();
	}
*/
	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;
//			tmpMax=-1.0; indexMax=0;
			for(j=2;j<=NP;j++){
				if(node[i].prob[siteNo][j]>tmpMax) { tmpMax=node[i].prob[siteNo][j]; indexMax=j;}
			}
//			if(indexMax==0) {printf("Warning: something wrong in optTree_Goldstein method\n"); }
			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;
//			node[i].seq_aa[siteNo]=aa[indexMax];

		}
//		else node[i].seq_aa[siteNo]=aa[0];
	}
	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;
//	printf("siteNo:%3d nodeNo:%3d root:%3d child:%3d\n",siteNo,nodeNo,root,child);
//	getchar();
	if((nodeNo==root)&&(gapTree[siteNo][nodeNo].rootFlag)){}
	else{
	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==26){
	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);
*/
	}
}

char *produceTree(int nodeNo)
{
	char treeLine[10000];
	if(gapTree[0][nodeNo].rootFlag) {
		strcpy(treeLine,"(");
		strcat(treeLine,produceTree(gapTree[0][nodeNo].childL));
		strcat(treeLine,",");
//		strcat(treeLine,produceTree(node[nodeNo].childR));
//		strcat(treeLine,",");
		strcat(treeLine,produceTree(gapTree[0][nodeNo].childR));
		strcat(treeLine,");");
		return treeLine;
	}
	else {
		if(nodeNo>seqNum) {
			strcpy(treeLine,"(");
			strcat(treeLine,produceTree(gapTree[0][nodeNo].childL));
			strcat(treeLine,",");
			strcat(treeLine,produceTree(gapTree[0][nodeNo].childR));
			strcat(treeLine,"):");
			sprintf(treeLine+strlen(treeLine),"%f",gapTree[0][nodeNo].dist);
			return treeLine;
		}
		else {
			sprintf(treeLine,"%s:%f",node[nodeNo].nodeName,gapTree[0][nodeNo].dist);
			return treeLine;
		}
	}
}

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=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);
	}
	branch=new BRANCH [DIM-1];
//	gapTree=topmatrix(0,seqLength+1,0,(seqNum*2-2));//seqLength);
	gapTree=new topology * [seqLength+2];
	for(i=0;i<seqLength+2;i++) gapTree[i]=new topology [seqNum*2-1];
	copyTree=new topology * [2];
	for(i=0;i<2;i++) copyTree[i]=new topology [seqNum*2-1];

	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=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");
	fprintf(resultFile,"%27s"," ");
	for(i=0;i<=seqLength/10;i++) fprintf(resultFile,"%-10d",i*10);
	fprintf(resultFile,"\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);
	//
	for(i=1;i<=seqNum;i++) for(j=1;j<=seqLength;j++) node[i].seq_aa[j]=toupper(node[i].seq_aa[j]);

	//assign seq_num
	for(i=1;i<=seqNum;i++){
		for(j=1;j<=seqLength;j++){
			bool unknownChar=true;
			for(k=0;k<=NP;k++){
				if(node[i].seq_aa[j]==aa[k]) { unknownChar=false; node[i].seq_num[j]=k; break; }
			}
			if(unknownChar) {node[i].seq_num[j]=0; node[i].seq_aa[j]='-';}
		}
	}
	//assign all ancestral sequence name
	strcpy(node[0].nodeName,"anc0");
	for(i=seqNum+1;i<=seqNum*2-2;i++) sprintf(node[i].nodeName,"anc%d",i);
/*	for(i=1;i<=seqNum;i++){
		for(j=1;j<=seqLength;j++) printf("%2d ",node[i].seq_num[j]);
		printf("\n");
	}
*/
//	for(i=1;i<=seqNum;i++) fprintf(resultFile,"%-20s %s\n",node[i].nodeName, node[i].seq_aa);
//        printf("sequence number before deleting duplication: %d\n",seqNum);
//        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);
//        printf("sequence number after deleting duplication: %d\n",seqNum);
	
	alfa= dvector(0,seqLength);
	//assign alfa vector
	for(i=1;i<=seqLength;i++) alfa[i]=1.0;

	//assign gapNumber vector
	gapNumber=ivector(1,seqLength);
	for(i=1;i<=seqLength;i++) {
		gapNumber[i]=0;
//		if(dashFlag[i]){
			for(j=1;j<=seqNum;j++) if(node[j].seq_aa[i]=='-') gapNumber[i]++;
//		}
//		fprintf(resultFile,"%4d\n",gapNumber[i]);
	}

	//calculate idMatrix
/*	idMatrix = dmatrix(1,seqNum,1,seqNum);
	for(i=1;i<seqNum;i++){
		for(j=i+1;j<=seqNum;j++){
			int count=0; int length=0;
			for(k=1;k<=seqLength;k++){
				if(node[i].seq_aa[k]!='-' && node[j].seq_aa[k]!='-'){
					length++;
					if(node[i].seq_aa[k]==node[j].seq_aa[k]) count++;
				}
			}
//			if(length>0) idMatrix[i][j]=1.0*count/length;
//			else idMatrix[i][j]=0.0;
			if(length>0) {
				idMatrix[i][j]=1.0*count/length; 
				if(idMatrix[i][j]>0.05) { idMatrix[i][j]=0.95/(idMatrix[i][j]-0.05)-1; if(idMatrix[i][j]>10.0) idMatrix[i][j]=10.0; }
				else idMatrix[i][j]=10.0;
			}
			else idMatrix[i][j]=10.0;
		}
	}
*/
//	for(i=2;i<=seqNum;i++) for(j=i-1;j<seqNum;j++) idMatrix[i][j]=idMatrix[j][i];
//	for(i=1;i<=seqNum;i++) idMatrix[i][i]=0.0;
//	for(i=1;i<seqNum;i++) for(j=i+1;j<=seqNum;j++) fprintf(resultFile,"row: %4d   col: %4d   id: %4.0f\n",i,j,idMatrix[i][j]*100);
//	printf("Finish reading alignment\n");
}


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

