source: branches/stable/GDE/PHYLIP/seq.c

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