source: trunk/GDE/PHYLIP/pars.c

Last change on this file was 19480, checked in by westram, 2 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 47.3 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, global_inseed0;
50double threshold;
51boolean jumble, usertree, thresh, weights, thorough, rearrfirst, trout,
52          progress, stepbox, ancseq, mulsets, justwts, firstset, mulf, global_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 **local_grbg, node *UNUSED_q, long UNUSED_len, long nodei,
438                        long *UNUSED_ntips, long *UNUSED_parens, initops whichinit,
439                        pointarray local_treenode, pointarray UNUSED_nodep, Char *str, Char *UNUSED_ch,
440                        FILE *UNUSED_intree)
441{
442    (void)UNUSED_q;
443    (void)UNUSED_len;
444    (void)UNUSED_ntips;
445    (void)UNUSED_parens;
446    (void)UNUSED_nodep;
447    (void)UNUSED_ch;
448    (void)UNUSED_intree;
449
450  /* initializes a node */
451
452  switch (whichinit) {
453  case bottom:
454    gnudisctreenode(local_grbg, p, nodei, endsite, zeros, zeros2);
455    local_treenode[nodei - 1] = *p;
456    break;
457  case nonbottom:
458    gnudisctreenode(local_grbg, p, nodei, endsite, zeros, zeros2);
459    break;
460  case tip:
461    match_names_to_data (str, local_treenode, p, spp);
462    break;
463  default:      /*cases hslength,hsnolength,treewt,unittrwt,iter,*/
464    break;      /*length should never occur                         */
465  }
466} /* initparsnode */
467
468
469void evaluate(node *r)
470{
471  /* determines the number of steps needed for a tree. this is
472     the minimum number of steps needed to evolve sequences on
473     this tree */
474  long i, steps;
475  long term;
476  double sum;
477
478  sum = 0.0;
479  for (i = 0; i < endsite; i++) {
480    steps = r->numsteps[i];
481    if ((long)steps <= threshwt[i])
482      term = steps;
483    else
484      term = threshwt[i];
485    sum += (double)term;
486    if (usertree && which <= maxuser)
487      fsteps[which - 1][i] = term;
488  }
489  if (usertree && which <= maxuser) {
490    nsteps[which - 1] = sum;
491    if (which == 1) {
492      minwhich = 1;
493      minsteps = sum;
494    } else if (sum < minsteps) {
495      minwhich = which;
496      minsteps = sum;
497    }
498  }
499  like = -sum;
500}  /* evaluate */
501
502
503void tryadd(node *p, node *item, node *nufork)
504{
505  /* temporarily adds one fork and one tip to the tree.
506     if the location where they are added yields greater
507     "likelihood" than other locations tested up to that
508     time, then keeps that location as there */
509  long pos;
510  double belowsum, parentsum;
511  boolean found, collapse, changethere, trysave;
512
513  if (!p->tip) {
514    memcpy(temp->discbase, p->discbase, endsite*sizeof(unsigned char));
515    memcpy(temp->numsteps, p->numsteps, endsite*sizeof(long));
516    memcpy(temp->discnumnuc, p->discnumnuc, endsite*sizeof(discnucarray));
517    temp->numdesc = p->numdesc + 1;
518    if (p->back) {
519      multifillin(temp, tempadd, 1);
520      sumnsteps2(tempsum, temp, p->back, 0, endsite, threshwt);
521    } else {
522      multisumnsteps(temp, tempadd, 0, endsite, threshwt);
523      tempsum->sumsteps = temp->sumsteps;
524    }
525    if (tempsum->sumsteps <= -bestyet) {
526      if (p->back)
527        sumnsteps2(tempsum, temp, p->back, endsite+1, endsite, threshwt);
528      else {
529        multisumnsteps(temp, temp1, endsite+1, endsite, threshwt);
530        tempsum->sumsteps = temp->sumsteps;
531      }
532    }
533    p->sumsteps = tempsum->sumsteps;
534  }
535  if (p == root)
536    sumnsteps2(temp, item, p, 0, endsite, threshwt);
537  else {
538    sumnsteps(temp1, item, p, 0, endsite);
539    sumnsteps2(temp, temp1, p->back, 0, endsite, threshwt);
540  }
541  if (temp->sumsteps <= -bestyet) {
542    if (p == root)
543      sumnsteps2(temp, item, p, endsite+1, endsite, threshwt);
544    else {
545      sumnsteps(temp1, item, p, endsite+1, endsite);
546      sumnsteps2(temp, temp1, p->back, endsite+1, endsite, threshwt);
547    }
548  }
549  belowsum = temp->sumsteps;
550  global_multf = false;
551  like = -belowsum;
552  if (!p->tip && belowsum >= p->sumsteps) {
553    global_multf = true;
554    like = -p->sumsteps;
555  }
556  trysave = true;
557  if (!global_multf && p != root) {
558    parentsum = treenode[p->back->index - 1]->sumsteps;
559    if (belowsum >= parentsum)
560      trysave = false;
561  }
562  if (lastrearr) {
563    changethere = true;
564    if (like >= bstlike2 && trysave) {
565      if (like > bstlike2)
566        found = false;
567      else {
568        addnsave(p, item, nufork, &root, &grbg, global_multf, treenode, place,
569                   zeros, zeros2);
570        pos = 0;
571        findtree(&found, &pos, nextree, place, bestrees);
572      }
573      if (!found) {
574        collapse = collapsible(item, p, temp, temp1, temp2, tempsum, temprm,
575                     tmpadd, global_multf, root, zeros, zeros2, treenode);
576        if (!thorough)
577          changethere = !collapse;
578        if (thorough || !collapse || like > bstlike2) {
579          if (like > bstlike2) {
580            addnsave(p, item, nufork, &root, &grbg, global_multf, treenode,
581                       place, zeros, zeros2);
582            bestlike = bstlike2 = like;
583            addbestever(&pos, &nextree, maxtrees, collapse, place, bestrees);
584          } else
585            addtiedtree(pos, &nextree, maxtrees, collapse, place, bestrees);
586        }
587      }
588    }
589    if (like >= bestyet) {
590      if (like > bstlike2)
591        bstlike2 = like;
592      if (changethere && trysave) {
593        bestyet = like;
594        there = p;
595        mulf = global_multf;
596      }
597    }
598  } else if ((like > bestyet) || (like >= bestyet && trysave)) {
599    bestyet = like;
600    there = p;
601    mulf = global_multf;
602  }
603}  /* tryadd */
604
605
606void addpreorder(node *p, node *item, node *nufork)
607{
608  /* traverses a n-ary tree, calling PROCEDURE tryadd
609     at a node before calling tryadd at its descendants */
610  node *q;
611
612  if (p == NULL)
613    return;
614  tryadd(p, item, nufork);
615  if (!p->tip) {
616    q = p->next;
617    while (q != p) {
618      addpreorder(q->back, item, nufork);
619      q = q->next;
620    }
621  }
622}  /* addpreorder */
623
624
625void trydescendants(node *item, node *forknode, node *parent,
626                        node *parentback, boolean trybelow)
627{
628  /* tries rearrangements at parent and below parent's descendants */
629  node *q, *tempblw;
630  boolean bestever=0, belowbetter, multf=0, saved, trysave;
631  double parentsum=0, belowsum;
632
633  memcpy(temp->discbase, parent->discbase, endsite*sizeof(unsigned char));
634  memcpy(temp->numsteps, parent->numsteps, endsite*sizeof(long));
635  memcpy(temp->discnumnuc, parent->discnumnuc, endsite*sizeof(discnucarray));
636  temp->numdesc = parent->numdesc + 1;
637  multifillin(temp, tempadd, 1);
638  sumnsteps2(tempsum, parentback, temp, 0, endsite, threshwt);
639  belowbetter = true;
640  if (lastrearr) {
641    parentsum = tempsum->sumsteps;
642    if (-tempsum->sumsteps >= bstlike2) {
643      belowbetter = false;
644      bestever = false;
645      multf = true;
646      if (-tempsum->sumsteps > bstlike2)
647        bestever = true;
648      savelocrearr(item, forknode, parent, tmp, tmp1, tmp2, tmp3, tmprm,
649                    tmpadd, &root, maxtrees, &nextree, multf, bestever,
650                    &saved, place, bestrees, treenode, &grbg, zeros, zeros2);
651      if (saved) {
652        like = bstlike2 = -tempsum->sumsteps;
653        there = parent;
654        mulf = true;
655      }
656    }
657  } else if (-tempsum->sumsteps >= like) {
658    there = parent;
659    mulf = true;
660    like = -tempsum->sumsteps;
661  }
662  if (trybelow) {
663    sumnsteps(temp, parent, tempadd, 0, endsite);
664    sumnsteps2(tempsum, temp, parentback, 0, endsite, threshwt);
665    if (lastrearr) {
666      belowsum = tempsum->sumsteps;
667      if (-tempsum->sumsteps >= bstlike2 && belowbetter && 
668            (forknode->numdesc > 2 ||
669               (forknode->numdesc == 2 && 
670                 parent->back->index != forknode->index))) {
671        trysave = false;
672        memcpy(temp->discbase, parentback->discbase, endsite*sizeof(unsigned char));
673        memcpy(temp->numsteps, parentback->numsteps, endsite*sizeof(long));
674        memcpy(temp->discnumnuc, parentback->discnumnuc, endsite*sizeof(discnucarray));
675        temp->numdesc = parentback->numdesc + 1;
676        multifillin(temp, tempadd, 1);
677        sumnsteps2(tempsum, parent, temp, 0, endsite, threshwt);
678        if (-tempsum->sumsteps < bstlike2) {
679          multf = false;
680          bestever = false;
681          trysave = true;
682        }
683        if (-belowsum > bstlike2) {
684          multf = false;
685          bestever = true;
686          trysave = true;
687        }
688        if (trysave) {
689          if (treenode[parent->index - 1] != parent)
690            tempblw = parent->back;
691          else
692            tempblw = parent;
693          savelocrearr(item, forknode, tempblw, tmp, tmp1, tmp2, tmp3, tmprm,
694                         tmpadd, &root, maxtrees, &nextree, multf, bestever,
695                         &saved, place, bestrees, treenode, &grbg,
696                         zeros, zeros2);
697          if (saved) {
698            like = bstlike2 = -belowsum;
699            there = tempblw;
700            mulf = false;
701          }
702        }
703      }
704    } else if (-tempsum->sumsteps > like) {
705      like = -tempsum->sumsteps;
706      if (-tempsum->sumsteps > bestyet) {
707        if (treenode[parent->index - 1] != parent)
708          tempblw = parent->back;
709        else
710          tempblw = parent;
711        there = tempblw;
712        mulf = false;
713      }
714    }
715  }
716  q = parent->next;
717  while (q != parent) {
718    if (q->back && q->back != item) {
719      memcpy(temp1->discbase, q->discbase, endsite*sizeof(unsigned char));
720      memcpy(temp1->numsteps, q->numsteps, endsite*sizeof(long));
721      memcpy(temp1->discnumnuc, q->discnumnuc, endsite*sizeof(discnucarray));
722      temp1->numdesc = q->numdesc;
723      multifillin(temp1, parentback, 0);
724      if (lastrearr)
725        belowbetter = (-parentsum < bstlike2);
726      if (!q->back->tip) {
727        memcpy(temp->discbase, q->back->discbase, endsite*sizeof(unsigned char));
728        memcpy(temp->numsteps, q->back->numsteps, endsite*sizeof(long));
729        memcpy(temp->discnumnuc, q->back->discnumnuc, endsite*sizeof(discnucarray));
730        temp->numdesc = q->back->numdesc + 1;
731        multifillin(temp, tempadd, 1);
732        sumnsteps2(tempsum, temp1, temp, 0, endsite, threshwt);
733        if (lastrearr) {
734          if (-tempsum->sumsteps >= bstlike2) {
735            belowbetter = false;
736            bestever = false;
737            multf = true;
738            if (-tempsum->sumsteps > bstlike2)
739              bestever = true;
740            savelocrearr(item, forknode, q->back, tmp, tmp1, tmp2, tmp3, tmprm,
741                          tmpadd, &root, maxtrees, &nextree, multf, bestever,
742                          &saved, place, bestrees, treenode, &grbg,
743                          zeros, zeros2);
744            if (saved) {
745              like = bstlike2 = -tempsum->sumsteps;
746              there = q->back;
747              mulf = true;
748            }
749          }
750        } else if (-tempsum->sumsteps >= like) {
751          like = -tempsum->sumsteps;
752          there = q->back;
753          mulf = true;
754        }
755      }
756      sumnsteps(temp, q->back, tempadd, 0, endsite);
757      sumnsteps2(tempsum, temp, temp1, 0, endsite, threshwt);
758      if (lastrearr) {
759        if (-tempsum->sumsteps >= bstlike2) {
760          trysave = false;
761          multf = false;
762          if (belowbetter) {
763            bestever = false;
764            trysave = true;
765          }
766          if (-tempsum->sumsteps > bstlike2) {
767            bestever = true;
768            trysave = true;
769          }
770          if (trysave) {
771            if (treenode[q->back->index - 1] != q->back)
772              tempblw = q;
773            else
774              tempblw = q->back;
775            savelocrearr(item, forknode, tempblw, tmp, tmp1, tmp2, tmp3, tmprm,
776                        tmpadd, &root, maxtrees, &nextree, multf, bestever,
777                        &saved, place, bestrees, treenode, &grbg,
778                        zeros, zeros2);
779            if (saved) {
780              like = bstlike2 = -tempsum->sumsteps;
781              there = tempblw;
782              mulf = false;
783            }
784          }
785        }
786      } else if (-tempsum->sumsteps > like) {
787        like = -tempsum->sumsteps;
788        if (-tempsum->sumsteps > bestyet) {
789          if (treenode[q->back->index - 1] != q->back)
790            tempblw = q;
791          else
792            tempblw = q->back;
793          there = tempblw;
794          mulf = false;
795        }
796      }
797    }
798    q = q->next;
799  }
800} /* trydescendants */
801
802
803void trylocal(node *item, node *forknode)
804{
805  /* rearranges below forknode, below descendants of forknode when
806     there are more than 2 descendants, then unroots the back of
807     forknode and rearranges on its descendants */
808  node *q;
809  boolean bestever, multf, saved;
810
811  memcpy(temprm->discbase, zeros2, endsite*sizeof(unsigned char));
812  memcpy(temprm->numsteps, zeros, endsite*sizeof(long));
813  memcpy(temprm->olddiscbase, item->discbase, endsite*sizeof(unsigned char));
814  memcpy(temprm->oldnumsteps, item->numsteps, endsite*sizeof(long));
815  memcpy(tempf->discbase, forknode->discbase, endsite*sizeof(unsigned char));
816  memcpy(tempf->numsteps, forknode->numsteps, endsite*sizeof(long));
817  memcpy(tempf->discnumnuc, forknode->discnumnuc, endsite*sizeof(discnucarray));
818  tempf->numdesc = forknode->numdesc - 1;
819  multifillin(tempf, temprm, -1);
820  if (!forknode->back) {
821    sumnsteps2(tempsum, tempf, tempadd, 0, endsite, threshwt);
822    if (lastrearr) {
823      if (-tempsum->sumsteps > bstlike2) {
824        bestever = true;
825        multf = false;
826        savelocrearr(item, forknode, forknode, tmp, tmp1, tmp2, tmp3, tmprm,
827                       tmpadd, &root, maxtrees, &nextree, multf, bestever,
828                       &saved, place, bestrees, treenode, &grbg,
829                       zeros, zeros2);
830        if (saved) {
831          like = bstlike2 = -tempsum->sumsteps;
832          there = forknode;
833          mulf = false;
834        }
835      }
836    } else if (-tempsum->sumsteps > like) {
837      like = -tempsum->sumsteps;
838      if (-tempsum->sumsteps > bestyet) {
839        there = forknode;
840        mulf = false;
841      }
842    }
843  } else {
844    sumnsteps(temp, tempf, tempadd, 0, endsite);
845    sumnsteps2(tempsum, temp, forknode->back, 0, endsite, threshwt);
846    if (lastrearr) {
847      if (-tempsum->sumsteps > bstlike2) {
848        bestever = true;
849        multf = false;
850        savelocrearr(item, forknode, forknode, tmp, tmp1, tmp2, tmp3, tmprm,
851                       tmpadd, &root, maxtrees, &nextree, multf, bestever,
852                       &saved, place, bestrees, treenode, &grbg,
853                       zeros, zeros2);
854        if (saved) {
855          like = bstlike2 = -tempsum->sumsteps;
856          there = forknode;
857          mulf = false;
858        }
859      }
860    } else if (-tempsum->sumsteps > like) {
861      like = -tempsum->sumsteps;
862      if (-tempsum->sumsteps > bestyet) {
863        there = forknode;
864        mulf = false;
865      }
866    }
867    trydescendants(item, forknode, forknode->back, tempf, false);
868  }
869  q = forknode->next;
870  while (q != forknode) {
871    if (q->back != item) {
872      memcpy(temp2->discbase, q->discbase, endsite*sizeof(unsigned char));
873      memcpy(temp2->numsteps, q->numsteps, endsite*sizeof(long));
874      memcpy(temp2->discnumnuc, q->discnumnuc, endsite*sizeof(discnucarray));
875      temp2->numdesc = q->numdesc - 1;
876      multifillin(temp2, temprm, -1);
877      if (!q->back->tip) {
878        trydescendants(item, forknode, q->back, temp2, true);
879      } else {
880        sumnsteps(temp1, q->back, tempadd, 0, endsite);
881        sumnsteps2(tempsum, temp1, temp2, 0, endsite, threshwt);
882        if (lastrearr) {
883          if (-tempsum->sumsteps > bstlike2) {
884            multf = false;
885            bestever = true;
886            savelocrearr(item, forknode, q->back, tmp, tmp1, tmp2, tmp3,
887                         tmprm, tmpadd, &root, maxtrees, &nextree, multf,
888                         bestever, &saved, place, bestrees, treenode,
889                         &grbg, zeros, zeros2);
890            if (saved) {
891              like = bstlike2 = -tempsum->sumsteps;
892              there = q->back;
893              mulf = false;
894            }
895          }
896        } else if (-tempsum->sumsteps > like) {
897          like = -tempsum->sumsteps;
898          if (-tempsum->sumsteps > bestyet) {
899            there = q->back;
900            mulf = false;
901          }
902        }
903      }
904    }
905    q = q->next;
906  }
907} /* trylocal */
908
909
910void trylocal2(node *item, node *forknode, node *other)
911{
912  /* rearranges below forknode, below descendants of forknode when
913     there are more than 2 descendants, then unroots the back of
914     forknode and rearranges on its descendants.  Used if forknode
915     has binary descendants */
916  node *q;
917  boolean bestever=0, multf, saved, belowbetter, trysave;
918
919  memcpy(tempf->discbase, other->discbase, endsite*sizeof(unsigned char));
920  memcpy(tempf->numsteps, other->numsteps, endsite*sizeof(long));
921  memcpy(tempf->olddiscbase, forknode->discbase, endsite*sizeof(unsigned char));
922  memcpy(tempf->oldnumsteps, forknode->numsteps, endsite*sizeof(long));
923  tempf->numdesc = other->numdesc;
924  if (forknode->back)
925    trydescendants(item, forknode, forknode->back, tempf, false);
926  if (!other->tip) {
927    memcpy(temp->discbase, other->discbase, endsite*sizeof(unsigned char));
928    memcpy(temp->numsteps, other->numsteps, endsite*sizeof(long));
929    memcpy(temp->discnumnuc, other->discnumnuc, endsite*sizeof(discnucarray));
930    temp->numdesc = other->numdesc + 1;
931    multifillin(temp, tempadd, 1);
932    if (forknode->back)
933      sumnsteps2(tempsum, forknode->back, temp, 0, endsite, threshwt);
934    else
935      sumnsteps2(tempsum, NULL, temp, 0, endsite, threshwt);
936    belowbetter = true;
937    if (lastrearr) {
938      if (-tempsum->sumsteps >= bstlike2) {
939        belowbetter = false;
940        bestever = false;
941        multf = true;
942        if (-tempsum->sumsteps > bstlike2)
943          bestever = true;
944        savelocrearr(item, forknode, other, tmp, tmp1, tmp2, tmp3, tmprm,
945                       tmpadd, &root, maxtrees, &nextree, multf, bestever,
946                       &saved, place, bestrees, treenode, &grbg,
947                       zeros, zeros2);
948        if (saved) {
949          like = bstlike2 = -tempsum->sumsteps;
950          there = other;
951          mulf = true;
952        }
953      }
954    } else if (-tempsum->sumsteps >= like) {
955      there = other;
956      mulf = true;
957      like = -tempsum->sumsteps;
958    }
959    if (forknode->back) {
960      memcpy(temprm->discbase, forknode->back->discbase, endsite*sizeof(unsigned char));
961      memcpy(temprm->numsteps, forknode->back->numsteps, endsite*sizeof(long));
962    } else {
963      memcpy(temprm->discbase, zeros2, endsite*sizeof(unsigned char));
964      memcpy(temprm->numsteps, zeros, endsite*sizeof(long));
965    }
966    memcpy(temprm->olddiscbase, other->back->discbase, endsite*sizeof(unsigned char));
967    memcpy(temprm->oldnumsteps, other->back->numsteps, endsite*sizeof(long));
968    q = other->next;
969    while (q != other) {
970      memcpy(temp2->discbase, q->discbase, endsite*sizeof(unsigned char));
971      memcpy(temp2->numsteps, q->numsteps, endsite*sizeof(long));
972      memcpy(temp2->discnumnuc, q->discnumnuc, endsite*sizeof(discnucarray));
973      if (forknode->back) {
974        temp2->numdesc = q->numdesc;
975        multifillin(temp2, temprm, 0);
976      } else {
977        temp2->numdesc = q->numdesc - 1;
978        multifillin(temp2, temprm, -1);
979      }
980      if (!q->back->tip)
981        trydescendants(item, forknode, q->back, temp2, true);
982      else {
983        sumnsteps(temp1, q->back, tempadd, 0, endsite);
984        sumnsteps2(tempsum, temp1, temp2, 0, endsite, threshwt);
985        if (lastrearr) {
986          if (-tempsum->sumsteps >= bstlike2) {
987            trysave = false;
988            multf = false;
989            if (belowbetter) {
990              bestever = false;
991              trysave = true;
992            }
993            if (-tempsum->sumsteps > bstlike2) {
994              bestever = true;
995              trysave = true;
996            }
997            if (trysave) {
998              savelocrearr(item, forknode, q->back, tmp, tmp1, tmp2, tmp3,
999                           tmprm, tmpadd, &root, maxtrees, &nextree, multf,
1000                           bestever, &saved, place, bestrees, treenode,
1001                           &grbg, zeros, zeros2);
1002              if (saved) {
1003                like = bstlike2 = -tempsum->sumsteps;
1004                there = q->back;
1005                mulf = false;
1006              }
1007            }
1008          }
1009        } else if (-tempsum->sumsteps > like) {
1010          like = -tempsum->sumsteps;
1011          if (-tempsum->sumsteps > bestyet) {
1012            there = q->back;
1013            mulf = false;
1014          }
1015        }
1016      }
1017      q = q->next;
1018    }
1019  }
1020} /* trylocal2 */
1021
1022
1023void tryrearr(node *p, boolean *success)
1024{
1025  /* evaluates one rearrangement of the tree.
1026     if the new tree has greater "likelihood" than the old
1027     one sets success = TRUE and keeps the new tree.
1028     otherwise, restores the old tree */
1029  node *forknode, *newfork, *other, *oldthere;
1030  double oldlike;
1031  boolean oldmulf;
1032
1033  if (p->back == NULL)
1034    return;
1035  forknode = treenode[p->back->index - 1]; 
1036  if (!forknode->back && forknode->numdesc <= 2 && alltips(forknode, p))
1037    return;
1038  oldlike = bestyet;
1039  like = -10.0 * spp * chars;
1040  memcpy(tempadd->discbase, p->discbase, endsite*sizeof(unsigned char));
1041  memcpy(tempadd->numsteps, p->numsteps, endsite*sizeof(long));
1042  memcpy(tempadd->olddiscbase, zeros2, endsite*sizeof(unsigned char));
1043  memcpy(tempadd->oldnumsteps, zeros, endsite*sizeof(long));
1044  if (forknode->numdesc > 2) {
1045    oldthere = there = forknode;
1046    oldmulf = mulf = true;
1047    trylocal(p, forknode);
1048  } else {
1049    findbelow(&other, p, forknode);
1050    oldthere = there = other;
1051    oldmulf = mulf = false;
1052    trylocal2(p, forknode, other);
1053  }
1054  if ((like <= oldlike) || (there == oldthere && mulf == oldmulf))
1055    return;
1056  recompute = true;
1057  re_move(p, &forknode, &root, recompute, treenode, &grbg, zeros, zeros2);
1058  if (mulf)
1059    add(there, p, NULL, &root, recompute, treenode, &grbg, zeros, zeros2);
1060  else {
1061    if (forknode->numdesc > 0)
1062      getnufork(&newfork, &grbg, treenode, zeros, zeros2);
1063    else
1064      newfork = forknode;
1065    add(there, p, newfork, &root, recompute, treenode, &grbg, zeros, zeros2);
1066  } 
1067  if (like > oldlike) {
1068    *success = true;
1069    bestyet = like;
1070  }
1071}  /* tryrearr */
1072
1073
1074void repreorder(node *p, boolean *success)
1075{
1076  /* traverses a binary tree, calling PROCEDURE tryrearr
1077     at a node before calling tryrearr at its descendants */
1078  node *q, *this;
1079
1080  if (p == NULL)
1081    return;
1082  if (!p->visited) {
1083    tryrearr(p, success);
1084    p->visited = true;
1085  }
1086  if (!p->tip) {
1087    q = p;
1088    while (q->next != p) {
1089      this = q->next->back;
1090      repreorder(q->next->back,success);
1091      if (q->next->back == this)
1092        q = q->next;
1093    }
1094  }
1095}  /* repreorder */
1096
1097
1098void rearrange(node **r)
1099{
1100  /* traverses the tree (preorder), finding any local
1101     rearrangement which decreases the number of steps.
1102     if traversal succeeds in increasing the tree's
1103     "likelihood", PROCEDURE rearrange runs traversal again */
1104  boolean success=true;
1105
1106  while (success) {
1107    success = false;
1108    clearvisited(treenode);
1109    repreorder(*r, &success);
1110  }
1111}  /* rearrange */
1112
1113
1114void describe()
1115{
1116  /* prints ancestors, steps and table of numbers of steps in
1117     each site */
1118
1119  if (treeprint) {
1120    fprintf(outfile, "\nrequires a total of %10.3f\n", like / -10.0);
1121    fprintf(outfile, "\n  between      and       length\n");
1122    fprintf(outfile, "  -------      ---       ------\n");
1123    printbranchlengths(root);
1124  }
1125  if (stepbox)
1126    writesteps(chars, weights, oldweight, root);
1127  if (ancseq) {
1128    hypstates(chars, root, treenode, &garbage);
1129    putc('\n', outfile);
1130  }
1131  putc('\n', outfile);
1132  if (trout) {
1133    col = 0;
1134    treeout3(root, nextree, &col, root);
1135  }
1136}  /* describe */
1137
1138
1139void pars_coordinates(node *p, double lengthsum, long *tipy,
1140                double *tipmax)
1141{
1142  /* establishes coordinates of nodes */
1143  node *q, *first, *last;
1144  double xx;
1145
1146  if (p == NULL)
1147    return;
1148  if (p->tip) {
1149    p->xcoord = (long)(over * lengthsum + 0.5);
1150    p->ycoord = (*tipy);
1151    p->ymin = (*tipy);
1152    p->ymax = (*tipy);
1153    (*tipy) += down;
1154    if (lengthsum > (*tipmax))
1155      (*tipmax) = lengthsum;
1156    return;
1157  }
1158  q = p->next;
1159  do {
1160    xx = q->v;
1161    if (xx > 100.0)
1162      xx = 100.0;
1163    pars_coordinates(q->back, lengthsum + xx, tipy,tipmax);
1164    q = q->next;
1165  } while (p != q);
1166  first = p->next->back;
1167  q = p;
1168  while (q->next != p)
1169    q = q->next;
1170  last = q->back;
1171  p->xcoord = (long)(over * lengthsum + 0.5);
1172  if ((p == root) || count_sibs(p) > 2)
1173    p->ycoord = p->next->next->back->ycoord;
1174  else
1175    p->ycoord = (first->ycoord + last->ycoord) / 2;
1176  p->ymin = first->ymin;
1177  p->ymax = last->ymax;
1178}  /* pars_coordinates */
1179
1180
1181void pars_printree()
1182{
1183  /* prints out diagram of the tree2 */
1184  long tipy;
1185  double scale, tipmax;
1186  long i;
1187 
1188  if (!treeprint)
1189    return;
1190  putc('\n', outfile);
1191  tipy = 1;
1192  tipmax = 0.0;
1193  pars_coordinates(root, 0.0, &tipy, &tipmax);
1194  scale = 1.0 / (long)(tipmax + 1.000);
1195  for (i = 1; i <= (tipy - down); i++)
1196    drawline3(i, scale, root);
1197  putc('\n', outfile);
1198}  /* pars_printree */
1199
1200
1201void globrearrange()
1202{
1203  /* does global rearrangements */
1204  long j;
1205  double gotlike;
1206  boolean frommulti;
1207  node *item, *nufork;
1208
1209  recompute = true;
1210  do {
1211    printf("   ");
1212    gotlike = bestlike;
1213    for (j = 0; j < nonodes; j++) {
1214      bestyet = -10.0 * spp * chars;
1215      if (j < spp)
1216        item = treenode[enterorder[j] -1];
1217      else 
1218        item = treenode[j];
1219      if ((item != root) && 
1220           ((j < spp) || ((j >= spp) && (item->numdesc > 0))) &&
1221           !((item->back->index == root->index) && (root->numdesc == 2)
1222              && alltips(root, item))) {
1223        re_move(item, &nufork, &root, recompute, treenode, &grbg, zeros, zeros2);
1224        frommulti = (nufork->numdesc > 0);
1225        clearcollapse(treenode);
1226        there = root;
1227        memcpy(tempadd->discbase, item->discbase, endsite*sizeof(unsigned char));
1228        memcpy(tempadd->numsteps, item->numsteps, endsite*sizeof(long));
1229        memcpy(tempadd->olddiscbase, zeros2, endsite*sizeof(unsigned char));
1230        memcpy(tempadd->oldnumsteps, zeros, endsite*sizeof(long));
1231        if (frommulti){
1232          oldnufork = nufork;
1233          getnufork(&nufork, &grbg, treenode, zeros, zeros2);
1234        }
1235        addpreorder(root, item, nufork);
1236        if (frommulti)
1237          oldnufork = NULL;
1238        if (!mulf)
1239          add(there, item, nufork, &root, recompute, treenode, &grbg,
1240                zeros, zeros2);
1241        else
1242          add(there, item, NULL, &root, recompute, treenode, &grbg,
1243                zeros, zeros2);
1244      }
1245      if (progress) {
1246        putchar('.');
1247        fflush(stdout);
1248      }
1249    }
1250    if (progress)
1251      putchar('\n');
1252  } while (bestlike > gotlike);
1253} /* globrearrange */
1254
1255void load_tree(long treei)
1256{ /* restores a tree from bestrees */
1257  long j, nextnode;
1258  boolean local_recompute = false;
1259  node *dummy;
1260
1261  for (j = spp - 1; j >= 1; j--)
1262    re_move(treenode[j], &dummy, &root, local_recompute, treenode, &grbg,
1263              zeros, zeros2);
1264  root = treenode[0];
1265  local_recompute = true;
1266  add(treenode[0], treenode[1], treenode[spp], &root, local_recompute,
1267    treenode, &grbg, zeros, zeros2);
1268  nextnode = spp + 2;
1269  for (j = 3; j <= spp; j++) {
1270    if (bestrees[treei].btree[j - 1] > 0)
1271      add(treenode[bestrees[treei].btree[j - 1] - 1], treenode[j - 1],
1272            treenode[nextnode++ - 1], &root, local_recompute, treenode, &grbg,
1273            zeros, zeros2);
1274    else
1275      add(treenode[treenode[-bestrees[treei].btree[j-1]-1]->back->index-1],
1276            treenode[j - 1], NULL, &root, local_recompute, treenode, &grbg,
1277            zeros, zeros2);
1278  }
1279}
1280
1281void grandrearr()
1282{
1283  /* calls either global rearrangement or local rearrangement on best trees */
1284  long treei;
1285  boolean done;
1286
1287  done = false;
1288  do {
1289    treei = findunrearranged(bestrees, nextree, true);
1290    if (treei < 0)
1291      done = true;
1292    else 
1293        bestrees[treei].gloreange = true;
1294   
1295    if (!done) {
1296      load_tree(treei);
1297      globrearrange();
1298      done = rearrfirst;
1299    }
1300  } while (!done);
1301} /* grandrearr */
1302
1303
1304void maketree()
1305{
1306  /* constructs a binary tree from the pointers in treenode.
1307     adds each node at location which yields highest "likelihood"
1308  then rearranges the tree for greatest "likelihood" */
1309  long i, j, numtrees, nextnode;
1310  boolean done, firsttree, goteof, haslengths;
1311  node *item, *nufork, *dummy;
1312  pointarray nodep;
1313
1314  if (!usertree) {
1315    for (i = 1; i <= spp; i++)
1316      enterorder[i - 1] = i;
1317    if (jumble)
1318      randumize(seed, enterorder);
1319    recompute = true;
1320    root = treenode[enterorder[0] - 1];
1321    add(treenode[enterorder[0] - 1], treenode[enterorder[1] - 1],
1322        treenode[spp], &root, recompute, treenode, &grbg, zeros, zeros2);
1323    if (progress) {
1324      printf("Adding species:\n");
1325      writename(0, 2, enterorder);
1326#ifdef WIN32
1327      phyFillScreenColor();
1328#endif
1329    }
1330    lastrearr = false;
1331    oldnufork = NULL;
1332    for (i = 3; i <= spp; i++) {
1333      bestyet = -10.0 * spp * chars;
1334      item = treenode[enterorder[i - 1] - 1];
1335      getnufork(&nufork, &grbg, treenode, zeros, zeros2);
1336      there = root;
1337      memcpy(tempadd->discbase, item->discbase, endsite*sizeof(unsigned char));
1338      memcpy(tempadd->numsteps, item->numsteps, endsite*sizeof(long));
1339      memcpy(tempadd->olddiscbase, zeros2, endsite*sizeof(unsigned char));
1340      memcpy(tempadd->oldnumsteps, zeros, endsite*sizeof(long));
1341      addpreorder(root, item, nufork);
1342      if (!mulf)
1343        add(there, item, nufork, &root, recompute, treenode, &grbg,
1344              zeros, zeros2);
1345      else
1346        add(there, item, NULL, &root, recompute, treenode, &grbg,
1347              zeros, zeros2);
1348      like = bestyet;
1349      rearrange(&root);
1350      if (progress) {
1351        writename(i - 1, 1, enterorder);
1352#ifdef WIN32
1353        phyFillScreenColor();
1354#endif
1355      }
1356      lastrearr = (i == spp);
1357      if (lastrearr) {
1358        bestlike = bestyet;
1359        if (jumb == 1) {
1360          bstlike2 = bestlike;
1361          nextree = 1;
1362          initbestrees(bestrees, maxtrees, true);
1363          initbestrees(bestrees, maxtrees, false);
1364        }
1365        if (progress) {
1366          printf("\nDoing global rearrangements");
1367          if (rearrfirst)
1368            printf(" on the first of the trees tied for best\n");
1369          else
1370            printf(" on all trees tied for best\n");
1371          printf("  !");
1372          for (j = 1; j <= nonodes; j++)
1373            putchar('-');
1374          printf("!\n");
1375#ifdef WIN32
1376          phyFillScreenColor();
1377#endif
1378    }
1379        globrearrange();
1380        rearrange(&root);
1381      }
1382    }
1383    done = false;
1384    while (!done && findunrearranged(bestrees, nextree, true) >= 0) {
1385      grandrearr();
1386      done = rearrfirst;
1387    }
1388    if (progress) {
1389      putchar('\n');
1390#ifdef WIN32
1391      phyFillScreenColor();
1392#endif
1393    }
1394    recompute = false;
1395    for (i = spp - 1; i >= 1; i--)
1396      re_move(treenode[i], &dummy, &root, recompute, treenode, &grbg,
1397                zeros, zeros2);
1398    if (jumb == njumble) {
1399      if (thorough && (nextree > 2))
1400        reducebestrees(bestrees, &nextree);
1401      if (treeprint) {
1402        putc('\n', outfile);
1403        if (nextree == 2)
1404          fprintf(outfile, "One most parsimonious tree found:\n");
1405        else
1406          fprintf(outfile, "%6ld trees in all found\n", nextree - 1);
1407      }
1408      if (nextree > maxtrees + 1) {
1409        if (treeprint)
1410          fprintf(outfile, "here are the first %4ld of them\n", (long)maxtrees);
1411        nextree = maxtrees + 1;
1412      }
1413      if (treeprint)
1414        putc('\n', outfile);
1415      for (i = 0; i <= (nextree - 2); i++) {
1416        root = treenode[0];
1417        add(treenode[0], treenode[1], treenode[spp], &root, recompute,
1418                   treenode, &grbg, zeros, zeros2);
1419        nextnode = spp + 2;
1420        for (j = 3; j <= spp; j++) {
1421          if (bestrees[i].btree[j - 1] > 0)
1422            add(treenode[bestrees[i].btree[j - 1] - 1], treenode[j - 1],
1423              treenode[nextnode++ - 1], &root, recompute, treenode, &grbg,
1424              zeros, zeros2);
1425          else
1426            add(treenode[treenode[-bestrees[i].btree[j - 1]-1]->back->index-1],
1427              treenode[j - 1], NULL, &root, recompute, treenode, &grbg,
1428              zeros, zeros2);
1429        }
1430        reroot(treenode[outgrno - 1], root);
1431        postorder(root);
1432        evaluate(root);
1433        treelength(root, chars, treenode);
1434        pars_printree();
1435        describe();
1436        for (j = 1; j < spp; j++)
1437          re_move(treenode[j], &dummy, &root, recompute, treenode,
1438                    &grbg, zeros, zeros2);
1439      }
1440    }
1441  } else {
1442    openfile(&intree,INTREE,"input tree", "r",progname,intreename);
1443    fscanf(intree, "%ld%*[^\n]", &numtrees);
1444    getc(intree);
1445    if (numtrees > MAXNUMTREES) {
1446      printf(
1447            "\n\nERROR: number of input trees is read incorrectly from %s\n\n",
1448             intreename);
1449      exxit(-1);
1450    }
1451    if (numtrees > 2)
1452      initseed(&inseed, &global_inseed0, seed);
1453    if (treeprint) {
1454      fprintf(outfile, "User-defined tree");
1455      if (numtrees > 1)
1456        putc('s', outfile);
1457      fprintf(outfile, ":\n");
1458    }
1459    fsteps = (long **)Malloc(maxuser*sizeof(long *));
1460    for (j = 1; j <= maxuser; j++)
1461      fsteps[j - 1] = (long *)Malloc(endsite*sizeof(long));
1462    nodep = NULL;
1463    which = 1;
1464    while (which <= numtrees) {
1465      firsttree = true;
1466      nextnode = 0;
1467      haslengths = true;
1468      treeread(intree, &root, treenode, &goteof, &firsttree,
1469                 nodep, &nextnode, &haslengths,
1470                 &grbg, initparsnode);
1471      if (treeprint)
1472        fprintf(outfile, "\n\n");
1473      if (outgropt)
1474        reroot(treenode[outgrno - 1], root);
1475      postorder(root);
1476      evaluate(root);
1477      treelength(root, chars, treenode);
1478      pars_printree();
1479      describe();
1480      if (which < numtrees)
1481        gdispose(root, &grbg, treenode);
1482      which++;
1483    }
1484    FClose(intree);
1485    putc('\n', outfile);
1486    if (numtrees > 1 && chars > 1 )
1487      standev(chars, numtrees, minwhich, minsteps, nsteps, fsteps, seed);
1488    for (j = 1; j <= maxuser; j++)
1489      free(fsteps[j - 1]);
1490    free(fsteps);
1491  }
1492  if (jumb == njumble) {
1493    if (progress) {
1494      printf("\nOutput written to file \"%s\"\n\n", outfilename);
1495      if (trout) {
1496        printf("Tree");
1497        if (numtrees > 1)
1498          printf("s");
1499        printf(" also written onto file \"%s\"\n\n", outtreename);
1500      }
1501    }
1502  }
1503}  /* maketree */
1504
1505
1506void freerest()
1507{
1508  if (!usertree) {
1509    freenode(&temp);
1510    freenode(&temp1);
1511    freenode(&temp2);
1512    freenode(&tempsum);
1513    freenode(&temprm);
1514    freenode(&tempadd);
1515    freenode(&tempf);
1516    freenode(&tmp);
1517    freenode(&tmp1);
1518    freenode(&tmp2);
1519    freenode(&tmp3);
1520    freenode(&tmprm);
1521    freenode(&tmpadd);
1522  }
1523  freegrbg(&grbg);
1524  if (ancseq)
1525    freegarbage(&garbage);
1526  free(threshwt);
1527  free(zeros);
1528  free(zeros2);
1529  freenodes(nonodes, treenode);
1530}  /* freerest*/
1531
1532
1533int main(int argc, Char *argv[])
1534{  /* Discrete character parsimony by uphill search */
1535
1536  /* reads in spp, chars, and the data. Then calls maketree to
1537     construct the tree */
1538#ifdef MAC
1539   argc = 1;                /* macsetup("Pars","");                */
1540   argv[0]="Pars";
1541#endif
1542  init(argc, argv);
1543  openfile(&infile,INFILE,"input file", "r",argv[0],infilename);
1544  openfile(&outfile,OUTFILE,"output file", "w",argv[0],outfilename);
1545
1546  ibmpc = IBMCRT;
1547  ansi = ANSICRT;
1548  msets = 1;
1549  firstset = true;
1550  garbage = NULL;
1551  grbg = NULL;
1552  doinit();
1553  if (weights || justwts)
1554    openfile(&weightfile,WEIGHTFILE,"weights file","r",argv[0],weightfilename);
1555  if (trout)
1556    openfile(&outtree,OUTTREE,"output tree file", "w",argv[0],outtreename);
1557  for (ith = 1; ith <= msets; ith++) {
1558    if (msets > 1 && !justwts) {
1559      fprintf(outfile, "\nData set # %ld:\n\n", ith);
1560      if (progress)
1561        printf("\nData set # %ld:\n\n", ith);
1562    }
1563    doinput();
1564    if (ith == 1)
1565      firstset = false;
1566    for (jumb = 1; jumb <= njumble; jumb++)
1567      maketree();
1568  freerest();   
1569  }
1570/*  freetree(nonodes, treenode);     debug */
1571  FClose(infile);
1572  FClose(outfile);
1573  if (weights || justwts)
1574    FClose(weightfile);
1575  if (trout)
1576    FClose(outtree);
1577  if (usertree)
1578    FClose(intree);
1579#ifdef MAC
1580  fixmacfile(outfilename);
1581  fixmacfile(outtreename);
1582#endif
1583  if (progress)
1584    printf("Done.\n\n");
1585#ifdef WIN32
1586  phyRestoreConsoleAttributes();
1587#endif
1588  return 0;
1589}  /* Discrete character parsimony by uphill search */
1590
Note: See TracBrowser for help on using the repository browser.