--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "distance_matrix.h"
+#include "utils.h"
+
+#define PUBLIC
+#define PRIVATE static
+
+#define MAXTAXA_FOR_LABELS 200
+
+typedef struct{
+ int set1;
+ int set2;
+ float distance;
+ float distance2;
+ } Union;
+
+typedef struct _Node_ {
+ float height;
+ float brr;
+ float brl;
+ int whoami;
+ int size;
+ struct _Node_ *father;
+ struct _Node_ *left;
+ struct _Node_ *right;
+ } Node;
+
+PUBLIC void PSplot_phylogeny(Union *cluster, char *filename, char *type);
+
+PRIVATE Node *W2Phylo(Union *cluster);
+PRIVATE Node *Nj2Phylo(Union *cluster);
+PRIVATE void free_phylo_tree(Node *root);
+PRIVATE void fill_height_into_tree(Node *root);
+PRIVATE void fill_br_from_height(Node *root);
+PRIVATE void plot_branch(Node *root, FILE *fp);
+PRIVATE void format_number(float x);
+
+PRIVATE char str[10];
+PRIVATE float threshold;
+PRIVATE int print_labels=1;
+
+/* --------------------------------------------------------------------------*/
+
+PUBLIC void PSplot_phylogeny(Union *cluster, char *filename, char *type)
+{
+ int n;
+ char outfile[50];
+ char *tmp, *tfont;
+ Node *root;
+ Node *tempnode;
+ FILE *fp;
+ float xsize, ysize, tfontsize, lfontsize, lwidth;
+
+ n = cluster[0].set1;
+ switch(type[0]) {
+ case 'W' :
+ root = W2Phylo(cluster);
+ break;
+ case 'N' :
+ root = Nj2Phylo(cluster);
+ break;
+ default :
+ return;
+ }
+
+ outfile[0] = '\0';
+ tmp = get_taxon_label(-1); /* retrieve data set identifier */
+ if(tmp) {strcat(outfile,tmp); strcat(outfile,"_"); free(tmp);}
+ strcat(outfile,filename);
+
+ fp = fopen(outfile,"w");
+ if (fp==NULL) {
+ fprintf(stderr,"couldn't open %s -- not doing treeplot\n", outfile);
+ return;
+ }
+ xsize = root->size;
+ tempnode = root;
+ while(tempnode->right)
+ tempnode = tempnode->right;
+ ysize = tempnode->height;
+ threshold = ysize*0.06;
+
+ lwidth = MIN2(1.5, 8/sqrt((double)n));
+ tfontsize = MIN2(15,550./(10+n));
+ lfontsize = 2./3.*tfontsize;
+ if (n>30) tfont = "Helvetica";
+ else tfont = "Times-Roman";
+ if(n > MAXTAXA_FOR_LABELS) print_labels=0;
+ else print_labels=1;
+ fprintf(fp,"%%!PS-Adobe-2.0 EPSF-1.2\n");
+ fprintf(fp,"%%%%Title: TreePlot (%s)\n",type);
+ fprintf(fp,"%%%%Creator: AnalyseDists\n");
+ fprintf(fp,"%%%%CreationDate: %s", time_stamp());
+ /* BoundingBox is only approximate */
+ fprintf(fp,"%%%%BoundingBox: 35 45 535 640\n");
+ fprintf(fp,"%%%%Pages: 1\n");
+ fprintf(fp,"%%%%EndComments: No comment!\n");
+ fprintf(fp,"288.5 50 translate\n");
+ fprintf(fp,"%3.1f setlinewidth\n"
+ "/cmtx matrix currentmatrix def\n", lwidth);
+ fprintf(fp,"500 %g div 360 %g div scale\n", xsize, ysize);
+ fprintf(fp,"/rotshow {gsave cmtx setmatrix\n");
+ fprintf(fp," 90 rotate 5 %4.1f rmoveto show grestore} def\n",
+ -tfontsize/4);
+ fprintf(fp,"/cshow {gsave cmtx setmatrix\n"
+ " /Helvetica findfont %4.1f scalefont setfont\n"
+ " 90 rotate 0 %3.1f rmoveto\n"
+ " dup stringwidth pop 2 div neg 0 rmoveto show\n"
+ " grestore} def\n", lfontsize, lfontsize/5);
+ fprintf(fp,"/%s findfont %4.1f scalefont setfont\n", tfont, tfontsize);
+
+ fprintf(fp,"0 0 moveto\n");
+ plot_branch(root,fp);
+
+ fprintf(fp,"cmtx setmatrix stroke\nshowpage\n");
+ free_phylo_tree(root);
+ fclose(fp);
+
+}
+
+/* --------------------------------------------------------------------------*/
+
+PRIVATE Node *W2Phylo(Union *cluster)
+{
+ int i,n;
+ float b;
+ Node **taxa;
+ Node *father;
+
+ n=cluster[0].set1;
+ taxa = (Node **) space(sizeof(Node*)*(n+1));
+
+ b = sqrt(MAX2(cluster[n-1].distance,0.));
+ for(i=1;i<=n;i++){
+ taxa[i] = (Node *) space(sizeof(Node));
+ taxa[i]->whoami = i;
+ taxa[i]->size = 1;
+ taxa[i]->height = b;
+ }
+
+ for(i=1;i<n;i++) {
+ father = (Node *) space(sizeof(Node));
+ father->whoami = 0;
+ father->left = taxa[cluster[i].set1];
+ father->right = taxa[cluster[i].set2];
+ father->size = father->left->size+father->right->size;
+ b = sqrt(MAX2(cluster[n-1].distance-cluster[i].distance,0.));
+ father->height = b;
+ father->left->father = father;
+ father->right->father = father;
+ taxa[cluster[i].set1] = father;
+ }
+ free(taxa);
+ father->whoami = -1;
+ fill_br_from_height(father);
+ return father;
+}
+
+
+/* --------------------------------------------------------------------------*/
+
+PRIVATE void fill_br_from_height(Node *root)
+{
+ root->brr = root->right->height - root->height;
+ root->brl = root->left->height - root->height;
+ if(!root->left->whoami) fill_br_from_height(root->left);
+ if(!root->right->whoami) fill_br_from_height(root->right);
+ return;
+}
+
+/* -------------------------------------------------------------------------- */
+#define ORDER(X) \
+if ((X)->brl>(X)->brr) { \
+ xnode=(X)->right; \
+ (X)->right=(X)->left; \
+ (X)->left=xnode; \
+ br = (X)->brr; \
+ (X)->brr=(X)->brl; \
+ (X)->brl=br; \
+}
+
+
+
+PRIVATE Node *Nj2Phylo(Union *cluster)
+{
+ int i,n, br;
+ float h1,h2,maxdist,dist;
+ Node **taxa;
+ Node *father;
+ Node *topnode; /* of the longest path in the tree */
+ Node *tempnode, *xnode;
+ Node *root;
+
+ maxdist = 0.;
+ n=cluster[0].set1;
+ taxa = (Node **) space(sizeof(Node*)*(n+1));
+
+ for(i=1;i<=n;i++){
+ taxa[i] = (Node *) space(sizeof(Node));
+ taxa[i]->whoami = i;
+ taxa[i]->size = 1;
+ }
+
+ for(i=1;i<n;i++) {
+ father = (Node *) space(sizeof(Node));
+ father->whoami = 0;
+ h1 = cluster[i].distance +taxa[cluster[i].set1]->height;
+ h2 = cluster[i].distance2+taxa[cluster[i].set2]->height;
+ dist = h1 + h2;
+ if(dist > maxdist) {
+ maxdist = dist;
+ topnode = father;
+ }
+ if(h1<=h2) {
+ father->brr = cluster[i].distance2;
+ father->brl = cluster[i].distance;
+ father->left = taxa[cluster[i].set1];
+ father->right = taxa[cluster[i].set2];
+ }
+ else {
+ father->brr = cluster[i].distance;
+ father->brl = cluster[i].distance2;
+ father->left = taxa[cluster[i].set2];
+ father->right = taxa[cluster[i].set1];
+ }
+ father->height = MAX2(h1,h2);
+ father->size = father->left->size + father->right->size;
+ father->left->father = father;
+ father->right->father = father;
+ taxa[cluster[i].set1] = father;
+ }
+
+ /* new root is inserted in the middle of the longest path */
+ dist = maxdist/2.;
+ tempnode = topnode;
+ while (tempnode->height > dist)
+ tempnode = tempnode->right;
+
+ topnode = tempnode->father;
+ /* new root must go between topnode and tempnode */
+ h1 = dist - tempnode->height;
+ h2 = topnode->height - dist ;
+
+ root = (Node *) space(sizeof(Node));
+ root->father = NULL;
+ root->whoami = -1;
+ root->right = tempnode;
+ root->left = topnode;
+ root->height = 0.0;
+ root->brr = h1;
+ root->brl = h2;
+ root->size = n;
+
+ topnode->size = n - tempnode->size;
+
+ tempnode = root;
+ while(topnode->father) {
+ topnode->right = topnode->father;
+ topnode->brr = topnode->father->brr;
+ topnode->father = tempnode;
+ tempnode = topnode;
+ topnode = topnode->right;
+ topnode->size = tempnode->size - tempnode->left->size;
+ ORDER(tempnode);
+ }
+ /* topnode is now old root, remove it */
+ if (tempnode->right==topnode) {
+ tempnode->right = topnode->left;
+ tempnode->brr = topnode->brr + topnode->brl;
+ tempnode->right->father = tempnode;
+ } else {
+ tempnode->left = topnode->left;
+ tempnode->brl = topnode->brr + topnode->brl;
+ tempnode->left->father = tempnode;
+ }
+ free(topnode);
+ ORDER(root);
+
+ fill_height_into_tree(root);
+
+ return root;
+}
+
+/* -------------------------------------------------------------------------- */
+
+PRIVATE void fill_height_into_tree(Node *root)
+{
+ if(root->whoami>0) return;
+ root->left->height = root->height + root->brl;
+ root->right->height = root->height + root->brr;
+ fill_height_into_tree(root->left);
+ fill_height_into_tree(root->right);
+ return;
+}
+
+/* -------------------------------------------------------------------------- */
+
+PRIVATE void free_phylo_tree(Node *root) {
+ if(root->left) free_phylo_tree(root->left);
+ if(root->right) free_phylo_tree(root->right);
+ free(root);
+}
+
+
+/* -------------------------------------------------------------------------- */
+
+PRIVATE void plot_branch(Node *root, FILE *fp)
+{
+ char *label;
+
+ fprintf(fp,"currentpoint %g 0 rlineto 0 %g rlineto \n",
+ -(float)(root->left->size)/2., root->brr);
+ if((print_labels)&&(root->brr > threshold)) {
+ format_number(root->brr);
+ fprintf(fp,"currentpoint 0 %g rlineto \n", -root->brr/2.);
+ fprintf(fp,"(%s) cshow moveto\n", str);
+ }
+ if(root->right->whoami==0) plot_branch(root->right,fp);
+ else if(print_labels) {
+ label = get_taxon_label(root->right->whoami);
+ fprintf(fp, "(%s) rotshow\n", label);
+ free(label);
+ }
+ fprintf(fp,"moveto\n");
+ fprintf(fp,"currentpoint %g 0 rlineto 0 %g rlineto \n",
+ +(float)(root->right->size)/2., root->brl);
+ if((print_labels)&&(root->brl > threshold)) {
+ format_number(root->brl);
+ fprintf(fp,"currentpoint 0 %g rlineto \n", -root->brl/2.);
+ fprintf(fp,"(%s) cshow moveto\n",str);
+ }
+ if(root->left->whoami==0) plot_branch(root->left,fp);
+ else if(print_labels) {
+ label = get_taxon_label(root->left->whoami);
+ fprintf(fp, "(%s) rotshow\n", label);
+ free(label);
+ }
+ fprintf(fp,"moveto\n");
+}
+
+/* -------------------------------------------------------------------------- */
+
+PRIVATE void format_number(float x)
+{
+ if(x >= 1.e5) { sprintf(str, "%4.1g",x ); return;}
+ if(x >= 100) { sprintf(str, "%.0f" ,x ); return;}
+ if(x >= 10) { sprintf(str, "%.1f", x ); return;}
+ if(x >= 1) { sprintf(str, "%.2f", x ); return;}
+ if(x >= 0.1) { sprintf(str, "%.2f", x ); return;}
+ if(x >= 0.01) { sprintf(str, "%.3f", x ); return;}
+ sprintf(str, "%4.1g",x );
+ return;
+}
+
+/* -------------------------------------------------------------------------- */
git://source.jalview.org / jabaws.git / blobdiff