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

   1
   2 #include "phylip.h"
   3 #include "seq.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 long nonodes, endsite, outgrno, nextree, which;
  11 boolean interleaved, printdata, outgropt, treeprint, dotdiff, transvp;
  12 steptr weight, category, alias, location, ally;
  13 sequence y;
  14
  15
  16 void fix_x(node* p,long site, double maxx, long rcategs)
  17 { /* dnaml dnamlk */
  18   long i,j;
  19   p->underflows[site] += log(maxx);
  20
  21   for ( i = 0 ; i < rcategs ; i++ ) {
  22     for ( j = 0 ; j < ((long)T - (long)A + 1) ; j++)
  23       p->x[site][i][j] /= maxx;
  24   }
  25 } /* fix_x */
  26
  27
  28 void fix_protx(node* p,long site, double maxx, long rcategs)
  29 { /* proml promlk */
  30   long i,m;
  31
  32   p->underflows[site] += log(maxx);
  33
  34   for ( i = 0 ; i < rcategs  ; i++ )
  35     for (m = 0; m <= 19; m++)
  36       p->protx[site][i][m] /= maxx;
  37 } /* fix_protx */
  38
  39
  40 void alloctemp(node **temp, long *zeros, long endsite)
  41 {
  42   /*used in dnacomp and dnapenny */
  43   *temp = (node *)Malloc(sizeof(node));
  44   (*temp)->numsteps = (steptr)Malloc(endsite*sizeof(long));
  45   (*temp)->base = (baseptr)Malloc(endsite*sizeof(long));
  46   (*temp)->numnuc = (nucarray *)Malloc(endsite*sizeof(nucarray));
  47   memcpy((*temp)->base, zeros, endsite*sizeof(long));
  48   memcpy((*temp)->numsteps, zeros, endsite*sizeof(long));
  49   zeronumnuc(*temp, endsite);
  50 }  /* alloctemp */
  51
  52
  53 void freetemp(node **temp)
  54 {
  55   /* used in dnacomp, dnapars, & dnapenny */
  56   free((*temp)->numsteps);
  57   free((*temp)->base);
  58   free((*temp)->numnuc);
  59   free(*temp);
  60 }  /* freetemp */
  61
  62
  63 void freetree2 (pointarray treenode, long nonodes)
  64 {
  65   /* The natural complement to alloctree2.  Free all elements of all
  66   the rings (normally triads) in treenode */
  67   long i;
  68   node *p, *q;
  69
  70   /* The first spp elements are just nodes, not rings */
  71   for (i = 0; i < spp; i++)
  72     free (treenode[i]);
  73
  74   /* The rest are rings */
  75   for (i = spp; i < nonodes; i++) {
  76     p = treenode[i]->next;
  77     while (p != treenode[i]) {
  78       q = p->next;
  79       free (p);
  80       p = q;
  81     }
  82     /* p should now point to treenode[i], which has yet to be freed */
  83     free (p);
  84   }
  85   free (treenode);
  86 }  /* freetree2 */
  87
  88
  89 void inputdata(long chars)
  90 {
  91   /* input the names and sequences for each species */
  92   /* used by dnacomp, dnadist, dnainvar, dnaml, dnamlk, dnapars, & dnapenny */
  93   long i, j, k, l, basesread, basesnew=0;
  94   Char charstate;
  95   boolean allread, done;
  96
  97   if (printdata)
  98     headings(chars, "Sequences", "---------");
  99   basesread = 0;
 100   allread = false;
 101   while (!(allread)) {
 102     /* eat white space -- if the separator line has spaces on it*/
 103     do {
 104       charstate = gettc(infile);
 105     } while (charstate == ' ' || charstate == '\t');
 106     ungetc(charstate, infile);
 107     if (eoln(infile))
 108       scan_eoln(infile);
 109     i = 1;
 110     while (i <= spp) {
 111       if ((interleaved && basesread == 0) || !interleaved)
 112         initname(i-1);
 113       j = (interleaved) ? basesread : 0;
 114       done = false;
 115       while (!done && !eoff(infile)) {
 116         if (interleaved)
 117           done = true;
 118         while (j < chars && !(eoln(infile) || eoff(infile))) {
 119           charstate = gettc(infile);
 120           if (charstate == '\n' || charstate == '\t')
 121             charstate = ' ';
 122           if (charstate == ' ' || (charstate >= '0' && charstate <= '9'))
 123             continue;
 124           uppercase(&charstate);
 125           if ((strchr("ABCDGHKMNRSTUVWXY?O-",charstate)) == NULL){
 126             printf("ERROR: bad base: %c at site %5ld of species %3ld\n",
 127                    charstate, j+1, i);
 128             if (charstate == '.') {
 129               printf("       Periods (.) may not be used as gap characters.\n");
 130               printf("       The correct gap character is (-)\n");
 131             }
 132             exxit(-1);
 133           }
 134           j++;
 135           y[i - 1][j - 1] = charstate;
 136         }
 137         if (interleaved)
 138           continue;
 139         if (j < chars)
 140           scan_eoln(infile);
 141         else if (j == chars)
 142           done = true;
 143       }
 144       if (interleaved && i == 1)
 145         basesnew = j;
 146
 147       scan_eoln(infile);
 148
 149       if ((interleaved && j != basesnew) ||
 150           (!interleaved && j != chars)) {
 151         printf("\nERROR: sequences out of alignment at position %ld", j+1);
 152         printf(" of species %ld\n\n", i);
 153         exxit(-1);
 154       }
 155       i++;
 156     }
 157
 158     if (interleaved) {
 159       basesread = basesnew;
 160       allread = (basesread == chars);
 161     } else
 162       allread = (i > spp);
 163   }
 164   if (!printdata)
 165     return;
 166   for (i = 1; i <= ((chars - 1) / 60 + 1); i++) {
 167     for (j = 1; j <= spp; j++) {
 168       for (k = 0; k < nmlngth; k++)
 169         putc(nayme[j - 1][k], outfile);
 170       fprintf(outfile, "   ");
 171       l = i * 60;
 172       if (l > chars)
 173         l = chars;
 174       for (k = (i - 1) * 60 + 1; k <= l; k++) {
 175         if (dotdiff && (j > 1 && y[j - 1][k - 1] == y[0][k - 1]))
 176           charstate = '.';
 177         else
 178           charstate = y[j - 1][k - 1];
 179         putc(charstate, outfile);
 180         if (k % 10 == 0 && k % 60 != 0)
 181           putc(' ', outfile);
 182       }
 183       putc('\n', outfile);
 184     }
 185     putc('\n', outfile);
 186   }
 187   putc('\n', outfile);
 188 }  /* inputdata */
 189
 190
 191 void alloctree(pointarray *treenode, long nonodes, boolean usertree)
 192 {
 193   /* allocate treenode dynamically */
 194   /* used in dnapars, dnacomp, dnapenny & dnamove */
 195   long i, j;
 196   node *p, *q;
 197
 198   *treenode = (pointarray)Malloc(nonodes*sizeof(node *));
 199   for (i = 0; i < spp; i++) {
 200     (*treenode)[i] = (node *)Malloc(sizeof(node));
 201     (*treenode)[i]->tip = true;
 202     (*treenode)[i]->index = i+1;
 203     (*treenode)[i]->iter = true;
 204     (*treenode)[i]->branchnum = 0;
 205     (*treenode)[i]->initialized = true;
 206   }
 207   if (!usertree)
 208     for (i = spp; i < nonodes; i++) {
 209       q = NULL;
 210       for (j = 1; j <= 3; j++) {
 211         p = (node *)Malloc(sizeof(node));
 212         p->tip = false;
 213         p->index = i+1;
 214         p->iter = true;
 215         p->branchnum = 0;
 216         p->initialized = false;
 217         p->next = q;
 218         q = p;
 219       }
 220       p->next->next->next = p;
 221       (*treenode)[i] = p;
 222     }
 223 } /* alloctree */
 224
 225
 226 void allocx(long nonodes, long rcategs, pointarray treenode, boolean usertree)
 227 {
 228   /* allocate x dynamically */
 229   /* used in dnaml & dnamlk */
 230   long i, j, k;
 231   node *p;
 232
 233   for (i = 0; i < spp; i++){
 234     treenode[i]->x = (phenotype)Malloc(endsite*sizeof(ratelike));
 235     treenode[i]->underflows = Malloc(endsite * sizeof (double));
 236     for (j = 0; j < endsite; j++)
 237       treenode[i]->x[j]  = (ratelike)Malloc(rcategs*sizeof(sitelike));
 238   }
 239   if (!usertree) {
 240     for (i = spp; i < nonodes; i++) {
 241       p = treenode[i];
 242       for (j = 1; j <= 3; j++) {
 243         p->underflows = Malloc (endsite * sizeof (double));
 244         p->x = (phenotype)Malloc(endsite*sizeof(ratelike));
 245         for (k = 0; k < endsite; k++)
 246           p->x[k] = (ratelike)Malloc(rcategs*sizeof(sitelike));
 247         p = p->next;
 248       }
 249     }
 250   }
 251 }  /* allocx */
 252
 253
 254 void prot_allocx(long nonodes, long rcategs, pointarray treenode,
 255                         boolean usertree)
 256 {
 257   /* allocate x dynamically */
 258   /* used in proml          */
 259   long i, j, k;
 260   node *p;
 261
 262   for (i = 0; i < spp; i++){
 263     treenode[i]->protx = (pphenotype)Malloc(endsite*sizeof(pratelike));
 264     treenode[i]->underflows = Malloc(endsite*sizeof(double));
 265     for (j = 0; j < endsite; j++)
 266       treenode[i]->protx[j]  = (pratelike)Malloc(rcategs*sizeof(psitelike));
 267   }
 268   if (!usertree) {
 269     for (i = spp; i < nonodes; i++) {
 270       p = treenode[i];
 271       for (j = 1; j <= 3; j++) {
 272         p->protx = (pphenotype)Malloc(endsite*sizeof(pratelike));
 273         p->underflows = Malloc(endsite*sizeof(double));
 274         for (k = 0; k < endsite; k++)
 275           p->protx[k] = (pratelike)Malloc(rcategs*sizeof(psitelike));
 276         p = p->next;
 277       }
 278     }
 279   }
 280 }  /* prot_allocx */
 281
 282
 283 void allocx2(long nonodes, long endsite, long sitelength, pointarray treenode,
 284    boolean usertree)
 285 {
 286   /* allocate x2 dynamically */
 287   /* used in restml */
 288   long i, j, k, l;
 289   node *p;
 290
 291   for (i = 0; i < spp; i++) {
 292     treenode[i]->x2 = (phenotype2)Malloc(endsite*sizeof(sitelike2));
 293     for ( j = 0 ; j < endsite ; j++ )
 294       treenode[i]->x2[j] = Malloc((sitelength + 1) * sizeof(double));
 295   }
 296   if (!usertree) {
 297     for (i = spp; i < nonodes; i++) {
 298       p = treenode[i];
 299       for (j = 1; j <= 3; j++) {
 300         p->x2 = (phenotype2)Malloc(endsite*sizeof(sitelike2));
 301         for (k = 0; k < endsite; k++) {
 302           p->x2[k] = Malloc((sitelength + 1) * sizeof(double));
 303           for (l = 0; l < sitelength; l++)
 304              p->x2[k][l] = 1.0;
 305         }
 306         p = p->next;
 307       }
 308     }
 309   }
 310 }  /* allocx2 */
 311
 312
 313 void setuptree(pointarray treenode, long nonodes, boolean usertree)
 314 {
 315   /* initialize treenodes */
 316   long i;
 317   node *p;
 318
 319   for (i = 1; i <= nonodes; i++) {
 320     if (i <= spp || !usertree) {
 321       treenode[i-1]->back = NULL;
 322       treenode[i-1]->tip = (i <= spp);
 323       treenode[i-1]->index = i;
 324       treenode[i-1]->numdesc = 0;
 325       treenode[i-1]->iter = true;
 326       treenode[i-1]->initialized = true;
 327       treenode[i-1]->tyme =  0.0;
 328     }
 329   }
 330   if (!usertree) {
 331     for (i = spp + 1; i <= nonodes; i++) {
 332       p = treenode[i-1]->next;
 333       while (p != treenode[i-1]) {
 334         p->back = NULL;
 335         p->tip = false;
 336         p->index = i;
 337         p->numdesc = 0;
 338         p->iter = true;
 339         p->initialized = false;
 340         p->tyme = 0.0;
 341         p = p->next;
 342       }
 343     }
 344   }
 345 } /* setuptree */
 346
 347
 348 void setuptree2(tree a)
 349 {
 350   /* initialize a tree */
 351   /* used in dnaml, dnamlk, & restml */
 352
 353   a.likelihood = -999999.0;
 354   a.start = a.nodep[0]->back;
 355   a.root = NULL;
 356 } /* setuptree2 */
 357
 358
 359 void alloctip(node *p, long *zeros)
 360 { /* allocate a tip node */
 361   /* used by dnacomp, dnapars, & dnapenny */
 362
 363   p->numsteps = (steptr)Malloc(endsite*sizeof(long));
 364   p->oldnumsteps = (steptr)Malloc(endsite*sizeof(long));
 365   p->base = (baseptr)Malloc(endsite*sizeof(long));
 366   p->oldbase = (baseptr)Malloc(endsite*sizeof(long));
 367   memcpy(p->base, zeros, endsite*sizeof(long));
 368   memcpy(p->numsteps, zeros, endsite*sizeof(long));
 369   memcpy(p->oldbase, zeros, endsite*sizeof(long));
 370   memcpy(p->oldnumsteps, zeros, endsite*sizeof(long));
 371 }  /* alloctip */
 372
 373
 374 void freetrans(transptr *trans, long nonodes,long sitelength)
 375 {
 376   long i ,j;
 377   for ( i = 0 ; i < nonodes ; i++ ) {
 378     for ( j = 0 ; j < sitelength + 1; j++) {
 379       free ((*trans)[i][j]);
 380     }
 381     free ((*trans)[i]);
 382   }
 383   free(*trans);
 384 }
 385
 386
 387 void getbasefreqs(double freqa, double freqc, double freqg, double freqt,
 388             double *freqr, double *freqy, double *freqar, double *freqcy,
 389             double *freqgr, double *freqty, double *ttratio, double *xi,
 390             double *xv, double *fracchange, boolean freqsfrom,
 391             boolean printdata)
 392 {
 393   /* used by dnadist, dnaml, & dnamlk */
 394   double aa, bb;
 395
 396   if (printdata) {
 397     putc('\n', outfile);
 398     if (freqsfrom)
 399       fprintf(outfile, "Empirical ");
 400     fprintf(outfile, "Base Frequencies:\n\n");
 401     fprintf(outfile, "   A    %10.5f\n", freqa);
 402     fprintf(outfile, "   C    %10.5f\n", freqc);
 403     fprintf(outfile, "   G    %10.5f\n", freqg);
 404     fprintf(outfile, "  T(U)  %10.5f\n", freqt);
 405   }
 406   *freqr = freqa + freqg;
 407   *freqy = freqc + freqt;
 408   *freqar = freqa / *freqr;
 409   *freqcy = freqc / *freqy;
 410   *freqgr = freqg / *freqr;
 411   *freqty = freqt / *freqy;
 412   aa = *ttratio * (*freqr) * (*freqy) - freqa * freqg - freqc * freqt;
 413   bb = freqa * (*freqgr) + freqc * (*freqty);
 414   *xi = aa / (aa + bb);
 415   *xv = 1.0 - *xi;
 416   if (*xi < 0.0) {
 417     printf("\n WARNING: This transition/transversion ratio\n");
 418     printf(" is impossible with these base frequencies!\n");
 419     *xi = 0.0;
 420     *xv = 1.0;
 421     (*ttratio) = (freqa*freqg+freqc*freqt)/((*freqr)*(*freqy));
 422     printf(" Transition/transversion parameter reset\n");
 423     printf("  so transition/transversion ratio is %10.6f\n\n", (*ttratio));
 424   }
 425   if (freqa <= 0.0)
 426     freqa = 0.000001;
 427   if (freqc <= 0.0)
 428     freqc = 0.000001;
 429   if (freqg <= 0.0)
 430     freqg = 0.000001;
 431   if (freqt <= 0.0)
 432     freqt = 0.000001;
 433   *fracchange = (*xi) * (2 * freqa * (*freqgr) + 2 * freqc * (*freqty)) +
 434       (*xv) * (1.0 - freqa * freqa - freqc * freqc - freqg * freqg
 435       - freqt * freqt);
 436 }  /* getbasefreqs */
 437
 438
 439 void empiricalfreqs(double *freqa, double *freqc, double *freqg,
 440                         double *freqt, steptr weight, pointarray treenode)
 441 {
 442   /* Get empirical base frequencies from the data */
 443   /* used in dnaml & dnamlk */
 444   long i, j, k;
 445   double sum, suma, sumc, sumg, sumt, w;
 446
 447   *freqa = 0.25;
 448   *freqc = 0.25;
 449   *freqg = 0.25;
 450   *freqt = 0.25;
 451   for (k = 1; k <= 8; k++) {
 452     suma = 0.0;
 453     sumc = 0.0;
 454     sumg = 0.0;
 455     sumt = 0.0;
 456     for (i = 0; i < spp; i++) {
 457       for (j = 0; j < endsite; j++) {
 458         w = weight[j];
 459         sum = (*freqa) * treenode[i]->x[j][0][0];
 460         sum += (*freqc) * treenode[i]->x[j][0][(long)C - (long)A];
 461         sum += (*freqg) * treenode[i]->x[j][0][(long)G - (long)A];
 462         sum += (*freqt) * treenode[i]->x[j][0][(long)T - (long)A];
 463         suma += w * (*freqa) * treenode[i]->x[j][0][0] / sum;
 464         sumc += w * (*freqc) * treenode[i]->x[j][0][(long)C - (long)A] / sum;
 465         sumg += w * (*freqg) * treenode[i]->x[j][0][(long)G - (long)A] / sum;
 466         sumt += w * (*freqt) * treenode[i]->x[j][0][(long)T - (long)A] / sum;
 467       }
 468     }
 469     sum = suma + sumc + sumg + sumt;
 470     *freqa = suma / sum;
 471     *freqc = sumc / sum;
 472     *freqg = sumg / sum;
 473     *freqt = sumt / sum;
 474   }
 475   if (*freqa <= 0.0)
 476     *freqa = 0.000001;
 477   if (*freqc <= 0.0)
 478     *freqc = 0.000001;
 479   if (*freqg <= 0.0)
 480     *freqg = 0.000001;
 481   if (*freqt <= 0.0)
 482     *freqt = 0.000001;
 483 }  /* empiricalfreqs */
 484
 485
 486 void sitesort(long chars, steptr weight)
 487 {
 488   /* Shell sort keeping sites, weights in same order */
 489   /* used in dnainvar, dnapars, dnacomp & dnapenny */
 490   long gap, i, j, jj, jg, k, itemp;
 491   boolean flip, tied;
 492
 493   gap = chars / 2;
 494   while (gap > 0) {
 495     for (i = gap + 1; i <= chars; i++) {
 496       j = i - gap;
 497       flip = true;
 498       while (j > 0 && flip) {
 499         jj = alias[j - 1];
 500         jg = alias[j + gap - 1];
 501         tied = true;
 502         k = 1;
 503         while (k <= spp && tied) {
 504           flip = (y[k - 1][jj - 1] > y[k - 1][jg - 1]);
 505           tied = (tied && y[k - 1][jj - 1] == y[k - 1][jg - 1]);
 506           k++;
 507         }
 508         if (!flip)
 509           break;
 510         itemp = alias[j - 1];
 511         alias[j - 1] = alias[j + gap - 1];
 512         alias[j + gap - 1] = itemp;
 513         itemp = weight[j - 1];
 514         weight[j - 1] = weight[j + gap - 1];
 515         weight[j + gap - 1] = itemp;
 516         j -= gap;
 517       }
 518     }
 519     gap /= 2;
 520   }
 521 }  /* sitesort */
 522
 523
 524 void sitecombine(long chars)
 525 {
 526   /* combine sites that have identical patterns */
 527   /* used in dnapars, dnapenny, & dnacomp */
 528   long i, j, k;
 529   boolean tied;
 530
 531   i = 1;
 532   while (i < chars) {
 533     j = i + 1;
 534     tied = true;
 535     while (j <= chars && tied) {
 536       k = 1;
 537       while (k <= spp && tied) {
 538         tied = (tied &&
 539             y[k - 1][alias[i - 1] - 1] == y[k - 1][alias[j - 1] - 1]);
 540         k++;
 541       }
 542       if (tied) {
 543         weight[i - 1] += weight[j - 1];
 544         weight[j - 1] = 0;
 545         ally[alias[j - 1] - 1] = alias[i - 1];
 546       }
 547       j++;
 548     }
 549     i = j - 1;
 550   }
 551 }  /* sitecombine */
 552
 553
 554 void sitescrunch(long chars)
 555 {
 556   /* move so one representative of each pattern of
 557      sites comes first */
 558   /* used in dnapars & dnacomp */
 559   long i, j, itemp;
 560   boolean done, found;
 561
 562   done = false;
 563   i = 1;
 564   j = 2;
 565   while (!done) {
 566     if (ally[alias[i - 1] - 1] != alias[i - 1]) {
 567       if (j <= i)
 568         j = i + 1;
 569       if (j <= chars) {
 570         do {
 571           found = (ally[alias[j - 1] - 1] == alias[j - 1]);
 572           j++;
 573         } while (!(found || j > chars));
 574         if (found) {
 575           j--;
 576           itemp = alias[i - 1];
 577           alias[i - 1] = alias[j - 1];
 578           alias[j - 1] = itemp;
 579           itemp = weight[i - 1];
 580           weight[i - 1] = weight[j - 1];
 581           weight[j - 1] = itemp;
 582         } else
 583           done = true;
 584       } else
 585         done = true;
 586     }
 587     i++;
 588     done = (done || i >= chars);
 589   }
 590 }  /* sitescrunch */
 591
 592
 593 void sitesort2(long sites, steptr aliasweight)
 594 {
 595   /* Shell sort keeping sites, weights in same order */
 596   /* used in dnaml & dnamnlk */
 597   long gap, i, j, jj, jg, k, itemp;
 598   boolean flip, tied, samewt;
 599
 600   gap = sites / 2;
 601   while (gap > 0) {
 602     for (i = gap + 1; i <= sites; i++) {
 603       j = i - gap;
 604       flip = true;
 605       while (j > 0 && flip) {
 606         jj = alias[j - 1];
 607         jg = alias[j + gap - 1];
 608         samewt = ((weight[jj - 1] != 0) && (weight[jg - 1] != 0))
 609                  || ((weight[jj - 1] == 0) && (weight[jg - 1] == 0));
 610         tied = samewt && (category[jj - 1] == category[jg - 1]);
 611         flip = ((!samewt) && (weight[jj - 1] == 0))
 612                || (samewt && (category[jj - 1] > category[jg - 1]));
 613         k = 1;
 614         while (k <= spp && tied) {
 615           flip = (y[k - 1][jj - 1] > y[k - 1][jg - 1]);
 616           tied = (tied && y[k - 1][jj - 1] == y[k - 1][jg - 1]);
 617           k++;
 618         }
 619         if (!flip)
 620           break;
 621         itemp = alias[j - 1];
 622         alias[j - 1] = alias[j + gap - 1];
 623         alias[j + gap - 1] = itemp;
 624         itemp = aliasweight[j - 1];
 625         aliasweight[j - 1] = aliasweight[j + gap - 1];
 626         aliasweight[j + gap - 1] = itemp;
 627         j -= gap;
 628       }
 629     }
 630     gap /= 2;
 631   }
 632 }  /* sitesort2 */
 633
 634
 635 void sitecombine2(long sites, steptr aliasweight)
 636 {
 637   /* combine sites that have identical patterns */
 638   /* used in dnaml & dnamlk */
 639   long i, j, k;
 640   boolean tied, samewt;
 641
 642   i = 1;
 643   while (i < sites) {
 644     j = i + 1;
 645     tied = true;
 646     while (j <= sites && tied) {
 647       samewt = ((aliasweight[i - 1] != 0) && (aliasweight[j - 1] != 0))
 648                || ((aliasweight[i - 1] == 0) && (aliasweight[j - 1] == 0));
 649       tied = samewt
 650              && (category[alias[i - 1] - 1] == category[alias[j - 1] - 1]);
 651       k = 1;
 652       while (k <= spp && tied) {
 653         tied = (tied &&
 654             y[k - 1][alias[i - 1] - 1] == y[k - 1][alias[j - 1] - 1]);
 655         k++;
 656       }
 657       if (!tied)
 658         break;
 659       aliasweight[i - 1] += aliasweight[j - 1];
 660       aliasweight[j - 1] = 0;
 661       ally[alias[j - 1] - 1] = alias[i - 1];
 662       j++;
 663     }
 664     i = j;
 665   }
 666 }  /* sitecombine2 */
 667
 668
 669 void sitescrunch2(long sites, long i, long j, steptr aliasweight)
 670 {
 671   /* move so positively weighted sites come first */
 672   /* used by dnainvar, dnaml, dnamlk, & restml */
 673   long itemp;
 674   boolean done, found;
 675
 676   done = false;
 677   while (!done) {
 678     if (aliasweight[i - 1] > 0)
 679       i++;
 680     else {
 681       if (j <= i)
 682         j = i + 1;
 683       if (j <= sites) {
 684         do {
 685           found = (aliasweight[j - 1] > 0);
 686           j++;
 687         } while (!(found || j > sites));
 688         if (found) {
 689           j--;
 690           itemp = alias[i - 1];
 691           alias[i - 1] = alias[j - 1];
 692           alias[j - 1] = itemp;
 693           itemp = aliasweight[i - 1];
 694           aliasweight[i - 1] = aliasweight[j - 1];
 695           aliasweight[j - 1] = itemp;
 696         } else
 697           done = true;
 698       } else
 699         done = true;
 700     }
 701     done = (done || i >= sites);
 702   }
 703 }  /* sitescrunch2 */
 704
 705
 706 void makevalues(pointarray treenode, long *zeros, boolean usertree)
 707 {
 708   /* set up fractional likelihoods at tips */
 709   /* used by dnacomp, dnapars, & dnapenny */
 710   long i, j;
 711   char ns = 0;
 712   node *p;
 713
 714   setuptree(treenode, nonodes, usertree);
 715   for (i = 0; i < spp; i++)
 716     alloctip(treenode[i], zeros);
 717   if (!usertree) {
 718     for (i = spp; i < nonodes; i++) {
 719       p = treenode[i];
 720       do {
 721         allocnontip(p, zeros, endsite);
 722         p = p->next;
 723       } while (p != treenode[i]);
 724     }
 725   }
 726   for (j = 0; j < endsite; j++) {
 727     for (i = 0; i < spp; i++) {
 728       switch (y[i][alias[j] - 1]) {
 729
 730       case 'A':
 731         ns = 1 << A;
 732         break;
 733
 734       case 'C':
 735         ns = 1 << C;
 736         break;
 737
 738       case 'G':
 739         ns = 1 << G;
 740         break;
 741
 742       case 'U':
 743         ns = 1 << T;
 744         break;
 745
 746       case 'T':
 747         ns = 1 << T;
 748         break;
 749
 750       case 'M':
 751         ns = (1 << A) | (1 << C);
 752         break;
 753
 754       case 'R':
 755         ns = (1 << A) | (1 << G);
 756         break;
 757
 758       case 'W':
 759         ns = (1 << A) | (1 << T);
 760         break;
 761
 762       case 'S':
 763         ns = (1 << C) | (1 << G);
 764         break;
 765
 766       case 'Y':
 767         ns = (1 << C) | (1 << T);
 768         break;
 769
 770       case 'K':
 771         ns = (1 << G) | (1 << T);
 772         break;
 773
 774       case 'B':
 775         ns = (1 << C) | (1 << G) | (1 << T);
 776         break;
 777
 778       case 'D':
 779         ns = (1 << A) | (1 << G) | (1 << T);
 780         break;
 781
 782       case 'H':
 783         ns = (1 << A) | (1 << C) | (1 << T);
 784         break;
 785
 786       case 'V':
 787         ns = (1 << A) | (1 << C) | (1 << G);
 788         break;
 789
 790       case 'N':
 791         ns = (1 << A) | (1 << C) | (1 << G) | (1 << T);
 792         break;
 793
 794       case 'X':
 795         ns = (1 << A) | (1 << C) | (1 << G) | (1 << T);
 796         break;
 797
 798       case '?':
 799         ns = (1 << A) | (1 << C) | (1 << G) | (1 << T) | (1 << O);
 800         break;
 801
 802       case 'O':
 803         ns = 1 << O;
 804         break;
 805
 806       case '-':
 807         ns = 1 << O;
 808         break;
 809       }
 810       treenode[i]->base[j] = ns;
 811       treenode[i]->numsteps[j] = 0;
 812     }
 813   }
 814 }  /* makevalues */
 815
 816
 817 void makevalues2(long categs, pointarray treenode, long endsite,
 818                         long spp, sequence y, steptr alias)
 819 {
 820   /* set up fractional likelihoods at tips */
 821   /* used by dnaml & dnamlk */
 822   long i, j, k, l;
 823   bases b;
 824
 825   for (k = 0; k < endsite; k++) {
 826     j = alias[k];
 827     for (i = 0; i < spp; i++) {
 828       for (l = 0; l < categs; l++) {
 829         for (b = A; (long)b <= (long)T; b = (bases)((long)b + 1))
 830           treenode[i]->x[k][l][(long)b - (long)A] = 0.0;
 831         switch (y[i][j - 1]) {
 832
 833         case 'A':
 834           treenode[i]->x[k][l][0] = 1.0;
 835           break;
 836
 837         case 'C':
 838           treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
 839           break;
 840
 841         case 'G':
 842           treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
 843           break;
 844
 845         case 'T':
 846           treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
 847           break;
 848
 849         case 'U':
 850           treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
 851           break;
 852
 853         case 'M':
 854           treenode[i]->x[k][l][0] = 1.0;
 855           treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
 856           break;
 857
 858         case 'R':
 859           treenode[i]->x[k][l][0] = 1.0;
 860           treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
 861           break;
 862
 863         case 'W':
 864           treenode[i]->x[k][l][0] = 1.0;
 865           treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
 866           break;
 867
 868         case 'S':
 869           treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
 870           treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
 871           break;
 872
 873         case 'Y':
 874           treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
 875           treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
 876           break;
 877
 878         case 'K':
 879           treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
 880           treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
 881           break;
 882
 883         case 'B':
 884           treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
 885           treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
 886           treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
 887           break;
 888
 889         case 'D':
 890           treenode[i]->x[k][l][0] = 1.0;
 891           treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
 892           treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
 893           break;
 894
 895         case 'H':
 896           treenode[i]->x[k][l][0] = 1.0;
 897           treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
 898           treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
 899           break;
 900
 901         case 'V':
 902           treenode[i]->x[k][l][0] = 1.0;
 903           treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
 904           treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
 905           break;
 906
 907         case 'N':
 908           for (b = A; (long)b <= (long)T; b = (bases)((long)b + 1))
 909             treenode[i]->x[k][l][(long)b - (long)A] = 1.0;
 910           break;
 911
 912         case 'X':
 913           for (b = A; (long)b <= (long)T; b = (bases)((long)b + 1))
 914             treenode[i]->x[k][l][(long)b - (long)A] = 1.0;
 915           break;
 916
 917         case '?':
 918           for (b = A; (long)b <= (long)T; b = (bases)((long)b + 1))
 919             treenode[i]->x[k][l][(long)b - (long)A] = 1.0;
 920           break;
 921
 922         case 'O':
 923           for (b = A; (long)b <= (long)T; b = (bases)((long)b + 1))
 924             treenode[i]->x[k][l][(long)b - (long)A] = 1.0;
 925           break;
 926
 927         case '-':
 928           for (b = A; (long)b <= (long)T; b = (bases)((long)b + 1))
 929             treenode[i]->x[k][l][(long)b - (long)A] = 1.0;
 930           break;
 931         }
 932       }
 933     }
 934   }
 935 }  /* makevalues2 */
 936
 937
 938 void fillin(node *p, node *left, node *rt)
 939 {
 940   /* sets up for each node in the tree the base sequence
 941      at that point and counts the changes.  */
 942   long i, j, k, n, purset, pyrset;
 943   node *q;
 944
 945   purset = (1 << (long)A) + (1 << (long)G);
 946   pyrset = (1 << (long)C) + (1 << (long)T);
 947   if (!left) {
 948     memcpy(p->base, rt->base, endsite*sizeof(long));
 949     memcpy(p->numsteps, rt->numsteps, endsite*sizeof(long));
 950     q = rt;
 951   } else if (!rt) {
 952     memcpy(p->base, left->base, endsite*sizeof(long));
 953     memcpy(p->numsteps, left->numsteps, endsite*sizeof(long));
 954     q = left;
 955   } else {
 956     for (i = 0; i < endsite; i++) {
 957       p->base[i] = left->base[i] & rt->base[i];
 958       p->numsteps[i] = left->numsteps[i] + rt->numsteps[i];
 959       if (p->base[i] == 0) {
 960         p->base[i] = left->base[i] | rt->base[i];
 961         if (transvp) {
 962           if (!((p->base[i] == purset) || (p->base[i] == pyrset)))
 963             p->numsteps[i] += weight[i];
 964         }
 965         else p->numsteps[i] += weight[i];
 966       }
 967     }
 968     q = rt;
 969   }
 970   if (left && rt) n = 2;
 971   else n = 1;
 972   for (i = 0; i < endsite; i++)
 973     for (j = (long)A; j <= (long)O; j++)
 974       p->numnuc[i][j] = 0;
 975   for (k = 1; k <= n; k++) {
 976     if (k == 2) q = left;
 977     for (i = 0; i < endsite; i++) {
 978       for (j = (long)A; j <= (long)O; j++) {
 979         if (q->base[i] & (1 << j))
 980           p->numnuc[i][j]++;
 981       }
 982     }
 983   }
 984 }  /* fillin */
 985
 986
 987 long getlargest(long *numnuc)
 988 {
 989   /* find the largest in array numnuc */
 990   long i, largest;
 991
 992   largest = 0;
 993   for (i = (long)A; i <= (long)O; i++)
 994     if (numnuc[i] > largest)
 995       largest = numnuc[i];
 996   return largest;
 997 } /* getlargest */
 998
 999
1000 void multifillin(node *p, node *q, long dnumdesc)
1001 {
1002   /* sets up for each node in the tree the base sequence
1003      at that point and counts the changes according to the
1004      changes in q's base */
1005   long i, j, b, largest, descsteps, purset, pyrset;
1006
1007   memcpy(p->oldbase, p->base, endsite*sizeof(long));
1008   memcpy(p->oldnumsteps, p->numsteps, endsite*sizeof(long));
1009   purset = (1 << (long)A) + (1 << (long)G);
1010   pyrset = (1 << (long)C) + (1 << (long)T);
1011   for (i = 0; i < endsite; i++) {
1012     descsteps = 0;
1013     for (j = (long)A; j <= (long)O; j++) {
1014       b = 1 << j;
1015       if ((descsteps == 0) && (p->base[i] & b))
1016         descsteps = p->numsteps[i]
1017                       - (p->numdesc - dnumdesc - p->numnuc[i][j]) * weight[i];
1018     }
1019     if (dnumdesc == -1)
1020       descsteps -= q->oldnumsteps[i];
1021     else if (dnumdesc == 0)
1022       descsteps += (q->numsteps[i] - q->oldnumsteps[i]);
1023     else
1024       descsteps += q->numsteps[i];
1025     if (q->oldbase[i] != q->base[i]) {
1026       for (j = (long)A; j <= (long)O; j++) {
1027         b = 1 << j;
1028         if (transvp) {
1029           if (b & purset) b = purset;
1030           if (b & pyrset) b = pyrset;
1031         }
1032         if ((q->oldbase[i] & b) && !(q->base[i] & b))
1033           p->numnuc[i][j]--;
1034         else if (!(q->oldbase[i] & b) && (q->base[i] & b))
1035           p->numnuc[i][j]++;
1036       }
1037     }
1038     largest = getlargest(p->numnuc[i]);
1039     if (q->oldbase[i] != q->base[i]) {
1040       p->base[i] = 0;
1041       for (j = (long)A; j <= (long)O; j++) {
1042         if (p->numnuc[i][j] == largest)
1043             p->base[i] |= (1 << j);
1044       }
1045     }
1046     p->numsteps[i] = (p->numdesc - largest) * weight[i] + descsteps;
1047   }
1048 } /* multifillin */
1049
1050
1051 void sumnsteps(node *p, node *left, node *rt, long a, long b)
1052 {
1053   /* sets up for each node in the tree the base sequence
1054      at that point and counts the changes. */
1055   long i;
1056   long ns, rs, ls, purset, pyrset;
1057
1058   if (!left) {
1059     memcpy(p->numsteps, rt->numsteps, endsite*sizeof(long));
1060     memcpy(p->base, rt->base, endsite*sizeof(long));
1061   } else if (!rt) {
1062     memcpy(p->numsteps, left->numsteps, endsite*sizeof(long));
1063     memcpy(p->base, left->base, endsite*sizeof(long));
1064   } else  {
1065     purset = (1 << (long)A) + (1 << (long)G);
1066     pyrset = (1 << (long)C) + (1 << (long)T);
1067     for (i = a; i < b; i++) {
1068       ls = left->base[i];
1069       rs = rt->base[i];
1070       ns = ls & rs;
1071       p->numsteps[i] = left->numsteps[i] + rt->numsteps[i];
1072       if (ns == 0) {
1073         ns = ls | rs;
1074         if (transvp) {
1075           if (!((ns == purset) || (ns == pyrset)))
1076             p->numsteps[i] += weight[i];
1077         }
1078         else p->numsteps[i] += weight[i];
1079       }
1080       p->base[i] = ns;
1081       }
1082     }
1083 }  /* sumnsteps */
1084
1085
1086 void sumnsteps2(node *p,node *left,node *rt,long a,long b,long *threshwt)
1087 {
1088   /* counts the changes at each node.  */
1089   long i, steps;
1090   long ns, rs, ls, purset, pyrset;
1091   long term;
1092
1093   if (a == 0) p->sumsteps = 0.0;
1094   if (!left)
1095     memcpy(p->numsteps, rt->numsteps, endsite*sizeof(long));
1096   else if (!rt)
1097     memcpy(p->numsteps, left->numsteps, endsite*sizeof(long));
1098   else {
1099     purset = (1 << (long)A) + (1 << (long)G);
1100     pyrset = (1 << (long)C) + (1 << (long)T);
1101     for (i = a; i < b; i++) {
1102       ls = left->base[i];
1103       rs = rt->base[i];
1104       ns = ls & rs;
1105       p->numsteps[i] = left->numsteps[i] + rt->numsteps[i];
1106       if (ns == 0) {
1107         ns = ls | rs;
1108         if (transvp) {
1109           if (!((ns == purset) || (ns == pyrset)))
1110             p->numsteps[i] += weight[i];
1111         }
1112         else p->numsteps[i] += weight[i];
1113       }
1114     }
1115   }
1116   for (i = a; i < b; i++) {
1117     steps = p->numsteps[i];
1118     if ((long)steps <= threshwt[i])
1119       term = steps;
1120     else
1121       term = threshwt[i];
1122     p->sumsteps += (double)term;
1123   }
1124 }  /* sumnsteps2 */
1125
1126
1127 void multisumnsteps(node *p, node *q, long a, long b, long *threshwt)
1128 {
1129   /* computes the number of steps between p and q */
1130   long i, j, steps, largest, descsteps, purset, pyrset, b1;
1131   long term;
1132
1133   if (a == 0) p->sumsteps = 0.0;
1134   purset = (1 << (long)A) + (1 << (long)G);
1135   pyrset = (1 << (long)C) + (1 << (long)T);
1136   for (i = a; i < b; i++) {
1137     descsteps = 0;
1138     for (j = (long)A; j <= (long)O; j++) {
1139       if ((descsteps == 0) && (p->base[i] & (1 << j)))
1140         descsteps = p->numsteps[i] -
1141                         (p->numdesc - 1 - p->numnuc[i][j]) * weight[i];
1142     }
1143     descsteps += q->numsteps[i];
1144     largest = 0;
1145     for (j = (long)A; j <= (long)O; j++) {
1146       b1 = (1 << j);
1147       if (transvp) {
1148         if (b1 & purset) b1 = purset;
1149         if (b1 & pyrset) b1 = pyrset;
1150       }
1151       if (q->base[i] & b1)
1152         p->numnuc[i][j]++;
1153       if (p->numnuc[i][j] > largest)
1154         largest = p->numnuc[i][j];
1155     }
1156     steps = (p->numdesc - largest) * weight[i] + descsteps;
1157     if ((long)steps <= threshwt[i])
1158       term = steps;
1159     else
1160       term = threshwt[i];
1161     p->sumsteps += (double)term;
1162   }
1163 } /* multisumnsteps */
1164
1165
1166 void multisumnsteps2(node *p)
1167 {
1168   /* counts the changes at each multi-way node. Sums up
1169      steps of all descendants */
1170   long i, j, largest, purset, pyrset, b1;
1171   node *q;
1172   baseptr b;
1173
1174   purset = (1 << (long)A) + (1 << (long)G);
1175   pyrset = (1 << (long)C) + (1 << (long)T);
1176   for (i = 0; i < endsite; i++) {
1177     p->numsteps[i] = 0;
1178     q = p->next;
1179     while (q != p) {
1180       if (q->back) {
1181         p->numsteps[i] += q->back->numsteps[i];
1182         b = q->back->base;
1183         for (j = (long)A; j <= (long)O; j++) {
1184           b1 = (1 << j);
1185           if (transvp) {
1186             if (b1 & purset) b1 = purset;
1187             if (b1 & pyrset) b1 = pyrset;
1188           }
1189           if (b[i] & b1) p->numnuc[i][j]++;
1190         }
1191       }
1192       q = q->next;
1193     }
1194     largest = getlargest(p->numnuc[i]);
1195     p->base[i] = 0;
1196     for (j = (long)A; j <= (long)O; j++) {
1197       if (p->numnuc[i][j] == largest)
1198         p->base[i] |= (1 << j);
1199     }
1200     p->numsteps[i] += ((p->numdesc - largest) * weight[i]);
1201   }
1202 }  /* multisumnsteps2 */
1203
1204 boolean alltips(node *forknode, node *p)
1205 {
1206   /* returns true if all descendants of forknode except p are tips;
1207      false otherwise.  */
1208   node *q, *r;
1209   boolean tips;
1210
1211   tips = true;
1212   r = forknode;
1213   q = forknode->next;
1214   do {
1215     if (q->back && q->back != p && !q->back->tip)
1216       tips = false;
1217     q = q->next;
1218   } while (tips && q != r);
1219   return tips;
1220 } /* alltips */
1221
1222
1223 void gdispose(node *p, node **grbg, pointarray treenode)
1224 {
1225   /* go through tree throwing away nodes */
1226   node *q, *r;
1227
1228   p->back = NULL;
1229   if (p->tip)
1230     return;
1231   treenode[p->index - 1] = NULL;
1232   q = p->next;
1233   while (q != p) {
1234     gdispose(q->back, grbg, treenode);
1235     q->back = NULL;
1236     r = q;
1237     q = q->next;
1238     chucktreenode(grbg, r);
1239   }
1240   chucktreenode(grbg, q);
1241 }  /* gdispose */
1242
1243
1244 void preorder(node *p, node *r, node *root, node *removing, node *adding,
1245                         node *changing, long dnumdesc)
1246 {
1247   /* recompute number of steps in preorder taking both ancestoral and
1248      descendent steps into account. removing points to a node being
1249      removed, if any */
1250   node *q, *p1, *p2;
1251
1252   if (p && !p->tip && p != adding) {
1253     q = p;
1254     do {
1255       if (p->back != r) {
1256         if (p->numdesc > 2) {
1257           if (changing)
1258             multifillin (p, r, dnumdesc);
1259           else
1260             multifillin (p, r, 0);
1261         } else {
1262           p1 = p->next;
1263           if (!removing)
1264             while (!p1->back)
1265               p1 = p1->next;
1266           else
1267             while (!p1->back || p1->back == removing)
1268               p1 = p1->next;
1269           p2 = p1->next;
1270           if (!removing)
1271             while (!p2->back)
1272               p2 = p2->next;
1273           else
1274             while (!p2->back || p2->back == removing)
1275               p2 = p2->next;
1276           p1 = p1->back;
1277           p2 = p2->back;
1278           if (p->back == p1) p1 = NULL;
1279           else if (p->back == p2) p2 = NULL;
1280           memcpy(p->oldbase, p->base, endsite*sizeof(long));
1281           memcpy(p->oldnumsteps, p->numsteps, endsite*sizeof(long));
1282           fillin(p, p1, p2);
1283         }
1284       }
1285       p = p->next;
1286     } while (p != q);
1287     q = p;
1288     do {
1289       preorder(p->next->back, p->next, root, removing, adding, NULL, 0);
1290       p = p->next;
1291     } while (p->next != q);
1292   }
1293 } /* preorder */
1294
1295
1296 void updatenumdesc(node *p, node *root, long n)
1297 {
1298   /* set p's numdesc to n.  If p is the root, numdesc of p's
1299   descendants are set to n-1. */
1300   node *q;
1301
1302   q = p;
1303   if (p == root && n > 0) {
1304     p->numdesc = n;
1305     n--;
1306     q = q->next;
1307   }
1308   do {
1309     q->numdesc = n;
1310     q = q->next;
1311   } while (q != p);
1312 }  /* updatenumdesc */
1313
1314
1315 void add(node *below,node *newtip,node *newfork,node **root,
1316         boolean recompute,pointarray treenode,node **grbg,long *zeros)
1317 {
1318   /* inserts the nodes newfork and its left descendant, newtip,
1319      to the tree.  below becomes newfork's right descendant.
1320      if newfork is NULL, newtip is added as below's sibling */
1321   /* used in dnacomp & dnapars */
1322   node *p;
1323
1324   if (below != treenode[below->index - 1])
1325     below = treenode[below->index - 1];
1326   if (newfork) {
1327     if (below->back != NULL)
1328       below->back->back = newfork;
1329     newfork->back = below->back;
1330     below->back = newfork->next->next;
1331     newfork->next->next->back = below;
1332     newfork->next->back = newtip;
1333     newtip->back = newfork->next;
1334     if (*root == below)
1335       *root = newfork;
1336     updatenumdesc(newfork, *root, 2);
1337   } else {
1338     gnutreenode(grbg, &p, below->index, endsite, zeros);
1339     p->back = newtip;
1340     newtip->back = p;
1341     p->next = below->next;
1342     below->next = p;
1343     updatenumdesc(below, *root, below->numdesc + 1);
1344   }
1345   if (!newtip->tip)
1346     updatenumdesc(newtip, *root, newtip->numdesc);
1347   (*root)->back = NULL;
1348   if (!recompute)
1349     return;
1350   if (!newfork) {
1351     memcpy(newtip->back->base, below->base, endsite*sizeof(long));
1352     memcpy(newtip->back->numsteps, below->numsteps, endsite*sizeof(long));
1353     memcpy(newtip->back->numnuc, below->numnuc, endsite*sizeof(nucarray));
1354     if (below != *root) {
1355       memcpy(below->back->oldbase, zeros, endsite*sizeof(long));
1356       memcpy(below->back->oldnumsteps, zeros, endsite*sizeof(long));
1357       multifillin(newtip->back, below->back, 1);
1358     }
1359     if (!newtip->tip) {
1360       memcpy(newtip->back->oldbase, zeros, endsite*sizeof(long));
1361       memcpy(newtip->back->oldnumsteps, zeros, endsite*sizeof(long));
1362       preorder(newtip, newtip->back, *root, NULL, NULL, below, 1);
1363     }
1364     memcpy(newtip->oldbase, zeros, endsite*sizeof(long));
1365     memcpy(newtip->oldnumsteps, zeros, endsite*sizeof(long));
1366     preorder(below, newtip, *root, NULL, newtip, below, 1);
1367     if (below != *root)
1368       preorder(below->back, below, *root, NULL, NULL, NULL, 0);
1369   } else {
1370     fillin(newtip->back, newtip->back->next->back,
1371              newtip->back->next->next->back);
1372     if (!newtip->tip) {
1373       memcpy(newtip->back->oldbase, zeros, endsite*sizeof(long));
1374       memcpy(newtip->back->oldnumsteps, zeros, endsite*sizeof(long));
1375       preorder(newtip, newtip->back, *root, NULL, NULL, newfork, 1);
1376     }
1377     if (newfork != *root) {
1378       memcpy(below->back->base, newfork->back->base, endsite*sizeof(long));
1379       memcpy(below->back->numsteps, newfork->back->numsteps, endsite*sizeof(long));
1380       preorder(newfork, newtip, *root, NULL, newtip, NULL, 0);
1381     } else {
1382       fillin(below->back, newtip, NULL);
1383       fillin(newfork, newtip, below);
1384       memcpy(below->back->oldbase, zeros, endsite*sizeof(long));
1385       memcpy(below->back->oldnumsteps, zeros, endsite*sizeof(long));
1386       preorder(below, below->back, *root, NULL, NULL, newfork, 1);
1387     }
1388     if (newfork != *root) {
1389       memcpy(newfork->oldbase, below->base, endsite*sizeof(long));
1390       memcpy(newfork->oldnumsteps, below->numsteps, endsite*sizeof(long));
1391       preorder(newfork->back, newfork, *root, NULL, NULL, NULL, 0);
1392     }
1393   }
1394 }  /* add */
1395
1396
1397 void findbelow(node **below, node *item, node *fork)
1398 {
1399   /* decide which of fork's binary children is below */
1400
1401   if (fork->next->back == item)
1402     *below = fork->next->next->back;
1403   else
1404     *below = fork->next->back;
1405 } /* findbelow */
1406
1407
1408 void re_move(node *item, node **fork, node **root, boolean recompute,
1409                         pointarray treenode, node **grbg, long *zeros)
1410 {
1411   /* removes nodes item and its ancestor, fork, from the tree.
1412      the new descendant of fork's ancestor is made to be
1413      fork's second descendant (other than item).  Also
1414      returns pointers to the deleted nodes, item and fork.
1415      If item belongs to a node with more than 2 descendants,
1416      fork will not be deleted */
1417   /* used in dnacomp & dnapars */
1418   node *p, *q, *other = NULL, *otherback = NULL;
1419
1420   if (item->back == NULL) {
1421     *fork = NULL;
1422     return;
1423   }
1424   *fork = treenode[item->back->index - 1];
1425   if ((*fork)->numdesc == 2) {
1426     updatenumdesc(*fork, *root, 0);
1427     findbelow(&other, item, *fork);
1428     otherback = other->back;
1429     if (*root == *fork) {
1430       *root = other;
1431       if (!other->tip)
1432         updatenumdesc(other, *root, other->numdesc);
1433     }
1434     p = item->back->next->back;
1435     q = item->back->next->next->back;
1436     if (p != NULL)
1437       p->back = q;
1438     if (q != NULL)
1439       q->back = p;
1440     (*fork)->back = NULL;
1441     p = (*fork)->next;
1442     while (p != *fork) {
1443       p->back = NULL;
1444       p = p->next;
1445     }
1446   } else {
1447     updatenumdesc(*fork, *root, (*fork)->numdesc - 1);
1448     p = *fork;
1449     while (p->next != item->back)
1450       p = p->next;
1451     p->next = item->back->next;
1452   }
1453   if (!item->tip) {
1454     updatenumdesc(item, item, item->numdesc);
1455     if (recompute) {
1456       memcpy(item->back->oldbase, item->back->base, endsite*sizeof(long));
1457       memcpy(item->back->oldnumsteps, item->back->numsteps, endsite*sizeof(long));
1458       memcpy(item->back->base, zeros, endsite*sizeof(long));
1459       memcpy(item->back->numsteps, zeros, endsite*sizeof(long));
1460       preorder(item, item->back, *root, item->back, NULL, item, -1);
1461     }
1462   }
1463   if ((*fork)->numdesc >= 2)
1464     chucktreenode(grbg, item->back);
1465   item->back = NULL;
1466   if (!recompute)
1467     return;
1468   if ((*fork)->numdesc == 0) {
1469     memcpy(otherback->oldbase, otherback->base, endsite*sizeof(long));
1470     memcpy(otherback->oldnumsteps, otherback->numsteps, endsite*sizeof(long));
1471     if (other == *root) {
1472       memcpy(otherback->base, zeros, endsite*sizeof(long));
1473       memcpy(otherback->numsteps, zeros, endsite*sizeof(long));
1474     } else {
1475       memcpy(otherback->base, other->back->base, endsite*sizeof(long));
1476       memcpy(otherback->numsteps, other->back->numsteps, endsite*sizeof(long));
1477     }
1478     p = other->back;
1479     other->back = otherback;
1480     if (other == *root)
1481       preorder(other, otherback, *root, otherback, NULL, other, -1);
1482     else
1483       preorder(other, otherback, *root, NULL, NULL, NULL, 0);
1484     other->back = p;
1485     if (other != *root) {
1486       memcpy(other->oldbase,(*fork)->base, endsite*sizeof(long));
1487       memcpy(other->oldnumsteps,(*fork)->numsteps, endsite*sizeof(long));
1488       preorder(other->back, other, *root, NULL, NULL, NULL, 0);
1489     }
1490   } else {
1491     memcpy(item->oldbase, item->base, endsite*sizeof(long));
1492     memcpy(item->oldnumsteps, item->numsteps, endsite*sizeof(long));
1493     memcpy(item->base, zeros, endsite*sizeof(long));
1494     memcpy(item->numsteps, zeros, endsite*sizeof(long));
1495     preorder(*fork, item, *root, NULL, NULL, *fork, -1);
1496     if (*fork != *root)
1497       preorder((*fork)->back, *fork, *root, NULL, NULL, NULL, 0);
1498     memcpy(item->base, item->oldbase, endsite*sizeof(long));
1499     memcpy(item->numsteps, item->oldnumsteps, endsite*sizeof(long));
1500   }
1501 }  /* remove */
1502
1503
1504 void postorder(node *p)
1505 {
1506   /* traverses an n-ary tree, suming up steps at a node's descendants */
1507   /* used in dnacomp, dnapars, & dnapenny */
1508   node *q;
1509
1510   if (p->tip)
1511     return;
1512   q = p->next;
1513   while (q != p) {
1514     postorder(q->back);
1515     q = q->next;
1516   }
1517   zeronumnuc(p, endsite);
1518   if (p->numdesc > 2)
1519     multisumnsteps2(p);
1520   else
1521     fillin(p, p->next->back, p->next->next->back);
1522 }  /* postorder */
1523
1524
1525 void getnufork(node **nufork,node **grbg,pointarray treenode,long *zeros)
1526 {
1527   /* find a fork not used currently */
1528   long i;
1529
1530   i = spp;
1531   while (treenode[i] && treenode[i]->numdesc > 0) i++;
1532   if (!treenode[i])
1533     gnutreenode(grbg, &treenode[i], i, endsite, zeros);
1534   *nufork = treenode[i];
1535 } /* getnufork */
1536
1537
1538 void reroot(node *outgroup, node *root)
1539 {
1540   /* reorients tree, putting outgroup in desired position. used if
1541      the root is binary. */
1542   /* used in dnacomp & dnapars */
1543   node *p, *q;
1544
1545   if (outgroup->back->index == root->index)
1546     return;
1547   p = root->next;
1548   q = root->next->next;
1549   p->back->back = q->back;
1550   q->back->back = p->back;
1551   p->back = outgroup;
1552   q->back = outgroup->back;
1553   outgroup->back->back = q;
1554   outgroup->back = p;
1555 }  /* reroot */
1556
1557
1558 void reroot2(node *outgroup, node *root)
1559 {
1560   /* reorients tree, putting outgroup in desired position. */
1561   /* used in dnacomp & dnapars */
1562   node *p;
1563
1564   p = outgroup->back->next;
1565   while (p->next != outgroup->back)
1566     p = p->next;
1567   root->next = outgroup->back;
1568   p->next = root;
1569 }  /* reroot2 */
1570
1571
1572 void reroot3(node *outgroup, node *root, node *root2, node *lastdesc,
1573                         node **grbg)
1574 {
1575   /* reorients tree, putting back outgroup in original position. */
1576   /* used in dnacomp & dnapars */
1577   node *p;
1578
1579   p = root->next;
1580   while (p->next != root)
1581     p = p->next;
1582   chucktreenode(grbg, root);
1583   p->next = outgroup->back;
1584   root2->next = lastdesc->next;
1585   lastdesc->next = root2;
1586 }  /* reroot3 */
1587
1588
1589 void savetraverse(node *p)
1590 {
1591   /* sets BOOLEANs that indicate which way is down */
1592   node *q;
1593
1594   p->bottom = true;
1595   if (p->tip)
1596     return;
1597   q = p->next;
1598   while (q != p) {
1599     q->bottom = false;
1600     savetraverse(q->back);
1601     q = q->next;
1602   }
1603 }  /* savetraverse */
1604
1605
1606 void newindex(long i, node *p)
1607 {
1608   /* assigns index i to node p */
1609
1610   while (p->index != i) {
1611     p->index = i;
1612     p = p->next;
1613   }
1614 } /* newindex */
1615
1616
1617 void flipindexes(long nextnode, pointarray treenode)
1618 {
1619   /* flips index of nodes between nextnode and last node.  */
1620   long last;
1621   node *temp;
1622
1623   last = nonodes;
1624   while (treenode[last - 1]->numdesc == 0)
1625     last--;
1626   if (last > nextnode) {
1627     temp = treenode[nextnode - 1];
1628     treenode[nextnode - 1] = treenode[last - 1];
1629     treenode[last - 1] = temp;
1630     newindex(nextnode, treenode[nextnode - 1]);
1631     newindex(last, treenode[last - 1]);
1632   }
1633 } /* flipindexes */
1634
1635
1636 boolean parentinmulti(node *anode)
1637 {
1638   /* sees if anode's parent has more than 2 children */
1639   node *p;
1640
1641   while (!anode->bottom) anode = anode->next;
1642   p = anode->back;
1643   while (!p->bottom)
1644     p = p->next;
1645   return (p->numdesc > 2);
1646 } /* parentinmulti */
1647
1648
1649 long sibsvisited(node *anode, long *place)
1650 {
1651   /* computes the number of nodes which are visited earlier than anode among
1652   its siblings */
1653   node *p;
1654   long nvisited;
1655
1656   while (!anode->bottom) anode = anode->next;
1657   p = anode->back->next;
1658   nvisited = 0;
1659   do {
1660     if (!p->bottom && place[p->back->index - 1] != 0)
1661       nvisited++;
1662     p = p->next;
1663   } while (p != anode->back);
1664   return nvisited;
1665 }  /* sibsvisited */
1666
1667
1668 long smallest(node *anode, long *place)
1669 {
1670   /* finds the smallest index of sibling of anode */
1671   node *p;
1672   long min;
1673
1674   while (!anode->bottom) anode = anode->next;
1675   p = anode->back->next;
1676   if (p->bottom) p = p->next;
1677   min = nonodes;
1678   do {
1679     if (p->back && place[p->back->index - 1] != 0) {
1680       if (p->back->index <= spp) {
1681         if (p->back->index < min)
1682           min = p->back->index;
1683       } else {
1684         if (place[p->back->index - 1] < min)
1685           min = place[p->back->index - 1];
1686       }
1687     }
1688     p = p->next;
1689     if (p->bottom) p = p->next;
1690   } while (p != anode->back);
1691   return min;
1692 }  /* smallest */
1693
1694
1695 void bintomulti(node **root, node **binroot, node **grbg, long *zeros)
1696 {  /* attaches root's left child to its right child and makes
1697       the right child new root */
1698   node *left, *right, *newnode, *temp;
1699
1700   right = (*root)->next->next->back;
1701   left = (*root)->next->back;
1702   if (right->tip) {
1703     (*root)->next = right->back;
1704     (*root)->next->next = left->back;
1705     temp = left;
1706     left = right;
1707     right = temp;
1708     right->back->next = *root;
1709   }
1710   gnutreenode(grbg, &newnode, right->index, endsite, zeros);
1711   newnode->next = right->next;
1712   newnode->back = left;
1713   left->back = newnode;
1714   right->next = newnode;
1715   (*root)->next->back = (*root)->next->next->back = NULL;
1716   *binroot = *root;
1717   (*binroot)->numdesc = 0;
1718   *root = right;
1719   (*root)->numdesc++;
1720   (*root)->back = NULL;
1721 } /* bintomulti */
1722
1723
1724 void backtobinary(node **root, node *binroot, node **grbg)
1725 { /* restores binary root */
1726   node *p;
1727
1728   binroot->next->back = (*root)->next->back;
1729   (*root)->next->back->back = binroot->next;
1730   p = (*root)->next;
1731   (*root)->next = p->next;
1732   binroot->next->next->back = *root;
1733   (*root)->back = binroot->next->next;
1734   chucktreenode(grbg, p);
1735   (*root)->numdesc--;
1736   *root = binroot;
1737   (*root)->numdesc = 2;
1738 } /* backtobinary */
1739
1740
1741 boolean outgrin(node *root, node *outgrnode)
1742 { /* checks if outgroup node is a child of root */
1743   node *p;
1744
1745   p = root->next;
1746   while (p != root) {
1747     if (p->back == outgrnode)
1748       return true;
1749     p = p->next;
1750   }
1751   return false;
1752 } /* outgrin */
1753
1754
1755 void flipnodes(node *nodea, node *nodeb)
1756 { /* flip nodes */
1757   node *backa, *backb;
1758
1759   backa = nodea->back;
1760   backb = nodeb->back;
1761   backa->back = nodeb;
1762   backb->back = nodea;
1763   nodea->back = backb;
1764   nodeb->back = backa;
1765 } /* flipnodes */
1766
1767
1768 void moveleft(node *root, node *outgrnode, node **flipback)
1769 { /* makes outgroup node to leftmost child of root */
1770   node *p;
1771   boolean done;
1772
1773   p = root->next;
1774   done = false;
1775   while (p != root && !done) {
1776     if (p->back == outgrnode) {
1777       *flipback = p;
1778       flipnodes(root->next->back, p->back);
1779       done = true;
1780     }
1781     p = p->next;
1782   }
1783 } /* moveleft */
1784
1785
1786 void savetree(node *root,  long *place, pointarray treenode,
1787                         node **grbg, long *zeros)
1788 { /* record in place where each species has to be
1789      added to reconstruct this tree */
1790   /* used by dnacomp & dnapars */
1791   long i, j, nextnode, nvisited;
1792   node *p, *q, *r = NULL, *root2, *lastdesc,
1793                 *outgrnode, *binroot, *flipback;
1794   boolean done, newfork;
1795
1796   binroot = NULL;
1797   lastdesc = NULL;
1798   root2 = NULL;
1799   flipback = NULL;
1800   outgrnode = treenode[outgrno - 1];
1801   if (root->numdesc == 2)
1802     bintomulti(&root, &binroot, grbg, zeros);
1803   if (outgrin(root, outgrnode)) {
1804     if (outgrnode != root->next->back)
1805       moveleft(root, outgrnode, &flipback);
1806   } else {
1807     root2 = root;
1808     lastdesc = root->next;
1809     while (lastdesc->next != root)
1810       lastdesc = lastdesc->next;
1811     lastdesc->next = root->next;
1812     gnutreenode(grbg, &root, outgrnode->back->index, endsite, zeros);
1813     root->numdesc = root2->numdesc;
1814     reroot2(outgrnode, root);
1815   }
1816   savetraverse(root);
1817   nextnode = spp + 1;
1818   for (i = nextnode; i <= nonodes; i++)
1819     if (treenode[i - 1]->numdesc == 0)
1820       flipindexes(i, treenode);
1821   for (i = 0; i < nonodes; i++)
1822     place[i] = 0;
1823   place[root->index - 1] = 1;
1824   for (i = 1; i <= spp; i++) {
1825     p = treenode[i - 1];
1826     while (place[p->index - 1] == 0) {
1827       place[p->index - 1] = i;
1828       while (!p->bottom)
1829         p = p->next;
1830       r = p;
1831       p = p->back;
1832     }
1833     if (i > 1) {
1834       q = treenode[i - 1];
1835       newfork = true;
1836       nvisited = sibsvisited(q, place);
1837       if (nvisited == 0) {
1838         if (parentinmulti(r)) {
1839           nvisited = sibsvisited(r, place);
1840           if (nvisited == 0)
1841             place[i - 1] = place[p->index - 1];
1842           else if (nvisited == 1)
1843             place[i - 1] = smallest(r, place);
1844           else {
1845             place[i - 1] = -smallest(r, place);
1846             newfork = false;
1847           }
1848         } else
1849           place[i - 1] = place[p->index - 1];
1850       } else if (nvisited == 1) {
1851         place[i - 1] = place[p->index - 1];
1852       } else {
1853         place[i - 1] = -smallest(q, place);
1854         newfork = false;
1855       }
1856       if (newfork) {
1857         j = place[p->index - 1];
1858         done = false;
1859         while (!done) {
1860           place[p->index - 1] = nextnode;
1861           while (!p->bottom)
1862             p = p->next;
1863           p = p->back;
1864           done = (p == NULL);
1865           if (!done)
1866             done = (place[p->index - 1] != j);
1867           if (done) {
1868             nextnode++;
1869           }
1870         }
1871       }
1872     }
1873   }
1874   if (flipback)
1875     flipnodes(outgrnode, flipback->back);
1876   else {
1877     if (root2) {
1878       reroot3(outgrnode, root, root2, lastdesc, grbg);
1879       root = root2;
1880     }
1881   }
1882   if (binroot)
1883     backtobinary(&root, binroot, grbg);
1884 }  /* savetree */
1885
1886
1887 void addnsave(node *p, node *item, node *nufork, node **root, node **grbg,
1888                 boolean multf, pointarray treenode, long *place, long *zeros)
1889 {  /* adds item to tree and save it.  Then removes item. */
1890   node *dummy;
1891
1892   if (!multf)
1893     add(p, item, nufork, root, false, treenode, grbg, zeros);
1894   else
1895     add(p, item, NULL, root, false, treenode, grbg, zeros);
1896   savetree(*root, place, treenode, grbg, zeros);
1897   if (!multf)
1898     re_move(item, &nufork, root, false, treenode, grbg, zeros);
1899   else
1900     re_move(item, &dummy, root, false, treenode, grbg, zeros);
1901 } /* addnsave */
1902
1903
1904 void addbestever(long *pos, long *nextree, long maxtrees, boolean collapse,
1905                         long *place, bestelm *bestrees)
1906 { /* adds first best tree */
1907
1908   *pos = 1;
1909   *nextree = 1;
1910   initbestrees(bestrees, maxtrees, true);
1911   initbestrees(bestrees, maxtrees, false);
1912   addtree(*pos, nextree, collapse, place, bestrees);
1913 } /* addbestever */
1914
1915
1916 void addtiedtree(long pos, long *nextree, long maxtrees, boolean collapse,
1917                         long *place, bestelm *bestrees)
1918 { /* add tied tree */
1919
1920   if (*nextree <= maxtrees)
1921     addtree(pos, nextree, collapse, place, bestrees);
1922 } /* addtiedtree */
1923
1924
1925 void clearcollapse(pointarray treenode)
1926 {
1927   /* clears collapse status at a node */
1928   long i;
1929   node *p;
1930
1931   for (i = 0; i < nonodes; i++) {
1932     treenode[i]->collapse = undefined;
1933     if (!treenode[i]->tip) {
1934       p = treenode[i]->next;
1935       while (p != treenode[i]) {
1936         p->collapse = undefined;
1937         p = p->next;
1938       }
1939     }
1940   }
1941 } /* clearcollapse */
1942
1943
1944 void clearbottom(pointarray treenode)
1945 {
1946   /* clears boolean bottom at a node */
1947   long i;
1948   node *p;
1949
1950   for (i = 0; i < nonodes; i++) {
1951     treenode[i]->bottom = false;
1952     if (!treenode[i]->tip) {
1953       p = treenode[i]->next;
1954       while (p != treenode[i]) {
1955         p->bottom = false;
1956         p = p->next;
1957       }
1958     }
1959   }
1960 } /* clearbottom */
1961
1962
1963 void collabranch(node *collapfrom, node *tempfrom, node *tempto)
1964 { /* collapse branch from collapfrom */
1965   long i, j, b, largest, descsteps;
1966   boolean done;
1967
1968   for (i = 0; i < endsite; i++) {
1969     descsteps = 0;
1970     for (j = (long)A; j <= (long)O; j++) {
1971       b = 1 << j;
1972       if ((descsteps == 0) && (collapfrom->base[i] & b))
1973         descsteps = tempfrom->oldnumsteps[i]
1974                      - (collapfrom->numdesc - collapfrom->numnuc[i][j])
1975                        * weight[i];
1976     }
1977     done = false;
1978     for (j = (long)A; j <= (long)O; j++) {
1979       b = 1 << j;
1980       if (!done && (tempto->base[i] & b)) {
1981         descsteps += (tempto->numsteps[i]
1982                       - (tempto->numdesc - collapfrom->numdesc
1983                         - tempto->numnuc[i][j]) * weight[i]);
1984         done = true;
1985       }
1986     }
1987     for (j = (long)A; j <= (long)O; j++)
1988       tempto->numnuc[i][j] += collapfrom->numnuc[i][j];
1989     largest = getlargest(tempto->numnuc[i]);
1990     tempto->base[i] = 0;
1991     for (j = (long)A; j <= (long)O; j++) {
1992       if (tempto->numnuc[i][j] == largest)
1993         tempto->base[i] |= (1 << j);
1994     }
1995     tempto->numsteps[i] = (tempto->numdesc - largest) * weight[i] + descsteps;
1996   }
1997 } /* collabranch */
1998
1999
2000 boolean allcommonbases(node *a, node *b, boolean *allsame)
2001 {  /* see if bases are common at all sites for nodes a and b */
2002   long i;
2003   boolean allcommon;
2004
2005   allcommon = true;
2006   *allsame = true;
2007   for (i = 0; i < endsite; i++) {
2008     if ((a->base[i] & b->base[i]) == 0)
2009       allcommon = false;
2010     else if (a->base[i] != b->base[i])
2011       *allsame = false;
2012   }
2013   return allcommon;
2014 } /* allcommonbases */
2015
2016
2017 void findbottom(node *p, node **bottom)
2018 { /* find a node with field bottom set at node p */
2019   node *q;
2020
2021   if (p->bottom)
2022     *bottom = p;
2023   else {
2024     q = p->next;
2025     while(!q->bottom && q != p)
2026       q = q->next;
2027     *bottom = q;
2028   }
2029 } /* findbottom */
2030
2031
2032 boolean moresteps(node *a, node *b)
2033 {  /* see if numsteps of node a exceeds those of node b */
2034   long i;
2035
2036   for (i = 0; i < endsite; i++)
2037     if (a->numsteps[i] > b->numsteps[i])
2038       return true;
2039   return false;
2040 } /* moresteps */
2041
2042
2043 boolean passdown(node *desc, node *parent, node *start, node *below,
2044                         node *item, node *added, node *total, node *tempdsc,
2045             node *tempprt, boolean multf)
2046 { /* track down to node start to see if an ancestor branch can be collapsed */
2047   node *temp;
2048   boolean done, allsame;
2049
2050   done = (parent == start);
2051   while (!done) {
2052     desc = parent;
2053     findbottom(parent->back, &parent);
2054     if (multf && start == below && parent == below)
2055       parent = added;
2056     memcpy(tempdsc->base, tempprt->base, endsite*sizeof(long));
2057     memcpy(tempdsc->numsteps, tempprt->numsteps, endsite*sizeof(long));
2058     memcpy(tempdsc->oldbase, desc->base, endsite*sizeof(long));
2059     memcpy(tempdsc->oldnumsteps, desc->numsteps, endsite*sizeof(long));
2060     memcpy(tempprt->base, parent->base, endsite*sizeof(long));
2061     memcpy(tempprt->numsteps, parent->numsteps, endsite*sizeof(long));
2062     memcpy(tempprt->numnuc, parent->numnuc, endsite*sizeof(nucarray));
2063     tempprt->numdesc = parent->numdesc;
2064     multifillin(tempprt, tempdsc, 0);
2065     if (!allcommonbases(tempprt, parent, &allsame))
2066       return false;
2067     else if (moresteps(tempprt, parent))
2068       return false;
2069     else if (allsame)
2070       return true;
2071     if (parent == added)
2072       parent = below;
2073     done = (parent == start);
2074     if (done && ((start == item) || (!multf && start == below))) {
2075       memcpy(tempdsc->base, tempprt->base, endsite*sizeof(long));
2076       memcpy(tempdsc->numsteps, tempprt->numsteps, endsite*sizeof(long));
2077       memcpy(tempdsc->oldbase, start->base, endsite*sizeof(long));
2078       memcpy(tempdsc->oldnumsteps, start->numsteps, endsite*sizeof(long));
2079       multifillin(added, tempdsc, 0);
2080       tempprt = added;
2081     }
2082   }
2083   temp = tempdsc;
2084   if (start == below || start == item)
2085     fillin(temp, tempprt, below->back);
2086   else
2087     fillin(temp, tempprt, added);
2088   return !moresteps(temp, total);
2089 } /* passdown */
2090
2091
2092 boolean trycollapdesc(node *desc, node *parent, node *start,
2093                         node *below, node *item, node *added, node *total,
2094             node *tempdsc, node *tempprt, boolean multf, long *zeros)
2095   { /* see if branch between nodes desc and parent can be collapsed */
2096   boolean allsame;
2097
2098   if (desc->numdesc == 1)
2099     return true;
2100   if (multf && start == below && parent == below)
2101     parent = added;
2102   memcpy(tempdsc->base, zeros, endsite*sizeof(long));
2103   memcpy(tempdsc->numsteps, zeros, endsite*sizeof(long));
2104   memcpy(tempdsc->oldbase, desc->base, endsite*sizeof(long));
2105   memcpy(tempdsc->oldnumsteps, desc->numsteps, endsite*sizeof(long));
2106   memcpy(tempprt->base, parent->base, endsite*sizeof(long));
2107   memcpy(tempprt->numsteps, parent->numsteps, endsite*sizeof(long));
2108   memcpy(tempprt->numnuc, parent->numnuc, endsite*sizeof(nucarray));
2109   tempprt->numdesc = parent->numdesc - 1;
2110   multifillin(tempprt, tempdsc, -1);
2111   tempprt->numdesc += desc->numdesc;
2112   collabranch(desc, tempdsc, tempprt);
2113   if (!allcommonbases(tempprt, parent, &allsame) ||
2114         moresteps(tempprt, parent)) {
2115     if (parent != added) {
2116       desc->collapse = nocollap;
2117       parent->collapse = nocollap;
2118     }
2119     return false;
2120   } else if (allsame) {
2121     if (parent != added) {
2122       desc->collapse = tocollap;
2123       parent->collapse = tocollap;
2124     }
2125     return true;
2126   }
2127   if (parent == added)
2128     parent = below;
2129   if ((start == item && parent == item) ||
2130         (!multf && start == below && parent == below)) {
2131     memcpy(tempdsc->base, tempprt->base, endsite*sizeof(long));
2132     memcpy(tempdsc->numsteps, tempprt->numsteps, endsite*sizeof(long));
2133     memcpy(tempdsc->oldbase, start->base, endsite*sizeof(long));
2134     memcpy(tempdsc->oldnumsteps, start->numsteps, endsite*sizeof(long));
2135     memcpy(tempprt->base, added->base, endsite*sizeof(long));
2136     memcpy(tempprt->numsteps, added->numsteps, endsite*sizeof(long));
2137     memcpy(tempprt->numnuc, added->numnuc, endsite*sizeof(nucarray));
2138     tempprt->numdesc = added->numdesc;
2139     multifillin(tempprt, tempdsc, 0);
2140     if (!allcommonbases(tempprt, added, &allsame))
2141       return false;
2142     else if (moresteps(tempprt, added))
2143       return false;
2144     else if (allsame)
2145       return true;
2146   }
2147   return passdown(desc, parent, start, below, item, added, total, tempdsc,
2148                     tempprt, multf);
2149 } /* trycollapdesc */
2150
2151
2152 void setbottom(node *p)
2153 { /* set field bottom at node p */
2154   node *q;
2155
2156   p->bottom = true;
2157   q = p->next;
2158   do {
2159     q->bottom = false;
2160     q = q->next;
2161   } while (q != p);
2162 } /* setbottom */
2163
2164 boolean zeroinsubtree(node *subtree, node *start, node *below, node *item,
2165                         node *added, node *total, node *tempdsc, node *tempprt,
2166                         boolean multf, node* root, long *zeros)
2167 { /* sees if subtree contains a zero length branch */
2168   node *p;
2169
2170   if (!subtree->tip) {
2171     setbottom(subtree);
2172     p = subtree->next;
2173     do {
2174       if (p->back && !p->back->tip &&
2175          !((p->back->collapse == nocollap) && (subtree->collapse == nocollap))
2176            && (subtree->numdesc != 1)) {
2177         if ((p->back->collapse == tocollap) && (subtree->collapse == tocollap)
2178             && multf && (subtree != below))
2179           return true;
2180         /* when root->numdesc == 2
2181          * there is no mandatory step at the root,
2182          * instead of checking at the root we check around it
2183          * we only need to check p because the first if
2184          * statement already gets rid of it for the subtree */
2185         else if ((p->back->index != root->index || root->numdesc > 2) &&
2186             trycollapdesc(p->back, subtree, start, below, item, added, total,
2187                 tempdsc, tempprt, multf, zeros))
2188           return true;
2189         else if ((p->back->index == root->index && root->numdesc == 2) &&
2190             !(root->next->back->tip) && !(root->next->next->back->tip) &&
2191             trycollapdesc(root->next->back, root->next->next->back, start,
2192                 below, item,added, total, tempdsc, tempprt, multf, zeros))
2193           return true;
2194       }
2195       p = p->next;
2196     } while (p != subtree);
2197     p = subtree->next;
2198     do {
2199       if (p->back && !p->back->tip) {
2200         if (zeroinsubtree(p->back, start, below, item, added, total,
2201                             tempdsc, tempprt, multf, root, zeros))
2202           return true;
2203       }
2204       p = p->next;
2205     } while (p != subtree);
2206   }
2207   return false;
2208 } /* zeroinsubtree */
2209
2210
2211 boolean collapsible(node *item, node *below, node *temp, node *temp1,
2212                         node *tempdsc, node *tempprt, node *added, node *total,
2213             boolean multf, node *root, long *zeros, pointarray treenode)
2214 {
2215   /* sees if any branch can be collapsed */
2216   node *belowbk;
2217   boolean allsame;
2218
2219   if (multf) {
2220     memcpy(tempdsc->base, item->base, endsite*sizeof(long));
2221     memcpy(tempdsc->numsteps, item->numsteps, endsite*sizeof(long));
2222     memcpy(tempdsc->oldbase, zeros, endsite*sizeof(long));
2223     memcpy(tempdsc->oldnumsteps, zeros, endsite*sizeof(long));
2224     memcpy(added->base, below->base, endsite*sizeof(long));
2225     memcpy(added->numsteps, below->numsteps, endsite*sizeof(long));
2226     memcpy(added->numnuc, below->numnuc, endsite*sizeof(nucarray));
2227     added->numdesc = below->numdesc + 1;
2228     multifillin(added, tempdsc, 1);
2229   } else {
2230     fillin(added, item, below);
2231     added->numdesc = 2;
2232   }
2233   fillin(total, added, below->back);
2234   clearbottom(treenode);
2235   if (below->back) {
2236     if (zeroinsubtree(below->back, below->back, below, item, added, total,
2237                         tempdsc, tempprt, multf, root, zeros))
2238       return true;
2239   }
2240   if (multf) {
2241     if (zeroinsubtree(below, below, below, item, added, total,
2242                        tempdsc, tempprt, multf, root, zeros))
2243       return true;
2244   } else if (!below->tip) {
2245     if (zeroinsubtree(below, below, below, item, added, total,
2246                         tempdsc, tempprt, multf, root, zeros))
2247       return true;
2248   }
2249   if (!item->tip) {
2250     if (zeroinsubtree(item, item, below, item, added, total,
2251                         tempdsc, tempprt, multf, root, zeros))
2252       return true;
2253   }
2254   if (multf && below->back && !below->back->tip) {
2255     memcpy(tempdsc->base, zeros, endsite*sizeof(long));
2256     memcpy(tempdsc->numsteps, zeros, endsite*sizeof(long));
2257     memcpy(tempdsc->oldbase, added->base, endsite*sizeof(long));
2258     memcpy(tempdsc->oldnumsteps, added->numsteps, endsite*sizeof(long));
2259     if (below->back == treenode[below->back->index - 1])
2260       belowbk = below->back->next;
2261     else
2262       belowbk = treenode[below->back->index - 1];
2263     memcpy(tempprt->base, belowbk->base, endsite*sizeof(long));
2264     memcpy(tempprt->numsteps, belowbk->numsteps, endsite*sizeof(long));
2265     memcpy(tempprt->numnuc, belowbk->numnuc, endsite*sizeof(nucarray));
2266     tempprt->numdesc = belowbk->numdesc - 1;
2267     multifillin(tempprt, tempdsc, -1);
2268     tempprt->numdesc += added->numdesc;
2269     collabranch(added, tempdsc, tempprt);
2270     if (!allcommonbases(tempprt, belowbk, &allsame))
2271       return false;
2272     else if (allsame && !moresteps(tempprt, belowbk))
2273       return true;
2274     else if (belowbk->back) {
2275       fillin(temp, tempprt, belowbk->back);
2276       fillin(temp1, belowbk, belowbk->back);
2277       return !moresteps(temp, temp1);
2278     }
2279   }
2280   return false;
2281 } /* collapsible */
2282
2283
2284 void replaceback(node **oldback, node *item, node *forknode,
2285                         node **grbg, long *zeros)
2286 { /* replaces back node of item with another */
2287   node *p;
2288
2289   p = forknode;
2290   while (p->next->back != item)
2291     p = p->next;
2292   *oldback = p->next;
2293   gnutreenode(grbg, &p->next, forknode->index, endsite, zeros);
2294   p->next->next = (*oldback)->next;
2295   p->next->back = (*oldback)->back;
2296   p->next->back->back = p->next;
2297   (*oldback)->next = (*oldback)->back = NULL;
2298 } /* replaceback */
2299
2300
2301 void putback(node *oldback, node *item, node *forknode, node **grbg)
2302 { /* restores node to back of item */
2303   node *p, *q;
2304
2305   p = forknode;
2306   while (p->next != item->back)
2307     p = p->next;
2308   q = p->next;
2309   oldback->next = p->next->next;
2310   p->next = oldback;
2311   oldback->back = item;
2312   item->back = oldback;
2313   oldback->index = forknode->index;
2314   chucktreenode(grbg, q);
2315 } /* putback */
2316
2317
2318 void savelocrearr(node *item, node *forknode, node *below, node *tmp,
2319         node *tmp1, node *tmp2, node *tmp3, node *tmprm, node *tmpadd,
2320         node **root, long maxtrees, long *nextree, boolean multf,
2321         boolean bestever, boolean *saved, long *place,
2322         bestelm *bestrees, pointarray treenode, node **grbg,
2323         long *zeros)
2324 { /* saves tied or better trees during local rearrangements by removing
2325      item from forknode and adding to below */
2326   node *other, *otherback = NULL, *oldfork, *nufork, *oldback;
2327   long pos;
2328   boolean found, collapse;
2329
2330   if (forknode->numdesc == 2) {
2331     findbelow(&other, item, forknode);
2332     otherback = other->back;
2333     oldback = NULL;
2334   } else {
2335     other = NULL;
2336     replaceback(&oldback, item, forknode, grbg, zeros);
2337   }
2338   re_move(item, &oldfork, root, false, treenode, grbg, zeros);
2339   if (!multf)
2340     getnufork(&nufork, grbg, treenode, zeros);
2341   else
2342     nufork = NULL;
2343   addnsave(below, item, nufork, root, grbg, multf, treenode, place, zeros);
2344   pos = 0;
2345   findtree(&found, &pos, *nextree, place, bestrees);
2346   if (other) {
2347     add(other, item, oldfork, root, false, treenode, grbg, zeros);
2348     if (otherback->back != other)
2349       flipnodes(item, other);
2350   } else
2351     add(forknode, item, NULL, root, false, treenode, grbg, zeros);
2352   *saved = false;
2353   if (found) {
2354     if (oldback)
2355       putback(oldback, item, forknode, grbg);
2356   } else {
2357     if (oldback)
2358       chucktreenode(grbg, oldback);
2359     re_move(item, &oldfork, root, true, treenode, grbg, zeros);
2360     collapse = collapsible(item, below, tmp, tmp1, tmp2, tmp3, tmprm,
2361                      tmpadd, multf, *root, zeros, treenode);
2362     if (!collapse) {
2363       if (bestever)
2364         addbestever(&pos, nextree, maxtrees, collapse, place, bestrees);
2365       else
2366         addtiedtree(pos, nextree, maxtrees, collapse, place, bestrees);
2367     }
2368     if (other)
2369       add(other, item, oldfork, root, true, treenode, grbg, zeros);
2370     else
2371       add(forknode, item, NULL, root, true, treenode, grbg, zeros);
2372     *saved = !collapse;
2373   }
2374 } /* savelocrearr */
2375
2376
2377 void clearvisited(pointarray treenode)
2378 {
2379   /* clears boolean visited at a node */
2380   long i;
2381   node *p;
2382
2383   for (i = 0; i < nonodes; i++) {
2384     treenode[i]->visited = false;
2385     if (!treenode[i]->tip) {
2386       p = treenode[i]->next;
2387       while (p != treenode[i]) {
2388         p->visited = false;
2389         p = p->next;
2390       }
2391     }
2392   }
2393 } /* clearvisited */
2394
2395
2396 void hyprint(long b1, long b2, struct LOC_hyptrav *htrav,
2397                         pointarray treenode, Char *basechar)
2398 {
2399   /* print out states in sites b1 through b2 at node */
2400   long i, j, k, n;
2401   boolean dot;
2402   bases b;
2403
2404   if (htrav->bottom) {
2405     if (!outgropt)
2406       fprintf(outfile, "       ");
2407     else
2408       fprintf(outfile, "root   ");
2409   } else
2410     fprintf(outfile, "%4ld   ", htrav->r->back->index - spp);
2411   if (htrav->r->tip) {
2412     for (i = 0; i < nmlngth; i++)
2413       putc(nayme[htrav->r->index - 1][i], outfile);
2414   } else
2415     fprintf(outfile, "%4ld      ", htrav->r->index - spp);
2416   if (htrav->bottom)
2417     fprintf(outfile, "          ");
2418   else if (htrav->nonzero)
2419     fprintf(outfile, "   yes    ");
2420   else if (htrav->maybe)
2421     fprintf(outfile, "  maybe   ");
2422   else
2423     fprintf(outfile, "   no     ");
2424   for (i = b1; i <= b2; i++) {
2425     j = location[ally[i - 1] - 1];
2426     htrav->tempset = htrav->r->base[j - 1];
2427     htrav->anc = htrav->hypset[j - 1];
2428     if (!htrav->bottom)
2429       htrav->anc = treenode[htrav->r->back->index - 1]->base[j - 1];
2430     dot = dotdiff && (htrav->tempset == htrav->anc && !htrav->bottom);
2431     if (dot)
2432       putc('.', outfile);
2433     else if (htrav->tempset == (1 << A))
2434       putc('A', outfile);
2435     else if (htrav->tempset == (1 << C))
2436       putc('C', outfile);
2437     else if (htrav->tempset == (1 << G))
2438       putc('G', outfile);
2439     else if (htrav->tempset == (1 << T))
2440       putc('T', outfile);
2441     else if (htrav->tempset == (1 << O))
2442       putc('-', outfile);
2443     else {
2444       k = 1;
2445       n = 0;
2446       for (b = A; b <= O; b = b + 1) {
2447         if (((1 << b) & htrav->tempset) != 0)
2448           n += k;
2449         k += k;
2450       }
2451       putc(basechar[n - 1], outfile);
2452     }
2453     if (i % 10 == 0)
2454       putc(' ', outfile);
2455   }
2456   putc('\n', outfile);
2457 }  /* hyprint */
2458
2459
2460 void gnubase(gbases **p, gbases **garbage, long endsite)
2461 {
2462   /* this and the following are do-it-yourself garbage collectors.
2463      Make a new node or pull one off the garbage list */
2464   if (*garbage != NULL) {
2465     *p = *garbage;
2466     *garbage = (*garbage)->next;
2467   } else {
2468     *p = (gbases *)Malloc(sizeof(gbases));
2469     (*p)->base = (baseptr)Malloc(endsite*sizeof(long));
2470   }
2471   (*p)->next = NULL;
2472 }  /* gnubase */
2473
2474
2475 void chuckbase(gbases *p, gbases **garbage)
2476 {
2477   /* collect garbage on p -- put it on front of garbage list */
2478   p->next = *garbage;
2479   *garbage = p;
2480 }  /* chuckbase */
2481
2482
2483 void hyptrav(node *r_, long *hypset_, long b1, long b2, boolean bottom_,
2484                         pointarray treenode, gbases **garbage, Char *basechar)
2485 {
2486   /*  compute, print out states at one interior node */
2487   struct LOC_hyptrav Vars;
2488   long i, j, k;
2489   long largest;
2490   gbases *ancset;
2491   nucarray *tempnuc;
2492   node *p, *q;
2493
2494   Vars.bottom = bottom_;
2495   Vars.r = r_;
2496   Vars.hypset = hypset_;
2497   gnubase(&ancset, garbage, endsite);
2498   tempnuc = (nucarray *)Malloc(endsite*sizeof(nucarray));
2499   Vars.maybe = false;
2500   Vars.nonzero = false;
2501   if (!Vars.r->tip)
2502     zeronumnuc(Vars.r, endsite);
2503   for (i = b1 - 1; i < b2; i++) {
2504     j = location[ally[i] - 1];
2505     Vars.anc = Vars.hypset[j - 1];
2506     if (!Vars.r->tip) {
2507       p = Vars.r->next;
2508       for (k = (long)A; k <= (long)O; k++)
2509         if (Vars.anc & (1 << k))
2510           Vars.r->numnuc[j - 1][k]++;
2511       do {
2512         for (k = (long)A; k <= (long)O; k++)
2513           if (p->back->base[j - 1] & (1 << k))
2514             Vars.r->numnuc[j - 1][k]++;
2515         p = p->next;
2516       } while (p != Vars.r);
2517       largest = getlargest(Vars.r->numnuc[j - 1]);
2518       Vars.tempset = 0;
2519       for (k = (long)A; k <= (long)O; k++) {
2520         if (Vars.r->numnuc[j - 1][k] == largest)
2521           Vars.tempset |= (1 << k);
2522       }
2523       Vars.r->base[j - 1] = Vars.tempset;
2524     }
2525     if (!Vars.bottom)
2526       Vars.anc = treenode[Vars.r->back->index - 1]->base[j - 1];
2527     Vars.nonzero = (Vars.nonzero || (Vars.r->base[j - 1] & Vars.anc) == 0);
2528     Vars.maybe = (Vars.maybe || Vars.r->base[j - 1] != Vars.anc);
2529   }
2530   hyprint(b1, b2, &Vars, treenode, basechar);
2531   Vars.bottom = false;
2532   if (!Vars.r->tip) {
2533     memcpy(tempnuc, Vars.r->numnuc, endsite*sizeof(nucarray));
2534     q = Vars.r->next;
2535     do {
2536       memcpy(Vars.r->numnuc, tempnuc, endsite*sizeof(nucarray));
2537       for (i = b1 - 1; i < b2; i++) {
2538         j = location[ally[i] - 1];
2539         for (k = (long)A; k <= (long)O; k++)
2540           if (q->back->base[j - 1] & (1 << k))
2541             Vars.r->numnuc[j - 1][k]--;
2542         largest = getlargest(Vars.r->numnuc[j - 1]);
2543         ancset->base[j - 1] = 0;
2544         for (k = (long)A; k <= (long)O; k++)
2545           if (Vars.r->numnuc[j - 1][k] == largest)
2546             ancset->base[j - 1] |= (1 << k);
2547         if (!Vars.bottom)
2548           Vars.anc = ancset->base[j - 1];
2549       }
2550       hyptrav(q->back, ancset->base, b1, b2, Vars.bottom,
2551                 treenode, garbage, basechar);
2552       q = q->next;
2553     } while (q != Vars.r);
2554   }
2555   chuckbase(ancset, garbage);
2556 }  /* hyptrav */
2557
2558
2559 void hypstates(long chars, node *root, pointarray treenode,
2560                         gbases **garbage, Char *basechar)
2561 {
2562   /* fill in and describe states at interior nodes */
2563   /* used in dnacomp, dnapars, & dnapenny */
2564   long i, n;
2565   baseptr nothing;
2566
2567   fprintf(outfile, "\nFrom    To     Any Steps?    State at upper node\n");
2568   fprintf(outfile, "                            ");
2569   if (dotdiff)
2570     fprintf(outfile, " ( . means same as in the node below it on tree)\n");
2571   nothing = (baseptr)Malloc(endsite*sizeof(long));
2572   for (i = 0; i < endsite; i++)
2573     nothing[i] = 0;
2574   for (i = 1; i <= ((chars - 1) / 40 + 1); i++) {
2575     putc('\n', outfile);
2576     n = i * 40;
2577     if (n > chars)
2578       n = chars;
2579     hyptrav(root, nothing, i * 40 - 39, n, true, treenode, garbage, basechar);
2580   }
2581   free(nothing);
2582 }  /* hypstates */
2583
2584
2585 void initbranchlen(node *p)
2586 {
2587   node *q;
2588
2589   p->v = 0.0;
2590   if (p->back)
2591     p->back->v = 0.0;
2592   if (p->tip)
2593     return;
2594   q = p->next;
2595   while (q != p) {
2596     initbranchlen(q->back);
2597     q = q->next;
2598   }
2599   q = p->next;
2600   while (q != p) {
2601     q->v = 0.0;
2602     q = q->next;
2603   }
2604 } /* initbranchlen */
2605
2606
2607 void initmin(node *p, long sitei, boolean internal)
2608 {
2609   long i;
2610
2611   if (internal) {
2612     for (i = (long)A; i <= (long)O; i++) {
2613       p->cumlengths[i] = 0;
2614       p->numreconst[i] = 1;
2615     }
2616   } else {
2617     for (i = (long)A; i <= (long)O; i++) {
2618       if (p->base[sitei - 1] & (1 << i)) {
2619         p->cumlengths[i] = 0;
2620         p->numreconst[i] = 1;
2621       } else {
2622         p->cumlengths[i] = -1;
2623         p->numreconst[i] = 0;
2624       }
2625     }
2626   }
2627 } /* initmin */
2628
2629
2630 void initbase(node *p, long sitei)
2631 {
2632   /* traverse tree to initialize base at internal nodes */
2633   node *q;
2634   long i, largest;
2635
2636   if (p->tip)
2637     return;
2638   q = p->next;
2639   while (q != p) {
2640     if (q->back) {
2641       memcpy(q->numnuc, p->numnuc, endsite*sizeof(nucarray));
2642       for (i = (long)A; i <= (long)O; i++) {
2643         if (q->back->base[sitei - 1] & (1 << i))
2644           q->numnuc[sitei - 1][i]--;
2645       }
2646       if (p->back) {
2647         for (i = (long)A; i <= (long)O; i++) {
2648           if (p->back->base[sitei - 1] & (1 << i))
2649             q->numnuc[sitei - 1][i]++;
2650         }
2651       }
2652       largest = getlargest(q->numnuc[sitei - 1]);
2653       q->base[sitei - 1] = 0;
2654       for (i = (long)A; i <= (long)O; i++) {
2655         if (q->numnuc[sitei - 1][i] == largest)
2656           q->base[sitei - 1] |= (1 << i);
2657       }
2658     }
2659     q = q->next;
2660   }
2661   q = p->next;
2662   while (q != p) {
2663     initbase(q->back, sitei);
2664     q = q->next;
2665   }
2666 } /* initbase */
2667
2668
2669 void inittreetrav(node *p, long sitei)
2670 {
2671   /* traverse tree to clear boolean initialized and set up base */
2672   node *q;
2673
2674   if (p->tip) {
2675     initmin(p, sitei, false);
2676     p->initialized = true;
2677     return;
2678   }
2679   q = p->next;
2680   while (q != p) {
2681     inittreetrav(q->back, sitei);
2682     q = q->next;
2683   }
2684   initmin(p, sitei, true);
2685   p->initialized = false;
2686   q = p->next;
2687   while (q != p) {
2688     initmin(q, sitei, true);
2689     q->initialized = false;
2690     q = q->next;
2691   }
2692 } /* inittreetrav */
2693
2694
2695 void compmin(node *p, node *desc)
2696 {
2697   /* computes minimum lengths up to p */
2698   long i, j, minn, cost, desclen, descrecon=0, maxx;
2699
2700   maxx = 10 * spp;
2701   for (i = (long)A; i <= (long)O; i++) {
2702     minn = maxx;
2703     for (j = (long)A; j <= (long)O; j++) {
2704       if (transvp) {
2705         if (
2706                (
2707                 ((i == (long)A) || (i == (long)G))
2708              && ((j == (long)A) || (j == (long)G))
2709                )
2710             || (
2711                 ((j == (long)C) || (j == (long)T))
2712              && ((i == (long)C) || (i == (long)T))
2713                )
2714            )
2715           cost = 0;
2716         else
2717           cost = 1;
2718       } else {
2719         if (i == j)
2720           cost = 0;
2721         else
2722           cost = 1;
2723       }
2724       if (desc->cumlengths[j] == -1) {
2725         desclen = maxx;
2726       } else {
2727         desclen = desc->cumlengths[j];
2728       }
2729       if (minn > cost + desclen) {
2730         minn = cost + desclen;
2731         descrecon = 0;
2732       }
2733       if (minn == cost + desclen) {
2734         descrecon += desc->numreconst[j];
2735       }
2736     }
2737     p->cumlengths[i] += minn;
2738     p->numreconst[i] *= descrecon;
2739   }
2740   p->initialized = true;
2741 } /* compmin */
2742
2743
2744 void minpostorder(node *p, pointarray treenode)
2745 {
2746   /* traverses an n-ary tree, computing minimum steps at each node */
2747   node *q;
2748
2749   if (p->tip) {
2750     return;
2751   }
2752   q = p->next;
2753   while (q != p) {
2754     if (q->back)
2755       minpostorder(q->back, treenode);
2756     q = q->next;
2757   }
2758   if (!p->initialized) {
2759     q = p->next;
2760     while (q != p) {
2761       if (q->back)
2762         compmin(p, q->back);
2763       q = q->next;
2764     }
2765   }
2766 }  /* minpostorder */
2767
2768
2769 void branchlength(node *subtr1, node *subtr2, double *brlen,
2770                         pointarray treenode)
2771 {
2772   /* computes a branch length between two subtrees for a given site */
2773   long i, j, minn, cost, nom, denom;
2774   node *temp;
2775
2776   if (subtr1->tip) {
2777     temp = subtr1;
2778     subtr1 = subtr2;
2779     subtr2 = temp;
2780   }
2781   if (subtr1->index == outgrno) {
2782     temp = subtr1;
2783     subtr1 = subtr2;
2784     subtr2 = temp;
2785   }
2786   minpostorder(subtr1, treenode);
2787   minpostorder(subtr2, treenode);
2788   minn = 10 * spp;
2789   nom = 0;
2790   denom = 0;
2791   for (i = (long)A; i <= (long)O; i++) {
2792     for (j = (long)A; j <= (long)O; j++) {
2793       if (transvp) {
2794         if (
2795                (
2796                 ((i == (long)A) || (i == (long)G))
2797              && ((j == (long)A) || (j == (long)G))
2798                )
2799             || (
2800                 ((j == (long)C) || (j == (long)T))
2801              && ((i == (long)C) || (i == (long)T))
2802                )
2803            )
2804           cost = 0;
2805         else
2806           cost = 1;
2807       } else {
2808         if (i == j)
2809           cost = 0;
2810         else
2811           cost = 1;
2812       }
2813       if (subtr1->cumlengths[i] != -1 && (subtr2->cumlengths[j] != -1)) {
2814         if (subtr1->cumlengths[i] + cost + subtr2->cumlengths[j] < minn) {
2815           minn = subtr1->cumlengths[i] + cost + subtr2->cumlengths[j];
2816           nom = 0;
2817           denom = 0;
2818         }
2819         if (subtr1->cumlengths[i] + cost + subtr2->cumlengths[j] == minn) {
2820           nom += subtr1->numreconst[i] * subtr2->numreconst[j] * cost;
2821           denom += subtr1->numreconst[i] * subtr2->numreconst[j];
2822         }
2823       }
2824     }
2825   }
2826   *brlen = (double)nom/(double)denom;
2827 } /* branchlength */
2828
2829
2830 void printbranchlengths(node *p)
2831 {
2832   node *q;
2833   long i;
2834
2835   if (p->tip)
2836     return;
2837   q = p->next;
2838   do {
2839     fprintf(outfile, "%6ld      ",q->index - spp);
2840     if (q->back->tip) {
2841       for (i = 0; i < nmlngth; i++)
2842         putc(nayme[q->back->index - 1][i], outfile);
2843     } else
2844       fprintf(outfile, "%6ld    ", q->back->index - spp);
2845     fprintf(outfile, "   %f\n",q->v);
2846     if (q->back)
2847       printbranchlengths(q->back);
2848     q = q->next;
2849   } while (q != p);
2850 } /* printbranchlengths */
2851
2852
2853 void branchlentrav(node *p, node *root, long sitei, long chars,
2854                         double *brlen, pointarray treenode)
2855   {
2856   /*  traverses the tree computing tree length at each branch */
2857   node *q;
2858
2859   if (p->tip)
2860     return;
2861   if (p->index == outgrno)
2862     p = p->back;
2863   q = p->next;
2864   do {
2865     if (q->back) {
2866       branchlength(q, q->back, brlen, treenode);
2867       q->v += ((weight[sitei - 1] / 10.0) * (*brlen)/chars);
2868       q->back->v += ((weight[sitei - 1] / 10.0) * (*brlen)/chars);
2869       if (!q->back->tip)
2870         branchlentrav(q->back, root, sitei, chars, brlen, treenode);
2871     }
2872     q = q->next;
2873   } while (q != p);
2874 }  /* branchlentrav */
2875
2876
2877 void treelength(node *root, long chars, pointarray treenode)
2878   {
2879   /*  calls branchlentrav at each site */
2880   long sitei;
2881   double trlen;
2882
2883   initbranchlen(root);
2884   for (sitei = 1; sitei <= endsite; sitei++) {
2885     trlen = 0.0;
2886     initbase(root, sitei);
2887     inittreetrav(root, sitei);
2888     branchlentrav(root, root, sitei, chars, &trlen, treenode);
2889   }
2890 } /* treelength */
2891
2892
2893 void coordinates(node *p, long *tipy, double f, long *fartemp)
2894 {
2895   /* establishes coordinates of nodes for display without lengths */
2896   node *q, *first, *last;
2897   node *mid1 = NULL, *mid2 = NULL;
2898   long numbranches, numb2;
2899
2900   if (p->tip) {
2901     p->xcoord = 0;
2902     p->ycoord = *tipy;
2903     p->ymin = *tipy;
2904     p->ymax = *tipy;
2905     (*tipy) += down;
2906     return;
2907   }
2908   numbranches = 0;
2909   q = p->next;
2910   do {
2911     coordinates(q->back, tipy, f, fartemp);
2912     numbranches += 1;
2913     q = q->next;
2914   } while (p != q);
2915   first = p->next->back;
2916   q = p->next;
2917   while (q->next != p)
2918     q = q->next;
2919   last = q->back;
2920   numb2 = 1;
2921   q = p->next;
2922   while (q != p) {
2923     if (numb2 == (long)(numbranches + 1)/2)
2924       mid1 = q->back;
2925     if (numb2 == (long)(numbranches/2 + 1))
2926       mid2 = q->back;
2927     numb2 += 1;
2928     q = q->next;
2929   }
2930   p->xcoord = (long)((double)(last->ymax - first->ymin) * f);
2931   p->ycoord = (long)((mid1->ycoord + mid2->ycoord) / 2);
2932   p->ymin = first->ymin;
2933   p->ymax = last->ymax;
2934   if (p->xcoord > *fartemp)
2935     *fartemp = p->xcoord;
2936 }  /* coordinates */
2937
2938
2939 void drawline(long i, double scale, node *root)
2940 {
2941   /* draws one row of the tree diagram by moving up tree */
2942   node *p, *q, *r, *first =NULL, *last =NULL;
2943   long n, j;
2944   boolean extra, done, noplus;
2945
2946   p = root;
2947   q = root;
2948   extra = false;
2949   noplus = false;
2950   if (i == (long)p->ycoord && p == root) {
2951     if (p->index - spp >= 10)
2952       fprintf(outfile, " %2ld", p->index - spp);
2953     else
2954       fprintf(outfile, "  %ld", p->index - spp);
2955     extra = true;
2956     noplus = true;
2957   } else
2958     fprintf(outfile, "  ");
2959   do {
2960     if (!p->tip) {
2961       r = p->next;
2962       done = false;
2963       do {
2964         if (i >= r->back->ymin && i <= r->back->ymax) {
2965           q = r->back;
2966           done = true;
2967         }
2968         r = r->next;
2969       } while (!(done || r == p));
2970       first = p->next->back;
2971       r = p->next;
2972       while (r->next != p)
2973         r = r->next;
2974       last = r->back;
2975     }
2976     done = (p == q);
2977     n = (long)(scale * (p->xcoord - q->xcoord) + 0.5);
2978     if (n < 3 && !q->tip)
2979       n = 3;
2980     if (extra) {
2981       n--;
2982       extra = false;
2983     }
2984     if ((long)q->ycoord == i && !done) {
2985       if (noplus) {
2986         putc('-', outfile);
2987         noplus = false;
2988       }
2989       else
2990         putc('+', outfile);
2991       if (!q->tip) {
2992         for (j = 1; j <= n - 2; j++)
2993           putc('-', outfile);
2994         if (q->index - spp >= 10)
2995           fprintf(outfile, "%2ld", q->index - spp);
2996         else
2997           fprintf(outfile, "-%ld", q->index - spp);
2998         extra = true;
2999         noplus = true;
3000       } else {
3001         for (j = 1; j < n; j++)
3002           putc('-', outfile);
3003       }
3004     } else if (!p->tip) {
3005       if ((long)last->ycoord > i && (long)first->ycoord < i
3006             && i != (long)p->ycoord) {
3007         putc('!', outfile);
3008         for (j = 1; j < n; j++)
3009           putc(' ', outfile);
3010       } else {
3011         for (j = 1; j <= n; j++)
3012           putc(' ', outfile);
3013       }
3014       noplus = false;
3015     } else {
3016       for (j = 1; j <= n; j++)
3017         putc(' ', outfile);
3018       noplus = false;
3019     }
3020     if (p != q)
3021       p = q;
3022   } while (!done);
3023   if ((long)p->ycoord == i && p->tip) {
3024     for (j = 0; j < nmlngth; j++)
3025       putc(nayme[p->index - 1][j], outfile);
3026   }
3027   putc('\n', outfile);
3028 }  /* drawline */
3029
3030
3031 void printree(node *root, double f)
3032 {
3033   /* prints out diagram of the tree */
3034   /* used in dnacomp, dnapars, & dnapenny */
3035   long i, tipy, dummy;
3036   double scale;
3037
3038   putc('\n', outfile);
3039   if (!treeprint)
3040     return;
3041   putc('\n', outfile);
3042   tipy = 1;
3043   dummy = 0;
3044   coordinates(root, &tipy, f, &dummy);
3045   scale = 1.5;
3046   putc('\n', outfile);
3047   for (i = 1; i <= (tipy - down); i++)
3048     drawline(i, scale, root);
3049   fprintf(outfile, "\n  remember:");
3050   if (outgropt)
3051     fprintf(outfile, " (although rooted by outgroup)");
3052   fprintf(outfile, " this is an unrooted tree!\n\n");
3053 }  /* printree */
3054
3055
3056 void writesteps(long chars, boolean weights, steptr oldweight, node *root)
3057 {
3058   /* used in dnacomp, dnapars, & dnapenny */
3059   long i, j, k, l;
3060
3061   putc('\n', outfile);
3062   if (weights)
3063     fprintf(outfile, "weighted ");
3064   fprintf(outfile, "steps in each site:\n");
3065   fprintf(outfile, "      ");
3066   for (i = 0; i <= 9; i++)
3067     fprintf(outfile, "%4ld", i);
3068   fprintf(outfile, "\n     *------------------------------------");
3069   fprintf(outfile, "-----\n");
3070   for (i = 0; i <= (chars / 10); i++) {
3071     fprintf(outfile, "%5ld", i * 10);
3072     putc('|', outfile);
3073     for (j = 0; j <= 9; j++) {
3074       k = i * 10 + j;
3075       if (k == 0 || k > chars)
3076         fprintf(outfile, "    ");
3077       else {
3078         l = location[ally[k - 1] - 1];
3079         if (oldweight[k - 1] > 0)
3080           fprintf(outfile, "%4ld",
3081                   oldweight[k - 1] *
3082                   (root->numsteps[l - 1] / weight[l - 1]));
3083         else
3084           fprintf(outfile, "   0");
3085       }
3086     }
3087     putc('\n', outfile);
3088   }
3089 } /* writesteps */
3090
3091
3092 void treeout(node *p, long nextree, long *col, node *root)
3093 {
3094   /* write out file with representation of final tree */
3095   /* used in dnacomp, dnamove, dnapars, & dnapenny */
3096   node *q;
3097   long i, n;
3098   Char c;
3099
3100   if (p->tip) {
3101     n = 0;
3102     for (i = 1; i <= nmlngth; i++) {
3103       if (nayme[p->index - 1][i - 1] != ' ')
3104         n = i;
3105     }
3106     for (i = 0; i < n; i++) {
3107       c = nayme[p->index - 1][i];
3108       if (c == ' ')
3109         c = '_';
3110       putc(c, outtree);
3111     }
3112     *col += n;
3113   } else {
3114     putc('(', outtree);
3115     (*col)++;
3116     q = p->next;
3117     while (q != p) {
3118       treeout(q->back, nextree, col, root);
3119       q = q->next;
3120       if (q == p)
3121         break;
3122       putc(',', outtree);
3123       (*col)++;
3124       if (*col > 60) {
3125         putc('\n', outtree);
3126         *col = 0;
3127       }
3128     }
3129     putc(')', outtree);
3130     (*col)++;
3131   }
3132   if (p != root)
3133     return;
3134   if (nextree > 2)
3135     fprintf(outtree, "[%6.4f];\n", 1.0 / (nextree - 1));
3136   else
3137     fprintf(outtree, ";\n");
3138 }  /* treeout */
3139
3140
3141 void treeout3(node *p, long nextree, long *col, node *root)
3142 {
3143   /* write out file with representation of final tree */
3144   /* used in dnapars -- writes branch lengths */
3145   node *q;
3146   long i, n, w;
3147   double x;
3148   Char c;
3149
3150   if (p->tip) {
3151     n = 0;
3152     for (i = 1; i <= nmlngth; i++) {
3153       if (nayme[p->index - 1][i - 1] != ' ')
3154         n = i;
3155     }
3156     for (i = 0; i < n; i++) {
3157       c = nayme[p->index - 1][i];
3158       if (c == ' ')
3159         c = '_';
3160       putc(c, outtree);
3161     }
3162     *col += n;
3163   } else {
3164     putc('(', outtree);
3165     (*col)++;
3166     q = p->next;
3167     while (q != p) {
3168       treeout3(q->back, nextree, col, root);
3169       q = q->next;
3170       if (q == p)
3171         break;
3172       putc(',', outtree);
3173       (*col)++;
3174       if (*col > 60) {
3175         putc('\n', outtree);
3176         *col = 0;
3177       }
3178     }
3179     putc(')', outtree);
3180     (*col)++;
3181   }
3182   x = p->v;
3183   if (x > 0.0)
3184     w = (long)(0.43429448222 * log(x));
3185   else if (x == 0.0)
3186     w = 0;
3187   else
3188     w = (long)(0.43429448222 * log(-x)) + 1;
3189   if (w < 0)
3190     w = 0;
3191   if (p != root) {
3192     fprintf(outtree, ":%*.5f", (int)(w + 7), x);
3193     *col += w + 8;
3194   }
3195   if (p != root)
3196     return;
3197   if (nextree > 2)
3198     fprintf(outtree, "[%6.4f];\n", 1.0 / (nextree - 1));
3199   else
3200     fprintf(outtree, ";\n");
3201 }  /* treeout3 */
3202
3203
3204 void drawline2(long i, double scale, tree curtree)
3205 {
3206   /* draws one row of the tree diagram by moving up tree */
3207   /* used in dnaml, proml, & restml */
3208   node *p, *q;
3209   long n, j;
3210   boolean extra;
3211   node *r, *first =NULL, *last =NULL;
3212   boolean done;
3213
3214   p = curtree.start;
3215   q = curtree.start;
3216   extra = false;
3217   if (i == (long)p->ycoord && p == curtree.start) {
3218     if (p->index - spp >= 10)
3219       fprintf(outfile, " %2ld", p->index - spp);
3220     else
3221       fprintf(outfile, "  %ld", p->index - spp);
3222     extra = true;
3223   } else
3224     fprintf(outfile, "  ");
3225   do {
3226     if (!p->tip) {
3227       r = p->next;
3228       done = false;
3229       do {
3230         if (i >= r->back->ymin && i <= r->back->ymax) {
3231           q = r->back;
3232           done = true;
3233         }
3234         r = r->next;
3235       } while (!(done || (p != curtree.start && r == p) ||
3236                  (p == curtree.start && r == p->next)));
3237       first = p->next->back;
3238       r = p;
3239       while (r->next != p)
3240         r = r->next;
3241       last = r->back;
3242       if (p == curtree.start)
3243         last = p->back;
3244     }
3245     done = (p->tip || p == q);
3246     n = (long)(scale * (q->xcoord - p->xcoord) + 0.5);
3247     if (n < 3 && !q->tip)
3248       n = 3;
3249     if (extra) {
3250       n--;
3251       extra = false;
3252     }
3253     if ((long)q->ycoord == i && !done) {
3254       if ((long)p->ycoord != (long)q->ycoord)
3255         putc('+', outfile);
3256       else
3257         putc('-', outfile);
3258       if (!q->tip) {
3259         for (j = 1; j <= n - 2; j++)
3260           putc('-', outfile);
3261         if (q->index - spp >= 10)
3262           fprintf(outfile, "%2ld", q->index - spp);
3263         else
3264           fprintf(outfile, "-%ld", q->index - spp);
3265         extra = true;
3266       } else {
3267         for (j = 1; j < n; j++)
3268           putc('-', outfile);
3269       }
3270     } else if (!p->tip) {
3271       if ((long)last->ycoord > i && (long)first->ycoord < i &&
3272           (i != (long)p->ycoord || p == curtree.start)) {
3273         putc('|', outfile);
3274         for (j = 1; j < n; j++)
3275           putc(' ', outfile);
3276       } else {
3277         for (j = 1; j <= n; j++)
3278           putc(' ', outfile);
3279       }
3280     } else {
3281       for (j = 1; j <= n; j++)
3282         putc(' ', outfile);
3283     }
3284     if (q != p)
3285       p = q;
3286   } while (!done);
3287   if ((long)p->ycoord == i && p->tip) {
3288     for (j = 0; j < nmlngth; j++)
3289       putc(nayme[p->index-1][j], outfile);
3290   }
3291   putc('\n', outfile);
3292 }  /* drawline2 */
3293
3294
3295 void drawline3(long i, double scale, node *start)
3296 {
3297   /* draws one row of the tree diagram by moving up tree */
3298   /* used in dnapars */
3299   node *p, *q;
3300   long n, j;
3301   boolean extra;
3302   node *r, *first =NULL, *last =NULL;
3303   boolean done;
3304
3305   p = start;
3306   q = start;
3307   extra = false;
3308   if (i == (long)p->ycoord) {
3309     if (p->index - spp >= 10)
3310       fprintf(outfile, " %2ld", p->index - spp);
3311     else
3312       fprintf(outfile, "  %ld", p->index - spp);
3313     extra = true;
3314   } else
3315     fprintf(outfile, "  ");
3316   do {
3317     if (!p->tip) {
3318       r = p->next;
3319       done = false;
3320       do {
3321         if (i >= r->back->ymin && i <= r->back->ymax) {
3322           q = r->back;
3323           done = true;
3324         }
3325         r = r->next;
3326       } while (!(done || (r == p)));
3327       first = p->next->back;
3328       r = p;
3329       while (r->next != p)
3330         r = r->next;
3331       last = r->back;
3332     }
3333     done = (p->tip || p == q);
3334     n = (long)(scale * (q->xcoord - p->xcoord) + 0.5);
3335     if (n < 3 && !q->tip)
3336       n = 3;
3337     if (extra) {
3338       n--;
3339       extra = false;
3340     }
3341     if ((long)q->ycoord == i && !done) {
3342       if ((long)p->ycoord != (long)q->ycoord)
3343         putc('+', outfile);
3344       else
3345         putc('-', outfile);
3346       if (!q->tip) {
3347         for (j = 1; j <= n - 2; j++)
3348           putc('-', outfile);
3349         if (q->index - spp >= 10)
3350           fprintf(outfile, "%2ld", q->index - spp);
3351         else
3352           fprintf(outfile, "-%ld", q->index - spp);
3353         extra = true;
3354       } else {
3355         for (j = 1; j < n; j++)
3356           putc('-', outfile);
3357       }
3358     } else if (!p->tip) {
3359       if ((long)last->ycoord > i && (long)first->ycoord < i &&
3360           (i != (long)p->ycoord || p == start)) {
3361         putc('|', outfile);
3362         for (j = 1; j < n; j++)
3363           putc(' ', outfile);
3364       } else {
3365         for (j = 1; j <= n; j++)
3366           putc(' ', outfile);
3367       }
3368     } else {
3369       for (j = 1; j <= n; j++)
3370         putc(' ', outfile);
3371     }
3372     if (q != p)
3373       p = q;
3374   } while (!done);
3375   if ((long)p->ycoord == i && p->tip) {
3376     for (j = 0; j < nmlngth; j++)
3377       putc(nayme[p->index-1][j], outfile);
3378   }
3379   putc('\n', outfile);
3380 }  /* drawline3 */
3381
3382
3383 void copynode(node *c, node *d, long categs)
3384 {
3385   long i, j;
3386
3387   for (i = 0; i < endsite; i++)
3388     for (j = 0; j < categs; j++)
3389       memcpy(d->x[i][j], c->x[i][j], sizeof(sitelike));
3390   memcpy(d->underflows,c->underflows,sizeof(double) * endsite);
3391   d->tyme = c->tyme;
3392   d->v = c->v;
3393   d->xcoord = c->xcoord;
3394   d->ycoord = c->ycoord;
3395   d->ymin = c->ymin;
3396   d->ymax = c->ymax;
3397   d->iter = c->iter;                   /* iter used in dnaml only */
3398   d->haslength = c->haslength;         /* haslength used in dnamlk only */
3399   d->initialized = c->initialized;     /* initialized used in dnamlk only */
3400 }  /* copynode */
3401
3402
3403 void prot_copynode(node *c, node *d, long categs)
3404 {
3405   /* a version of copynode for proml */
3406   long i, j;
3407
3408   for (i = 0; i < endsite; i++)
3409     for (j = 0; j < categs; j++)
3410       memcpy(d->protx[i][j], c->protx[i][j], sizeof(psitelike));
3411   memcpy(d->underflows,c->underflows,sizeof(double) * endsite);
3412   d->tyme = c->tyme;
3413   d->v = c->v;
3414   d->xcoord = c->xcoord;
3415   d->ycoord = c->ycoord;
3416   d->ymin = c->ymin;
3417   d->ymax = c->ymax;
3418   d->iter = c->iter;                   /* iter used in dnaml only */
3419   d->haslength = c->haslength;         /* haslength used in dnamlk only */
3420   d->initialized = c->initialized;     /* initialized used in dnamlk only */
3421 }  /* prot_copynode */
3422
3423
3424 void copy_(tree *a, tree *b, long nonodes, long categs)
3425 {
3426   /* used in dnamlk */
3427   long i;
3428   node *p, *q, *r, *s, *t;
3429
3430   for (i = 0; i < spp; i++) {
3431     copynode(a->nodep[i], b->nodep[i], categs);
3432     if (a->nodep[i]->back) {
3433       if (a->nodep[i]->back == a->nodep[a->nodep[i]->back->index - 1])
3434         b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1];
3435       else if (a->nodep[i]->back == a->nodep[a->nodep[i]->back->index - 1]->next)
3436         b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1]->next;
3437       else
3438         b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1]->next->next;
3439     }
3440     else b->nodep[i]->back = NULL;
3441   }
3442   for (i = spp; i < nonodes; i++) {
3443     if (a->nodep[i]) {
3444       p = a->nodep[i];
3445       q = b->nodep[i];
3446           r = p;
3447       do {
3448         copynode(p, q, categs);
3449         if (p->back) {
3450           s = a->nodep[p->back->index - 1];
3451           t = b->nodep[p->back->index - 1];
3452           if (s->tip) {
3453             if(p->back == s)
3454               q->back = t;
3455           } else {
3456             do {
3457               if (p->back == s)
3458                 q->back = t;
3459               s = s->next;
3460               t = t->next;
3461             } while (s != a->nodep[p->back->index - 1]);
3462           }
3463         }
3464         else
3465           q->back = NULL;
3466         p = p->next;
3467         q = q->next;
3468       } while (p != r);
3469     }
3470   }
3471   b->likelihood = a->likelihood;
3472   b->start = a->start;               /* start used in dnaml only */
3473   b->root = a->root;                 /* root used in dnamlk only */
3474 }  /* copy_ */
3475
3476
3477 void prot_copy_(tree *a, tree *b, long nonodes, long categs)
3478 {
3479   /* used in promlk */
3480   /* identical to copy_() except for calls to prot_copynode rather */
3481   /* than copynode.                                                */
3482   long i;
3483   node *p, *q, *r, *s, *t;
3484
3485   for (i = 0; i < spp; i++) {
3486     prot_copynode(a->nodep[i], b->nodep[i], categs);
3487     if (a->nodep[i]->back) {
3488       if (a->nodep[i]->back == a->nodep[a->nodep[i]->back->index - 1])
3489         b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1];
3490       else if (a->nodep[i]->back == a->nodep[a->nodep[i]->back->index - 1]->next
3491 )
3492         b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1]->next;
3493       else
3494         b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1]->next->next;
3495     }
3496     else b->nodep[i]->back = NULL;
3497   }
3498   for (i = spp; i < nonodes; i++) {
3499     if (a->nodep[i]) {
3500       p = a->nodep[i];
3501       q = b->nodep[i];
3502           r = p;
3503       do {
3504         prot_copynode(p, q, categs);
3505         if (p->back) {
3506           s = a->nodep[p->back->index - 1];
3507           t = b->nodep[p->back->index - 1];
3508           if (s->tip)
3509             {
3510                 if(p->back == s)
3511                   q->back = t;
3512           } else {
3513             do {
3514               if (p->back == s)
3515                 q->back = t;
3516               s = s->next;
3517               t = t->next;
3518             } while (s != a->nodep[p->back->index - 1]);
3519           }
3520         }
3521         else
3522           q->back = NULL;
3523         p = p->next;
3524         q = q->next;
3525       } while (p != r);
3526     }
3527   }
3528   b->likelihood = a->likelihood;
3529   b->start = a->start;               /* start used in dnaml only */
3530   b->root = a->root;                 /* root used in dnamlk only */
3531 }  /* prot_copy_ */
3532
3533
3534 void standev(long chars, long numtrees, long minwhich, double minsteps,
3535                         double *nsteps, long **fsteps, longer seed)
3536 {  /* do paired sites test (KHT or SH test) on user-defined trees */
3537    /* used in dnapars & protpars */
3538   long i, j, k;
3539   double wt, sumw, sum, sum2, sd;
3540   double temp;
3541   double **covar, *P, *f, *r;
3542
3543 #define SAMPLES 1000
3544   if (numtrees == 2) {
3545     fprintf(outfile, "Kishino-Hasegawa-Templeton test\n\n");
3546     fprintf(outfile, "Tree    Steps   Diff Steps   Its S.D.");
3547     fprintf(outfile, "   Significantly worse?\n\n");
3548     which = 1;
3549     while (which <= numtrees) {
3550       fprintf(outfile, "%3ld%10.1f", which, nsteps[which - 1] / 10);
3551       if (minwhich == which)
3552         fprintf(outfile, "  <------ best\n");
3553       else {
3554         sumw = 0.0;
3555         sum = 0.0;
3556         sum2 = 0.0;
3557         for (i = 0; i < endsite; i++) {
3558           if (weight[i] > 0) {
3559             wt = weight[i] / 10.0;
3560             sumw += wt;
3561             temp = (fsteps[which - 1][i] - fsteps[minwhich - 1][i]) / 10.0;
3562             sum += wt * temp;
3563             sum2 += wt * temp * temp;
3564           }
3565         }
3566         sd = sqrt(sumw / (sumw - 1.0) * (sum2 - sum * sum / sumw));
3567         fprintf(outfile, "%10.1f%12.4f",
3568                 (nsteps[which - 1] - minsteps) / 10, sd);
3569         if ((sum > 0.0) && (sum > 1.95996 * sd))
3570           fprintf(outfile, "           Yes\n");
3571         else
3572           fprintf(outfile, "           No\n");
3573       }
3574       which++;
3575     }
3576     fprintf(outfile, "\n\n");
3577   } else {           /* Shimodaira-Hasegawa test using normal approximation */
3578     if(numtrees > MAXSHIMOTREES){
3579       fprintf(outfile, "Shimodaira-Hasegawa test on first %d of %ld trees\n\n"
3580               , MAXSHIMOTREES, numtrees);
3581       numtrees = MAXSHIMOTREES;
3582     } else {
3583       fprintf(outfile, "Shimodaira-Hasegawa test\n\n");
3584     }
3585     covar = (double **)Malloc(numtrees*sizeof(double *));
3586     sumw = 0.0;
3587     for (i = 0; i < endsite; i++)
3588       sumw += weight[i] / 10.0;
3589     for (i = 0; i < numtrees; i++)
3590       covar[i] = (double *)Malloc(numtrees*sizeof(double));
3591     for (i = 0; i < numtrees; i++) {        /* compute covariances of trees */
3592       sum = nsteps[i]/(10.0*sumw);
3593       for (j = 0; j <=i; j++) {
3594         sum2 = nsteps[j]/(10.0*sumw);
3595         temp = 0.0;
3596         for (k = 0; k < endsite; k++) {
3597           if (weight[k] > 0) {
3598             wt = weight[k]/10.0;
3599             temp = temp + wt*(fsteps[i][k]/10.0-sum)
3600                             *(fsteps[j][k]/10.0-sum2);
3601           }
3602         }
3603         covar[i][j] = temp;
3604         if (i != j)
3605           covar[j][i] = temp;
3606       }
3607     }
3608     for (i = 0; i < numtrees; i++) { /* in-place Cholesky decomposition
3609                                         of trees x trees covariance matrix */
3610       sum = 0.0;
3611       for (j = 0; j <= i-1; j++)
3612         sum = sum + covar[i][j] * covar[i][j];
3613       if (covar[i][i] <= sum)
3614         temp = 0.0;
3615       else
3616         temp = sqrt(covar[i][i] - sum);
3617       covar[i][i] = temp;
3618       for (j = i+1; j < numtrees; j++) {
3619         sum = 0.0;
3620         for (k = 0; k < i; k++)
3621           sum = sum + covar[i][k] * covar[j][k];
3622         if (fabs(temp) < 1.0E-12)
3623           covar[j][i] = 0.0;
3624         else
3625           covar[j][i] = (covar[j][i] - sum)/temp;
3626       }
3627     }
3628     f = (double *)Malloc(numtrees*sizeof(double)); /* resampled sums */
3629     P = (double *)Malloc(numtrees*sizeof(double)); /* vector of P's of trees */
3630     r = (double *)Malloc(numtrees*sizeof(double)); /* store Normal variates */
3631     for (i = 0; i < numtrees; i++)
3632       P[i] = 0.0;
3633     sum2 = nsteps[0]/10.0;               /* sum2 will be smallest # of steps */
3634     for (i = 1; i < numtrees; i++)
3635       if (sum2 > nsteps[i]/10.0)
3636         sum2 = nsteps[i]/10.0;
3637     for (i = 1; i <= SAMPLES; i++) {          /* loop over resampled trees */
3638       for (j = 0; j < numtrees; j++)          /* draw Normal variates */
3639         r[j] = normrand(seed);
3640       for (j = 0; j < numtrees; j++) {        /* compute vectors */
3641         sum = 0.0;
3642         for (k = 0; k <= j; k++)
3643           sum += covar[j][k]*r[k];
3644         f[j] = sum;
3645       }
3646       sum = f[1];
3647       for (j = 1; j < numtrees; j++)          /* get min of vector */
3648         if (f[j] < sum)
3649           sum = f[j];
3650       for (j = 0; j < numtrees; j++)          /* accumulate P's */
3651         if (nsteps[j]/10.0-sum2 <= f[j] - sum)
3652           P[j] += 1.0/SAMPLES;
3653     }
3654     fprintf(outfile, "Tree    Steps   Diff Steps   P value");
3655     fprintf(outfile, "   Significantly worse?\n\n");
3656     for (i = 0; i < numtrees; i++) {
3657       fprintf(outfile, "%3ld%10.1f", i+1, nsteps[i]/10);
3658       if ((minwhich-1) == i)
3659         fprintf(outfile, "  <------ best\n");
3660       else {
3661         fprintf(outfile, " %9.1f  %10.3f", nsteps[i]/10.0-sum2, P[i]);
3662         if (P[i] < 0.05)
3663           fprintf(outfile, "           Yes\n");
3664         else
3665           fprintf(outfile, "           No\n");
3666       }
3667     }
3668   fprintf(outfile, "\n");
3669   free(P);             /* free the variables we Malloc'ed */
3670   free(f);
3671   free(r);
3672   for (i = 0; i < numtrees; i++)
3673     free(covar[i]);
3674   free(covar);
3675   }
3676 }  /* standev */
3677
3678
3679 void standev2(long numtrees, long maxwhich, long a, long b, double maxlogl,
3680               double *l0gl, double **l0gf, steptr aliasweight, longer seed)
3681 {  /* do paired sites test (KHT or SH) for user-defined trees */
3682   /* used in dnaml, dnamlk, proml, promlk, and restml */
3683   double **covar, *P, *f, *r;
3684   long i, j, k;
3685   double wt, sumw, sum, sum2, sd;
3686   double temp;
3687
3688 #define SAMPLES 1000
3689   if (numtrees == 2) {
3690     fprintf(outfile, "Kishino-Hasegawa-Templeton test\n\n");
3691     fprintf(outfile, "Tree    logL    Diff logL    Its S.D.");
3692     fprintf(outfile, "   Significantly worse?\n\n");
3693     which = 1;
3694     while (which <= numtrees) {
3695       fprintf(outfile, "%3ld %9.1f", which, l0gl[which - 1]);
3696       if (maxwhich == which)
3697         fprintf(outfile, "  <------ best\n");
3698       else {
3699         sumw = 0.0;
3700         sum = 0.0;
3701         sum2 = 0.0;
3702         for (i = a; i <= b; i++) {
3703           if (aliasweight[i] > 0) {
3704             wt = aliasweight[i];
3705             sumw += wt;
3706             temp = l0gf[which - 1][i] - l0gf[maxwhich - 1][i];
3707             sum += temp * wt;
3708             sum2 += wt * temp * temp;
3709           }
3710         }
3711         temp = sum / sumw;
3712         sd = sqrt(sumw / (sumw - 1.0) * (sum2 - sum * sum / sumw ));
3713         fprintf(outfile, "%10.1f %11.4f", (l0gl[which - 1])-maxlogl, sd);
3714         if ((sum < 0.0) && ((-sum) > 1.95996 * sd))
3715           fprintf(outfile, "           Yes\n");
3716         else
3717           fprintf(outfile, "           No\n");
3718       }
3719       which++;
3720     }
3721     fprintf(outfile, "\n\n");
3722   } else {           /* Shimodaira-Hasegawa test using normal approximation */
3723     if(numtrees > MAXSHIMOTREES){
3724       fprintf(outfile, "Shimodaira-Hasegawa test on first %d of %ld trees\n\n"
3725               , MAXSHIMOTREES, numtrees);
3726       numtrees = MAXSHIMOTREES;
3727     } else {
3728       fprintf(outfile, "Shimodaira-Hasegawa test\n\n");
3729     }
3730     covar = (double **)Malloc(numtrees*sizeof(double *));
3731     sumw = 0.0;
3732     for (i = a; i <= b; i++)
3733       sumw += aliasweight[i];
3734     for (i = 0; i < numtrees; i++)
3735       covar[i] = (double *)Malloc(numtrees*sizeof(double));
3736     for (i = 0; i < numtrees; i++) {        /* compute covariances of trees */
3737       sum = l0gl[i]/sumw;
3738       for (j = 0; j <=i; j++) {
3739         sum2 = l0gl[j]/sumw;
3740         temp = 0.0;
3741         for (k = a; k <= b ; k++) {
3742           if (aliasweight[k] > 0) {
3743             wt = aliasweight[k];
3744             temp = temp + wt*(l0gf[i][k]-sum)
3745                             *(l0gf[j][k]-sum2);
3746           }
3747         }
3748         covar[i][j] = temp;
3749         if (i != j)
3750           covar[j][i] = temp;
3751       }
3752     }
3753     for (i = 0; i < numtrees; i++) { /* in-place Cholesky decomposition
3754                                         of trees x trees covariance matrix */
3755       sum = 0.0;
3756       for (j = 0; j <= i-1; j++)
3757         sum = sum + covar[i][j] * covar[i][j];
3758       if (covar[i][i] <= sum)
3759         temp = 0.0;
3760       else
3761         temp = sqrt(covar[i][i] - sum);
3762       covar[i][i] = temp;
3763       for (j = i+1; j < numtrees; j++) {
3764         sum = 0.0;
3765         for (k = 0; k < i; k++)
3766           sum = sum + covar[i][k] * covar[j][k];
3767         if (fabs(temp) < 1.0E-12)
3768           covar[j][i] = 0.0;
3769         else
3770           covar[j][i] = (covar[j][i] - sum)/temp;
3771       }
3772     }
3773     f = (double *)Malloc(numtrees*sizeof(double)); /* resampled likelihoods */
3774     P = (double *)Malloc(numtrees*sizeof(double)); /* vector of P's of trees */
3775     r = (double *)Malloc(numtrees*sizeof(double)); /* store Normal variates */
3776     for (i = 0; i < numtrees; i++)
3777       P[i] = 0.0;
3778     for (i = 1; i <= SAMPLES; i++) {          /* loop over resampled trees */
3779       for (j = 0; j < numtrees; j++)          /* draw Normal variates */
3780         r[j] = normrand(seed);
3781       for (j = 0; j < numtrees; j++) {        /* compute vectors */
3782         sum = 0.0;
3783         for (k = 0; k <= j; k++)
3784           sum += covar[j][k]*r[k];
3785         f[j] = sum;
3786       }
3787       sum = f[1];
3788       for (j = 1; j < numtrees; j++)          /* get max of vector */
3789         if (f[j] > sum)
3790           sum = f[j];
3791       for (j = 0; j < numtrees; j++)          /* accumulate P's */
3792         if (maxlogl-l0gl[j] <= sum-f[j])
3793           P[j] += 1.0/SAMPLES;
3794     }
3795     fprintf(outfile, "Tree    logL    Diff logL    P value");
3796     fprintf(outfile, "   Significantly worse?\n\n");
3797     for (i = 0; i < numtrees; i++) {
3798       fprintf(outfile, "%3ld%10.1f", i+1, l0gl[i]);
3799       if ((maxwhich-1) == i)
3800         fprintf(outfile, "  <------ best\n");
3801       else {
3802         fprintf(outfile, " %9.1f  %10.3f", l0gl[i]-maxlogl, P[i]);
3803         if (P[i] < 0.05)
3804           fprintf(outfile, "           Yes\n");
3805         else
3806           fprintf(outfile, "           No\n");
3807       }
3808     }
3809   fprintf(outfile, "\n");
3810   free(P);             /* free the variables we Malloc'ed */
3811   free(f);
3812   free(r);
3813   for (i = 0; i < numtrees; i++)
3814     free(covar[i]);
3815   free(covar);
3816   }
3817 }  /* standev */
3818
3819
3820 void freetip(node *anode)
3821 {
3822   /* used in dnacomp, dnapars, & dnapenny */
3823
3824   free(anode->numsteps);
3825   free(anode->oldnumsteps);
3826   free(anode->base);
3827   free(anode->oldbase);
3828 }  /* freetip */
3829
3830
3831 void freenontip(node *anode)
3832 {
3833   /* used in dnacomp, dnapars, & dnapenny */
3834
3835   free(anode->numsteps);
3836   free(anode->oldnumsteps);
3837   free(anode->base);
3838   free(anode->oldbase);
3839   free(anode->numnuc);
3840 }  /* freenontip */
3841
3842
3843 void freenodes(long nonodes, pointarray treenode)
3844 {
3845   /* used in dnacomp, dnapars, & dnapenny */
3846   long i;
3847   node *p;
3848
3849   for (i = 0; i < spp; i++)
3850     freetip(treenode[i]);
3851   for (i = spp; i < nonodes; i++) {
3852     if (treenode[i] != NULL) {
3853       p = treenode[i]->next;
3854       do {
3855         freenontip(p);
3856         p = p->next;
3857       } while (p != treenode[i]);
3858       freenontip(p);
3859     }
3860   }
3861 }  /* freenodes */
3862
3863
3864 void freenode(node **anode)
3865 {
3866   /* used in dnacomp, dnapars, & dnapenny */
3867
3868   freenontip(*anode);
3869   free(*anode);
3870 }  /* freenode */
3871
3872
3873 void freetree(long nonodes, pointarray treenode)
3874 {
3875   /* used in dnacomp, dnapars, & dnapenny */
3876   long i;
3877   node *p, *q;
3878
3879   for (i = 0; i < spp; i++)
3880     free(treenode[i]);
3881   for (i = spp; i < nonodes; i++) {
3882     if (treenode[i] != NULL) {
3883       p = treenode[i]->next;
3884       do {
3885         q = p->next;
3886         free(p);
3887         p = q;
3888       } while (p != treenode[i]);
3889       free(p);
3890     }
3891   }
3892   free(treenode);
3893 }  /* freetree */
3894
3895
3896 void prot_freex_notip(long nonodes, pointarray treenode)
3897 {
3898   /* used in proml */
3899   long i, j;
3900   node *p;
3901
3902   for (i = spp; i < nonodes; i++) {
3903     p = treenode[i];
3904     if ( p == NULL ) continue;
3905     do {
3906       for (j = 0; j < endsite; j++){
3907         free(p->protx[j]);
3908         p->protx[j] = NULL;
3909       }
3910       free(p->protx);
3911       p->protx = NULL;
3912       p = p->next;
3913     } while (p != treenode[i]);
3914   }
3915 }  /* prot_freex_notip */
3916
3917
3918 void prot_freex(long nonodes, pointarray treenode)
3919 {
3920   /* used in proml */
3921   long i, j;
3922   node *p;
3923
3924   for (i = 0; i < spp; i++) {
3925     for (j = 0; j < endsite; j++)
3926       free(treenode[i]->protx[j]);
3927     free(treenode[i]->protx);
3928   }
3929   for (i = spp; i < nonodes; i++) {
3930     p = treenode[i];
3931     do {
3932       for (j = 0; j < endsite; j++)
3933         free(p->protx[j]);
3934       free(p->protx);
3935       p = p->next;
3936     } while (p != treenode[i]);
3937   }
3938 }  /* prot_freex */
3939
3940
3941 void freex_notip(long nonodes, pointarray treenode)
3942 {
3943   /* used in dnaml & dnamlk */
3944   long i, j;
3945   node *p;
3946
3947   for (i = spp; i < nonodes; i++) {
3948     p = treenode[i];
3949     if ( p == NULL ) continue;
3950     do {
3951       for (j = 0; j < endsite; j++)
3952         free(p->x[j]);
3953       free(p->x);
3954       p = p->next;
3955     } while (p != treenode[i]);
3956   }
3957 }  /* freex_notip */
3958
3959
3960 void freex(long nonodes, pointarray treenode)
3961 {
3962   /* used in dnaml & dnamlk */
3963   long i, j;
3964   node *p;
3965
3966   for (i = 0; i < spp; i++) {
3967     for (j = 0; j < endsite; j++)
3968       free(treenode[i]->x[j]);
3969     free(treenode[i]->x);
3970   }
3971   for (i = spp; i < nonodes; i++) {
3972     if(treenode[i]){
3973       p = treenode[i];
3974       do {
3975         for (j = 0; j < endsite; j++)
3976           free(p->x[j]);
3977         free(p->x);
3978         p = p->next;
3979       } while (p != treenode[i]);
3980     }
3981   }
3982 }  /* freex */
3983
3984
3985 void freex2(long nonodes, pointarray treenode)
3986 {
3987   /* used in restml */
3988   long i, j;
3989   node *p;
3990
3991   for (i = 0; i < spp; i++)
3992     free(treenode[i]->x2);
3993   for (i = spp; i < nonodes; i++) {
3994     p = treenode[i];
3995     for (j = 1; j <= 3; j++) {
3996       free(p->x2);
3997       p = p->next;
3998     }
3999   }
4000 }  /* freex2 */
4001
4002
4003 void freegarbage(gbases **garbage)
4004 {
4005   /* used in dnacomp, dnapars, & dnapenny */
4006   gbases *p;
4007
4008   while (*garbage) {
4009     p = *garbage;
4010     *garbage = (*garbage)->next;
4011     free(p->base);
4012     free(p);
4013   }
4014 }  /*freegarbage */
4015
4016
4017 void freegrbg(node **grbg)
4018 {
4019   /* used in dnacomp, dnapars, & dnapenny */
4020   node *p;
4021
4022   while (*grbg) {
4023     p = *grbg;
4024     *grbg = (*grbg)->next;
4025     freenontip(p);
4026     free(p);
4027   }
4028 } /*freegrbg */
4029
4030
4031 void collapsetree(node *p, node *root, node **grbg, pointarray treenode,
4032                   long *zeros)
4033 {
4034   /*  Recurse through tree searching for zero length brances between */
4035   /*  nodes (not to tips).  If one exists, collapse the nodes together, */
4036   /*  removing the branch. */
4037   node *q, *x1, *y1, *x2, *y2;
4038   long i, j, index, index2, numd;
4039   if (p->tip)
4040     return;
4041   q = p->next;
4042   do {
4043     if (!q->back->tip && q->v == 0.000000) {
4044       /* merge the two nodes. */
4045       x1 = y2 = q->next;
4046       x2 = y1 = q->back->next;
4047       while(x1->next != q)
4048         x1 = x1-> next;
4049       while(y1->next != q->back)
4050         y1 = y1-> next;
4051       x1->next = x2;
4052       y1->next = y2;
4053
4054       index = q->index;
4055       index2 = q->back->index;
4056       numd = treenode[index-1]->numdesc + q->back->numdesc -1;
4057       chucktreenode(grbg, q->back);
4058       chucktreenode(grbg, q);
4059       q = x2;
4060
4061       /* update the indicies around the node circle */
4062       do{
4063         if(q->index != index){
4064           q->index = index;
4065         }
4066         q = q-> next;
4067       }while(x2 != q);
4068       updatenumdesc(treenode[index-1], root, numd);
4069
4070       /* Alter treenode to point to real nodes, and update indicies */
4071       /* acordingly. */
4072        j = 0; i=0;
4073       for(i = (index2-1); i < nonodes-1 && treenode[i+1]; i++){
4074         treenode[i]=treenode[i+1];
4075         treenode[i+1] = NULL;
4076         x1=x2=treenode[i];
4077         do{
4078           x1->index = i+1;
4079           x1 = x1 -> next;
4080         } while(x1 != x2);
4081       }
4082
4083       /* Create a new empty fork in the blank spot of treenode */
4084       x1=NULL;
4085       for(i=1; i <=3 ; i++){
4086         gnutreenode(grbg, &x2, index2, endsite, zeros);
4087         x2->next = x1;
4088         x1 = x2;
4089       }
4090       x2->next->next->next = x2;
4091       treenode[nonodes-1]=x2;
4092       if (q->back)
4093         collapsetree(q->back, root, grbg, treenode, zeros);
4094     } else {
4095       if (q->back)
4096         collapsetree(q->back, root, grbg, treenode, zeros);
4097       q = q->next;
4098     }
4099   } while (q != p);
4100 } /* collapsetree */
4101
4102
4103 void collapsebestrees(node **root, node **grbg, pointarray treenode,
4104                       bestelm *bestrees, long *place, long *zeros,
4105                       long chars, boolean recompute, boolean progress)
4106
4107 {
4108   /* Goes through all best trees, collapsing trees where possible, and  */
4109   /* deleting trees that are not unique.    */
4110   long i,j, k, pos, nextnode, oldnextree;
4111   boolean found;
4112   node *dummy;
4113
4114   oldnextree = nextree;
4115   for(i = 0 ; i < (oldnextree - 1) ; i++){
4116     bestrees[i].collapse = true;
4117   }
4118
4119   if(progress)
4120     printf("Collapsing best trees\n   ");
4121   k = 0;
4122   for(i = 0 ; i < (oldnextree - 1) ; i++){
4123     if(progress){
4124       if(i % (((oldnextree-1) / 72) + 1) == 0)
4125         putchar('.');
4126       fflush(stdout);
4127     }
4128     while(!bestrees[k].collapse)
4129       k++;
4130     /* Reconstruct tree. */
4131     *root = treenode[0];
4132     add(treenode[0], treenode[1], treenode[spp], root, recompute,
4133         treenode, grbg, zeros);
4134     nextnode = spp + 2;
4135     for (j = 3; j <= spp; j++) {
4136       if (bestrees[k].btree[j - 1] > 0)
4137         add(treenode[bestrees[k].btree[j - 1] - 1], treenode[j - 1],
4138             treenode[nextnode++ - 1], root, recompute, treenode, grbg,
4139             zeros);
4140       else
4141           add(treenode[treenode[-bestrees[k].btree[j - 1]-1]->back->index-1],
4142               treenode[j - 1], NULL, root, recompute, treenode, grbg, zeros);
4143     }
4144     reroot(treenode[outgrno - 1], *root);
4145
4146     treelength(*root, chars, treenode);
4147     collapsetree(*root, *root, grbg, treenode, zeros);
4148     savetree(*root, place, treenode, grbg, zeros);
4149     /* move everything down in the bestree list */
4150     for(j = k ; j < (nextree - 2) ; j++){
4151       memcpy(bestrees[j].btree, bestrees[j + 1].btree, spp * sizeof(long));
4152       bestrees[j].gloreange = bestrees[j + 1].gloreange;
4153       bestrees[j + 1].gloreange = false;
4154       bestrees[j].locreange = bestrees[j + 1].locreange;
4155       bestrees[j + 1].locreange = false;
4156       bestrees[j].collapse = bestrees[j + 1].collapse;
4157     }
4158     pos=0;
4159     findtree(&found, &pos, nextree-1, place, bestrees);
4160
4161     /* put the new tree at the end of the list if it wasn't found */
4162     nextree--;
4163     if(!found)
4164       addtree(pos, &nextree, false, place, bestrees);
4165
4166     /* Deconstruct the tree */
4167     for (j = 1; j < spp; j++){
4168       re_move(treenode[j], &dummy, root, recompute, treenode,
4169               grbg, zeros);
4170     }
4171   }
4172   if (progress) {
4173     putchar('\n');
4174 #ifdef WIN32
4175     phyFillScreenColor();
4176 #endif
4177   }
4178 }