/*
By B.Y. Wu, Feb, 2001.

This program roots a tree to minimize its ultrametric size.
The time complexity of the algorithm is O(n^2) and is optimal in time.
This is an academic program and is freely available at
http://www.personal.stu.edu.tw/bangye/mutroot.htm


The input file should be in the following format:

<tree in phylip format, ending with ";">
<The upper triangle of the distance matrix. the distance are
nonnegative floating point numbers separated by white space
(space or carrage return). for example:

m(1,2) m(1,3) ....m(1,n)
m(2,3) m(2,4) ....,m(2,n)
....
m(n-1,n)

the labels of the input tree should be between 0 and n.
The tree should be binary.
num<MAX_N is the number of species.

The output is in phylip format */

/**********/
#include <stdio.h>
#include <io.h>
#include <time.h>
#include <stdlib.h>
#include <float.h>
/* MAX_N is the maximum number of species */
#define MAX_N 80
#define INF 999999.9

struct UTNODE { /* node of ultrametric tree */
		  short label;
		  struct UTNODE *parent,*lson,*rson;
		  float height;
		  short ref;
		  } ;
/* prototype of subroutines */
void badin(void);
struct UTNODE *read_tree(void);
struct UTNODE *getnode(void);
struct UTNODE *findnode(struct UTNODE *,short);
void getinput(void);
struct UTNODE *reroot(struct UTNODE *,struct UTNODE *);
void max_dist(void);
void outree(struct UTNODE *);
void calc_h(struct UTNODE *);
void calc_hi(struct UTNODE *);
float treesize(struct UTNODE *);
float tot_hi(struct UTNODE *);

FILE *inf,*outf;
struct UTNODE *intree,*p1p,*p2p;
/* p1 and p2 are the nodes of longest input distance*/
int n_spe,p1,p2;
float dist[MAX_N][MAX_N],mat[MAX_N][MAX_N],max_d;


main(int argc, char *argv[])
{
int i,num,j;
float cost[MAX_N],cost2[MAX_N],c2;
struct UTNODE *u,*v;
if (argc<2) {
	printf("use root inputfile outputfile\n");
	exit(0);}
if ((inf=fopen(argv[1],"r"))==NULL)
	 {printf("cannot open input file");exit(0);}
getinput();
fclose(inf);
p1p=findnode(intree,p1);
p2p=findnode(intree,p2);
intree=reroot(p1p,p1p->parent);
/* the tree now is rooted at edge incident to the node p1 */
calc_hi(intree);
i=0;  u=p2p;v=u->parent;
cost[i++]=0.0;
while (v->parent!=NULL) {
	c2=(v->lson==u)? tot_hi(v->rson) : tot_hi(v->lson);
	cost[i]=cost[i-1]+c2+v->height;
	u=v;v=v->parent;i++;}

intree=reroot(p2p,p2p->parent);
calc_hi(intree);
u=p1p;v=u->parent;
num=i;
cost2[i-1]=0.0;
i=i-2;
while (v->parent!=NULL) {
	c2=(v->lson==u)? tot_hi(v->rson) : tot_hi(v->lson);
	cost2[i]=cost2[i+1]+c2+v->height;
	u=v;v=v->parent;i--;}
for (i=0;i<num;i++) cost[i]+=cost2[i];
c2=INF;
for (i=0;i<num;i++)
	if (cost[i]<c2) {
		c2=cost[i];j=i;}
u=p1p;v=u->parent;
i=j;
while (i<num-1) {u=v;v=v->parent;i++;}
intree=reroot(u,v);
calc_hi(intree);

if ((outf=fopen(argv[2],"w"))==NULL)
	{printf("cannot open file");exit(0);}
outree(intree);
fprintf(outf,"\nThe minimum tree size is %f\n",c2+max_d);
fclose(outf);
}

struct UTNODE *findnode(struct UTNODE *tree,short ch)
/* find the node with specified label */
{struct UTNODE *t;
if (tree==NULL) return(NULL);
if (tree->label==ch) return(tree);
else if ((t=findnode(tree->lson,ch))!=NULL)
	return(t);
else return(findnode(tree->rson,ch));

}

float treesize(struct UTNODE *root)
{/* the tree size is the total hight of all node plus the
height of root */
return(root->height+tot_hi(root));
}

float tot_hi(struct UTNODE *root)
{
if (root->lson==NULL)
	return(root->height);
else
	return(root->height+tot_hi(root->lson)+tot_hi(root->rson));
}

void calc_hi(struct UTNODE *root)
{/* recursively call calc_h to find the total height of all node */
int i,j;
/* reserve the distances */
for (i=0;i<n_spe-1;i++)
	for (j=i;j<n_spe;j++)
		 mat[i][j]=mat[j][i]=dist[i][j];
calc_h(root);
}

void calc_h(struct UTNODE *root)
/* calculate the height of nodes in the tree rooted at root */
{
int i;
float d1,d2;
if (root->lson!=NULL) {
	calc_h(root->lson);
	calc_h(root->rson);
	root->height=(root->lson->height+root->rson->height
		+mat[root->lson->ref][root->rson->ref])/2;
	if (root->height<root->lson->height)
		root->height=root->lson->height;
	if (root->height<root->rson->height)
		root->height=root->rson->height;
	root->ref=root->lson->ref;
	for (i=0;i<n_spe;i++) {
		d1=mat[i][root->lson->ref]-root->height+root->lson->height;
		d2=mat[i][root->rson->ref]-root->height+root->rson->height;
		d1=(d1>d2)? d1 : d2;
		d1=(d1>0)? d1 : 0.0;
		mat[i][root->ref]=mat[root->ref][i]=d1;}
	}
else {
	root->ref=root->label;
	root->height=0.0;}
}

struct UTNODE *reroot(struct UTNODE *r1, struct UTNODE *r2)
{/* reroot a tree at an edge. */
struct UTNODE *t1,*t2,*t3,*t4;

t3= (r1->parent==r2)? r1:r2;
t2=t3->parent;
/* insert a node between r1,r2 */
t4=t1=getnode();
t1->lson=t3;t1->rson=t2;t1->parent=NULL;
t3->parent=t1;
if (t2->lson==t3) t2->lson=t1;
else t2->rson=t1;
/* reroot the tree at the new node */
while (t2!=NULL) {
	if (t1->lson==NULL) t1->lson=t2;
	else t1->rson=t2;
	t3=t2->parent;
	t2->parent=t1;
	if (t2->lson==t1) t2->lson=NULL;
	else t2->rson=NULL;
	t1=t2;t2=t3;
	}
/* t1 is the old root */
t2=(t1->lson!=NULL) ? t1->lson : t1->rson;
if (t1->parent->lson==t1) t1->parent->lson=t2;
else t1->parent->rson=t2;
t2->parent=t1->parent;
return(t4);
}

struct UTNODE *getnode()
{struct UTNODE *t;
if ((t=malloc(sizeof(struct UTNODE)))==NULL) {
	printf("insufficient memory\n");exit(-1);}
return(t);
}

struct UTNODE *read_tree()
{/* read a binary tree in phylip format. branch lengths will be ignored */
int ch;
struct UTNODE *stack[2*MAX_N];
int i=0,newnode=MAX_N,j;
struct UTNODE *t1;
ch=getc(inf);
while ((ch!='(')&&(!feof(inf)))
	 ch=getc(inf);
if (ch!='(') badin();
while ((ch!=';')&& !feof(inf)) {
	switch (ch) {
	case ',': break;
	case '(': stack[++i]=NULL; break;
	case ')': /* pop two node from stack and merge the two subtrees */
				 if ((i<3)||(stack[i]==NULL)||(stack[i-1]==NULL)||
					  (stack[i-2]!=NULL))
					  badin();
				 t1=getnode();
				 t1->label=newnode++;
				 t1->parent=NULL;
				 t1->rson=stack[i--];
				 t1->lson=stack[i--];
				 stack[i]=t1;
				 t1->rson->parent=t1->lson->parent=t1;
				 t1->height=0;
				 break;
	default : t1=getnode();
				 t1->parent=t1->lson=t1->rson=NULL;
				 t1->label=ch;t1->height=0;
				 stack[++i]=t1;
	}
	/* read next token from file*/
	j=getc(inf);
	while (((j==0x20)||(j==0x0d)||(j==0x0a))&&(!feof(inf)))
		 j=getc(inf);
	if ((j=='(') || (j==')') ||(j==';')||(j==',')) ch=j;
	else {
		if ((j<48)||(j>57)) badin();
		ch=0;
		while ((j>=48)&&(j<58)) {ch=ch*10+j-48;j=getc(inf);}
		if (j==':') while ((j!=',')&&(j!=')')&&(!feof(inf))) j=getc(inf);
		if (j==')') ungetc(j,inf);
		if (ch>=n_spe) badin();
		}
}
if (i!=1) badin();
return(stack[1]);
}

void badin()
{printf("bad input format\n");fclose(inf);exit(-1);}

/*outree output a tree file in phylip format */
void outree(struct UTNODE *root)
{
if ((root->lson==NULL)&&(root->rson==NULL))
	{fprintf(outf,"%u", root->label);
	fprintf(outf,":");
	fprintf(outf,"%.3f",root->parent->height-root->height);
	 }
else {
fprintf(outf,"(");
if (root->lson!=NULL) outree(root->lson);
fprintf(outf,",\n");
if (root->rson!=NULL) outree(root->rson);
fprintf(outf,")");
if (root->parent!=NULL) {
	fprintf(outf,":");
	fprintf(outf,"%.3f",root->parent->height-root->height);}
else fprintf(outf,";");
}
}

void getinput(void)
/* get input data from input file; */
{
int i,j;
float value;
fscanf(inf,"%u",&i);
n_spe=i;
if ((n_spe<2)||(n_spe>MAX_N)) badin();
intree=read_tree();
max_d=0.0;p1=1;p2=2;
for (i=0;i<n_spe;i++)
	 for (j=i+1;j<n_spe;j++) {
		if (fscanf(inf,"%f",&value)==NULL)
			badin();
		dist[i][j]=value;
		if (value>max_d) {
			max_d=value;
			p1=i;p2=j;
			}
		}
for (i=0;i<n_spe;i++) {
	for (j=0;j<i;j++)
		dist[i][j]=dist[j][i];
	dist[i][i]=0.0;}

}