source: trunk/GDE/PHYLIP/mix.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: 29.5 KB
Line 
1
2#include "phylip.h"
3#include "disc.h"
4#include "wagner.h"
5
6/* version 3.6. (c) Copyright 1993-2002 by the University of Washington.
7   Written by Joseph Felsenstein, Akiko Fuseki, Sean Lamont, and Andrew Keeffe.
8   Permission is granted to copy and use this program provided no fee is
9   charged for it and provided that this copyright notice is not removed. */
10
11#define maxtrees        100  /* maximum number of tied trees stored     */
12
13typedef long *placeptr;
14
15#ifndef OLDC
16/* function prototypes */
17void   getoptions(void);
18void   allocrest(void);
19void   doinit(void);
20void   inputoptions(void);
21void   doinput(void);
22void   evaluate(node2 *);
23void   reroot(node2 *);
24void   savetraverse(node2 *);
25void   savetree(void);
26void   mix_addtree(long *pos);
27
28void   mix_findtree(boolean *, long *, long, long *, long **);
29void   tryadd(node2 *, node2 **, node2 **);
30void   addpreorder(node2 *, node2 *, node2 *);
31void   tryrearr(node2 *, node2 **, boolean *);
32void   repreorder(node2 *, node2 **, boolean *);
33void   rearrange(node2 **r);
34void   mix_addelement(node2 **, long *, long *, boolean *);
35void   mix_treeread(void);
36void   describe(void);
37void   maketree(void);
38void   reallocchars(void);
39/* function prototypes */
40#endif
41
42
43
44Char infilename[FNMLNGTH], outfilename[FNMLNGTH], intreename[FNMLNGTH], outtreename[FNMLNGTH], weightfilename[FNMLNGTH], ancfilename[FNMLNGTH], mixfilename[FNMLNGTH];
45node2 *root;
46long outgrno, msets, ith, njumble, jumb;
47/*  outgrno indicates outgroup */
48long inseed, inseed0;
49boolean jumble, usertree, weights, thresh, ancvar, questions, allsokal,
50               allwagner, mixture, trout, noroot, outgropt, didreroot,
51               progress, treeprint, stepbox, ancseq, mulsets, firstset,
52               justwts;
53boolean *ancone, *anczero, *ancone0, *anczero0;
54pointptr2 treenode;   /* pointers to all nodes in tree */
55double threshold;
56double *threshwt;
57bitptr wagner, wagner0;
58longer seed;
59long *enterorder;
60double **fsteps;
61char *guess;
62long **bestrees;
63steptr numsteps, numsone, numszero;
64gbit *garbage;
65char global_ch;
66char *progname;
67
68/* Local variables for maketree: */
69long minwhich;
70double like, bestyet, bestlike, bstlike2, minsteps;
71boolean lastrearr,full;
72double nsteps[maxuser];
73node2 *there;
74long fullset;
75bitptr steps, zeroanc, oneanc, fulzeroanc, empzeroanc;
76long *place, col;
77
78void getoptions()
79{
80  /* interactively set options */
81  long loopcount, loopcount2;
82  Char ch, ch2;
83
84  fprintf(outfile, "\nMixed parsimony algorithm, version %s\n\n",VERSION);
85  putchar('\n');
86  jumble = false;
87  njumble = 1;
88  outgrno = 1;
89  outgropt = false;
90  thresh = false;
91  threshold = spp;
92  trout = true;
93  usertree = false;
94  weights = false;
95  justwts = false;
96  ancvar = false;
97  allsokal = false;
98  allwagner = true;
99  mixture = false;
100  printdata = false;
101  progress = true;
102  treeprint = true;
103  stepbox = false;  ancseq = false;
104  loopcount = 0;
105  for (;;) {
106    cleerhome();
107    printf("\nMixed parsimony algorithm, version %s\n\n",VERSION);
108    printf("Settings for this run:\n");
109    printf("  U                 Search for best tree?  %s\n",
110           (usertree ? "No, use user trees in input file" : "Yes"));
111    printf("  X                     Use Mixed method?  %s\n",
112           (mixture ? "Yes" : "No"));
113    printf("  P                     Parsimony method?");
114    if (!mixture) {
115      printf("  %s\n",(allwagner ? "Wagner" : "Camin-Sokal"));
116    } else
117      printf("  (methods in mixture)\n");
118    if (!usertree) {
119      printf("  J     Randomize input order of species?");
120      if (jumble)
121        printf("  Yes (seed =%8ld,%3ld times)\n", inseed0, njumble);
122      else
123        printf("  No. Use input order\n");
124    }
125    printf("  O                        Outgroup root?");
126    if (outgropt)
127      printf("  Yes, at species number%3ld\n", outgrno);
128    else
129      printf("  No, use as outgroup species%3ld\n", outgrno);
130    printf("  T              Use Threshold parsimony?");
131    if (thresh)
132      printf("  Yes, count steps up to%4.1f per char.\n", threshold);
133    else
134      printf("  No, use ordinary parsimony\n");
135    printf("  A   Use ancestral states in input file?  %s\n",
136    (ancvar ? "Yes" : "No"));
137    printf("  W                       Sites weighted?  %s\n",
138           (weights ? "Yes" : "No"));
139    printf("  M           Analyze multiple data sets?");
140    if (mulsets)
141    printf("  Yes, %2ld %s\n", msets,
142               (justwts ? "sets of weights" : "data sets"));
143    else
144      printf("  No\n");
145    printf("  0   Terminal type (IBM PC, ANSI, none)?  %s\n",
146           (ibmpc ? "IBM PC" : ansi  ? "ANSI" : "(none)"));
147    printf("  1    Print out the data at start of run  %s\n",
148           (printdata ? "Yes" : "No"));
149    printf("  2  Print indications of progress of run  %s\n",
150           (progress ? "Yes" : "No"));
151    printf("  3                        Print out tree  %s\n",
152          (treeprint ? "Yes" : "No"));
153    printf("  4     Print out steps in each character  %s\n",
154           (stepbox ? "Yes" : "No"));
155    printf("  5     Print states at all nodes of tree  %s\n",
156           (ancseq ? "Yes" : "No"));
157    printf("  6       Write out trees onto tree file?  %s\n",
158           (trout ? "Yes" : "No"));
159    if(weights && justwts){
160        printf(
161         "WARNING:  W option and Multiple Weights options are both on.  ");
162        printf(
163         "The W menu option is unnecessary and has no additional effect. \n");
164    }
165    printf("\nAre these settings correct? ");
166    printf("(type Y or 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("WJOTUMPAX1234560",ch) != NULL){
178      switch (ch) {
179
180      case 'W':
181        weights = !weights;
182        break;
183       
184      case 'U':
185        usertree = !usertree;
186        break;
187       
188      case 'X':
189        mixture = !mixture;
190        break;
191       
192      case 'P':
193        allwagner = !allwagner;
194        break;
195       
196      case 'A':
197        ancvar = !ancvar;
198        break;
199       
200      case 'J':
201        jumble = !jumble;
202        if (jumble)
203          initjumble(&inseed, &inseed0, seed, &njumble);
204        else njumble = 1;
205        break;
206       
207      case 'O':
208        outgropt = !outgropt;
209        if (outgropt)
210          initoutgroup(&outgrno, spp);
211        break;
212       
213      case 'T':
214        thresh = !thresh;
215        if (thresh)
216          initthreshold(&threshold);
217        break;
218       
219      case 'M':
220        mulsets = !mulsets;
221        if (mulsets){
222            printf("Multiple data sets or multiple weights?");
223          loopcount2 = 0;
224          do {
225            printf(" (type D or W)\n");
226#ifdef WIN32
227            phyFillScreenColor();
228#endif
229            scanf("%c%*[^\n]", &ch2);
230            getchar();
231            if (ch2 == '\n')
232              ch2 = ' ';
233            uppercase(&ch2);
234            countup(&loopcount2, 10);
235          } while ((ch2 != 'W') && (ch2 != 'D'));
236          justwts = (ch2 == 'W');
237          if (justwts)
238            justweights(&msets);
239          else
240            initdatasets(&msets);
241          if (!jumble) {
242            jumble = true;
243            initjumble(&inseed, &inseed0, seed, &njumble);
244          }
245        }
246        break;
247
248      case '0':
249        initterminal(&ibmpc, &ansi);
250        break;
251       
252      case '1':
253        printdata = !printdata;
254        break;
255       
256      case '2':
257        progress = !progress;
258        break;
259       
260      case '3':
261        treeprint = !treeprint;
262        break;
263       
264      case '4':
265        stepbox = !stepbox;
266        break;
267       
268      case '5':
269        ancseq = !ancseq;
270        break;
271       
272      case '6':
273        trout = !trout;
274        break;
275      }
276    } else
277      printf("Not a possible option!\n");
278    countup(&loopcount, 100);
279  }
280  allsokal = (!allwagner && !mixture);
281}  /* getoptions */
282
283void reallocchars() 
284{
285  long i;
286
287  if (usertree) {
288    for (i = 0; i < maxuser; i++) {
289      free (fsteps[i]);
290      fsteps[i] = (double *)Malloc(chars*sizeof(double));
291    }
292  }
293  free(extras);
294  free(weight);
295  free(threshwt);
296  free(numsteps);
297  free(numszero);
298  free(numsone);
299  free(guess);
300  free(ancone);
301  free(anczero);
302  free(ancone0);
303  free(anczero0);
304
305  extras = (steptr)Malloc(chars*sizeof(long));
306  weight = (steptr)Malloc(chars*sizeof(long));
307  threshwt = (double *)Malloc(chars*sizeof(double));
308  numsteps = (steptr)Malloc(chars*sizeof(long));
309  numszero = (steptr)Malloc(chars*sizeof(long));
310  numsone = (steptr)Malloc(chars*sizeof(long));
311  guess = (Char *)Malloc(chars*sizeof(Char));
312  ancone = (boolean *)Malloc(chars*sizeof(boolean));
313  anczero = (boolean *)Malloc(chars*sizeof(boolean));
314  ancone0 = (boolean *)Malloc(chars*sizeof(boolean));
315  anczero0 = (boolean *)Malloc(chars*sizeof(boolean));
316}
317
318void allocrest()
319{
320  long i;
321
322  if (usertree) {
323    fsteps = (double **)Malloc(maxuser*sizeof(double *));
324    for (i = 0; i < maxuser; i++)
325      fsteps[i] = (double *)Malloc(chars*sizeof(double));
326  }
327  bestrees = (long **)Malloc(maxtrees*sizeof(long *));
328  for (i = 1; i <= maxtrees; i++)
329    bestrees[i - 1] = (long *)Malloc(spp*sizeof(long));
330  extras = (steptr)Malloc(chars*sizeof(long));
331  weight = (steptr)Malloc(chars*sizeof(long));
332  threshwt = (double *)Malloc(chars*sizeof(double));
333  numsteps = (steptr)Malloc(chars*sizeof(long));
334  numszero = (steptr)Malloc(chars*sizeof(long));
335  numsone = (steptr)Malloc(chars*sizeof(long));
336  guess = (Char *)Malloc(chars*sizeof(Char));
337  nayme = (naym *)Malloc(spp*sizeof(naym));
338  enterorder = (long *)Malloc(spp*sizeof(long));
339  ancone = (boolean *)Malloc(chars*sizeof(boolean));
340  anczero = (boolean *)Malloc(chars*sizeof(boolean));
341  ancone0 = (boolean *)Malloc(chars*sizeof(boolean));
342  anczero0 = (boolean *)Malloc(chars*sizeof(boolean));
343  wagner = (bitptr)Malloc(words*sizeof(long));
344  wagner0 = (bitptr)Malloc(words*sizeof(long));
345  place = (long *)Malloc(nonodes*sizeof(long));
346  steps = (bitptr)Malloc(words*sizeof(long));
347  zeroanc = (bitptr)Malloc(words*sizeof(long));
348  oneanc = (bitptr)Malloc(words*sizeof(long));
349  fulzeroanc = (bitptr)Malloc(words*sizeof(long));
350  empzeroanc = (bitptr)Malloc(words*sizeof(long));
351}  /* allocrest */
352
353
354void doinit()
355{
356  /* initializes variables */
357
358  inputnumbers(&spp, &chars, &nonodes, 1);
359  words = chars / bits + 1;
360  getoptions();
361  if (printdata)
362      fprintf(outfile, "%ld species, %ld characters\n\n", spp, chars);
363  alloctree2(&treenode);
364  setuptree2(treenode);
365  allocrest();
366}  /* doinit */
367
368
369void inputoptions()
370{
371  /* input the information on the options */
372  long i;
373
374  if(justwts){
375    if(firstset){
376      scan_eoln(infile);
377      if (ancvar)
378        inputancestorsnew(anczero0, ancone0);
379      if (mixture)
380        inputmixturenew(wagner0);
381      }
382    for (i = 0; i < (chars); i++)
383      weight[i] = 1;
384    inputweights(chars, weight, &weights);
385    for (i = 0; i < (words); i++) {
386      if (mixture)
387        wagner[i] = wagner0[i];
388      else if (allsokal)
389        wagner[i] = 0;
390      else
391        wagner[i] = (1L << (bits + 1)) - (1L << 1);
392    }
393  }
394  else {
395    if (!firstset) {
396      samenumsp(&chars, ith);
397      reallocchars();
398    }
399    scan_eoln(infile);
400    for (i = 0; i < (chars); i++)
401      weight[i] = 1;
402    if (ancvar)
403      inputancestorsnew(anczero0, ancone0);
404    if (mixture)
405      inputmixturenew(wagner0);
406    if (weights)
407      inputweights(chars, weight, &weights);
408    for (i = 0; i < (words); i++) {
409      if (mixture)
410        wagner[i] = wagner0[i];
411      else if (allsokal)
412        wagner[i] = 0;
413      else
414        wagner[i] = (1L << (bits + 1)) - (1L << 1);
415    }
416  }
417  for (i = 0; i < (chars); i++) {
418    if (!ancvar) {
419      anczero[i] = true;
420      ancone[i] = (((1L << (i % bits + 1)) & wagner[i / bits]) != 0);
421    } else {
422      anczero[i] = anczero0[i];
423      ancone[i] = ancone0[i];
424    }
425  }
426  noroot = true;
427  questions = false;
428  for (i = 0; i < (chars); i++) {
429    if (weight[i] > 0) {
430      noroot = (noroot && ancone[i] && anczero[i] &&
431                ((((1L << (i % bits + 1)) & wagner[i / bits]) != 0)
432                 || threshold <= 2.0));
433    }
434    questions = (questions || (ancone[i] && anczero[i]));
435    threshwt[i] = threshold * weight[i];
436  }
437}  /* inputoptions */
438
439
440void doinput()
441{
442  /* reads the input data */
443  inputoptions();
444  if(!justwts || firstset)
445      inputdata2(treenode);
446}  /* doinput */
447
448
449void evaluate(node2 *r)
450{
451  /* Determines the number of steps needed for a tree.
452    This is the minimum number needed to evolve chars on
453    this tree */
454  long i, stepnum, smaller;
455  double sum, term;
456
457  sum = 0.0;
458  for (i = 0; i < (chars); i++) {
459    numszero[i] = 0;
460    numsone[i] = 0;
461  }
462  full = true;
463  for (i = 0; i < (words); i++)
464    zeroanc[i] = fullset;
465  postorder(r, fullset, full, wagner, zeroanc);
466  cpostorder(r, full, zeroanc, numszero, numsone);
467  count(r->fulstte1, zeroanc, numszero, numsone);
468  for (i = 0; i < (words); i++)
469    zeroanc[i] = 0;
470  full = false;
471  postorder(r, fullset, full, wagner, zeroanc);
472  cpostorder(r, full, zeroanc, numszero, numsone);
473  count(r->empstte0, zeroanc, numszero, numsone);
474  for (i = 0; i < (chars); i++) {
475    smaller = spp * weight[i];
476    numsteps[i] = smaller;
477    if (anczero[i]) {
478      numsteps[i] = numszero[i];
479      smaller = numszero[i];
480    }
481    if (ancone[i] && numsone[i] < smaller)
482      numsteps[i] = numsone[i];
483    stepnum = numsteps[i] + extras[i];
484    if (stepnum <= threshwt[i])
485      term = stepnum;
486    else
487      term = threshwt[i];
488    sum += term;
489    if (usertree && which <= maxuser)
490      fsteps[which - 1][i] = term;
491    guess[i] = '?';
492    if (!ancone[i] || (anczero[i] && numszero[i] < numsone[i]))
493      guess[i] = '0';
494    else if (!anczero[i] || (ancone[i] && numsone[i] < numszero[i]))
495      guess[i] = '1';
496  }
497  if (usertree && which <= maxuser) {
498    nsteps[which - 1] = sum;
499    if (which == 1) {
500      minwhich = 1;
501      minsteps = sum;
502    } else if (sum < minsteps) {
503      minwhich = which;
504      minsteps = sum;
505    }
506  }
507  like = -sum;
508}  /* evaluate */
509
510
511void reroot(node2 *outgroup)
512{
513  /* reorients tree, putting outgroup in desired position. */
514  node2 *p, *q;
515
516  if (outgroup->back->index == root->index)
517    return;
518  p = root->next;
519  q = root->next->next;
520  p->back->back = q->back;
521  q->back->back = p->back;
522  p->back = outgroup;
523  q->back = outgroup->back;
524  outgroup->back->back = q;
525  outgroup->back = p;
526}  /* reroot */
527
528
529void savetraverse(node2 *p)
530{
531  /* sets BOOLEANs that indicate which way is down */
532  p->bottom = true;
533  if (p->tip)
534    return;
535  p->next->bottom = false;
536  savetraverse(p->next->back);
537  p->next->next->bottom = false;
538  savetraverse(p->next->next->back);
539}  /* savetraverse */
540
541
542void savetree()
543{
544  /* record in place where each species has to be
545    added to reconstruct this tree */
546  long i, j;
547  node2 *p;
548  boolean done;
549
550  if (noroot)
551    reroot(treenode[outgrno - 1]);
552  savetraverse(root);
553  for (i = 0; i < (nonodes); i++)
554   place[i] = 0;
555  place[root->index - 1] = 1;
556  for (i = 1; i <= (spp); i++) {
557    p = treenode[i - 1];
558    while (place[p->index - 1] == 0) {
559      place[p->index - 1] = i;
560      while (!p->bottom)
561        p = p->next;
562      p = p->back;
563    }
564    if (i > 1) {
565      place[i - 1] = place[p->index - 1];
566      j = place[p->index - 1];
567      done = false;
568      while (!done) {
569        place[p->index - 1] = spp + i - 1;
570        while (!p->bottom)
571          p = p->next;
572        p = p->back;
573        done = (p == NULL);
574        if (!done)
575          done = (place[p->index - 1] != j);
576      }
577    }
578  }
579}  /* savetree */
580
581
582void mix_addtree(long *pos)
583{
584  /* puts tree from ARRAY place in its proper position
585    in ARRAY bestrees */
586  long i;
587  for (i =nextree - 1; i >= (*pos); i--)
588    memcpy(bestrees[i], bestrees[i - 1], spp*sizeof(long));
589  for (i = 0; i < (spp); i++)
590    bestrees[(*pos) - 1][i] = place[i];
591  nextree++;
592}  /* mix_addtree */
593
594
595void mix_findtree(boolean *found, long *pos, long local_nextree,
596                        long *local_place, long **local_bestrees)
597{
598  /* finds tree given by ARRAY local_place in ARRAY local_bestrees by binary search */
599  /* used by dnacomp, dnapars, dollop, mix, & protpars */
600  long i, lower, upper;
601  boolean below, done;
602
603  below = false;
604  lower = 1;
605  upper = local_nextree - 1;
606  (*found) = false;
607  while (!(*found) && lower <= upper) {
608    (*pos) = (lower + upper) / 2;
609    i = 3;
610    done = false;
611    while (!done) {
612      done = (i > spp);
613      if (!done)
614        done = (local_place[i - 1] != local_bestrees[(*pos) - 1][i - 1]);
615      if (!done)
616        i++;
617    }
618    (*found) = (i > spp);
619    below = (local_place[i - 1] <  local_bestrees[(*pos )- 1][i - 1]);
620    if (*found)
621      break;
622    if (below)
623      upper = (*pos) - 1;
624    else
625      lower = (*pos) + 1;
626  }
627  if (!(*found) && !below)
628    (*pos)++;
629}  /* mix_findtree */
630
631
632void tryadd(node2 *p, node2 **item, node2 **nufork)
633{
634  /* temporarily adds one fork and one tip to the tree.
635    if the location where they are added yields greater
636    "likelihood" than other locations tested up to that
637    time, then keeps that location as there */
638
639  long pos;
640  boolean found; 
641  node2 *rute;
642
643  add3(p, *item, *nufork, &root, treenode);
644  evaluate(root);
645  if (lastrearr) {
646    if (like >= bstlike2) {
647      rute = root->next->back;
648      savetree();
649      reroot(rute);
650      if (like > bstlike2) {
651        bestlike = bstlike2 = like;
652        pos = 1;
653        nextree = 1;
654        mix_addtree(&pos);
655      } else {
656        pos = 0;
657        mix_findtree(&found, &pos, nextree, place, bestrees);
658        if (!found) { 
659          if (nextree <= maxtrees)
660            mix_addtree(&pos);
661        }
662      }
663    }
664  }
665  if (like > bestyet) {
666    bestyet = like;
667    there = p;
668  }
669  re_move3(item, nufork, &root, treenode);
670}  /* tryadd */
671
672
673void addpreorder(node2 *p, node2 *item, node2 *nufork)
674{
675  /* traverses a binary tree, calling PROCEDURE tryadd
676    at a node before calling tryadd at its descendants */
677
678  if (p == NULL)
679    return;
680  tryadd(p, &item, &nufork);
681  if (!p->tip) {
682    addpreorder(p->next->back, item, nufork);
683    addpreorder(p->next->next->back, item, nufork);
684  }
685}  /* addpreorder */
686
687
688void tryrearr(node2 *p, node2 **r, boolean *success)
689{
690  /* evaluates one rearrangement of the tree.
691    if the new tree has greater "likelihood" than the old
692    one sets success := TRUE and keeps the new tree.
693    otherwise, restores the old tree */
694  node2 *frombelow, *whereto, *forknode;
695  double oldlike;
696
697  if (p->back == NULL)
698    return;
699  forknode = treenode[p->back->index - 1];
700  if (forknode->back == NULL)
701    return;
702  oldlike = bestyet;
703  if (p->back->next->next == forknode)
704    frombelow = forknode->next->next->back;
705  else
706    frombelow = forknode->next->back;
707  whereto = treenode[forknode->back->index - 1];
708  re_move3(&p, &forknode, &root, treenode);
709  add3(whereto, p, forknode, &root, treenode);
710  evaluate(*r);
711  if (like <= oldlike) {
712    re_move3(&p, &forknode, &root, treenode);
713    add3(frombelow, p, forknode, &root, treenode);
714  } else {
715    *success = true;
716    bestyet = like;
717  }
718}  /* tryrearr */
719
720
721void repreorder(node2 *p, node2 **r, boolean *success)
722{
723  /* traverses a binary tree, calling PROCEDURE tryrearr
724    at a node before calling tryrearr at its descendants */
725  if (p == NULL)
726    return;
727  tryrearr(p, r, success);
728  if (!p->tip) {
729    repreorder(p->next->back, r,success);
730    repreorder(p->next->next->back, r,success);
731  }
732}  /* repreorder */
733
734
735void rearrange(node2 **r)
736{
737  /* traverses the tree (preorder), finding any local
738    rearrangement which decreases the number of steps.
739    if traversal succeeds in increasing the tree's
740    "likelihood", PROCEDURE rearrange runs traversal again */
741  boolean success=true;
742  while (success) {
743    success = false;
744    repreorder(*r,r,&success);
745  }
746}  /* rearrange */
747
748
749void mix_addelement(node2 **p, long *nextnode, long *lparens,
750                        boolean *names)
751{
752  /* recursive procedure adds nodes to user-defined tree */
753  node2 *q;
754  long i, n;
755  boolean found;
756  Char str[nmlngth];
757
758  getch(&global_ch, lparens, intree);
759  if (global_ch == '(' ) {
760    if ((*lparens) >= spp) {
761      printf("\n\nERROR IN USER TREE: Too many left parentheses\n\n");
762      exxit(-1);
763    }
764    (*nextnode)++;
765    q = treenode[(*nextnode) - 1];
766    mix_addelement(&q->next->back, nextnode, lparens, names);
767    q->next->back->back = q->next;
768    findch(',', &global_ch, which);
769    mix_addelement(&q->next->next->back, nextnode, lparens, names);
770    q->next->next->back->back = q->next->next;
771    findch(')', &global_ch, which);
772    *p = q;
773    return;
774  }
775  for (i = 0; i < nmlngth; i++)
776    str[i] = ' ';
777  n = 1;
778  do {
779    if (global_ch == '_')
780      global_ch = ' ';
781    str[n - 1] =global_ch;
782    if (eoln(intree)) 
783      scan_eoln(intree);
784    global_ch = gettc(intree);
785    if (global_ch == '\n')
786      global_ch = ' ';
787    n++;
788  } while (global_ch != ',' && global_ch != ')' && global_ch != ':' && n <= nmlngth);
789  n = 1;
790  do {
791    found = true;
792    for (i = 0; i < nmlngth; i++)
793      found = (found && ((str[i] == nayme[n - 1][i]) ||
794                         ((nayme[n - 1][i] == '_') && (str[i] == ' '))));
795    if (found) {
796      if (names[n - 1] == false) {
797        *p = treenode[n - 1];
798        names[n - 1] = true;
799      } else {
800        printf("\n\nERROR IN USER TREE: Duplicate name found: ");
801        for (i = 0; i < nmlngth; i++)
802          putchar(nayme[n - 1][i]);
803        printf("\n\n");
804        exxit(-1);
805      }
806    } else
807      n++;
808  } while (!(n > spp || found ));
809  if (n <= spp)
810    return;
811  printf("CANNOT FIND SPECIES: ");
812  for (i = 0; i < nmlngth; i++)
813    putchar(str[i]);
814  putchar('\n');
815}  /* mix_addelement */
816
817
818void mix_treeread()
819{
820  /* read in user-defined tree and set it up */
821  long nextnode, lparens, i;
822  boolean *names;
823
824  root = treenode[spp];
825  nextnode = spp;
826  root->back = NULL;
827  names = (boolean *)Malloc(spp*sizeof(boolean));
828  for (i = 0; i < (spp); i++)
829    names[i] = false;
830  lparens = 0;
831  mix_addelement(&root, &nextnode, &lparens, names);
832  if (global_ch == '[') {
833    do
834      global_ch = gettc(intree);
835    while (global_ch != ']');
836    global_ch = gettc(intree);
837  }
838  findch(';', &global_ch, which);
839  if (progress)
840    printf(".");
841  scan_eoln(intree);
842  free(names);
843}  /* mix_treeread */
844
845
846void describe()
847{
848  /* prints ancestors, steps and table of numbers of steps in
849    each character */
850
851  if (treeprint)
852    fprintf(outfile, "\nrequires a total of %10.3f\n", -like);
853  putc('\n', outfile);
854  if (stepbox)
855    writesteps(weights, numsteps);
856  if (questions && (!noroot || didreroot))
857    guesstates(guess);
858  if (ancseq) {
859    hypstates(fullset, full, noroot, didreroot, root, wagner,
860                zeroanc, oneanc, treenode, guess, garbage);
861    putc('\n', outfile);
862  }
863  putc('\n', outfile);
864  if (trout) {
865    col = 0;
866    treeout2(root, &col, root);
867  }
868}  /* describe */
869
870
871void maketree()
872{
873  /* constructs a binary tree from the pointers in treenode.
874    adds each node at location which yields highest "likelihood"
875    then rearranges the tree for greatest "likelihood" */
876  long i, j, numtrees;
877  double gotlike;
878  node2 *item, *nufork, *dummy;
879
880  fullset = (1L << (bits + 1)) - (1L << 1);
881  for (i=0 ; i<words ; ++i){
882    fulzeroanc[i]=fullset;
883    empzeroanc[i]= 0;}
884  if (!usertree) {
885    for (i = 1; i <= (spp); i++)
886      enterorder[i - 1] = i;
887    if (jumble)
888      randumize(seed, enterorder);
889    root = treenode[enterorder[0] - 1];
890    add3(treenode[enterorder[0] - 1], treenode[enterorder[1] - 1],
891        treenode[spp], &root, treenode);
892    if (progress) {
893      printf("Adding species:\n");
894      writename(0, 2, enterorder);
895#ifdef WIN32
896      phyFillScreenColor();
897#endif
898    }
899    lastrearr = false;
900    for (i = 3; i <= (spp); i++) {
901      bestyet = -350.0 * spp * chars;
902      item = treenode[enterorder[i - 1] - 1];
903      nufork = treenode[spp + i - 2];
904      addpreorder(root, item, nufork);
905      add3(there, item, nufork, &root, treenode);
906      like = bestyet;
907      rearrange(&root);
908      if (progress) {
909        writename(i - 1, 1, enterorder);
910#ifdef WIN32
911        phyFillScreenColor();
912#endif
913      }
914      lastrearr = (i == spp);
915      if (lastrearr) {
916        if (progress) {
917          printf("\nDoing global rearrangements\n");
918          printf("  !");
919          for (j = 1; j <= (nonodes); j++)
920            putchar('-');
921          printf("!\n");
922#ifdef WIN32
923          phyFillScreenColor();
924#endif
925        }
926        bestlike = bestyet;
927        if (jumb == 1) {
928          bstlike2 = bestlike;
929          nextree = 1;
930        }
931        do {
932          if (progress)
933            printf("   ");
934          gotlike = bestlike;
935          for (j = 0; j < (nonodes); j++) {
936            bestyet = -350.0  * spp * chars;
937            item = treenode[j];
938            if (item != root) {
939              re_move3(&item, &nufork, &root, treenode);
940              there = root;
941              addpreorder(root, item, nufork);
942              add3(there, item, nufork, &root, treenode);
943            }
944            if (progress) {
945              putchar('.');
946              fflush(stdout);
947            }
948          }
949          if (progress) {
950            putchar('\n');
951#ifdef WIN32
952            phyFillScreenColor();
953#endif
954          }
955        } while (bestlike > gotlike);
956      }
957    }
958    if (progress)
959      putchar('\n');
960    for (i = spp - 1; i >= 1; i--)
961      re_move3(&treenode[i], &dummy, &root, treenode);
962    if (jumb == njumble) {
963      if (treeprint) {
964        putc('\n', outfile);
965        if (nextree == 2)
966          fprintf(outfile, "One most parsimonious tree found:\n");
967        else
968          fprintf(outfile, "%6ld trees in all found\n", nextree - 1);
969      }
970      if (nextree > maxtrees + 1) {
971        if (treeprint)
972          fprintf(outfile, "here are the first%4ld of them\n",(long)maxtrees);
973        nextree = maxtrees + 1;
974      }
975      if (treeprint)
976        putc('\n', outfile);
977      for (i = 0; i <= (nextree - 2); i++) {
978        root = treenode[0];
979        add3(treenode[0], treenode[1], treenode[spp], &root, treenode);
980        for (j = 3; j <= (spp); j++)
981            add3(treenode[bestrees[i][j - 1] - 1], treenode[j - 1],
982            treenode[spp + j - 2], &root, treenode);
983        if (noroot)
984          reroot(treenode[outgrno - 1]);
985        didreroot = (outgropt && noroot);
986        evaluate(root);
987        printree(treeprint, noroot, didreroot, root);
988        describe();
989        for (j = 1; j < (spp); j++)
990          re_move3(&treenode[j], &dummy, &root, treenode);
991      }
992    }
993  } else {
994    openfile(&intree,INTREE,"input tree file", "r",progname,intreename);
995    numtrees = countsemic(&intree);
996    if (numtrees > 2)
997      initseed(&inseed, &inseed0, seed);
998    if (treeprint) {
999      fprintf(outfile, "User-defined tree");
1000      if (numtrees > 1)
1001        putc('s', outfile);
1002      fprintf(outfile, ":\n\n");
1003    }
1004    which = 1;
1005    if (progress)
1006      printf("   ");
1007    while (which <= numtrees ) {
1008      mix_treeread();
1009      didreroot = (outgropt && noroot);
1010      if (noroot)
1011        reroot(treenode[outgrno - 1]);
1012      evaluate(root);
1013      printree(treeprint, noroot, didreroot, root);
1014      describe();
1015      which++;
1016    }
1017    if (progress)
1018      printf("\n");
1019    FClose(intree);
1020    fprintf(outfile, "\n\n");
1021    if (numtrees > 2 && chars > 1 ) {
1022      if (progress)
1023        printf("   sampling for SH test\n");
1024      standev(numtrees, minwhich, minsteps, nsteps, fsteps, seed);
1025    }
1026  }
1027  if (jumb == njumble) {
1028    if (progress) {
1029      printf("\nOutput written to file \"%s\"\n\n", outfilename);
1030      if (trout)
1031        printf("Trees also written onto file \"%s\"\n", outtreename);
1032      putchar('\n');
1033    }
1034  }
1035  if (ancseq)
1036    freegarbage(&garbage);
1037}  /* maketree */
1038
1039int main(int argc, Char *argv[])
1040{  /* Mixed parsimony by uphill search */
1041#ifdef MAC
1042  argc = 1;         /* macsetup("Mix","");                */
1043  argv[0] = "Mix";
1044#endif
1045  init(argc, argv);
1046  progname = argv[0];
1047  openfile(&infile,INFILE,"input file", "r",argv[0],infilename);
1048  openfile(&outfile,OUTFILE,"output file", "w",argv[0],outfilename);
1049
1050  ibmpc = IBMCRT;
1051  ansi = ANSICRT;
1052  mulsets = false;
1053  msets = 1;
1054  firstset = true;
1055  garbage = NULL;
1056  bits = 8*sizeof(long) - 1;
1057  doinit();
1058  if (weights || justwts)
1059    openfile(&weightfile,WEIGHTFILE,"weights file","r",argv[0],weightfilename);
1060  if (trout)
1061    openfile(&outtree,OUTTREE,"output tree file", "w",argv[0],outtreename);
1062  if(ancvar)
1063      openfile(&ancfile,ANCFILE,"ancestors file", "r",argv[0],ancfilename);
1064  if(mixture)
1065      openfile(&mixfile,MIXFILE,"mixture file", "r",argv[0],mixfilename);
1066 
1067  for (ith = 1; ith <= msets; ith++) {
1068    if(firstset){
1069        if (allsokal && !mixture)
1070            fprintf(outfile, "Camin-Sokal parsimony method\n\n");
1071        if (allwagner && !mixture)
1072            fprintf(outfile, "Wagner parsimony method\n\n");
1073    }
1074    doinput();
1075    if (msets > 1 && !justwts) {
1076      fprintf(outfile, "Data set # %ld:\n\n",ith);
1077      if (progress)
1078        printf("\nData set # %ld:\n",ith);
1079    }
1080    if (justwts){
1081        if(firstset && mixture && printdata)
1082            printmixture(outfile, wagner);
1083        fprintf(outfile, "Weights set # %ld:\n\n", ith);
1084        if (progress)
1085            printf("\nWeights set # %ld:\n\n", ith);
1086    }
1087    else if (mixture && printdata)
1088        printmixture(outfile, wagner);
1089    if (printdata){
1090        if (weights || justwts)
1091            printweights(outfile, 0, chars, weight, "Characters");
1092        if (ancvar)
1093            printancestors(outfile, anczero, ancone);
1094    }
1095   
1096
1097    if (ith == 1)
1098      firstset = false;
1099    for (jumb = 1; jumb <= njumble; jumb++)
1100      maketree();
1101  }
1102  free(place);
1103  free(steps);
1104  free(zeroanc);
1105  free(oneanc);
1106  free(fulzeroanc);
1107  free(empzeroanc);
1108  FClose(outfile);
1109  FClose(infile);
1110  FClose(outtree);
1111#ifdef MAC
1112  fixmacfile(outtreename);
1113  fixmacfile(outfilename);
1114#endif
1115#ifdef WIN32
1116  phyRestoreConsoleAttributes();
1117#endif
1118  return 0;
1119}  /* Mixed parsimony by uphill search */
Note: See TracBrowser for help on using the repository browser.