source: trunk/GDE/PHYLIP/dnamove.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: 46.2 KB
Line 
1
2#include "phylip.h"
3#include "moves.h"
4#include "seq.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 overr           4
12#define which           1
13
14typedef enum {
15  horiz, vert, up, overt, upcorner, downcorner, aa, cc, gg, tt, question
16  } chartype;
17
18typedef enum {
19  rearr, flipp, reroott, none
20  } rearrtype;
21
22typedef struct gbase2 {
23  baseptr2 base2;
24  struct gbase2 *next;
25} gbase2;
26
27typedef enum {
28  arb, use, spec
29  } howtree;
30
31typedef node **pointptr;
32
33#ifndef OLDC
34/* function prototypes */
35void dnamove_gnu(gbase2 **);
36void dnamove_chuck(gbase2 *);
37void getoptions(void);
38void inputoptions(void);
39void dnamove_inputdata(void);
40void allocrest(void);
41void doinput(void);
42void configure(void);
43void prefix(chartype);
44void postfix(chartype);
45       
46void makechar(chartype);
47void dnamove_add(node *, node *, node *);
48void dnamove_re_move(node **, node **);
49void dnamove_fillin(node *);
50void dnamove_postorder(node *);
51void evaluate(node *);
52void dnamove_reroot(node *);
53void dnamove_ancestset(long, long, long *);
54void firstrav(node *, long);
55void dnamove_hyptrav(node *, long *, long, boolean *);
56
57void dnamove_hypstates(void);
58void grwrite(chartype, long, long *);
59void dnamove_drawline(long);
60void dnamove_printree(void);
61void arbitree(void);
62void yourtree(void);
63void initdnamovenode(node **, node **, node *, long, long, long *,
64        long *, initops, pointarray, pointarray, Char *, Char *,
65        FILE *);
66void buildtree(void);
67void setorder(void);
68void mincomp(void);
69       
70void rearrange(void);
71void dnamove_nextinc(void);
72void dnamove_nextchar(void);
73void dnamove_prevchar(void);
74void dnamove_show(void);
75void tryadd(node *, node **, node **, double *);
76void addpreorder(node *, node *, node *, double *);
77void try(void);
78void undo(void);
79void treewrite(boolean);
80       
81void clade(void);
82void flip(void);
83void changeoutgroup(void);
84void redisplay(void);
85void treeconstruct(void);
86/* function prototypes */
87#endif
88
89
90
91char infilename[FNMLNGTH],intreename[FNMLNGTH],outtreename[FNMLNGTH], weightfilename[FNMLNGTH];
92node *root;
93long chars, screenlines, col, treelines, leftedge, topedge, vmargin,
94   hscroll, vscroll, scrollinc, screenwidth, farthest;
95boolean weights, thresh, waswritten;
96boolean usertree, goteof, firsttree, haslengths;  /*treeread variables*/
97pointarray nodep;                                  /*treeread variables*/
98node *grbg = NULL;                                  /*treeread variables*/
99long *zeros;                                          /*treeread variables*/
100pointptr treenode;   /* pointers to all nodes in tree */
101double threshold;
102double *threshwt;
103boolean reversed[11];
104boolean graphic[11];
105unsigned char chh[11];
106howtree how;
107gbase2 *garbage;
108char *progname;
109
110/* Local variables for treeconstruct, propogated global for C version: */
111
112long dispchar, atwhat, what, fromwhere, towhere, oldoutgrno, compatible;
113double like, bestyet, gotlike;
114boolean display, newtree, changed, subtree, written, oldwritten, restoring,
115  wasleft, oldleft, earlytree;
116steptr necsteps;
117boolean *in_tree;
118long sett[31];
119steptr numsteps;
120node *nuroot;
121rearrtype lastop;
122Char  global_ch;
123boolean *names;
124
125
126void dnamove_gnu(gbase2 **p)
127{
128  /* this and the following are do-it-yourself garbage collectors.
129     Make a new node or pull one off the garbage list */
130  if (garbage != NULL) {
131    *p = garbage;
132    garbage = garbage->next;
133  } else {
134    *p = (gbase2 *)Malloc(sizeof(gbase2));
135    (*p)->base2 = (baseptr2)Malloc(chars*sizeof(long));
136  }
137  (*p)->next = NULL;
138}  /* dnamove_gnu */
139
140
141void dnamove_chuck(gbase2 *p)
142{
143  /* collect garbage on p -- put it on front of garbage list */
144  p->next = garbage;
145  garbage = p;
146}  /* dnamove_chuck */
147
148void getoptions()
149{
150  /* interactively set options */
151  Char ch;
152  boolean done, gotopt;
153  long loopcount;
154
155  how = arb;
156  usertree = false;
157  goteof = false;
158  outgrno = 1;
159  outgropt = false;
160  thresh = false;
161  weights = false;
162  interleaved = true;
163  loopcount = 0;
164  do {
165    cleerhome();
166    printf("\nInteractive DNA parsimony, version %s\n\n",VERSION);
167    printf("Settings for this run:\n");
168    printf("  O                             Outgroup root?");
169    if (outgropt)
170      printf("  Yes, at sequence number%3ld\n", outgrno);
171    else
172      printf("  No, use as outgroup species%3ld\n", outgrno);
173    printf("  W                            Sites weighted?  %s\n",
174           (weights ? "Yes" : "No"));
175    printf("  T                   Use Threshold parsimony?");
176    if (thresh)
177      printf("  Yes, count up to%4.1f per site\n", threshold);
178    else
179      printf("  No, use ordinary parsimony\n");
180    printf("  I               Input sequences interleaved?  %s\n",
181           (interleaved ? "Yes" : "No, sequential"));
182
183    printf("  U   Initial tree (arbitrary, user, specify)?  %s\n",
184           (how == arb) ? "Arbitrary"                :
185           (how == use) ? "User tree from tree file" : "Tree you specify");
186    printf("  0        Graphics type (IBM PC, ANSI, none)?  %s\n",
187           ibmpc ? "IBM PC\n" : ansi  ? "ANSI"     : "(none)");
188    printf("  S                  Width of terminal screen?");
189    printf("%4ld\n", screenwidth);
190    printf("  L                 Number of lines on screen?%4ld\n",screenlines);
191    printf("\nAre these settings correct? ");
192    printf("(type Y or the letter for one to change)\n");
193#ifdef WIN32
194    phyFillScreenColor();
195#endif
196    scanf("%c%*[^\n]", &ch);
197    getchar();
198    uppercase(&ch);
199    done = (ch == 'Y');
200    gotopt = (strchr("SOTIU0WL",ch) != NULL) ? true : false;
201    if (gotopt) {
202      switch (ch) {
203       
204      case 'O':
205        outgropt = !outgropt;
206        if (outgropt)
207          initoutgroup(&outgrno, spp);
208        break;
209       
210        case 'T':
211          thresh = !thresh;
212          if (thresh)
213            initthreshold(&threshold);
214          break;
215       
216        case 'I':
217          interleaved = !interleaved;
218          break;
219
220        case 'W':
221          weights = !weights;
222          break;
223
224        case 'U':
225          if (how == arb){
226            how = use;
227            usertree = 1;}
228          else if (how == use){
229            how = spec;
230            usertree = 0;}
231          else
232            how = arb;
233          break;
234       
235        case '0':
236          initterminal(&ibmpc, &ansi);
237          break;
238       
239        case 'S':
240          screenwidth= readlong("Width of terminal screen (in characters)?\n");
241          break;
242
243        case 'L':
244          initnumlines(&screenlines);
245          break;
246        }
247      }
248    if (!(gotopt || done))
249      printf("Not a possible option!\n");
250    countup(&loopcount, 100);   
251  } while (!done);
252  if (scrollinc < screenwidth / 2.0)
253    hscroll = scrollinc;
254  else
255    hscroll = screenwidth / 2;
256  if (scrollinc < screenlines / 2.0)
257    vscroll = scrollinc;
258  else
259    vscroll = screenlines / 2;
260}  /* getoptions */
261
262
263void inputoptions()
264{
265  /* input the information on the options */
266  long i;
267
268  for (i = 0; i < (chars); i++)
269    weight[i] = 1;
270  if (weights){
271      inputweights(chars, weight, &weights);
272      printweights(stdout, 0, chars, weight, "Sites");
273  }
274  if (!thresh)
275    threshold = spp;
276  for (i = 0; i < (chars); i++)
277    threshwt[i] = threshold * weight[i];
278}  /* inputoptions */
279
280
281void dnamove_inputdata()
282{
283  /* input the names and sequences for each species */
284  long i,j, basesread, basesnew=0;
285  Char charstate;
286  boolean allread, done;
287  long ns = 0;   /* temporary base set for input */
288
289  basesread = 0;
290  allread = false;
291  while (!(allread)) {
292    if (eoln(infile)) 
293      scan_eoln(infile);
294    i = 1;
295    while (i <= spp ) {
296      if ((interleaved && basesread == 0) || !interleaved)
297        initname(i - 1);
298      if (interleaved)
299        j = basesread;
300      else
301        j = 0;
302      done = false;
303      while (!done && !eoff(infile)) {
304        if (interleaved)
305          done = true;
306        while (j < chars && !(eoln(infile) || eoff(infile))) {
307          charstate = gettc(infile);
308          if (global_ch == '\n')
309            global_ch = ' ';
310          if (charstate == ' ' || (charstate >= '0' && charstate <= '9'))
311            continue;
312          uppercase(&charstate);
313          if ((strchr("ABCDGHKMNRSTUVWXY?O-",charstate)) == NULL){
314            printf("ERROR: bad base: %c at site %5ld of species %3ld\n",
315                   charstate, j+1, i);
316            if (charstate == '.') {
317              printf("       Periods (.) may not be used as gap characters.\n");
318              printf("       The correct gap character is (-)\n");
319            }
320            exxit(-1);
321          }
322          j++;
323          switch (charstate) {
324       
325          case 'A':
326            ns = 1L << ((long)A);
327            break;
328       
329          case 'C':
330            ns = 1L << ((long)C);
331            break;
332       
333          case 'G':
334            ns = 1L << ((long)G);
335            break;
336       
337          case 'U':
338            ns = 1L << ((long)T);
339            break;
340       
341          case 'T':
342            ns = 1L << ((long)T);
343            break;
344       
345          case 'M':
346            ns = (1L << ((long)A)) | (1L << ((long)C));
347            break;
348       
349          case 'R':
350            ns = (1L << ((long)A)) | (1L << ((long)G));
351            break;
352       
353          case 'W':
354            ns = (1L << ((long)A)) | (1L << ((long)T));
355            break;
356       
357          case 'S':
358            ns = (1L << ((long)C)) | (1L << ((long)G));
359            break;
360       
361          case 'Y':
362            ns = (1L << ((long)C)) | (1L << ((long)T));
363            break;
364       
365          case 'K':
366            ns = (1L << ((long)G)) | (1L << ((long)T));
367            break;
368       
369          case 'B':
370            ns = (1L << ((long)C)) | (1L << ((long)G)) | (1L << ((long)T));
371            break;
372       
373          case 'D':
374            ns = (1L << ((long)A)) | (1L << ((long)G)) | (1L << ((long)T));
375            break;
376       
377          case 'H':
378            ns = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)T));
379            break;
380       
381          case 'V':
382            ns = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)G));
383            break;
384       
385          case 'N':
386            ns = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)G)) |
387              (1L << ((long)T));
388            break;
389       
390          case 'X':
391            ns = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)G)) |
392              (1L << ((long)T));
393            break;
394       
395          case '?':
396            ns = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)G)) |
397                    (1L << ((long)T)) | (1L << ((long)O));
398            break;
399       
400          case 'O':
401            ns = 1L << ((long)O);
402            break;
403       
404          case '.':
405            ns = treenode[0]->base2[j - 1];
406            break;
407       
408          case '-':
409            ns = 1L << ((long)O);
410            break;
411          }
412          treenode[i - 1]->base2[j - 1] = ns;
413        }
414        if (interleaved)
415          continue;
416        if (j < chars)
417          scan_eoln(infile);
418        else if (j == chars)
419          done = true;
420      }
421      if (interleaved && i == 1)
422        basesnew = j;
423      scan_eoln(infile);
424      if ((interleaved && j != basesnew) || ((!interleaved) && j != chars)){
425        printf("ERROR: sequences out of alignment\n");
426        exxit(-1);}
427      i++;
428    }
429    if (interleaved) {
430      basesread = basesnew;
431      allread = (basesread == chars);
432    } else
433      allread = (i > spp);
434  }
435  root = NULL;
436  printf("\n\n");
437}  /* dnamove_inputdata */
438
439
440void allocrest()
441{
442  nayme = (naym *)Malloc(spp*sizeof(naym));
443  in_tree = (boolean *)Malloc(nonodes*sizeof(boolean));
444  weight = (steptr)Malloc(chars*sizeof(long));
445  numsteps = (steptr)Malloc(chars*sizeof(long));
446  necsteps = (steptr)Malloc(chars*sizeof(long));
447  threshwt = (double *)Malloc(chars*sizeof(double));
448}  /* allocrest */
449
450
451void doinput()
452{
453  /* reads the input data */
454  long i, j;
455  node *p;
456
457  inputnumbers(&spp, &chars, &nonodes, 1);
458  printf("%2ld species, %3ld  sites\n", spp, chars);
459  getoptions();
460  printf("\nReading input file ...\n\n");
461  if (weights)
462    openfile(&weightfile,WEIGHTFILE,"weights file","r",progname,weightfilename);
463  allocrest();
464  inputoptions();
465  alloctree(&treenode, nonodes, usertree);
466  for (i = 0; i < (spp); i++)
467    treenode[i]->base2 = (baseptr2)Malloc(chars*sizeof(long));
468  if(!usertree) {
469        for (i = spp; i < (nonodes); i++) {
470          p = treenode[i];
471          for (j = 1; j <= 3; j++) {
472            p->base2 = (baseptr2)Malloc(chars*sizeof(long));
473            p = p->next;
474          }
475        }
476  } 
477  setuptree(treenode, nonodes, usertree);
478  dnamove_inputdata();
479}  /* doinput */
480
481
482void configure()
483{
484  /* configure to machine -- set up special characters */
485  chartype a;
486
487  for (a = horiz; (long)a <= (long)question; a = (chartype)((long)a + 1))
488    reversed[(long)a] = false;
489  for (a = horiz; (long)a <= (long)question; a = (chartype)((long)a + 1))
490    graphic[(long)a] = false;
491  if (ibmpc) {
492    chh[(long)horiz] = 205;
493    graphic[(long)horiz] = true;
494    chh[(long)vert] = 186;
495    graphic[(long)vert] = true;
496    chh[(long)up] = 186;
497    graphic[(long)up] = true;
498    chh[(long)overt] = 205;
499    graphic[(long)overt] = true;
500    chh[(long)upcorner] = 200;
501    graphic[(long)upcorner] = true;
502    chh[(long)downcorner] = 201;
503    graphic[(long)downcorner] = true;
504    chh[(long)aa] = 176;
505    chh[(long)cc] = 178;
506    chh[(long)gg] = 177;
507    chh[(long)tt] = 219;
508    chh[(long)question] = '\001';
509    return;
510  }
511  if (ansi) {
512    chh[(long)horiz] = ' ';
513    reversed[(long)horiz] = true;
514    chh[(long)vert] = chh[(long)horiz];
515    reversed[(long)vert] = true;
516    chh[(long)up] = 'x';
517    graphic[(long)up] = true;
518    chh[(long)overt] = 'q';
519    graphic[(long)overt] = true;
520    chh[(long)upcorner] = 'm';
521    graphic[(long)upcorner] = true;
522    chh[(long)downcorner] = 'l';
523    graphic[(long)downcorner] = true;
524    chh[(long)aa] = 'a';
525    reversed[(long)aa] = true;
526    chh[(long)cc] = 'c';
527    reversed[(long)cc] = true;
528    chh[(long)gg] = 'g';
529    reversed[(long)gg] = true;
530    chh[(long)tt] = 't';
531    reversed[(long)tt] = true;
532    chh[(long)question] = '?';
533    reversed[(long)question] = true;
534    return;
535  }
536  chh[(long)horiz] = '=';
537  chh[(long)vert] = ' ';
538  chh[(long)up] = '!';
539  chh[(long)upcorner] = '`';
540  chh[(long)downcorner] = ',';
541  chh[(long)overt] = '-';
542  chh[(long)aa] = 'a';
543  chh[(long)cc] = 'c';
544  chh[(long)gg] = 'g';
545  chh[(long)tt] = 't';
546  chh[(long)question] = '.';
547}  /* configure */
548
549
550void prefix(chartype a)
551{
552  /* give prefix appropriate for this character */
553  if (reversed[(long)a])
554    prereverse(ansi);
555  if (graphic[(long)a])
556    pregraph2(ansi);
557}  /* prefix */
558
559
560void postfix(chartype a)
561{
562  /* give postfix appropriate for this character */
563  if (reversed[(long)a])
564    postreverse(ansi);
565  if (graphic[(long)a])
566    postgraph2(ansi);
567}  /* postfix */
568
569
570void makechar(chartype a)
571{
572  /* print out a character with appropriate prefix or postfix */
573  prefix(a);
574  putchar(chh[(long)a]);
575  postfix(a);
576}  /* makechar */
577
578
579void dnamove_add(node *below, node *newtip, node *newfork)
580{
581  /* inserts the nodes newfork and its left descendant, newtip,
582     to the tree.  below becomes newfork's right descendant */
583  boolean putleft;
584  node *leftdesc, *rtdesc;
585
586  if (below != treenode[below->index - 1])
587    below = treenode[below->index - 1];
588  if (below->back != NULL)
589    below->back->back = newfork;
590  newfork->back = below->back;
591  putleft = true;
592  if (restoring)
593    putleft = wasleft;
594  if (putleft) {
595    leftdesc = newtip;
596    rtdesc = below;
597  } else {
598    leftdesc = below;
599    rtdesc = newtip;
600  }
601  rtdesc->back = newfork->next->next;
602  newfork->next->next->back = rtdesc;
603  newfork->next->back = leftdesc;
604  leftdesc->back = newfork->next;
605  if (root == below)
606    root = newfork;
607  root->back = NULL;
608}  /* dnamove_add */
609
610
611void dnamove_re_move(node **item, node **fork)
612{
613  /* removes nodes item and its ancestor, fork, from the tree.
614     the new descendant of fork's ancestor is made to be
615     fork's second descendant (other than item).  Also
616     returns pointers to the deleted nodes, item and fork */
617  node *p, *q;
618
619  if ((*item)->back == NULL) {
620    *fork = NULL;
621    return;
622  }
623  *fork = treenode[(*item)->back->index - 1];
624  if (*item == (*fork)->next->back) {
625    if (root == *fork)
626      root = (*fork)->next->next->back;
627    wasleft = true;
628  } else {
629    if (root == *fork)
630      root = (*fork)->next->back;
631    wasleft = false;
632  }
633  p = (*item)->back->next->back;
634  q = (*item)->back->next->next->back;
635  if (p != NULL)
636    p->back = q;
637  if (q != NULL)
638    q->back = p;
639  (*fork)->back = NULL;
640  p = (*fork)->next;
641  while (p != *fork) {
642    p->back = NULL;
643    p = p->next;
644  }
645  (*item)->back = NULL;
646}  /* dnamove_re_move */
647
648
649void dnamove_fillin(node *p)
650{
651  /* sets up for each node in the tree the base sequence
652     at that point and counts the changes.  The program
653     spends much of its time in this PROCEDURE */
654  long i;
655  long ns, rs, ls;
656
657  for (i = 0; i < (chars); i++) {
658    ls = p->next->back->base2[i];
659    rs = p->next->next->back->base2[i];
660    ns = ls & rs;
661    if (ns == 0) {
662      ns = ls | rs;
663      numsteps[i] += weight[i];
664    }
665    p->base2[i] = ns;
666  }
667}  /* dnamove_fillin */
668
669
670void dnamove_postorder(node *p)
671{
672  /* traverses a binary tree, calling PROCEDURE fillin at a
673     node's descendants before calling fillin at the node */
674
675  if (p->tip)
676    return;
677  dnamove_postorder(p->next->back);
678  dnamove_postorder(p->next->next->back);
679  dnamove_fillin(p);
680}  /* dnamove_postorder */
681
682
683void evaluate(node *r)
684{
685  /* determines the number of steps needed for a tree. this is
686     the minimum number of steps needed to evolve sequences on
687     this tree */
688  long i, steps;
689  double sum;
690
691  compatible = 0;
692  sum = 0.0;
693  for (i = 0; i < (chars); i++)
694    numsteps[i] = 0;
695  dnamove_postorder(r);
696  for (i = 0; i < (chars); i++) {
697    steps = numsteps[i];
698    if (steps <= threshwt[i])
699      sum += steps;
700    else
701      sum += threshwt[i];
702    if (steps <= necsteps[i] && !earlytree)
703      compatible += weight[i];
704  }
705  like = -sum;
706  /*printf("like: %f\n",like);*/
707}  /* evaluate */
708
709
710void dnamove_reroot(node *outgroup)
711{
712  /* reorients tree, putting outgroup in desired position. */
713  node *p, *q, *newbottom, *oldbottom;
714  boolean onleft;
715
716  if (outgroup->back->index == root->index)
717    return;
718  newbottom = outgroup->back;
719  p = treenode[newbottom->index - 1]->back;
720  while (p->index != root->index) {
721    oldbottom = treenode[p->index - 1];
722    treenode[p->index - 1] = p;
723    p = oldbottom->back;
724  }
725  onleft = (p == root->next);
726  if (restoring)
727    if (!onleft && wasleft){
728      p = root->next->next;
729      q = root->next;
730    } else {
731      p = root->next;
732      q = root->next->next;
733    }
734  else {
735    if (onleft)
736      oldoutgrno = root->next->next->back->index;
737    else
738      oldoutgrno = root->next->back->index;
739    wasleft = onleft;
740    p = root->next;
741    q = root->next->next;
742  }
743  p->back->back = q->back;
744  q->back->back = p->back;
745  p->back = outgroup;
746  q->back = outgroup->back;
747  if (restoring) {
748    if (!onleft && wasleft) {
749      outgroup->back->back = root->next;
750      outgroup->back = root->next->next;
751    } else {
752      outgroup->back->back = root->next->next;
753      outgroup->back = root->next;
754    }
755  } else {
756    outgroup->back->back = root->next->next;
757    outgroup->back = root->next;
758  }
759  treenode[newbottom->index - 1] = newbottom;
760}  /* dnamove_reroot */
761
762
763void dnamove_ancestset(long a, long b, long *c)
764{
765  /* make the set of ancestral states below nodes
766     whose base sets are a and b */
767  *c = a & b;
768  if (*c == 0)
769    *c = a | b;
770}  /* dnamove_ancestset */
771
772
773void firstrav(node *r, long i)
774{
775  /* initial traverse for hypothetical states */
776  if (r->tip)
777    return;
778  firstrav(r->next->back, i);
779  firstrav(r->next->next->back, i);
780  dnamove_ancestset(r->next->back->base2[i - 1],
781            r->next->next->back->base2[i - 1], &r->base2[i - 1]);
782}  /* firstrav */
783
784
785void dnamove_hyptrav(node *r, long *hypset, long i, boolean *local_bottom)
786{
787  /*  compute, print out state at one interior node */
788  long tempset, local_left, rt, anc;
789  gbase2 *temparray, *ancset;
790
791  dnamove_gnu(&ancset);
792  dnamove_gnu(&temparray);
793  anc = hypset[i - 1];
794  if (!r->tip) {
795    local_left = r->next->back->base2[i - 1];
796    rt = r->next->next->back->base2[i - 1];
797    tempset = local_left & rt & anc;
798    if (tempset == 0) {
799      tempset = (local_left & rt) | (local_left & anc) | (rt & anc);
800      if (tempset == 0)
801        tempset = local_left | rt | anc;
802    }
803    r->base2[i - 1] = tempset;
804  }
805  r->state = '?';
806  if (r->base2[dispchar - 1] == 1L << ((long)A))
807    r->state = 'A';
808  if (r->base2[dispchar - 1] == 1L << ((long)C))
809    r->state = 'C';
810  if (r->base2[dispchar - 1] == 1L << ((long)G))
811    r->state = 'G';
812  if (r->base2[dispchar - 1] == 1L << ((long)T))
813    r->state = 'T';
814  *local_bottom = false;
815  if (!r->tip) {
816    memcpy(temparray->base2, r->next->back->base2, chars*sizeof(long));
817    dnamove_ancestset(hypset[i - 1], r->next->next->back->base2[i - 1],
818              &ancset->base2[i - 1]);
819    dnamove_hyptrav(r->next->back, ancset->base2, i,local_bottom);
820    dnamove_ancestset(hypset[i - 1], temparray->base2[i - 1],
821              &ancset->base2[i - 1]);
822    dnamove_hyptrav(r->next->next->back, ancset->base2, i,local_bottom);
823  }
824  dnamove_chuck(temparray);
825  dnamove_chuck(ancset);
826}  /* dnamove_hyptrav */
827
828void dnamove_hypstates()
829{
830  /* fill in and describe states at interior nodes */
831  long i;
832  boolean local_bottom;
833  baseptr2 nothing;
834
835  i = dispchar;
836  nothing = (baseptr2)Malloc(chars*sizeof(long));
837  nothing[i - 1] = 0;
838  local_bottom = true;
839  firstrav(root, i);
840  dnamove_hyptrav(root, nothing, i,&local_bottom);
841  free(nothing);
842}  /* dnamove_hypstates */
843
844
845void grwrite(chartype c, long num, long *pos)
846{
847  long i;
848
849  prefix(c);
850  for (i = 1; i <= num; i++) {
851    if ((*pos) >= leftedge && (*pos) - leftedge + 1 < screenwidth)
852      putchar(chh[(long)c]);
853    (*pos)++;
854  }
855  postfix(c);
856}  /* grwrite */
857
858
859void dnamove_drawline(long i)
860{
861  /* draws one row of the tree diagram by moving up tree */
862  node *p, *q, *r, *first =NULL, *last =NULL;
863  long n, j, pos;
864  boolean extra, done;
865  Char st;
866  chartype c, d;
867
868  pos = 1;
869  p = nuroot;
870  q = nuroot;
871  extra = false;
872  if (i == p->ycoord && (p == root || subtree)) {
873    extra = true;
874    c = overt;
875    if (display) {
876      switch (p->state) {
877       
878      case 'A':
879        c = aa;
880        break;
881       
882      case 'C':
883        c = cc;
884        break;
885       
886      case 'G':
887        c = gg;
888        break;
889       
890      case 'T':
891        c = tt;
892        break;
893       
894      case '?':
895        c = question;
896        break;
897      }
898    }
899    if ((subtree))
900      stwrite("Subtree:", 8, &pos, leftedge, screenwidth);
901    if (p->index >= 100)
902      nnwrite(p->index, 3, &pos, leftedge, screenwidth);
903    else if (p->index >= 10) {
904      grwrite(c, 1, &pos);
905      nnwrite(p->index, 2, &pos, leftedge, screenwidth);
906    } else {
907      grwrite(c, 2, &pos);
908      nnwrite(p->index, 1, &pos, leftedge, screenwidth);
909    }
910  } else {
911    if ((subtree))
912      stwrite("          ", 10, &pos, leftedge, screenwidth);
913    else
914      stwrite("  ", 2, &pos, leftedge, screenwidth);
915  }
916  do {
917    if (!p->tip) {
918      r = p->next;
919      done = false;
920      do {
921        if (i >= r->back->ymin && i <= r->back->ymax) {
922          q = r->back;
923          done = true;
924        }
925        r = r->next;
926      } while (!(done || r == p));
927      first = p->next->back;
928      r = p->next;
929      while (r->next != p)
930        r = r->next;
931      last = r->back;
932    }
933    done = (p == q);
934    n = p->xcoord - q->xcoord;
935    if (n < 3 && !q->tip)
936      n = 3;
937    if (extra) {
938      n--;
939      extra = false;
940    }
941    if (q->ycoord == i && !done) {
942      if (q->ycoord > p->ycoord)
943        d = upcorner;
944      else
945        d = downcorner;
946      c = overt;
947      if (display) {
948        switch (q->state) {
949       
950        case 'A':
951          c = aa;
952          break;
953       
954        case 'C':
955          c = cc;
956          break;
957       
958        case 'G':
959          c = gg;
960          break;
961       
962        case 'T':
963          c = tt;
964          break;
965       
966        case '?':
967          c = question;
968          break;
969        }
970        d = c;
971      }
972      if (n > 1) {
973        grwrite(d, 1, &pos);
974        grwrite(c, n - 3, &pos);
975      }
976      if (q->index >= 100)
977        nnwrite(q->index, 3, &pos, leftedge, screenwidth);
978      else if (q->index >= 10) {
979        grwrite(c, 1, &pos);
980        nnwrite(q->index, 2, &pos, leftedge, screenwidth);
981      } else {
982        grwrite(c, 2, &pos);
983        nnwrite(q->index, 1, &pos, leftedge, screenwidth);
984      }
985      extra = true;
986    } else if (!q->tip) {
987      if (last->ycoord > i && first->ycoord < i && i != p->ycoord) {
988        c = up;
989        if (i < p->ycoord)
990          st = p->next->back->state;
991        else
992          st = p->next->next->back->state;
993        if (display) {
994          switch (st) {
995       
996          case 'A':
997            c = aa;
998            break;
999       
1000          case 'C':
1001            c = cc;
1002            break;
1003       
1004          case 'G':
1005            c = gg;
1006            break;
1007       
1008          case 'T':
1009            c = tt;
1010            break;
1011       
1012          case '?':
1013            c = question;
1014            break;
1015          }
1016        }
1017        grwrite(c, 1, &pos);
1018        chwrite(' ', n - 1, &pos, leftedge, screenwidth);
1019      } else
1020        chwrite(' ', n, &pos, leftedge, screenwidth);
1021    } else
1022      chwrite(' ', n, &pos, leftedge, screenwidth);
1023    if (p != q)
1024      p = q;
1025  } while (!done);
1026  if (p->ycoord == i && p->tip) {
1027    n = 0;
1028    for (j = 1; j <= nmlngth; j++) {
1029      if (nayme[p->index - 1][j - 1] != '\0')
1030        n = j;
1031    }
1032    chwrite(':', 1, &pos, leftedge, screenwidth);
1033    for (j = 0; j < n; j++)
1034      chwrite(nayme[p->index - 1][j], 1, &pos, leftedge, screenwidth);
1035  }
1036  putchar('\n');
1037}  /* dnamove_drawline */
1038
1039void dnamove_printree()
1040{
1041  /* prints out diagram of the tree */
1042  long tipy;
1043  long i, dow;
1044
1045  if (!subtree)
1046    nuroot = root;
1047  if (changed || newtree)
1048    evaluate(root);
1049  if (display)
1050    dnamove_hypstates();
1051  printf((ansi || ibmpc) ? "\033[2J\033[H" : "\n");
1052  tipy = 1;
1053  dow = down;
1054  if (spp * dow > screenlines && !subtree)
1055    dow--;
1056
1057  printf("(unrooted)");
1058  if (display) {
1059    printf(" ");
1060    makechar(aa);
1061    printf(":A ");
1062    makechar(cc);
1063    printf(":C ");
1064    makechar(gg);
1065    printf(":G ");
1066    makechar(tt);
1067    printf(":T ");
1068    makechar(question);
1069    printf(":?");
1070  } else
1071    printf("                    ");
1072  if (!earlytree) {
1073    printf("%10.1f Steps", -like);
1074  }
1075  if (display)
1076    printf("  SITE%4ld", dispchar);
1077  else
1078    printf("         ");
1079  if (!earlytree) {
1080    printf("  %3ld sites compatible\n", compatible);
1081  }
1082
1083  printf("                            ");
1084  if (changed && !earlytree) {
1085    if (-like < bestyet) {
1086      printf("     BEST YET!");
1087      bestyet = -like;
1088    } else if (fabs(-like - bestyet) < 0.000001)
1089      printf("     (as good as best)");
1090    else {
1091      if (-like < gotlike)
1092        printf("     better");
1093      else if (-like > gotlike)
1094        printf("     worse!");
1095    }
1096  }
1097  printf("\n");
1098
1099  farthest = 0;
1100  coordinates(nuroot, &tipy, 1.5, &farthest);
1101  vmargin = 4;
1102  treelines = tipy - dow;
1103  if (topedge != 1) {
1104    printf("** %ld lines above screen **\n", topedge - 1);
1105    vmargin++;
1106  }
1107  if ((treelines - topedge + 1) > (screenlines - vmargin))
1108    vmargin++;
1109  for (i = 1; i <= treelines; i++) {
1110    if (i >= topedge && i < topedge + screenlines - vmargin)
1111      dnamove_drawline(i);
1112  }
1113  if ((treelines - topedge + 1) > (screenlines - vmargin)) {
1114    printf("** %ld", treelines - (topedge - 1 + screenlines - vmargin));
1115    printf(" lines below screen **\n");
1116  }
1117  if (treelines - topedge + vmargin + 1 < screenlines)
1118    putchar('\n');
1119  gotlike = -like;
1120  changed = false;
1121}  /* dnamove_printree */
1122
1123
1124void arbitree()
1125{
1126  long i;
1127
1128  root = treenode[0];
1129  dnamove_add(treenode[0], treenode[1], treenode[spp]);
1130  for (i = 3; i <= (spp); i++)
1131    dnamove_add(treenode[spp + i - 3], treenode[i - 1], treenode[spp + i - 2]);
1132  for (i = 0; i < (nonodes); i++)
1133    in_tree[i] = true;
1134}  /* arbitree */
1135
1136
1137void yourtree()
1138{
1139  long i, j;
1140  boolean ok;
1141
1142  root = treenode[0];
1143  dnamove_add(treenode[0], treenode[1], treenode[spp]);
1144  i = 2;
1145  do {
1146    i++;
1147    dnamove_printree();
1148    printf("Add species%3ld: ", i);
1149    for (j = 0; j < nmlngth; j++)
1150      putchar(nayme[i - 1][j]);
1151    do {
1152      printf("\n before node (type number): ");
1153      inpnum(&j, &ok);
1154      ok = (ok && ((j >= 1 && j < i) || (j > spp && j < spp + i - 1)));
1155      if (!ok)
1156        printf("Impossible number. Please try again:\n");
1157    } while (!ok);
1158    dnamove_add(treenode[j - 1], treenode[i - 1], treenode[spp + i - 2]);
1159  } while (i != spp);
1160  for (i = 0; i < (nonodes); i++)
1161    in_tree[i] = true;
1162}  /* yourtree */
1163
1164
1165void initdnamovenode(node **p, node **local_grbg, node *UNUSED_q, long UNUSED_len, long nodei,
1166                        long *UNUSED_ntips, long *parens, initops whichinit,
1167                        pointarray local_treenode, pointarray UNUSED_nodep, Char *str, Char *ch,
1168                        FILE *fp_intree)
1169{
1170    (void)UNUSED_q;
1171    (void)UNUSED_len;
1172    (void)UNUSED_ntips;
1173    (void)UNUSED_nodep;
1174
1175  /* initializes a node */
1176  /* LM 7/27  I added this function and the commented lines around */
1177  /* treeread() to get the program running, but all 4 move programs*/
1178  /* are improperly integrated into the v4.0 support files.  As is */
1179  /* endsite = chars and this is a patchwork function                   */
1180  boolean minusread;
1181  double valyew, divisor;
1182
1183  switch (whichinit) {
1184  case bottom:
1185    gnutreenode(local_grbg, p, nodei, endsite, zeros);
1186    local_treenode[nodei - 1] = *p;
1187    break;
1188  case nonbottom:
1189    gnutreenode(local_grbg, p, nodei, endsite, zeros);
1190    break;
1191  case tip:
1192    match_names_to_data (str, local_treenode, p, spp);
1193    break;
1194  case length:
1195    processlength(&valyew, &divisor, ch, &minusread, fp_intree, parens);
1196    /* process lengths and discard */
1197  default:      /*cases hslength,hsnolength,treewt,unittrwt,iter,*/
1198    break;
1199  }
1200} /* initdnamovenode */
1201
1202
1203void buildtree()
1204{
1205  long i, nextnode;
1206  node *p;
1207  long j;
1208
1209  changed = false;
1210  newtree = false;
1211  switch (how) {
1212
1213  case arb:
1214    arbitree();
1215    break;
1216
1217  case use:
1218    openfile(&intree,intreename,"input tree file", "r",progname,intreename);
1219    names = (boolean *)Malloc(spp*sizeof(boolean));
1220    firsttree = true;                                        /**/
1221    nodep = NULL;                                        /**/
1222    nextnode = 0;                                        /**/
1223    haslengths = 0;                                        /**/
1224    endsite = chars;                        /*debug*/
1225    zeros = (long *)Malloc(endsite*sizeof(long));        /**/
1226    for (i = 0; i < endsite; i++)                        /**/
1227      zeros[i] = 0;                                        /**/
1228    treeread(intree, &root, treenode, &goteof, &firsttree,
1229                nodep, &nextnode, &haslengths,
1230                &grbg, initdnamovenode); /*debug*/
1231    for (i = spp; i < (nonodes); i++) {
1232      p = treenode[i];
1233      for (j = 1; j <= 3; j++) {
1234        p->base2 = (baseptr2)Malloc(chars*sizeof(long));
1235        p = p->next;
1236      } 
1237    } /* debug: see comment at initdnamovenode() */
1238
1239    for (i = 0; i < (spp); i++)
1240      in_tree[i] = names[i];
1241    free(names);
1242    FClose(intree);
1243    break;
1244
1245  case spec:
1246    yourtree();
1247    break;
1248  }
1249  if (!outgropt)
1250    outgrno = root->next->back->index;
1251  if (outgropt && in_tree[outgrno - 1])
1252    dnamove_reroot(treenode[outgrno - 1]);
1253}  /* buildtree */
1254
1255
1256void setorder()
1257{
1258  /* sets in order of number of members */
1259  sett[0] = 1L << ((long)A);
1260  sett[1] = 1L << ((long)C);
1261  sett[2] = 1L << ((long)G);
1262  sett[3] = 1L << ((long)T);
1263  sett[4] = 1L << ((long)O);
1264  sett[5] = (1L << ((long)A)) | (1L << ((long)C));
1265  sett[6] = (1L << ((long)A)) | (1L << ((long)G));
1266  sett[7] = (1L << ((long)A)) | (1L << ((long)T));
1267  sett[8] = (1L << ((long)A)) | (1L << ((long)O));
1268  sett[9] = (1L << ((long)C)) | (1L << ((long)G));
1269  sett[10] = (1L << ((long)C)) | (1L << ((long)T));
1270  sett[11] = (1L << ((long)C)) | (1L << ((long)O));
1271  sett[12] = (1L << ((long)G)) | (1L << ((long)T));
1272  sett[13] = (1L << ((long)G)) | (1L << ((long)O));
1273  sett[14] = (1L << ((long)T)) | (1L << ((long)O));
1274  sett[15] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)G));
1275  sett[16] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)T));
1276  sett[17] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)O));
1277  sett[18] = (1L << ((long)A)) | (1L << ((long)G)) | (1L << ((long)T));
1278  sett[19] = (1L << ((long)A)) | (1L << ((long)G)) | (1L << ((long)O));
1279  sett[20] = (1L << ((long)A)) | (1L << ((long)T)) | (1L << ((long)O));
1280  sett[21] = (1L << ((long)C)) | (1L << ((long)G)) | (1L << ((long)T));
1281  sett[22] = (1L << ((long)C)) | (1L << ((long)G)) | (1L << ((long)O));
1282  sett[23] = (1L << ((long)C)) | (1L << ((long)T)) | (1L << ((long)O));
1283  sett[24] = (1L << ((long)G)) | (1L << ((long)T)) | (1L << ((long)O));
1284  sett[25] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)G)) |
1285    (1L << ((long)T));
1286  sett[26] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)G)) |
1287    (1L << ((long)O));
1288  sett[27] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)T)) |
1289    (1L << ((long)O));
1290  sett[28] = (1L << ((long)A)) | (1L << ((long)G)) | (1L << ((long)T)) |
1291    (1L << ((long)O));
1292  sett[29] = (1L << ((long)C)) | (1L << ((long)G)) | (1L << ((long)T)) |
1293    (1L << ((long)O));
1294  sett[30] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)G)) |
1295    (1L << ((long)T)) | (1L << ((long)O));
1296}  /* setorder */
1297
1298
1299void mincomp()
1300{
1301  /* computes for each site the minimum number of steps
1302     necessary to accomodate those species already
1303     in the analysis */
1304  long i, j, k;
1305  boolean done;
1306
1307  for (i = 0; i < (chars); i++) {
1308    done = false;
1309    j = 0;
1310    while (!done) {
1311      j++;
1312      done = true;
1313      k = 1;
1314      do {
1315        if (in_tree[k - 1])
1316          done = (done && (treenode[k - 1]->base2[i] & sett[j - 1]) != 0);
1317        k++;
1318      } while (k <= spp && done);
1319    }
1320    if (j == 31)
1321      necsteps[i] = 4;
1322    if (j <= 30)
1323      necsteps[i] = 3;
1324    if (j <= 25)
1325      necsteps[i] = 2;
1326    if (j <= 15)
1327      necsteps[i] = 1;
1328    if (j <= 5)
1329      necsteps[i] = 0;
1330    necsteps[i] *= weight[i];
1331  }
1332}  /* mincomp */
1333
1334
1335void rearrange()
1336{
1337  long i, j;
1338  boolean ok1, ok2;
1339  node *p, *q;
1340
1341  printf("Remove everything to the right of which node? ");
1342  inpnum(&i, &ok1);
1343  ok1 = (ok1 && i >= 1 && i < spp * 2 && i != root->index);
1344  if (ok1) {
1345    printf("Add before which node? ");
1346    inpnum(&j, &ok2);
1347    ok2 = (ok2 && j >= 1 && j < spp * 2);
1348    if (ok2) {
1349      ok2 = (treenode[j - 1] != treenode[treenode[i - 1]->back->index - 1]);
1350      p = treenode[j - 1];
1351      while (p != root) {
1352        ok2 = (ok2 && p != treenode[i - 1]);
1353        p = treenode[p->back->index - 1];
1354      }
1355      if (ok1 && ok2) {
1356        what = i;
1357        q = treenode[treenode[i - 1]->back->index - 1];
1358        if (q->next->back->index == i)
1359          fromwhere = q->next->next->back->index;
1360        else
1361          fromwhere = q->next->back->index;
1362        towhere = j;
1363        dnamove_re_move(&treenode[i - 1], &q);
1364        dnamove_add(treenode[j - 1], treenode[i - 1], q);
1365      }
1366      lastop = rearr;
1367    }
1368  }
1369  changed = (ok1 && ok2);
1370  dnamove_printree();
1371  if (!(ok1 && ok2))
1372    printf("Not a possible rearrangement.  Try again: \n");
1373  else {
1374    oldwritten = written;
1375    written = false;
1376  }
1377}  /* rearrange */
1378
1379void dnamove_nextinc()
1380{
1381  /* show next incompatible site */
1382  long disp0;
1383  boolean done;
1384
1385  display = true;
1386  disp0 = dispchar;
1387  done = false;
1388  do {
1389    dispchar++;
1390    if (dispchar > chars) {
1391      dispchar = 1;
1392      done = (disp0 == 0);
1393    }
1394  } while (!(necsteps[dispchar - 1] != numsteps[dispchar - 1] ||
1395             dispchar == disp0 || done));
1396  dnamove_printree();
1397}  /* dnamove_nextinc */
1398
1399void dnamove_nextchar()
1400{
1401  /* show next site */
1402  display = true;
1403  dispchar++;
1404  if (dispchar > chars)
1405    dispchar = 1;
1406  dnamove_printree();
1407}  /* dnamove_nextchar */
1408
1409void dnamove_prevchar()
1410{
1411  /* show previous site */
1412  display = true;
1413  dispchar--;
1414  if (dispchar < 1)
1415    dispchar = chars;
1416  dnamove_printree();
1417}  /* dnamove_prevchar */
1418
1419void dnamove_show()
1420{
1421  long i;
1422  boolean ok;
1423
1424  do {
1425    printf("SHOW: (Character number or 0 to see none)? ");
1426    inpnum(&i, &ok);
1427    ok = (ok && (i == 0 || (i >= 1 && i <= chars)));
1428    if (ok && i != 0) {
1429      display = true;
1430      dispchar = i;
1431    }
1432    if (ok && i == 0)
1433      display = false;
1434  } while (!ok);
1435  dnamove_printree();
1436}  /* dnamove_show */
1437
1438
1439void tryadd(node *p, node **item, node **nufork, double *place)
1440{
1441  /* temporarily adds one fork and one tip to the tree.
1442     Records scores in ARRAY place */
1443  dnamove_add(p, *item, *nufork);
1444  evaluate(root);
1445  place[p->index - 1] = -like;
1446  dnamove_re_move(item, nufork);
1447}  /* tryadd */
1448
1449
1450void addpreorder(node *p, node *item_, node *nufork_, double *place)
1451{
1452  /* traverses a binary tree, calling PROCEDURE tryadd
1453     at a node before calling tryadd at its descendants */
1454  node *item, *nufork;
1455
1456  item = item_;
1457  nufork = nufork_;
1458  if (p == NULL)
1459    return;
1460  tryadd(p,&item,&nufork,place);
1461  if (!p->tip) {
1462    addpreorder(p->next->back, item,nufork,place);
1463    addpreorder(p->next->next->back,item,nufork,place);
1464  }
1465}  /* addpreorder */
1466
1467
1468void try()
1469{
1470  /* Remove node, try it in all possible places */
1471  double *place;
1472  long i, j, oldcompat;
1473  double current;
1474  node *q, *dummy, *rute;
1475  boolean tied, better, ok;
1476
1477  printf("Try other positions for which node? ");
1478  inpnum(&i, &ok);
1479  if (!(ok && i >= 1 && i <= nonodes && i != root->index)) {
1480    printf("Not a possible choice! ");
1481    return;
1482  }
1483  printf("WAIT ...\n");
1484  place = (double *)Malloc(nonodes*sizeof(double));
1485  for (j = 0; j < (nonodes); j++)
1486    place[j] = -1.0;
1487  evaluate(root);
1488  current = -like;
1489  oldcompat = compatible;
1490  what = i;
1491  q = treenode[treenode[i - 1]->back->index - 1];
1492  if (q->next->back->index == i)
1493    fromwhere = q->next->next->back->index;
1494  else
1495    fromwhere = q->next->back->index;
1496  rute = root;
1497  if (root == treenode[treenode[i - 1]->back->index - 1]) {
1498    if (treenode[treenode[i - 1]->back->index - 1]->next->back == treenode[i - 1])
1499      rute = treenode[treenode[i - 1]->back->index - 1]->next->next->back;
1500    else
1501      rute = treenode[treenode[i - 1]->back->index - 1]->next->back;
1502  }
1503  dnamove_re_move(&treenode[i - 1], &dummy);
1504  oldleft = wasleft;
1505  root = rute;
1506  addpreorder(root, treenode[i - 1], dummy, place);
1507  wasleft = oldleft;
1508  restoring = true;
1509  dnamove_add(treenode[fromwhere - 1], treenode[what - 1], q);
1510  like = -current;
1511  compatible = oldcompat;
1512  restoring = false;
1513  better = false;
1514  printf("       BETTER: ");
1515  for (j = 1; j <= (nonodes); j++) {
1516    if (place[j - 1] < current && place[j - 1] >= 0.0) {
1517      printf("%3ld:%6.2f", j, place[j - 1]);
1518      better = true;
1519    }
1520  }
1521  if (!better)
1522    printf(" NONE");
1523  printf("\n       TIED:    ");
1524  tied = false;
1525  for (j = 1; j <= (nonodes); j++) {
1526    if (fabs(place[j - 1] - current) < 1.0e-6 && j != fromwhere) {
1527      if (j < 10)
1528        printf("%2ld", j);
1529      else
1530        printf("%3ld", j);
1531      tied = true;
1532    }
1533  }
1534  if (tied)
1535    printf(":%6.2f\n", current);
1536  else
1537    printf("NONE\n");
1538  changed = true;
1539  free(place);
1540}  /* try */
1541
1542
1543void undo()
1544{
1545  /* restore to tree before last rearrangement */
1546  long temp;
1547  boolean btemp;
1548  node *q;
1549
1550  switch (lastop) {
1551
1552  case rearr:
1553    restoring = true;
1554    oldleft = wasleft;
1555    dnamove_re_move(&treenode[what - 1], &q);
1556    btemp = wasleft;
1557    wasleft = oldleft;
1558    dnamove_add(treenode[fromwhere - 1], treenode[what - 1],q);
1559    wasleft = btemp;
1560    restoring = false;
1561    temp = fromwhere;
1562    fromwhere = towhere;
1563    towhere = temp;
1564    changed = true;
1565    break;
1566
1567  case flipp:
1568    q = treenode[atwhat - 1]->next->back;
1569    treenode[atwhat - 1]->next->back =treenode[atwhat - 1]->next->next->back;
1570    treenode[atwhat - 1]->next->next->back = q;
1571    treenode[atwhat - 1]->next->back->back = treenode[atwhat - 1]->next;
1572    treenode[atwhat - 1]->next->next->back->back =
1573      treenode[atwhat - 1]->next->next;
1574    break;
1575
1576  case reroott:
1577    restoring = true;
1578    temp =oldoutgrno;
1579    oldoutgrno = outgrno;
1580    outgrno = temp;
1581    dnamove_reroot(treenode[outgrno - 1]);
1582    restoring = false;
1583    break;
1584
1585  case none:
1586    /* blank case */
1587    break;
1588  }
1589  dnamove_printree();
1590  if (lastop == none) {
1591    printf("No operation to undo! ");
1592    return;
1593  }
1594  btemp = oldwritten;
1595  oldwritten = written;
1596  written = btemp;
1597}  /* undo */
1598
1599
1600void treewrite(boolean done)
1601{
1602  /* write out tree to a file */
1603  Char ch;
1604
1605  treeoptions(waswritten, &ch, &outtree, outtreename, progname);
1606  if (!done)
1607    dnamove_printree();
1608  if (waswritten && ch != 'A' && ch != 'R')
1609    return;
1610  col = 0;
1611  treeout(root, 1, &col, root);
1612  printf("\nTree written to file \"%s\"\n\n", outtreename);
1613  waswritten = true;
1614  written = true;
1615  FClose(outtree);
1616#ifdef MAC
1617  fixmacfile(outtreename);
1618#endif
1619}
1620/* treewrite */
1621
1622
1623void clade()
1624{
1625  /* pick a subtree and show only that on screen */
1626  long i;
1627  boolean ok;
1628
1629  printf("Select subtree rooted at which node (0 for whole tree)? ");
1630  inpnum(&i, &ok);
1631  ok = (ok && ((unsigned)i) <= ((unsigned)nonodes));
1632  if (ok) {
1633    subtree = (i > 0);
1634    if (subtree)
1635      nuroot = treenode[i - 1];
1636    else
1637      nuroot = root;
1638  }
1639  dnamove_printree();
1640  if (!ok)
1641    printf("Not possible to use this node. ");
1642}  /* clade */
1643
1644
1645void flip()
1646{
1647  /* flip at a node left-right */
1648  long i;
1649  boolean ok;
1650  node *p;
1651
1652  printf("Flip branches at which node? ");
1653  inpnum(&i, &ok);
1654  ok = (ok && i > spp && i <= nonodes);
1655  if (ok) {
1656    p = treenode[i - 1]->next->back;
1657    treenode[i - 1]->next->back = treenode[i - 1]->next->next->back;
1658    treenode[i - 1]->next->next->back = p;
1659    treenode[i - 1]->next->back->back = treenode[i - 1]->next;
1660    treenode[i - 1]->next->next->back->back = treenode[i - 1]->next->next;
1661    atwhat = i;
1662    lastop = flipp;
1663  }
1664  dnamove_printree();
1665  if (ok) {
1666    oldwritten = written;
1667    written = false;
1668    return;
1669  }
1670  if (i >= 1 && i <= spp)
1671    printf("Can't flip there. ");
1672  else
1673    printf("No such node. ");
1674}  /* flip */
1675
1676void changeoutgroup()
1677{
1678  long i;
1679  boolean ok;
1680
1681  oldoutgrno = outgrno;
1682  do {
1683    printf("Which node should be the new outgroup? ");
1684    inpnum(&i, &ok);
1685    ok = (ok && in_tree[i - 1] && i >= 1 && i <= nonodes &&
1686          i != root->index);
1687    if (ok)
1688      outgrno = i;
1689  } while (!ok);
1690  if (in_tree[outgrno - 1])
1691    dnamove_reroot(treenode[outgrno - 1]);
1692  changed = true;
1693  lastop = reroott;
1694  dnamove_printree();
1695  oldwritten = written;
1696  written = false;
1697}  /* changeoutgroup */
1698
1699
1700void redisplay()
1701{
1702  boolean done = false;
1703  waswritten = false;
1704  do {
1705    printf("NEXT? (Options: R # + - S . T U W O F H J K L C ? X Q) ");
1706    printf("(? for Help) ");
1707#ifdef WIN32
1708    phyFillScreenColor();
1709#endif
1710    scanf("%c%*[^\n]", &global_ch);
1711    getchar();
1712    uppercase(&global_ch);
1713    if (strchr("HJKLCFORSTUXQ+#-.W?",global_ch) != NULL){
1714      switch (global_ch) {
1715       
1716      case 'R':
1717        rearrange();
1718        break;
1719       
1720      case '#':
1721        dnamove_nextinc();
1722        break;
1723       
1724      case '+':
1725        dnamove_nextchar();
1726        break;
1727       
1728      case '-':
1729        dnamove_prevchar();
1730        break;
1731       
1732      case 'S':
1733        dnamove_show();
1734        break;
1735       
1736      case '.':
1737        dnamove_printree();
1738        break;
1739       
1740      case 'T':
1741        try();
1742        break;
1743       
1744      case 'U':
1745        undo();
1746        break;
1747       
1748      case 'W':
1749        treewrite(done);
1750        break;
1751       
1752      case 'O':
1753        changeoutgroup();
1754        break;
1755       
1756      case 'F':
1757        flip();
1758        break;
1759       
1760      case 'H':
1761        window(left, &leftedge, &topedge, hscroll, vscroll, treelines,
1762                 screenlines, screenwidth, farthest, subtree);
1763        dnamove_printree();
1764        break;
1765
1766      case 'J':
1767        window(downn, &leftedge, &topedge, hscroll, vscroll, treelines,
1768                 screenlines, screenwidth, farthest, subtree);
1769        dnamove_printree();
1770        break;
1771
1772      case 'K':
1773        window(upp, &leftedge, &topedge, hscroll, vscroll, treelines,
1774                 screenlines, screenwidth, farthest, subtree);
1775        dnamove_printree();
1776        break;
1777
1778      case 'L':
1779        window(right, &leftedge, &topedge, hscroll, vscroll, treelines,
1780                 screenlines, screenwidth, farthest, subtree);
1781        dnamove_printree();
1782        break;
1783
1784      case 'C':
1785        clade();
1786        break;
1787       
1788      case '?':
1789        help("site");
1790        dnamove_printree();
1791        break;
1792       
1793      case 'X':
1794        done = true;
1795        break;
1796       
1797      case 'Q':
1798        done = true;
1799        break;
1800      }
1801    }
1802  } while (!done);
1803  if (written)
1804    return;
1805  do {
1806    printf("Do you want to write out the tree to a file? (Y or N) ");
1807#ifdef WIN32
1808    phyFillScreenColor();
1809#endif
1810    scanf("%c%*[^\n]", &global_ch);
1811    getchar();
1812    if (global_ch == 'Y' || global_ch == 'y')
1813      treewrite(done);
1814  } while (global_ch != 'Y' && global_ch != 'y' && global_ch != 'N' && global_ch != 'n');
1815}  /* redisplay */
1816
1817
1818void treeconstruct()
1819{
1820  /* constructs a binary tree from the pointers in treenode. */
1821
1822  restoring = false;
1823  subtree = false;
1824  display = false;
1825  dispchar = 0;
1826  earlytree = true;
1827  waswritten = false;
1828  buildtree();
1829  printf("\nComputing steps needed for compatibility in sites ...\n\n");
1830  setorder();
1831  mincomp();
1832  newtree = true;
1833  earlytree = false;
1834  dnamove_printree();
1835  bestyet = -like;
1836  gotlike = -like;
1837  lastop = none;
1838  newtree = false;
1839  written = false;
1840  redisplay();
1841}  /* treeconstruct */
1842
1843
1844int main(int argc, Char *argv[])
1845{ /* Interactive DNA parsimony */
1846  /* reads in spp, chars, and the data. Then calls treeconstruct to
1847     construct the tree and query the user */
1848#ifdef MAC
1849  argc = 1;                /* macsetup("Dnamove","");        */
1850  argv[0] = "Dnamove";
1851#endif
1852  init(argc, argv);
1853  progname = argv[0];
1854  strcpy(infilename,INFILE);
1855  strcpy(intreename,INTREE);
1856  strcpy(outtreename,OUTTREE);
1857
1858  openfile(&infile,infilename,"input file", "r",argv[0],infilename);
1859  openfile(&outtree,outtreename,"output file", "w",argv[0],outtreename);
1860
1861  garbage = NULL;
1862  screenlines = 24;
1863  scrollinc = 20;
1864  screenwidth = 80;
1865  topedge = 1;
1866  leftedge = 1;
1867  ibmpc = IBMCRT;
1868  ansi = ANSICRT;
1869  doinput();
1870  configure();
1871  treeconstruct();
1872  FClose(infile);
1873  FClose(outtree);
1874#ifdef MAC
1875  fixmacfile(outtreename);
1876#endif
1877#ifdef WIN32
1878  phyRestoreConsoleAttributes();
1879#endif
1880  return 0;
1881}  /* Interactive DNA parsimony */
Note: See TracBrowser for help on using the repository browser.