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

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

upgrade to PHYLIP 3.6
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 45.9 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  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 (ch == '\n')
309            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 *bottom)
786{
787  /*  compute, print out state at one interior node */
788  long tempset, 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    left = r->next->back->base2[i - 1];
796    rt = r->next->next->back->base2[i - 1];
797    tempset = left & rt & anc;
798    if (tempset == 0) {
799      tempset = (left & rt) | (left & anc) | (rt & anc);
800      if (tempset == 0)
801        tempset = 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  *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,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,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 bottom;
833  baseptr2 nothing;
834
835  i = dispchar;
836  nothing = (baseptr2)Malloc(chars*sizeof(long));
837  nothing[i - 1] = 0;
838  bottom = true;
839  firstrav(root, i);
840  dnamove_hyptrav(root, nothing, i,&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 **grbg, node *q, long len, long nodei,
1166                        long *ntips, long *parens, initops whichinit,
1167                        pointarray treenode, pointarray nodep, Char *str, Char *ch,
1168                        FILE *intree)
1169{
1170  /* initializes a node */
1171  /* LM 7/27  I added this function and the commented lines around */
1172  /* treeread() to get the program running, but all 4 move programs*/
1173  /* are improperly integrated into the v4.0 support files.  As is */
1174  /* endsite = chars and this is a patchwork function                   */
1175  boolean minusread;
1176  double valyew, divisor;
1177
1178  switch (whichinit) {
1179  case bottom:
1180    gnutreenode(grbg, p, nodei, endsite, zeros);
1181    treenode[nodei - 1] = *p;
1182    break;
1183  case nonbottom:
1184    gnutreenode(grbg, p, nodei, endsite, zeros);
1185    break;
1186  case tip:
1187    match_names_to_data (str, treenode, p, spp);
1188    break;
1189  case length:
1190    processlength(&valyew, &divisor, ch, &minusread, intree, parens);
1191    /* process lengths and discard */
1192  default:      /*cases hslength,hsnolength,treewt,unittrwt,iter,*/
1193    break;
1194  }
1195} /* initdnamovenode */
1196
1197
1198void buildtree()
1199{
1200  long i, nextnode;
1201  node *p;
1202  long j;
1203
1204  changed = false;
1205  newtree = false;
1206  switch (how) {
1207
1208  case arb:
1209    arbitree();
1210    break;
1211
1212  case use:
1213    openfile(&intree,intreename,"input tree file", "r",progname,intreename);
1214    names = (boolean *)Malloc(spp*sizeof(boolean));
1215    firsttree = true;                                        /**/
1216    nodep = NULL;                                        /**/
1217    nextnode = 0;                                        /**/
1218    haslengths = 0;                                        /**/
1219    endsite = chars;                        /*debug*/
1220    zeros = (long *)Malloc(endsite*sizeof(long));        /**/
1221    for (i = 0; i < endsite; i++)                        /**/
1222      zeros[i] = 0;                                        /**/
1223    treeread(intree, &root, treenode, &goteof, &firsttree,
1224                nodep, &nextnode, &haslengths,
1225                &grbg, initdnamovenode); /*debug*/
1226    for (i = spp; i < (nonodes); i++) {
1227      p = treenode[i];
1228      for (j = 1; j <= 3; j++) {
1229        p->base2 = (baseptr2)Malloc(chars*sizeof(long));
1230        p = p->next;
1231      } 
1232    } /* debug: see comment at initdnamovenode() */
1233
1234    for (i = 0; i < (spp); i++)
1235      in_tree[i] = names[i];
1236    free(names);
1237    FClose(intree);
1238    break;
1239
1240  case spec:
1241    yourtree();
1242    break;
1243  }
1244  if (!outgropt)
1245    outgrno = root->next->back->index;
1246  if (outgropt && in_tree[outgrno - 1])
1247    dnamove_reroot(treenode[outgrno - 1]);
1248}  /* buildtree */
1249
1250
1251void setorder()
1252{
1253  /* sets in order of number of members */
1254  sett[0] = 1L << ((long)A);
1255  sett[1] = 1L << ((long)C);
1256  sett[2] = 1L << ((long)G);
1257  sett[3] = 1L << ((long)T);
1258  sett[4] = 1L << ((long)O);
1259  sett[5] = (1L << ((long)A)) | (1L << ((long)C));
1260  sett[6] = (1L << ((long)A)) | (1L << ((long)G));
1261  sett[7] = (1L << ((long)A)) | (1L << ((long)T));
1262  sett[8] = (1L << ((long)A)) | (1L << ((long)O));
1263  sett[9] = (1L << ((long)C)) | (1L << ((long)G));
1264  sett[10] = (1L << ((long)C)) | (1L << ((long)T));
1265  sett[11] = (1L << ((long)C)) | (1L << ((long)O));
1266  sett[12] = (1L << ((long)G)) | (1L << ((long)T));
1267  sett[13] = (1L << ((long)G)) | (1L << ((long)O));
1268  sett[14] = (1L << ((long)T)) | (1L << ((long)O));
1269  sett[15] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)G));
1270  sett[16] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)T));
1271  sett[17] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)O));
1272  sett[18] = (1L << ((long)A)) | (1L << ((long)G)) | (1L << ((long)T));
1273  sett[19] = (1L << ((long)A)) | (1L << ((long)G)) | (1L << ((long)O));
1274  sett[20] = (1L << ((long)A)) | (1L << ((long)T)) | (1L << ((long)O));
1275  sett[21] = (1L << ((long)C)) | (1L << ((long)G)) | (1L << ((long)T));
1276  sett[22] = (1L << ((long)C)) | (1L << ((long)G)) | (1L << ((long)O));
1277  sett[23] = (1L << ((long)C)) | (1L << ((long)T)) | (1L << ((long)O));
1278  sett[24] = (1L << ((long)G)) | (1L << ((long)T)) | (1L << ((long)O));
1279  sett[25] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)G)) |
1280    (1L << ((long)T));
1281  sett[26] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)G)) |
1282    (1L << ((long)O));
1283  sett[27] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)T)) |
1284    (1L << ((long)O));
1285  sett[28] = (1L << ((long)A)) | (1L << ((long)G)) | (1L << ((long)T)) |
1286    (1L << ((long)O));
1287  sett[29] = (1L << ((long)C)) | (1L << ((long)G)) | (1L << ((long)T)) |
1288    (1L << ((long)O));
1289  sett[30] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)G)) |
1290    (1L << ((long)T)) | (1L << ((long)O));
1291}  /* setorder */
1292
1293
1294void mincomp()
1295{
1296  /* computes for each site the minimum number of steps
1297     necessary to accomodate those species already
1298     in the analysis */
1299  long i, j, k;
1300  boolean done;
1301
1302  for (i = 0; i < (chars); i++) {
1303    done = false;
1304    j = 0;
1305    while (!done) {
1306      j++;
1307      done = true;
1308      k = 1;
1309      do {
1310        if (in_tree[k - 1])
1311          done = (done && (treenode[k - 1]->base2[i] & sett[j - 1]) != 0);
1312        k++;
1313      } while (k <= spp && done);
1314    }
1315    if (j == 31)
1316      necsteps[i] = 4;
1317    if (j <= 30)
1318      necsteps[i] = 3;
1319    if (j <= 25)
1320      necsteps[i] = 2;
1321    if (j <= 15)
1322      necsteps[i] = 1;
1323    if (j <= 5)
1324      necsteps[i] = 0;
1325    necsteps[i] *= weight[i];
1326  }
1327}  /* mincomp */
1328
1329
1330void rearrange()
1331{
1332  long i, j;
1333  boolean ok1, ok2;
1334  node *p, *q;
1335
1336  printf("Remove everything to the right of which node? ");
1337  inpnum(&i, &ok1);
1338  ok1 = (ok1 && i >= 1 && i < spp * 2 && i != root->index);
1339  if (ok1) {
1340    printf("Add before which node? ");
1341    inpnum(&j, &ok2);
1342    ok2 = (ok2 && j >= 1 && j < spp * 2);
1343    if (ok2) {
1344      ok2 = (treenode[j - 1] != treenode[treenode[i - 1]->back->index - 1]);
1345      p = treenode[j - 1];
1346      while (p != root) {
1347        ok2 = (ok2 && p != treenode[i - 1]);
1348        p = treenode[p->back->index - 1];
1349      }
1350      if (ok1 && ok2) {
1351        what = i;
1352        q = treenode[treenode[i - 1]->back->index - 1];
1353        if (q->next->back->index == i)
1354          fromwhere = q->next->next->back->index;
1355        else
1356          fromwhere = q->next->back->index;
1357        towhere = j;
1358        dnamove_re_move(&treenode[i - 1], &q);
1359        dnamove_add(treenode[j - 1], treenode[i - 1], q);
1360      }
1361      lastop = rearr;
1362    }
1363  }
1364  changed = (ok1 && ok2);
1365  dnamove_printree();
1366  if (!(ok1 && ok2))
1367    printf("Not a possible rearrangement.  Try again: \n");
1368  else {
1369    oldwritten = written;
1370    written = false;
1371  }
1372}  /* rearrange */
1373
1374void dnamove_nextinc()
1375{
1376  /* show next incompatible site */
1377  long disp0;
1378  boolean done;
1379
1380  display = true;
1381  disp0 = dispchar;
1382  done = false;
1383  do {
1384    dispchar++;
1385    if (dispchar > chars) {
1386      dispchar = 1;
1387      done = (disp0 == 0);
1388    }
1389  } while (!(necsteps[dispchar - 1] != numsteps[dispchar - 1] ||
1390             dispchar == disp0 || done));
1391  dnamove_printree();
1392}  /* dnamove_nextinc */
1393
1394void dnamove_nextchar()
1395{
1396  /* show next site */
1397  display = true;
1398  dispchar++;
1399  if (dispchar > chars)
1400    dispchar = 1;
1401  dnamove_printree();
1402}  /* dnamove_nextchar */
1403
1404void dnamove_prevchar()
1405{
1406  /* show previous site */
1407  display = true;
1408  dispchar--;
1409  if (dispchar < 1)
1410    dispchar = chars;
1411  dnamove_printree();
1412}  /* dnamove_prevchar */
1413
1414void dnamove_show()
1415{
1416  long i;
1417  boolean ok;
1418
1419  do {
1420    printf("SHOW: (Character number or 0 to see none)? ");
1421    inpnum(&i, &ok);
1422    ok = (ok && (i == 0 || (i >= 1 && i <= chars)));
1423    if (ok && i != 0) {
1424      display = true;
1425      dispchar = i;
1426    }
1427    if (ok && i == 0)
1428      display = false;
1429  } while (!ok);
1430  dnamove_printree();
1431}  /* dnamove_show */
1432
1433
1434void tryadd(node *p, node **item, node **nufork, double *place)
1435{
1436  /* temporarily adds one fork and one tip to the tree.
1437     Records scores in ARRAY place */
1438  dnamove_add(p, *item, *nufork);
1439  evaluate(root);
1440  place[p->index - 1] = -like;
1441  dnamove_re_move(item, nufork);
1442}  /* tryadd */
1443
1444
1445void addpreorder(node *p, node *item_, node *nufork_, double *place)
1446{
1447  /* traverses a binary tree, calling PROCEDURE tryadd
1448     at a node before calling tryadd at its descendants */
1449  node *item, *nufork;
1450
1451  item = item_;
1452  nufork = nufork_;
1453  if (p == NULL)
1454    return;
1455  tryadd(p,&item,&nufork,place);
1456  if (!p->tip) {
1457    addpreorder(p->next->back, item,nufork,place);
1458    addpreorder(p->next->next->back,item,nufork,place);
1459  }
1460}  /* addpreorder */
1461
1462
1463void try()
1464{
1465  /* Remove node, try it in all possible places */
1466  double *place;
1467  long i, j, oldcompat;
1468  double current;
1469  node *q, *dummy, *rute;
1470  boolean tied, better, ok;
1471
1472  printf("Try other positions for which node? ");
1473  inpnum(&i, &ok);
1474  if (!(ok && i >= 1 && i <= nonodes && i != root->index)) {
1475    printf("Not a possible choice! ");
1476    return;
1477  }
1478  printf("WAIT ...\n");
1479  place = (double *)Malloc(nonodes*sizeof(double));
1480  for (j = 0; j < (nonodes); j++)
1481    place[j] = -1.0;
1482  evaluate(root);
1483  current = -like;
1484  oldcompat = compatible;
1485  what = i;
1486  q = treenode[treenode[i - 1]->back->index - 1];
1487  if (q->next->back->index == i)
1488    fromwhere = q->next->next->back->index;
1489  else
1490    fromwhere = q->next->back->index;
1491  rute = root;
1492  if (root == treenode[treenode[i - 1]->back->index - 1]) {
1493    if (treenode[treenode[i - 1]->back->index - 1]->next->back == treenode[i - 1])
1494      rute = treenode[treenode[i - 1]->back->index - 1]->next->next->back;
1495    else
1496      rute = treenode[treenode[i - 1]->back->index - 1]->next->back;
1497  }
1498  dnamove_re_move(&treenode[i - 1], &dummy);
1499  oldleft = wasleft;
1500  root = rute;
1501  addpreorder(root, treenode[i - 1], dummy, place);
1502  wasleft = oldleft;
1503  restoring = true;
1504  dnamove_add(treenode[fromwhere - 1], treenode[what - 1], q);
1505  like = -current;
1506  compatible = oldcompat;
1507  restoring = false;
1508  better = false;
1509  printf("       BETTER: ");
1510  for (j = 1; j <= (nonodes); j++) {
1511    if (place[j - 1] < current && place[j - 1] >= 0.0) {
1512      printf("%3ld:%6.2f", j, place[j - 1]);
1513      better = true;
1514    }
1515  }
1516  if (!better)
1517    printf(" NONE");
1518  printf("\n       TIED:    ");
1519  tied = false;
1520  for (j = 1; j <= (nonodes); j++) {
1521    if (fabs(place[j - 1] - current) < 1.0e-6 && j != fromwhere) {
1522      if (j < 10)
1523        printf("%2ld", j);
1524      else
1525        printf("%3ld", j);
1526      tied = true;
1527    }
1528  }
1529  if (tied)
1530    printf(":%6.2f\n", current);
1531  else
1532    printf("NONE\n");
1533  changed = true;
1534  free(place);
1535}  /* try */
1536
1537
1538void undo()
1539{
1540  /* restore to tree before last rearrangement */
1541  long temp;
1542  boolean btemp;
1543  node *q;
1544
1545  switch (lastop) {
1546
1547  case rearr:
1548    restoring = true;
1549    oldleft = wasleft;
1550    dnamove_re_move(&treenode[what - 1], &q);
1551    btemp = wasleft;
1552    wasleft = oldleft;
1553    dnamove_add(treenode[fromwhere - 1], treenode[what - 1],q);
1554    wasleft = btemp;
1555    restoring = false;
1556    temp = fromwhere;
1557    fromwhere = towhere;
1558    towhere = temp;
1559    changed = true;
1560    break;
1561
1562  case flipp:
1563    q = treenode[atwhat - 1]->next->back;
1564    treenode[atwhat - 1]->next->back =treenode[atwhat - 1]->next->next->back;
1565    treenode[atwhat - 1]->next->next->back = q;
1566    treenode[atwhat - 1]->next->back->back = treenode[atwhat - 1]->next;
1567    treenode[atwhat - 1]->next->next->back->back =
1568      treenode[atwhat - 1]->next->next;
1569    break;
1570
1571  case reroott:
1572    restoring = true;
1573    temp =oldoutgrno;
1574    oldoutgrno = outgrno;
1575    outgrno = temp;
1576    dnamove_reroot(treenode[outgrno - 1]);
1577    restoring = false;
1578    break;
1579
1580  case none:
1581    /* blank case */
1582    break;
1583  }
1584  dnamove_printree();
1585  if (lastop == none) {
1586    printf("No operation to undo! ");
1587    return;
1588  }
1589  btemp = oldwritten;
1590  oldwritten = written;
1591  written = btemp;
1592}  /* undo */
1593
1594
1595void treewrite(boolean done)
1596{
1597  /* write out tree to a file */
1598  Char ch;
1599
1600  treeoptions(waswritten, &ch, &outtree, outtreename, progname);
1601  if (!done)
1602    dnamove_printree();
1603  if (waswritten && ch != 'A' && ch != 'R')
1604    return;
1605  col = 0;
1606  treeout(root, 1, &col, root);
1607  printf("\nTree written to file \"%s\"\n\n", outtreename);
1608  waswritten = true;
1609  written = true;
1610  FClose(outtree);
1611#ifdef MAC
1612  fixmacfile(outtreename);
1613#endif
1614}
1615/* treewrite */
1616
1617
1618void clade()
1619{
1620  /* pick a subtree and show only that on screen */
1621  long i;
1622  boolean ok;
1623
1624  printf("Select subtree rooted at which node (0 for whole tree)? ");
1625  inpnum(&i, &ok);
1626  ok = (ok && ((unsigned)i) <= ((unsigned)nonodes));
1627  if (ok) {
1628    subtree = (i > 0);
1629    if (subtree)
1630      nuroot = treenode[i - 1];
1631    else
1632      nuroot = root;
1633  }
1634  dnamove_printree();
1635  if (!ok)
1636    printf("Not possible to use this node. ");
1637}  /* clade */
1638
1639
1640void flip()
1641{
1642  /* flip at a node left-right */
1643  long i;
1644  boolean ok;
1645  node *p;
1646
1647  printf("Flip branches at which node? ");
1648  inpnum(&i, &ok);
1649  ok = (ok && i > spp && i <= nonodes);
1650  if (ok) {
1651    p = treenode[i - 1]->next->back;
1652    treenode[i - 1]->next->back = treenode[i - 1]->next->next->back;
1653    treenode[i - 1]->next->next->back = p;
1654    treenode[i - 1]->next->back->back = treenode[i - 1]->next;
1655    treenode[i - 1]->next->next->back->back = treenode[i - 1]->next->next;
1656    atwhat = i;
1657    lastop = flipp;
1658  }
1659  dnamove_printree();
1660  if (ok) {
1661    oldwritten = written;
1662    written = false;
1663    return;
1664  }
1665  if (i >= 1 && i <= spp)
1666    printf("Can't flip there. ");
1667  else
1668    printf("No such node. ");
1669}  /* flip */
1670
1671void changeoutgroup()
1672{
1673  long i;
1674  boolean ok;
1675
1676  oldoutgrno = outgrno;
1677  do {
1678    printf("Which node should be the new outgroup? ");
1679    inpnum(&i, &ok);
1680    ok = (ok && in_tree[i - 1] && i >= 1 && i <= nonodes &&
1681          i != root->index);
1682    if (ok)
1683      outgrno = i;
1684  } while (!ok);
1685  if (in_tree[outgrno - 1])
1686    dnamove_reroot(treenode[outgrno - 1]);
1687  changed = true;
1688  lastop = reroott;
1689  dnamove_printree();
1690  oldwritten = written;
1691  written = false;
1692}  /* changeoutgroup */
1693
1694
1695void redisplay()
1696{
1697  boolean done = false;
1698  waswritten = false;
1699  do {
1700    printf("NEXT? (Options: R # + - S . T U W O F H J K L C ? X Q) ");
1701    printf("(? for Help) ");
1702#ifdef WIN32
1703    phyFillScreenColor();
1704#endif
1705    scanf("%c%*[^\n]", &ch);
1706    getchar();
1707    uppercase(&ch); 
1708    if (strchr("HJKLCFORSTUXQ+#-.W?",ch) != NULL){
1709      switch (ch) {
1710       
1711      case 'R':
1712        rearrange();
1713        break;
1714       
1715      case '#':
1716        dnamove_nextinc();
1717        break;
1718       
1719      case '+':
1720        dnamove_nextchar();
1721        break;
1722       
1723      case '-':
1724        dnamove_prevchar();
1725        break;
1726       
1727      case 'S':
1728        dnamove_show();
1729        break;
1730       
1731      case '.':
1732        dnamove_printree();
1733        break;
1734       
1735      case 'T':
1736        try();
1737        break;
1738       
1739      case 'U':
1740        undo();
1741        break;
1742       
1743      case 'W':
1744        treewrite(done);
1745        break;
1746       
1747      case 'O':
1748        changeoutgroup();
1749        break;
1750       
1751      case 'F':
1752        flip();
1753        break;
1754       
1755      case 'H':
1756        window(left, &leftedge, &topedge, hscroll, vscroll, treelines,
1757                 screenlines, screenwidth, farthest, subtree);
1758        dnamove_printree();
1759        break;
1760
1761      case 'J':
1762        window(downn, &leftedge, &topedge, hscroll, vscroll, treelines,
1763                 screenlines, screenwidth, farthest, subtree);
1764        dnamove_printree();
1765        break;
1766
1767      case 'K':
1768        window(upp, &leftedge, &topedge, hscroll, vscroll, treelines,
1769                 screenlines, screenwidth, farthest, subtree);
1770        dnamove_printree();
1771        break;
1772
1773      case 'L':
1774        window(right, &leftedge, &topedge, hscroll, vscroll, treelines,
1775                 screenlines, screenwidth, farthest, subtree);
1776        dnamove_printree();
1777        break;
1778
1779      case 'C':
1780        clade();
1781        break;
1782       
1783      case '?':
1784        help("site");
1785        dnamove_printree();
1786        break;
1787       
1788      case 'X':
1789        done = true;
1790        break;
1791       
1792      case 'Q':
1793        done = true;
1794        break;
1795      }
1796    }
1797  } while (!done);
1798  if (written)
1799    return;
1800  do {
1801    printf("Do you want to write out the tree to a file? (Y or N) ");
1802#ifdef WIN32
1803    phyFillScreenColor();
1804#endif
1805    scanf("%c%*[^\n]", &ch);
1806    getchar();
1807    if (ch == 'Y' || ch == 'y')
1808      treewrite(done);
1809  } while (ch != 'Y' && ch != 'y' && ch != 'N' && ch != 'n');
1810}  /* redisplay */
1811
1812
1813void treeconstruct()
1814{
1815  /* constructs a binary tree from the pointers in treenode. */
1816
1817  restoring = false;
1818  subtree = false;
1819  display = false;
1820  dispchar = 0;
1821  earlytree = true;
1822  waswritten = false;
1823  buildtree();
1824  printf("\nComputing steps needed for compatibility in sites ...\n\n");
1825  setorder();
1826  mincomp();
1827  newtree = true;
1828  earlytree = false;
1829  dnamove_printree();
1830  bestyet = -like;
1831  gotlike = -like;
1832  lastop = none;
1833  newtree = false;
1834  written = false;
1835  redisplay();
1836}  /* treeconstruct */
1837
1838
1839int main(int argc, Char *argv[])
1840{ /* Interactive DNA parsimony */
1841  /* reads in spp, chars, and the data. Then calls treeconstruct to
1842     construct the tree and query the user */
1843#ifdef MAC
1844  argc = 1;                /* macsetup("Dnamove","");        */
1845  argv[0] = "Dnamove";
1846#endif
1847  init(argc, argv);
1848  progname = argv[0];
1849  strcpy(infilename,INFILE);
1850  strcpy(intreename,INTREE);
1851  strcpy(outtreename,OUTTREE);
1852
1853  openfile(&infile,infilename,"input file", "r",argv[0],infilename);
1854  openfile(&outtree,outtreename,"output file", "w",argv[0],outtreename);
1855
1856  garbage = NULL;
1857  screenlines = 24;
1858  scrollinc = 20;
1859  screenwidth = 80;
1860  topedge = 1;
1861  leftedge = 1;
1862  ibmpc = IBMCRT;
1863  ansi = ANSICRT;
1864  doinput();
1865  configure();
1866  treeconstruct();
1867  FClose(infile);
1868  FClose(outtree);
1869#ifdef MAC
1870  fixmacfile(outtreename);
1871#endif
1872#ifdef WIN32
1873  phyRestoreConsoleAttributes();
1874#endif
1875  return 0;
1876}  /* Interactive DNA parsimony */
Note: See TracBrowser for help on using the repository browser.