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

Last change on this file was 2175, checked in by westram, 20 years ago

upgrade to PHYLIP 3.6
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 46.9 KB
Line 
1
2#include "phylip.h"
3#include "discrete.h"
4
5/* version 3.6 (c) Copyright 1993-2002 by the University of Washington.
6   Written by Joseph Felsenstein, Akiko Fuseki, Sean Lamont, and Andrew Keeffe.
7   Permission is granted to copy and use this program provided no fee is
8   charged for it and provided that this copyright notice is not removed. */
9
10#define MAXNUMTREES     1000000  /* bigger than number of user trees can be */
11
12#ifndef OLDC
13/* function prototypes */
14void   getoptions(void);
15void   allocrest(void);
16void   doinit(void);
17void   makeweights(void);
18void   doinput(void);
19void   initparsnode(node **, node **, node *, long, long, long *, long *,
20                initops, pointarray, pointarray, Char *, Char *, FILE *);
21void   evaluate(node *);
22void   tryadd(node *, node *, node *);
23void   addpreorder(node *, node *, node *);
24void   trydescendants(node *, node *, node *, node *, boolean);
25
26void   trylocal(node *, node *);
27void   trylocal2(node *, node *, node *);
28void   tryrearr(node *p, boolean *);
29void   repreorder(node *p, boolean *);
30void   rearrange(node **);
31void   describe(void);
32void   pars_coordinates(node *, double, long *, double *);
33void   pars_printree(void);
34void   globrearrange(void);
35void   grandrearr(void);
36
37void   maketree(void);
38void   freerest(void);
39void   load_tree(long treei);
40/* function prototypes */
41#endif
42
43
44Char infilename[FNMLNGTH], outfilename[FNMLNGTH], intreename[FNMLNGTH], outtreename[FNMLNGTH],
45     weightfilename[FNMLNGTH];
46node *root;
47long chars, col, msets, ith, njumble, jumb, maxtrees;
48/*   chars = number of sites in actual sequences */
49long inseed, inseed0;
50double threshold;
51boolean jumble, usertree, thresh, weights, thorough, rearrfirst, trout,
52          progress, stepbox, ancseq, mulsets, justwts, firstset, mulf, multf;
53steptr oldweight;
54longer seed;
55pointarray treenode;            /* pointers to all nodes in tree */
56long *enterorder;
57Char progname[20];
58long *zeros;
59unsigned char *zeros2;
60
61/* local variables for Pascal maketree, propagated globally for C version: */
62
63long minwhich;
64double like, minsteps, bestyet, bestlike, bstlike2;
65boolean lastrearr, recompute;
66double nsteps[maxuser];
67long **fsteps;
68node *there, *oldnufork;
69long *place;
70bestelm *bestrees;
71long *threshwt;
72discbaseptr nothing;
73gbases *garbage;
74node *temp, *temp1, *temp2, *tempsum, *temprm, *tempadd, *tempf, *tmp, *tmp1,
75       *tmp2, *tmp3, *tmprm, *tmpadd;
76boolean *names;
77node *grbg;
78
79
80void getoptions()
81{
82  /* interactively set options */
83  long inseed0, loopcount, loopcount2;
84  Char ch, ch2;
85
86  fprintf(outfile, "\nDiscrete character parsimony algorithm, version %s\n\n",
87            VERSION);
88  jumble = false;
89  njumble = 1;
90  outgrno = 1;
91  outgropt = false;
92  thresh = false;
93  thorough = true;
94  rearrfirst = false;
95  maxtrees = 100;
96  trout = true;
97  usertree = false;
98  weights = false;
99  mulsets = false;
100  printdata = false;
101  progress = true;
102  treeprint = true;
103  stepbox = false;
104  ancseq = false;
105  dotdiff = true;
106  interleaved = true;
107  loopcount = 0;
108  for (;;) {
109    cleerhome();
110    printf("\nDiscrete character parsimony algorithm, version %s\n\n",VERSION);
111    printf("Setting for this run:\n");
112    printf("  U                 Search for best tree?  %s\n",
113           (usertree ? "No, use user trees in input file" : "Yes"));
114    if (!usertree) {
115      printf("  S                        Search option?  ");
116      if (thorough)
117        printf("More thorough search\n");
118      else if (rearrfirst)
119        printf("Rearrange on one best tree\n");
120      else
121        printf("Less thorough\n");
122      printf("  V              Number of trees to save?  %ld\n", maxtrees);
123      printf("  J     Randomize input order of species?");
124      if (jumble)
125        printf("  Yes (seed =%8ld,%3ld times)\n", inseed0, njumble);
126      else
127        printf("  No. Use input order\n");
128    }
129    printf("  O                        Outgroup root?");
130    if (outgropt)
131      printf("  Yes, at species number %ld\n", outgrno);
132    else
133      printf("  No, use as outgroup species %ld\n", outgrno);
134    printf("  T              Use Threshold parsimony?");
135    if (thresh)
136      printf("  Yes, count steps up to%4.1f per site\n", threshold);
137    else
138      printf("  No, use ordinary parsimony\n");
139    printf("  W                       Sites weighted?  %s\n",
140           (weights ? "Yes" : "No"));
141    printf("  M           Analyze multiple data sets?");
142    if (mulsets)
143      printf("  Yes, %2ld %s\n", msets,
144               (justwts ? "sets of weights" : "data sets"));
145    else
146      printf("  No\n");
147    printf("  I            Input species interleaved?  %s\n",
148           (interleaved ? "Yes" : "No, sequential"));
149    printf("  0   Terminal type (IBM PC, ANSI, none)?  %s\n",
150           ibmpc ? "IBM PC" : ansi  ? "ANSI"  : "(none)");
151    printf("  1    Print out the data at start of run  %s\n",
152           (printdata ? "Yes" : "No"));
153    printf("  2  Print indications of progress of run  %s\n",
154           progress ? "Yes" : "No");
155    printf("  3                        Print out tree  %s\n",
156           treeprint ? "Yes" : "No");
157    printf("  4          Print out steps in each site  %s\n",
158           stepbox ? "Yes" : "No");
159    printf("  5  Print character at all nodes of tree  %s\n",
160           ancseq ? "Yes" : "No");
161    if (ancseq || printdata)
162      printf("  .  Use dot-differencing to display them  %s\n",
163           dotdiff ? "Yes" : "No");
164    printf("  6       Write out trees onto tree file?  %s\n",
165           trout ? "Yes" : "No");
166    printf("\n  Y to accept these or type the letter for one to change\n");
167#ifdef WIN32
168    phyFillScreenColor();
169#endif
170    scanf("%c%*[^\n]", &ch);
171    getchar();
172    if (ch == '\n')
173      ch = ' ';
174    uppercase(&ch);
175    if (ch == 'Y')
176      break;
177    if (strchr("WSVJOTUMI12345.60",ch) != NULL){
178      switch (ch) {
179       
180      case 'J':
181        jumble = !jumble;
182        if (jumble)
183          initjumble(&inseed, &inseed0, seed, &njumble);
184        else njumble = 1;
185        break;
186       
187      case 'O':
188        outgropt = !outgropt;
189        if (outgropt)
190          initoutgroup(&outgrno, spp);
191        break;
192       
193      case 'T':
194        thresh = !thresh;
195        if (thresh)
196          initthreshold(&threshold);
197        break;
198       
199      case 'W':
200        weights = !weights;
201        break;
202
203      case 'M':
204        mulsets = !mulsets;
205        if (mulsets) {
206          printf("Multiple data sets or multiple weights?");
207          loopcount2 = 0;
208          do {
209            printf(" (type D or W)\n");
210#ifdef WIN32
211            phyFillScreenColor();
212#endif
213            scanf("%c%*[^\n]", &ch2);
214            getchar();
215            if (ch2 == '\n')
216              ch2 = ' ';
217            uppercase(&ch2);
218            countup(&loopcount2, 10);
219          } while ((ch2 != 'W') && (ch2 != 'D'));
220          justwts = (ch2 == 'W');
221          if (justwts)
222            justweights(&msets);
223          else
224            initdatasets(&msets);
225          if (!jumble) {
226            jumble = true;
227            initjumble(&inseed, &inseed0, seed, &njumble);
228          }
229        }
230        break;
231       
232      case 'U':
233        usertree = !usertree;
234        break;
235       
236      case 'S':
237        thorough = !thorough;
238        if (!thorough)
239          printf("Rearrange on just one best tree?");
240          loopcount2 = 0;
241          do {
242            printf(" (type Y or N)\n");
243#ifdef WIN32
244            phyFillScreenColor();
245#endif
246            scanf("%c%*[^\n]", &ch2);
247            getchar();
248            if (ch2 == '\n')
249              ch2 = ' ';
250            uppercase(&ch2);
251            countup(&loopcount2, 10);
252          } while ((ch2 != 'Y') && (ch2 != 'N'));
253          rearrfirst = (ch2 == 'Y');
254        break;
255
256      case 'V':
257        loopcount2 = 0;
258        do {
259          printf("type the number of trees to save\n");
260#ifdef WIN32
261          phyFillScreenColor();
262#endif
263          scanf("%ld%*[^\n]", &maxtrees);
264          if (maxtrees  > MAXNUMTREES)
265            maxtrees = MAXNUMTREES;
266          getchar();
267          countup(&loopcount2, 10);
268        } while (maxtrees < 1);
269        break;
270
271      case 'I':
272        interleaved = !interleaved;
273        break;
274       
275      case '0':
276        initterminal(&ibmpc, &ansi);
277        break;
278       
279      case '1':
280        printdata = !printdata;
281        break;
282       
283      case '2':
284        progress = !progress;
285        break;
286       
287      case '3':
288        treeprint = !treeprint;
289        break;
290       
291      case '4':
292        stepbox = !stepbox;
293        break;
294       
295      case '5':
296        ancseq = !ancseq;
297        break;
298       
299      case '.':
300        dotdiff = !dotdiff;
301        break;
302       
303      case '6':
304        trout = !trout;
305        break;
306      }
307    } else
308      printf("Not a possible option!\n");
309    countup(&loopcount, 100);
310  }
311}  /* getoptions */
312
313
314void allocrest()
315{
316  long i;
317
318  y = (Char **)Malloc(spp*sizeof(Char *));
319  for (i = 0; i < spp; i++)
320    y[i] = (Char *)Malloc(chars*sizeof(Char));
321  convtab = (Char **)Malloc(spp*sizeof(Char *));
322  for (i = 0; i < spp; i++)
323    convtab[i] = (Char *)Malloc(chars*sizeof(Char));
324  bestrees = (bestelm *)Malloc(maxtrees*sizeof(bestelm));
325  for (i = 1; i <= maxtrees; i++)
326    bestrees[i - 1].btree = (long *)Malloc(nonodes*sizeof(long));
327  nayme = (naym *)Malloc(spp*sizeof(naym));
328  enterorder = (long *)Malloc(spp*sizeof(long));
329  place = (long *)Malloc(nonodes*sizeof(long));
330  weight = (long *)Malloc(chars*sizeof(long));
331  oldweight = (long *)Malloc(chars*sizeof(long));
332  alias = (long *)Malloc(chars*sizeof(long));
333  ally = (long *)Malloc(chars*sizeof(long));
334  location = (long *)Malloc(chars*sizeof(long));
335}  /* alocrest */
336
337
338void doinit()
339{
340  /* initializes variables */
341
342  inputnumbers(&spp, &chars, &nonodes, 1);
343  getoptions();
344  if (printdata)
345    fprintf(outfile, "%2ld species, %3ld  sites\n\n", spp, chars);
346  alloctree(&treenode, nonodes, usertree);
347  allocrest();
348}  /* doinit */
349
350
351void makeweights()
352{
353  /* make up weights vector to avoid duplicate computations */
354  long i;
355
356  for (i = 1; i <= chars; i++) {
357    alias[i - 1] = i;
358    oldweight[i - 1] = weight[i - 1];
359    ally[i - 1] = i;
360  }
361  sitesort(chars, weight);
362  sitecombine(chars);
363  sitescrunch(chars);
364  endsite = 0;
365  for (i = 1; i <= chars; i++) {
366    if (ally[i - 1] == i)
367      endsite++;
368  }
369  for (i = 1; i <= endsite; i++)
370    location[alias[i - 1] - 1] = i;
371  if (!thresh)
372    threshold = spp;
373  threshwt = (long *)Malloc(endsite*sizeof(long));
374  for (i = 0; i < endsite; i++) {
375    weight[i] *= 10;
376    threshwt[i] = (long)(threshold * weight[i] + 0.5);
377  }
378  zeros = (long *)Malloc(endsite*sizeof(long));
379  for (i = 0; i < endsite; i++)
380    zeros[i] = 0;
381  zeros2 = (unsigned char *)Malloc(endsite*sizeof(unsigned char));
382  for (i = 0; i < endsite; i++)
383    zeros2[i] = 0;
384}  /* makeweights */
385
386
387void doinput()
388{
389  /* reads the input data */
390  long i;
391
392  if (justwts) {
393    if (firstset)
394      inputdata(chars);
395    for (i = 0; i < chars; i++)
396      weight[i] = 1;
397    inputweights(chars, weight, &weights);
398    if (justwts) {
399      fprintf(outfile, "\n\nWeights set # %ld:\n\n", ith);
400      if (progress)
401        printf("\nWeights set # %ld:\n\n", ith);
402    }
403    if (printdata)
404      printweights(outfile, 0, chars, weight, "Sites");
405  } else {
406    if (!firstset)
407      samenumsp(&chars, ith);
408    inputdata(chars);
409    for (i = 0; i < chars; i++)
410      weight[i] = 1;
411    if (weights) {
412      inputweights(chars, weight, &weights);
413      if (printdata)
414        printweights(outfile, 0, chars, weight, "Sites");
415    }
416  }
417  makeweights();
418  makevalues(treenode, zeros, zeros2, usertree);
419  if (!usertree) {
420    allocdiscnode(&temp, zeros, zeros2, endsite);
421    allocdiscnode(&temp1, zeros, zeros2, endsite);
422    allocdiscnode(&temp2, zeros, zeros2, endsite);
423    allocdiscnode(&tempsum, zeros, zeros2, endsite);
424    allocdiscnode(&temprm, zeros, zeros2, endsite);
425    allocdiscnode(&tempadd, zeros, zeros2, endsite);
426    allocdiscnode(&tempf, zeros, zeros2, endsite);
427    allocdiscnode(&tmp, zeros, zeros2, endsite);
428    allocdiscnode(&tmp1, zeros, zeros2, endsite);
429    allocdiscnode(&tmp2, zeros, zeros2, endsite);
430    allocdiscnode(&tmp3, zeros, zeros2, endsite);
431    allocdiscnode(&tmprm, zeros, zeros2, endsite);
432    allocdiscnode(&tmpadd, zeros, zeros2, endsite);
433  }
434}  /* doinput */
435
436
437void initparsnode(node **p, node **grbg, node *q, long len, long nodei,
438                        long *ntips, long *parens, initops whichinit,
439                        pointarray treenode, pointarray nodep, Char *str, Char *ch,
440                        FILE *intree)
441{
442  /* initializes a node */
443
444  switch (whichinit) {
445  case bottom:
446    gnudisctreenode(grbg, p, nodei, endsite, zeros, zeros2);
447    treenode[nodei - 1] = *p;
448    break;
449  case nonbottom:
450    gnudisctreenode(grbg, p, nodei, endsite, zeros, zeros2);
451    break;
452  case tip:
453    match_names_to_data (str, treenode, p, spp);
454    break;
455  default:      /*cases hslength,hsnolength,treewt,unittrwt,iter,*/
456    break;      /*length should never occur                         */
457  }
458} /* initparsnode */
459
460
461void evaluate(node *r)
462{
463  /* determines the number of steps needed for a tree. this is
464     the minimum number of steps needed to evolve sequences on
465     this tree */
466  long i, steps;
467  long term;
468  double sum;
469
470  sum = 0.0;
471  for (i = 0; i < endsite; i++) {
472    steps = r->numsteps[i];
473    if ((long)steps <= threshwt[i])
474      term = steps;
475    else
476      term = threshwt[i];
477    sum += (double)term;
478    if (usertree && which <= maxuser)
479      fsteps[which - 1][i] = term;
480  }
481  if (usertree && which <= maxuser) {
482    nsteps[which - 1] = sum;
483    if (which == 1) {
484      minwhich = 1;
485      minsteps = sum;
486    } else if (sum < minsteps) {
487      minwhich = which;
488      minsteps = sum;
489    }
490  }
491  like = -sum;
492}  /* evaluate */
493
494
495void tryadd(node *p, node *item, node *nufork)
496{
497  /* temporarily adds one fork and one tip to the tree.
498     if the location where they are added yields greater
499     "likelihood" than other locations tested up to that
500     time, then keeps that location as there */
501  long pos;
502  double belowsum, parentsum;
503  boolean found, collapse, changethere, trysave;
504
505  if (!p->tip) {
506    memcpy(temp->discbase, p->discbase, endsite*sizeof(unsigned char));
507    memcpy(temp->numsteps, p->numsteps, endsite*sizeof(long));
508    memcpy(temp->discnumnuc, p->discnumnuc, endsite*sizeof(discnucarray));
509    temp->numdesc = p->numdesc + 1;
510    if (p->back) {
511      multifillin(temp, tempadd, 1);
512      sumnsteps2(tempsum, temp, p->back, 0, endsite, threshwt);
513    } else {
514      multisumnsteps(temp, tempadd, 0, endsite, threshwt);
515      tempsum->sumsteps = temp->sumsteps;
516    }
517    if (tempsum->sumsteps <= -bestyet) {
518      if (p->back)
519        sumnsteps2(tempsum, temp, p->back, endsite+1, endsite, threshwt);
520      else {
521        multisumnsteps(temp, temp1, endsite+1, endsite, threshwt);
522        tempsum->sumsteps = temp->sumsteps;
523      }
524    }
525    p->sumsteps = tempsum->sumsteps;
526  }
527  if (p == root)
528    sumnsteps2(temp, item, p, 0, endsite, threshwt);
529  else {
530    sumnsteps(temp1, item, p, 0, endsite);
531    sumnsteps2(temp, temp1, p->back, 0, endsite, threshwt);
532  }
533  if (temp->sumsteps <= -bestyet) {
534    if (p == root)
535      sumnsteps2(temp, item, p, endsite+1, endsite, threshwt);
536    else {
537      sumnsteps(temp1, item, p, endsite+1, endsite);
538      sumnsteps2(temp, temp1, p->back, endsite+1, endsite, threshwt);
539    }
540  }
541  belowsum = temp->sumsteps;
542  multf = false;
543  like = -belowsum;
544  if (!p->tip && belowsum >= p->sumsteps) {
545    multf = true;
546    like = -p->sumsteps;
547  }
548  trysave = true;
549  if (!multf && p != root) {
550    parentsum = treenode[p->back->index - 1]->sumsteps;
551    if (belowsum >= parentsum)
552      trysave = false;
553  }
554  if (lastrearr) {
555    changethere = true;
556    if (like >= bstlike2 && trysave) {
557      if (like > bstlike2)
558        found = false;
559      else {
560        addnsave(p, item, nufork, &root, &grbg, multf, treenode, place,
561                   zeros, zeros2);
562        pos = 0;
563        findtree(&found, &pos, nextree, place, bestrees);
564      }
565      if (!found) {
566        collapse = collapsible(item, p, temp, temp1, temp2, tempsum, temprm,
567                     tmpadd, multf, root, zeros, zeros2, treenode);
568        if (!thorough)
569          changethere = !collapse;
570        if (thorough || !collapse || like > bstlike2) {
571          if (like > bstlike2) {
572            addnsave(p, item, nufork, &root, &grbg, multf, treenode, 
573                       place, zeros, zeros2);
574            bestlike = bstlike2 = like;
575            addbestever(&pos, &nextree, maxtrees, collapse, place, bestrees);
576          } else
577            addtiedtree(pos, &nextree, maxtrees, collapse, place, bestrees);
578        }
579      }
580    }
581    if (like >= bestyet) {
582      if (like > bstlike2)
583        bstlike2 = like;
584      if (changethere && trysave) {
585        bestyet = like;
586        there = p;
587        mulf = multf;
588      }
589    }
590  } else if ((like > bestyet) || (like >= bestyet && trysave)) {
591    bestyet = like;
592    there = p;
593    mulf = multf;
594  }
595}  /* tryadd */
596
597
598void addpreorder(node *p, node *item, node *nufork)
599{
600  /* traverses a n-ary tree, calling PROCEDURE tryadd
601     at a node before calling tryadd at its descendants */
602  node *q;
603
604  if (p == NULL)
605    return;
606  tryadd(p, item, nufork);
607  if (!p->tip) {
608    q = p->next;
609    while (q != p) {
610      addpreorder(q->back, item, nufork);
611      q = q->next;
612    }
613  }
614}  /* addpreorder */
615
616
617void trydescendants(node *item, node *forknode, node *parent,
618                        node *parentback, boolean trybelow)
619{
620  /* tries rearrangements at parent and below parent's descendants */
621  node *q, *tempblw;
622  boolean bestever=0, belowbetter, multf=0, saved, trysave;
623  double parentsum=0, belowsum;
624
625  memcpy(temp->discbase, parent->discbase, endsite*sizeof(unsigned char));
626  memcpy(temp->numsteps, parent->numsteps, endsite*sizeof(long));
627  memcpy(temp->discnumnuc, parent->discnumnuc, endsite*sizeof(discnucarray));
628  temp->numdesc = parent->numdesc + 1;
629  multifillin(temp, tempadd, 1);
630  sumnsteps2(tempsum, parentback, temp, 0, endsite, threshwt);
631  belowbetter = true;
632  if (lastrearr) {
633    parentsum = tempsum->sumsteps;
634    if (-tempsum->sumsteps >= bstlike2) {
635      belowbetter = false;
636      bestever = false;
637      multf = true;
638      if (-tempsum->sumsteps > bstlike2)
639        bestever = true;
640      savelocrearr(item, forknode, parent, tmp, tmp1, tmp2, tmp3, tmprm,
641                    tmpadd, &root, maxtrees, &nextree, multf, bestever,
642                    &saved, place, bestrees, treenode, &grbg, zeros, zeros2);
643      if (saved) {
644        like = bstlike2 = -tempsum->sumsteps;
645        there = parent;
646        mulf = true;
647      }
648    }
649  } else if (-tempsum->sumsteps >= like) {
650    there = parent;
651    mulf = true;
652    like = -tempsum->sumsteps;
653  }
654  if (trybelow) {
655    sumnsteps(temp, parent, tempadd, 0, endsite);
656    sumnsteps2(tempsum, temp, parentback, 0, endsite, threshwt);
657    if (lastrearr) {
658      belowsum = tempsum->sumsteps;
659      if (-tempsum->sumsteps >= bstlike2 && belowbetter && 
660            (forknode->numdesc > 2 ||
661               (forknode->numdesc == 2 && 
662                 parent->back->index != forknode->index))) {
663        trysave = false;
664        memcpy(temp->discbase, parentback->discbase, endsite*sizeof(unsigned char));
665        memcpy(temp->numsteps, parentback->numsteps, endsite*sizeof(long));
666        memcpy(temp->discnumnuc, parentback->discnumnuc, endsite*sizeof(discnucarray));
667        temp->numdesc = parentback->numdesc + 1;
668        multifillin(temp, tempadd, 1);
669        sumnsteps2(tempsum, parent, temp, 0, endsite, threshwt);
670        if (-tempsum->sumsteps < bstlike2) {
671          multf = false;
672          bestever = false;
673          trysave = true;
674        }
675        if (-belowsum > bstlike2) {
676          multf = false;
677          bestever = true;
678          trysave = true;
679        }
680        if (trysave) {
681          if (treenode[parent->index - 1] != parent)
682            tempblw = parent->back;
683          else
684            tempblw = parent;
685          savelocrearr(item, forknode, tempblw, tmp, tmp1, tmp2, tmp3, tmprm,
686                         tmpadd, &root, maxtrees, &nextree, multf, bestever,
687                         &saved, place, bestrees, treenode, &grbg,
688                         zeros, zeros2);
689          if (saved) {
690            like = bstlike2 = -belowsum;
691            there = tempblw;
692            mulf = false;
693          }
694        }
695      }
696    } else if (-tempsum->sumsteps > like) {
697      like = -tempsum->sumsteps;
698      if (-tempsum->sumsteps > bestyet) {
699        if (treenode[parent->index - 1] != parent)
700          tempblw = parent->back;
701        else
702          tempblw = parent;
703        there = tempblw;
704        mulf = false;
705      }
706    }
707  }
708  q = parent->next;
709  while (q != parent) {
710    if (q->back && q->back != item) {
711      memcpy(temp1->discbase, q->discbase, endsite*sizeof(unsigned char));
712      memcpy(temp1->numsteps, q->numsteps, endsite*sizeof(long));
713      memcpy(temp1->discnumnuc, q->discnumnuc, endsite*sizeof(discnucarray));
714      temp1->numdesc = q->numdesc;
715      multifillin(temp1, parentback, 0);
716      if (lastrearr)
717        belowbetter = (-parentsum < bstlike2);
718      if (!q->back->tip) {
719        memcpy(temp->discbase, q->back->discbase, endsite*sizeof(unsigned char));
720        memcpy(temp->numsteps, q->back->numsteps, endsite*sizeof(long));
721        memcpy(temp->discnumnuc, q->back->discnumnuc, endsite*sizeof(discnucarray));
722        temp->numdesc = q->back->numdesc + 1;
723        multifillin(temp, tempadd, 1);
724        sumnsteps2(tempsum, temp1, temp, 0, endsite, threshwt);
725        if (lastrearr) {
726          if (-tempsum->sumsteps >= bstlike2) {
727            belowbetter = false;
728            bestever = false;
729            multf = true;
730            if (-tempsum->sumsteps > bstlike2)
731              bestever = true;
732            savelocrearr(item, forknode, q->back, tmp, tmp1, tmp2, tmp3, tmprm,
733                          tmpadd, &root, maxtrees, &nextree, multf, bestever,
734                          &saved, place, bestrees, treenode, &grbg,
735                          zeros, zeros2);
736            if (saved) {
737              like = bstlike2 = -tempsum->sumsteps;
738              there = q->back;
739              mulf = true;
740            }
741          }
742        } else if (-tempsum->sumsteps >= like) {
743          like = -tempsum->sumsteps;
744          there = q->back;
745          mulf = true;
746        }
747      }
748      sumnsteps(temp, q->back, tempadd, 0, endsite);
749      sumnsteps2(tempsum, temp, temp1, 0, endsite, threshwt);
750      if (lastrearr) {
751        if (-tempsum->sumsteps >= bstlike2) {
752          trysave = false;
753          multf = false;
754          if (belowbetter) {
755            bestever = false;
756            trysave = true;
757          }
758          if (-tempsum->sumsteps > bstlike2) {
759            bestever = true;
760            trysave = true;
761          }
762          if (trysave) {
763            if (treenode[q->back->index - 1] != q->back)
764              tempblw = q;
765            else
766              tempblw = q->back;
767            savelocrearr(item, forknode, tempblw, tmp, tmp1, tmp2, tmp3, tmprm,
768                        tmpadd, &root, maxtrees, &nextree, multf, bestever,
769                        &saved, place, bestrees, treenode, &grbg,
770                        zeros, zeros2);
771            if (saved) {
772              like = bstlike2 = -tempsum->sumsteps;
773              there = tempblw;
774              mulf = false;
775            }
776          }
777        }
778      } else if (-tempsum->sumsteps > like) {
779        like = -tempsum->sumsteps;
780        if (-tempsum->sumsteps > bestyet) {
781          if (treenode[q->back->index - 1] != q->back)
782            tempblw = q;
783          else
784            tempblw = q->back;
785          there = tempblw;
786          mulf = false;
787        }
788      }
789    }
790    q = q->next;
791  }
792} /* trydescendants */
793
794
795void trylocal(node *item, node *forknode)
796{
797  /* rearranges below forknode, below descendants of forknode when
798     there are more than 2 descendants, then unroots the back of
799     forknode and rearranges on its descendants */
800  node *q;
801  boolean bestever, multf, saved;
802
803  memcpy(temprm->discbase, zeros2, endsite*sizeof(unsigned char));
804  memcpy(temprm->numsteps, zeros, endsite*sizeof(long));
805  memcpy(temprm->olddiscbase, item->discbase, endsite*sizeof(unsigned char));
806  memcpy(temprm->oldnumsteps, item->numsteps, endsite*sizeof(long));
807  memcpy(tempf->discbase, forknode->discbase, endsite*sizeof(unsigned char));
808  memcpy(tempf->numsteps, forknode->numsteps, endsite*sizeof(long));
809  memcpy(tempf->discnumnuc, forknode->discnumnuc, endsite*sizeof(discnucarray));
810  tempf->numdesc = forknode->numdesc - 1;
811  multifillin(tempf, temprm, -1);
812  if (!forknode->back) {
813    sumnsteps2(tempsum, tempf, tempadd, 0, endsite, threshwt);
814    if (lastrearr) {
815      if (-tempsum->sumsteps > bstlike2) {
816        bestever = true;
817        multf = false;
818        savelocrearr(item, forknode, forknode, tmp, tmp1, tmp2, tmp3, tmprm,
819                       tmpadd, &root, maxtrees, &nextree, multf, bestever,
820                       &saved, place, bestrees, treenode, &grbg,
821                       zeros, zeros2);
822        if (saved) {
823          like = bstlike2 = -tempsum->sumsteps;
824          there = forknode;
825          mulf = false;
826        }
827      }
828    } else if (-tempsum->sumsteps > like) {
829      like = -tempsum->sumsteps;
830      if (-tempsum->sumsteps > bestyet) {
831        there = forknode;
832        mulf = false;
833      }
834    }
835  } else {
836    sumnsteps(temp, tempf, tempadd, 0, endsite);
837    sumnsteps2(tempsum, temp, forknode->back, 0, endsite, threshwt);
838    if (lastrearr) {
839      if (-tempsum->sumsteps > bstlike2) {
840        bestever = true;
841        multf = false;
842        savelocrearr(item, forknode, forknode, tmp, tmp1, tmp2, tmp3, tmprm,
843                       tmpadd, &root, maxtrees, &nextree, multf, bestever,
844                       &saved, place, bestrees, treenode, &grbg,
845                       zeros, zeros2);
846        if (saved) {
847          like = bstlike2 = -tempsum->sumsteps;
848          there = forknode;
849          mulf = false;
850        }
851      }
852    } else if (-tempsum->sumsteps > like) {
853      like = -tempsum->sumsteps;
854      if (-tempsum->sumsteps > bestyet) {
855        there = forknode;
856        mulf = false;
857      }
858    }
859    trydescendants(item, forknode, forknode->back, tempf, false);
860  }
861  q = forknode->next;
862  while (q != forknode) {
863    if (q->back != item) {
864      memcpy(temp2->discbase, q->discbase, endsite*sizeof(unsigned char));
865      memcpy(temp2->numsteps, q->numsteps, endsite*sizeof(long));
866      memcpy(temp2->discnumnuc, q->discnumnuc, endsite*sizeof(discnucarray));
867      temp2->numdesc = q->numdesc - 1;
868      multifillin(temp2, temprm, -1);
869      if (!q->back->tip) {
870        trydescendants(item, forknode, q->back, temp2, true);
871      } else {
872        sumnsteps(temp1, q->back, tempadd, 0, endsite);
873        sumnsteps2(tempsum, temp1, temp2, 0, endsite, threshwt);
874        if (lastrearr) {
875          if (-tempsum->sumsteps > bstlike2) {
876            multf = false;
877            bestever = true;
878            savelocrearr(item, forknode, q->back, tmp, tmp1, tmp2, tmp3,
879                         tmprm, tmpadd, &root, maxtrees, &nextree, multf,
880                         bestever, &saved, place, bestrees, treenode,
881                         &grbg, zeros, zeros2);
882            if (saved) {
883              like = bstlike2 = -tempsum->sumsteps;
884              there = q->back;
885              mulf = false;
886            }
887          }
888        } else if (-tempsum->sumsteps > like) {
889          like = -tempsum->sumsteps;
890          if (-tempsum->sumsteps > bestyet) {
891            there = q->back;
892            mulf = false;
893          }
894        }
895      }
896    }
897    q = q->next;
898  }
899} /* trylocal */
900
901
902void trylocal2(node *item, node *forknode, node *other)
903{
904  /* rearranges below forknode, below descendants of forknode when
905     there are more than 2 descendants, then unroots the back of
906     forknode and rearranges on its descendants.  Used if forknode
907     has binary descendants */
908  node *q;
909  boolean bestever=0, multf, saved, belowbetter, trysave;
910
911  memcpy(tempf->discbase, other->discbase, endsite*sizeof(unsigned char));
912  memcpy(tempf->numsteps, other->numsteps, endsite*sizeof(long));
913  memcpy(tempf->olddiscbase, forknode->discbase, endsite*sizeof(unsigned char));
914  memcpy(tempf->oldnumsteps, forknode->numsteps, endsite*sizeof(long));
915  tempf->numdesc = other->numdesc;
916  if (forknode->back)
917    trydescendants(item, forknode, forknode->back, tempf, false);
918  if (!other->tip) {
919    memcpy(temp->discbase, other->discbase, endsite*sizeof(unsigned char));
920    memcpy(temp->numsteps, other->numsteps, endsite*sizeof(long));
921    memcpy(temp->discnumnuc, other->discnumnuc, endsite*sizeof(discnucarray));
922    temp->numdesc = other->numdesc + 1;
923    multifillin(temp, tempadd, 1);
924    if (forknode->back)
925      sumnsteps2(tempsum, forknode->back, temp, 0, endsite, threshwt);
926    else
927      sumnsteps2(tempsum, NULL, temp, 0, endsite, threshwt);
928    belowbetter = true;
929    if (lastrearr) {
930      if (-tempsum->sumsteps >= bstlike2) {
931        belowbetter = false;
932        bestever = false;
933        multf = true;
934        if (-tempsum->sumsteps > bstlike2)
935          bestever = true;
936        savelocrearr(item, forknode, other, tmp, tmp1, tmp2, tmp3, tmprm,
937                       tmpadd, &root, maxtrees, &nextree, multf, bestever,
938                       &saved, place, bestrees, treenode, &grbg,
939                       zeros, zeros2);
940        if (saved) {
941          like = bstlike2 = -tempsum->sumsteps;
942          there = other;
943          mulf = true;
944        }
945      }
946    } else if (-tempsum->sumsteps >= like) {
947      there = other;
948      mulf = true;
949      like = -tempsum->sumsteps;
950    }
951    if (forknode->back) {
952      memcpy(temprm->discbase, forknode->back->discbase, endsite*sizeof(unsigned char));
953      memcpy(temprm->numsteps, forknode->back->numsteps, endsite*sizeof(long));
954    } else {
955      memcpy(temprm->discbase, zeros2, endsite*sizeof(unsigned char));
956      memcpy(temprm->numsteps, zeros, endsite*sizeof(long));
957    }
958    memcpy(temprm->olddiscbase, other->back->discbase, endsite*sizeof(unsigned char));
959    memcpy(temprm->oldnumsteps, other->back->numsteps, endsite*sizeof(long));
960    q = other->next;
961    while (q != other) {
962      memcpy(temp2->discbase, q->discbase, endsite*sizeof(unsigned char));
963      memcpy(temp2->numsteps, q->numsteps, endsite*sizeof(long));
964      memcpy(temp2->discnumnuc, q->discnumnuc, endsite*sizeof(discnucarray));
965      if (forknode->back) {
966        temp2->numdesc = q->numdesc;
967        multifillin(temp2, temprm, 0);
968      } else {
969        temp2->numdesc = q->numdesc - 1;
970        multifillin(temp2, temprm, -1);
971      }
972      if (!q->back->tip)
973        trydescendants(item, forknode, q->back, temp2, true);
974      else {
975        sumnsteps(temp1, q->back, tempadd, 0, endsite);
976        sumnsteps2(tempsum, temp1, temp2, 0, endsite, threshwt);
977        if (lastrearr) {
978          if (-tempsum->sumsteps >= bstlike2) {
979            trysave = false;
980            multf = false;
981            if (belowbetter) {
982              bestever = false;
983              trysave = true;
984            }
985            if (-tempsum->sumsteps > bstlike2) {
986              bestever = true;
987              trysave = true;
988            }
989            if (trysave) {
990              savelocrearr(item, forknode, q->back, tmp, tmp1, tmp2, tmp3,
991                           tmprm, tmpadd, &root, maxtrees, &nextree, multf,
992                           bestever, &saved, place, bestrees, treenode,
993                           &grbg, zeros, zeros2);
994              if (saved) {
995                like = bstlike2 = -tempsum->sumsteps;
996                there = q->back;
997                mulf = false;
998              }
999            }
1000          }
1001        } else if (-tempsum->sumsteps > like) {
1002          like = -tempsum->sumsteps;
1003          if (-tempsum->sumsteps > bestyet) {
1004            there = q->back;
1005            mulf = false;
1006          }
1007        }
1008      }
1009      q = q->next;
1010    }
1011  }
1012} /* trylocal2 */
1013
1014
1015void tryrearr(node *p, boolean *success)
1016{
1017  /* evaluates one rearrangement of the tree.
1018     if the new tree has greater "likelihood" than the old
1019     one sets success = TRUE and keeps the new tree.
1020     otherwise, restores the old tree */
1021  node *forknode, *newfork, *other, *oldthere;
1022  double oldlike;
1023  boolean oldmulf;
1024
1025  if (p->back == NULL)
1026    return;
1027  forknode = treenode[p->back->index - 1]; 
1028  if (!forknode->back && forknode->numdesc <= 2 && alltips(forknode, p))
1029    return;
1030  oldlike = bestyet;
1031  like = -10.0 * spp * chars;
1032  memcpy(tempadd->discbase, p->discbase, endsite*sizeof(unsigned char));
1033  memcpy(tempadd->numsteps, p->numsteps, endsite*sizeof(long));
1034  memcpy(tempadd->olddiscbase, zeros2, endsite*sizeof(unsigned char));
1035  memcpy(tempadd->oldnumsteps, zeros, endsite*sizeof(long));
1036  if (forknode->numdesc > 2) {
1037    oldthere = there = forknode;
1038    oldmulf = mulf = true;
1039    trylocal(p, forknode);
1040  } else {
1041    findbelow(&other, p, forknode);
1042    oldthere = there = other;
1043    oldmulf = mulf = false;
1044    trylocal2(p, forknode, other);
1045  }
1046  if ((like <= oldlike) || (there == oldthere && mulf == oldmulf))
1047    return;
1048  recompute = true;
1049  re_move(p, &forknode, &root, recompute, treenode, &grbg, zeros, zeros2);
1050  if (mulf)
1051    add(there, p, NULL, &root, recompute, treenode, &grbg, zeros, zeros2);
1052  else {
1053    if (forknode->numdesc > 0)
1054      getnufork(&newfork, &grbg, treenode, zeros, zeros2);
1055    else
1056      newfork = forknode;
1057    add(there, p, newfork, &root, recompute, treenode, &grbg, zeros, zeros2);
1058  } 
1059  if (like > oldlike) {
1060    *success = true;
1061    bestyet = like;
1062  }
1063}  /* tryrearr */
1064
1065
1066void repreorder(node *p, boolean *success)
1067{
1068  /* traverses a binary tree, calling PROCEDURE tryrearr
1069     at a node before calling tryrearr at its descendants */
1070  node *q, *this;
1071
1072  if (p == NULL)
1073    return;
1074  if (!p->visited) {
1075    tryrearr(p, success);
1076    p->visited = true;
1077  }
1078  if (!p->tip) {
1079    q = p;
1080    while (q->next != p) {
1081      this = q->next->back;
1082      repreorder(q->next->back,success);
1083      if (q->next->back == this)
1084        q = q->next;
1085    }
1086  }
1087}  /* repreorder */
1088
1089
1090void rearrange(node **r)
1091{
1092  /* traverses the tree (preorder), finding any local
1093     rearrangement which decreases the number of steps.
1094     if traversal succeeds in increasing the tree's
1095     "likelihood", PROCEDURE rearrange runs traversal again */
1096  boolean success=true;
1097
1098  while (success) {
1099    success = false;
1100    clearvisited(treenode);
1101    repreorder(*r, &success);
1102  }
1103}  /* rearrange */
1104
1105
1106void describe()
1107{
1108  /* prints ancestors, steps and table of numbers of steps in
1109     each site */
1110
1111  if (treeprint) {
1112    fprintf(outfile, "\nrequires a total of %10.3f\n", like / -10.0);
1113    fprintf(outfile, "\n  between      and       length\n");
1114    fprintf(outfile, "  -------      ---       ------\n");
1115    printbranchlengths(root);
1116  }
1117  if (stepbox)
1118    writesteps(chars, weights, oldweight, root);
1119  if (ancseq) {
1120    hypstates(chars, root, treenode, &garbage);
1121    putc('\n', outfile);
1122  }
1123  putc('\n', outfile);
1124  if (trout) {
1125    col = 0;
1126    treeout3(root, nextree, &col, root);
1127  }
1128}  /* describe */
1129
1130
1131void pars_coordinates(node *p, double lengthsum, long *tipy,
1132                double *tipmax)
1133{
1134  /* establishes coordinates of nodes */
1135  node *q, *first, *last;
1136  double xx;
1137
1138  if (p == NULL)
1139    return;
1140  if (p->tip) {
1141    p->xcoord = (long)(over * lengthsum + 0.5);
1142    p->ycoord = (*tipy);
1143    p->ymin = (*tipy);
1144    p->ymax = (*tipy);
1145    (*tipy) += down;
1146    if (lengthsum > (*tipmax))
1147      (*tipmax) = lengthsum;
1148    return;
1149  }
1150  q = p->next;
1151  do {
1152    xx = q->v;
1153    if (xx > 100.0)
1154      xx = 100.0;
1155    pars_coordinates(q->back, lengthsum + xx, tipy,tipmax);
1156    q = q->next;
1157  } while (p != q);
1158  first = p->next->back;
1159  q = p;
1160  while (q->next != p)
1161    q = q->next;
1162  last = q->back;
1163  p->xcoord = (long)(over * lengthsum + 0.5);
1164  if ((p == root) || count_sibs(p) > 2)
1165    p->ycoord = p->next->next->back->ycoord;
1166  else
1167    p->ycoord = (first->ycoord + last->ycoord) / 2;
1168  p->ymin = first->ymin;
1169  p->ymax = last->ymax;
1170}  /* pars_coordinates */
1171
1172
1173void pars_printree()
1174{
1175  /* prints out diagram of the tree2 */
1176  long tipy;
1177  double scale, tipmax;
1178  long i;
1179 
1180  if (!treeprint)
1181    return;
1182  putc('\n', outfile);
1183  tipy = 1;
1184  tipmax = 0.0;
1185  pars_coordinates(root, 0.0, &tipy, &tipmax);
1186  scale = 1.0 / (long)(tipmax + 1.000);
1187  for (i = 1; i <= (tipy - down); i++)
1188    drawline3(i, scale, root);
1189  putc('\n', outfile);
1190}  /* pars_printree */
1191
1192
1193void globrearrange()
1194{
1195  /* does global rearrangements */
1196  long j;
1197  double gotlike;
1198  boolean frommulti;
1199  node *item, *nufork;
1200
1201  recompute = true;
1202  do {
1203    printf("   ");
1204    gotlike = bestlike;
1205    for (j = 0; j < nonodes; j++) {
1206      bestyet = -10.0 * spp * chars;
1207      if (j < spp)
1208        item = treenode[enterorder[j] -1];
1209      else 
1210        item = treenode[j];
1211      if ((item != root) && 
1212           ((j < spp) || ((j >= spp) && (item->numdesc > 0))) &&
1213           !((item->back->index == root->index) && (root->numdesc == 2)
1214              && alltips(root, item))) {
1215        re_move(item, &nufork, &root, recompute, treenode, &grbg, zeros, zeros2);
1216        frommulti = (nufork->numdesc > 0);
1217        clearcollapse(treenode);
1218        there = root;
1219        memcpy(tempadd->discbase, item->discbase, endsite*sizeof(unsigned char));
1220        memcpy(tempadd->numsteps, item->numsteps, endsite*sizeof(long));
1221        memcpy(tempadd->olddiscbase, zeros2, endsite*sizeof(unsigned char));
1222        memcpy(tempadd->oldnumsteps, zeros, endsite*sizeof(long));
1223        if (frommulti){
1224          oldnufork = nufork;
1225          getnufork(&nufork, &grbg, treenode, zeros, zeros2);
1226        }
1227        addpreorder(root, item, nufork);
1228        if (frommulti)
1229          oldnufork = NULL;
1230        if (!mulf)
1231          add(there, item, nufork, &root, recompute, treenode, &grbg,
1232                zeros, zeros2);
1233        else
1234          add(there, item, NULL, &root, recompute, treenode, &grbg,
1235                zeros, zeros2);
1236      }
1237      if (progress) {
1238        putchar('.');
1239        fflush(stdout);
1240      }
1241    }
1242    if (progress)
1243      putchar('\n');
1244  } while (bestlike > gotlike);
1245} /* globrearrange */
1246
1247void load_tree(long treei)
1248{ /* restores a tree from bestrees */
1249  long j, nextnode;
1250  boolean recompute = false;
1251  node *dummy;
1252
1253  for (j = spp - 1; j >= 1; j--)
1254    re_move(treenode[j], &dummy, &root, recompute, treenode, &grbg,
1255              zeros, zeros2);
1256  root = treenode[0];
1257  recompute = true;
1258  add(treenode[0], treenode[1], treenode[spp], &root, recompute,
1259    treenode, &grbg, zeros, zeros2);
1260  nextnode = spp + 2;
1261  for (j = 3; j <= spp; j++) {
1262    if (bestrees[treei].btree[j - 1] > 0)
1263      add(treenode[bestrees[treei].btree[j - 1] - 1], treenode[j - 1],
1264            treenode[nextnode++ - 1], &root, recompute, treenode, &grbg,
1265            zeros, zeros2);
1266    else
1267      add(treenode[treenode[-bestrees[treei].btree[j-1]-1]->back->index-1],
1268            treenode[j - 1], NULL, &root, recompute, treenode, &grbg,
1269            zeros, zeros2);
1270  }
1271}
1272
1273void grandrearr()
1274{
1275  /* calls either global rearrangement or local rearrangement on best trees */
1276  long treei;
1277  boolean done;
1278
1279  done = false;
1280  do {
1281    treei = findunrearranged(bestrees, nextree, true);
1282    if (treei < 0)
1283      done = true;
1284    else 
1285        bestrees[treei].gloreange = true;
1286   
1287    if (!done) {
1288      load_tree(treei);
1289      globrearrange();
1290      done = rearrfirst;
1291    }
1292  } while (!done);
1293} /* grandrearr */
1294
1295
1296void maketree()
1297{
1298  /* constructs a binary tree from the pointers in treenode.
1299     adds each node at location which yields highest "likelihood"
1300  then rearranges the tree for greatest "likelihood" */
1301  long i, j, numtrees, nextnode;
1302  boolean done, firsttree, goteof, haslengths;
1303  node *item, *nufork, *dummy;
1304  pointarray nodep;
1305
1306  if (!usertree) {
1307    for (i = 1; i <= spp; i++)
1308      enterorder[i - 1] = i;
1309    if (jumble)
1310      randumize(seed, enterorder);
1311    recompute = true;
1312    root = treenode[enterorder[0] - 1];
1313    add(treenode[enterorder[0] - 1], treenode[enterorder[1] - 1],
1314        treenode[spp], &root, recompute, treenode, &grbg, zeros, zeros2);
1315    if (progress) {
1316      printf("Adding species:\n");
1317      writename(0, 2, enterorder);
1318#ifdef WIN32
1319      phyFillScreenColor();
1320#endif
1321    }
1322    lastrearr = false;
1323    oldnufork = NULL;
1324    for (i = 3; i <= spp; i++) {
1325      bestyet = -10.0 * spp * chars;
1326      item = treenode[enterorder[i - 1] - 1];
1327      getnufork(&nufork, &grbg, treenode, zeros, zeros2);
1328      there = root;
1329      memcpy(tempadd->discbase, item->discbase, endsite*sizeof(unsigned char));
1330      memcpy(tempadd->numsteps, item->numsteps, endsite*sizeof(long));
1331      memcpy(tempadd->olddiscbase, zeros2, endsite*sizeof(unsigned char));
1332      memcpy(tempadd->oldnumsteps, zeros, endsite*sizeof(long));
1333      addpreorder(root, item, nufork);
1334      if (!mulf)
1335        add(there, item, nufork, &root, recompute, treenode, &grbg,
1336              zeros, zeros2);
1337      else
1338        add(there, item, NULL, &root, recompute, treenode, &grbg,
1339              zeros, zeros2);
1340      like = bestyet;
1341      rearrange(&root);
1342      if (progress) {
1343        writename(i - 1, 1, enterorder);
1344#ifdef WIN32
1345        phyFillScreenColor();
1346#endif
1347      }
1348      lastrearr = (i == spp);
1349      if (lastrearr) {
1350        bestlike = bestyet;
1351        if (jumb == 1) {
1352          bstlike2 = bestlike;
1353          nextree = 1;
1354          initbestrees(bestrees, maxtrees, true);
1355          initbestrees(bestrees, maxtrees, false);
1356        }
1357        if (progress) {
1358          printf("\nDoing global rearrangements");
1359          if (rearrfirst)
1360            printf(" on the first of the trees tied for best\n");
1361          else
1362            printf(" on all trees tied for best\n");
1363          printf("  !");
1364          for (j = 1; j <= nonodes; j++)
1365            putchar('-');
1366          printf("!\n");
1367#ifdef WIN32
1368          phyFillScreenColor();
1369#endif
1370    }
1371        globrearrange();
1372        rearrange(&root);
1373      }
1374    }
1375    done = false;
1376    while (!done && findunrearranged(bestrees, nextree, true) >= 0) {
1377      grandrearr();
1378      done = rearrfirst;
1379    }
1380    if (progress) {
1381      putchar('\n');
1382#ifdef WIN32
1383      phyFillScreenColor();
1384#endif
1385    }
1386    recompute = false;
1387    for (i = spp - 1; i >= 1; i--)
1388      re_move(treenode[i], &dummy, &root, recompute, treenode, &grbg,
1389                zeros, zeros2);
1390    if (jumb == njumble) {
1391      if (thorough && (nextree > 2))
1392        reducebestrees(bestrees, &nextree);
1393      if (treeprint) {
1394        putc('\n', outfile);
1395        if (nextree == 2)
1396          fprintf(outfile, "One most parsimonious tree found:\n");
1397        else
1398          fprintf(outfile, "%6ld trees in all found\n", nextree - 1);
1399      }
1400      if (nextree > maxtrees + 1) {
1401        if (treeprint)
1402          fprintf(outfile, "here are the first %4ld of them\n", (long)maxtrees);
1403        nextree = maxtrees + 1;
1404      }
1405      if (treeprint)
1406        putc('\n', outfile);
1407      for (i = 0; i <= (nextree - 2); i++) {
1408        root = treenode[0];
1409        add(treenode[0], treenode[1], treenode[spp], &root, recompute,
1410                   treenode, &grbg, zeros, zeros2);
1411        nextnode = spp + 2;
1412        for (j = 3; j <= spp; j++) {
1413          if (bestrees[i].btree[j - 1] > 0)
1414            add(treenode[bestrees[i].btree[j - 1] - 1], treenode[j - 1],
1415              treenode[nextnode++ - 1], &root, recompute, treenode, &grbg,
1416              zeros, zeros2);
1417          else
1418            add(treenode[treenode[-bestrees[i].btree[j - 1]-1]->back->index-1],
1419              treenode[j - 1], NULL, &root, recompute, treenode, &grbg,
1420              zeros, zeros2);
1421        }
1422        reroot(treenode[outgrno - 1], root);
1423        postorder(root);
1424        evaluate(root);
1425        treelength(root, chars, treenode);
1426        pars_printree();
1427        describe();
1428        for (j = 1; j < spp; j++)
1429          re_move(treenode[j], &dummy, &root, recompute, treenode,
1430                    &grbg, zeros, zeros2);
1431      }
1432    }
1433  } else {
1434    openfile(&intree,INTREE,"input tree", "r",progname,intreename);
1435    fscanf(intree, "%ld%*[^\n]", &numtrees);
1436    getc(intree);
1437    if (numtrees > MAXNUMTREES) {
1438      printf(
1439            "\n\nERROR: number of input trees is read incorrectly from %s\n\n",
1440             intreename);
1441      exxit(-1);
1442    }
1443    if (numtrees > 2)
1444      initseed(&inseed, &inseed0, seed);
1445    if (treeprint) {
1446      fprintf(outfile, "User-defined tree");
1447      if (numtrees > 1)
1448        putc('s', outfile);
1449      fprintf(outfile, ":\n");
1450    }
1451    fsteps = (long **)Malloc(maxuser*sizeof(long *));
1452    for (j = 1; j <= maxuser; j++)
1453      fsteps[j - 1] = (long *)Malloc(endsite*sizeof(long));
1454    nodep = NULL;
1455    which = 1;
1456    while (which <= numtrees) {
1457      firsttree = true;
1458      nextnode = 0;
1459      haslengths = true;
1460      treeread(intree, &root, treenode, &goteof, &firsttree,
1461                 nodep, &nextnode, &haslengths,
1462                 &grbg, initparsnode);
1463      if (treeprint)
1464        fprintf(outfile, "\n\n");
1465      if (outgropt)
1466        reroot(treenode[outgrno - 1], root);
1467      postorder(root);
1468      evaluate(root);
1469      treelength(root, chars, treenode);
1470      pars_printree();
1471      describe();
1472      if (which < numtrees)
1473        gdispose(root, &grbg, treenode);
1474      which++;
1475    }
1476    FClose(intree);
1477    putc('\n', outfile);
1478    if (numtrees > 1 && chars > 1 )
1479      standev(chars, numtrees, minwhich, minsteps, nsteps, fsteps, seed);
1480    for (j = 1; j <= maxuser; j++)
1481      free(fsteps[j - 1]);
1482    free(fsteps);
1483  }
1484  if (jumb == njumble) {
1485    if (progress) {
1486      printf("\nOutput written to file \"%s\"\n\n", outfilename);
1487      if (trout) {
1488        printf("Tree");
1489        if (numtrees > 1)
1490          printf("s");
1491        printf(" also written onto file \"%s\"\n\n", outtreename);
1492      }
1493    }
1494  }
1495}  /* maketree */
1496
1497
1498void freerest()
1499{
1500  if (!usertree) {
1501    freenode(&temp);
1502    freenode(&temp1);
1503    freenode(&temp2);
1504    freenode(&tempsum);
1505    freenode(&temprm);
1506    freenode(&tempadd);
1507    freenode(&tempf);
1508    freenode(&tmp);
1509    freenode(&tmp1);
1510    freenode(&tmp2);
1511    freenode(&tmp3);
1512    freenode(&tmprm);
1513    freenode(&tmpadd);
1514  }
1515  freegrbg(&grbg);
1516  if (ancseq)
1517    freegarbage(&garbage);
1518  free(threshwt);
1519  free(zeros);
1520  free(zeros2);
1521  freenodes(nonodes, treenode);
1522}  /* freerest*/
1523
1524
1525int main(int argc, Char *argv[])
1526{  /* Discrete character parsimony by uphill search */
1527
1528  /* reads in spp, chars, and the data. Then calls maketree to
1529     construct the tree */
1530#ifdef MAC
1531   argc = 1;                /* macsetup("Pars","");                */
1532   argv[0]="Pars";
1533#endif
1534  init(argc, argv);
1535  openfile(&infile,INFILE,"input file", "r",argv[0],infilename);
1536  openfile(&outfile,OUTFILE,"output file", "w",argv[0],outfilename);
1537
1538  ibmpc = IBMCRT;
1539  ansi = ANSICRT;
1540  msets = 1;
1541  firstset = true;
1542  garbage = NULL;
1543  grbg = NULL;
1544  doinit();
1545  if (weights || justwts)
1546    openfile(&weightfile,WEIGHTFILE,"weights file","r",argv[0],weightfilename);
1547  if (trout)
1548    openfile(&outtree,OUTTREE,"output tree file", "w",argv[0],outtreename);
1549  for (ith = 1; ith <= msets; ith++) {
1550    if (msets > 1 && !justwts) {
1551      fprintf(outfile, "\nData set # %ld:\n\n", ith);
1552      if (progress)
1553        printf("\nData set # %ld:\n\n", ith);
1554    }
1555    doinput();
1556    if (ith == 1)
1557      firstset = false;
1558    for (jumb = 1; jumb <= njumble; jumb++)
1559      maketree();
1560  freerest();   
1561  }
1562/*  freetree(nonodes, treenode);     debug */
1563  FClose(infile);
1564  FClose(outfile);
1565  if (weights || justwts)
1566    FClose(weightfile);
1567  if (trout)
1568    FClose(outtree);
1569  if (usertree)
1570    FClose(intree);
1571#ifdef MAC
1572  fixmacfile(outfilename);
1573  fixmacfile(outtreename);
1574#endif
1575  if (progress)
1576    printf("Done.\n\n");
1577#ifdef WIN32
1578  phyRestoreConsoleAttributes();
1579#endif
1580  return 0;
1581}  /* Discrete character parsimony by uphill search */
1582
Note: See TracBrowser for help on using the repository browser.