/*
 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>
int NumStrings;
char InputString[50][60];

main(ac,av)
int ac;
char *av[];
{
	FILE *B1andD, *B2andD, *fopen();
	int c,i,length;
	char line[400],paren[200];

	if(ac<2)
		{
		fprintf(stderr,"\nPurpose: Determine Branch Lengths");
		fprintf(stderr," Based on Distance Matrix.\n");
		fprintf(stderr,"Usage: ");
		fprintf(stderr,"%s rf [>wf]\n",av[0]);
		fprintf(stderr,"\trf: `albr' (an output file from tree)\n");
		fprintf(stderr,"\twf: Branch Lengths\n");
		fprintf(stderr,"** Use `tomac' for tree drawing. **\n\n");
		exit(0);
		}

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

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

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

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

	paren[0] = '(';
	for(i=1; i<3; ++i)
		paren[i] = start[0][i-1];
	paren[i++] = ')';
	paren[i] = '\0';

	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[200][200],res1[1],res2[1];
char parameter[200], FoundLetters[200];
int ParameterIndex, CurrentBranch, LastBranch;
int bnum[200],or[60],nseq;

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

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

	fd1=fopen("zecret","r");
	fscanf(fd1,"%s",parameter);

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

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

	fd2=fopen("cecret","w");
	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");
			}
		}

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

	LeftBranch[0] = parameter[ParameterIndex++];
	LeftBranch[1] = '\0';
	if(LeftBranch[0]=='(')
		GetBranch(LeftBranch);
	strcpy(branch[CurrentBranch],LeftBranch);
	CurrentBranch++;
	LastBranch++;

	RightBranch[0] = parameter[ParameterIndex++];
	RightBranch[1] = '\0';
	if(RightBranch[0]=='(')
		GetBranch(RightBranch);
	strcpy(branch[CurrentBranch],RightBranch);
	CurrentBranch++;
	LastBranch++;

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

#define NSEQ 52
#define tnsq 2*NSEQ-3
#define nsq NSEQ*(NSEQ-1)/2
#define tpnsq 2*(2*NSEQ-3)

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

brchlen()
{
	FILE *fd1, *fd2, *fd3;
	int i, j, k, kk, nrow, ncol, kont, kmax, pt[30];
	char c, orderstring[300];	
	double A[nsq][tnsq], cum, R[tnsq], Y[nsq], temp;
	double sum, dum, x;

	fd1 = fopen("cecret","r");
	fd2 = fopen("zecret","r");
	fd3 = fopen("tomac","w");
	fscanf(fd2,"%s",orderstring);
	fprintf(fd3,"%s\n",orderstring);

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

	for(i=0; i<nrow; ++i) 
		{
		k=0;
		while((c=fgetc(fd1))!='\n')
			{
			ungetc(c,fd1);
			fscanf(fd1,"%d",&pt[k]);
			k++;
			}
		sort(pt,k);
		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(fd1,"%lf",&x);
		Y[i]=x;
		}
	fclose(fd1);

   	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];
		R[j]=cum;
		}
	printf("\n");

	cum=0;
	printf("Branch lengths are (r1,r2,r3,r4,r5):\n");
	printf("                   (r6,r7,........):\n");
	for(j=0; j<ncol; )
		{
		cum += R[j];
		printf("%8.2f ",R[j]); 
		fprintf(fd3,"%8.2f ",R[j++]);
		if(j%5==0)
			{
			printf("\n");
			fprintf(fd3,"\n");
			}
		}
	printf("\n");
	fprintf(fd3,"\n");

	sum=0;
	kont=0;
	kmax=nseq-1;
	for(i=0; i<nrow; ++i)
		{
		cum=0;
		for(j=0; j<ncol; ++j)
			cum += A[i][j]*R[j];
		dum = fabs(cum - Y[i])*100/Y[i];
		sum += dum*dum;
		}
	printf("\n");

	printf("%% s.d. (observed vs calculated) = %7.2f\n",sqrt(sum/(nrow-1)));
}

inverse(ns,tns)
int ns, tns;
{
	FILE *fd;
	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;
				}
}