#include <stdio.h>
#include <string.h>

#define VERSION 2.11

#define  MEASURE_EUCLDIST  0x01
#define  MEASURE_VECTANGLE 0x02
#define  MEASURE_VECTLEN   0x04
#define  MEASURE_CITYBLOCK 0x08
#define  MERGE_MEAN        0x01
#define  MERGE_MIN         0x02

#define  STARTX            20    /* The minimum x-margin of the figure  */
#define  STRLEN           200    /* The number of points used for the label */
#define  STARTY            80    /* The y-margin for the graph */

#define  DEFAULT_DIST      20    /* Distance (in points) between 2 entries */
#define  DEFAULT_FIGWIDTH 600    /* The width (in points) of the figure */
#define  DEFAULT_FONTSIZE  12    /* The size of all characters in the fig. */
#define  DEFAULT_RESOLUTION 0.5  /* scale resolution of axis */
#define  DEFAULT_FILENAME "noname.fig" /* Filename for fig-file */
#define  DEFAULT_MERGEFLAG MERGE_MIN

#define  VOID              -1
#define  ABS(x)      ((x)<0.0?(-(x)):(x))
#define  NOTVOID(x)  (x>-0.1)
#define  MIN(x,y)    ((x)<(y)?(x):(y))
#define  MEAN(x,y)   (((x)+(y))/2.0)

extern double atof(char *str);
extern double sqrt(double x);
extern int    atoi(char *str);
extern double acos(double x);

int MERGEFLAG;
int debugflag=0;

/* structures. */
  
typedef struct treenode {
      int x, y;             /* the components */
      double distance;      /* the distance between x and y */
      struct treenode *ptr; /* pointer to the next "collapse" */
} TNODE;

typedef struct position {
      int x, y;
} POINT;


/* matrix manipulation routines. 
   note that we only use the upper right corner. 
   */

double get(double **m, int r, int c) /* get wanted distance in matrix */
{
   if (NOTVOID(m[r][c]))
      return m[r][c];
   else return m[c][r];
}

void set(double **m, int r, int c, double value) /* set matrix element */
{
   if (NOTVOID(m[r][c]))
      m[r][c]=value;
   else m[c][r]=value;
}

/* the recursive function for merging distances. */

TNODE *collapse(double **m, int dim)
{
   double min=99999.0; 
   int minr=0, minc=0;
   int r, c;
   TNODE *node;
   
   for (r=0; r<dim-1; r++) {
      for (c=r+1; c<dim; c++) {
	 if (NOTVOID(m[r][c])) {
	   if (m[r][c]<min) {
	     min=m[r][c];
	     minr=r;
	     minc=c;
	   }
	 }
      }
   }

   if (minc==minr) {
      return (TNODE *)NULL;
   }
   
   if (debugflag) {
      printf("-----------------------------------------\n");
      for (r=0; r<dim; r++) {
	 for (c=0; c<dim; c++) {
	    if (NOTVOID(m[r][c]))
	       printf("%+5.3lf",m[r][c]);
	    else
	       printf("      ");
	    if ((r==minr) && (c==minc))
	       printf("* ");
	    else
	       printf("  ");
	 }
	 printf("\n");
      }
   }
		   
   for (r=0; r<dim; r++) {
      if (NOTVOID(get(m,r,minr)) && NOTVOID(get(m,r,minc))) {
	if (MERGEFLAG==MERGE_MIN) {
	  /* update distance: DISTxy=MIN(DISTx,DISTy) */
	  set(m,minr,r,MIN(get(m,r,minr),get(m,r,minc)));
	} else 	if (MERGEFLAG==MERGE_MEAN) {
	  /* update distance: DISTxy=(DISTx+DISTy)/2 */
	  set(m,minr,r,MEAN(get(m,r,minr),get(m,r,minc)));
	} else {
	  fprintf(stderr,"Merge flag not set\n");
	  exit(22);
	}
      }
   }
   
   for (r=0; r<dim-1; r++) {
      for (c=r+1; c<dim; c++) {
	 if (r==minc)
	    m[minc][c]=VOID;
	 else if (c==minc)
	    m[r][minc]=VOID;
      }
   } 


   node=(TNODE *)malloc(sizeof(TNODE));
   node->ptr=collapse(m, dim);
   node->x=minr; node->y=minc;
   node->distance=min;
   return node;
}

/* help function(s) for printing the map. */

int how_many_smaller(int *list, int no, int value)
{
   int i,a=0;
   for (i=0; i<no; i++) {
      if (list[i]<value)
	 a++;
   }
   return a;
}
      
/* calculate the euclidean distance. */

double eucl(double *v1, double *v2, int k)
{
  int i;
  double sum=0.0;

  for (i=0; i<k; i++)
    sum+=(v1[i]-v2[i])*(v1[i]-v2[i]);
  return sqrt(sum);
}

double vectormult(double *v1, double *v2,int k)
{
  int i,j;
  double sum=0.0;

  for (i=0; i<k; i++) {
    sum+=v1[i]*v2[i];
  }
  return sum;
}

double vlength(double *v1, int k)
{
  return (sqrt(vectormult(v1,v1,k)));
}

/* calculate the city-block distance between two vectors */
double cityblock(double *v1, double *v2, int k)
{
  double sum;
  int i;

  sum=0.0;
  for (i=0; i<k; i++)
    sum+=ABS(v1[i]-v2[i]);
  return sum;
}

/* calculate the angle between two vectors. */

double angle(double *v1, double *v2, int k)
{
   return acos(vectormult(v1,v2,k)/(vlength(v1,k)*vlength(v2,k)));
}

/* calculate the absolute difference in length between two vectors. */

double length(double *v1, double *v2, int k)
{
   return (ABS((vlength(v1,k)-vlength(v2,k))));
}

/* help function for reading input correctly. */

char *removeNL(char *tok)
{
  int i;
  for (i=0; tok[i]!='\n' && tok[i]!='\0'; i++) ;
  tok[i]='\0';
  return tok;
}  

/* the main program: 
   (1) read std input, 
   (2) create matrix with euclidean distances, 
   (3) create a sequential list of merged pairs,
   (4) decide the optimal order, and
   (5) print them out.
   */

int main(int argc, char **argv)
{
   int a,i,j,d,k,mode=MEASURE_EUCLDIST;
   double **v, **m;
   double max;
   char buf[10001], str[80], *tok;
   TNODE *head, *ptr;
   char **labels;
   POINT *points;
   POINT *pairs;
   int   *order;
   FILE  *fp;
   char   filename[80];
   
   int DIST=          DEFAULT_DIST;
   int FIGWIDTH=      DEFAULT_FIGWIDTH;
   int FONTSIZE=      DEFAULT_FONTSIZE;
   double RESOLUTION= DEFAULT_RESOLUTION;

   MERGEFLAG=         DEFAULT_MERGEFLAG;

   strcpy(filename,DEFAULT_FILENAME);

   for (a=1; a<argc; a++) {
      if (argv[a][0]=='-') {
	 for (i=1; argv[a][i]!='\0'; i++) {
	    switch (argv[a][i]) {
	       case 'h': 
	       fprintf(stderr,"-Hierarchical-cluster-analysis-utility-------------------------v%3.1f--\nThis program generates a hierarchical cluster map in FIG v2.1 format.\n",VERSION);
	       fprintf(stderr,"It takes as input rows of text on the following format:\n\t<textlabel> <value> <value> <value> ... <value>\nIf no label appears the system will generate one which identifies the\nrow by number.\n\nYou will need the fig2dev utility in the Transfig-package  (available\nfrom export.lcs.mit.edu) to convert the generated file to printer and \nword-processing languages (e.g. postscript, latex etc.). In addition,\nthe xfig-facility  (available from  ftp.cs.cornell.edu) allows you to\nedit the generated file.\nThis program was written by Mikael Boden.\n---------------------------------------------------------------------\n"); 
	       fprintf(stderr,"Usage: %s [-ealp] [-f string] [-s int] [-r float] [-w int] [-d int]\n\t-e use euclidean distance (default)\n\t-c use city-block distance\n\t-a use vector angle\n\t-l use vector length\n\t-m use merge method -m1 DISTxy=min(DISTx,DISTy) (default) -m2 DISTxy=(DISTx+DISTy)/2\n\t-p produce debug-info\n\t-f write to specified file (default \"%s\")\n\t-s specify fontsize (default 12)\n\t-r specify scale of axis (default 0.5)\n\t-w specify width of produced figure (default 600, min 220)\n\t-d specify y-distance between printed items (default 20)\n",argv[0],DEFAULT_FILENAME);
	       exit(0);
	       case 'e': mode=MEASURE_EUCLDIST; break; /* Default */
	       case 'c': mode=MEASURE_CITYBLOCK; break;
	       case 'a': mode=MEASURE_VECTANGLE; break;
	       case 'l': mode=MEASURE_VECTLEN; break; 
	       case 'm': switch (argv[a][i+1]) {
	                   case '1': MERGEFLAG=MERGE_MIN; i++; break;
			   case '2': MERGEFLAG=MERGE_MEAN; i++; break;
			   default: fprintf(stderr,"Unknown merge method %c\n",argv[a][i+1]);
			 }  break;
	       case 'p': debugflag=1; break;
	       case 'f': if ((a+1)<argc) strcpy(filename,argv[++a]); goto nextarg;
	       case 's': if ((a+1)<argc) FONTSIZE=atoi(argv[++a]); goto nextarg;
	       case 'r': if ((a+1)<argc) RESOLUTION=atof(argv[++a]); goto nextarg;
	       case 'w': if ((a+1)<argc) FIGWIDTH=atoi(argv[++a]); goto nextarg;
	       case 'd': if ((a+1)<argc) DIST=atoi(argv[++a]); goto nextarg;
	       default: 
	       fprintf(stderr,"Unknown option: -%c.\n",argv[a][i]);
	       fprintf(stderr,"Usage: %s [-ealp] [-f string] [-s int] [-r float] [-w int] [-d int]\n\t-e use euclidean distance (default)\n\t-c use city-block distance\n\t-a use vector angle\n\t-l use vector length\n\t-m use merge method -m1 DISTxy=min(DISTx,DISTy) -m2 DISTxy=(DISTx+DISTy)/2\n\t-d produce debug-info\n\t-f write to specified file (default \"%s\")\n\t-s specify fontsize (default 12)\n\t-r specify scale of axis (default 0.5)\n\t-w specify width of produced figure (default 600, min 220)\n\t-d specify y-distance between printed items (default 20)\n\t-h Help\n",argv[0],DEFAULT_FILENAME); exit(1);
	    }
	 }
      }
     nextarg:;
   }
   v=(double **)malloc(sizeof(double *));
   labels=(char **)malloc(sizeof(char *));
   for (i=0; fgets(buf,10000,stdin); i++) {
      v=(double **)realloc(v,sizeof(double *)*(i+1));
      labels=(char **)realloc(labels,sizeof(char *)*(i+1));
      sprintf(str,"[%d]",i);
      labels[i]=(char *)malloc(sizeof(char)*(strlen(str)+1));
      strcpy(labels[i],str);
      v[i]=(double *)malloc(sizeof(double));
      tok=strtok(buf," ");
      j=0;
      while (tok) {
	 if ((isdigit(tok[0]) || (tok[0]=='-') || (tok[0]=='+') || (tok[0]=='.'))) {
	    v[i]=(double *)realloc(v[i],sizeof(double)*(j+1));
	    v[i][j++]=atof(tok);
	 } else if (isalpha(tok[0])) {
	    labels[i]=(char *)realloc(labels[i],sizeof(char)*(strlen(tok)+1));
	    strcpy(labels[i],removeNL(tok));
	 }
	 tok=strtok(NULL," ");
      }
      if (i==0) 
	 d=j;
      else if (d!=j) {
	 fprintf(stderr,"Vectors should have similar dimension\n");
	 fprintf(stderr,"Error in row %d. Skipping.\n",i);
	 i--;
      }
   }
   k=i; /* no of vectors */

   /* vectors in v; v[vector][elem], labels in labels; labels[vector] */

   /* calculate the distance between each vector and all others */

   m=(double **)malloc(sizeof(double *)*k);
   for (i=0; i<k; i++) {
      m[i]=(double *)malloc(sizeof(double)*k);
      for (j=0; j<k; j++)
	 if (j>i) {
	    if (mode==MEASURE_EUCLDIST)
	       m[i][j]=eucl(v[i],v[j],d);
	    else if (mode==MEASURE_CITYBLOCK)
	       m[i][j]=cityblock(v[i],v[j],d);
	    else if (mode==MEASURE_VECTLEN)
	       m[i][j]=length(v[i],v[j],d);
	    else if (mode==MEASURE_VECTANGLE)
	       m[i][j]=angle(v[i],v[j],d);
	 } else
	    m[i][j]=VOID;
   }

   head=ptr=collapse(m,k);  /* make the list */

   /* organise the list(s) */

   pairs=(POINT *)malloc(sizeof(POINT)*(k));
   for (i=0; ptr; i++) {
      max=ptr->distance;
      pairs[i].x=ptr->x;
      pairs[i].y=ptr->y;
      ptr=ptr->ptr;
   }

   order=(int *)calloc(sizeof(int),k);
   order[pairs[i-1].y]=1;
   order[pairs[i-1].x]=2;
   for (i=i-2; i>=0; i--) { 
      for (j=0; j<k; j++) {
	 if (order[j]>order[pairs[i].x])
	    order[j]++;
      }
      if (i%2) {
	 order[pairs[i].y]=order[pairs[i].y]=order[pairs[i].x]+1;
      } else {
	 order[pairs[i].y]=order[pairs[i].x];
	 order[pairs[i].x]++;
      }
   }

   points=(POINT *)malloc(sizeof(POINT)*k);
   for (i=0; i<k; i++) {
     /* points[i].x=STARTX+STRLEN; */
      points[i].x=0;
      points[i].y=STARTY+(DIST*(how_many_smaller(order,k,order[i])+1));
   }
   

   /*+++++++++++ print it as a fig-file +++++++++++++*/

   printf("Generating FIG2.1 file: %s.\n",filename, FONTSIZE, RESOLUTION);

   fp=fopen(filename,"w");
   fprintf(fp,"#FIG 2.1\n80 2\n");

   /*   for (i=0; i<k; i++)
      fprintf(fp,"4 2 -1 %d 0 -1 0 0.00000 4 12 %d %d %d %s\1\n",FONTSIZE,strlen(labels[i])*6,points[i].x,points[i].y,labels[i]); */
   
   ptr=head;
   while (ptr) {
      if (points[ptr->x].x==0) {
	fprintf(fp,"4 2 -1 %d 0 -1 0 0.00000 4 12 %d %d %d %s\1\n",FONTSIZE,
		strlen(labels[ptr->x])*6,
		STARTX+STRLEN+(int)((double)((FIGWIDTH-STARTX-STRLEN)/max)*((double)ptr->distance))-DEFAULT_DIST,
		points[ptr->x].y,
		labels[ptr->x]);
	fprintf(fp,"2 1 0 1 -1 0 0 0 0.000 -1 0 0\n\t%d %d %d %d 9999 9999\n",
		STARTX+STRLEN+(int)((double)((FIGWIDTH-STARTX-STRLEN)/max)*((double)ptr->distance))-DEFAULT_DIST,
		points[ptr->x].y,
		STARTX+STRLEN+(int)((double)((FIGWIDTH-STARTX-STRLEN)/max)*((double)ptr->distance)),
		points[ptr->x].y);
      } else {
	fprintf(fp,"2 1 0 1 -1 0 0 0 0.000 -1 0 0\n\t%d %d %d %d 9999 9999\n",
		points[ptr->x].x, 
		points[ptr->x].y,
		STARTX+STRLEN+(int)((double)((FIGWIDTH-STARTX-STRLEN)/max)*((double)ptr->distance)),
		points[ptr->x].y);
      }
      points[ptr->x].x=STARTX+STRLEN+(int)((double)((FIGWIDTH-STARTX-STRLEN)/max)*((double)ptr->distance));
      if (points[ptr->y].x==0) {
	fprintf(fp,"4 2 -1 %d 0 -1 0 0.00000 4 12 %d %d %d %s\1\n",FONTSIZE,
		strlen(labels[ptr->y])*6,
		STARTX+STRLEN+(int)((double)((FIGWIDTH-STARTX-STRLEN)/max)*((double)ptr->distance))-DEFAULT_DIST,
		points[ptr->y].y,
		labels[ptr->y]);
	fprintf(fp,"2 1 0 1 -1 0 0 0 0.000 -1 0 0\n\t%d %d %d %d 9999 9999\n",
		STARTX+STRLEN+(int)((double)((FIGWIDTH-STARTX-STRLEN)/max)*((double)ptr->distance))-DEFAULT_DIST,
		points[ptr->y].y,
		STARTX+STRLEN+(int)((double)((FIGWIDTH-STARTX-STRLEN)/max)*((double)ptr->distance)),
		points[ptr->y].y);
      } else {
	fprintf(fp,"2 1 0 1 -1 0 0 0 0.000 -1 0 0\n\t%d %d %d %d 9999 9999\n",
		points[ptr->y].x,points[ptr->y].y,
		STARTX+STRLEN+(int)((double)((FIGWIDTH-STARTX-STRLEN)/max)*((double)ptr->distance)),
		points[ptr->y].y);
      }
      points[ptr->y].x=STARTX+STRLEN+(int)((double)((FIGWIDTH-STARTX-STRLEN)/max)*((double)ptr->distance));
      fprintf(fp,"2 1 0 1 -1 0 0 0 0.000 -1 0 0\n\t%d %d %d %d 9999 9999\n",
	      points[ptr->x].x,
	      points[ptr->x].y,
	      points[ptr->y].x,
	      points[ptr->y].y);
      points[ptr->x].y=(points[ptr->x].y+points[ptr->y].y)/2;
      points[ptr->y].y=points[ptr->x].y;
      ptr=ptr->ptr;
   }
   fprintf(fp,"2 1 0 1 -1 0 0 0 0.000 -1 0 0\n\t%d %d %d %d 9999 9999\n",
	   STARTX+STRLEN,STARTY-50,FIGWIDTH,STARTY-50);
   for (i=1; i<((int)(max/RESOLUTION)+1); i++) {
      fprintf(fp,"2 1 0 1 -1 0 0 0 0.000 -1 0 0\n\t%d %d %d %d 9999 9999\n",
	      STARTX+STRLEN+(int)((double)((FIGWIDTH-STARTX-STRLEN)/max)*((double)i*RESOLUTION)),
	      STARTY-55,
	      STARTX+STRLEN+(int)((double)((FIGWIDTH-STARTX-STRLEN)/max)*((double)i*RESOLUTION)),
	      STARTY-45);
      sprintf(str,"%3.1lf",i*RESOLUTION);
      fprintf(fp,"4 1 -1 %d 0 -1 0 0.00000 4 12 %d %d %d %s\1\n",
	      FONTSIZE,
	      strlen(str)*6, 
	      STARTX+STRLEN+(int)((double)((FIGWIDTH-STARTX-STRLEN)/max)*((double)i*RESOLUTION)), 
	      STARTY-60, str);
   }
   fprintf(fp,"2 1 0 1 -1 0 0 0 0.000 -1 0 0\n\t%d %d %d %d 9999 9999\n",
	   STARTX+STRLEN, STARTY-50, STARTX+STRLEN, STARTY-40);
   sprintf(str,"%3.1lf",0.0);
   fprintf(fp,"4 1 -1 %d 0 -1 0 0.00000 4 12 %d %d %d %s\1\n",
	   FONTSIZE,
	   strlen(str)*6, STARTX+STRLEN, STARTY-20, str);
   fprintf(fp,"2 1 0 1 -1 0 0 0 0.000 -1 0 0\n\t%d %d %d %d 9999 9999\n",
	   FIGWIDTH, STARTY-50, FIGWIDTH, STARTY-40);
   sprintf(str,"%5.3lf",max);
   fprintf(fp,"4 1 -1 %d 0 -1 0 0.00000 4 12 %d %d %d %s\1\n",
	   FONTSIZE,
	   strlen(str)*6, FIGWIDTH, STARTY-20, str);
   fclose(fp);
}