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

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

void cal_dist_matrix();
void assignTree(char *treeFile1,char *treeFile2);
void clusterTree(char *treeFile1, char *treeFile2, char *branchRank);

struct BRANCH{
//	int No;
	double length;
	int LPoint;
	int RPoint;
	int LPointNumber;
	int RPointNumber;
};

//Global Variables
NODE *node;
// NODE *secondBest;
int seqLength,seqNum,DIM;
char aa[NP+1]={'-','A','R','N','D','C','Q','E','G','H','I','L','K','M','F','P','S','T','W','Y','V'};
//for output file
FILE *resultFile=stdout;
topology **gapTree;
double **distMatrix1, **distMatrix2;

void usage()
{
	printf("Usage: compareTree -i sequencefile -t1 treefile1 -t2 treefile2 [-o outputfile]\n");
	printf("       -i   name: read sequences from file\n");
	printf("       -t1  name: original tree file\n");
	printf("       -t2  name: new tree file but with same tree structure but different branch length as original tree\n");
	printf("       -t3  name: rank file which ranks branches based on its statistical significance, used with -C option\n");
	printf("       -o   name: write the result into file or standard output device\n");
	printf("       -C   name: cluster tree into subgroups base on optimized tree and the branch rank list\n");
	exit(1);
}
	

int main(int argc, char **argv)
{
	int i;
	FILE *input; int number;
	char seqFile[50], treeFile1[50],treeFile2[50],rankFile[50];
	strcpy(seqFile,"empty"); strcpy(treeFile1,"empty"); strcpy(treeFile2,"empty"); strcpy(rankFile,"empty");
	bool clusterFlag=false;

	if(argc<=1) usage();
	for(i=1;i<argc;++i){
		char option,optidx;
		if(argv[i][0]=='-') option=argv[i][1];
		else usage();
		switch(option){
			case 'C':
				clusterFlag=true;
				break;
			case 'i':
				++i;
				optidx=0;
				if(i<argc) {
					strcpy(seqFile,argv[i]);
//					printf("%s\n",seqFile);
					input=fopen(seqFile,"r");
					if(!input) {
						printf("can not open file:%s\n",seqFile);
						exit(1);
					}
					fclose(input);
				}
				else usage();
				break;
			case 't':
				if(argv[i][2]=='1') {
					++i;
					optidx=0;
					if(i<argc) {
						strcpy(treeFile1,argv[i]);
					}
					else usage();
				}
				else if(argv[i][2]=='2') {
					++i;
					optidx=0;
					if(i<argc) {
						strcpy(treeFile2,argv[i]);
					}
					else usage();
				}
				else if(argv[i][2]=='3') {
					++i;
					optidx=0;
					if(i<argc) {
						strcpy(rankFile,argv[i]);
					}
					else usage();
				}

				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();
		}
	}
	if(strcmp(seqFile,"empty")==0 || strcmp(treeFile1,"empty")==0 || strcmp(treeFile2,"empty")==0 ) usage();
	if(clusterFlag) { if(strcmp(rankFile,"empty")==0) usage(); }
	readali(seqFile);
	if(clusterFlag) clusterTree(treeFile1, treeFile2, rankFile);
    	else assignTree(treeFile1,treeFile2);//assign tree structure to gapTree
}

void cal_dist_matrix()
{
	int i,j,k,m;
	distMatrix1=dmatrix(1,seqNum,1,seqNum);
	distMatrix2=dmatrix(1,seqNum,1,seqNum);
	int **parents;
	parents=imatrix(1,seqNum,1,seqNum-2);
	double **parentsDist1,**parentsDist2;
	parentsDist1=dmatrix(1,seqNum,1,seqNum-2);
	parentsDist2=dmatrix(1,seqNum,1,seqNum-2);
	int *parentsIndex;
	parentsIndex=ivector(1,seqNum);
	//initiate matrix and vector
	for(i=1;i<=seqNum;i++) {
		parentsIndex[i]=0;
		for(j=1;j<=seqNum-2;j++) {
			parentsDist1[i][j]=0.0;
			parentsDist2[i][j]=0.0;
		}
	}
	//find all the parent nodes and dist from leaf to root
	for(i=1;i<=seqNum;i++) {
		int parent=gapTree[1][i].parent;
                while(!gapTree[1][parent].rootFlag) {
                        parents[i][++parentsIndex[i]]=parent;
			if(parentsIndex[i]==1) {
				parentsDist1[i][parentsIndex[i]]=gapTree[1][i].dist;
				parentsDist2[i][parentsIndex[i]]=gapTree[2][i].dist;
			}
			else {
				parentsDist1[i][parentsIndex[i]]=parentsDist1[i][parentsIndex[i]-1]+gapTree[1][parents[i][parentsIndex[i]-1]].dist;
				parentsDist2[i][parentsIndex[i]]=parentsDist2[i][parentsIndex[i]-1]+gapTree[2][parents[i][parentsIndex[i]-1]].dist;
			}
                        parent=gapTree[1][parent].parent;
                }
		parents[i][++parentsIndex[i]]=parent;
		if(parentsIndex[i]==1) {
			parentsDist1[i][parentsIndex[i]]=gapTree[1][i].dist;
			parentsDist2[i][parentsIndex[i]]=gapTree[2][i].dist;
		}
		else {
			parentsDist1[i][parentsIndex[i]]=parentsDist1[i][parentsIndex[i]-1]+gapTree[1][parents[i][parentsIndex[i]-1]].dist;
			parentsDist2[i][parentsIndex[i]]=parentsDist2[i][parentsIndex[i]-1]+gapTree[2][parents[i][parentsIndex[i]-1]].dist;
		}
	}
	for(i=1;i<seqNum;i++) {
		for(j=i+1;j<=seqNum;j++) {
			for(k=1;k<=parentsIndex[i];k++) {
				for(m=1;m<=parentsIndex[j];m++) {
					if(parents[i][k]==parents[j][m]) {
						distMatrix1[i][j]=parentsDist1[i][k]+parentsDist1[j][m];
						distMatrix2[i][j]=parentsDist2[i][k]+parentsDist2[j][m];
						k=parentsIndex[i]+1; m=parentsIndex[j]+1;
						//break;
						//break;
					}
				}
			}
		}
	}
	//output test result
/*	for(i=1;i<=seqNum;i++) {
		fprintf(resultFile,"node:%4d\n",i);
		for(j=1;j<=parentsIndex[i];j++) 
			fprintf(resultFile,"tree1: parent: %4d dist: %8.3f, tree2: dist: %8.3f\n",parents[i][j], parentsDist1[i][j],parentsDist2[i][j]);
	}
	for(i=1;i<seqNum;i++) {
		for(j=i+1;j<=seqNum;j++) fprintf(resultFile,"%8.3f %8.3f\n",distMatrix1[i][j],distMatrix2[i][j]);
	}
*/
}

void clusterTree(char *treeFile1, char *treeFile2, char *branchRank)
{
	int i,j,k,m;
	//read tree information
	readTree(treeFile1);
	//store tree information into gapTree[1];
	for(i=1;i<=(seqNum*2-2);i++) {
		gapTree[1][i].childL=node[i].childL;
		gapTree[1][i].childR=node[i].childR;
		gapTree[1][i].childE=node[i].childE;
		gapTree[1][i].parent=node[i].parent;
		gapTree[1][i].rootFlag=node[i].rootFlag;
		gapTree[1][i].dist=node[i].dist;
	}
	//read new tree structure
	readTree(treeFile2);
	//store new tree information into gapTree[2]
	for(i=1;i<=(seqNum*2-2);i++) {
		gapTree[2][i].childL=node[i].childL;
		gapTree[2][i].childR=node[i].childR;
		gapTree[2][i].childE=node[i].childE;
		gapTree[2][i].parent=node[i].parent;
		gapTree[2][i].rootFlag=node[i].rootFlag;
		gapTree[2][i].dist=node[i].dist;
	}

	//calculate the distance matrix for both trees
	cal_dist_matrix();

	//read branch rank file
	FILE *rankFile;
	rankFile=fopen(branchRank,"r");
	if(!rankFile) {printf("Can not open file: %s\n",branchRank); exit(1); }
	char line[1000];
	int *rank;
	rank=ivector(1,seqNum*2-2);
	int rankIndex=0;
        while(fgets(line,1000,rankFile)!=NULL){
                bool flag=false;
                for(i=0;i<strlen(line);i++) if(isalnum(line[i])) { flag=true; break; }
                if(strlen(line)>1 && flag) {
			++rankIndex;
			if(rankIndex>seqNum*2-2) { printf("Too many branches! Check your rank list!\n"); exit(1); }
                        rank[rankIndex]=atoi(line);
                }
        }
        fclose(rankFile);
	if(rankIndex<seqNum*2-2) printf("Warning: your total branch number is less than expected!\n");

	//find left/right children for each node
        int **children, *childIndex, **children1;
        children=imatrix(1,seqNum*2-2,1,seqNum*2-2);
        childIndex=ivector(1,seqNum*2-2);
        children1=imatrix(1,seqNum*2-2,1,seqNum*2-2);
	//store all children including internal nodes
	int **children2, **children3,*childIndex1;
	childIndex1=ivector(1,seqNum*2-2);
	children2=imatrix(1,seqNum*2-2,1,seqNum*2-2);
	children3=imatrix(1,seqNum*2-2,1,seqNum*2-2);

        //search children for each node
        for(i=1;i<=seqNum*2-2;i++) {
//		if(i<=seqNum) { childIndex[i]=1; children[i][1]=i; childIndex1[i]=1; children2[i][1]=i; }
//		else { childIndex[i]=0; childIndex1[i]=0; }
		childIndex[i]=childIndex1[i]=0;
	}
        for(i=1;i<=seqNum;i++) {
            if(i!=seqNum+1) {
                int parent=gapTree[1][i].parent;
                while(!gapTree[1][parent].rootFlag) {
                        children[parent][++childIndex[parent]]=i;
                        parent=gapTree[1][parent].parent;
                }
                children[seqNum+1][++childIndex[seqNum+1]]=i;
            }
        }
        for(i=1;i<=seqNum*2-2;i++) {
            if(i!=seqNum+1) {
                int parent=gapTree[1][i].parent;
                while(!gapTree[1][parent].rootFlag) {
                        children2[parent][++childIndex1[parent]]=i;
                        parent=gapTree[1][parent].parent;
                }
                children2[seqNum+1][++childIndex1[seqNum+1]]=i;
            }
        }

	//add leaf node to children array
	for(i=1;i<=seqNum;i++) {
		childIndex[i]++;
		children[i][childIndex[i]]=i;
		childIndex1[i]++;
		children2[i][childIndex1[i]]=i;
	}

	//assign children1 with child not in children
	for(i=1;i<=seqNum*2-2;i++) {
		int tmp_childIndex=0;
		for(j=1;j<=seqNum;j++) {
			if(j!=i) {
				bool childFlag=false;
				for(k=1;k<=childIndex[i];k++) {
					if(children[i][k]==j) { childFlag=true; break; }
				}
				if(!childFlag) {
					children1[i][++tmp_childIndex]=j;
				}
			}
		}
		tmp_childIndex=0;
		for(j=1;j<=seqNum*2-2;j++) {
			if(j!=i) {
				bool childFlag=false;
				for(k=1;k<=childIndex1[i];k++) {
					if(children2[i][k]==j) { childFlag=true; break; }
				}
				if(!childFlag) {
					children3[i][++tmp_childIndex]=j;
				}
			}
		}

	}

	//test
/*	for(i=1;i<=seqNum*2-2;i++) {
		fprintf(resultFile, "branch: %4d\n Its right children:",i);
		for(j=1; j<=childIndex[i]; j++) fprintf(resultFile," %4d",children[i][j]);
		fprintf(resultFile, "\nIts left children: ");
		for(j=1; j<=seqNum-childIndex[i]; j++) fprintf(resultFile," %4d",children1[i][j]);
		fprintf(resultFile,"\n");
	}
*/
	//left, right children number for each node
	int *left_child, *right_child;
	left_child=ivector(1,seqNum*2-2);
	right_child=ivector(1,seqNum*2-2);
	for(i=1;i<=seqNum*2-2;i++) {
		right_child[i]=childIndex[i];
		left_child[i]=seqNum-right_child[i];//-1;
	}

	double total_length1=0.0,total_length2=0.0;
	double *left_length1, *right_length1;
	double *left_length2, *right_length2;
//	int **children, *childIndex;
	double *diff_left_length, *diff_right_length;
	left_length1=dvector(1,seqNum*2-2);
	right_length1=dvector(1,seqNum*2-2);
	left_length2=dvector(1,seqNum*2-2);
	right_length2=dvector(1,seqNum*2-2);

	//calculate the total length
	for(i=1;i<=(seqNum*2-2);i++) {
		total_length1+=gapTree[1][i].dist;
		total_length2+=gapTree[2][i].dist;
	}
	fprintf(resultFile, "length for tree1: %8.3f and length for tree2: %8.3f\n",total_length1,total_length2);

	//scale both tree, total_length=10.0
	for(i=1;i<=seqNum*2-2;i++) {
		gapTree[1][i].dist*=10.0/total_length1;
		gapTree[2][i].dist*=10.0/total_length2;
	}
	//assign total length to 10.0
	total_length1=total_length2=10.0;


        //calculate the right_length and left_length
       for(i=1;i<=seqNum*2-2;i++) {
                left_length1[i]=right_length1[i]=0.0;
		left_length2[i]=right_length2[i]=0.0;
        }
        for(i=1;i<=seqNum*2-2;i++) {
                for(j=1;j<=childIndex1[i];j++) {
                        right_length1[i]+=gapTree[1][children2[i][j]].dist;
			right_length2[i]+=gapTree[2][children2[i][j]].dist;
               }
                left_length1[i]=total_length1-gapTree[1][i].dist-right_length1[i];
		left_length2[i]=total_length2-gapTree[2][i].dist-right_length2[i];
        }

	//test: output children information and left & right distance
/*	for(i=1;i<=seqNum*2-2;i++) {
		fprintf(resultFile, "branch: %4d, it's right children: ",i);
		for(j=1;j<=childIndex1[i];j++) fprintf(resultFile, "%4d ",children2[i][j]);
		fprintf(resultFile, "\nIt's left children: ");
		for(j=1;j<=seqNum*2-2-childIndex1[i];j++) fprintf(resultFile, "%4d ",children3[i][j]);
		fprintf(resultFile, "\n");
	}
*/
	for(i=1;i<=seqNum*2-2;i++) 
		fprintf(resultFile, "branch: %4d left dist1: %8.3f right dist1: %8.3f left dist2: %8.3f right dist2: %8.3f \n", i,left_length1[i], right_length1[i],left_length2[i],right_length2[i]);

	//define groups
	int **group, *groupIndex;
	group=imatrix(1,seqNum*2-2,1,seqNum);
	groupIndex=ivector(1,seqNum*2-2);
	for(i=1;i<=seqNum*2-2;i++) {
		groupIndex[i]=0;
		for(j=1;j<=seqNum;j++) group[i][j]=0;
	}

	int unclusterNum=seqNum;
	int *unclusterNode;
	unclusterNode=ivector(1,seqNum);
	for(i=1;i<=seqNum;i++) unclusterNode[i]=1;

/*
**	big partition does not supersize small partition if they overlap
*/
	//try to break tree into groups (maybe one group will contain another group)
	for(int n=1;n<=rankIndex;n++) {
		i=rank[n];
		//test
//		fprintf(resultFile, "Cluster tree based on branch: %4d\n", i);

		if(i==seqNum+1) continue;
		if(unclusterNum<=0) break;

		//calculate T1/T2 and T1'/T2' and choose site with the bigger one
		double ratio_left=0.0, ratio_right=0.0;
		//calculate average distance in subtree
/*		ratio_right=right_length2[i]/(childIndex[i]*2-2);
		ratio_left=left_length2[i]/((seqNum-childIndex[i])*2-2);
*/
		//calculate average pair dist of left and right
		for(j=1;j<childIndex[i];j++) {
			for(k=j+1;k<=childIndex[i];k++) {
				ratio_right+=distMatrix2[children[i][j]][children[i][k]]/(childIndex[i]*(childIndex[i]-1)/2);
			}
		}
		for(j=1;j<seqNum-childIndex[i];j++) {
			for(k=j+1;k<=seqNum-childIndex[i];k++) {
				ratio_left+=distMatrix2[children1[i][j]][children1[i][k]]/((seqNum-childIndex[i])*(seqNum-childIndex[i]-1)/2);
			}
		}
		//test: average branch length
		fprintf(resultFile,"branch: %4d, left: %8.3f num: %4d left_avg: %8.3f right: %8.3f num: %4d right_avg: %8.3f\n", i, left_length2[i],seqNum-childIndex[i], ratio_left, right_length2[i], childIndex[i], ratio_right);

		//check if a group is already a subgroup of big one
		int isSubgroup=0;
		if(ratio_right<ratio_left) {
			for(j=1;j<=childIndex[i];j++) isSubgroup+=unclusterNode[children[i][j]];
			if(isSubgroup<childIndex[i]) { /*fprintf(resultFile, "Branch: %4d is skipped\n",i);*/ continue; }
		}
		else {
			for(j=1;j<=seqNum-childIndex[i];j++) isSubgroup+=unclusterNode[children1[i][j]];
			if(isSubgroup<(seqNum-childIndex[i])) { /*fprintf(resultFile, "Branch: %4d is skipped\n",i);*/ continue; }
		}
/*		if(isSubgroup==0) continue;
		else if(isSubgroup==1) {  //if a new subgroup only one more leaf node than old subgroup(s), make the leaf node singleton
			//find the No. of the singleton
			if(ratio_right<ratio_left) {
				for(j=1;j<=childIndex[i];j++) {
					if(unclusterNode[children[i][j]]==1) {
						int singleton=children[i][j];
						group[singleton][++groupIndex[singleton]]=singleton;
						unclusterNode[singleton]=0;
						unclusterNum--;
						fprintf(resultFile,"Add Singleton: branch: %4d\n",singleton); 
					}
				}
			}
			else if(ratio_right>ratio_left) {
				for(j=1;j<=seqNum-childIndex[i];j++) {
					if(unclusterNode[children1[i][j]]==1) {
						int singleton=children1[i][j];
						group[singleton][++groupIndex[singleton]]=singleton;
						unclusterNode[singleton]=0;
						unclusterNum--;
						fprintf(resultFile,"Add Singleton: branch: %4d\n",singleton); 
					}
				}
			}
			continue;
		}
*/
		if(ratio_right<ratio_left) {
			for(j=1;j<=childIndex[i];j++) {
				group[i][j]=children[i][j];
				if(unclusterNode[group[i][j]]==1) {
					unclusterNode[group[i][j]]=0;
					unclusterNum--;
				}
			}
			groupIndex[i]=childIndex[i];
			
		}
		else if(ratio_right>ratio_left) {
			for(j=1;j<=seqNum-childIndex[i];j++) {
				group[i][j]=children1[i][j];
				if(unclusterNode[group[i][j]]==1) {
					unclusterNode[group[i][j]]=0;
					unclusterNum--;
				}
			}
			groupIndex[i]=seqNum-childIndex[i];
			
		}
		else {
			printf("Warning: Something wrong with clustering\n");
		}
	}

	//output the group information
	for(k=1;k<=rankIndex;k++) {
		i=rank[k];
		if(groupIndex[i]>0) {
			fprintf(resultFile, "group %4d, contains %4d member, ",i,groupIndex[i]);
			for(j=1;j<=groupIndex[i];j++) {
				fprintf(resultFile, "%4d ",group[i][j]);
			}
			fprintf(resultFile, "\n");
		}
	}

	//check if one subgroup is contained by another subgroup
	int *subFlag; subFlag=ivector(1,seqNum*2-2);
	for(i=1;i<=seqNum*2-2;i++) {
		if(groupIndex[i]==0) subFlag[i]=0;
		else subFlag[i]=1;
	}
	for(i=1;i<=seqNum*2-2;i++) {
		bool subgroupFlag=false;
		for(j=1;j<=seqNum*2-2;j++) {
			if(groupIndex[i]>0 && groupIndex[j]>groupIndex[i]) {
				for(k=1;k<=groupIndex[j];k++) {
					if(group[j][k]==group[i][1]) {subgroupFlag=true; break;break;}
				}
			}
		}
		if(subgroupFlag) subFlag[i]=0;
	}

	//output the group information
	fprintf(resultFile,"\n\n\n");
	for(k=1;k<=rankIndex;k++) {
		i=rank[k];
		if(groupIndex[i]>0 && subFlag[i]>0) {
			fprintf(resultFile, "group %4d, contains %4d member, ",i,groupIndex[i]);
			for(j=1;j<=groupIndex[i];j++) {
				fprintf(resultFile, "%4d ",group[i][j]);
			}
			fprintf(resultFile, "\n");
		}
	}
	//check
	int total_group_member=0;
	for(i=1;i<=seqNum*2-2;i++) {
		if(groupIndex[i]>0 && subFlag[i]>0) total_group_member+=groupIndex[i];
	}
	if(total_group_member!=seqNum) fprintf(resultFile,"Warning: total group member is not correct\n");


/*
**	big partition does supersize small partition if they overlap
*/
	fprintf(resultFile,"\n\nBig partition supersizes small partition if they overlap\n");
	//initialize some array
	for(i=1;i<=seqNum*2-2;i++) {
		groupIndex[i]=0;
		for(j=1;j<=seqNum;j++) group[i][j]=0;
	}
	unclusterNum=seqNum;
	for(i=1;i<=seqNum;i++) unclusterNode[i]=1;

	//try to break tree into groups (maybe one group will contain another group)
	for(int n=1;n<=rankIndex;n++) {
		i=rank[n];

		if(i==seqNum+1) continue;
		if(unclusterNum<=0) break;

		//calculate T1/T2 and T1'/T2' and choose site with the bigger one
		double ratio_left=0.0, ratio_right=0.0;
		//calculate average distance in subtree
/*		ratio_right=right_length2[i]/(childIndex[i]*2-2);
		ratio_left=left_length2[i]/((seqNum-childIndex[i])*2-2);
*/
		//calculate average pair dist of left and right
		for(j=1;j<childIndex[i];j++) {
			for(k=j+1;k<=childIndex[i];k++) {
				ratio_right+=distMatrix2[children[i][j]][children[i][k]]/(childIndex[i]*(childIndex[i]-1)/2);
			}
		}
		for(j=1;j<seqNum-childIndex[i];j++) {
			for(k=j+1;k<=seqNum-childIndex[i];k++) {
				ratio_left+=distMatrix2[children1[i][j]][children1[i][k]]/((seqNum-childIndex[i])*(seqNum-childIndex[i]-1)/2);
			}
		}
		//test: average branch length
		fprintf(resultFile,"branch: %4d, left: %8.3f num: %4d left_avg: %8.3f right: %8.3f num: %4d right_avg: %8.3f\n", i, left_length2[i],seqNum-childIndex[i], ratio_left, right_length2[i], childIndex[i], ratio_right);

		//check if a group is already a subgroup of big one
		int isSubgroup=0;
		if(ratio_right<ratio_left) {
			for(j=1;j<=childIndex[i];j++) isSubgroup+=unclusterNode[children[i][j]];
			if(isSubgroup<=childIndex[i] && isSubgroup>0) { 
				for(j=1;j<=childIndex[i];j++) {
					group[i][j]=children[i][j];
					if(unclusterNode[group[i][j]]==1) {
						unclusterNode[group[i][j]]=0;
						unclusterNum--;
					}
				}
				groupIndex[i]=childIndex[i];
			}
			else if(isSubgroup == 0) {
				/*fprintf(resultFile, "Branch: %4d is skipped\n",i);*/ continue; }
		}
		else {
			for(j=1;j<=seqNum-childIndex[i];j++) isSubgroup+=unclusterNode[children1[i][j]];
			if(isSubgroup<=(seqNum-childIndex[i]) && isSubgroup>0) {
				for(j=1;j<=seqNum-childIndex[i];j++) {
					group[i][j]=children1[i][j];
					if(unclusterNode[group[i][j]]==1) {
						unclusterNode[group[i][j]]=0;
						unclusterNum--;
					}
				}
				groupIndex[i]=seqNum-childIndex[i];
			} else if(isSubgroup == 0) {
				/*fprintf(resultFile, "Branch: %4d is skipped\n",i);*/ continue; 
			}
		}
/*
		if(ratio_right<ratio_left) {
			for(j=1;j<=childIndex[i];j++) {
				group[i][j]=children[i][j];
				if(unclusterNode[group[i][j]]==1) {
					unclusterNode[group[i][j]]=0;
					unclusterNum--;
				}
			}
			groupIndex[i]=childIndex[i];
			
		}
		else if(ratio_right>ratio_left) {
			for(j=1;j<=seqNum-childIndex[i];j++) {
				group[i][j]=children1[i][j];
				if(unclusterNode[group[i][j]]==1) {
					unclusterNode[group[i][j]]=0;
					unclusterNum--;
				}
			}
			groupIndex[i]=seqNum-childIndex[i];
			
		}
		else {
			printf("Warning: Something wrong with clustering\n");
		}
*/
	}

	//output the group information
	for(k=1;k<=rankIndex;k++) {
		i=rank[k];
		if(groupIndex[i]>0) {
			fprintf(resultFile, "group %4d, contains %4d member, ",i,groupIndex[i]);
			for(j=1;j<=groupIndex[i];j++) {
				fprintf(resultFile, "%4d ",group[i][j]);
			}
			fprintf(resultFile, "\n");
		}
	}

	//check if one subgroup is contained by another subgroup
//	int *subFlag; subFlag=ivector(1,seqNum*2-2);
	for(i=1;i<=seqNum*2-2;i++) {
		if(groupIndex[i]==0) subFlag[i]=0;
		else subFlag[i]=1;
	}
	for(i=1;i<=seqNum*2-2;i++) {
		bool subgroupFlag=false;
		for(j=1;j<=seqNum*2-2;j++) {
			if(groupIndex[i]>0 && groupIndex[j]>groupIndex[i]) {
				for(k=1;k<=groupIndex[j];k++) {
					if(group[j][k]==group[i][1]) {subgroupFlag=true; break;break;}
				}
			}
		}
		if(subgroupFlag) subFlag[i]=0;
	}

	//output the group information
	fprintf(resultFile,"\n\n\n");
	for(k=1;k<=rankIndex;k++) {
		i=rank[k];
		if(groupIndex[i]>0 && subFlag[i]>0) {
			fprintf(resultFile, "group %4d, contains %4d member, ",i,groupIndex[i]);
			for(j=1;j<=groupIndex[i];j++) {
				fprintf(resultFile, "%4d ",group[i][j]);
			}
			fprintf(resultFile, "\n");
		}
	}
	//check
//	int total_group_member=0;
	total_group_member=0;
	for(i=1;i<=seqNum*2-2;i++) {
		if(groupIndex[i]>0 && subFlag[i]>0) total_group_member+=groupIndex[i];
	}
	if(total_group_member!=seqNum) fprintf(resultFile,"Warning: total group member is not correct\n");


/*
**	new partition will be generated by big partition minus small partition if they overlap
*/
	fprintf(resultFile,"\n\nBig partition minus small partition if they overlap\n");
	//initialize some array
	for(i=1;i<=seqNum*2-2;i++) {
		groupIndex[i]=0;
		for(j=1;j<=seqNum;j++) group[i][j]=0;
	}
	unclusterNum=seqNum;
	for(i=1;i<=seqNum;i++) unclusterNode[i]=1;

	//try to break tree into groups (maybe one group will contain another group)
	for(int n=1;n<=rankIndex;n++) {
		i=rank[n];

		if(i==seqNum+1) continue;
		if(unclusterNum<=0) break;

		//calculate T1/T2 and T1'/T2' and choose site with the bigger one
		double ratio_left=0.0, ratio_right=0.0;
		//calculate average distance in subtree
/*		ratio_right=right_length2[i]/(childIndex[i]*2-2);
		ratio_left=left_length2[i]/((seqNum-childIndex[i])*2-2);
*/
		//calculate average pair dist of left and right
		for(j=1;j<childIndex[i];j++) {
			for(k=j+1;k<=childIndex[i];k++) {
				ratio_right+=distMatrix2[children[i][j]][children[i][k]]/(childIndex[i]*(childIndex[i]-1)/2);
			}
		}
		for(j=1;j<seqNum-childIndex[i];j++) {
			for(k=j+1;k<=seqNum-childIndex[i];k++) {
				ratio_left+=distMatrix2[children1[i][j]][children1[i][k]]/((seqNum-childIndex[i])*(seqNum-childIndex[i]-1)/2);
			}
		}
		//test: average branch length
		fprintf(resultFile,"branch: %4d, left: %8.3f num: %4d left_avg: %8.3f right: %8.3f num: %4d right_avg: %8.3f\n", i, left_length2[i],seqNum-childIndex[i], ratio_left, right_length2[i], childIndex[i], ratio_right);

		//check if a group is already a subgroup of big one
		int isSubgroup=0;
		if(ratio_right<ratio_left) {
			for(j=1;j<=childIndex[i];j++) isSubgroup+=unclusterNode[children[i][j]];
//			if(i==68) printf("branch 68 has %d leaf nodes remaining and its childIndex: %d\n",isSubgroup,childIndex[i]);
			if(isSubgroup<=childIndex[i] && isSubgroup>0) { 
				groupIndex[i]=0;
				for(j=1;j<=childIndex[i];j++) {
					if(unclusterNode[children[i][j]]==1) {
						group[i][++groupIndex[i]]=children[i][j];
						unclusterNode[group[i][groupIndex[i]]]=0;
						unclusterNum--;
					}
				}
			}
			else if(isSubgroup == 0) {
				/*fprintf(resultFile, "Branch: %4d is skipped\n",i);*/ continue; }
		}
		else {
			for(j=1;j<=seqNum-childIndex[i];j++) isSubgroup+=unclusterNode[children1[i][j]];
			if(isSubgroup<=(seqNum-childIndex[i]) && isSubgroup>0) {
				groupIndex[i]=0;
				for(j=1;j<=seqNum-childIndex[i];j++) {
					if(unclusterNode[children1[i][j]]==1) {
						group[i][++groupIndex[i]]=children1[i][j];
						unclusterNode[group[i][groupIndex[i]]]=0;
						unclusterNum--;
					}
				}
			} else if(isSubgroup == 0) {
				/*fprintf(resultFile, "Branch: %4d is skipped\n",i);*/ continue; 
			}
		}

	}

	//output the group information
	for(k=1;k<=rankIndex;k++) {
		i=rank[k];
		if(groupIndex[i]>0) {
			fprintf(resultFile, "group %4d, contains %4d member, ",i,groupIndex[i]);
			for(j=1;j<=groupIndex[i];j++) {
				fprintf(resultFile, "%4d ",group[i][j]);
			}
			fprintf(resultFile, "\n");
		}
	}

	//check if one subgroup is contained by another subgroup
//	int *subFlag; subFlag=ivector(1,seqNum*2-2);
	for(i=1;i<=seqNum*2-2;i++) {
		if(groupIndex[i]==0) subFlag[i]=0;
		else subFlag[i]=1;
	}
	for(i=1;i<=seqNum*2-2;i++) {
		bool subgroupFlag=false;
		for(j=1;j<=seqNum*2-2;j++) {
			if(groupIndex[i]>0 && groupIndex[j]>groupIndex[i]) {
				for(k=1;k<=groupIndex[j];k++) {
					if(group[j][k]==group[i][1]) {subgroupFlag=true; break;break;}
				}
			}
		}
		if(subgroupFlag) subFlag[i]=0;
	}

	//output the group information
	fprintf(resultFile,"\n\n\n");
	for(k=1;k<=rankIndex;k++) {
		i=rank[k];
		if(groupIndex[i]>0 && subFlag[i]>0) {
			fprintf(resultFile, "group %4d, contains %4d member, ",i,groupIndex[i]);
			for(j=1;j<=groupIndex[i];j++) {
				fprintf(resultFile, "%4d ",group[i][j]);
			}
			fprintf(resultFile, "\n");
		}
	}
	//check
//	int total_group_member=0;
	total_group_member=0;
	for(i=1;i<=seqNum*2-2;i++) {
		if(groupIndex[i]>0 && subFlag[i]>0) total_group_member+=groupIndex[i];
	}
	if(total_group_member!=seqNum) fprintf(resultFile,"Warning: total group member is not correct\n");

}

			

void assignTree(char *treeFile1,char *treeFile2)
{
	int i,j,k;
	//read original tree structure
	readTree(treeFile1);
	//store original tree information into gapTree[1];
	for(i=1;i<=(seqNum*2-2);i++) {
		gapTree[1][i].childL=node[i].childL;
		gapTree[1][i].childR=node[i].childR;
		gapTree[1][i].childE=node[i].childE;
		gapTree[1][i].parent=node[i].parent;
		gapTree[1][i].rootFlag=node[i].rootFlag;
		gapTree[1][i].dist=node[i].dist;
	}
	//read new tree structure
	readTree(treeFile2);
	//store new tree information into gapTree[2]
	for(i=1;i<=(seqNum*2-2);i++) {
		gapTree[2][i].childL=node[i].childL;
		gapTree[2][i].childR=node[i].childR;
		gapTree[2][i].childE=node[i].childE;
		gapTree[2][i].parent=node[i].parent;
		gapTree[2][i].rootFlag=node[i].rootFlag;
		gapTree[2][i].dist=node[i].dist;
	}

	//calculate Z score for each branch in original tree ; added by Wei 03/17/2005
	double *Z_score, avg_length=0.0, stdev_length=0.0;
	Z_score = dvector(1,seqNum*2-2);
	for(i=1;i<=(seqNum*2-2);i++) avg_length+=gapTree[1][i].dist/(seqNum*2-3);
	for(i=1;i<=(seqNum*2-2);i++) if(i!=seqNum+1) stdev_length+=(gapTree[1][i].dist-avg_length)*(gapTree[1][i].dist-avg_length);
	stdev_length = sqrt(stdev_length/(seqNum*2-4));
	for(i=1;i<=seqNum*2-2;i++) {
		if(i!=seqNum+1) Z_score[i]=(gapTree[1][i].dist-avg_length)/stdev_length;
		else Z_score[i]=0.0;
	}
	printf("average length: %12.6f and stdev length: %12.6f\n",avg_length, stdev_length);

        //calculate Z score for each branch in optimized tree ; added by Wei 03/17/2005
        double *Z_score1;
	avg_length=0.0, stdev_length=0.0;
        Z_score1 = dvector(1,seqNum*2-2);
        for(i=1;i<=(seqNum*2-2);i++) avg_length+=gapTree[2][i].dist/(seqNum*2-3);
        for(i=1;i<=(seqNum*2-2);i++) if(i!=seqNum+1) stdev_length+=(gapTree[2][i].dist-avg_length)*(gapTree[2][i].dist-avg_length);
        stdev_length = sqrt(stdev_length/(seqNum*2-4));
        for(i=1;i<=seqNum*2-2;i++) {
                if(i!=seqNum+1) Z_score1[i]=(gapTree[2][i].dist-avg_length)/stdev_length;
                else Z_score1[i]=0.0;
        }
        printf("average length: %12.6f and stdev length: %12.6f\n",avg_length, stdev_length);


	//generate branch rank list for optimized tree used for control
	int *nodeIndex; double *lengthIndex;
	nodeIndex = ivector(1,seqNum*2-2);
	lengthIndex = dvector(1,seqNum*2-2);
	for(i=1;i<=seqNum*2-2;i++) {
		nodeIndex[i] = i;
		lengthIndex[i]=gapTree[2][i].dist;
	}
	for(i=1;i<seqNum*2-2;i++) {
		for(j=i+1;j<=seqNum*2-2;j++) {
			if(lengthIndex[i]<lengthIndex[j]) {
				int tmp_nodeIndex = nodeIndex[i];
				nodeIndex[i] = nodeIndex[j];
				nodeIndex[j] = tmp_nodeIndex;
				double tmp_lengthIndex = lengthIndex[i];
				lengthIndex[i] = lengthIndex[j];
				lengthIndex[j] = tmp_lengthIndex;
			}
		}
	}
	fprintf(resultFile,"Branch rank list based on branch length of optimized tree:\n");
	for(i=1;i<=seqNum*2-2;i++) {
		fprintf(resultFile,"%4d %8.4f\n",nodeIndex[i],lengthIndex[i]);
	}

	//calculate the distance matrix for both trees
	cal_dist_matrix();

	double total_length1=0.0,total_length2=0.0;
	double *left_length1, *right_length1;
	double *left_length2, *right_length2;
	int **children, *childIndex;
	double *diff_left_length, *diff_right_length;
	left_length1=dvector(1,seqNum*2-2);
	right_length1=dvector(1,seqNum*2-2);
	left_length2=dvector(1,seqNum*2-2);
	right_length2=dvector(1,seqNum*2-2);
	diff_left_length=dvector(1,seqNum*2-2);
	diff_right_length=dvector(1,seqNum*2-2);
	children=imatrix(1,seqNum*2-2,1,seqNum*2-2);
	childIndex=ivector(1,seqNum*2-2);

	int **children1;
	children1=imatrix(1,seqNum*2-2,1,seqNum*2-2);

	//calculate the total length
	for(i=1;i<=(seqNum*2-2);i++) {
		total_length1+=gapTree[1][i].dist;
		total_length2+=gapTree[2][i].dist;
	}
	//scale both tree, total_length=10.0
/*	for(i=1;i<=seqNum*2-2;i++) {
		gapTree[1][i].dist*=10.0/total_length1;
		gapTree[2][i].dist*=10.0/total_length2;
	}
	//assign total length to 10.0
	total_length1=total_length2=10.0;
*/
	//search children for each node
	for(i=1;i<=seqNum*2-2;i++) childIndex[i]=0;
	for(i=1;i<=seqNum;i++) {
	    if(i!=seqNum+1) {
		int parent=gapTree[1][i].parent;
		while(!gapTree[1][parent].rootFlag) {
			children[parent][++childIndex[parent]]=i;
			parent=gapTree[1][parent].parent;
		}
		children[seqNum+1][++childIndex[seqNum+1]]=i;
	    }
	}

	//add leaf node to its neighbor branch's children set
	//added at Nov.2, 2004
	for(i=1;i<=seqNum;i++) {
		childIndex[i]++;
		children[i][childIndex[i]]=i;
	}

		//assign children1 with child not in children
	for(i=1;i<=seqNum*2-2;i++) {
		int tmp_childIndex=0;
		for(j=1;j<=seqNum;j++) {
			if(j!=i) {
				bool childFlag=false;
				for(k=1;k<=childIndex[i];k++) {
					if(children[i][k]==j) { childFlag=true; break; }
				}
				if(!childFlag) {
					children1[i][++tmp_childIndex]=j;
				}
			}
		}
	}

	//calculate the right_length and left_length
	for(i=1;i<=seqNum*2-2;i++) {
		left_length1[i]=right_length1[i]=left_length2[i]=right_length2[i]=diff_left_length[i]=diff_right_length[i]=0.0;
	}

	//calculate the left and right subtree total branch length
/*	for(i=1;i<=seqNum*2-2;i++) {
		if(i>seqNum) {
			for(j=1;j<=childIndex[i];j++) {
				right_length1[i]+=gapTree[1][children[i][j]].dist;
				right_length2[i]+=gapTree[2][children[i][j]].dist;
			}
		}
		left_length1[i]=total_length1-gapTree[1][i].dist-right_length1[i];
		left_length2[i]=total_length2-gapTree[2][i].dist-right_length2[i];
	}

	for(i=1;i<=seqNum*2-2;i++) {
		printf("branch: %4d it's right children: ",i);
		for(j=1;j<=childIndex[i];j++) printf("%4d ",children[i][j]);
		printf("\nits left children: ");
		for(j=1;j<=seqNum-childIndex[i];j++) printf("%4d ",children1[i][j]);
		printf("\n");
	}
*/
	printf("before calculate the pairwise distance\n");
	//calculate the average pair wise distance for left and right subtree
	for(i=1;i<=seqNum*2-2;i++) {
		for(j=1;j<childIndex[i];j++) {
			for(k=j+1;k<=childIndex[i];k++) {
				right_length1[i]+=distMatrix1[children[i][j]][children[i][k]]/(childIndex[i]*(childIndex[i]-1)/2);
				right_length2[i]+=distMatrix2[children[i][j]][children[i][k]]/(childIndex[i]*(childIndex[i]-1)/2);
			}
		}
		for(j=1;j<seqNum-childIndex[i];j++) {
			for(k=j+1;k<=seqNum-childIndex[i];k++) {
				left_length1[i]+=distMatrix1[children1[i][j]][children1[i][k]]/((seqNum-childIndex[i])*(seqNum-childIndex[i]-1)/2);
				left_length2[i]+=distMatrix2[children1[i][j]][children1[i][k]]/((seqNum-childIndex[i])*(seqNum-childIndex[i]-1)/2);
			}
		}
	}
	printf("finish calculate the pairwise distance\n");

    //output the left_length and right_length, used for test
	fprintf(resultFile,"Branch Length:    tree1      tree2\n");
	for(i=1;i<=seqNum*2-2;i++) {
		fprintf(resultFile,"           %10.4f      %10.4f\n",gapTree[1][i].dist,gapTree[2][i].dist);
	}
	//print out format |Lo |To1 |To2 |Ln |Tn1 |Tn2 (original branch length, original left subtree length, original right subtree length)
	fprintf(resultFile,"No  branch_length1 left_length1   right_length1  branch_length2 left_length2   right_length2\n");
	for(i=1;i<=seqNum*2-2;i++) {
		fprintf(resultFile,"%4d %12.6f %12.6f   %12.6f  %12.6f  %12.6f   %12.6f\n",i,gapTree[1][i].dist,left_length1[i],right_length1[i],gapTree[2][i].dist,left_length2[i],right_length2[i]);
	}
	//for every branch in gapTree[2], make it equal to corresponding branch in gapTree[1], then calculate the diff in left and righ length
	//if branch length in gapTree[2] equals to 0.0, then assign a very bigger value to diff_left and diff_right length
	for(i=1;i<=seqNum*2-2;i++) {
		if(gapTree[2][i].dist<1.0e-8) { diff_left_length[i]=1.0e6; diff_right_length[i]=1.0e6; }
		else {
			double ratio=gapTree[1][i].dist/gapTree[2][i].dist;
			diff_left_length[i]=left_length2[i]*ratio-left_length1[i];
			diff_right_length[i]=right_length2[i]*ratio-right_length1[i];
		}
	}
	//output result
	fprintf(resultFile,"  No  diff_left_length     diff_right_length     children\n");
	for(i=1;i<=seqNum*2-2;i++) {
		fprintf(resultFile,"%4d  %16.6f    %16.6f    ",i,diff_left_length[i],diff_right_length[i]);
		for(j=1;j<=childIndex[i];j++) fprintf(resultFile,"%3d ",children[i][j]);//node[children[i][j]].nodeName);
		fprintf(resultFile,"\n");
	}
	//calculate the sum of diff_left_length and diff_right_length, then sort them by ascend
	double *total_diff;
	int *index_total_diff;
	total_diff=dvector(1,seqNum*2-2);
	index_total_diff=ivector(1,seqNum*2-2);
	for(i=1;i<=seqNum*2-2;i++) {
		index_total_diff[i]=i;
		total_diff[i]=diff_left_length[i]+diff_right_length[i];
	}
	for(i=1;i<seqNum*2-2;i++) {
		for(j=i+1;j<=seqNum*2-2;j++) {
			if(total_diff[i]>total_diff[j]) {
				double tmp_swap=total_diff[i];
				total_diff[i]=total_diff[j]; total_diff[j]=tmp_swap;
				int tmp_index=index_total_diff[i];
				index_total_diff[i]=index_total_diff[j];
				index_total_diff[j]=tmp_index;
			}
		}
	}
	fprintf(resultFile,"\n\n  No oiginLength newLength diff_left_length     diff_right_length   total_diff   children\n");
	for(i=1;i<=seqNum*2-2;i++) {
		int tmp_index=index_total_diff[i];
		fprintf(resultFile,"%4d %10.6f %10.6f %16.6f  %16.6f %16.6f   ",tmp_index,gapTree[1][tmp_index].dist,gapTree[2][tmp_index].dist*total_length1/total_length2,diff_left_length[tmp_index],diff_right_length[tmp_index],total_diff[i]);
		for(j=1;j<=childIndex[tmp_index];j++) fprintf(resultFile,"%s ",node[children[tmp_index][j]].nodeName);
		fprintf(resultFile,"\n");
	}

	//calculate the target function max(T1/T2,T1'/T2')*(Ln-Lo)
	double *target_value; target_value=dvector(1,seqNum*2-2);
	int *index_target; index_target=ivector(1,seqNum*2-2);
	for(i=1;i<=seqNum*2-2;i++) {
		index_target[i]=i;
		target_value[i]=gapTree[2][i].dist-gapTree[1][i].dist;
		//Nov.11, 2004  modify target function to max(T1/T2,T1'/T2')*(Ln-Lo)*Lo
//		target_value[i]=(gapTree[2][i].dist-gapTree[1][i].dist)*gapTree[1][i].dist;
		//Mar.17, 2005  modify target function to max(T1/T2,T1'/T2')*(Ln-Lo)*Zo
//		target_value[i]=(gapTree[2][i].dist-gapTree[1][i].dist)*Z_score[i]*Z_score[i];
//		target_value[i]=(gapTree[2][i].dist-gapTree[1][i].dist)*Z_score[i]*Z_score1[i];
		double left=0.0,right=0.0;
		if(left_length1[i]<1.0e-8 && left_length2[i]<1.0e-8) left=1.0;
		else if(left_length2[i]<1.0e-8) left=1.0e10;
		else left=left_length1[i]/left_length2[i];
		if(right_length1[i]<1.0e-8 && right_length2[i]<1.0e-8) right=1.0;
		else if(right_length2[i]<1.0e-8) right=1.0e10;
		else right=right_length1[i]/right_length2[i];
		if(left>=right) target_value[i]*=left;
		else target_value[i]*=right;
	}
	fprintf(resultFile,"Target1: ");
	for(i=1;i<=seqNum*2-2;i++) fprintf(resultFile,"%12.6f ",target_value[i]);
	fprintf(resultFile,"\n");

	for(i=1;i<seqNum*2-2;i++) {
		for(j=i+1;j<=seqNum*2-2;j++) {
			if(target_value[i]<target_value[j]) {
				double tmp_target=target_value[j];
				target_value[j]=target_value[i]; target_value[i]=tmp_target;
				int tmp_index_target=index_target[j];
				index_target[j]=index_target[i]; index_target[i]=tmp_index_target;
			}
		}
	}
	fprintf(resultFile,"No  branch_length1 left_length1   right_length1  branch_length2 left_length2   right_length2 Z_score Z_score1 target_value\n");
	for(j=1;j<=seqNum*2-2;j++) {
		i=index_target[j];
		fprintf(resultFile,"%4d %12.6f %12.6f   %12.6f  %12.6f  %12.6f   %12.6f  %12.6f  %12.6f  %12.6f\n",i,gapTree[1][i].dist,left_length1[i],right_length1[i],gapTree[2][i].dist,left_length2[i],right_length2[i], Z_score[i], Z_score1[i], target_value[j]);
	}
	

	//calculate the target function max((T1-T2),(T1'-T2'))*(Ln-Lo)
	double *target_value1; target_value1=dvector(1,seqNum*2-2);
	for(i=1;i<=seqNum*2-2;i++) {
		if((left_length1[i]-left_length2[i])>(right_length1[i]-right_length2[i])) target_value1[i]=left_length1[i]-left_length2[i];
		else target_value1[i]=right_length1[i]-right_length2[i];
		target_value1[i]*=(gapTree[2][i].dist-gapTree[1][i].dist);
	}
	fprintf(resultFile,"Target2: ");
	for(i=1;i<=seqNum*2-2;i++) fprintf(resultFile,"%12.6f ",target_value1[i]);
	fprintf(resultFile,"\n");

}

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)
		printf("Data file %s not found\n",treeFile);
	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");

	//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");
	printf("Finish reading tree\n");
}


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=new topology * [seqLength+2];
	for(i=0;i<seqLength+2;i++) gapTree[i]=new topology [seqNum*2-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]='-';}
		}
	}
	//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);
	
	printf("Finish reading alignment\n");
}



