source: trunk/GDE/PHYLIP/dnapars.c

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