#include "ancestral.h"

void chooseSeq(char *listFile, int number);
void chooseSeq1(char *listFile, int number);
void histgram(char *seqFile, int lowBoundItem, int pickItem, char *alnFile);

struct IDENTITY {
	int LPoint;
	int RPoint;
	double identity;
	bool LClusterFlag;
	bool RClusterFlag;
};


NODE *node;
FILE *resultFile=stdout;
int seqNum, seqLength;
char aa[NP+1]={'-','A','R','N','D','C','Q','E','G','H','I','L','K','M','F','P','S','T','W','Y','V'};
long idum=(-64);
void usage()
{
//	printf("Usage: chooseSeq -i sequencefile -l listfile -n sequenceNumber\n");
	printf("Usage: chooseSeq [-H] [-S] -i sequencefile -L (1~7) -C (1~4) -o output alignment\n");
	exit(1);
}


int main(int argc, char **argv)
{
	int i;
	int site;
	FILE *input;
	char seqFile[20], treeFile[20], alnFile[40];
	bool histFlag=false,seqFlag=false;
	int lowBoundItem=0, pickItem=0;

	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 'H':
				histFlag=true;
				break;
			case 'S':
				seqFlag=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 'L':
				++i;
				optidx=0;
				if(i<argc) {
					lowBoundItem=atoi(argv[i]);
					if(lowBoundItem<1 || lowBoundItem>7) {printf("L parameter is wrong\n"); usage(); }
				}
				else usage();
				break;
			case 'C':
				++i;
				optidx=0;
				if(i<argc) {
					pickItem=atoi(argv[i]);
					if(pickItem<1 || pickItem>4) {printf("C parameter is wrong\n"); usage(); }
				}
				else usage();
				break;
			case 'o':
				++i;
				optidx=0;
				if(i<argc){
					strcpy(alnFile,argv[i]);
				}
				else usage();
				break;
			case 'l':
				++i;
				optidx=0;
				if(i<argc) {
					strcpy(treeFile,argv[i]);
					input=fopen(treeFile,"r");
					if(!input) {
						printf("can not open file:%s\n",treeFile);
						exit(1);
					}
					fclose(input);
				}
				else usage();
				break;
			case 'n':
				++i;
				optidx=0;
				if(i<argc) {
					site=atoi(argv[i]);
					if(site<=0) {
						printf("site should be greater than 0\n");
						exit(1);
					}
				}
				else usage();
				break;
			default:
				usage();
		}
	}

	//read S_ij matrix
	readali(seqFile);
	if(histFlag) histgram(seqFile,lowBoundItem,pickItem,alnFile);
	if(seqFlag) chooseSeq1(treeFile,site);
}


void histgram(char *seqFile, int lowBoundItem, int pickItem, char *alnFile)
{
	int i,j,k;
	double **idMatrix;
	idMatrix=dmatrix(1,seqNum,1,seqNum);
//	char outFile[20],outFile1[20];
	char recFile[50];
	strcpy(recFile,alnFile);
//	recFile[strlen(recFile)-3]='\0';
	for(i=strlen(recFile);i>=0;i--) if(recFile[i]=='.') { recFile[i+1]='\0'; break; }
	strcat(recFile,"rec");
	FILE *output,*output1;
/*
	for(i=0;i<20;i++){
		if(seqFile[i]!='.') {
			outFile[i]=seqFile[i];
		}
		else { outFile[i]='\0'; break;}
	}
	strcpy(outFile1,outFile);
	strcat(outFile,"_1.aln");
	strcat(outFile1,"_2.aln");
*/
	output1=fopen(recFile,"w");

	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[i].seq_aa[k]!='.')&&(node[j].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;
		}
	}
	for(i=1;i<=seqNum;i++) idMatrix[i][i]=0.0;
	for(i=1;i<seqNum;i++) for(j=i+1;j<=seqNum;j++) idMatrix[j][i]=idMatrix[i][j];
	int *binCount;
	binCount=ivector(1,20);
	for(i=1;i<=20;i++) binCount[i]=0;
	int tmpCount=0;
	for(i=1;i<seqNum;i++){
		for(j=i+1;j<=seqNum;j++){
//			printf("%3.0f\n",idMatrix[i][j]*100.0);
			tmpCount++;
			for(k=0;k<=19;k++){
				if((idMatrix[i][j]>=k*0.05)&&(idMatrix[i][j]<(k+1)*0.05)) {
					binCount[k+1]++;
					break;
				}
				else if(idMatrix[i][j]>0.95) {
					binCount[20]++;
					break;
				}
//				if(idMatrix[i][j]<0.0)  printf("Warning\n");
			}
		}
	}
	fprintf(output1,"seqNum: %d  total number of pairs: %d  tmpCount:%d\n",seqNum,seqNum*(seqNum-1)/2,tmpCount);
	for(i=1;i<=20;i++) fprintf(output1,"%5d\n",binCount[i]);

	//define struct IDENTITY
/*	IDENTITY *identity;
	identity = new IDENTITY [seqNum*(seqNum-1)/2+1];
	int *cluster,clusterIndex=1;
	cluster=ivector(1,seqNum);
	for(i=1;i<=seqNum;i++) cluster[i]=0;
	//get all identity pair
	tmpCount=1;
	for(i=1;i<seqNum;i++){
		for(j=i+1;j<=seqNum;j++){
			identity[tmpCount].LPoint=i;
			identity[tmpCount].RPoint=j;
			identity[tmpCount].identity=idMatrix[i][j];
			identity[tmpCount].LClusterFlag=false;
			identity[tmpCount].RClusterFlag=false;
			tmpCount++;
		}
	}
	//sort the identity vector
	for(i=1;i<seqNum*(seqNum-1)/2;i++){
		for(j=i+1;j<=seqNum*(seqNum-1)/2;j++){
			if(identity[i].identity<identity[j].identity){
				IDENTITY tmpID;
				tmpID=identity[j]; identity[j]=identity[i]; identity[i]=tmpID;
			}
		}
	}
*/
	//find the highBound and lowBound
//	double highBound=identity[(int) (0.1*seqNum*(seqNum-1)/2)].identity;
//	double lowBound=identity[(int) (0.15*seqNum*(seqNum-1)/2)].identity;
	double highBound=0.8;
	double lowBound=0.0;
	for(i=1;i<seqNum;i++) for(j=i+1;j<=seqNum;j++) lowBound+=idMatrix[i][j];
	lowBound=lowBound/(seqNum*(seqNum-1)/2);
	if(lowBoundItem==1) lowBound=0.25;
	else if(lowBoundItem==2) lowBound=0.30;
	else if(lowBoundItem==3) lowBound=0.35;
	else if(lowBoundItem==4) lowBound=0.40;
	else if(lowBoundItem==6) lowBound+=0.05;
	else if(lowBoundItem==7) lowBound-=0.05;
//	if(lowBound<0.3) lowBound=0.3;
//	double biggest=identity[1].identity;
//	int row=identity[1].LPoint;
//	int col=identity[1].RPoint;
	double biggest=0.0; int row,col;
	for(i=1;i<seqNum;i++) for(j=i+1;j<=seqNum;j++) if(biggest<idMatrix[i][j]) {biggest=idMatrix[i][j]; row=i; col=j; }
	int *cluster,clusterIndex=1;
	cluster=ivector(1,seqNum);
	for(i=1;i<=seqNum;i++) cluster[i]=0;

	//prepare matrix of identity
	double **idCopy;
	idCopy=dmatrix(1,seqNum,1,seqNum);
	for(i=1;i<seqNum;i++) for(j=i+1;j<=seqNum;j++) idCopy[i][j]=idMatrix[i][j];
	for(i=1;i<seqNum;i++) for(j=i+1;j<=seqNum;j++) idCopy[j][i]=idCopy[i][j];
	for(i=1;i<=seqNum;i++) idCopy[i][i]=0.0;

	//test
/*	for(i=1;i<=seqNum;i++){
		for(j=1;j<=seqNum;j++) printf("%3d ",(int) (idCopy[i][j]*1000.0));
		printf("\n");
	}
*/
	fprintf(output1,"highBound: %3d lowBound: %3d\n",(int) (1000.0*highBound),(int) (1000.0*lowBound));
	//cluster with highBound, method: single linkage
	while(biggest>=highBound){
//		printf("biggest: %3d  clusterIndex:%d\n",(int) (biggest*1000.0),clusterIndex);
		if((cluster[row]==0)&&(cluster[col]==0)){
			cluster[row]=cluster[col]=clusterIndex;
			clusterIndex++;
		}
		else if(cluster[row]==0) cluster[row]=cluster[col];
		else if(cluster[col]==0) cluster[col]=cluster[row];
		else {
			int changeFrom,changeTo;
			if(cluster[row]<cluster[col]) {changeFrom=cluster[col]; changeTo=cluster[row]; }
			else { changeFrom=cluster[row]; changeTo=cluster[col]; }
			if(changeFrom!=changeTo){
//				printf("cluster:%4d and cluster:%4d combine together\n",changeFrom,changeTo);

				for(i=1;i<=seqNum;i++){
					if(cluster[i]==changeFrom) cluster[i]=changeTo;
				}

				for(i=1;i<=seqNum;i++) if(cluster[i]>changeFrom) cluster[i]--;
				clusterIndex--;
			}

		}
		//change the matrix
		idCopy[row][col]=idCopy[col][row]=0.0;
/*
		for(i=1;i<row;i++){
			if(idCopy[i][row]<idCopy[i][col]) idCopy[i][col]=0.0;
			else { idCopy[i][row]=idCopy[i][col]; idCopy[i][col]=0.0; }
		}
		for(i=col+1;i<=seqNum;i++){
			if(idCopy[row][i]<idCopy[col][i]) idCopy[col][i]=0.0;
			else { idCopy[row][i]=idCopy[col][i]; idCopy[col][i]=0.0; }
		}
		for(i=row+1;i<=col-1;i++){
			if(idCopy[row][i]<idCopy[i][col]) idCopy[i][col]=0.0;
			else { idCopy[row][i]=idCopy[i][col]; idCopy[i][col]=0.0; }
		}
*/
		//test
/*		printf("row:%3d  col:%3d  biggest:%3d\n",row,col,(int)(biggest*1000.0));
		for(i=1;i<=seqNum;i++){
			for(j=1;j<=seqNum;j++) printf("%3d ",(int) (idCopy[i][j]*1000.0));
			printf("\n");
		}
		printf("\n");
*/
		//find the biggest
		biggest=0.0; row=col=0;
		for(i=1;i<seqNum;i++){
			for(j=i+1;j<=seqNum;j++){
				if(idCopy[i][j]>biggest) {
					biggest=idCopy[i][j];
					row=i; col=j;
				}
			}
		}
	}
	//test
	fprintf(output1,"clusterIndex:%d\n",clusterIndex);
	for(i=0;i<clusterIndex;i++){
		fprintf(output1,"Cluster %d : ",i);
		for(j=1;j<=seqNum;j++)	if(cluster[j]==i) fprintf(output1,"%d ",j);
		fprintf(output1,"\n");
	}
	int clusterNumber=clusterIndex-1;
	for(j=1;j<=seqNum;j++) if(cluster[j]==0) clusterNumber++;
	fprintf(output1,"Total cluster number: %d\n",clusterNumber);

	//pick sequences with longest non-gap length
	//define cluster
	int **group;
	if(clusterIndex<500) group=imatrix(1,500,1,seqNum);
	else group=imatrix(1,clusterIndex-1,1,seqNum);
	int *groupIndex;
	if(clusterIndex<500) groupIndex=ivector(1,500);
	else groupIndex=ivector(1,clusterIndex-1);
	for(i=1;i<clusterIndex;i++) groupIndex[i]=1;
	for(i=1;i<clusterIndex;i++){
		for(j=1;j<=seqNum;j++){
			if(cluster[j]==i) group[i][groupIndex[i]++]=j;
		}
	}
	//find the site gap flag;
	bool *siteFlag;
	siteFlag=new bool [seqLength+1];
	for(i=1;i<=seqLength;i++){
		int count=0;
		for(j=1;j<=seqNum;j++){
			if(node[j].seq_aa[i]=='.'||node[j].seq_aa[i]=='-') count++;
		}
		if(2*count>seqNum) siteFlag[i]=true;
		else siteFlag[i]=false;
	}
	//pick up sequence
	for(i=1;i<clusterIndex;i++){
		int *groupCount;
		groupCount=ivector(1,groupIndex[i]-1);
		for(j=1;j<groupIndex[i];j++){
			groupCount[j]=0;
			for(k=1;k<=seqLength;k++){
				if(!siteFlag[k]){
					if(node[group[i][j]].seq_aa[k]!='.' && node[group[i][j]].seq_aa[k]!='-')
						groupCount[j]++;
				}
			}
//			printf("%d ",groupCount[j]);
		}
		//find the longest
		int longest=0; int longestIndex=0;
		for(j=1;j<groupIndex[i];j++) if(longest<groupCount[j]) {longest=groupCount[j]; longestIndex=j; }
		if(longestIndex==0) { printf("warning: there is wrong with picking sequence from group: %d\n",i); exit(0); }
		group[i][groupIndex[i]]=group[i][longestIndex];
//		printf("\n longest:%d\n",group[i][groupIndex[i]]);
		free_ivector(groupCount,1,groupIndex[i]-1);
	}
	//find maximum, minimum and average of identity for each group
	double *highMax,*highMin,*highAvg;
	highMax=dvector(1,clusterIndex-1);
	highMin=dvector(1,clusterIndex-1);
	highAvg=dvector(1,clusterIndex-1);
	for(i=1;i<clusterIndex;i++) { highMax[i]=0.0; highMin[i]=100.0; highAvg[i]=0.0; }
	for(i=1;i<clusterIndex;i++){
		for(j=1;j<groupIndex[i]-1;j++){
			for(k=j+1;k<groupIndex[i];k++){
				highAvg[i]+=idMatrix[group[i][j]][group[i][k]];
				if(highMax[i]<idMatrix[group[i][j]][group[i][k]]) highMax[i]=idMatrix[group[i][j]][group[i][k]];
				if(highMin[i]>idMatrix[group[i][j]][group[i][k]]) highMin[i]=idMatrix[group[i][j]][group[i][k]];
			}
		}
		highAvg[i]=highAvg[i]/((groupIndex[i]-1)*(groupIndex[i]-2)/2.0);
	}
	for(i=1;i<clusterIndex;i++){
		fprintf(output1,"cluster:%3d Max ID:%4.0f  Min ID:%4.0f Avg ID:%4.0f\n",i,highMax[i]*1000.0,highMin[i]*1000.0,highAvg[i]*1000.0);
	}

	//copy identity matrix
	for(i=1;i<seqNum;i++) for(j=i+1;j<=seqNum;j++) idCopy[i][j]=idMatrix[i][j];
	//change cluster vector and idCopy
	for(i=1;i<clusterIndex;i++){
		for(j=1;j<=groupIndex[i]-1;j++){
			if(group[i][j]!=group[i][groupIndex[i]]){
				cluster[group[i][j]]=-1;
				for(k=1;k<=group[i][j]-1;k++) idCopy[k][group[i][j]]=0.0;
				for(k=group[i][j]+1;k<=seqNum;k++) idCopy[group[i][j]][k]=0.0;
			}
			else cluster[group[i][j]]=0;
		}
	}
	//cluster with low bound
	//find the biggest
	biggest=0.0; row=col=0;
	for(i=1;i<seqNum;i++){
		for(j=i+1;j<=seqNum;j++) if(idCopy[i][j]>biggest){ biggest=idCopy[i][j]; row=i;col=j; }
	}
	clusterIndex=1;
	//cluster with lowBound, method: complete linkage
	while(biggest>=lowBound){
		if((cluster[row]==0)&&(cluster[col]==0)){
			cluster[row]=cluster[col]=clusterIndex;
			clusterIndex++;
		}
		else if(cluster[row]==0) cluster[row]=cluster[col];
		else if(cluster[col]==0) cluster[col]=cluster[row];
		else {
			int changeFrom,changeTo;
			if(cluster[row]<cluster[col]) {changeFrom=cluster[col]; changeTo=cluster[row]; }
			else { changeFrom=cluster[row]; changeTo=cluster[col]; }
			if(changeFrom!=changeTo){
//				printf("cluster:%4d and cluster:%4d combine together\n",changeFrom,changeTo);

				for(i=1;i<=seqNum;i++){
					if(cluster[i]==changeFrom) cluster[i]=changeTo;
				}

				for(i=1;i<=seqNum;i++) if(cluster[i]>changeFrom) cluster[i]--;
				clusterIndex--;
			}

		}
		//change the matrix
		idCopy[row][col]=0.0;

		for(i=1;i<row;i++){
			if(idCopy[i][row]<idCopy[i][col]) idCopy[i][col]=0.0;
			else { idCopy[i][row]=idCopy[i][col]; idCopy[i][col]=0.0; }
		}
		for(i=col+1;i<=seqNum;i++){
			if(idCopy[row][i]<idCopy[col][i]) idCopy[col][i]=0.0;
			else { idCopy[row][i]=idCopy[col][i]; idCopy[col][i]=0.0; }
		}
		for(i=row+1;i<=col-1;i++){
			if(idCopy[row][i]<idCopy[i][col]) idCopy[i][col]=0.0;
			else { idCopy[row][i]=idCopy[i][col]; idCopy[i][col]=0.0; }
		}

		//test
/*		printf("row:%3d  col:%3d  biggest:%3d\n",row,col,(int)(biggest*1000.0));
		for(i=1;i<=seqNum;i++){
			for(j=1;j<=seqNum;j++) printf("%3d ",(int) (idCopy[i][j]*1000.0));
			printf("\n");
		}
		printf("\n");
*/
		//find the biggest
		biggest=0.0; row=col=0;
		for(i=1;i<seqNum;i++){
			for(j=i+1;j<=seqNum;j++){
				if(idCopy[i][j]>biggest) {
					biggest=idCopy[i][j];
					row=i; col=j;
				}
			}
		}
	}

	//statistics analysis
	for(i=1;i<clusterIndex;i++) groupIndex[i]=1;
	for(i=1;i<clusterIndex;i++){
		for(j=1;j<=seqNum;j++){
			if(cluster[j]==i) group[i][groupIndex[i]++]=j;
		}
	}
	for(i=1;i<clusterIndex;i++) { highMax[i]=0.0; highMin[i]=100.0; highAvg[i]=0.0; }
	for(i=1;i<clusterIndex;i++){
		for(j=1;j<groupIndex[i]-1;j++){
			for(k=j+1;k<groupIndex[i];k++){
				highAvg[i]+=idMatrix[group[i][j]][group[i][k]];
				if(highMax[i]<idMatrix[group[i][j]][group[i][k]]) highMax[i]=idMatrix[group[i][j]][group[i][k]];
				if(highMin[i]>idMatrix[group[i][j]][group[i][k]]) highMin[i]=idMatrix[group[i][j]][group[i][k]];
			}
		}
		highAvg[i]=highAvg[i]/((groupIndex[i]-1)*(groupIndex[i]-2)/2.0);
	}
	for(i=1;i<clusterIndex;i++){
		fprintf(output1,"cluster:%3d Max ID:%4.0f  Min ID:%4.0f Avg ID:%4.0f\n",i,highMax[i]*1000.0,highMin[i]*1000.0,highAvg[i]*1000.0);
	}

	//calculate the max & avg identity between each two groups
	double **maxIDGroups, **avgIDGroups;
	maxIDGroups=dmatrix(1,clusterIndex-1,1,clusterIndex-1);
	avgIDGroups=dmatrix(1,clusterIndex-1,1,clusterIndex-1);
	for(i=1;i<clusterIndex-1;i++) {
		for(j=i+1;j<=clusterIndex-1;j++) {
			maxIDGroups[i][j]=0.0; avgIDGroups[i][j]=0.0;
			for(k=1;k<groupIndex[i];k++) {
				for(int m=1;m<groupIndex[j];m++) {
					if(maxIDGroups[i][j]<idMatrix[group[i][k]][group[j][m]]) maxIDGroups[i][j]=idMatrix[group[i][k]][group[j][m]];
					avgIDGroups[i][j]+=idMatrix[group[i][k]][group[j][m]]/(groupIndex[i]-1)/(groupIndex[j]-1);
				}
			}
		}
	}
	for(i=1;i<clusterIndex-1;i++) {
		maxIDGroups[i][i]=avgIDGroups[i][i]=0.0;
		for(j=i+1;j<=clusterIndex-1;j++) {
			maxIDGroups[j][i]=maxIDGroups[i][j];
			avgIDGroups[j][i]=avgIDGroups[i][j];
		}
	}
	fprintf(output1,"MAX identities between each two clusters\n");
	fprintf(output1,"Cluster: ");
	for(i=1;i<=clusterIndex-1;i++) fprintf(output1,"%4d ",i);
	fprintf(output1,"\n");
	for(i=1;i<=clusterIndex-1;i++) {
		fprintf(output1,"%8d ",i);
		for(j=1;j<=clusterIndex-1;j++) fprintf(output1,"%4.0f ",maxIDGroups[i][j]*1000.0);
		fprintf(output1,"\n");
	}
	fprintf(output1,"AVG identities between each two clusters\n");
	fprintf(output1,"Cluster: ");
	for(i=1;i<=clusterIndex-1;i++) fprintf(output1,"%4d ",i);
	fprintf(output1,"\n");
	for(i=1;i<=clusterIndex-1;i++) {
		fprintf(output1,"%8d ",i);
		for(j=1;j<=clusterIndex-1;j++) fprintf(output1,"%4.0f ",avgIDGroups[i][j]*1000.0);
		fprintf(output1,"\n");
	}

	//histogram of number of members in each cluster
	//find the largest number
	int largest=0,largest0;
	for(i=1;i<=seqNum;i++) if(cluster[i]==0) largest++;
	largest0=largest;
	for(i=1;i<clusterIndex;i++){
		if(largest<groupIndex[i]-1) largest=groupIndex[i]-1;
	}
	int *lowBin;
	lowBin=ivector(1,largest);
	lowBin[1]=largest0;
	for(i=2;i<=largest;i++) lowBin[i]=0;
	for(j=1;j<clusterIndex;j++) lowBin[groupIndex[j]-1]++;
//	for(i=1;i<=largest;i++) printf("%d\n",lowBin[i]*i);

	//test
	int *clusterCount; clusterCount=ivector(1,clusterIndex-1);
	int *clusterCountIndex; clusterCountIndex=ivector(1,clusterIndex-1);
	for(i=1;i<clusterIndex;i++) {clusterCount[i]=0; clusterCountIndex[i]=i;}
	for(i=1;i<=seqNum;i++) if(cluster[i]>=0) clusterCount[cluster[i]]++;
	for(i=1;i<clusterIndex-1;i++) for(j=i+1;j<clusterIndex;j++) if(clusterCount[i]<clusterCount[j])
	{	int tmpCluster=clusterCount[i]; int tmpClusterIndex=clusterCountIndex[i];
		clusterCount[i]=clusterCount[j]; clusterCountIndex[i]=clusterCountIndex[j];
		clusterCount[j]=tmpCluster; clusterCountIndex[j]=tmpClusterIndex;
	}

	for(i=1;i<clusterIndex;i++) fprintf(output1,"%d\n",clusterCount[i]);
	int outputScope;
	if(pickItem==1) outputScope=2;
	else if(pickItem==2) outputScope=3;
	else if(pickItem==3) outputScope=4;
	else if(pickItem==4) outputScope=5;
	int outputScopeCopy=outputScope;
	if(outputScope>=clusterIndex-1) outputScope=clusterIndex-1;
	else {
		while(clusterCount[outputScope]==clusterCount[outputScope+1]) {
			outputScope++;
			if(outputScope==clusterIndex-1) break;
		}
	}

	fprintf(output1,"clusterIndex:%d\n",clusterIndex);
	for(i=1;i<clusterIndex;i++){
		fprintf(output1,"Cluster %3d : ",clusterCountIndex[i]);
		for(j=1;j<=seqNum;j++)	if(cluster[j]==clusterCountIndex[i]) fprintf(output1,"%3d ",j);
		fprintf(output1,"\n");
	}
	for(j=1;j<=seqNum;j++) if(cluster[j]==0) fprintf(output1,"Cluster   0: %3d\n",j);
//	clusterNumber=clusterIndex-1;
	if(outputScope>=outputScopeCopy) for(j=1;j<=seqNum;j++) if(cluster[j]==0) clusterNumber--;
	fprintf(output1,"Total cluster number: %d\n",clusterNumber);

	//redefine siteFlag
	for(i=1;i<=seqLength;i++) siteFlag[i]=false;
	for(j=1;j<=seqLength;j++){
		int count=0;
		for(i=1;i<=seqNum;i++){
			if(cluster[i]>0){
				if(node[i].seq_aa[j]=='.' || node[i].seq_aa[j]=='-') count++;
			}
		}
		if(2*count>clusterNumber) siteFlag[j]=true;
	}
//	strcat(outFile1,"_2.aln");
//	output1=fopen(outFile1,"w");
	fprintf(output1,"%20s ",node[0].nodeName);
	for(j=1;j<=seqLength;j++)
		fprintf(output1,"%c",node[0].seq_aa[j]);
	fprintf(output1,"\n");

	for(i=1;i<=clusterIndex-1;i++){
//		for(j=1;j<groupIndex[clusterCountIndex[i]];j++){
		for(j=1;j<groupIndex[i];j++) {
			fprintf(output1,"%20s ",node[group[i][j]].nodeName);
			for(k=1;k<=seqLength;k++)
				fprintf(output1,"%c",node[group[i][j]].seq_aa[k]);
			fprintf(output1,"\n");
		}
		fprintf(output1,"\n");
	}
/*	for(i=1;i<=seqNum;i++){
		if(cluster[i]>0){
			fprintf(output1,"%20s ",node[i].nodeName);
			for(j=1;j<=seqLength;j++)
				fprintf(output1,"%c",node[i].seq_aa[j]);
			fprintf(output1,"\n");
		}
	}
*/
	fclose(output1);
	//output *_1.aln
	output=fopen(alnFile,"w");
	for(i=1;i<=outputScope;i++){
		for(j=1;j<groupIndex[clusterCountIndex[i]];j++){
			char tmpChar[20];
			if(strlen(node[group[clusterCountIndex[i]][j]].nodeName)>14) node[group[clusterCountIndex[i]][j]].nodeName[15]='\0';
			strcat(node[group[clusterCountIndex[i]][j]].nodeName,"$");
			sprintf(tmpChar,"%d",clusterCountIndex[i]);
			strcat(node[group[clusterCountIndex[i]][j]].nodeName,tmpChar);
			fprintf(output,"%20s ",node[group[clusterCountIndex[i]][j]].nodeName);
			for(k=1;k<=seqLength;k++){
				if(node[0].seq_aa[k]!='-'){
					if(node[group[clusterCountIndex[i]][j]].seq_aa[k]=='.')
						fprintf(output,"%c",'-');
					else fprintf(output,"%c",node[group[clusterCountIndex[i]][j]].seq_aa[k]);
				}
			}
			fprintf(output,"\n");
		}
	}
	if(outputScope<outputScopeCopy){
		for(i=1;i<=seqNum;i++){
			if(cluster[i]==0){
				if(strlen(node[i].nodeName)>14) node[i].nodeName[15]='\0';
				strcat(node[i].nodeName,"$0");
				fprintf(output,"%20s ",node[i].nodeName);
				for(k=1;k<=seqLength;k++){
					if(node[0].seq_aa[k]!='-'){
						if(node[i].seq_aa[k]=='.')
							fprintf(output,"%c",'-');
						else fprintf(output,"%c",node[i].seq_aa[k]);
					}
				}
				fprintf(output,"\n");
			}
		}
	}

/*
	for(i=1;i<=seqNum;i++){
		if(cluster[i]>0){
			if(strlen(node[i].nodeName)>14) node[i].nodeName[15]='\0';
			strcat(node[i].nodeName,"$");
			strcat(node[i].nodeName,itoa(cluster[i]));
			fprintf(output,"%20s ",node[i].nodeName);
			for(j=1;j<=seqLength;j++){
//				if(!siteFlag[j]) {
				if(node[0].seq_aa[j]!='-'){
					if(node[i].seq_aa[j]=='.')
						fprintf(output,"%c",'-');
					else fprintf(output,"%c",node[i].seq_aa[j]);
				}
			}
			fprintf(output,"\n");
		}
	}
*/
	fclose(output);

}


void chooseSeq1(char *listFile, int number)
{
	int i,j,k;
	int listIndex=1;
	char **listName;
	int *accordIndex;
	char outFile[20],outFile1[20],outFile2[20];
	listName=cmatrix(1,seqNum,1,25);
	FILE *output,*output1;
	for(i=0;i<20;i++){
		if(listFile[i]!='.') {
			outFile[i]=listFile[i];
		}
		else break;
	}
	strcpy(outFile1,outFile);
	strcpy(outFile2,outFile);

	ifstream input;
	input.open(listFile,ios::in);
	for(i=1;i<=seqNum;i++){
		listIndex=i;
		input>>listName[i];
		if(strcmp(listName[i],"")==0) break;
	}
	input.close();
	cout<<"finishing reading list name"<<endl;
	accordIndex=ivector(1,listIndex-1);
	for(i=1;i<listIndex;i++){
		bool tmpFlag=false;
		for(j=1;j<=seqNum;j++){
			if(strcmp(listName[i],node[j].nodeName)==0){
				accordIndex[i]=j;
				tmpFlag=true;
				break;
			}
		}
		if(!tmpFlag) { printf("sequence: %s can not find according sequence\n",listName[i]); exit(0); }
	}

	//calculate the identity matrix
	double **idMatrix;
	idMatrix=dmatrix(1,listIndex-1,1,listIndex-1);
	for(i=1;i<listIndex-1;i++){
		for(j=i+1;j<=listIndex-1;j++){
			int count=0; int length=0;
			for(k=1;k<=seqLength;k++){
				if(((node[accordIndex[i]].seq_aa[k]!='-')&&(node[accordIndex[i]].seq_aa[k]!='.'))&&((node[accordIndex[j]].seq_aa[k]!='-')&&(node[accordIndex[j]].seq_aa[k]!='.')))
				{
					length++;
					if(node[accordIndex[i]].seq_aa[k]==node[accordIndex[j]].seq_aa[k]) count++;
				}
			}
			if(length>0) idMatrix[i][j]=1.0*count/length;
			else idMatrix[i][j]=0.0;
		}
	}
	for(i=1;i<listIndex-1;i++) for(j=i+1;j<=listIndex-1;j++) idMatrix[j][i]=idMatrix[i][j];
	for(i=1;i<=listIndex-1;i++) idMatrix[i][i]=0.0;
	//test
/*	for(i=1;i<=listIndex-1;i++){
		printf("%20s ",listName[i]);
		for(j=1;j<=listIndex-1;j++) printf("%3.0f",idMatrix[i][j]*100.0);
		printf("\n");
	}
*/
	double *idvector;
	idvector=dvector(1,listIndex-1);
	int *idIndex;
	idIndex=ivector(1,listIndex-1);
	for(i=1;i<=listIndex-1;i++){
		idvector[i]=0.0; idIndex[i]=i;
		for(j=1;j<=listIndex-1;j++){
			if(j!=i) idvector[i]+=idMatrix[i][j];
		}
		idvector[i]/=(listIndex-2);
//		printf("%-20s  %d\n", listName[i],(int) (idvector[i]*100.0));
	}

	//sort idvector
	for(i=1;i<listIndex-1;i++){
		for(j=i+1;j<=listIndex-1;j++){
			if(idvector[i]<idvector[j]){
				double tmpvector=idvector[i]; int tmpIndex=idIndex[i];
				idvector[i]=idvector[j]; idIndex[i]=idIndex[j];
				idvector[j]=tmpvector; idIndex[j]=tmpIndex;
			}
		}
	}
	if(number>listIndex-1) {
		printf("WARNING: sample base is smaller than your expectation\n");
		number=listIndex-1;
	}
	if(idvector[number]<0.26) {
		while(idvector[number]<0.26){
			number--;
		}
	}
	bool *siteFlag;
	siteFlag=new bool [seqLength+1];
	for(i=1;i<=seqLength;i++){
		int count=0; siteFlag[i]=false;
		for(j=1;j<=number;j++){
			if((node[accordIndex[idIndex[j]]].seq_aa[i]=='-')||(node[accordIndex[idIndex[j]]].seq_aa[i]=='.')) count++;
		}
		if(2*count>number) siteFlag[i]=true;
	}

	//output
	strcat(outFile,"_1.aln");
	output=fopen(outFile,"w");
	for(i=1;i<=number;i++){
		fprintf(output,"%20s ",node[accordIndex[idIndex[i]]].nodeName);
		for(j=1;j<=seqLength;j++){
			if(!siteFlag[j]) {
				if(node[accordIndex[idIndex[i]]].seq_aa[j]=='.')
					fprintf(output,"%c",'-');
				else fprintf(output,"%c",node[accordIndex[idIndex[i]]].seq_aa[j]);
			}
		}
		fprintf(output,"\n");
	}
	fclose(output);
	strcat(outFile1,"_2.aln");
	output1=fopen(outFile1,"w");
	for(i=1;i<=number;i++){
		fprintf(output1,"%20s ",node[accordIndex[idIndex[i]]].nodeName);
		for(j=1;j<=seqLength;j++)
			fprintf(output1,"%c",node[accordIndex[idIndex[i]]].seq_aa[j]);
		fprintf(output1,"\n");
	}
	fclose(output1);

	//statistics analysis for mean, var, max, min, histograph
	int *binCount;
	double *bin;
	bin=dvector(1,11);
	binCount=ivector(1,10);
	for(i=1;i<=10;i++) binCount[i]=0;
	double mean=0.0, var=0.0, max=-0.1, min=101.0;
	for(i=1;i<number-1;i++){
		for(j=i+1;j<=number;j++){
			if(idMatrix[idIndex[i]][idIndex[j]]>max) max=idMatrix[idIndex[i]][idIndex[j]];
			if(idMatrix[idIndex[i]][idIndex[j]]<min) min=idMatrix[idIndex[i]][idIndex[j]];
			mean+=idMatrix[idIndex[i]][idIndex[j]];
		}
	}
	mean=mean/(number*(number-1)/2);
	bin[1]=min; bin[11]=max;
	for(i=2;i<=10;i++) bin[i]=bin[i-1]+(max-min)/10.0;
	for(i=1;i<number-1;i++){
		for(j=i+1;j<=number;j++){
			var+=(idMatrix[idIndex[i]][idIndex[j]]-mean)*(idMatrix[idIndex[i]][idIndex[j]]-mean);
			for(k=1;k<=10;k++){
				if((idMatrix[idIndex[i]][idIndex[j]]>bin[k])&&(idMatrix[idIndex[i]][idIndex[j]]<bin[k+1])) {
					binCount[k]++;
					break;
				}
			}
		}
	}
	var=sqrt(var/(number*(number-1)/2-1));
	//find the biggest binCount then normalize
	int binMax=binCount[1];
	for(i=2;i<=10;i++) if(binCount[i]>binMax) binMax=binCount[i];
	for(i=1;i<=10;i++) binCount[i]=binCount[i]*100/binMax;
	//output
	strcat(outFile2,".his");
	output=fopen(outFile2,"w");
	fprintf(output,"%s\n",outFile2);
	fprintf(output,"Mean: %3.0f\nVariance: %3.0f\nMax Value: %3.0f\nMin Value: %3.0f\n",mean*100,var*100,max*100,min*100);
	fprintf(output,"Histograph:\n");
	for(i=1;i<=10;i++){
		fprintf(output,"bin%2d(%2.0f) ",i,(100*(bin[i+1]+bin[i])/2));
		for(j=1;j<=binCount[i];j++) fprintf(output,"=");
		fprintf(output,"\n");
	}
	fprintf(output,"\n\n");
	for(i=1;i<=number;i++)	fprintf(output,"%-20s  %d\n", listName[idIndex[i]],(int) (idvector[i]*100.0));

	fclose(output);

	free_dvector(bin,1,11);
	free_ivector(binCount,1,10);

}


void chooseSeq(char *listFile, int number)
{
	double ran1(long *idum);

	int i,j;
	char **listName;
	char **chosenName;
	int *chosenIndex;
	int listIndex=1;
	bool *siteFlag;
	char outFile[20],outFile1[20];
	siteFlag=new bool [seqLength+1];
	listName=cmatrix(1,seqNum,1,25);
	chosenName=cmatrix(1,number,1,25);
	chosenIndex=ivector(1,number);
	ifstream input;
	FILE *output,*output1;
	input.open(listFile,ios::in);
	for(i=1;i<=seqNum;i++){
		listIndex=i;
		input>>listName[i];
		if(strcmp(listName[i],"")==0) break;
	}
	input.close();
	cout<<"finishing reading list name"<<endl;
	cout<<"the number you choose:"<<number<<endl;
	int ranIndex=1;
	if(number<listIndex-1){
		while(ranIndex<=number){
			int tmpIndex=(int) (ran1(&idum)*(listIndex-1)+1.0);
	//		cout<<"tmpIndex:"<<tmpIndex<<"   ranIndex:"<<ranIndex<<endl;
	//		getchar();
			bool tmpFlag=false;
			for(i=1;i<ranIndex;i++) {
				if(strcmp(listName[tmpIndex],chosenName[i])==0) { tmpFlag=true; break;}
			}
			if(!tmpFlag){
				strcpy(chosenName[ranIndex],listName[tmpIndex]);
				ranIndex++;
			}
		}
	}
	else{
		number=listIndex-1;
		for(i=1;i<=number;i++) strcpy(chosenName[i],listName[i]);
	}
	cout<<"actual number:"<<number<<endl;
	cout<<"finish choosing sequence\n";
	for(i=1;i<=number;i++){
		bool testFlag=false;
		for(j=1;j<=seqNum;j++){
			if(strcmp(chosenName[i],node[j].nodeName)==0){
				chosenIndex[i]=j;
				testFlag=true;
				break;
			}
		}
		if(!testFlag) printf("%s chosenName[i] did not find its sequence\n");
	}
	for(i=1;i<=seqLength;i++){
		int count=0;
		for(j=1;j<=number;j++){
			if((node[chosenIndex[j]].seq_aa[i]=='-')||(node[chosenIndex[j]].seq_aa[i]=='.')) count++;
		}
		if(2*count>=number) siteFlag[i]=true;
	}
	for(i=0;i<20;i++){
		if(listFile[i]!='.') {
			outFile[i]=listFile[i];
		}
		else break;
	}
	strcpy(outFile1,outFile);
	strcat(outFile,"_1.aln");
//	cout<<"outFile:"<<outFile<<endl;
	output=fopen(outFile,"w");
	for(i=1;i<=number;i++){
		fprintf(output,"%20s ",node[chosenIndex[i]].nodeName);
		for(j=1;j<=seqLength;j++){
			if(!siteFlag[j]) {
				if(node[chosenIndex[i]].seq_aa[j]=='.')
					fprintf(output,"%c",'-');
				else fprintf(output,"%c",node[chosenIndex[i]].seq_aa[j]);
			}
		}
		fprintf(output,"\n");
	}
	fclose(output);
//	strcpy(outFile," ");
	strcat(outFile1,"_2.aln");
//	cout<<"outFile1:"<<outFile1<<endl;
	output1=fopen(outFile1,"w");
	for(i=1;i<=number;i++){
		fprintf(output1,"%20s ",node[chosenIndex[i]].nodeName);
		for(j=1;j<=seqLength;j++)
			fprintf(output1,"%c",node[chosenIndex[i]].seq_aa[j]);
		fprintf(output1,"\n");
	}
	fclose(output1);

//	cout<<"listIndex:"<<listIndex<<endl;
//	for(i=1;i<listIndex;i++) cout<<listName[i]<<endl;
//	cout<<"ranIndex:"<<ranIndex<<endl;
//	for(i=1;i<ranIndex;i++) cout<<chosenName[i]<<endl;
}


static char str[MAXSTR+1];

char **aname, **aseq;
int nal, alilen, *astart, *alen;


static char *mymalloc(int size)
//int size;
{
	char *buf;

	if ((buf = (char *)malloc((size_t) (size*sizeof(char)))) == NULL) {
		fprintf(stderr, "Not enough memory: %d\n", size);
		exit(1);
	}
	return buf;
}

char *ccvector(long nl, long nh)
/* allocate an char vector with subscript range v[nl..nh] */
{
	char *v;

	v=(char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(char)));
	if (!v) nrerror("allocation failure in cvector()");
	return v-nl+NR_END;
}

char *strsave(char *str)
//char *str;
{
	char *buf;
	int l;

	l = strlen(str);
	buf = mymalloc(l + 1);
//	buf=ccvector(0,l+1);
	strcpy(buf, str);
	return buf;
}

char *strnsave(char *str, int l)
//char *str;
//int l;
{
	char *buf;

	buf = mymalloc(l + 1);
//	buf=ccvector(0,l+1);
	memcpy(buf, str, l);
	buf[l] = '\0';
	return buf;
}

static char **incbuf(int n, char **was)
{
	char **buf;
	buf = (char **) malloc((size_t) (n+1) * sizeof(buf[0]));
	if(!buf) {fprintf(stderr, "Not enough memory: %d\n",n+1); exit(1); }
//	buf=cmatrix(0,(n+1),0,(n+1));
	if (n > 0) {
		memcpy(buf, was, n * sizeof(was[0]));
		free(was);
	}
	buf[n] = NULL;
	return buf;
}

static int *incibuf(int n, int *was)
//int n, *was;
{
	int *ibuf;

	ibuf = (int *) malloc((size_t) (n+1) * sizeof(int));
	if(!ibuf) {fprintf(stderr, "Not enough memory: %d\n",n+1); exit(1); }
//	ibuf=ivector(0,(n+1)*sizeof(ibuf[0]));
	if (n > 0) {
		memcpy(ibuf, was, n * sizeof(was[0]));
		free(was);
	}
	ibuf[n] = 0;
	return ibuf;
}


void readali(char *filename)
/* fills alignment from the file *filename to
nal - number of sequences in the alignment;
alilen - alignment length;
aname[0..nal-1][0..until '\0'] - names of sequences;
aseq[0..nal-1][0..alilen-1] - sequences;
astart[0..nal-1] - starting residue number;
alen[0..nal-1] - sequence length;
*/

{
	FILE *fp;
	char *s, *ss, *seqbuf;
	int i,j,k, n, l, len, len0;
	int ii,mark=1;

	printf("Reading sequences\n");
	if ((fp = fopen(filename, "r")) == NULL) {
		fprintf(stderr, "No such file: \"%s\"\n", filename);
//		err_readali(1);
		exit(1);
	}
	alilen = 0;
	nal = 0;
	n = 0;
	if(fgets(str, MAXSTR, fp) != NULL) {
		if(strncmp(str,"CLUSTAL W",9)!=0){rewind(fp);}
					   }
	while (fgets(str, MAXSTR, fp) != NULL) {
		for (ss = str; isspace(*ss); ss++) ;
//		printf("SS1:%s\n",ss);
		if ((ii<=ss-str)&&(mark==0)) {continue;}
		if (*ss == '\0') {
			if (n == 0) {
				continue;
			}
			if (nal == 0) {
				if (n == 0) {
					fprintf(stderr, "No alignments read\n");
//					err_readali(2);
					exit(1);
				}
				nal = n;
			} else if (n != nal) {
				fprintf(stderr, "Wrong nal, was: %d, now: %d\n", nal, n);
//				err_readali(3);
				exit(1);
			}
			n = 0;
			continue;
		}
		for (s = ss; *s != '\0' && !isspace(*s); s++) ;
//		printf("S1:%s\n",s);
		*s++ = '\0';
		if (nal == 0) {
			astart = incibuf(n, astart);
			alen = incibuf(n, alen);
			aseq = incbuf(n, aseq);
			aname = incbuf(n, aname);
			aname[n] = strsave(ss);
		} else {
			if (n < 0 || n >= nal) {
				fprintf(stderr, "Bad sequence number: %d of %d\n", n, nal);
//				err_readali(4);
				exit(1);
			}
			if (strcmp(ss, aname[n]) != 0) {
				fprintf(stderr, "Names do not match");
				fprintf(stderr, ", was: %s, now: %s\n", aname[n], ss);
//				err_readali(5);
				exit(1);
			}
		}
		for (ss = s; isspace(*ss); ss++);
//		printf("SS2:%s\n",ss);
		if(mark==1){
		ii = ss-str;
		mark=0;}

		for (s = ss; isdigit(*s); s++) ;
//		printf("S2:%s\n",s);
		if (isspace(*s)) {
			if (nal == 0) {
				astart[n] = atoi(ss);
			}
			for (ss = s; isspace(*ss); ss++);
		}
		for (s = ss, l = 0; *s != '\0' && !isspace(*s); s++) {
			if (isalpha(*s)) {
				l++;
			}
		}
		len = s - ss;
//		printf("n:%d  nal:%d\n",n,nal);
		if (n == 0) {
			len0 = len;
			alilen += len;
		} else if (len != len0) {
			fprintf(stderr, "wrong len for %s", aname[n]);
			fprintf(stderr, ", was: %d, now: %d\n", len0, len);
//			err_readali(6);
			exit(1);
		}
		alen[n] += l;
		if (aseq[n] == NULL) {
			aseq[n] = strnsave(ss, len);
		} else {
//			printf("save sequence\n");
//			seqbuf=ccvector(0,alilen);
			seqbuf=mymalloc(alilen+1);
			memcpy(seqbuf, aseq[n], alilen-len);
			free(aseq[n]);
			aseq[n] = seqbuf;
			memcpy(seqbuf+alilen-len, ss, len);
			seqbuf[alilen] = '\0';
//			printf("after save sequence\n");
		}
		n++;
	}
	if (nal == 0) {
		if (n == 0) {
			fprintf(stderr, "No alignments read\n");
//			err_readali(7);
			exit(1);
		}
		nal = n;
	} else if (n != 0 && n != nal) {
		fprintf(stderr, "Wrong nal, was: %d, now: %d\n", nal, n);
//		err_readali(8);
		exit(1);
	}
	fclose(fp);
//	nal=n;

//	printf("n=%d alilen=%d nal=%d\n",n,alilen,nal);
	//define new structures
	seqNum=nal; seqLength=alilen;
//	int DIM=(seqNum-1)*2+1;
//	node=new NODE [DIM];
	int DIM=seqNum+1;
	node=new NODE [seqNum+1];
	for(i=0;i<DIM;i++) {
		node[i].seq_aa=new char [seqLength+2];
		node[i].seq_num=new int [seqLength+1];
		node[i].prob= dmatrix(1,seqLength+1,1,NP+1);
	}
	//transfer sequence to NODE structure
	for(i=1;i<=(seqNum*2-2);i++) node[i].nodeNum=i;
	for(i=1;i<=seqNum;i++) strcpy(node[i].nodeName,aname[i-1]);
	for(i=0;i<=seqNum/*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];
	for(i=1;i<=seqNum;i++) for(j=1;j<=seqLength;j++) node[i].seq_aa[j]=toupper(node[i].seq_aa[j]);
	//deal with non-amino acid character
	for(i=1;i<=seqNum;i++){
		for(j=1;j<=seqLength;j++){
			bool flag=false;
			for(int k=0;k<=NP;k++) if(node[i].seq_aa[j]==aa[k]) flag=true;
			if(!flag) 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);

	//put template sequence into node[0]
	for(i=0;i<=seqNum-1;i++) node[i]=node[i+1];
	seqNum--;
//	printf("seqNum:%d\n",seqNum);
	//find the template sequence non gap length
	int stdLength=0;
	for(i=1;i<=seqLength;i++) if(node[0].seq_aa[i]!='-') stdLength++;
	//find those sequences that can not cover 80% of template sequence non-gap position
	int *lengthDist; lengthDist=ivector(1,seqNum);
	for(i=1;i<=seqNum;i++){
		lengthDist[i]=0;
		for(j=1;j<=seqLength;j++){
			if(node[0].seq_aa[j]!='-'){
				if(node[i].seq_aa[j]!='-') lengthDist[i]++;
			}
		}
		lengthDist[i]=lengthDist[i]*100/stdLength;
	}
	//histogram
	int *binCount;
	binCount=ivector(1,20);
	for(i=1;i<=20;i++) binCount[i]=0;
	for(i=1;i<=seqNum;i++){
			for(k=0;k<=19;k++){
				if((lengthDist[i]>=k*5)&&(lengthDist[i]<(k+1)*5)) {
					binCount[k+1]++;
					break;
				}
				else if(lengthDist[i]>95) {
					binCount[20]++;
					break;
				}
			}

	}
//	for(i=1;i<=20;i++) printf("%5d\n",binCount[i]);
	//delete sequence
	for(i=1;i<=seqNum;i++){
		if(lengthDist[i]<75){
/*			for(j=i;j<seqNum;j++) { node[j]=node[j+1]; lengthDist[j]=lengthDist[j+1]; }
			i--;
			seqNum--;
*/			while(lengthDist[seqNum]<75) { seqNum--;}
			node[i]=node[seqNum];
			seqNum--;

		}
	}
//	printf("seqNum:%d\n",seqNum);
	//print the sequence info
	fprintf(resultFile,"The input sequences:\n");
	for(i=0;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");
	}
	fprintf(resultFile,"\n\n");

}


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