/*This is the papa3.c program.
Reference: R.F. Doolittle & D.F. Feng (1990) Meth. Enzymol., 183:659-669. 


Note to compile this program on the Iris Indigo, type
   cc -cckr -signed papa3.c						*/



#define MAXS 500   /*Maximum size for sequence*/
#define NSEQ 50    /*Maximum number of sequences*/

#include <stdio.h>
int nseq;
char seq[NSEQ][MAXS];
#define nse NSEQ*(NSEQ-1)/2
float mat[nse];
int kind[14];
int xtable[5][14] =	{
		{30, 0, 0, 0, 0,15,15, 0, 0, 0,20,30,30,18},
		{ 0,30, 0, 0, 0,15,15, 0,30,30,20, 0, 0,18},
		{ 0, 0, 0, 0,30, 0, 0,20, 0, 0,20, 0, 0,18},
		{ 0, 0,30, 0, 0,15,15,20, 0,30, 0,30, 0,18}, 
		{ 0, 0, 0,30, 0,15,15,20,30, 0, 0, 0,30,18}, 
			};
int ytable[5][14] =	{
		{30, 0, 0, 0,15, 0,15, 0, 0, 0,30,20,30,18},
		{ 0, 0,30, 0,15, 0,15,30, 0,30, 0,20, 0,18}, 
		{ 0, 0, 0, 0, 0,30, 0, 0,20, 0, 0,20, 0,18},
		{ 0,30, 0, 0,15, 0,15, 0,20,30,30, 0, 0,18}, 
		{ 0, 0, 0,30,15, 0,15,30,20, 0, 0, 0,30,18}, 
			};
int ztable[5][14] =	{
		{30, 0, 0, 0,15,15, 0, 0, 0, 0,30,30,20,18}, 
		{ 0, 0, 0,30,15,15, 0,30,30, 0, 0, 0,20,18}, 
		{ 0, 0, 0, 0, 0, 0,30, 0, 0,20, 0, 0,20,18},
		{ 0,30, 0, 0,15,15, 0, 0,30,20,30, 0, 0,18}, 
		{ 0, 0,30, 0,15,15, 0,30, 0,20, 0,30, 0,18}, 
			};
float xdtable[5][14] =	{
		{1,0,0,0,0,.5,.5, 0,0,0,.7,1,1,.6},
		{0,1,0,0,0,.5,.5, 0,1,1,.7,0,0,.6},
		{0,0,0,0,1, 0, 0,.7,0,0,.7,0,0,.6},
		{0,0,1,0,0,.5,.5,.7,0,1, 0,1,0,.6}, 
		{0,0,0,1,0,.5,.5,.7,1,0, 0,0,1,.6}, 
			};
float ydtable[5][14] =	{
		{1,0,0,0,.5,0,.5,0, 0,0,1,.7,1,.6},
		{0,0,1,0,.5,0,.5,1, 0,1,0,.7,0,.6}, 
		{0,0,0,0, 0,1, 0,0,.7,0,0,.7,0,.6},
		{0,1,0,0,.5,0,.5,0,.7,1,1, 0,0,.6}, 
		{0,0,0,1,.5,0,.5,1,.7,0,0, 0,1,.6}, 
			};
float zdtable[5][14] =	{
		{1,0,0,0,.5,.5,0,0,0, 0,1,1,.7,.6}, 
		{0,0,0,1,.5,.5,0,1,1, 0,0,0,.7,.6}, 
		{0,0,0,0, 0, 0,1,0,0,.7,0,0,.7,.6},
		{0,1,0,0,.5,.5,0,0,1,.7,1,0, 0,.6}, 
		{0,0,1,0,.5,.5,0,1,0,.7,0,1, 0,.6}, 
			};

main(ac,av)
int ac;
char *av[];
{
	FILE *fd, *fsd;	
	int flag, i, ii, j, jj, k, l, m, n;
	int O,Q,sum,pprt;
	int item,nterm,index;
	int ip,jp,kp,lp,ipp,ipr,iprr;
	float max,cm,X,Y,Z;
	float length[5];

 	char c, stitle[NSEQ][100], code[NSEQ][5];
	char s1,s2,s3,s4;
	
	if(ac<2)
		{
		fprintf(stderr,"Parsimony After Progressiv Alignment (PAPA)\n");
		fprintf(stderr,"usage: %s rf [opt] [>wf]\n",av[0]);
		fprintf(stderr,"\trf: Aligned sequences in NEWAT format.\n");
		fprintf(stderr,"\topt:Set opt=1 for intermed steps printout\n");
		fprintf(stderr,"\twf: Summary Results and Branching Order!\n");
		fprintf(stderr,"Current limits: NSEQ=<%d, length=<%d\n",
				NSEQ,MAXS);
		fprintf(stderr,"X,Y,or Z tree selected on max c-length\n");
		fprintf(stderr,"Branch order: FM-method using diff matrix.\n");
		exit();
		}

	if((fd=fopen(av[1],"r"))==NULL)
		exit("Can't open %s\n",av[1]);

	if((pprt = atoi(av[2]))!=0)
		pprt = 1;

	nseq=0;
	for(i=0; ; ++i)
		{
		if((c=getc(fd))==EOF)
			break;
		ungetc(c,fd);
		nseq++;
		while((c=getc(fd))==' ')
			;
		fgets(stitle[i], sizeof stitle[i], fd);
		fgetc(fd);
		fgets(code[i],5,fd);
		code[i][5] = '\0';
		for(k=0; k<6; ++k)
			fgetc(fd);
	       	for(j=0;(c=getc(fd)) && c != '*'; ++j)
			{
			if(c=='\n')
				{
				for(k=0; k<11; ++k)
					fgetc(fd);
				--j;
				continue;
				}
			if(c==' ')
				{
				--j;
				continue;
				}
			seq[i][j]=c;
		       }
		seq[i][j]='*';
		seq[i][j+1]='\0';
		getc(fd);
		}

	for(i=0; i<nseq; ++i)
		{
		if(i<26)
			printf("%c) %s",65+i,stitle[i]);
		else if(i<52)
			printf("%c) %s",71+i,stitle[i]);
		else
			exit(printf(stderr,"Too many sequences\n"));
		}

	for(i=0; seq[0][i]!='\0'; ++i)
		{
		flag=0;
		for(j=0; j<nseq; ++j)
			if(seq[j][i]!='X')
				flag=1;
		if(flag==0)
			{
			for(n=0; n<nseq; ++n)
				{
				for(m=i; seq[n][m]!='\0'; ++m)
					seq[n][m]=seq[n][m+1];
				seq[n][m]='\0';
				}
			i--;
			}
		}

	ipp = item = 0;
	nterm = nseq*(nseq-1)/2;
	for(i=0; i<nterm+1; ++i)
		mat[i] = 0;
	if(pprt != 0)
		printf("Comb Tree     a        b        c        d        e\n");
	for(i=0; i<nseq; ++i)
	for(j=i+1; j<nseq; ++j)
	for(k=j+1; k<nseq; ++k)
	for(l=k+1; l<nseq; ++l)
		{
		ipp++;
		for(ii=0; ii<14; ++ii)
			kind[ii]=0;
		O=Q=0;
		for(m=0; seq[0][m]!='*'; )
			{
			s1 = seq[i][m];
			s2 = seq[j][m];
			s3 = seq[k][m];
			s4 = seq[l][m];
			if(s1=='X'||s2=='X'||s3=='X'||s4=='X')
				Q++;
			else if(s1==s2 && s1==s3 && s1==s4)
				{
				if(s1=='X')
					;
				else
					O++;
				}
			else if(s1==s2 && s3==s4 && s1!=s3)
				kind[4]++;
			else if(s1==s3 && s2==s4 && s1!=s2)
				kind[5]++;
			else if(s1==s4 && s2==s3 && s1!=s2)
				kind[6]++;
			else if(s2==s3 && s2==s4 && s1!=s2)
				kind[0]++;
			else if(s1==s2 && s1==s3 && s1!=s4)
				kind[3]++;
			else if(s1==s2 && s1==s4 && s3!=s4)
				kind[2]++;
			else if(s1==s3 && s1==s4 && s1!=s2)
				kind[1]++;
			else if(s1==s2 && s1!=s3 && s1!=s4 && s3!=s4)
				kind[7]++;
			else if(s1==s3 && s1!=s2 && s1!=s4 && s2!=s4)
				kind[8]++;
			else if(s1==s4 && s1!=s2 && s1!=s3 && s2!=s3)
				kind[9]++;
			else if(s3==s4 && s1!=s2 && s1!=s3 && s2!=s3)
				kind[10]++;
			else if(s2==s4 && s1!=s2 && s1!=s3 && s1!=s3)
				kind[11]++;
			else if(s2==s3 && s2!=s1 && s2!=s4 && s1!=s4)
				kind[12]++;
			else 
				kind[13]++;
			m++;
			}
		X=Y=Z=0;
		for(jj=0; jj<14; ++jj)
			X += xdtable[2][jj]*kind[jj];
		for(jj=0; jj<14; ++jj)
			Y += ydtable[2][jj]*kind[jj];
		for(jj=0; jj<14; ++jj)
			Z += zdtable[2][jj]*kind[jj];

		ip = i+1;
		jp = j+1;
		kp = k+1;
		lp = l+1;
		max = X;
		max = (max>Y)?max:Y;
		max = (max>Z)?max:Z;
		if(max==X)
			{
			if(pprt != 0)
		   	printf("%2d,%2d,%2d,%2d     X    ",i+1,j+1,k+1,l+1);
			for(ii=0; ii<5; ++ii)
				{
				cm=0;
				for(jj=0; jj<14; ++jj)
					cm += xdtable[ii][jj]*kind[jj];
				length[ii] = cm;
				if(pprt != 0)
					printf("%7.1f  ",cm);
				}
			if(pprt != 0)
				printf("\n");
			index=(ip-1)*(2*nseq-ip)/2+jp-ip;
			if(index==1)
				item++;
			mat[index] += length[0]+length[1];
			index=(ip-1)*(2*nseq-ip)/2+kp-ip;
			mat[index] += length[0]+length[2]+length[3];
			index=(ip-1)*(2*nseq-ip)/2+lp-ip;
			mat[index] += length[0]+length[2]+length[4];
			index=(jp-1)*(2*nseq-jp)/2+kp-jp;
			mat[index] += length[1]+length[2]+length[3];
			index=(jp-1)*(2*nseq-jp)/2+lp-jp;
			mat[index] += length[1]+length[2]+length[4];
			index=(kp-1)*(2*nseq-kp)/2+lp-kp;
			mat[index] += length[3]+length[4];
			}
		else if(max==Y)
			{
			if(pprt != 0)
			printf("%2d,%2d,%2d,%2d     Y    ",i+1,j+1,k+1,l+1);
			for(ii=0; ii<5; ++ii)
				{
				cm=0;
				for(jj=0; jj<14; ++jj)
					cm += ydtable[ii][jj]*kind[jj];
				length[ii] = cm;
				if(pprt != 0)
					printf("%7.1f  ",cm);
				}
			if(pprt != 0)
				printf("\n");
			index=(ip-1)*(2*nseq-ip)/2+jp-ip;
			if(index==1)
				item++;
			mat[index] += length[0]+length[2]+length[3];
			index=(ip-1)*(2*nseq-ip)/2+kp-ip;
			mat[index] += length[0]+length[1];
			index=(ip-1)*(2*nseq-ip)/2+lp-ip;
			mat[index] += length[0]+length[2]+length[4];
			index=(jp-1)*(2*nseq-jp)/2+kp-jp;
			mat[index] += length[1]+length[2]+length[3];
			index=(jp-1)*(2*nseq-jp)/2+lp-jp;
			mat[index] += length[3]+length[4];
			index=(kp-1)*(2*nseq-kp)/2+lp-kp;
			mat[index] += length[1]+length[2]+length[4];
			}
		else
			{
			if(pprt != 0)
			printf("%2d,%2d,%2d,%2d     Z    ",i+1,j+1,k+1,l+1);
			for(ii=0; ii<5; ++ii)
				{
				cm=0;
				for(jj=0; jj<14; ++jj)
					cm += zdtable[ii][jj]*kind[jj];
				length[ii] = cm;
				if(pprt != 0)
					printf("%7.1f  ",cm);
				}
			if(pprt != 0)	
				printf("\n");
			index=(ip-1)*(2*nseq-ip)/2+jp-ip;
			if(index==1)
				item++;
			mat[index] += length[0]+length[2]+length[3];
			index=(ip-1)*(2*nseq-ip)/2+kp-ip;
			mat[index] += length[0]+length[2]+length[4];
			index=(ip-1)*(2*nseq-ip)/2+lp-ip;
			mat[index] += length[0]+length[1];
			index=(jp-1)*(2*nseq-jp)/2+kp-jp;
			mat[index] += length[3]+length[4];
			index=(jp-1)*(2*nseq-jp)/2+lp-jp;
			mat[index] += length[1]+length[2]+length[3];
			index=(kp-1)*(2*nseq-kp)/2+lp-kp;
			mat[index] += length[1]+length[2]+length[4];
			}
		}
	fsd = fopen("sscore","w");
	printf("\n\nSummary Table:\n\n");
	ipr=0;
	iprr = 1;
	printf("     ");
	for(i=0; i<nseq; ++i)
		printf("%s ",code[i]);
	printf("\n");
	printf("%s",code[0]);
	for(i=0; i<nterm+1; ++i)
		{
		ipr++;
		printf(" %4.0f",mat[i]/item);
		if(i==0)
			;
		else
			fprintf(fsd," %4.0f",mat[i]/item);
		if(ipr==nseq)
			{
			ipr = 0;
			fprintf(fsd,"\n");
			printf("\n%s",code[iprr]);
			for(j=0; j<iprr; ++j)
				{
				ipr++;
				printf(" ....");
				}
			printf("    0");
			ipr++;
			iprr++;
			}
		}
	fprintf(fsd,"\n");
	printf("\n");
	fclose(fsd);
	brchord();
}

float score[NSEQ][NSEQ],nmatrix[NSEQ][NSEQ];
char code[NSEQ][200];

brchord()
{
	FILE *fd, *fopen();
	float sc;
	int i,j,item,newnseq,test;

	fd = fopen("sscore","r");

	newnseq = nseq;

	for(i=0; i<nseq-1; ++i)
		{
		if(i<26)
			code[i][0] = 65+i;
		else if(i<52)
			code[i][0] = 71+i;
		else
			exit(printf(stderr,"Too many sequences!!\n"));
		code[i][1] = '\0';
		for(j=i+1; j<nseq; ++j)
			{
			fscanf(fd,"%f",&sc);
			score[i][j] = score[j][i] = sc;
			}
		}
	if(i<26)
		code[i][0] = 65+i;
	else
		code[i][0] = 71+i;
	code[i][1] = '\0';

	for(i=0; i<nseq; ++i)
		for(j=0; j<nseq; ++j)
			if(i==j)
				score[i][j] = 0;
	fclose(fd);

	for(i=0; i<nseq; ++i)
		for(j=i; j<nseq; ++j)
				nmatrix[i][j] = nmatrix[j][i] = score[i][j];
	
	printf("\n\nSequential Branching Order:\n");
	while(newnseq>2)
		{
		grouping(newnseq);
		newnseq--;
		}
	strcat(code[1],code[0]);
	printf("\n\nOverall Branching Order:\n");
	printf("%s\n",code[1]);
}

grouping(newnseq)
int newnseq;
{
	int imin,jmin,i,j,k;
	int ntimes,j1,j2,row,col;
	float min,sum;
	char c1,c2;

	imin = jmin = 0;
	min = 100000;
	for(i=0; i<newnseq-1; ++i)
		for(j=i+1; j<newnseq; ++j)
			if(nmatrix[i][j] < min)
				{
				min = nmatrix[i][j];
				imin = i;
				jmin = j;
				}

	row = col = 0;
	for(i=0; i<newnseq; ++i)
		{
		if(i==imin || i==jmin)
			;
		else
			{
			for(j=0; j<newnseq; ++j)
				nmatrix[row][j] = nmatrix[i][j];
			row++;
			}
		}

	for(i=0; i<newnseq; ++i)
		{
		if(i==imin || i==jmin)
			;
		else
			{
			for(j=0; j<newnseq; ++j)
				nmatrix[j][col] = nmatrix[j][i];
			col++;
			}
		}
	
	if(strlen(code[imin])==1 && strlen(code[jmin])==1)
		{
		strcpy(code[newnseq],"(");
		strcat(code[newnseq],code[imin]);
		strcat(code[newnseq],code[jmin]);
		strcat(code[newnseq],")");
		}
	else
		{
		strcpy(code[newnseq],"(");
		strcat(code[newnseq],code[jmin]);
		strcat(code[newnseq],code[imin]);
		strcat(code[newnseq],")");
		}
	printf("%s\n",code[newnseq]);
	j=0;
	for(i=0; i<newnseq; ++i)
		{
		if(i==imin || i==jmin)
			;
		else
			{
			strcpy(code[j],code[i]);
			j++;
			if(j>=newnseq-2)	
				break;
			}
		}
	strcpy(code[newnseq-2],code[newnseq]);

	nmatrix[newnseq-2][newnseq-2] = 0.0;
	k=0;
	for(i=0; i<newnseq-2; ++i)
		{
		sum=0.0;
		j1=ntimes=0;
		while((c1=code[newnseq-2][j1])!='\0')
			{
			if((c1>='A' && c1<='Z') || (c1>='a' && c1<='z'))
				{
				j2=0;
				while((c2=code[i][j2])!='\0')
					{
					if((c2>='A' && c2<='Z') ||
					   (c2>='a' && c2<='z'))
						{
						if((c1>='A' && c1<='Z') &&
						   (c2>='A' && c2<='Z'))
							sum += score[c2-'A']
								    [c1-'A'];
						else if((c1>='A' && c1<='Z') &&
						        (c2>='a' && c2<='z'))
							sum += score[c2-'a'+26]
								    [c1-'A'];
						else if((c1>='a' && c1<='z') &&
						        (c2>='A' && c2<='Z'))
							sum += score[c2-'A']
								    [c1-'a'+26];
						else if((c1>='a' && c1<='z') &&
						        (c2>='a' && c2<='z'))
							sum += score[c2-'a'+26]
								    [c1-'a'+26];
						ntimes++;
						}
					j2++;
					}
				}
			j1++;
			}
		nmatrix[newnseq-2][k] = nmatrix[k][newnseq-2]
				      = sum/(float)ntimes;
		k++;
		}	
}