// 08/19/02 add -G parameter to correct estimated alpha calcualation

#include "ancestral.h"
//#include "alpha.c"


void assignBranch(int nodeNo);
void findMidpoint();
int reconstructTree(int siteNo, int root, bool rootedFlag);
void reconstruct(int siteNo, int nodeNo);
void readSpecies(char *predFile);
char *produceTree(int nodeNo);

//Global Variables
NODE *node;
//double S_ij[NP+1][NP+1];
//double Q_ij[NP+1][NP+1];
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);

void usage()
{
	printf("Usage: likely2_produceTree -i sequencefile -t treefile [-o outputfile] -T species file\n");
	exit(1);
}
	

int main(int argc, char **argv)
{
	bool optPIFlag=false;  //used for optimize PI
	bool rootOnly=false;
	int i;
	int refineFlag=1, conditionalProb=0;
	bool flag=true,treefileFlag=false,bioRoot=false,ZFlag=false;
	int refineSeqItem=1,cutDistItem=1;
	FILE *input; int number;
	char seqFile[50], treeFile[50],predFile[50];
	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 'i':
				++i;
				optidx=0;
				if(i<argc) {
					strcpy(seqFile,argv[i]);
//					printf("%s\n",seqFile);
					input=fopen(seqFile,"r");
					if(!input) {
						printf("can not open file:%s\n",seqFile);
						exit(1);
					}
					fclose(input);
				}
				else usage();
				break;
			case 't':
				++i;
				optidx=0;
				if(i<argc) {
					treefileFlag=true;
					strcpy(treeFile,argv[i]);
				}
				else usage();
				break;
			case 'o':
				++i;
				optidx=0;
				if(i<argc) {
					resultFile=fopen(argv[i],"w");
					if(!resultFile) {
						printf("can not open file:%s\n",seqFile);
						exit(1);
					}
				}
				else usage();
				break;
			case 'T':
				++i;
				optidx=0;
				if(i<argc) {
					strcpy(predFile,argv[i]);
//					printf("%s\n",seqFile);
					input=fopen(predFile,"r");
					if(!input) {
						printf("can not open file:%s\n",predFile);
						exit(1);
					}
					fclose(input);
				}
				else usage();
				break;
			default:
				usage();
		}
	}

	//read S_ij matrix
	calTime1=time(NULL);
	init_S();
	readali(seqFile);
	if(!treefileFlag) {
		usage();
	}
	else readTree(treeFile);
    assignTree();
	readSpecies(predFile);

}

void readSpecies(char *predFile)
{
	FILE *input;
	int i,j;
	char **estSites;
	char line[2000];
	int siteIndex=0;
	input=fopen(predFile,"r");
	if(!input) { printf("can not open file: %s\n",predFile); exit(1); }
        while(fgets(line,2000,input)!=NULL){
                if(strlen(line)>1) {
                        ++siteIndex;
                }
        }
        rewind(input);
	//define data array
	estSites = new char* [siteIndex+1];
	for(i=0;i<=siteIndex;i++) estSites[i]= new char [50];

    siteIndex=0;
    while(fgets(line,2000,input)!=NULL){
            if(strlen(line)>1) {
                    strcpy(estSites[++siteIndex], line);
            }
    }
    fclose(input);

	for(i=1;i<=siteIndex;i++) {
		estSites[i][strlen(estSites[i])-1]='\0';
//		printf("%s\n",estSites[i]);
	}

	//change the gapTree[0] 
	for(i=1;i<=seqNum;i++) {
		bool estSiteFlag = false;
		for(j=1;j<=siteIndex;j++) {
			if(memcmp(node[i].nodeName, estSites[j],strlen(estSites[j]))==0) {
				gapTree[0][i].gapFlag=false;
				estSiteFlag=true;
				break;
			}
		}
		if(!estSiteFlag) gapTree[0][i].gapFlag=true;
	}
/*
	for(i=1;i<=seqNum*2-2;i++) {
		printf("seq: %4d name: %-12s parent: %4d childL: %4d childR: %4d childE:%4d ",i,node[i].nodeName,gapTree[0][i].parent,gapTree[0][i].childL,gapTree[0][i].childR,gapTree[0][i].childE);
		if(gapTree[0][i].gapFlag) printf("gap: Yes\n");
		else printf("gap: No\n");
	}
	printf("\n\n");
*/
	int root0=rebuildTree(0);
	//test
/*	for(i=1;i<=seqNum*2-2;i++) {
		printf("seq: %4d name: %-12s parent: %4d childL: %4d childR: %4d childE:%4d ",i,node[i].nodeName,gapTree[0][i].parent,gapTree[0][i].childL,gapTree[0][i].childR,gapTree[0][i].childE);
		if(gapTree[0][i].gapFlag) printf("gap: Yes\n");
		else printf("gap: No\n");
	}
*/
	printf("%s\n",produceTree(root0));


}

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(gapTree[0][nodeNo].childR));
		strcat(treeLine,",");
		strcat(treeLine,produceTree(gapTree[0][nodeNo].childE));
		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;
		}
	}
}


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 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=0;i<=seqLength;i++){
		if(!dashFlag[i]) treeRoot[i]=seqNum+1;
		else if(calFlag[i]) treeRoot[i]=rebuildTree(i);
		else {
			treeRoot[i]=0;
			for(j=1;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++;}  
						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;
//								printf("node: %d dist: %f\n",j,node[j].dist);
//								printf("in tree Left: i=%d, j=%d, nodeNum=%d, dist=%f\n",i,j,node[j].nodeNum,node[j].dist);
								break;							}
//							else printf("No\n");
						}
					}
					else {//3
						if(treeR[i]==0) {
							tmpIndex=0;
							while((c=fgetc(input))!=':') { tmpName[tmpIndex]=c; tmpIndex++;}
//							printf("i=%d\n",i);
							for(j=1;j<=index;j++) {
								if(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;
//									printf("in tree Right: i=%d, j=%d, nodeNum=%d, dist=%f\n",i,j,node[j].nodeNum,node[j].dist);
									break;
								}
							}
						}
						else {//4
							tmpIndex=0;
							while((c=fgetc(input))!=':') { tmpName[tmpIndex]=c; tmpIndex++;}
							for(j=1;j<=index;j++) {
								if(tmpIndex==strlen(node[j].nodeName) && memcmp(node[j].nodeName,tmpName,tmpIndex)==0) {
									treeE=node[j].nodeNum;
									fscanf(input,"%f",&tmpValue); node[j].dist=tmpValue;
//									printf("in tree E: i=%d, j=%d, nodeNum=%d, dist=%f\n",i,j,node[j].nodeNum,node[j].dist);
									break;
								}
							}
						}//end else 4
					}//end else 3
				}//end else 2
			}//end else 1
		}//end if
	}//end while
	fclose(input);
//	printf("Finish reading tree file\n");
/*	for(i=index;i<=(index-1)*2-2;i++){
		printf("i: %d  treeLeft: %d  treeRight: %d\n",i,treeL[i],treeR[i]);
	}
	printf("treeE: %d\n",treeE);
	for(i=1;i<=(index-1)*2-2;i++){
		printf("node: %d, nodeName: %10s distance: %f\n", node[i].nodeNum, node[i].nodeName, node[i].dist);
	}
*/
	//input child, parent, level information into node[]
	//root node
	node[index].rootFlag=true;
	node[index].childL=treeL[index];node[index].childR=treeR[index];node[index].childE=treeE;
	node[index].parent=node[index].nodeNum;
	node[treeL[index]].parent=index;node[treeR[index]].parent=index;node[treeE].parent=index;
	//other nodes
	for(i=index+1;i<=(index-1)*2-2;i++){
		node[i].childL=treeL[i]; node[i].childR=treeR[i];
		node[treeL[i]].parent=i; node[treeR[i]].parent=i;
	}
//	printf("Finish assigning node, before assign level\n");
	setLevel(index);
/*
	for(i=1;i<=(index-1)*2-2;i++){
		levelIndex=0;
		levelP=i;
		while(node[levelP].parent!=index){
			levelIndex++;
			levelP=node[levelP].parent;
		}
		if(i!=index) levelIndex++;
		node[i].nodeLevel=levelIndex;
	}

	for(i=1;i<=(index-1)*2-2;i++){
		printf("nodeNum:%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;
}		


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

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);
*/
	}
}


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);
	}
//	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("%-20s %d\n",node[i].nodeName,strlen(node[i].nodeName));
	}
//	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]='-';}
		}
	}
/*	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);
/*	secondBest=new NODE [DIM];
	for(i=0;i<DIM;i++) { 
		secondBest[i].seq_aa=new char [seqLength+2];
		secondBest[i].seq_num=new int [seqLength+1];
		secondBest[i].prob= dmatrix(1,seqLength+1,1,NP+1);
	}
	for(i=0;i<DIM;i++) { secondBest[i].seq_aa[0]=' ';secondBest[i].seq_aa[seqLength+1]='\0';}
	for(i=0;i<DIM;i++) for(j=0;j<=seqLength;j++) secondBest[i].seq_num[j]=0;
*/	//compare sequence name to find duplicates
//	printf("try to find the duplicate name of different sequences\n");
	bool dupFlag=false;
	for(i=1;i<seqNum;i++){
		for(j=i+1;j<=seqNum;j++){
			if(strcmp(node[i].nodeName, node[j].nodeName)==0){ 
				dupFlag=true;
				printf("duplication: seqNo:%d seqName:%s  seqNo:%d seqName:%s\n",i,node[i].nodeName,j,node[j].nodeName);
			}
		}
	}
	if(dupFlag) exit(1);
	
	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;
}

