/*
 Program From Progressive Alignment of Amino Acid
 Sequences and Construction of Phylogenetic Trees
 from Them, by D. F. Feng and R. F. Doolittle,
 Methods of Enzymology, v266, p368 (1996)
				2/7/1996
 Read "readme.doc" before using.
*/

#include <stdio.h>
#include <math.h>
#define NSEQ  52         	/*Maximum number of sequences*/
#define MAXSE 501        	/*Maximum length for sequences*/ 
#define MAXAS MAXSE+100  	/*Maximum length for aligned sequences*/
#define tnsq 2*NSEQ-3		/*Total number of branches*/
#define nsq NSEQ*(NSEQ-1)/2	/*Upper triangle*/
#define tpnsq 2*(2*NSEQ-3)	/*Total number of matrix elements*/

int pam250[21][21] =	{
{20, 8, 6, 5, 6, 5, 4, 3, 3, 3, 5, 4, 3, 3, 6, 2, 6, 4, 8, 0, 0},
{ 8,10, 9, 9, 9, 9, 9, 8, 8, 7, 7, 8, 8, 6, 7, 5, 7, 5, 5, 6, 0},
{ 6, 9,11, 8, 9, 8, 8, 8, 8, 7, 7, 7, 8, 7, 8, 6, 8, 5, 5, 3, 0},
{ 5, 9, 8,14, 9, 7, 7, 7, 7, 8, 8, 8, 7, 6, 6, 5, 7, 3, 3, 2, 0},
{ 6, 9, 9, 9,10, 9, 8, 8, 8, 8, 7, 6, 7, 7, 7, 6, 8, 4, 5, 2, 0},
{ 5, 9, 8, 7, 9,13, 8, 9, 8, 7, 6, 5, 6, 5, 5, 4, 7, 3, 3, 1, 0},
{ 4, 9, 8, 7, 8, 8,10,10, 9, 9,10, 8, 9, 6, 6, 5, 6, 4, 6, 4, 0},
{ 3, 8, 8, 7, 8, 9,10,12,11,10, 9, 7, 8, 5, 6, 4, 6, 2, 4, 1, 0},
{ 3, 8, 8, 7, 8, 8, 9,11,12,10, 9, 7, 8, 6, 6, 5, 6, 3, 4, 1, 0},
{ 3, 7, 7, 8, 8, 7, 9,10,10,12,11, 9, 9, 7, 6, 6, 6, 3, 4, 3, 0},
{ 5, 7, 7, 8, 7, 6,10, 9, 9,11,14,10, 8, 6, 6, 6, 6, 6, 8, 5, 0},
{ 4, 8, 7, 8, 6, 5, 8, 7, 7, 9,10,14,11, 8, 6, 5, 6, 4, 4,10, 0},
{ 3, 8, 8, 7, 7, 6, 9, 8, 8, 9, 8,11,13, 8, 6, 5, 6, 3, 4, 5, 0},
{ 3, 6, 7, 6, 7, 5, 6, 5, 6, 7, 6, 8, 8,14,10,12,10, 8, 6, 4, 0},
{ 6, 7, 8, 6, 7, 5, 6, 6, 6, 6, 6, 6, 6,10,13,10,12, 9, 7, 3, 0},
{ 2, 5, 6, 5, 6, 4, 5, 4, 5, 6, 6, 5, 5,12,10,14,10,10, 7, 6, 0},
{ 6, 7, 8, 7, 8, 7, 6, 6, 6, 6, 6, 6, 6,10,12,10,12, 7, 6, 2, 0},
{ 4, 5, 5, 3, 4, 3, 4, 2, 3, 3, 6, 4, 3, 8, 9,10, 7,17,15, 8, 0},
{ 8, 5, 5, 3, 5, 3, 6, 4, 4, 4, 8, 4, 4, 6, 7, 7, 6,15,18, 8, 0},
{ 0, 6, 3, 2, 2, 1, 4, 1, 1, 3, 5,10, 5, 4, 3, 6, 2, 8, 8,25, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
			};

int Gonnet[21][21] =	{
{20, 8, 7, 5, 9, 6, 6, 5, 5, 6, 7, 6, 5, 7, 7, 6, 8, 7, 8, 7, 0},
{ 8,10,10, 8, 9, 8, 9, 9, 8, 8, 8, 8, 8, 7, 6, 6, 7, 5, 6, 5, 0},
{ 7,10,11, 8, 9, 7,10, 8, 8, 8, 8, 8, 8, 7, 7, 7, 8, 6, 6, 4, 0},
{ 5, 8, 8,15, 8, 6, 7, 7, 7, 8, 7, 7, 7, 6, 5, 6, 6, 4, 5, 3, 0},
{ 9, 9, 9, 8,10, 9, 8, 8, 8, 8, 7, 7, 8, 7, 7, 7, 8, 6, 6, 4, 0},
{ 6, 8, 7, 6, 9,15, 8, 8, 7, 7, 7, 7, 7, 4, 3, 4, 4, 3, 4, 4, 0},
{ 6, 9,10, 7, 8, 8,12,10, 9, 9, 9, 8, 9, 6, 5, 5, 6, 5, 7, 4, 0},
{ 5, 9, 8, 7, 8, 8,10,13,11, 9, 8, 8, 9, 5, 4, 4, 5, 3, 5, 3, 0},
{ 5, 8, 8, 7, 8, 7, 9,11,12,10, 8, 8, 9, 6, 5, 5, 6, 4, 5, 4, 0},
{ 6, 8, 8, 8, 8, 7, 9, 9,10,11, 9,10,10, 7, 6, 6, 6, 5, 6, 5, 0},
{ 7, 8, 8, 7, 7, 7, 9, 8, 8, 9,14, 9, 9, 7, 6, 6, 6, 8,10, 7, 0},
{ 6, 8, 8, 7, 7, 7, 8, 8, 8,10, 9,13,11, 6, 6, 6, 6, 5, 6, 6, 0},
{ 5, 8, 8, 7, 8, 7, 9, 9, 9,10, 9,11,11, 7, 6, 6, 6, 5, 6, 4, 0},
{ 7, 7, 7, 6, 7, 4, 6, 5, 6, 7, 7, 6, 7,12,11,11,10,10, 8, 7, 0},
{ 7, 6, 7, 5, 7, 3, 5, 4, 5, 6, 6, 6, 6,11,12,11,11, 9, 7, 6, 0},
{ 6, 6, 7, 6, 7, 4, 5, 4, 5, 6, 6, 6, 6,11,11,12,10,10, 8, 7, 0},
{ 8, 7, 8, 6, 8, 4, 6, 5, 6, 6, 6, 6, 6,10,11,10,11, 8, 7, 5, 0},
{ 7, 5, 6, 4, 6, 3, 5, 3, 4, 5, 8, 5, 5,10, 9,10, 8,15,13,12, 0},
{ 7, 6, 6, 5, 6, 4, 7, 5, 5, 6,10, 6, 6, 8, 7, 8, 7,13,15,12, 0},
{ 7, 5, 4, 3, 4, 4, 4, 3, 4, 5, 7, 6, 4, 7, 6, 7, 5,12,12,22, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
			};

int Blosum62[21][21] =	{
{13, 3, 3, 1, 4, 1, 1, 1, 0, 1, 1, 1, 1, 3, 3, 3, 3, 2, 2, 2, 0},
{ 3, 8, 5, 3, 5, 4, 5, 4, 4, 4, 3, 3, 4, 3, 2, 2, 2, 2, 2, 1, 0},
{ 3, 5, 9, 3, 4, 2, 4, 3, 3, 3, 2, 3, 3, 3, 3, 3, 4, 2, 2, 2, 0},
{ 1, 3, 3,11, 3, 2, 2, 3, 3, 3, 2, 2, 3, 2, 1, 1, 2, 0, 1, 0, 0},
{ 4, 5, 4, 3, 8, 4, 2, 2, 3, 3, 2, 3, 3, 3, 3, 3, 4, 2, 2, 1, 0},
{ 1, 4, 2, 2, 4,10, 4, 3, 2, 2, 2, 2, 2, 1, 0, 0, 1, 1, 1, 2, 0},
{ 1, 5, 4, 2, 2, 4,10, 5, 4, 4, 5, 4, 4, 2, 1, 1, 1, 1, 2, 0, 0},
{ 1, 4, 3, 3, 2, 3, 5,10, 6, 4, 3, 2, 3, 1, 1, 0, 1, 1, 1, 0, 0},
{ 0, 4, 3, 3, 3, 2, 4, 6, 9, 6, 4, 4, 5, 2, 1, 1, 2, 1, 2, 1, 0},
{ 1, 4, 3, 3, 3, 2, 4, 4, 6, 9, 4, 5, 5, 4, 1, 2, 2, 1, 3, 2, 0},
{ 1, 3, 2, 2, 2, 2, 5, 3, 4, 4,12, 4, 3, 2, 1, 1, 1, 3, 6, 2, 0},
{ 1, 3, 3, 2, 3, 2, 4, 2, 4, 5, 4, 9, 6, 3, 1, 2, 1, 1, 2, 1, 0},
{ 1, 4, 3, 3, 3, 2, 4, 3, 5, 5, 3, 6, 9, 3, 1, 2, 2, 1, 2, 1, 0},
{ 3, 3, 3, 2, 3, 1, 2, 1, 2, 4, 2, 3, 3, 9, 5, 6, 5, 4, 3, 3, 0},
{ 3, 2, 3, 1, 3, 0, 1, 1, 1, 1, 1, 1, 1, 5, 8, 6, 7, 4, 3, 1, 0},
{ 3, 2, 3, 1, 3, 0, 1, 0, 1, 2, 1, 2, 2, 6, 6, 8, 5, 4, 3, 2, 0},
{ 3, 2, 4, 2, 4, 1, 1, 1, 2, 2, 1, 1, 2, 5, 7, 5, 8, 3, 3, 1, 0},
{ 2, 2, 2, 0, 2, 1, 1, 1, 1, 1, 3, 1, 1, 4, 4, 4, 3,10, 7, 5, 0},
{ 2, 2, 2, 1, 2, 1, 2, 1, 2, 3, 6, 2, 2, 3, 3, 3, 3, 7,11, 6, 0},
{ 2, 1, 2, 0, 1, 2, 0, 0, 1, 2, 2, 1, 1, 3, 1, 2, 1, 5, 6,15, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
			};

int mptable[21][21];

int index[26]={  4, 7, 0, 7, 8,17, 5,10,
		14,20,12,15,13, 6,20, 3,
		 9,11, 1, 2,20,16,19,20,
		18, 8,
		 };

char seq1[MAXAS], seq2[MAXAS];  /*Used in shuffle, dummy arrays*/
char seqt[NSEQ][100],code[NSEQ][8];
char seq[NSEQ][MAXAS],oseq[NSEQ][MAXAS],wseq[NSEQ][MAXAS];
char aseq1[MAXAS],aseq2[MAXAS];
char clus[NSEQ];

int L[MAXAS][MAXAS], R[MAXAS][MAXAS], M[MAXAS][MAXAS];
int lseq[NSEQ], nseq, cnum[NSEQ];
int pen;     /*None negative gap penalty*/
int wgseq, wgt;
int iseq, family[NSEQ][NSEQ],ZGROUP,print_flag;
int loop,nrow,ncol;

float rndmax = {2147483647.};

double sreal[NSEQ][NSEQ],shmean[NSEQ][NSEQ],siden[NSEQ][NSEQ];
double D[NSEQ][NSEQ],pro[NSEQ][NSEQ];
double RB[tnsq],A[nsq][tnsq],Y[nsq];

main(ac,av)
int ac;
char *av[];
{
	FILE *fd1, *fd2, *fd3, *fd4, *fd5, *fd6, *fd7, *fd8, *fd9, *fopen();

	char dseq[MAXAS], dummy_array[100];  /*dummy sequence array*/
	char cdum;
	char segment[NSEQ];
	char ori_branch[400],rnd_branch[400];
	char orderstring[400],temp_string[400];

	int i,j,si,sj,NSDL,nax,cu_nseq,tmp_nseq,tmp_si,tmp_sj;
	int cflag1,cflag2,ik;
	int nsub_cl,kk,clus_cnt,dum;
	int pchange,nchange,ntimes,hits,ixx,achange,ichange;
	int ipp,kind,same,pc_low,pc_high;
	int k,t,v,q,n,ta[NSEQ],side,rn,flseq,ifront,iend;
	int count_x,cutoff_x,drop_x[200],ndrop_x;
	int adj;

	float bR,psd,psdmin,sum_psd,num_noneg;
	double sim_score1(),score1,score2,sum,ddum,cum;

   /* Check for correct command line */
	if(ac<4)
		{
		fprintf(stderr,"\nPurpose: Sample Multiple Alignments and");
		fprintf(stderr," Generate New Trees for Boostrapping\n\t");
		fprintf(stderr," Procedure.\n");
		fprintf(stderr,"Usage: ");
		fprintf(stderr,"%s rf wf1 M [>wf2]\n",av[0]);
		fprintf(stderr,"\trf:  All sequences arranged in ");
		fprintf(stderr,"binary alignment order.\n");
		fprintf(stderr,"\twf1: Aligned sequences.\n");
		fprintf(stderr,"\tM:   Scoring matrix.\n");
		fprintf(stderr,"\t\tType D(for PAM250), ");
		fprintf(stderr,"G(for Gonnet), or B(for Blosum62).\n");
		fprintf(stderr,"\twf2: Alignment results.\n");
		fprintf(stderr,"Current limit: Sequence length =< %d,",MAXSE-1);
		fprintf(stderr," NSEQ =< %d.\n\n",NSEQ);
		exit(0);
		}

   /* Open files */
	if((fd1 = fopen(av[1],"r")) == NULL)   /*error check*/
		exit(fprintf(stderr,"%s: can't open %s\n",av[1]) );
	if((fd2 = fopen(av[2],"w")) == NULL)   /*error check*/
		exit(fprintf(stderr,"%s: can't open %s\n",av[2]));
	if((fd9 = fopen("rboot","w")) == NULL)   /*error check*/
		exit(fprintf(stderr,"%s: can't open rboot\n"));

	srandom(6);

   /* Set weighting, shuffle times and gap penalty */
	wgseq=50;
	wgt=1;
	if(*av[3]=='D')
		{
		pen=8;
		for(i=0; i<21; i++)
			for(j=0; j<21; ++j)
				mptable[i][j] = pam250[i][j];
		}
	else if(*av[3]=='G')
		{
		pen=8;
		for(i=0; i<21; i++)
			for(j=0; j<21; ++j)
				mptable[i][j] = Gonnet[i][j];
		}
	else if(*av[3]=='B')
		{
		pen=6;
		for(i=0; i<21; i++)
			for(j=0; j<21; ++j)
				mptable[i][j] = Blosum62[i][j];
		}
	else
		{
		printf("\n*** Matrix-%s is not available. ***\n",av[3]);
		exit(0);
		}

	print_flag = 0; /*Control printing D matrix */
	psdmin = 10000.0; /*Minimum %s.d.*/
	for(i=0; i<200; ++i)
		drop_x[i] = -1;

   /* Read sequences */
	Read_seq(fd1);
	fclose(fd1);

	printf("\t\t-- PROGRESSIVE (similarity) ALIGNMENT --\n");
	if(*av[3]=='D')
		printf("\t\t\t\tPAM250\n\n");
	else if(*av[3]=='B')
		printf("\t\t\t\tBlosum62\n\n");
	else if(*av[3]=='G')
		printf("\t\t\t\tGonnet\n\n");
	printf("Sequences are entered in the following order:\n");
	for(i=0; i<nseq; ++i)
		printf("%2d)%s%s",i+1,code[i],seqt[i]);
	printf("\n");

   /* Set output parameters */
	NSDL=nseq;
	nax=35;
	
   /* Align sequences recursively*/
	/*Aligning the first two or three sequences*/
	si = 0;
	sj = 1;
	cu_nseq = 2;
	nsub_cl = 0;
	cflag1=cflag2=0;
	Align(si,sj,nsub_cl,cflag1,cflag2);
	BkUtoX(0);
	if(nseq == 2)
		goto skip; /*Print results if only two sequences*/

	   /*A & B are aligned, align C and realign A & B*/
	control(&si,&sj,&cu_nseq,3,nsub_cl,cflag1,cflag2);
	score1 = sim_score1(3);  /*score for (A,B)C arrangement*/

	for(i=0; i<3; ++i)
		{
		strcpy(wseq[i],seq[i]);
		strcpy(seq[i],oseq[i]);
		}
	switch_array(2,seq);  /*Switching A & B*/
	cu_nseq = 2;
	Align(si,sj,nsub_cl,cflag1,cflag2);
	BkUtoX(0);
	control(&si,&sj,&cu_nseq,3,nsub_cl,cflag1,cflag2);

	score2 = sim_score1(3);   /*score for (B,A)C arrangement*/
	if(score1 >= score2)
		{
		for(i=0; i<2; ++i)
			strcpy(seq[i],wseq[i]);
		if(nseq == 3)
			{
			strcpy(seq[2],wseq[2]);
			rm_cX1(3,seq);
			goto skip;
			}
		}
	else
		{
		strcpy(dummy_array,code[0]);
		strcpy(code[0],code[1]);
		strcpy(code[1],dummy_array);
		strcpy(dummy_array,seqt[0]);
		strcpy(seqt[0],seqt[1]);
		strcpy(seqt[1],dummy_array);
		switch_array(2,oseq);
		if(nseq == 3)
			goto skip;
		}
	strcpy(seq[2],oseq[2]);

	tmp_nseq = 3;
	si=0;
	sj=1;
	while(tmp_nseq <= nseq)
		{
		nsub_cl = 0;
		if(clus[sj+1]=='J')
			cflag1 = 1;
		else
			cflag1= 0;
		if((sj+2)<nseq)
			{
			if(clus[sj+2]=='J')
				cflag2 = 1;
			else
				cflag2 = 0;
			}
		else cflag2 = 2;

		tmp_si = si;
		tmp_sj = sj;
		cu_nseq = tmp_nseq-1;

		if(cflag1==1 && cflag2==1)
			{
			kk=2;
			nsub_cl=0;
			clus_cnt=1;
			while(clus[sj+kk]=='J' && cnum[sj+kk]==clus_cnt)
				{
				nsub_cl++;
				kk++;
				clus_cnt++;
				}
			control(&tmp_si,&tmp_sj,&cu_nseq,tmp_nseq,nsub_cl,
				cflag1,cflag2);
			while(nsub_cl>0)
				{
				realign_cl(sj+1,sj+2);
				si++;
				sj++;
				tmp_nseq++;
				nsub_cl--;
				}
			goto NORM;
			}

		control(&tmp_si,&tmp_sj,&cu_nseq,tmp_nseq,nsub_cl,
			cflag1,cflag2);
		rm_cX1(tmp_nseq,seq);
		score1 = sim_score1(tmp_nseq);


		if(tmp_nseq==nseq) /*When end is reached, no exchange*/
			goto skip;

		for(i=0; i<tmp_nseq; ++i)
			strcpy(wseq[i],seq[i]);
		strcpy(seq[tmp_nseq-1],oseq[tmp_nseq]);
		strcpy(seq[tmp_nseq],oseq[tmp_nseq-1]);
		if(cflag1==0 && cflag2==1)
			{
			kk=3;
			clus_cnt=1;
			while(clus[sj+kk]=='J' && cnum[sj+kk]==clus_cnt)
				{
				nsub_cl++;
				kk++;
				clus_cnt++;
				}
			}
		cu_nseq = tmp_nseq-1;
		tmp_si = si;
		tmp_sj = sj;
		control(&tmp_si,&tmp_sj,&cu_nseq,tmp_nseq,nsub_cl,
			cflag1,cflag2);
		rm_cX1(tmp_nseq,seq);
		score2 = sim_score1(tmp_nseq);

		if(score1>=score2)
			{
			for(i=0; i<tmp_nseq; ++i)
				strcpy(seq[i],wseq[i]);
			strcpy(seq[tmp_nseq],oseq[tmp_nseq]);
			}
		else
			{
			switch_array(tmp_nseq+1,oseq);

			strcpy(dummy_array,code[tmp_nseq]);
			strcpy(code[tmp_nseq],code[tmp_nseq-1]);
			strcpy(code[tmp_nseq-1],dummy_array);

			strcpy(dummy_array,seqt[tmp_nseq]);
			strcpy(seqt[tmp_nseq],seqt[tmp_nseq-1]);
			strcpy(seqt[tmp_nseq-1],dummy_array);

			cdum=clus[tmp_nseq];
			clus[tmp_nseq]=clus[tmp_nseq-1];
			clus[tmp_nseq-1]=cdum;

			dum=cnum[tmp_nseq];
			cnum[tmp_nseq]=cnum[tmp_nseq-1];
			cnum[tmp_nseq-1]=dum;

			if(cflag1==0 && cflag2==1)
				{
				ik=sj+3;
				si++;
				sj++;
				tmp_nseq++;
				while(nsub_cl>0)
					{
					switch_array(ik+1,seq);
					switch_array(ik+1,oseq);

					strcpy(dummy_array,code[ik]);
					strcpy(code[ik],code[ik-1]);
					strcpy(code[ik-1],dummy_array);

					strcpy(dummy_array,seqt[ik]);
					strcpy(seqt[ik],seqt[ik-1]);
					strcpy(seqt[ik-1],dummy_array);

					cdum=clus[ik];
					clus[ik]=clus[ik-1];
					clus[ik-1]=cdum;

					dum=cnum[tmp_nseq];
					cnum[tmp_nseq]=cnum[tmp_nseq-1];
					cnum[tmp_nseq-1]=dum;
					
					realign_cl(ik-2,ik-1);
					si++;
					sj++;
					tmp_nseq++;
					nsub_cl--;
					ik++;
					}
				goto NORM1;
				}
			}
		NORM:
		si++;
		sj++;
		tmp_nseq++;
		NORM1:
			;
		}
	skip:
	rpJwX();
	rm_cX1(nseq,seq);

   /* Re-format the original sequences: remove the blank space at the 
      beginning and append a '*' at the end of the sequences */
	for(i=0; i<nseq; ++i)
		format_seq(i);

   /* Print the multiple alignment */
	mulalign(fd2,NSDL,nax);
	fclose(fd2);

   /* Remove excess overhangs at both ends */
	ifront=0;
	iend=10000;  /*If aligned seq >10000, this number has to be increased*/
	flseq = strlen(seq[0])-1;
	for(i=0; i<nseq; ++i)
		{
		for(j=0; j<flseq; )
			{
			if(seq[i][j]=='X')
				j++;
			else if(j>ifront)
				ifront = j;
			else
				break;
			}
		}
	for(i=0; i<nseq; ++i)
		{
		for(j=flseq-1; j>0; )
			{
			if(seq[i][j]=='X')
				j--;
			else if(j<iend)
				iend = j;
			else
				break;
			}
		}
	iend++;
	for(i=0; i<nseq; ++i)
		{
		k=0;
		for(j=ifront; j<iend; ++j)
			wseq[i][k++]=seq[i][j];
		wseq[i][k]='*';	
		wseq[i][k+1]='\0';	
		}
	
	flseq = strlen(wseq[0])-1;

   /* Setup exclusion numbers for bootstrap method */
	cutoff_x = nseq*0.6;
	k = 0;
	for(i=0; i<flseq; ++i)
		{
		count_x = 0;
		for(j=0; j<nseq; ++j)
			if(wseq[j][i]=='X')
				count_x++;
		if(count_x >= cutoff_x)
			drop_x[k++] = i;
		}
	ndrop_x = k;

	sum_psd = num_noneg = 0;
	for(loop=0; loop<18000; loop++)
	{
   /* Print computed parameters */
	para();

   /* Calculate Srand */
	shuffle();

   /* Compute difference score */ 
	adj = 0;
	for(i=0; i<nseq+1; ++i)
		D[i][i] = 0.0;
	for(i=0; i<nseq-1; ++i)
		{
		for(j=i+1; j<nseq; ++j)
			{
			if(sreal[i][j] <= shmean[i][j])
				{
				adj++;
			/*	printf("**%s,%s %8.3f %8.3f ",
				code[i],code[j],sreal[i][j],shmean[i][j]); */
				sreal[i][j]=shmean[i][j]+0.001;
			/*	printf("new_sreal=%8.3f\n",sreal[i][j]); */	
				}
			D[i+1][j+1] = D[j+1][i+1] =-(log((sreal[i][j]-
			shmean[i][j])/(siden[i][j]-shmean[i][j])))
			*100.0;
			}
		}
	if(loop<0)
		{
		printf("Siden:\n");
		printf("      ");
		for(i=0; i<nseq; ++i)
			printf("  %s  ",code[i]);
		for(i=0; i<nseq; ++i)
			{
			printf("\n%s ",code[i]);
			for(j=0; j<nseq; ++j)
				printf("%7.2f ",siden[i][j]);
			}
		printf("\n\n");
		printf("Srand:\n");
		printf("      ");
		for(i=0; i<nseq; ++i)
			printf("  %s  ",code[i]);
		for(i=0; i<nseq; ++i)
			{
			printf("\n%s ",code[i]);
			for(j=0; j<nseq; ++j)
				printf("%7.2f ",shmean[i][j]);
			}
		printf("\n\n");
		}
	if(loop==0)
		{
		printf("Distance:\n");
		printf("      ");
		for(i=0; i<nseq; ++i)
			printf("  %s  ",code[i]);
		for(i=0; i<nseq; ++i)
			{
			printf("\n%s ",code[i]);
			for(j=0; j<nseq; ++j)
				printf("%7.2f ",D[i+1][j+1]);
			}
		printf("\n\n");
		}

   /*Branching order & branch lengths*/
	bord();

	if((fd3 = fopen("albr1","w")) == NULL )   /*error check*/
		exit(fprintf(stderr,"%s: can't open albr1 in main.\n"));

	for(i=0;i<=ZGROUP;++i)
		{
		if(i==0)
			{
			j=0;
			for(v=0;family[i][v]!=0;++v)
				{
				if(v<26)
					ta[family[i][v]] = j+65;
				else if(v<52)
					ta[family[i][v]] = j+71;
				else
					exit(printf("Too many sequences!\n"));
				j++;
				}
			}
		for(n=0;family[i][n]!=0;++n)
			fprintf(fd3,"%c",ta[family[i][n]]);
		fprintf(fd3,"\n");
		}
	fprintf(fd3,"*\n");
	for(i=0;i<nseq-1;++i)
		{
		for(j=i+1;j<nseq;++j)
			fprintf(fd3," %7.2f",D[family[0][i]][family[0][j]]);
		fprintf(fd3,"\n");
		}
	fclose(fd3);
	blen();

	for(i=0; i<ncol; ++i)
		{
		if(RB[i]<0)
			goto jump;
		}
	sum=0;
	for(i=0; i<nrow; ++i)
		{
		if(Y[i]<=0.0)
			goto jump;
		cum=0;
		for(j=0; j<ncol; ++j)
			cum += A[i][j]*RB[j];
		ddum = fabs(cum - Y[i])*100/Y[i];
		sum += ddum*ddum;
		}
	psd=sqrt(sum/(nrow-1));
	sum_psd += psd;
	num_noneg++;
	if(num_noneg==101)
		goto out;
	if((fd8 = fopen("zecret","r")) == NULL)   /*error check*/
		exit(fprintf(stderr,"%s: can't open zecret for reading in main.\n"));
	fscanf(fd8,"%s",temp_string);
	fclose(fd8);
	i=j=0;
	fprintf(fd9,"(");
	while(temp_string[i]!='\0')
		{
		if(temp_string[i]>='A' && temp_string[i]<='Z')
			fprintf(fd9,"%s",code[family[0][j++]-1]);	
		else if(temp_string[i]>='a' && temp_string[i]<='z')
			fprintf(fd9,"%s",code[family[0][j++]-1]);	
		else
			fprintf(fd9,"%c",temp_string[i]);		
		i++;
		}
	fprintf(fd9,")\n");

	if((fd4 = fopen("albr","w")) == NULL)   /*error check*/
		exit(fprintf(stderr,"%s: can't open albr for writing in main.\n"));


	for(i=0;i<=ZGROUP;++i)
		{
		if(i==0)
			{
			j=0;
			for(v=0;family[i][v]!=0;++v)
				{
				if(v<26)
					ta[family[i][v]] = j+65;
				else if(v<52)
					ta[family[i][v]] = j+71;
				else
					exit(printf("Too many sequences!\n"));
				j++;
				}
			}
		for(n=0;family[i][n]!=0;++n)
			fprintf(fd4,"%c",ta[family[i][n]]);
		fprintf(fd4,"\n");
		}
	fprintf(fd4,"*\n");

	if((fd5 = fopen("algr","w")) == NULL)   /*error check*/
		exit(fprintf(stderr,"%s: can't open algr for writing in main.\n"));
	for(i=0; i<nseq; ++i)
		fprintf(fd5,"%s\n",code[family[0][i]-1]);
	fclose(fd5);
	for(i=0;i<nseq-1;++i)
		{
		for(j=i+1;j<nseq;++j)
			fprintf(fd4," %7.2f",D[family[0][i]][family[0][j]]);
		fprintf(fd4,"\n");
		}
	fclose(fd4);
	if((fd6 = fopen("zecret","r")) == NULL)   /*error check*/
		exit(fprintf(stderr,"%s: can't open zecret for reading in main.\n"));
	fscanf(fd6,"%s",orderstring);
	fclose(fd6);
	if((fd7 = fopen("duup","w")) == NULL)   /*error check*/
		exit(fprintf(stderr,"%s: can't open duup for writing in main.\n"));
	fprintf(fd7,"%s\n",orderstring);

	psdmin = psd;
	printf("\nBranching order: ");
	printf("Trial # = %d ",loop);
	if(adj>0)
		printf("**Certain Sreal are adjusted!**");
	printf("\nMain Branch: ");
	side = 0;
	for(i=0;i<=ZGROUP;++i)
		{
		if(side==1)
			printf("\n   Clusters: ");
		else if(side>1)
			printf("\n             ");
		printf("(%s,%s)",code[family[i][0]-1],code[family[i][1]-1]);
		for(v=2;family[i][v]!=0;++v)
			{
			printf("%s)",code[family[i][v]-1]);
			if(v%11==0 && v==nseq-1)
				;
			else if(v%11==0)
				printf("\n             ");
			}
		side++;
		}
	printf("\n");

	printf("Branch lengths: (r1,r2,r3,r4,r5)\n");
	printf("                (r6,r7,........)\n");
	for(j=0; j<ncol; )
		{
		printf("%8.2f ",RB[j]); 
		fprintf(fd7,"%8.2f ",RB[j++]); 
		if(j%5==0)
			{
			printf("\n");
			fprintf(fd7,"\n");
			}
		}
	fclose(fd7);
	printf("\n%% s.d. (observed vs calculated) = %7.2f\n\n",psdmin);

/*	if(loop==0)
		goto out; */

	jump:
	for(i=0; i<flseq+100; ++i)
		{
		get_new_rand:
		rn = (random()/rndmax)*flseq;
	/* If 60% of any column are occupied by X, this column is dropped*/
	/*	for(k=0; k<ndrop_x; ++k)
			if(rn == drop_x[k])
				goto get_new_rand; */
		for(j=0; j<nseq; ++j)
			oseq[j][i] = wseq[j][rn];
		}
	for(j=0; j<nseq; ++j)
		{
		oseq[j][i] = '*';
		oseq[j][i+1] = '\0';
		}
	for(i=0; i<nseq; ++i)
		strcpy(seq[i],oseq[i]);
	}
out:
printf("Total # trial=%d, #_noneg=%3.0f, ",loop,num_noneg); 
if(num_noneg!=0)
	printf("Average %%sd = %6.2f\n",sum_psd/num_noneg);
else
	printf("\n");
fclose(fd9);
remove("albr1");
remove("cecret");
remove("zecret");
remove("albr");
remove("duup");
remove("algr");
}

Read_seq(fd)  /* Enter sequences */
FILE *fd;
{
	int i, ii, m, it;
	char s, c;

	m=nseq=0;
	for( ; ; )
		{
		if(fgetc(fd)==EOF)
			return;
		fgets(code[m],5L,fd); /*Get the 1st 4 letters as code*/
		fgets(seqt[m], sizeof seqt[m], fd);  /*Get sequence title*/

 		seq[m][0]=' ';
		oseq[m][0]=' ';
		clus[m]=fgetc(fd);
		cnum[m]=fgetc(fd)-48;
		it=9;
		if(cnum[m-1]>=9)
			{
			c = fgetc(fd);
			it--;
			if(c>='0' && c<='9')
				cnum[m] = 10*cnum[m] + (c - 48);
			}
		for(i=0; i<it; ++i)
			fgetc(fd);
       		for(i=1; (s=fgetc(fd)) && (s != '*'); ++i)
			{
			if(s=='\n')
				{
				--i;
				for(ii=0; ii<11; ++ii)
					fgetc(fd);
				}
			else if(s==' ')
				--i;
			else
				{
				seq[m][i]=s;
				oseq[m][i]=s;
				}
	       		}
		fgetc(fd);
		lseq[m]=i;
		if(i > MAXSE-1)
			{
			printf("%s length exceeded %d,",code[m],MAXSE);
			printf(" execution terminated!\n");
			printf("Shorten %s or increase MAXAS in program.\n",code[m]);
			exit();
			}
		seq[m][i]='\0';
		oseq[m][i]='\0';
		m++;
		nseq++;
		}
}

Align(si,sj,nsc,cf1,cf2)
int si,sj,nsc,cf1,cf2;
{
	int sigma;   /*Similarity constant between pairs of letters*/
	int Sim;   /*Aligned sequences optimal similarity score*/
	int ssigma;   /*Sum of sigmas*/
	int Mmax,dum1,dum2,jj,kk,kp,m,n,i,j,ij,k,c,t,dcd;
	int mtch,coeff;

	float fidty,lsmall,fsim,NAS;
	float fsigma,sigmaf,test;

   /* Construct the first row and column of the 3 matrices */

	for( j=0; seq[si][j]!='\0'; ++j )
		L[j][0] = R[j][0] = M[j][0] = 0;

	for( k=0; seq[sj][k]!='\0'; ++k )
		L[0][k] = R[0][k] = M[0][k] = 0;


   /* Construct the remaining elements of the 3 matrices */
	for( j=1; seq[si][j]!='\0'; ++j )
		{
		jj = j - 1;
		for( k=1; seq[sj][k]!='\0'; ++k )
			{
			kk = k - 1;
			coeff=1;
			dcd=sigma=mtch=0;

		   /*Scoring current sequence against all earlier sequences*/
			for(i=0; i<sj; ++i)
				{
				sigma += mptable[index[seq[i][j]-'A']]
					        [index[seq[sj][k]-'A']];
				dcd++;
				if(seq[i][j]==seq[sj][k])
					mtch++;
				}

			if(cf1==1 && cf2==1)
				kp = sj + 1;
			else
				kp=sj+2;
			for(ij=0; ij<nsc; ++ij)
				{
				for(i=0; i<sj; ++i)
					{
					dcd++;
					sigma += mptable[index[seq[i][j]-'A']]
					        	[index[seq[kp][k]-'A']];
					if(seq[i][j]==seq[kp][k])
						mtch++;
					}
				kp++;
				}
			if( sj > (wgseq-2) && mtch == dcd )
				coeff = wgt;

			fsigma = (float)sigma/(float)dcd;
			sigma = fsigma;
			sigmaf = sigma;
			test = fsigma - sigmaf;
			if(test>0.5)
				sigma++;
			sigma = coeff*sigma;

			Mmax = M[jj][kk];     /*Determine M[j][k]*/
			dum1 = L[jj][kk] - pen;
			dum2 = R[jj][kk] - pen;
			if( Mmax < dum1 )
				Mmax = dum1;
			if( Mmax < dum2 )
				Mmax = dum2;
			M[j][k] = Mmax + sigma;

			L[j][k] = L[jj][k];   /*Determine L[j][k]*/
			if( L[jj][k] < M[jj][k] )
				L[j][k] = M[jj][k];

			R[j][k] = R[j][kk];    /*Determine R[j][k]*/
			if( R[j][kk] < M[j][kk] )
				R[j][k] = M[j][kk];
			}
		}

   /* Determine the optimal score of the aligned sequences */

	m = strlen(seq[si])-1;
	n = strlen(seq[sj])-1;
	Sim = ( R[m][n] > L[m][n] ) ? R[m][n] : L[m][n];
	Sim = ( Sim > M[m][n] ) ? Sim : M[m][n];

   /* Tracing out the aligned sequences */
	j = strlen(seq[si])-1;
	k = strlen(seq[sj])-1;
	c = 0;
	if( Sim == M[j][k] )
		{
		aseq1[c] = seq[si][j];
		aseq2[c] = seq[sj][k];
		c++;
		branchM(si,sj,&j,&k,&c);
		}

	else if( Sim == L[j][k] )
		{
		aseq1[c] = seq[si][j];
		aseq2[c] = ' ';
		c++;
		branchL(si,sj,&j,&k,&c);
		}

	else
		{
		aseq1[c] = ' ';
		aseq2[c] = seq[sj][k];
		c++;
		branchR(si,sj,&j,&k,&c);
		}

   /* Catch the N-terminal overhangs */
	if( j==0 )
		{
		while( k>-1 )
			{
			aseq1[c] = ' ';
			aseq2[c++] = seq[sj][--k];
			}
		}

	if( k==0 )
		{
		while( j>-1 )
			{
			aseq1[c] = seq[si][--j];
			aseq2[c++] = ' ';
			}
		}

	if( c>MAXAS )   /*error check*/
		exit(printf("job aborted, please increase MAXAS"));

   /* Replace the sequences in the aligned form */
	j=c-3;
	i=1;
	if(si>0)
		{
		while(j>-1)
			{
			if(aseq1[j]==' ')
				seq[si][i]='U';
			else
				seq[si][i]=aseq1[j];
			i++;
			j--;
			}
		seq[si][i]='\0';
		}
	else
		{
		while(j>-1)
			{
			if(aseq1[j]==' ')
				seq[si][i]='X';
			else
				seq[si][i]=aseq1[j];
			i++;
			j--;
			}
		seq[si][i]='\0';
		}

	j=c-3;
	i=1;  /*First element must be blank*/
	while(j>-1)
		{
		if(aseq2[j]==' ')
			seq[sj][i]='X';
		else
			seq[sj][i]=aseq2[j];
		i++;
		j--;
		}
	seq[sj][i]='\0';
}

branchM(si,sj,pj,pk,pc)
int si,sj, *pj, *pk, *pc;
{
	int Mmax;

	if( (*pj==0) || (*pk==0) )
		return;
	(*pj)--;
	(*pk)--;
	Mmax = M[*pj][*pk];
	if( Mmax < L[*pj][*pk] - pen )
		Mmax = L[*pj][*pk] - pen;
	if( Mmax < R[*pj][*pk] - pen )
		Mmax = R[*pj][*pk] -pen;

	if( Mmax == M[*pj][*pk] )
		{
		aseq1[*pc] = seq[si][*pj];
		aseq2[*pc] = seq[sj][*pk];
		(*pc)++;
		branchM(si,sj,pj,pk,pc);
		return;
		}

	else if( Mmax == (L[*pj][*pk] -pen) )
		{
		aseq1[*pc] = seq[si][*pj];
		aseq2[*pc] = ' ';
		(*pc)++;
		branchL(si,sj,pj,pk,pc);
		return;
		}

	else
		{
		aseq1[*pc] = ' ';
		aseq2[*pc] = seq[sj][*pk];
		(*pc)++;
		branchR(si,sj,pj,pk,pc);
		return;
		}
}

branchL(si,sj,pj,pk,pc)
int si,sj, *pj, *pk, *pc;
{
	if( *pj==0 )
		return;
	(*pj)--;
	aseq1[*pc] = seq[si][*pj];

	if( L[*pj][*pk] > M[*pj][*pk] )
		{
		aseq2[*pc] = ' ';
		(*pc)++;
		branchL(si,sj,pj,pk,pc);
		return;
		}

	else
		{
		aseq2[*pc] = seq[sj][*pk];
		(*pc)++;
		branchM(si,sj,pj,pk,pc);
		return;
		}
}

branchR(si,sj,pj,pk,pc)
int si,sj, *pj, *pk, *pc;
{
	if( *pk==0 )
		return;
	(*pk)--;
	aseq2[*pc] = seq[sj][*pk];

	if( R[*pj][*pk] > M[*pj][*pk] )
		{
		aseq1[*pc] = ' ';
		(*pc)++;
		branchR(si,sj,pj,pk,pc);
		return;
		}

	else
		{
		aseq1[*pc] = seq[si][*pj];
		(*pc)++;
		branchM(si,sj,pj,pk,pc);
		return;
		}
}

BkUtoX(si)  /*Change Blank or 'U' to 'X'*/
int si;
{
	int i;

	for(i=1; ;++i)
		{
		if(seq[si][i]=='\0')
			break;
		else if(seq[si][i]==' '||seq[si][i]=='U')
			seq[si][i]='X';
		}
	seq[si][i]='\0';
}

Re_align(si,sj)  /*sj contains 'U' */
int si,sj;
{
	int i, ii;
	char d_asi[MAXAS];  /* Dummy array */

	ii=1;
	d_asi[0]=' ';
	for(i=1; ;++i) /*i=1 because 1st element is blank*/
		{
		if(seq[sj][i]=='\0')
			break;
		else if(seq[sj][i]!='U')
			{
			d_asi[i]=seq[si][ii];
			ii++;
			}
		else if(seq[sj][i]=='U')
			d_asi[i]='U';
		}
	d_asi[i]='\0';
	strcpy(seq[si],d_asi);
	return;
}

format_seq(si)   /*Remove first blank and append '*' in seq[]*/
int si;
{
	int i;
	char d_as[MAXAS];

	for(i=1; ;++i)
		{
		if(seq[si][i]=='\0')
			break;
		else
			d_as[i-1]=seq[si][i];
		}
	d_as[i-1]='*';
	d_as[i]='\0';
	strcpy(seq[si],d_as);
	return;
}

mulalign(fd,NSDL,nax)
int NSDL, nax;
FILE *fd;
{
	int i, j, k, nline;

   /* Copy all sequences into a file */	
	for(i=0; i<nseq; ++i)
		{
		j = k = 0;
		fprintf(fd,"T%s%s",code[i],seqt[i]);
		fprintf(fd,"P%s  %3d",code[i],j+1);
		while( seq[i][j]!='*')
			{
			fprintf(fd,"%2c",seq[i][j++]);
			k++;
			if(k==30)
				{
				fprintf(fd,"\nP%s  %3d",code[i],j+1);
				k = 0;
				}
			}
		fprintf(fd,"%2c\n",'*');
		}

   /*Print progressvely aligned sequences*/
	printf("\t\t\t\t---- Alignment ----\n");
	printf("\t\t\t\t (gap penalty = %d)",pen);
	nline=strlen(seq[0])-1;
	if(nline%nax)
		nline=nline/nax+1;
	else 
		nline=nline/nax;
	
	for(i=0; i<nline; ++i)
		{
		printf("\n");
		for(j=0; j<nseq; ++j)
			{
			printf("\n%s  ",code[j]);
			for(k=i*nax; k<(nax+i*nax) && (seq[j][k]!='*'); ++k)
				{
				if(seq[j][k]=='X'||seq[j][k]=='U')
					printf("%-2c",' ');
				else if( match(j,k)>=NSDL )
					printf("%c_",seq[j][k]);
				else 
					printf("%-2c",seq[j][k]); 
				}
			}
		}
/*	printf("\n\nSequence order after progressive alignment:\n");
	for(i=0; i<nseq; ++i)
		printf("%2d)%s%s",i+1,code[i],seqt[i]);  */
	return;
}

match(j,k)
int j,k;
{
	int i, mtch;

	mtch=0;
	for(i=0; i<nseq; ++i)
		{
		if((seq[j][k]==seq[i][k])&&(seq[j][k]!=' '))
			mtch++;
		}
	return(mtch);
}

control(psi,psj,pcu_nseq,tmp_nseq,nsc,cf1,cf2)  /*Controls recursive loop*/
int *psi,*psj,*pcu_nseq,tmp_nseq,nsc,cf1,cf2;
{
	int i,sip,sjp;
	if( *pcu_nseq<tmp_nseq )
		{
	 	(*psi)++;
		(*psj)++;
		Align(*psi,*psj,nsc,cf1,cf2);
		palign(psi,psj);
		}
	else
		return;
	sip = *pcu_nseq - 1;
	sjp = *pcu_nseq;
	for(i=0; i< *pcu_nseq; ++i)
		BkUtoX(i);
	(*pcu_nseq)++;
	control(&sip,&sjp,pcu_nseq,tmp_nseq,nsc,cf1,cf2);
}

palign(psi,psj)  /*Making sure "once a gap, always a gap"*/
int *psi,*psj;
{
	int i;

	while( *psi>0 )
		{
		(*psi)--;
		(*psj)--;
		Re_align(*psi,*psj);
		palign(psi,psj);
		}
	return;
}

para()
{
	int asl,flag,i,iden,nar,j,k,l,ls,shtseq,rsl[NSEQ];
	double fns,fls,piden[NSEQ][NSEQ],Sreal,Stotal,sum,sum1,sum2;

	for(i=0; i<nseq; ++i)
		{
		lseq[i] = strlen(seq[i])-1;
		asl=ls=0;   /*Actual sequence length*/
		for(j=0; ; ++j)
			{
			if(seq[i][j]=='*')
				break;
			else if(seq[i][j]!='X')
				asl++;
			}
		rsl[i]=asl;
		}

	if(loop==0)
		{
		printf("\n\t\t\t\t-- Results --\n\n");		
		printf("Length of sequences with gaps = %d\n",lseq[0]);
		printf("Average sequence length: ");
		sum=0;
		for(i=0; i<nseq; )
			sum += rsl[i++];
		printf("%4.0f\n\n",sum/nseq);
		}

   /* Zero sreal and piden arrays */
	for(i=0; i<nseq; ++i)
		for(j=0; j<nseq; ++j)
			{
			piden[i][j]=0.0;
			sreal[i][j]=0.0;
			}
	
   /*Compute the similarity (length indep. overhange ignored)
     between each pair*/
	Stotal=0;
	for(i=0; i<nseq-1; ++i)
		{
		for(j=i+1; j<nseq; ++j)
			{
			flag=1;   /*Assuming overhang present at N-term*/
			Sreal=iden=nar=sum1=sum2=0;
			for(k=0; k<lseq[j]; ++k)
				{
				if((seq[i][k]!='X')&&(seq[j][k]!='X'))
					{
					Sreal += mptable[index[seq[i][k]-'A']]
					 	        [index[seq[j][k]-'A']];
					nar++;
					flag=0;
					if(seq[i][k]==seq[j][k])
						iden++;
					}
				else if((flag==0)&&((seq[i][k]!='X')||
						    (seq[j][k]!='X')))
					{
					Sreal -= pen;
					flag=1;
					}
				else if((flag>0)&&((seq[i][k]=='X')&&
						   (seq[j][k]!='X')))
					;
				else if(flag>0)
					;
				sum1 += mptable[index[seq[i][k]-'A']]
					       [index[seq[i][k]-'A']];
				sum2 += mptable[index[seq[j][k]-'A']]
					       [index[seq[j][k]-'A']];
				}
			if(flag>0)  /*True if overhang present at C-term*/
				Sreal += pen;
			Sreal = Sreal/10;
			Stotal += Sreal;
			sreal[i][j] = sreal[j][i] = Sreal;
			siden[i][j] = siden[j][i] = (sum1+sum2)/20;

			piden[i][j] = (float)iden*100.0/(float)nar;
			piden[j][i] = piden[i][j];
			}
		}

	if(loop==0)
		{
		printf("\nSimilarity Scores:\n");
		printf("      ");
		for(i=0; i<nseq; ++i)
			printf("  %s  ",code[i]);
		for(i=0; i<nseq; ++i)
			{
			printf("\n%s ",code[i]);
			for(j=0; j<nseq; ++j)
				printf("%7.2f ",sreal[i][j]);
			}
		printf("\n\t\t*** Sreal(total) = %7.2f ***\n\n",Stotal);
	   
		printf("Percent identity (based on aligned regions):\n");
		printf("      ");
		for(i=0; i<nseq; ++i)
			printf("  %s  ",code[i]);
		for(i=0; i<nseq; ++i)
			{
			printf("\n%s ",code[i]);
			for(j=0; j<nseq; ++j)
				printf("%7.2f ",piden[i][j]);
			}
		printf("\n\n");
		printf("\n\t--- Branching Order and Branch Lengths ");
		printf("Determined by Bootstrap Method ---\n");
		}
	return;   
}

double sim_score1(tmp_nseq)
int tmp_nseq;
{
	int tmp_lseq,flag,i,j,k,nar;
	float Sreal,Stotal;

	tmp_lseq = strlen(seq[0]); /*Length of seq[0] = tmp_lseq+1*/

   /*Compute the similarity (length indep. overhange ignored)
     between each pair*/
	Stotal=nar=0;
	for(i=0; i<tmp_nseq-1; ++i)
		{
		for(j=i+1; j<tmp_nseq; ++j)
			{
			flag=1;   /*Assuming overhang present at N-term*/
			Sreal=0;
			for(k=1; k<tmp_lseq; ++k)
				{
				if((seq[i][k]!='X')&&(seq[j][k]!='X'))
					{
					Sreal += mptable[index[seq[i][k]-'A']]
					 	        [index[seq[j][k]-'A']];
					nar++;
					flag=0;
					}
				else if((flag==0)&&((seq[i][k]!='X')||
						    (seq[j][k]!='X')))
					{
					Sreal -= pen;
					flag=1;
					}
				else if((flag>0)&&((seq[i][k]=='X')&&
						   (seq[j][k]!='X')))
					;
				else if(flag>0)
					;
				}
			if(flag>0)  /*True if overhang present at C-term*/
				Sreal += pen;
			Sreal = Sreal/10;
			Stotal += Sreal;
			}
		}
	Stotal = Stotal/nar;
	return(Stotal);		
}

switch_array(tmp_nseq,array)
int tmp_nseq;
char (*array)[MAXAS];
{
	int i;
	char dum_array[MAXAS];

	strcpy(dum_array,array[tmp_nseq-1]);
	strcpy(array[tmp_nseq-1],array[tmp_nseq-2]);
	strcpy(array[tmp_nseq-2],dum_array);
}

rm_cX1(tmp_nseq,array)
int tmp_nseq;
char (*array)[MAXAS];
{
	int flag,i,j,m,n;

   /*Removing X's that appear in all sequences at the same location*/
	for(i=1 ; array[0][i]!='\0'; ++i)
		{
		flag=0;  /*means all 'X'*/
		for(j=0; j<tmp_nseq; ++j)
			if(array[j][i]!='X')
				flag=1;  /*means not all 'X'*/
		if(flag==0)
			{
			for(n=0; n<tmp_nseq; ++n)
				{
				for(m=i; array[n][m]!='\0'; ++m)
					array[n][m]=array[n][m+1];
				array[n][m]='\0';
				}
			i--;
			}
		}
}

shuffle() /*Calculate Srand*/
{
	int i,j,tlseq1,tlseq2,dum;
	int ii,jj;
	double rand_score,rand_score1;
	char aa[27] = "ABCDEFGHIJKLMNOPQRSTUVWXYZ";
	int h1[26],h2[26];
	int ttlseq1,ttlseq2;
	int gaps1[20],gaps2[20],gps1,gps2,gflag,tngap,tlgap;

	for(i=0; i<nseq; ++i)
		{
		for(j=0; seq[i][j] != '*'; ++j)
			;
		seq[i][j] = '\0';
		lseq[i] = strlen(seq[i]);
		}

	for(i=0; j<nseq; ++i)
		for(j=0; j<nseq; ++j)
			shmean[i][j] = 0.0;

	for(i=0; i<nseq-1; ++i)
		{
		for(j=i+1; j<nseq; ++j)
			{
			strcpy(seq1,seq[i]);
			strcpy(seq2,seq[j]);
			tlseq1 = lseq[i];
			tlseq2 = lseq[j];
			dum = 0;
		/*Common gaps in seq1 and seq2 are removed*/
			while(seq1[dum]!='\0')
				{
				if(seq1[dum]=='X' && seq2[dum]=='X')
					{
					rm_cX2(&tlseq1,&tlseq2);
					break;
					}
				dum++;
				}
			for(ii=0; ii<20; ++ii)
				gaps1[ii] = gaps2[ii] = 0;
			gps1=0;
			gflag=0;
			for(ii=0; ii<tlseq1;)
				{
				while(seq1[ii]=='X')
					{
					gaps1[gps1]++;
					gflag++;
					ii++;
					}
				if(gflag>0)
					{
					gps1++;
					gflag=0;
					}
				ii++;
				}
			gps2=0;
			gflag=0;
			for(ii=0; ii<tlseq2;)
				{
				while(seq2[ii]=='X')
					{
					gaps2[gps2]++;
					gflag++;
					ii++;
					}
				if(gflag>0)
					{
					gps2++;
					gflag=0;
					}
				ii++;
				}
			rand_score = 0;
			for(ii=0; ii<26; ++ii)
				h1[ii] = h2[ii] = 0;
			ttlseq1 = 0;
			for(ii=0; ii<tlseq1; ++ii)
				{
				if(seq1[ii] != 'X' && seq2[ii] != 'X')
					{
					h1[seq1[ii]-'A']++;
					ttlseq1++;
					}
				if(seq1[ii] != 'X' && seq2[ii] != 'X')
					{
					h2[seq2[ii]-'A']++;
					}
				}
			tlgap=tngap=0;
			for(ii=0; ii<20; ++ii)
				if(gaps1[ii]>0)
					{
					tlgap += gaps1[ii];
					tngap++;
					}
			for(ii=0; ii<20; ++ii)
				if(gaps2[ii]>0)
					{
					tlgap += gaps2[ii];
					tngap++;
					}	
			if(seq1[0]=='X')
				tngap--;
			if(seq1[tlseq1-1]=='X')
				tngap--;
			if(seq2[0]=='X')
				tngap--;
			if(seq2[tlseq1-1]=='X')
				tngap--;
			for(ii=0; ii<26; ++ii)
				for(jj=0; jj<26; ++jj)
					rand_score += (float)mptable[index[aa
					[ii]-'A']][index[aa[jj]-'A']]*h1[ii]
					*h2[jj];
			rand_score1 = ((rand_score/ttlseq1)-pen*tngap)/10.0;
			shmean[i][j] = shmean[j][i] = rand_score1;
			}
		}
}

rm_cX2(plseq1,plseq2)
int *plseq1, *plseq2;
{
	int i,k,kp;

	for(i=0; seq1[i]!='\0'; ++i)
		{
		if(seq1[i]=='X' && seq1[i]==seq2[i])
			{
			k=i;
			while(seq1[k+1]!='\0')
				{
				seq1[k]=seq1[k+1];
				k++;
				}
			seq1[k]='\0';
			kp=i;
			while(seq2[kp+1]!='\0')
				{
				seq2[kp]=seq2[kp+1];
				kp++;
				}
			seq2[kp]='\0';
			i--;
			}
	/* If seq have no X's, have to take into account here */
		}
	*plseq1 = k;
	*plseq2 = kp;
	return;
}

rpJwX()   /*Replace J in sequences with X*/
{
	int i,j;
	for(i=0; i<nseq; ++i)
		{
		for(j=0; seq[i][j]!='\0'; ++j)
			if(seq[i][j]=='J')
				seq[i][j]='X';
		}
}

realign_cl(si,sj)  /* seq[si] contains 'X' & `J`; seq[sj] contains 'J'*/
int si,sj;
{
	int i, ii;
	char d_asi[MAXAS];  /* Dummy array */

	ii=1;
	d_asi[0] = ' ';
	for(i=1;seq[si][i]!='\0';++i) /*i=1 because 1st element is blank*/
		{
		if(seq[si][i]=='X')
			d_asi[i] = 'X';
		else
			d_asi[i] = seq[sj][ii++];
		}
	d_asi[i]='\0';
	strcpy(seq[sj],d_asi);
	return;
}

bord()
{
	int i,ii,ii1,j,k,q,qz,ss1,v,vv,w,z,zt,zz;
	int minj1,mink1,n,nj,r,s,t;

   /*NULLified the array family for repeated usage*/
	for(i=0; i<10; ++i)
		{
		for(j=0; j<NSEQ; ++j)
			family[i][j] = 0;
		}	

	i=qz=1;
	i=s=ZGROUP=0;
	for(j=1;j<nseq;++j)
		{
		for(k=j+1;k<=nseq;++k)
			pro[j][k]=pro[k][j]=D[j][k];
		}

	FINDMIN(&minj1, &mink1);
	zt=1;
	ss1=nseq-1;
	pro[mink1][minj1]=pro[minj1][mink1]=1000;
	family[0][0]=minj1;
	family[0][1]=mink1;
	for(zz=1;zz<=ss1;++zz)
		{
		NEWMATRIX();
		FINDMIN(&minj1, &mink1);
		FAMFIND(&ZGROUP, &ss1, minj1, mink1);
		if(family[qz][0]!=0)
			{
			zz--;
			qz++;
			}
		if(zz== ss1-1 && zt==1)
			return;
		}
}

FINDMIN(pminj1,pmink1)
int *pminj1, *pmink1;
{
	int j,k;

	j=1;
	k=2;
	*pminj1 = j;
	*pmink1 = k;
	while(j<nseq)
		{
		while(k<nseq)
			{
			if(pro[*pminj1][*pmink1] > pro[j][++k])
				{
				*pminj1 = j;
				*pmink1 = k;
				}
			else 
				continue;
			}
		++j;
		k=j;
		}
	return;
}

NEWMATRIX()
{
	int i,j,k,q,qq1,qq2;

	j=1;
	k=2;
	while(j<nseq)
		{
		while(k<=nseq)
			{
			for(i=0;family[i][0]!=0;++i)
				{
				for(q=0;family[i][q]!=0;q++)
					{
					if(pro[j][k]==1000)
						goto NEWK;
					else if(j==family[i][q])
						{
						qq1=k;
						qq2=j;
						AV(qq1,qq2,i);
						goto NEWK;
						}
					else if(k==family[i][q])
						{
						qq1=j;
						qq2=k;
						AV(qq1,qq2,i);
						goto NEWK;
						}
					else 
						continue;
					}
				}
	 		NEWK: ++k;
			}
		j++;
		k=j+1;
		}
	return;
}

FAMFIND(pZ,pss1,minj1,mink1)
int *pZ, *pss1, minj1, mink1;
{
	int i,q,q1,W;

	for(i=0;family[i][0]!=0;++i)
		{
  		for(q=0;family[i][q]!=0;++q)
			{
    			if(family[i][q]==minj1)
				{
				q1=mink1;
				OTHER(q1,i);
				goto tree;
				}
    			if(family[i][q]==mink1)
				{
				q1=minj1;
				OTHER(q1,i);
				goto tree;
				}
    			}
		}
	(*pZ)++;
	if(W=1)
		W=0;
		/*A new pair to be placed into a new family*/

	pro[minj1][mink1]=pro[mink1][minj1]=1000;
	family[*pZ][0]=minj1;
	family[*pZ][1]=mink1;
	W++;
	if(W==1)
		(*pss1)--;
	tree: return;
}

OTHER(q1,i)
int q1,i;
{
		/*Checks to see if the other min. value belongs*/
		/*to another family*/
	int ii,ii1,v,vv;

	for(ii=0;family[ii][0]!=0;ii++)
		{
		for(ii1=0;family[ii][ii1]!=0;++ii1)
			{
			if(q1==family[ii][ii1])
				{
				ANOTHER(&ii,i);
				return;
				}
			}
		}

	/*Does not belong to a new family*/

	for(v=0;v<MAXAS;++v)
		{
		if(family[i][v]==0)
			{
			family[i][v]=q1;
			break;
			}
		}
	for(v=0;family[i][v+1]!=0;++v)
		{
		for(vv=1;family[i][vv]!=0;++vv)
			{
			if(family[i][v]!=family[i][vv])
				pro[family[i][v]][family[i][vv]]=
				pro[family[i][vv]][family[i][v]]=1000;
			}
		}
	return;
}

ANOTHER(pii,i)
int *pii, i;
{
	      /*Min. does belong to another family...add to*/
	      /*family of first min. value*/
	int v,vv;

	for(v=0;v<MAXAS;v++)
		{
		if(family[i][v]==0)
			{
			for(vv=0;family[*pii][vv]!=0;++vv)
				family[i][v++]=family[*pii][vv];
			break;
			}
		}

	for(v=0;family[i][v+1]!=0;v++)
		{
		for(vv=1;family[i][vv]!=0;++vv)
			{
			if(family[i][v]!=family[i][vv])
				pro[family[i][v]][family[i][vv]]=
				pro[family[i][vv]][family[i][v]]=1000;
			}
		}
	return;
}

AV(qq1,qq2,i)
int qq1,qq2,i;
{
	 /*Determines whether qq1 also belongs to a*/
	 /*family; if so, average throughout both families*/

	int ii,ii1,m,n,q,v;

	m=n=v=0;
	pro[m][n]=0;
	for(ii=0;family[ii][0]!=0;++ii)
		{
		for(ii1=0;family[ii][ii1]!=0;++ii1)
			{
			if(qq1==family[ii][ii1])
				goto heman;
			else 
				continue;
			}
		}
	for(q=0;family[i][q]!=0;++q)
		{
		++v;
		pro[m][n] += D[family[i][q]][qq1];
		}
	pro[m][n]=pro[m][n]/v;
	pro[qq1][qq2]=pro[qq2][qq1]=pro[m][n];
	goto heman2;

	    /*qq1 belongs to a separate family average*/
	    /*thoughout both families*/

	heman:	for(ii1=0;family[ii][ii1]!=0;ii1++)
			{
			for(q=0;family[i][q]!=0;++q)
				{
				v++;
				pro[m][n] += D[family[i][q]]
						 [family[ii][ii1]];
				}
			}
		pro[m][n]=pro[m][n]/v;
		for(ii1=0;family[ii][ii1]!=0;++ii1)
			{
			for(q=0;family[i][q]!=0;++q)
				pro[family[ii][ii1]][family[i][q]]=
				pro[family[i][q]][family[ii][ii1]]=pro[m][n];
			}

	heman2: return;
}

int NumStrings;
char InputString[50][60];
blen()
{
	FILE *B1andD, *B2andD, *fopen();
	int i,length,c;
	char line[400],paren[200];

	if((B1andD = fopen("albr1","r")) == NULL)   /*error check*/
		exit(fprintf(stderr,"%s: can't open albr1 for reading in blen.\n"));

	i=0;
	while((c=getc(B1andD)) != '*')
		{
		ungetc(c,B1andD);
		fscanf(B1andD,"%s",InputString[i]);
		i++;
		getc(B1andD);
		}
	NumStrings = i;
	GetParen(InputString,paren);
	length = strlen(paren);
	for(i=1; i<(length-1); ++i)
		paren[i-1] = paren[i];
	paren[i-1] = '\0';

	if((B2andD = fopen("zecret","w")) == NULL)   /*error check*/
		exit(fprintf(stderr,"%s: can't open zecret for writing in blen.\n"));

	fprintf(B2andD,"%s\n",paren);
	while((c=fgetc(B1andD)) != EOF)
		{
		fgets(line, sizeof line, B1andD);
		fprintf(B2andD,"%s",line);
		}
	fclose(B2andD);
	fclose(B1andD);
	branchlength();
	return;
}

GetParen(start,paren)
char start[50][60], paren[];
{
	int i,j,match;
	char c[2],SubParen[100],oldparen[200];

/*printf("\tstart: %s\n",start[0]);*/
	paren[0] = '(';
	for(i=1; i<3; ++i)
		paren[i] = start[0][i-1];
	paren[i++] = ')';
	paren[i] = '\0';
/*	printf("**%s\n",paren);*/

	i=2;
	while(i<strlen(start[0]))
		{
		c[0] = start[0][i];
		c[1] = NULL;
		match = 0;
		if(NumStrings>1)
			{
			for(j=1; j<NumStrings; ++j)
				{
				if(c[0] == InputString[j][0])
					{
					match = 1;
					break;
					}
				}
			}
		if(match==1)
			{
			GetParen((InputString+j),SubParen);
			strcpy(oldparen,"(");
			strcat(oldparen,paren);
			strcat(oldparen,SubParen);
			strcat(oldparen,")");
			strcpy(paren,oldparen);
			i += strlen(InputString[j]);
			}
		else
			{
			strcpy(oldparen,"(");
			strcat(oldparen,paren);
			strcat(oldparen,c);
			strcat(oldparen,")");
			strcpy(paren,oldparen);
			i++;
			}
		}
}

char branch[100][100],res1[1],res2[1];
char parameter[100], FoundLetters[100];
int ParameterIndex, CurrentBranch, LastBranch;
int bnum[100],or[60],lastor,nseq;

branchlength()
{
	FILE *fd1,*fd2,*fopen();
	char FoundBranch[100],*pts;
	int i,ilimit,j,k,flag;
	double D[50][50],x;

	ParameterIndex = 0;
	CurrentBranch = 0;
	LastBranch = 0;

   /*Read branching order*/
	if((fd1 = fopen("zecret","r")) == NULL)   /*error check*/
		exit(fprintf(stderr,"%s: can't open zecret for reading in branchlength.\n"));

	fscanf(fd1,"%s",parameter);

	nseq = 0;
	pts = parameter;
	while(*pts)
		{
		if(*pts>='A' && *pts<='Z')
			{
			or[nseq] = *pts - 65;
			lastor = or[nseq] + 1;
			nseq++;
			}
		else if(*pts>='a' && *pts<='z')
			{
			or[nseq] = *pts - 97 + lastor;
			nseq++;
			}
		pts++;
		}

   /*Read distance scores*/
	for(i=0; i<nseq-1; ++i)
		for(j=i+1; j<nseq; ++j)
			{
			fscanf(fd1,"%lf",&x);
			D[i][j] = D[j][i] = x;
			}

	GetBranch(FoundBranch);
	LastBranch--;
	for(i=0; i<LastBranch; ++i)
		bnum[i] = i+1;

	if((fd2 = fopen("cecret","w")) == NULL)   /*error check*/
		exit(fprintf(stderr,"%s: can't open cecret for writing in branchlength\n"));
	ilimit = strlen(FoundBranch);	
	for(i=0; i<ilimit-1; ++i)
		{
		for(j=i+1; j<ilimit; ++j)
			{
			for(k=0; k<LastBranch; ++k)
				{
				res1[0] = FoundBranch[i];
				res2[0] = FoundBranch[j];
				flag = compare(k);
				if(flag==2)
					break;
				else if(flag==0)
					;
				else
					fprintf(fd2," %d",bnum[k]);
				}
			fprintf(fd2,"\n");
			}
		}

   /*Put distance scores in the right order*/	
	for(i=0; i<nseq-1; ++i)
		{
		for(j=i+1; j<nseq; ++j)
			fprintf(fd2,"%7.2f",D[or[i]][or[j]]);
		fprintf(fd2,"\n");
		}
	fclose(fd1);
	fclose(fd2);
	brchlen();
}

GetBranch(FoundBranch)
char FoundBranch[];
{
	char LeftBranch[100], RightBranch[100];

   /*Find LEFT side of branch*/
	LeftBranch[0] = parameter[ParameterIndex++];
	LeftBranch[1] = '\0';
	if(LeftBranch[0]=='(')
		GetBranch(LeftBranch);
	strcpy(branch[CurrentBranch],LeftBranch);
	CurrentBranch++;
	LastBranch++;

   /*Find RIGHT side of branch*/
	RightBranch[0] = parameter[ParameterIndex++];
	RightBranch[1] = '\0';
	if(RightBranch[0]=='(')
		GetBranch(RightBranch);
	strcpy(branch[CurrentBranch],RightBranch);
	CurrentBranch++;
	LastBranch++;

   /*Remove trailing ")"*/
	if(parameter[ParameterIndex]==')')
		ParameterIndex++;
	
	strcpy(FoundBranch,LeftBranch);
	strcat(FoundBranch,RightBranch);
	return;
}

compare(n)
int n;
{
	int match,i;

	match = 0;	
	for(i=0; branch[n][i]!='\0'; ++i)
		if(res1[0] == branch[n][i])
			{
			match++;
			break;
			}

	for(i=0; branch[n][i]!='\0'; ++i)
		if(res2[0] == branch[n][i])
			{
			match++;
			break;
			}
	return(match);
}


double a[nsq][tpnsq], iA[tnsq][tnsq];

brchlen()
{
	FILE *fd;
	int i, j, k, kk, pt[30];
	char c;	
	double temp;
	double cum,x;

	if((fd = fopen("cecret","r")) == NULL)   /*error check*/
		exit(fprintf(stderr,"%s: can't open cecret for reading in brchlen.\n"));

	nrow=nseq;
	ncol=2*nrow-3;
	nrow=nrow*(nrow-1)/2;

	for(i=0; i<nrow; ++i) 
		{
		k=0;
		while((c=fgetc(fd))!='\n')
			{
			ungetc(c,fd);
			fscanf(fd,"%d",&pt[k]);
			k++;
			}
		sort(pt,k);   /*Arrange numbers in pt in ascending order*/
		for(j=0; j<ncol; ++j)
			{
			A[i][j]=0;
			for(kk=0; kk<k; ++kk)
				if((j+1)==pt[kk])
					A[i][j]=1;
			}
		}

	for(i=0; i<nrow; ++i)
		{
		fscanf(fd,"%lf",&x);
		Y[i]=x;
		}
	fclose(fd);

   	for(i=0; i<ncol; ++i)
		for(j=0; j<ncol; ++j)
			{
			temp=0;
			for(k=0; k<nrow; ++k)
				temp += A[k][i]*A[k][j];
			a[i][j]=temp;
			}

	inverse(ncol, 2*ncol);
	for(i=0; i<nrow; ++i)
		for(j=0; j<ncol; ++j)
			{
			cum=0;
			for(k=0; k<ncol; ++k)
				cum += A[i][k]*iA[k][j];
			a[i][j]=cum;
			}

	for(j=0; j<ncol; ++j)
		{
		cum=0;
		for(k=0; k<nrow; ++k)
			cum += Y[k]*a[k][j];
		RB[j]=cum;
		}
}

inverse(ns,tns)
int ns, tns;
{
	int i, ii, ip, j, jj, k, l, n, np;
	double cum, den, dum, m, s;
	n=ns;

	for(i=0; i<n; i++)
		{
		ip=i+n;
		for(j=n; j<2*n; j++)
			{
			if(ip==j)
				a[i][j]=1;
			else
				a[i][j]=0;
			}
		}
	for(k=0; k<n-1; ++k)
		{
		l=k;
		for(i=k+1; i<n; i++)
			{
			if(fabs(a[i][k])>fabs(a[l][k]))
				l=i;
			}
		if(l!=k)
			{
			for(j=k; j<2*n; j++)
				{
				dum=a[k][j];
				a[k][j]=a[l][j];
				a[l][j]=dum;
				}
			}

		for(i=k+1; i<n; ++i)
			{
			m=a[i][k]/a[k][k];
			a[i][k]=0.0;
			j=k+1;
			for(; j<2*n; ++j)
				a[i][j] -= m*a[k][j];
			}
		}


	np=n-1;
	den=a[np][np];
	for(i=n; i<tns; ++i)
		{
		a[np][i]=a[np][i]/a[np][np];	
		iA[np][i-n]=a[np][i];
		}

	for(i=np-1; i>=0; --i)
		for(k=n; k<tns; ++k)
			{
			s=0;
			for(j=i+1; j<=np; ++j)
				s += a[i][j]*a[j][k];
			a[i][k]=(a[i][k]-s)/a[i][i];
			iA[i][k-n]=a[i][k];
			}
}

sort(a,n)
int a[], n;
{
	int i,j,temp;

	for(i=0; i<n-1; ++i)
		for(j=i+1; j<n; ++j)
			if(a[i]>a[j])
				{
				temp = a[i];
				a[i] = a[j];
				a[j] = temp;
				}
}