/*
 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,ndiff;
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: Find Best Branch Order With");
		fprintf(stderr," No Negative Branch Lengths.\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=ndiff=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);
	for(i=0; i<60; i++)
		{
		branchlength();
		setree();
		flip();
		}
}

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;

   /*Read branching order*/
	fd1=fopen("zecret","r");
	fscanf(fd1,"%s",parameter);
	printf("%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++;
		}

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

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

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

#define NSEQ 55
#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];
	int negno;
	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);
	fclose(fd2);
	fprintf(fd3,"%s\n",orderstring);

   /*Enter A_matrix*/
	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);   /*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;
			}
		}

   /*Enter Y_matrix*/
	for(i=0; i<nrow; ++i)
		{
		fscanf(fd1,"%lf",&x);
		Y[i]=x;
		}
	fclose(fd1);

   /*Take the product of transpose_of_A and A*/
   	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;
			}

   /*Take the inverse of a_matrix*/
	inverse(ncol, 2*ncol);

   /*Take the product matrix of A and ia*/
	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;
			}

   /*Product matrix of transpose of Y and a*/
	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");

   /*Solution*/
	negno=0;
	for(j=0; j<ncol; ++j)
		{
		if(R[j]<0)
			{
			negno=1;
			break;
			}
		}
	

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

   /*Recompute difference matrix*/
	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;
		}

   /*% standard deviation*/
	printf("\t%% s.d. (obs. vs calc.) = %7.2f\n",sqrt(sum/(nrow-1)));
}

inverse(ns,tns)  /*Take the inverse of a square matrix*/
int ns, tns;
{
	FILE *fd;
	int i, ii, ip, j, jj, k, l, n, np;
	double cum, den, dum, m, s;

	n=ns;
   /*Attach an identity matrix to a-matrix*/
	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;
			}
		}

   /*Use Gaussian elimination method to convert the a_matrix into an 
     upper triangular matrix*/

	for(k=0; k<n-1; ++k)  /*Begin: Gauss elimination procedure*/
		{             /*k is column, stop at the matrix*/
		l=k;   /*Begin: Minimize roundoff error*/
		for(i=k+1; i<n; i++)  /*i is row*/
			{
			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;
				}
			}   /*End: Minimize roundoff error*/

		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];
			}
		}   /*End: Gauss elimination procedure*/	


   /*Back elimination*/
	np=n-1;   /*Begin: Back elimination procedure*/
	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;
				}
}

setree()
{
	FILE *fd, *fd1, *fopen();
	int i,j,k;
	char branch_order[300], branch_length[1300][8];
	char bostring[500];
	char c,c1,c2,c3,sign;

	fd = fopen("tomac","r");
	fscanf(fd,"%s",branch_order);

	i=0;
	while(fscanf(fd,"%s",branch_length[i++]) != EOF)
		;
	fclose(fd);

	fd1 = fopen("duup1","w");
	i=k=0;
	while(branch_order[i] != '\0')
		{
		c1 = branch_order[i];
		c2 = branch_order[i+1];
		c3 = branch_order[i+2];
		if(c3 == '\0' && ((c2>='A' && c2<='Z')||(c2>='a'&&c2<='z')))
			{
			if(branch_length[k][0]=='-')
				fprintf(fd1,")*%c",c2);
			else
				fprintf(fd1,")%c",c2);
			goto out;
			}
		else
			{
	/* (  */	if(c1=='(')
				fprintf(fd1,"%c",c1);
	/* XY */	else if((c1>='A' && c1<='Z') && (c2>='A' && c2<='Z'))
				{
				sign = branch_length[k][0];
				if(sign=='-')
					fprintf(fd1,"%c*",c1);
				else
					fprintf(fd1,"%c",c1);
				k++;
				}
	/* xy */	else if((c1>='a' && c1<='z') && (c2>='a' && c2<='z'))
				{
				sign = branch_length[k][0];
				if(sign=='-')
					fprintf(fd1,"%c*",c1);
				else
					fprintf(fd1,"%c",c1);
				k++;
				}
	/* Xy */	else if((c1>='A' && c1<='Z') && (c2>='a' && c2<='z'))
				{
				sign = branch_length[k][0];
				if(sign=='-')
					fprintf(fd1,"%c*",c1);
				else
					fprintf(fd1,"%c",c1);
				k++;
				}
	/* xY */	else if((c1>='a' && c1<='z') && (c2>='A' && c2<='Z'))
				{
				sign = branch_length[k][0];
				if(sign=='-')
					fprintf(fd1,"%c*",c1);
				else
					fprintf(fd1,"%c",c1);
				k++;
				}
	/* X) */	else if((c1>='A' && c1<='Z') && c2==')')
				{
				sign = branch_length[k][0];
				if(sign=='-')
					fprintf(fd1,"%c*",c1);
				else
					fprintf(fd1,"%c",c1);
				k++;
				}
	/* x) */	else if((c1>='a' && c1<='z') && c2==')')
				{
				sign = branch_length[k][0];
				if(sign=='-')
					fprintf(fd1,"%c*",c1);
				else
					fprintf(fd1,"%c",c1);
				k++;
				}
	/* )Y */	else if(c1==')' && (c2>='A' && c2<='Z'))
				{
				sign = branch_length[k][0];
				if(sign=='-')
					fprintf(fd1,"%c*",c1);
				else
					fprintf(fd1,"%c",c1);
				k++;
				}
	/* Y( */	else if((c1>='A' && c1<='Z') && c2=='(')
				{
				sign = branch_length[k][0];
				if(sign=='-')
					fprintf(fd1,"%c*",c1);
				else
					fprintf(fd1,"%c",c1);
				k++;
				}
	/* y( */	else if((c1>='a' && c1<='z') && c2=='(')
				{
				sign = branch_length[k][0];
				if(sign=='-')
					fprintf(fd1,"%c*",c1);
				else
					fprintf(fd1,"%c",c1);
				k++;
				}
	/* )z */	else if(c1==')' && (c2>='a' && c2<='z'))
				{
				sign = branch_length[k][0];
				if(sign=='-')
					fprintf(fd1,"%c*",c1);
				else
					fprintf(fd1,"%c",c1);
				k++;
				}
	/* )) */	else if(c1==')' && c2==')')
				{
				sign = branch_length[k][0];
				if(sign=='-')
					fprintf(fd1,"%c*",c1);
				else
					fprintf(fd1,"%c",c1);
				k++;
				}
	/* )( */	else if(c1==')' && c2=='(')
				{
				sign = branch_length[k][0];
				if(sign=='-')
					fprintf(fd1,"%c*",c1);
				else
					fprintf(fd1,"%c",c1);
				k++;
				}
	/* )\0 */	else if(c1==')' && c2=='\0')
				fprintf(fd1,")");
	/* X\0 */	else if((c1>='A' && c1<='Z') && c2=='\0')
				fprintf(fd1,"%c",c1);
	/* x\0 */	else if((c1>='a' && c1<='z') && c2=='\0')
				fprintf(fd1,"%c",c1);
			else
				exit(printf("Problem in setree.\n"));
			}
		i++;
		}
	out:
	fprintf(fd1,"\n");
	fclose(fd1);
}

flip()
{
	FILE *fd, *fdz, *fdd, *fopen();
	int i,j,star_pos,left,ibrkt,isw,mleft,mmleft,mright;
	int li,ri,nstar,old_star_pos,lcluster,rcluster;
	int olddiff,newdiff,k,l,cleft,cright;
	char c,swstring[400],newstring[400];
	char left_entry[100],right_entry[100];
	char line[500],sub_left[50],sub_right[50];

	fd=fopen("duup1","r");
	fgets(swstring, sizeof swstring, fd);
	fclose(fd);

	nstar=0;
	for(i=0; swstring[i]!='\0'; ++i)
		if(swstring[i]=='*')
			{
			c = swstring[i-1];
			if(c>='A'&&c<='Z' || c>='a'&&c<='z')
				{
				printf("Negative terminal branch lengths");
				printf(" detected.  Job aborted.\n");
				printf("You may have two very similar ");
				printf("sequences.\n");
				exit(0);
				}
			else
				nstar++;
			}
	if(nstar==0)
		{
		remove("duup1");
		exit(0);
		}

	star_pos = 0;
	while(swstring[star_pos]!='*')
		star_pos++;		/*Locate * position*/

	newdiff = abs(star_pos - old_star_pos);
	if(newdiff==olddiff)
		ndiff++;
	else
		ndiff=0;

	/*Encountered deep node*/
	if(swstring[star_pos-1]==')' && swstring[star_pos+1]==')')
		{
		/*Move left, find (*/
		i=star_pos-2;
		ibrkt=1;
		while(ibrkt!=0)
			{
			if(swstring[i]==')')
				ibrkt++;
			else if(swstring[i]=='(')
				ibrkt--;
			i--;
			}
		left = i; /*1st outside the (*/
		/*End: Move left, find (*/

		/*Move left to corresponding (*/
		mleft=left;
		ibrkt=0;
		for(;;)
			{
			if(swstring[mleft]==')')
				ibrkt++;
			else if(swstring[mleft]=='(')
				ibrkt--;
			if(ibrkt==0)
				break;
			mleft--;
			}
		/*End: Move left to corresponding (*/

		i=0;
		for(isw=0; isw<mleft+1; ++isw)
			newstring[i++] = swstring[isw];
		newstring[i++] = '('; /*Insert a (*/

		for(isw=mleft+1; isw<left+1; ++isw)
			newstring[i++] = swstring[isw];

		isw++;  /*skip one in swstring*/
		while(isw < star_pos-1)
			{
			/*No cluster found*/
			if(swstring[isw]!='(')
				newstring[i++] = swstring[isw++];

			else		/*Cluster found*/
				{
				ibrkt=0;
				for(;;)
					{
					if(swstring[isw]=='(')
						{
						newstring[i++]=swstring[isw++];
						ibrkt++;
						}
					else if(swstring[isw]==')')
						{
						newstring[i++]=swstring[isw++];
						ibrkt--;
						}
					else
						newstring[i++]=swstring[isw++];
					if(ibrkt==0)
						break;
					}
				}
				/*End: Cluster found*/
			newstring[i++] = ')';
			}
			/*No cluster found*/
		for(j=star_pos+2; swstring[j]!='\0'; ++j)
			newstring[i++]=swstring[j];
		newstring[i] = '\0';
		}
	/*End: Encountered deep node*/

	else if(ndiff>5 && star_pos>old_star_pos)
		{
		/*Encounter type 2 duplicated shallow nodes*/
		lcluster=li=0;
		mleft=star_pos-2;
		mmleft=star_pos-1;
		if(swstring[mleft]!=')') /*No cluster*/
			{
			lcluster++;
			}
		else
			{
			ibrkt=0;
			for(;;)
				{
				if(swstring[mleft]==')')
					ibrkt++;
				else if(swstring[mleft]=='(')
					ibrkt--;
				if(ibrkt==0)
					break;
				mleft--;
				}
			ibrkt=0;
			for(;;)
				{
				if(swstring[mmleft]==')')
					ibrkt++;
				else if(swstring[mmleft]=='(')
					ibrkt--;
				if(ibrkt==0)
					break;
				mmleft--;
				}
			}
	
		rcluster=ri=0;
		mright=star_pos+1;
		if(swstring[star_pos+1]!='(') /*No cluster*/
			{
			right_entry[ri++] = swstring[star_pos+1];
			rcluster++;
			}
		else  /*Cluster found*/
			{
			ibrkt=0;
			for(;;)
				{
				if(swstring[mright]=='(')
					ibrkt++;
				else if(swstring[mright]==')')
					ibrkt--;
				if(ibrkt==0)
					break;
				mright++;
				}
			for(i=star_pos+1; i<=mright; ++i)
				right_entry[ri++] = swstring[i];
			}
		right_entry[ri] = '\0';

		/* cluster on the left and no cluster on the right */
		if(lcluster==0 && rcluster!=0)
			{
			j=0;
			for(i=0; i<mmleft; ++i)
				newstring[j++] = swstring[i];
			newstring[j] = '\0';
			for(i=mmleft+1; i<mleft; ++i)
				newstring[j++] = swstring[i];
			newstring[j++] = '(';
			newstring[j] = '\0';
			for(i=mleft; i<star_pos-1; ++i)
				newstring[j++] = swstring[i];
			newstring[j] = '\0';
			strcat(newstring,right_entry);
			j=strlen(newstring);
			for(i=star_pos-1; swstring[i]!='*'; ++i)
				newstring[j++]=swstring[i];
			newstring[j] = '\0';
			for(i=mright+1; swstring[i]!='\n'; ++i)
				newstring[j++]=swstring[i];
			newstring[j] = '\0';
			}

		/* no cluster on the left and cluster on the right */
		else if(lcluster!=0 && rcluster==0)
			{
			while(swstring[mleft]!='(')
				mleft--;
			j=0;
			for(i=0; i<mleft; ++i)
				newstring[j++] = swstring[i];
			newstring[j++] = '(';
			newstring[j] = '\0';
			for(i=mleft; i<star_pos-1; ++i)
				newstring[j++] = swstring[i];
			newstring[j] = '\0';

			k=1;
			l=0;
			if(right_entry[k]!='(')
				sub_left[l++] = right_entry[k++];
			else
				{
				ibrkt=0;
				for(;;)
					{
					if(right_entry[k]=='(')
						{
						sub_left[l++]=right_entry[k];
						ibrkt++;
						}
					else if(right_entry[k]==')')
						{
						sub_left[l++]=right_entry[k];
						ibrkt--;
						}
					if(ibrkt==0)
						break;
					k++;
					}
				}
			sub_left[l] = '\0';

			l=0;
			if(right_entry[k]!='(')
				sub_right[l++] = right_entry[k++];
			else
				{
				ibrkt=0;
				for(;;)
					{
					if(right_entry[k]=='(')
						{
						sub_right[l++]=right_entry[k];
						ibrkt++;
						}
					else if(right_entry[k]==')')
						{
						sub_right[l++]=right_entry[k];
						ibrkt--;
						}
					if(ibrkt==0)
						break;
					k++;
					}
				}
			sub_right[l] = '\0';

			newstring[j++]=')';
			newstring[j] = '\0';
			strcat(newstring,sub_right);
			j=strlen(newstring);
			newstring[j] = '\0';

			newstring[j++] = ')';
			newstring[j] = '\0';
			strcat(newstring,sub_left);
			j=strlen(newstring);
			newstring[j] = '\0';
			for(i=mright+1; swstring[i]!='\0'; ++i)
				newstring[j++] = swstring[i];
			newstring[j]='\0';
			}
		else
			{
			li=0;
			mleft=star_pos-2;
			mv_pointer_left(&mleft,swstring);
			for(i=mleft; i<star_pos-1; ++i)
				left_entry[li++] = swstring[i];
			left_entry[li++] = '\0';

			/*mmleft is the position of the leftmost (*/
			mmleft=star_pos-1;
			mv_pointer_left(&mmleft,swstring);

			/*set up the new string*/
			j=0;
			for(i=0; i<mmleft; ++i)
				newstring[j++] = swstring[i];
			for(i=mmleft+1; i<mleft; ++i)
				newstring[j++] = swstring[i];
			newstring[j] = '\0';
			j=strlen(newstring);
			newstring[j++] = '(';
			newstring[j] = '\0';
			strcat(newstring,left_entry);
			strcat(newstring,right_entry);
			j=strlen(newstring);
			newstring[j++] = ')';
			for(i=mright+1; swstring[i]!='\0'; ++i)
				newstring[j++] = swstring[i];
			newstring[j] = '\0';
			cleft = cright = 0;
			}
		ndiff=0;
		/*End: Encounter type 2 duplicated shallow nodes*/
		}

	else if(old_star_pos == star_pos && cright==0 && cleft==0)
		{
		/*Encounter type 1 duplicated shallow nodes*/
		li=0;
		mleft=star_pos-3;
		if(swstring[mleft]!=')') /*No cluster*/
			{
			left_entry[li++] = swstring[mleft];
			}
		else  /*Cluster found*/
			{
			ibrkt=0;
			for(;;)
				{
				if(swstring[mleft]==')')
					ibrkt++;
				else if(swstring[mleft]=='(')
					ibrkt--;
				if(ibrkt==0)
					break;
				mleft--;
				}
			for(i=mleft; i<star_pos-2; ++i)
				left_entry[li++] = swstring[i];
			}
		left_entry[li] = '\0';
	
		ri=0;
		mright=star_pos+1;
		if(swstring[star_pos+1]!='(') /*No cluster*/
			{
			right_entry[ri++] = swstring[star_pos+1];
			}
		else  /*Cluster found*/
			{
			ibrkt=0;
			for(;;)
				{
				if(swstring[mright]=='(')
					ibrkt++;
				else if(swstring[mright]==')')
					ibrkt--;
				if(ibrkt==0)
					break;
				mright++;
				}
			for(i=star_pos+1; i<=mright; ++i)
				right_entry[ri++] = swstring[i];
			}
		right_entry[ri] = '\0';

		j=0;
		for(i=0; i<mleft; ++i)
			newstring[j++] = swstring[i];
		newstring[j] = '\0';
		strcat(newstring,right_entry);
		j=strlen(newstring);
		for(i=star_pos-2; swstring[i]!='*'; ++i)
			newstring[j++]=swstring[i];
		newstring[j] = '\0';
		strcat(newstring,left_entry);
		j=strlen(newstring);
		for(i=mright+1; swstring[i]!='\0'; ++i)
			newstring[j++] = swstring[i];
		newstring[j] = '\0';
		/*End: Encounter type 2 duplicated shallow nodes*/
		}

	else if(old_star_pos == star_pos)
		{
	/*Encounter type 1 duplicated shallow nodes*/
		li=0;
		mleft=star_pos-2;
		if(swstring[mleft]!=')') /*No cluster*/
			{
			left_entry[li++] = swstring[mleft];
			}
		else  /*Cluster found*/
			{
			ibrkt=0;
			for(;;)
				{
				if(swstring[mleft]==')')
					ibrkt++;
				else if(swstring[mleft]=='(')
					ibrkt--;
				if(ibrkt==0)
					break;
				mleft--;
				}
			for(i=mleft; i<star_pos-1; ++i)
				left_entry[li++] = swstring[i];
			}
		left_entry[li] = '\0';
		
		/*Move further left*/
		ibrkt=0;
		mmleft=star_pos-1;
		for(;;)
			{
			if(swstring[mmleft]==')')
				ibrkt++;
			else if(swstring[mmleft]=='(')
				ibrkt--;
			if(ibrkt==0)
				break;
			mmleft--;
			}
		/*The ( at mmleft must be eliminated*/

		ri=0;
		mright=star_pos+1;
		if(swstring[star_pos+1]!='(') /*No cluster*/
			{
			right_entry[ri++] = swstring[star_pos+1];
			}
		else  /*Cluster found*/
			{
			ibrkt=0;
			for(;;)
				{
				if(swstring[mright]=='(')
					ibrkt++;
				else if(swstring[mright]==')')
					ibrkt--;
				if(ibrkt==0)
					break;
				mright++;
				}
			for(i=star_pos+1; i<=mright; ++i)
				right_entry[ri++] = swstring[i];
			}
		right_entry[ri] = '\0';

		j=0;
		for(i=0; i<mmleft; ++i)
			newstring[j++] = swstring[i];
		for(i=mmleft+1; i<mleft; ++i)
			newstring[j++] = swstring[i];
		newstring[j] = '\0';
		j=strlen(newstring);
		newstring[j++] = '(';
		newstring[j] = '\0';
		strcat(newstring,left_entry);
		strcat(newstring,right_entry);
		j=strlen(newstring);
		newstring[j++] = ')';
		for(i=mright+1; swstring[i]!='\0'; ++i)
			newstring[j++] = swstring[i];
		newstring[j] = '\0';
		cleft = cright = 0;
		/*End: Encounter type 2 duplicated shallow nodes*/
		}
	else
		{
		/*Encounter shallow nodes*/
		li=cleft=cright=0;
		mleft=star_pos-2;
		if(swstring[mleft]!=')') /*No cluster*/
			{
			left_entry[li++] = swstring[mleft];
			}
		else  /*Cluster found*/
			{
			cleft=1;
			mv_pointer_left(&mleft,swstring);
			for(i=mleft; i<star_pos-1; ++i)
				left_entry[li++] = swstring[i];
			}
		left_entry[li] = '\0';
	
		ri=0;
		mright=star_pos+1;
		if(swstring[mright]!='(') /*No cluster*/
			{
			right_entry[ri++] = swstring[mright];
			}
		else  /*Cluster found*/
			{
			cright=1;
			ibrkt=0;
			for(;;)
				{
				if(swstring[mright]=='(')
					ibrkt++;
				else if(swstring[mright]==')')
					ibrkt--;
				if(ibrkt==0)
					break;
				mright++;
				}
			for(i=star_pos+1; i<=mright; ++i)
				right_entry[ri++] = swstring[i];
			}
		right_entry[ri] = '\0';

		j=0;
		for(i=0; i<mleft; ++i)
			newstring[j++] = swstring[i];
		newstring[j] = '\0';
		strcat(newstring,right_entry);
		j=strlen(newstring);
		newstring[j++]=')';
		newstring[j] = '\0';
		strcat(newstring,left_entry);
		j=strlen(newstring);
		for(i=mright+1; swstring[i]!='\n'; ++i)
			newstring[j++] = swstring[i];
		newstring[j] = '\0';
		/*End: Encounter shallow nodes*/
		}

	/*Remove * from newstring*/
	j=0;
	for(i=0; newstring[i]!='\0'; ++i)
		{
		if(newstring[i]>='A' && newstring[i]<='Z')
			newstring[j++] = newstring[i];
		else if(newstring[i]>='a' && newstring[i]<='z')
			newstring[j++] = newstring[i];
		else if(newstring[i]=='(' || newstring[i]==')')
			newstring[j++] = newstring[i];
		else
			;	
		}
	newstring[j] = '\0';

	fdz=fopen("zecret","r");
	fdd=fopen("duup1","w");
	fprintf(fdd,"%s\n",newstring);
	fgets(line, sizeof line, fdz);
	while(fgets(line, sizeof line, fdz)!=NULL)
		{
		fprintf(fdd,"%s",line);
		}
	fclose(fdz);
	fclose(fdd);
	fdz=fopen("zecret","w");
	fdd=fopen("duup1","r");
	while(fgets(line, sizeof line, fdd)!=NULL)
		{
		fprintf(fdz,"%s",line);
		}
	fclose(fdz);
	fclose(fdd);
	old_star_pos = star_pos;
	olddiff = newdiff;
	remove("duup1");
}

checking(s)
char *s;
{
	int i,nc,noq,ncq,nn,ns,sc;

	nc=noq=ncq=nn=ns=sc=0;
	for(i=0; s[i]!='\0'; ++i)
		{
		if(s[i]>='A' && s[i]<='Z')
			nc++;
		else if(s[i]=='(')
			noq++;
		else if(s[i]==')')
			ncq++;
		else if(s[i]=='\n')
			nn++;
		else if(s[i]=='*')
			ns++;
		else
			sc++;
		}
}

mv_pointer_left(mleft,swstring)
int *mleft;
char *swstring;
{
	int ibrkt;

	ibrkt=0;
	for(;;)
		{
		if(swstring[*mleft]==')')
			ibrkt++;
		else if(swstring[*mleft]=='(')
			ibrkt--;
		if(ibrkt==0)
			break;
		(*mleft)--;
		}
	return;
}