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

   1
   2 #include "phylip.h"
   3 #include "dist.h"
   4
   5 /* version 3.6. (c) Copyright 1993-2004 by the University of Washington.
   6    Written by Mary Kuhner, Jon Yamato, Joseph Felsenstein, Akiko Fuseki,
   7    Sean Lamont, and Andrew Keeffe.
   8    Permission is granted to copy and use this program provided no fee is
   9    charged for it and provided that this copyright notice is not removed. */
  10
  11
  12 #ifndef OLDC
  13 /* function prototypes */
  14 void getoptions(void);
  15 void allocrest(void);
  16 void doinit(void);
  17 void inputoptions(void);
  18 void getinput(void);
  19 void describe(node *, double);
  20 void summarize(void);
  21 void nodelabel(boolean);
  22 void jointree(void);
  23 void maketree(void);
  24 void freerest(void);
  25 /* function prototypes */
  26 #endif
  27
  28
  29 Char infilename[FNMLNGTH], outfilename[FNMLNGTH], outtreename[FNMLNGTH];
  30 long nonodes2, outgrno, col, datasets, ith;
  31 long inseed;
  32 vector *x;
  33 intvector *reps;
  34 boolean jumble, lower, upper, outgropt, replicates, trout,
  35                printdata, progress, treeprint, mulsets, njoin;
  36 tree curtree;
  37 longer seed;
  38 long *enterorder;
  39 Char progname[20];
  40
  41 /* variables for maketree, propagated globally for C version: */
  42 node **cluster;
  43
  44
  45 void getoptions()
  46 {
  47   /* interactively set options */
  48   long inseed0 = 0, loopcount;
  49   Char ch;
  50
  51   fprintf(outfile, "\nNeighbor-Joining/UPGMA method version %s\n\n",VERSION);
  52   putchar('\n');
  53   jumble = false;
  54   lower = false;
  55   outgrno = 1;
  56   outgropt = false;
  57   replicates = false;
  58   trout = true;
  59   upper = false;
  60   printdata = false;
  61   progress = true;
  62   treeprint = true;
  63   njoin = true;
  64   loopcount = 0;
  65   for(;;) {
  66     cleerhome();
  67     printf("\nNeighbor-Joining/UPGMA method version %s\n\n",VERSION);
  68     printf("Settings for this run:\n");
  69     printf("  N       Neighbor-joining or UPGMA tree?  %s\n",
  70            (njoin ? "Neighbor-joining" : "UPGMA"));
  71     if (njoin) {
  72       printf("  O                        Outgroup root?");
  73       if (outgropt)
  74         printf("  Yes, at species number%3ld\n", outgrno);
  75       else
  76         printf("  No, use as outgroup species%3ld\n", outgrno);
  77     }
  78     printf("  L         Lower-triangular data matrix?  %s\n",
  79            (lower ? "Yes" : "No"));
  80     printf("  R         Upper-triangular data matrix?  %s\n",
  81            (upper ? "Yes" : "No"));
  82     printf("  S                        Subreplicates?  %s\n",
  83            (replicates ? "Yes" : "No"));
  84     printf("  J     Randomize input order of species?");
  85     if (jumble)
  86       printf("  Yes (random number seed =%8ld)\n", inseed0);
  87     else
  88       printf("  No. Use input order\n");
  89     printf("  M           Analyze multiple data sets?");
  90     if (mulsets)
  91       printf("  Yes, %2ld sets\n", datasets);
  92     else
  93       printf("  No\n");
  94     printf("  0   Terminal type (IBM PC, ANSI, none)?  %s\n",
  95            (ibmpc ? "IBM PC" : ansi ? "ANSI" : "(none)"));
  96     printf("  1    Print out the data at start of run  %s\n",
  97            (printdata ? "Yes" : "No"));
  98     printf("  2  Print indications of progress of run  %s\n",
  99            (progress ? "Yes" : "No"));
 100     printf("  3                        Print out tree  %s\n",
 101            (treeprint ? "Yes" : "No"));
 102     printf("  4       Write out trees onto tree file?  %s\n",
 103            (trout ? "Yes" : "No"));
 104     printf("\n\n  Y to accept these or type the letter for one to change\n");
 105 #ifdef WIN32
 106     phyFillScreenColor();
 107 #endif
 108     scanf("%c%*[^\n]", &ch);
 109     getchar();
 110     if (ch == '\n')
 111       ch = ' ';
 112     uppercase(&ch);
 113     if  (ch == 'Y')
 114       break;
 115     if (strchr("NJOULRSM01234",ch) != NULL){
 116       switch (ch) {
 117
 118       case 'J':
 119         jumble = !jumble;
 120          if (jumble)
 121           initseed(&inseed, &inseed0, seed);
 122         break;
 123
 124       case 'L':
 125         lower = !lower;
 126         break;
 127
 128       case 'O':
 129         outgropt = !outgropt;
 130         if (outgropt)
 131           initoutgroup(&outgrno, spp);
 132         else
 133           outgrno = 1;
 134         break;
 135
 136       case 'R':
 137         upper = !upper;
 138         break;
 139
 140       case 'S':
 141         replicates = !replicates;
 142         break;
 143
 144       case 'N':
 145         njoin = !njoin;
 146         break;
 147
 148       case 'M':
 149         mulsets = !mulsets;
 150         if (mulsets)
 151           initdatasets(&datasets);
 152         jumble = true;
 153          if (jumble)
 154           initseed(&inseed, &inseed0, seed);
 155         break;
 156
 157       case '0':
 158         initterminal(&ibmpc, &ansi);
 159         break;
 160
 161       case '1':
 162         printdata = !printdata;
 163         break;
 164
 165       case '2':
 166         progress = !progress;
 167         break;
 168
 169       case '3':
 170         treeprint = !treeprint;
 171         break;
 172
 173       case '4':
 174         trout = !trout;
 175         break;
 176       }
 177     } else
 178       printf("Not a possible option!\n");
 179     countup(&loopcount, 100);
 180   }
 181 }  /* getoptions */
 182
 183
 184 void allocrest()
 185 {
 186   long i;
 187
 188   x = (vector *)Malloc(spp*sizeof(vector));
 189   for (i = 0; i < spp; i++)
 190     x[i] = (vector)Malloc(spp*sizeof(double));
 191   reps = (intvector *)Malloc(spp*sizeof(intvector));
 192   for (i = 0; i < spp; i++)
 193     reps[i] = (intvector)Malloc(spp*sizeof(long));
 194   nayme = (naym *)Malloc(spp*sizeof(naym));
 195   enterorder = (long *)Malloc(spp*sizeof(long));
 196   cluster = (node **)Malloc(spp*sizeof(node *));
 197 }  /* allocrest */
 198
 199
 200 void freerest()
 201 {
 202   long i;
 203
 204   for (i = 0; i < spp; i++)
 205     free(x[i]);
 206   free(x);
 207   for (i = 0; i < spp; i++)
 208     free(reps[i]);
 209   free(reps);
 210   free(nayme);
 211   free(enterorder);
 212   free(cluster);
 213 }  /* freerest */
 214
 215
 216 void doinit()
 217 {
 218   /* initializes variables */
 219   node *p;
 220
 221   inputnumbers2(&spp, &nonodes2, 2);
 222   nonodes2 += (njoin ? 0 : 1);
 223   getoptions();
 224   alloctree(&curtree.nodep, nonodes2+1);
 225   p = curtree.nodep[nonodes2]->next->next;
 226   curtree.nodep[nonodes2]->next = curtree.nodep[nonodes2];
 227   free(p);
 228   allocrest();
 229
 230 }  /* doinit */
 231
 232
 233 void inputoptions()
 234 {
 235   /* read options information */
 236
 237   if (ith != 1)
 238     samenumsp2(ith);
 239   putc('\n', outfile);
 240   if (njoin)
 241     fprintf(outfile, " Neighbor-joining method\n");
 242   else
 243     fprintf(outfile, " UPGMA method\n");
 244   fprintf(outfile, "\n Negative branch lengths allowed\n\n");
 245 }  /* inputoptions */
 246
 247
 248 void describe(node *p, double height)
 249 {
 250   /* print out information for one branch */
 251   long i;
 252   node *q;
 253
 254   q = p->back;
 255   if (njoin)
 256     fprintf(outfile, "%4ld          ", q->index - spp);
 257   else
 258     fprintf(outfile, "%4ld     ", q->index - spp);
 259   if (p->tip) {
 260     for (i = 0; i < nmlngth; i++)
 261       putc(nayme[p->index - 1][i], outfile);
 262     putc(' ', outfile);
 263   } else {
 264     if (njoin)
 265       fprintf(outfile, "%4ld       ", p->index - spp);
 266     else {
 267       fprintf(outfile, "%4ld       ", p->index - spp);
 268     }
 269   }
 270   if (njoin)
 271     fprintf(outfile, "%12.5f\n", q->v);
 272   else
 273     fprintf(outfile, "%10.5f      %10.5f\n", q->v, q->v+height);
 274   if (!p->tip) {
 275     describe(p->next->back, height+q->v);
 276     describe(p->next->next->back, height+q->v);
 277   }
 278 }  /* describe */
 279
 280
 281 void summarize()
 282 {
 283   /* print out branch lengths etc. */
 284   putc('\n', outfile);
 285   if (njoin) {
 286     fprintf(outfile, "remember:");
 287     if (outgropt)
 288       fprintf(outfile, " (although rooted by outgroup)");
 289     fprintf(outfile, " this is an unrooted tree!\n");
 290   }
 291   if (njoin) {
 292     fprintf(outfile, "\nBetween        And            Length\n");
 293     fprintf(outfile, "-------        ---            ------\n");
 294   } else {
 295     fprintf(outfile, "From     To            Length          Height\n");
 296     fprintf(outfile, "----     --            ------          ------\n");
 297   }
 298   describe(curtree.start->next->back, 0.0);
 299   describe(curtree.start->next->next->back, 0.0);
 300   if (njoin)
 301     describe(curtree.start->back, 0.0);
 302   fprintf(outfile, "\n\n");
 303 }  /* summarize */
 304
 305
 306 void nodelabel(boolean isnode)
 307 {
 308   if (isnode)
 309     printf("node");
 310   else
 311     printf("species");
 312 }  /* nodelabel */
 313
 314
 315 void jointree()
 316 {
 317   /* calculate the tree */
 318   long nc, nextnode, mini=0, minj=0, i, j, ia, ja, ii, jj, nude, iter;
 319   double fotu2, total, tmin, dio, djo, bi, bj, bk, dmin=0, da;
 320   long el[3];
 321   vector av;
 322   intvector oc;
 323
 324   double *R;   /* added in revisions by Y. Ina */
 325   R = (double *)Malloc(spp * sizeof(double));
 326
 327   for (i = 0; i <= spp - 2; i++) {
 328     for (j = i + 1; j < spp; j++) {
 329       da = (x[i][j] + x[j][i]) / 2.0;
 330       x[i][j] = da;
 331       x[j][i] = da;
 332     }
 333   }
 334   /* First initialization */
 335   fotu2 = spp - 2.0;
 336   nextnode = spp + 1;
 337   av = (vector)Malloc(spp*sizeof(double));
 338   oc = (intvector)Malloc(spp*sizeof(long));
 339   for (i = 0; i < spp; i++) {
 340     av[i] = 0.0;
 341     oc[i] = 1;
 342   }
 343   /* Enter the main cycle */
 344   if (njoin)
 345     iter = spp - 3;
 346   else
 347     iter = spp - 1;
 348   for (nc = 1; nc <= iter; nc++) {
 349     for (j = 2; j <= spp; j++) {
 350       for (i = 0; i <= j - 2; i++)
 351         x[j - 1][i] = x[i][j - 1];
 352     }
 353     tmin = 99999.0;
 354     /* Compute sij and minimize */
 355     if (njoin) {     /* many revisions by Y. Ina from here ... */
 356       for (i = 0; i < spp; i++)
 357         R[i] = 0.0;
 358       for (ja = 2; ja <= spp; ja++) {
 359         jj = enterorder[ja - 1];
 360         if (cluster[jj - 1] != NULL) {
 361           for (ia = 0; ia <= ja - 2; ia++) {
 362             ii = enterorder[ia];
 363             if (cluster[ii - 1] != NULL) {
 364               R[ii - 1] += x[ii - 1][jj - 1];
 365               R[jj - 1] += x[ii - 1][jj - 1];
 366             }
 367           }
 368         }
 369       }
 370     } /* ... to here */
 371     for (ja = 2; ja <= spp; ja++) {
 372       jj = enterorder[ja - 1];
 373       if (cluster[jj - 1] != NULL) {
 374         for (ia = 0; ia <= ja - 2; ia++) {
 375           ii = enterorder[ia];
 376           if (cluster[ii - 1] != NULL) {
 377             if (njoin) {
 378               total = fotu2 * x[ii - 1][jj - 1] - R[ii - 1] - R[jj - 1];
 379                /* this statement part of revisions by Y. Ina */
 380             } else
 381               total = x[ii - 1][jj - 1];
 382             if (total < tmin) {
 383               tmin = total;
 384               mini = ii;
 385               minj = jj;
 386             }
 387           }
 388         }
 389       }
 390     }
 391     /* compute lengths and print */
 392     if (njoin) {
 393       dio = 0.0;
 394       djo = 0.0;
 395       for (i = 0; i < spp; i++) {
 396         dio += x[i][mini - 1];
 397         djo += x[i][minj - 1];
 398       }
 399       dmin = x[mini - 1][minj - 1];
 400       dio = (dio - dmin) / fotu2;
 401       djo = (djo - dmin) / fotu2;
 402       bi = (dmin + dio - djo) * 0.5;
 403       bj = dmin - bi;
 404       bi -= av[mini - 1];
 405       bj -= av[minj - 1];
 406     } else {
 407       bi = x[mini - 1][minj - 1] / 2.0 - av[mini - 1];
 408       bj = x[mini - 1][minj - 1] / 2.0 - av[minj - 1];
 409       av[mini - 1] += bi;
 410     }
 411     if (progress) {
 412       printf("Cycle %3ld: ", iter - nc + 1);
 413       if (njoin)
 414         nodelabel((boolean)(av[mini - 1] > 0.0));
 415       else
 416         nodelabel((boolean)(oc[mini - 1] > 1.0));
 417       printf(" %ld (%10.5f) joins ", mini, bi);
 418       if (njoin)
 419         nodelabel((boolean)(av[minj - 1] > 0.0));
 420       else
 421         nodelabel((boolean)(oc[minj - 1] > 1.0));
 422       printf(" %ld (%10.5f)\n", minj, bj);
 423 #ifdef WIN32
 424       phyFillScreenColor();
 425 #endif
 426     }
 427     hookup(curtree.nodep[nextnode - 1]->next, cluster[mini - 1]);
 428     hookup(curtree.nodep[nextnode - 1]->next->next, cluster[minj - 1]);
 429     cluster[mini - 1]->v = bi;
 430     cluster[minj - 1]->v = bj;
 431     cluster[mini - 1]->back->v = bi;
 432     cluster[minj - 1]->back->v = bj;
 433     cluster[mini - 1] = curtree.nodep[nextnode - 1];
 434     cluster[minj - 1] = NULL;
 435     nextnode++;
 436     if (njoin)
 437       av[mini - 1] = dmin * 0.5;
 438     /* re-initialization */
 439     fotu2 -= 1.0;
 440     for (j = 0; j < spp; j++) {
 441       if (cluster[j] != NULL) {
 442         if (njoin) {
 443           da = (x[mini - 1][j] + x[minj - 1][j]) * 0.5;
 444           if (mini - j - 1 < 0)
 445             x[mini - 1][j] = da;
 446           if (mini - j - 1 > 0)
 447             x[j][mini - 1] = da;
 448         } else {
 449           da = x[mini - 1][j] * oc[mini - 1] + x[minj - 1][j] * oc[minj - 1];
 450           da /= oc[mini - 1] + oc[minj - 1];
 451           x[mini - 1][j] = da;
 452           x[j][mini - 1] = da;
 453         }
 454       }
 455     }
 456     for (j = 0; j < spp; j++) {
 457       x[minj - 1][j] = 0.0;
 458       x[j][minj - 1] = 0.0;
 459     }
 460     oc[mini - 1] += oc[minj - 1];
 461   }
 462   /* the last cycle */
 463   nude = 1;
 464   for (i = 1; i <= spp; i++) {
 465     if (cluster[i - 1] != NULL) {
 466       el[nude - 1] = i;
 467       nude++;
 468     }
 469   }
 470   if (!njoin) {
 471     curtree.start = cluster[el[0] - 1];
 472     curtree.start->back = NULL;
 473     free(av);
 474     free(oc);
 475     return;
 476   }
 477   bi = (x[el[0] - 1][el[1] - 1] + x[el[0] - 1][el[2] - 1] - x[el[1] - 1]
 478         [el[2] - 1]) * 0.5;
 479   bj = x[el[0] - 1][el[1] - 1] - bi;
 480   bk = x[el[0] - 1][el[2] - 1] - bi;
 481   bi -= av[el[0] - 1];
 482   bj -= av[el[1] - 1];
 483   bk -= av[el[2] - 1];
 484   if (progress) {
 485     printf("last cycle:\n");
 486     putchar(' ');
 487     nodelabel((boolean)(av[el[0] - 1] > 0.0));
 488     printf(" %ld  (%10.5f) joins ", el[0], bi);
 489     nodelabel((boolean)(av[el[1] - 1] > 0.0));
 490     printf(" %ld  (%10.5f) joins ", el[1], bj);
 491     nodelabel((boolean)(av[el[2] - 1] > 0.0));
 492     printf(" %ld  (%10.5f)\n", el[2], bk);
 493 #ifdef WIN32
 494     phyFillScreenColor();
 495 #endif
 496   }
 497   hookup(curtree.nodep[nextnode - 1], cluster[el[0] - 1]);
 498   hookup(curtree.nodep[nextnode - 1]->next, cluster[el[1] - 1]);
 499   hookup(curtree.nodep[nextnode - 1]->next->next, cluster[el[2] - 1]);
 500   cluster[el[0] - 1]->v = bi;
 501   cluster[el[1] - 1]->v = bj;
 502   cluster[el[2] - 1]->v = bk;
 503   cluster[el[0] - 1]->back->v = bi;
 504   cluster[el[1] - 1]->back->v = bj;
 505   cluster[el[2] - 1]->back->v = bk;
 506   curtree.start = cluster[el[0] - 1]->back;
 507   free(av);
 508   free(oc);
 509   free(R);
 510 }  /* jointree */
 511
 512
 513 void maketree()
 514 {
 515   /* construct the tree */
 516   long i ;
 517
 518   inputdata(replicates, printdata, lower, upper, x, reps);
 519   if (njoin && (spp < 3)) {
 520     printf("\nERROR: Neighbor-Joining runs must have at least 3 species\n\n");
 521     exxit(-1);
 522   }
 523   if (progress)
 524     putchar('\n');
 525   if (ith == 1)
 526     setuptree(&curtree, nonodes2 + 1);
 527   for (i = 1; i <= spp; i++)
 528     enterorder[i - 1] = i;
 529   if (jumble)
 530     randumize(seed, enterorder);
 531   for (i = 0; i < spp; i++)
 532     cluster[i] = curtree.nodep[i];
 533   jointree();
 534   if (njoin)
 535     curtree.start = curtree.nodep[outgrno - 1]->back;
 536   printree(curtree.start, treeprint, njoin, (boolean)(!njoin));
 537   if (treeprint)
 538     summarize();
 539   if (trout) {
 540     col = 0;
 541     if (njoin)
 542       treeout(curtree.start, &col, 0.43429448222, njoin, curtree.start);
 543     else
 544       curtree.root = curtree.start,
 545       treeoutr(curtree.start,&col,&curtree);
 546   }
 547   if (progress) {
 548     printf("\nOutput written on file \"%s\"\n\n", outfilename);
 549     if (trout)
 550       printf("Tree written on file \"%s\"\n\n", outtreename);
 551   }
 552 }  /* maketree */
 553
 554
 555 int main(int argc, Char *argv[])
 556 {  /* main program */
 557 #ifdef MAC
 558   argc = 1;                /* macsetup("Neighbor","");                */
 559   argv[0] = "Neighbor";
 560 #endif
 561   init(argc, argv);
 562   openfile(&infile,INFILE,"input file", "r",argv[0],infilename);
 563   openfile(&outfile,OUTFILE,"output file", "w",argv[0],outfilename);
 564   ibmpc = IBMCRT;
 565   ansi = ANSICRT;
 566   mulsets = false;
 567   datasets = 1;
 568   doinit();
 569   if (trout)
 570     openfile(&outtree,OUTTREE,"output tree file", "w",argv[0],outtreename);
 571   ith = 1;
 572   while (ith <= datasets) {
 573     if (datasets > 1) {
 574       fprintf(outfile, "Data set # %ld:\n",ith);
 575       if (progress)
 576         printf("Data set # %ld:\n",ith);
 577     }
 578     inputoptions();
 579     maketree();
 580     if (eoln(infile) && (ith < datasets))
 581       scan_eoln(infile);
 582     ith++;
 583   }
 584   FClose(infile);
 585   FClose(outfile);
 586   FClose(outtree);
 587   freerest();
 588 #ifdef MAC
 589   fixmacfile(outfilename);
 590   fixmacfile(outtreename);
 591 #endif
 592   printf("Done.\n\n");
 593 #ifdef WIN32
 594   phyRestoreConsoleAttributes();
 595 #endif
 596   return 0;
 597 }
 598
 599
 600
 601
 602