forester/archive/RIO/others/phylip_mod/src/fitch.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 Joseph Felsenstein, Akiko Fuseki, Sean Lamont, and Andrew Keeffe.
   7    Permission is granted to copy and use this program provided no fee is
   8    charged for it and provided that this copyright notice is not removed. */
   9
  10 #define zsmoothings     10    /* number of zero-branch correction iterations */
  11 #define epsilonf        0.000001   /* a very small but not too small number  */
  12 #define delta           0.0001      /* a not quite so small number */
  13 #define MAXNUMTREES   100000000 /* a number bigger than conceivable numtrees */
  14
  15
  16 #ifndef OLDC
  17 /* function prototypes */
  18 void   getoptions(void);
  19 void   allocrest(void);
  20 void   doinit(void);
  21 void   inputoptions(void);
  22 void   fitch_getinput(void);
  23 void   secondtraverse(node *, double , long *, double *);
  24 void   firsttraverse(node *, long *, double *);
  25 double evaluate(tree *);
  26 void   nudists(node *, node *);
  27 void   makedists(node *);
  28
  29 void   makebigv(node *);
  30 void   correctv(node *);
  31 void   alter(node *, node *);
  32 void   nuview(node *);
  33 void   update(node *);
  34 void   smooth(node *);
  35 void   filltraverse(node *, node *, boolean);
  36 void   fillin(node *, node *, boolean);
  37 void   insert_(node *, node *, boolean);
  38 void   copynode(node *, node *);
  39
  40 void   copy_(tree *, tree *);
  41 void   setuptipf(long, tree *);
  42 void   buildnewtip(long , tree *, long);
  43 void   buildsimpletree(tree *, long);
  44 void   addtraverse(node *, node *, boolean, long *, boolean *);
  45 void   re_move(node **, node **);
  46 void   rearrange(node *, long *, long *, boolean *);
  47 void   describe(node *);
  48 void   summarize(long);
  49 void   nodeinit(node *);
  50 void   initrav(node *);
  51 void   treevaluate(void);
  52 void   maketree(void);
  53 void   globrearrange(long* numtrees,boolean* succeeded);
  54 /* function prototypes */
  55 #endif
  56
  57
  58
  59 Char infilename[FNMLNGTH], outfilename[FNMLNGTH], intreename[FNMLNGTH], outtreename[FNMLNGTH];
  60 long nonodes2, outgrno, nums, col, datasets, ith, njumble, jumb=0;
  61 long inseed;
  62 vector *x;
  63 intvector *reps;
  64 boolean minev, global, jumble, lengths, usertree, lower, upper, negallowed,
  65         outgropt, replicates, trout, printdata, progress, treeprint,
  66         mulsets, firstset;
  67 double power;
  68 double trweight; /* to make treeread happy */
  69 boolean goteof, haslengths;  /* ditto ... */
  70 boolean first; /* ditto ... */
  71 node *addwhere;
  72
  73 longer seed;
  74 long *enterorder;
  75 tree curtree, priortree, bestree, bestree2;
  76 Char ch;
  77 char *progname;
  78
  79
  80
  81 void getoptions()
  82 {
  83   /* interactively set options */
  84   long inseed0=0, loopcount;
  85   Char ch;
  86   boolean done=false;
  87
  88   putchar('\n');
  89   minev = false;
  90   global = false;
  91   jumble = false;
  92   njumble = 1;
  93   lengths = false;
  94   lower = false;
  95   negallowed = false;
  96   outgrno = 1;
  97   outgropt = false;
  98   power = 2.0;
  99   replicates = false;
 100   trout = true;
 101   upper = false;
 102   usertree = false;
 103   printdata = false;
 104   progress = true;
 105   treeprint = true;
 106   loopcount = 0;
 107   do {
 108     cleerhome();
 109     printf("\nFitch-Margoliash method version %s\n\n",VERSION);
 110     printf("Settings for this run:\n");
 111     printf("  D      Method (F-M, Minimum Evolution)?  %s\n",
 112              (minev ? "Minimum Evolution" : "Fitch-Margoliash"));
 113     printf("  U                 Search for best tree?  %s\n",
 114            (usertree ? "No, use user trees in input file" : "Yes"));
 115     if (usertree) {
 116       printf("  N          Use lengths from user trees?  %s\n",
 117              (lengths ? "Yes" : "No"));
 118     }
 119     printf("  P                                Power?%9.5f\n",power);
 120     printf("  -      Negative branch lengths allowed?  %s\n",
 121            negallowed ? "Yes" : "No");
 122     printf("  O                        Outgroup root?");
 123     if (outgropt)
 124       printf("  Yes, at species number%3ld\n", outgrno);
 125     else
 126       printf("  No, use as outgroup species%3ld\n", outgrno);
 127     printf("  L         Lower-triangular data matrix?");
 128     if (lower)
 129       printf("  Yes\n");
 130     else
 131       printf("  No\n");
 132     printf("  R         Upper-triangular data matrix?");
 133     if (upper)
 134       printf("  Yes\n");
 135     else
 136       printf("  No\n");
 137     printf("  S                        Subreplicates?");
 138     if (replicates)
 139       printf("  Yes\n");
 140     else
 141       printf("  No\n");
 142     if (!usertree) {
 143       printf("  G                Global rearrangements?");
 144       if (global)
 145         printf("  Yes\n");
 146       else
 147         printf("  No\n");
 148       printf("  J     Randomize input order of species?");
 149       if (jumble)
 150         printf("  Yes (seed =%8ld,%3ld times)\n", inseed0, njumble);
 151       else
 152         printf("  No. Use input order\n");
 153     }
 154     printf("  M           Analyze multiple data sets?");
 155     if (mulsets)
 156       printf("  Yes, %2ld sets\n", datasets);
 157     else
 158       printf("  No\n");
 159     printf("  0   Terminal type (IBM PC, ANSI, none)?");
 160     if (ibmpc)
 161       printf("  IBM PC\n");
 162     if (ansi)
 163       printf("  ANSI\n");
 164     if (!(ibmpc || ansi))
 165       printf("  (none)\n");
 166     printf("  1    Print out the data at start of run");
 167     if (printdata)
 168       printf("  Yes\n");
 169     else
 170       printf("  No\n");
 171     printf("  2  Print indications of progress of run");
 172     if (progress)
 173       printf("  Yes\n");
 174     else
 175       printf("  No\n");
 176     printf("  3                        Print out tree");
 177     if (treeprint)
 178       printf("  Yes\n");
 179     else
 180       printf("  No\n");
 181     printf("  4       Write out trees onto tree file?");
 182     if (trout)
 183       printf("  Yes\n");
 184     else
 185       printf("  No\n");
 186     printf(
 187    "\n  Y to accept these or type the letter for one to change\n");
 188 #ifdef WIN32
 189     phyFillScreenColor();
 190 #endif
 191     scanf("%c%*[^\n]", &ch);
 192     getchar();
 193     uppercase(&ch);
 194     done = (ch == 'Y');
 195    if (!done) {
 196       if (strchr("DJOUNPG-LRSM01234",ch) != NULL) {
 197         switch (ch) {
 198
 199         case 'D':
 200           minev = !minev;
 201           if (minev && (!negallowed))
 202             negallowed = true;
 203           break;
 204
 205         case '-':
 206           negallowed = !negallowed;
 207           break;
 208
 209         case 'G':
 210           global = !global;
 211           break;
 212
 213         case 'J':
 214           jumble = !jumble;
 215           if (jumble)
 216             initjumble(&inseed, &inseed0, seed, &njumble);
 217           else njumble = 1;
 218           break;
 219
 220         case 'L':
 221           lower = !lower;
 222           break;
 223
 224          case 'N':
 225            lengths = !lengths;
 226            break;
 227
 228         case 'O':
 229           outgropt = !outgropt;
 230           if (outgropt)
 231             initoutgroup(&outgrno, spp);
 232           break;
 233
 234         case 'P':
 235           initpower(&power);
 236           break;
 237
 238         case 'R':
 239           upper = !upper;
 240           break;
 241
 242         case 'S':
 243           replicates = !replicates;
 244           break;
 245
 246         case 'U':
 247           usertree = !usertree;
 248           break;
 249
 250         case 'M':
 251           mulsets = !mulsets;
 252           if (mulsets)
 253             initdatasets(&datasets);
 254           jumble = true;
 255           if (jumble)
 256             initseed(&inseed, &inseed0, seed);
 257           break;
 258
 259         case '0':
 260           initterminal(&ibmpc, &ansi);
 261           break;
 262
 263         case '1':
 264           printdata = !printdata;
 265           break;
 266
 267         case '2':
 268           progress = !progress;
 269           break;
 270
 271         case '3':
 272           treeprint = !treeprint;
 273           break;
 274
 275         case '4':
 276           trout = !trout;
 277           break;
 278         }
 279       } else
 280         printf("Not a possible option!\n");
 281     }
 282     countup(&loopcount, 100);
 283   } while (!done);
 284   if (lower && upper) {
 285     printf("ERROR: Data matrix cannot be both uppeR and Lower triangular\n");
 286     exxit(-1);
 287   }
 288 }  /* getoptions */
 289
 290
 291 void allocrest()
 292 {
 293   long i;
 294
 295   x = (vector *)Malloc(spp*sizeof(vector));
 296   reps = (intvector *)Malloc(spp*sizeof(intvector));
 297   for (i=0;i<spp;++i){
 298     x[i]=(vector)Malloc(nonodes2 * sizeof(double));
 299     reps[i]=(intvector)Malloc(spp * sizeof(long));
 300   }
 301   nayme = (naym *)Malloc(spp*sizeof(naym));
 302   enterorder = (long *)Malloc(spp*sizeof(long));
 303 }
 304
 305
 306 void doinit()
 307 {
 308   /* initializes variables */
 309
 310   inputnumbers2(&spp, &nonodes2, 1);
 311   getoptions();
 312   if ( !usertree )
 313     nonodes2--;
 314   alloctree(&curtree.nodep, nonodes2);
 315   allocd(nonodes2, curtree.nodep);
 316   allocw(nonodes2, curtree.nodep);
 317   if (!usertree) {
 318     alloctree(&bestree.nodep, nonodes2);
 319     allocd(nonodes2, bestree.nodep);
 320     allocw(nonodes2, bestree.nodep);
 321     alloctree(&priortree.nodep, nonodes2);
 322     allocd(nonodes2, priortree.nodep);
 323     allocw(nonodes2, priortree.nodep);
 324     if (njumble > 1) {
 325       alloctree(&bestree2.nodep, nonodes2);
 326       allocd(nonodes2, bestree2.nodep);
 327       allocw(nonodes2, bestree2.nodep);
 328     }
 329   }
 330   allocrest();
 331 }  /* doinit */
 332
 333
 334 void inputoptions()
 335 {
 336   /* print options information */
 337   if (!firstset)
 338     samenumsp2(ith);
 339   fprintf(outfile, "\nFitch-Margoliash method version %s\n\n",VERSION);
 340   if (minev)
 341     fprintf(outfile, "Minimum evolution method option\n\n");
 342   fprintf(outfile, "                  __ __             2\n");
 343   fprintf(outfile, "                  \\  \\   (Obs - Exp)\n");
 344   fprintf(outfile, "Sum of squares =  /_ /_  ------------\n");
 345   fprintf(outfile, "                               ");
 346   if (power == (long)power)
 347     fprintf(outfile, "%2ld\n", (long)power);
 348   else
 349     fprintf(outfile, "%4.1f\n", power);
 350   fprintf(outfile, "                   i  j      Obs\n\n");
 351   fprintf(outfile, "Negative branch lengths ");
 352   if (!negallowed)
 353     fprintf(outfile, "not ");
 354   fprintf(outfile, "allowed\n\n");
 355   if (global)
 356     fprintf(outfile, "global optimization\n\n");
 357 }  /* inputoptions */
 358
 359
 360 void fitch_getinput()
 361 {
 362   /* reads the input data */
 363   inputoptions();
 364 }  /* fitch_getinput */
 365
 366
 367 void secondtraverse(node *q, double y, long *nx, double *sum)
 368 {
 369   /* from each of those places go back to all others */
 370    /* nx comes from firsttraverse */
 371    /* sum comes from evaluate via firsttraverse */
 372   double z=0.0, TEMP=0.0;
 373
 374   z = y + q->v;
 375   if (q->tip) {
 376     TEMP = q->d[(*nx) - 1] - z;
 377     *sum += q->w[(*nx) - 1] * (TEMP * TEMP);
 378   } else {
 379     secondtraverse(q->next->back, z, nx, sum);
 380     secondtraverse(q->next->next->back, z, nx,sum);
 381   }
 382 }  /* secondtraverse */
 383
 384
 385 void firsttraverse(node *p, long *nx, double *sum)
 386 {
 387   /* go through tree calculating branch lengths */
 388   if (minev && (p != curtree.start))
 389     *sum += p->v;
 390   if (p->tip) {
 391     if (!minev) {
 392       *nx = p->index;
 393       secondtraverse(p->back, 0.0, nx, sum);
 394       }
 395   } else {
 396     firsttraverse(p->next->back, nx,sum);
 397     firsttraverse(p->next->next->back, nx,sum);
 398   }
 399 }  /* firsttraverse */
 400
 401
 402 double evaluate(tree *t)
 403 {
 404   double sum=0.0;
 405   long nx=0;
 406   /* evaluate likelihood of a tree */
 407   firsttraverse(t->start->back ,&nx, &sum);
 408   firsttraverse(t->start, &nx, &sum);
 409   if ((!minev) && replicates && (lower || upper))
 410     sum /= 2;
 411   t->likelihood = -sum;
 412   return (-sum);
 413 }  /* evaluate */
 414
 415
 416 void nudists(node *x, node *y)
 417 {
 418   /* compute distance between an interior node and tips */
 419   long nq=0, nr=0, nx=0, ny=0;
 420   double dil=0, djl=0, wil=0, wjl=0, vi=0, vj=0;
 421   node *qprime, *rprime;
 422
 423   qprime = x->next;
 424   rprime = qprime->next->back;
 425   qprime = qprime->back;
 426   ny = y->index;
 427   dil = qprime->d[ny - 1];
 428   djl = rprime->d[ny - 1];
 429   wil = qprime->w[ny - 1];
 430   wjl = rprime->w[ny - 1];
 431   vi = qprime->v;
 432   vj = rprime->v;
 433   x->w[ny - 1] = wil + wjl;
 434   if (wil + wjl <= 0.0)
 435     x->d[ny - 1] = 0.0;
 436   else
 437     x->d[ny - 1] = ((dil - vi) * wil + (djl - vj) * wjl) / (wil + wjl);
 438   nx = x->index;
 439   nq = qprime->index;
 440   nr = rprime->index;
 441   dil = y->d[nq - 1];
 442   djl = y->d[nr - 1];
 443   wil = y->w[nq - 1];
 444   wjl = y->w[nr - 1];
 445   y->w[nx - 1] = wil + wjl;
 446   if (wil + wjl <= 0.0)
 447     y->d[nx - 1] = 0.0;
 448   else
 449     y->d[nx - 1] = ((dil - vi) * wil + (djl - vj) * wjl) / (wil + wjl);
 450 }  /* nudists */
 451
 452
 453 void makedists(node *p)
 454 {
 455   /* compute distances among three neighbors of a node */
 456   long i=0, nr=0, ns=0;
 457   node *q, *r, *s;
 458
 459   r = p->back;
 460   nr = r->index;
 461   for (i = 1; i <= 3; i++) {
 462     q = p->next;
 463     s = q->back;
 464     ns = s->index;
 465     if (s->w[nr - 1] + r->w[ns - 1] <= 0.0)
 466       p->dist = 0.0;
 467     else
 468       p->dist = (s->w[nr - 1] * s->d[nr - 1] + r->w[ns - 1] * r->d[ns - 1]) /
 469                 (s->w[nr - 1] + r->w[ns - 1]);
 470     p = q;
 471     r = s;
 472     nr = ns;
 473   }
 474 }  /* makedists */
 475
 476
 477 void makebigv(node *p)
 478 {
 479   /* make new branch length */
 480   long i=0;
 481   node *temp, *q, *r;
 482
 483   q = p->next;
 484   r = q->next;
 485   for (i = 1; i <= 3; i++) {
 486     if (p->iter) {
 487       p->v = (p->dist + r->dist - q->dist) / 2.0;
 488       p->back->v = p->v;
 489     }
 490     temp = p;
 491     p = q;
 492     q = r;
 493     r = temp;
 494   }
 495 }  /* makebigv */
 496
 497
 498 void correctv(node *p)
 499 {
 500   /* iterate branch lengths if some are to be zero */
 501   node *q, *r, *temp;
 502   long i=0, j=0, n=0, nq=0, nr=0, ntemp=0;
 503   double wq=0.0, wr=0.0;
 504
 505   q = p->next;
 506   r = q->next;
 507   n = p->back->index;
 508   nq = q->back->index;
 509   nr = r->back->index;
 510   for (i = 1; i <= zsmoothings; i++) {
 511     for (j = 1; j <= 3; j++) {
 512       if (p->iter) {
 513         wr = r->back->w[n - 1] + p->back->w[nr - 1];
 514         wq = q->back->w[n - 1] + p->back->w[nq - 1];
 515         if (wr + wq <= 0.0 && !negallowed)
 516           p->v = 0.0;
 517         else
 518           p->v = ((p->dist - q->v) * wq + (r->dist - r->v) * wr) / (wr + wq);
 519         if (p->v < 0 && !negallowed)
 520           p->v = 0.0;
 521         p->back->v = p->v;
 522       }
 523       temp = p;
 524       p = q;
 525       q = r;
 526       r = temp;
 527       ntemp = n;
 528       n = nq;
 529       nq = nr;
 530       nr = ntemp;
 531     }
 532   }
 533 }  /* correctv */
 534
 535
 536 void alter(node *x, node *y)
 537 {
 538   /* traverse updating these views */
 539   nudists(x, y);
 540   if (!y->tip) {
 541     alter(x, y->next->back);
 542     alter(x, y->next->next->back);
 543   }
 544 }  /* alter */
 545
 546
 547 void nuview(node *p)
 548 {
 549   /* renew information about subtrees */
 550   long i=0;
 551   node *q, *r, *pprime, *temp;
 552
 553   q = p->next;
 554   r = q->next;
 555   for (i = 1; i <= 3; i++) {
 556     temp = p;
 557     pprime = p->back;
 558     alter(p, pprime);
 559     p = q;
 560     q = r;
 561     r = temp;
 562   }
 563 }  /* nuview */
 564
 565
 566 void update(node *p)
 567 {
 568   /* update branch lengths around a node */
 569
 570   if (p->tip)
 571     return;
 572   makedists(p);
 573   if (p->iter || p->next->iter || p->next->next->iter) {
 574     makebigv(p);
 575     correctv(p);
 576   }
 577   nuview(p);
 578 }  /* update */
 579
 580
 581 void smooth(node *p)
 582 {
 583   /* go through tree getting new branch lengths and views */
 584   if (p->tip)
 585     return;
 586   update(p);
 587   smooth(p->next->back);
 588   smooth(p->next->next->back);
 589 }  /* smooth */
 590
 591
 592 void filltraverse(node *pb, node *qb, boolean contin)
 593 {
 594   if (qb->tip)
 595     return;
 596   if (contin) {
 597     filltraverse(pb, qb->next->back,contin);
 598     filltraverse(pb, qb->next->next->back,contin);
 599     nudists(qb, pb);
 600     return;
 601   }
 602   if (!qb->next->back->tip)
 603     nudists(qb->next->back, pb);
 604   if (!qb->next->next->back->tip)
 605     nudists(qb->next->next->back, pb);
 606 }  /* filltraverse */
 607
 608
 609 void fillin(node *pa, node *qa, boolean contin)
 610 {
 611   if (!pa->tip) {
 612     fillin(pa->next->back, qa, contin);
 613     fillin(pa->next->next->back, qa, contin);
 614   }
 615   filltraverse(pa, qa, contin);
 616 }  /* fillin */
 617
 618
 619 void insert_(node *p, node *q, boolean contin_)
 620 {
 621   /* put p and q together and iterate info. on resulting tree */
 622   double x=0.0, oldlike;
 623   hookup(p->next->next, q->back);
 624   hookup(p->next, q);
 625   x = q->v / 2.0;
 626   p->v = 0.0;
 627   p->back->v = 0.0;
 628   p->next->v = x;
 629   p->next->back->v = x;
 630   p->next->next->back->v = x;
 631   p->next->next->v = x;
 632   fillin(p->back, p, contin_);
 633   evaluate(&curtree);
 634   do {
 635     oldlike = curtree.likelihood;
 636     smooth(p);
 637     smooth(p->back);
 638     evaluate(&curtree);
 639   } while (fabs(curtree.likelihood - oldlike) > delta);
 640 }  /* insert_ */
 641
 642
 643 void copynode(node *c, node *d)
 644 {
 645   /* make a copy of a node */
 646
 647   memcpy(d->d, c->d, nonodes2*sizeof(double));
 648   memcpy(d->w, c->w, nonodes2*sizeof(double));
 649   d->v = c->v;
 650   d->iter = c->iter;
 651   d->dist = c->dist;
 652   d->xcoord = c->xcoord;
 653   d->ycoord = c->ycoord;
 654   d->ymin = c->ymin;
 655   d->ymax = c->ymax;
 656 }  /* copynode */
 657
 658
 659 void copy_(tree *a, tree *b)
 660 {
 661   /* make copy of a tree a to tree b */
 662   long i, j=0;
 663   node *p, *q;
 664
 665   for (i = 0; i < spp; i++) {
 666     copynode(a->nodep[i], b->nodep[i]);
 667     if (a->nodep[i]->back) {
 668       if (a->nodep[i]->back == a->nodep[a->nodep[i]->back->index - 1])
 669         b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1];
 670       else if (a->nodep[i]->back
 671                  == a->nodep[a->nodep[i]->back->index - 1]->next)
 672         b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1]->next;
 673       else
 674         b->nodep[i]->back
 675           = b->nodep[a->nodep[i]->back->index - 1]->next->next;
 676     }
 677     else b->nodep[i]->back = NULL;
 678   }
 679   for (i = spp; i < nonodes2; i++) {
 680     p = a->nodep[i];
 681     q = b->nodep[i];
 682     for (j = 1; j <= 3; j++) {
 683       copynode(p, q);
 684       if (p->back) {
 685         if (p->back == a->nodep[p->back->index - 1])
 686           q->back = b->nodep[p->back->index - 1];
 687         else if (p->back == a->nodep[p->back->index - 1]->next)
 688           q->back = b->nodep[p->back->index - 1]->next;
 689         else
 690           q->back = b->nodep[p->back->index - 1]->next->next;
 691       }
 692       else
 693         q->back = NULL;
 694       p = p->next;
 695       q = q->next;
 696     }
 697   }
 698   b->likelihood = a->likelihood;
 699   b->start = a->start;
 700 }  /* copy_ */
 701
 702
 703 void setuptipf(long m, tree *t)
 704 {
 705   /* initialize branch lengths and views in a tip */
 706   long i=0;
 707   intvector n=(long *)Malloc(spp * sizeof(long));
 708   node *WITH;
 709
 710   WITH = t->nodep[m - 1];
 711   memcpy(WITH->d, x[m - 1], (nonodes2 * sizeof(double)));
 712   memcpy(n, reps[m - 1], (spp * sizeof(long)));
 713   for (i = 0; i < spp; i++) {
 714     if (i + 1 != m && n[i] > 0) {
 715       if (WITH->d[i] < epsilonf)
 716         WITH->d[i] = epsilonf;
 717       WITH->w[i] = n[i] / exp(power * log(WITH->d[i]));
 718     } else {
 719       WITH->w[i] = 0.0;
 720       WITH->d[i] = 0.0;
 721     }
 722   }
 723   for (i = spp; i < nonodes2; i++) {
 724     WITH->w[i] = 1.0;
 725     WITH->d[i] = 0.0;
 726   }
 727   WITH->index = m;
 728   if (WITH->iter) WITH->v = 0.0;
 729   free(n);
 730 }  /* setuptipf */
 731
 732
 733 void buildnewtip(long m, tree *t, long nextsp)
 734 {
 735   /* initialize and hook up a new tip */
 736   node *p;
 737   setuptipf(m, t);
 738   p = t->nodep[nextsp + spp - 3];
 739   hookup(t->nodep[m - 1], p);
 740 }  /* buildnewtip */
 741
 742
 743 void buildsimpletree(tree *t, long nextsp)
 744 {
 745   /* make and initialize a three-species tree */
 746   curtree.start=curtree.nodep[enterorder[0] - 1];
 747   setuptipf(enterorder[0], t);
 748   setuptipf(enterorder[1], t);
 749   hookup(t->nodep[enterorder[0] - 1], t->nodep[enterorder[1] - 1]);
 750   buildnewtip(enterorder[2], t, nextsp);
 751   insert_(t->nodep[enterorder[2] - 1]->back, t->nodep[enterorder[0] - 1],
 752           false);
 753 }  /* buildsimpletree */
 754
 755
 756 void addtraverse(node *p, node *q, boolean contin, long *numtrees,
 757                         boolean *succeeded)
 758 {
 759  /* traverse through a tree, finding best place to add p */
 760   insert_(p, q, true);
 761   (*numtrees)++;
 762   if (evaluate(&curtree) > (bestree.likelihood +
 763                             epsilonf * fabs(bestree.likelihood))){
 764     copy_(&curtree, &bestree);
 765     addwhere = q;
 766     (*succeeded)=true;
 767   }
 768   copy_(&priortree, &curtree);
 769   if (!q->tip && contin) {
 770     addtraverse(p, q->next->back, contin,numtrees,succeeded);
 771     addtraverse(p, q->next->next->back, contin,numtrees,succeeded);
 772   }
 773 }  /* addtraverse */
 774
 775
 776 void re_move(node **p, node **q)
 777 {
 778   /* re_move p and record in q where it was */
 779   *q = (*p)->next->back;
 780   hookup(*q, (*p)->next->next->back);
 781   (*p)->next->back = NULL;
 782   (*p)->next->next->back = NULL;
 783   update(*q);
 784   update((*q)->back);
 785 }  /* re_move */
 786
 787
 788 void globrearrange(long* numtrees,boolean* succeeded)
 789 {
 790   /* does global rearrangements */
 791   tree globtree;
 792   tree oldtree;
 793   int i,j,k,num_sibs,num_sibs2;
 794   node *where,*sib_ptr,*sib_ptr2;
 795   double oldbestyet = curtree.likelihood;
 796   int success = false;
 797
 798   alloctree(&globtree.nodep,nonodes2);
 799   alloctree(&oldtree.nodep,nonodes2);
 800   setuptree(&globtree,nonodes2);
 801   setuptree(&oldtree,nonodes2);
 802   allocd(nonodes2, globtree.nodep);
 803   allocd(nonodes2, oldtree.nodep);
 804   allocw(nonodes2, globtree.nodep);
 805   allocw(nonodes2, oldtree.nodep);
 806   copy_(&curtree,&globtree);
 807   copy_(&curtree,&oldtree);
 808   for ( i = spp ; i < nonodes2 ; i++ ) {
 809     num_sibs = count_sibs(curtree.nodep[i]);
 810     sib_ptr  = curtree.nodep[i];
 811     if ( (i - spp) % (( nonodes2 / 72 ) + 1 ) == 0 )
 812       putchar('.');
 813     fflush(stdout);
 814     for ( j = 0 ; j <= num_sibs ; j++ ) {
 815       re_move(&sib_ptr,&where);
 816       copy_(&curtree,&priortree);
 817
 818       if (where->tip) {
 819         copy_(&oldtree,&curtree);
 820         copy_(&oldtree,&bestree);
 821         sib_ptr=sib_ptr->next;
 822         continue;
 823       }
 824       else num_sibs2 = count_sibs(where);
 825       sib_ptr2 = where;
 826       for ( k = 0 ; k < num_sibs2 ; k++ ) {
 827         addwhere = NULL;
 828         addtraverse(sib_ptr,sib_ptr2->back,true,numtrees,succeeded);
 829         if ( addwhere && where != addwhere && where->back != addwhere
 830               && bestree.likelihood > globtree.likelihood) {
 831             copy_(&bestree,&globtree);
 832             success = true;
 833         }
 834         sib_ptr2 = sib_ptr2->next;
 835       }
 836       copy_(&oldtree,&curtree);
 837       copy_(&oldtree,&bestree);
 838       sib_ptr = sib_ptr->next;
 839     }
 840   }
 841   copy_(&globtree,&curtree);
 842   copy_(&globtree,&bestree);
 843   if (success && globtree.likelihood > oldbestyet)  {
 844     *succeeded = true;
 845   }
 846   else  {
 847     *succeeded = false;
 848   }
 849   freed(nonodes2, globtree.nodep);
 850   freed(nonodes2, oldtree.nodep);
 851   freew(nonodes2, globtree.nodep);
 852   freew(nonodes2, oldtree.nodep);
 853   freetree(&globtree.nodep,nonodes2);
 854   freetree(&oldtree.nodep,nonodes2);
 855 }
 856
 857
 858 void rearrange(node *p, long *numtrees, long *nextsp, boolean *succeeded)
 859 {
 860   node *q, *r;
 861   if (!p->tip && !p->back->tip) {
 862     r = p->next->next;
 863     re_move(&r, &q);
 864     copy_(&curtree, &priortree);
 865     addtraverse(r, q->next->back, false, numtrees,succeeded);
 866     addtraverse(r, q->next->next->back, false, numtrees,succeeded);
 867     copy_(&bestree, &curtree);
 868     if (global && ((*nextsp) == spp)) {
 869       putchar('.');
 870       fflush(stdout);
 871     }
 872   }
 873   if (!p->tip) {
 874     rearrange(p->next->back, numtrees,nextsp,succeeded);
 875     rearrange(p->next->next->back, numtrees,nextsp,succeeded);
 876   }
 877 }  /* rearrange */
 878
 879
 880 void describe(node *p)
 881 {
 882   /* print out information for one branch */
 883   long i=0;
 884   node *q;
 885
 886   q = p->back;
 887   fprintf(outfile, "%4ld          ", q->index - spp);
 888   if (p->tip) {
 889     for (i = 0; i < nmlngth; i++)
 890       putc(nayme[p->index - 1][i], outfile);
 891   } else
 892     fprintf(outfile, "%4ld      ", p->index - spp);
 893   fprintf(outfile, "%15.5f\n", q->v);
 894   if (!p->tip) {
 895     describe(p->next->back);
 896     describe(p->next->next->back);
 897   }
 898 }  /* describe */
 899
 900
 901 void summarize(long numtrees)
 902 {
 903   /* print out branch lengths etc. */
 904   long i, j, totalnum;
 905
 906   fprintf(outfile, "\nremember:");
 907   if (outgropt)
 908     fprintf(outfile, " (although rooted by outgroup)");
 909   fprintf(outfile, " this is an unrooted tree!\n\n");
 910   if (!minev)
 911     fprintf(outfile, "Sum of squares = %11.5f\n\n", -curtree.likelihood);
 912   else
 913     fprintf(outfile, "Sum of branch lengths = %11.5f\n\n", -curtree.likelihood);
 914   if ((power == 2.0) && !minev) {
 915     totalnum = 0;
 916     for (i = 1; i <= nums; i++) {
 917       for (j = 1; j <= nums; j++) {
 918         if (i != j)
 919           totalnum += reps[i - 1][j - 1];
 920       }
 921     }
 922     fprintf(outfile, "Average percent standard deviation = ");
 923     fprintf(outfile, "%11.5f\n\n",
 924             100 * sqrt(-curtree.likelihood / (totalnum - 2)));
 925   }
 926   fprintf(outfile, "Between        And            Length\n");
 927   fprintf(outfile, "-------        ---            ------\n");
 928   describe(curtree.start->next->back);
 929   describe(curtree.start->next->next->back);
 930   describe(curtree.start->back);
 931   fprintf(outfile, "\n\n");
 932   if (trout) {
 933     col = 0;
 934     treeout(curtree.start, &col, 0.43429445222, true,
 935               curtree.start);
 936   }
 937 }  /* summarize */
 938
 939
 940 void nodeinit(node *p)
 941 {
 942   /* initialize a node */
 943   long i, j;
 944
 945   for (i = 1; i <= 3; i++) {
 946     for (j = 0; j < nonodes2; j++) {
 947       p->w[j] = 1.0;
 948       p->d[j] = 0.0;
 949     }
 950     p = p->next;
 951   }
 952   if ((!lengths) || p->iter)
 953     p->v = 1.0;
 954   if ((!lengths) || p->back->iter)
 955     p->back->v = 1.0;
 956 }  /* nodeinit */
 957
 958
 959 void initrav(node *p)
 960 {
 961   /* traverse to initialize */
 962   if (p->tip)
 963     return;
 964   nodeinit(p);
 965   initrav(p->next->back);
 966   initrav(p->next->next->back);
 967 }  /* initrav */
 968
 969 void treevaluate()
 970 {
 971   /* evaluate user-defined tree, iterating branch lengths */
 972   long i;
 973   double oldlike;
 974
 975   for (i = 1; i <= spp; i++)
 976     setuptipf(i, &curtree);
 977   unroot(&curtree,nonodes2);
 978
 979   initrav(curtree.start);
 980   if (curtree.start->back != NULL) {
 981     initrav(curtree.start->back);
 982     evaluate(&curtree);
 983     do {
 984       oldlike = curtree.likelihood;
 985       smooth(curtree.start);
 986       evaluate(&curtree);
 987     } while (fabs(curtree.likelihood - oldlike) > delta);
 988   }
 989   evaluate(&curtree);
 990 }  /* treevaluate */
 991
 992
 993 void maketree()
 994 {
 995   /* contruct the tree */
 996   long nextsp,numtrees;
 997   boolean succeeded=false;
 998   long i, j, which;
 999
1000   if (usertree) {
1001     inputdata(replicates, printdata, lower, upper, x, reps);
1002     setuptree(&curtree, nonodes2);
1003     for (which = 1; which <= spp; which++)
1004       setuptipf(which, &curtree);
1005     if (eoln(infile))
1006       scan_eoln(infile);
1007     openfile(&intree,INTREE,"input tree file","r",progname,intreename);
1008     numtrees = countsemic(&intree);
1009     if (numtrees > MAXNUMTREES) {
1010       printf("\nERROR: number of input trees is read incorrectly from %s\n",
1011         intreename);
1012       exxit(-1);
1013     }
1014     if (treeprint) {
1015       fprintf(outfile, "User-defined tree");
1016       if (numtrees > 1)
1017         putc('s', outfile);
1018       fprintf(outfile, ":\n\n");
1019     }
1020     first = true;
1021     which = 1;
1022     while (which <= numtrees) {
1023       treeread2 (intree, &curtree.start, curtree.nodep,
1024         lengths, &trweight, &goteof, &haslengths, &spp,false,nonodes2);
1025       nums = spp;
1026       curtree.start = curtree.nodep[outgrno - 1]->back;
1027       treevaluate();
1028       printree(curtree.start, treeprint, false, false);
1029       summarize(numtrees);
1030       clear_connections(&curtree,nonodes2);
1031       which++;
1032     }
1033     FClose(intree);
1034   } else {
1035     if (jumb == 1) {
1036       inputdata(replicates, printdata, lower, upper, x, reps);
1037       setuptree(&curtree, nonodes2);
1038       setuptree(&priortree, nonodes2);
1039       setuptree(&bestree, nonodes2);
1040       if (njumble > 1) setuptree(&bestree2, nonodes2);
1041     }
1042     for (i = 1; i <= spp; i++)
1043       enterorder[i - 1] = i;
1044     if (jumble)
1045       randumize(seed, enterorder);
1046     nextsp = 3;
1047     buildsimpletree(&curtree, nextsp);
1048     curtree.start = curtree.nodep[enterorder[0] - 1]->back;
1049     if (jumb == 1) numtrees = 1;
1050     nextsp = 4;
1051     if (progress) {
1052       printf("Adding species:\n");
1053       writename(0, 3, enterorder);
1054 #ifdef WIN32
1055       phyFillScreenColor();
1056 #endif
1057     }
1058     while (nextsp <= spp) {
1059       nums = nextsp;
1060       buildnewtip(enterorder[nextsp - 1], &curtree, nextsp);
1061       copy_(&curtree, &priortree);
1062       bestree.likelihood = -99999.0;
1063       curtree.start = curtree.nodep[enterorder[0] - 1]->back;
1064       addtraverse(curtree.nodep[enterorder[nextsp - 1] - 1]->back,
1065                   curtree.start, true, &numtrees,&succeeded);
1066       copy_(&bestree, &curtree);
1067       if (progress) {
1068         writename(nextsp  - 1, 1, enterorder);
1069 #ifdef WIN32
1070         phyFillScreenColor();
1071 #endif
1072       }
1073       if (global && nextsp == spp) {
1074         if (progress) {
1075           printf("Doing global rearrangements\n");
1076           printf("  !");
1077           for (j = spp; j < nonodes2; j++)
1078             if ( (j - spp) % (( nonodes2 / 72 ) + 1 ) == 0 )
1079               putchar('-');
1080           printf("!\n");
1081           printf("   ");
1082         }
1083       }
1084       succeeded = true;
1085       while (succeeded) {
1086         succeeded = false;
1087         curtree.start = curtree.nodep[enterorder[0] - 1]->back;
1088         if (nextsp == spp  && global)
1089           globrearrange (&numtrees,&succeeded);
1090         else{
1091           rearrange(curtree.start,&numtrees,&nextsp,&succeeded);
1092         }
1093         if (global && ((nextsp) == spp) && progress)
1094           printf("\n   ");
1095       }
1096       if (global && nextsp == spp) {
1097         putc('\n', outfile);
1098         if (progress)
1099           putchar('\n');
1100       }
1101       if (njumble > 1) {
1102         if (jumb == 1 && nextsp == spp)
1103           copy_(&bestree, &bestree2);
1104         else if (nextsp == spp) {
1105           if (bestree2.likelihood < bestree.likelihood)
1106             copy_(&bestree, &bestree2);
1107         }
1108       }
1109       if (nextsp == spp && jumb == njumble) {
1110         if (njumble > 1) copy_(&bestree2, &curtree);
1111         curtree.start = curtree.nodep[outgrno - 1]->back;
1112         printree(curtree.start, treeprint, true, false);
1113         summarize(numtrees);
1114       }
1115       nextsp++;
1116     }
1117   }
1118   if (jumb == njumble && progress) {
1119     printf("\nOutput written to file \"%s\"\n\n", outfilename);
1120     if (trout) {
1121       printf("Tree also written onto file \"%s\"\n", outtreename);
1122       putchar('\n');
1123     }
1124   }
1125 }  /* maketree */
1126
1127
1128 int main(int argc, Char *argv[])
1129 {
1130   int i;
1131 #ifdef MAC
1132   argc = 1;                /* macsetup("Fitch","");        */
1133   argv[0]="Fitch";
1134 #endif
1135   init(argc,argv);
1136   progname = argv[0];
1137   openfile(&infile,INFILE,"input file","r",argv[0],infilename);
1138   openfile(&outfile,OUTFILE,"output file","w",argv[0],outfilename);
1139
1140   ibmpc = IBMCRT;
1141   ansi = ANSICRT;
1142   mulsets = false;
1143   datasets = 1;
1144   firstset = true;
1145   doinit();
1146   if (trout)
1147     openfile(&outtree,OUTTREE,"output tree file","w",argv[0],outtreename);
1148   for (i=0;i<spp;++i){
1149     enterorder[i]=0;}
1150   for (ith = 1; ith <= datasets; ith++) {
1151     if (datasets > 1) {
1152       fprintf(outfile, "Data set # %ld:\n\n",ith);
1153       if (progress)
1154         printf("\nData set # %ld:\n\n",ith);
1155     }
1156     fitch_getinput();
1157     for (jumb = 1; jumb <= njumble; jumb++)
1158         maketree();
1159     firstset = false;
1160     if (eoln(infile) && (ith < datasets))
1161       scan_eoln(infile);
1162   }
1163   if (trout)
1164     FClose(outtree);
1165   FClose(outfile);
1166   FClose(infile);
1167 #ifdef MAC
1168   fixmacfile(outfilename);
1169   fixmacfile(outtreename);
1170 #endif
1171   printf("Done.\n\n");
1172 #ifdef WIN32
1173   phyRestoreConsoleAttributes();
1174 #endif
1175   return 0;
1176 }