forester/archive/RIO/others/phylip_mod/src/dist.c

   1 #include "phylip.h"
   2 #include "dist.h"
   3
   4 /* version 3.6. (c) Copyright 1993-2004 by the University of Washington.
   5    Written by Joseph Felsenstein, Akiko Fuseki, Sean Lamont, and Andrew Keeffe.
   6    Permission is granted to copy and use this program provided no fee is
   7    charged for it and provided that this copyright notice is not removed. */
   8
   9 void alloctree(pointptr *treenode, long nonodes)
  10 {
  11   /* allocate treenode dynamically */
  12   /* used in fitch, kitsch & neighbor */
  13   long i, j;
  14   node *p, *q;
  15
  16   *treenode = (pointptr)Malloc(nonodes*sizeof(node *));
  17   for (i = 0; i < spp; i++)
  18     (*treenode)[i] = (node *)Malloc(sizeof(node));
  19   for (i = spp; i < nonodes; i++) {
  20     q = NULL;
  21     for (j = 1; j <= 3; j++) {
  22       p = (node *)Malloc(sizeof(node));
  23       p->next = q;
  24       q = p;
  25     }
  26     p->next->next->next = p;
  27     (*treenode)[i] = p;
  28   }
  29 } /* alloctree */
  30
  31
  32 void freetree(pointptr *treenode, long nonodes)
  33 {
  34   long i, j;
  35   node *p, *q;
  36
  37   for (i = 0; i < spp; i++)
  38     free((*treenode)[i]);
  39   for (i = spp; i < nonodes; i++) {
  40     p = (*treenode)[i];
  41     for (j = 1; j <= 3; j++) {
  42       q = p;
  43       p = p->next;
  44       free(q);
  45     }
  46   }
  47   free(*treenode);
  48 } /* freetree */
  49
  50
  51 void allocd(long nonodes, pointptr treenode)
  52 {
  53   /* used in fitch & kitsch */
  54   long i, j;
  55   node *p;
  56
  57   for (i = 0; i < (spp); i++) {
  58     treenode[i]->d = (vector)Malloc(nonodes*sizeof(double));
  59   }
  60   for (i = spp; i < nonodes; i++) {
  61     p = treenode[i];
  62     for (j = 1; j <= 3; j++) {
  63       p->d = (vector)Malloc(nonodes*sizeof(double));
  64       p = p->next;
  65     }
  66   }
  67 }
  68
  69
  70 void freed(long nonodes, pointptr treenode)
  71 {
  72   /* used in fitch */
  73   long i, j;
  74   node *p;
  75
  76   for (i = 0; i < (spp); i++) {
  77     free(treenode[i]->d);
  78   }
  79   for (i = spp; i < nonodes; i++) {
  80     p = treenode[i];
  81     for (j = 1; j <= 3; j++) {
  82       free(p->d);
  83       p = p->next;
  84     }
  85   }
  86 }
  87
  88
  89 void allocw(long nonodes, pointptr treenode)
  90 {
  91   /* used in fitch & kitsch */
  92   long i, j;
  93   node *p;
  94
  95   for (i = 0; i < (spp); i++) {
  96       treenode[i]->w = (vector)Malloc(nonodes*sizeof(double));
  97   }
  98   for (i = spp; i < nonodes; i++) {
  99     p = treenode[i];
 100     for (j = 1; j <= 3; j++) {
 101       p->w = (vector)Malloc(nonodes*sizeof(double));
 102       p = p->next;
 103     }
 104   }
 105 }
 106
 107
 108 void freew(long nonodes, pointptr treenode)
 109 {
 110   /* used in fitch */
 111   long i, j;
 112   node *p;
 113
 114   for (i = 0; i < (spp); i++) {
 115     free(treenode[i]->w);
 116   }
 117   for (i = spp; i < nonodes; i++) {
 118     p = treenode[i];
 119     for (j = 1; j <= 3; j++) {
 120       free(p->w);
 121       p = p->next;
 122     }
 123   }
 124 }
 125
 126
 127 void setuptree(tree *a, long nonodes)
 128 {
 129   /* initialize a tree */
 130   /* used in fitch, kitsch, & neighbor */
 131   long i=0;
 132   node *p;
 133
 134   for (i = 1; i <= nonodes; i++) {
 135     a->nodep[i - 1]->back = NULL;
 136     a->nodep[i - 1]->tip = (i <= spp);
 137     a->nodep[i - 1]->iter = true;
 138     a->nodep[i - 1]->index = i;
 139     a->nodep[i - 1]->t = 0.0;
 140     a->nodep[i - 1]->sametime = false;
 141     a->nodep[i - 1]->v = 0.0;
 142     if (i > spp) {
 143       p = a->nodep[i-1]->next;
 144       while (p != a->nodep[i-1]) {
 145         p->back = NULL;
 146         p->tip = false;
 147         p->iter = true;
 148         p->index = i;
 149         p->t = 0.0;
 150         p->sametime = false;
 151         p = p->next;
 152       }
 153     }
 154   }
 155   a->likelihood = -1.0;
 156   a->start = a->nodep[0];
 157   a->root = NULL;
 158 }  /* setuptree */
 159
 160
 161 void inputdata(boolean replicates, boolean printdata, boolean lower,
 162                         boolean upper, vector *x, intvector *reps)
 163 {
 164   /* read in distance matrix */
 165   /* used in fitch & neighbor */
 166   long i=0, j=0, k=0, columns=0;
 167   boolean skipit=false, skipother=false;
 168
 169   if (replicates)
 170     columns = 4;
 171   else
 172     columns = 6;
 173   if (printdata) {
 174     fprintf(outfile, "\nName                       Distances");
 175     if (replicates)
 176       fprintf(outfile, " (replicates)");
 177     fprintf(outfile, "\n----                       ---------");
 178     if (replicates)
 179       fprintf(outfile, "-------------");
 180     fprintf(outfile, "\n\n");
 181   }
 182   for (i = 0; i < spp; i++) {
 183     x[i][i] = 0.0;
 184     scan_eoln(infile);
 185     initname(i);
 186     for (j = 0; j < spp; j++) {
 187       skipit = ((lower && j + 1 >= i + 1) || (upper && j + 1 <= i + 1));
 188       skipother = ((lower && i + 1 >= j + 1) || (upper && i + 1 <= j + 1));
 189       if (!skipit) {
 190         if (eoln(infile))
 191           scan_eoln(infile);
 192         if (fscanf(infile, "%lf", &x[i][j]) != 1)  {
 193           printf("The infile is of the wrong type\n");
 194           exxit(-1);
 195         }
 196         if (replicates) {
 197           if (eoln(infile))
 198             scan_eoln(infile);
 199           fscanf(infile, "%ld", &reps[i][j]);
 200         } else
 201           reps[i][j] = 1;
 202       }
 203       if (!skipit && skipother) {
 204           x[j][i] = x[i][j];
 205           reps[j][i] = reps[i][j];
 206       }
 207       if ((i == j) && (fabs(x[i][j]) > 0.000000001)) {
 208        printf("\nERROR: diagonal element of row %ld of distance matrix ", i+1);
 209         printf("is not zero.\n");
 210         printf("       Is it a distance matrix?\n\n");
 211         exxit(-1);
 212       }
 213       if ((j < i) && (fabs(x[i][j]-x[j][i]) > 0.000000001)) {
 214         printf("ERROR: distance matrix is not symmetric:\n");
 215         printf("       (%ld,%ld) element and (%ld,%ld) element are unequal.\n",
 216                i+1, j+1, j+1, i+1);
 217         printf("       They are %10.6f and %10.6f, respectively.\n",
 218                x[i][j], x[j][i]);
 219         printf("       Is it a distance matrix?\n\n");
 220         exxit(-1);
 221       }
 222     }
 223   }
 224   scan_eoln(infile);
 225   if (!printdata)
 226     return;
 227   for (i = 0; i < spp; i++) {
 228     for (j = 0; j < nmlngth; j++)
 229       putc(nayme[i][j], outfile);
 230     putc(' ', outfile);
 231     for (j = 1; j <= spp; j++) {
 232       fprintf(outfile, "%10.5f", x[i][j - 1]);
 233       if (replicates)
 234         fprintf(outfile, " (%3ld)", reps[i][j - 1]);
 235       if (j % columns == 0 && j < spp) {
 236         putc('\n', outfile);
 237         for (k = 1; k <= nmlngth + 1; k++)
 238           putc(' ', outfile);
 239       }
 240     }
 241     putc('\n', outfile);
 242   }
 243   putc('\n', outfile);
 244 }  /* inputdata */
 245
 246
 247 void coordinates(node *p, double lengthsum, long *tipy, double *tipmax,
 248                         node *start, boolean njoin)
 249 {
 250   /* establishes coordinates of nodes */
 251   node *q, *first, *last;
 252
 253   if (p->tip) {
 254     p->xcoord = (long)(over * lengthsum + 0.5);
 255     p->ycoord = *tipy;
 256     p->ymin = *tipy;
 257     p->ymax = *tipy;
 258     (*tipy) += down;
 259     if (lengthsum > *tipmax)
 260       *tipmax = lengthsum;
 261     return;
 262   }
 263   q = p->next;
 264   do {
 265     if (q->back)
 266       coordinates(q->back, lengthsum + q->v, tipy,tipmax, start, njoin);
 267     q = q->next;
 268   } while ((p == start || p != q) && (p != start || p->next != q));
 269   first = p->next->back;
 270   q = p;
 271   while (q->next != p && q->next->back)  /* is this right ? */
 272     q = q->next;
 273   last = q->back;
 274   p->xcoord = (long)(over * lengthsum + 0.5);
 275   if (p == start && p->back)
 276     p->ycoord = p->next->next->back->ycoord;
 277   else
 278     p->ycoord = (first->ycoord + last->ycoord) / 2;
 279   p->ymin = first->ymin;
 280   p->ymax = last->ymax;
 281 }  /* coordinates */
 282
 283
 284 void drawline(long i, double scale, node *start, boolean rooted)
 285 {
 286   /* draws one row of the tree diagram by moving up tree */
 287   node *p, *q;
 288   long n=0, j=0;
 289   boolean extra=false, trif=false;
 290   node *r, *first =NULL, *last =NULL;
 291   boolean done=false;
 292
 293   p = start;
 294   q = start;
 295   extra = false;
 296   trif = false;
 297   if (i == (long)p->ycoord && p == start) {  /* display the root */
 298     if (rooted) {
 299       if (p->index - spp >= 10)
 300         fprintf(outfile, "-");
 301       else
 302         fprintf(outfile, "--");
 303     }
 304     else {
 305       if (p->index - spp >= 10)
 306         fprintf(outfile, " ");
 307       else
 308         fprintf(outfile, "  ");
 309     }
 310     if (p->index - spp >= 10)
 311       fprintf(outfile, "%2ld", p->index - spp);
 312     else
 313       fprintf(outfile, "%ld", p->index - spp);
 314     extra = true;
 315     trif = true;
 316   } else
 317     fprintf(outfile, "  ");
 318   do {
 319     if (!p->tip) { /* internal nodes */
 320       r = p->next;
 321       /* r->back here is going to the same node. */
 322       do {
 323         if (!r->back) {
 324           r = r->next;
 325           continue;
 326         }
 327         if (i >= r->back->ymin && i <= r->back->ymax) {
 328           q = r->back;
 329           break;
 330         }
 331         r = r->next;
 332       } while (!((p != start && r == p) || (p == start && r == p->next)));
 333       first = p->next->back;
 334       r = p;
 335       while (r->next != p)
 336         r = r->next;
 337       last = r->back;
 338       if (!rooted && (p == start))
 339         last = p->back;
 340     } /* end internal node case... */
 341     /* draw the line: */
 342     done = (p->tip || p == q);
 343     n = (long)(scale * (q->xcoord - p->xcoord) + 0.5);
 344     if (!q->tip) {
 345       if ((n < 3) && (q->index - spp >= 10))
 346         n = 3;
 347       if ((n < 2) && (q->index - spp < 10))
 348         n = 2;
 349     }
 350     if (extra) {
 351       n--;
 352       extra = false;
 353     }
 354     if ((long)q->ycoord == i && !done) {
 355       if (p->ycoord != q->ycoord)
 356         putc('+', outfile);
 357       if (trif) {
 358         n++;
 359         trif = false;
 360       }
 361       if (!q->tip) {
 362         for (j = 1; j <= n - 2; j++)
 363           putc('-', outfile);
 364         if (q->index - spp >= 10)
 365           fprintf(outfile, "%2ld", q->index - spp);
 366         else
 367           fprintf(outfile, "-%ld", q->index - spp);
 368         extra = true;
 369       } else {
 370         for (j = 1; j < n; j++)
 371           putc('-', outfile);
 372       }
 373     } else if (!p->tip) {
 374       if ((long)last->ycoord > i && (long)first->ycoord < i
 375            && i != (long)p->ycoord) {
 376         putc('!', outfile);
 377         for (j = 1; j < n; j++)
 378           putc(' ', outfile);
 379       } else {
 380         for (j = 1; j <= n; j++)
 381           putc(' ', outfile);
 382         trif = false;
 383       }
 384     }
 385     if (q != p)
 386       p = q;
 387   } while (!done);
 388   if ((long)p->ycoord == i && p->tip) {
 389     for (j = 0; j < nmlngth; j++)
 390       putc(nayme[p->index - 1][j], outfile);
 391   }
 392   putc('\n', outfile);
 393 }  /* drawline */
 394
 395
 396 void printree(node *start, boolean treeprint,
 397                         boolean        njoin, boolean rooted)
 398 {
 399   /* prints out diagram of the tree */
 400   /* used in fitch & neighbor */
 401   long i;
 402   long tipy;
 403   double scale,tipmax;
 404
 405   if (!treeprint)
 406     return;
 407   putc('\n', outfile);
 408   tipy = 1;
 409   tipmax = 0.0;
 410   coordinates(start, 0.0, &tipy, &tipmax, start, njoin);
 411   scale = 1.0 / (long)(tipmax + 1.000);
 412   for (i = 1; i <= (tipy - down); i++)
 413     drawline(i, scale, start, rooted);
 414   putc('\n', outfile);
 415 }  /* printree */
 416
 417
 418 void treeoutr(node *p, long *col, tree *curtree)
 419 {
 420   /* write out file with representation of final tree.
 421    * Rooted case. Used in kitsch and neighbor. */
 422   long i, n, w;
 423   Char c;
 424   double x;
 425
 426   if (p->tip) {
 427     n = 0;
 428     for (i = 1; i <= nmlngth; i++) {
 429       if (nayme[p->index - 1][i - 1] != ' ')
 430         n = i;
 431     }
 432     for (i = 0; i < n; i++) {
 433       c = nayme[p->index - 1][i];
 434       if (c == ' ')
 435         c = '_';
 436       putc(c, outtree);
 437     }
 438     (*col) += n;
 439   } else {
 440     putc('(', outtree);
 441     (*col)++;
 442     treeoutr(p->next->back,col,curtree);
 443     putc(',', outtree);
 444     (*col)++;
 445     if ((*col) > 55) {
 446       putc('\n', outtree);
 447       (*col) = 0;
 448     }
 449     treeoutr(p->next->next->back,col,curtree);
 450     putc(')', outtree);
 451     (*col)++;
 452   }
 453   x = p->v;
 454   if (x > 0.0)
 455     w = (long)(0.43429448222 * log(x));
 456   else if (x == 0.0)
 457     w = 0;
 458   else
 459     w = (long)(0.43429448222 * log(-x)) + 1;
 460   if (w < 0)
 461     w = 0;
 462   if (p == curtree->root)
 463     fprintf(outtree, ";\n");
 464   else {
 465     fprintf(outtree, ":%*.5f", (int)(w + 7), x);
 466     (*col) += w + 8;
 467   }
 468 }  /* treeoutr */
 469
 470
 471 void treeout(node *p, long *col, double m, boolean njoin, node *start)
 472 {
 473   /* write out file with representation of final tree */
 474   /* used in fitch & neighbor */
 475   long i=0, n=0, w=0;
 476   Char c;
 477   double x=0.0;
 478
 479   if (p->tip) {
 480     n = 0;
 481     for (i = 1; i <= nmlngth; i++) {
 482       if (nayme[p->index - 1][i - 1] != ' ')
 483         n = i;
 484     }
 485     for (i = 0; i < n; i++) {
 486       c = nayme[p->index - 1][i];
 487       if (c == ' ')
 488         c = '_';
 489       putc(c, outtree);
 490     }
 491     *col += n;
 492   } else {
 493     putc('(', outtree);
 494     (*col)++;
 495     treeout(p->next->back, col, m, njoin, start);
 496     putc(',', outtree);
 497     (*col)++;
 498     if (*col > 55) {
 499       putc('\n', outtree);
 500       *col = 0;
 501     }
 502     treeout(p->next->next->back, col, m, njoin, start);
 503     if (p == start && njoin) {
 504       putc(',', outtree);
 505       treeout(p->back, col, m, njoin, start);
 506     }
 507     putc(')', outtree);
 508     (*col)++;
 509   }
 510   x = p->v;
 511   if (x > 0.0)
 512     w = (long)(m * log(x));
 513   else if (x == 0.0)
 514     w = 0;
 515   else
 516     w = (long)(m * log(-x)) + 1;
 517   if (w < 0)
 518     w = 0;
 519   if (p == start)
 520     fprintf(outtree, ";\n");
 521   else {
 522     fprintf(outtree, ":%*.5f", (int) w + 7, x);
 523     *col += w + 8;
 524   }
 525 }  /* treeout */
 526