/* DNATREE version 1.2. (c) Copyright 1998-2004 by The University of
   Washington.  Written by Joseph Felsenstein, Akiko Fuseki, Sean Lamont,
   and Andrew Keeffe.
   Permission is granted to copy and use this program provided no fee is
   charged for it and provided that this copyright notice is not removed. */

#include "dnatree.h"

#define nmlngth        10    /* number of characters in species name    */
#define maxsp         100    /* maximum number of species */
#define maxsz         200    /* twice the previous value */
#define maxsites    20000    /* maximum number of sites */
#define downn           2    /* this has to do with spacing on tree diagram */
#define overr           4    /* ... and so does this */

#define INPUT stdin

typedef long *steptr;
typedef enum {  A, C, G, U, O } bases;
typedef long *baseptr;
typedef enum {
  horiz, vert, up, over, upcorner, downcorner, aa, cc, gg, tt, quest
  } chartype;
/* nodes will form a binary tree */
typedef enum {
  rearr, flipp, reroott, none
  } rearrtype;

typedef struct node {         /* describes a tip species or an ancestor */
  struct node *next, *back;   /* pointers to nodes                      */
  long index;                 /* number of the node                     */
  boolean tip;                /* present species are tips of tree       */
  double tyme;                /* time in simulation                     */
  double cumprob;
  double probxi, probxv;      /* for treedna part                       */
  double v;                   /* for treedna part                       */
  baseptr base;               /* the sequence                           */
  Char state;                 /* stores state for printing              */
  long xcoord, ycoord;        /* used by printree                       */
  long ymin, ymax;
} node;

typedef Char naym[nmlngth];
typedef struct gbase {
  baseptr base;
  struct gbase *next;
} gbase;

typedef enum {
  ran, arb, use, spec
  } howtree;

typedef long vector[maxsites];   /* a data vector. values 0 to 3 for A to T */
typedef long longer[6];          /* array stores a random number */

char outfilename[100],trfilename[100];
node *root;
FILE *outfile, *treefile;
long spp, nonodes, nsites, outgrno, screenlines, col;
/* spp = number of species
   nonodes = number of nodes in tree
   nsites = number of sites in actual sequences
   outgrno indicates outgroup */
boolean treesim, outgropt, weights, thresh,  waswritten, interleaved,
  ibmpc, ansi;
steptr weight;
naym *nayme;   /* names of species */
double threshold;
double *threshwt;
boolean reversed[11];
boolean graphic[11];
Char chh[11];
howtree how;
gbase *garbage;
Char progname[20];
Char  ch;
/* Local variables for treeconstruct, propogated global for C version: */

long dispchar, atwhat, what, fromwhere, towhere, oldoutgrno, compatible;
double like, bestyet, gotlike;
boolean display, newtree, changed, subtree, written, oldwritten, restoring,
  wasleft, oldleft, earlytree;
steptr necsteps;
boolean *intree;
long sett[31];
steptr numsteps;
node *nuroot;
rearrtype lastop;

double ttratio, aaa, bbb;
double pchange;
long i, theseed, ntarget;
longer seed;
node *root;

long j, spp, nonodes;
node **treenode;               /* will be array of nodeptrs to tips */
vector data[maxsp];            /* the data matrix */
long sporder[maxsp];


void uppercase(ch)
     Char *ch;
{  /* convert a character to upper case -- either ASCII or EBCDIC */
  *ch = (islower(*ch) ?  toupper(*ch) : (*ch));
}  /* uppercase */


int filexists(filename)
     char *filename;
{
  /* check whether this file already exists */
  FILE *fp;
  fp = fopen(filename,"r");
  if (fp) {
    fclose(fp);
    return 1;
  } else
    return 0;
} /* filexists */


void openfile(fp, filename, filedesc, mode)
FILE **fp;
char *filename;
char *filedesc;
char *mode;
{
  /* open a file, testing whether it exists etc. */
  FILE *of;
  char file[100];
  char input[100];
  char ch;
  char *p;
  boolean done;

  p = strchr(filename, '\n');
  if (p)
    *p = '\0';
  strcpy(file,filename);
  done = false;
  while (!done) {
    if (mode[0] == 'w' && filexists(file)){
      printf("\nDnatree: the file  %s  that you wanted to\n", filename);
      printf("     use as %s already exists.\n", filedesc);
      printf("     Do you want to Replace it, Append to it,\n");
      printf("     write to a new File, or Quit?\n");
      do {
        printf("     (please type R, A, F, or Q) \n");
        fgets(input, 100, INPUT);
        ch  = input[0];
        uppercase(&ch);
      } while (ch != 'A' && ch != 'R' && ch != 'F' && ch != 'Q');
      if (ch == 'Q')
        exit(-1);
      if (ch == 'A') {
        strcpy(mode, "a");
        continue;
      }
      else if (ch == 'F') {
          file[0] = '\0';
          while (file[0] =='\0') {
            printf("Please enter a new file name> ");
            fgets(file, 100, INPUT);
          }
          strcpy(mode,"w");
          continue;
        }
      }
    of = fopen(file, mode);
    if (of) {
      done = true;
    }
    else {
      switch (mode[0]){

      case 'r':
        printf("Dnatree: can't find %s  %s\n", filedesc, file);
        file[0] = '\0';
        while (file[0] =='\0'){
          printf("Please enter a new file name> ");
          fgets(file, 100, INPUT);
        }
        break;

      case 'w':
      case 'a':
        printf("Dnatree: can't write %s file %s\n", filedesc, file);
        file[0] = '\0';
        while (file[0] =='\0') {
          printf("Please enter a new file name> ");
          fgets(file, 100, INPUT);
        }
        continue;
      default:
        printf("There is some error in the call of openfile. Unknown mode.\n");
        exit(-1);
      }
    }
  }
  *fp = of;
} /* openfile */


void gnu(p)
     gbase **p;
{
  /* this and the following are do-it-yourself garbage collectors.
     Make a new node or pull one off the garbage list */
  if (garbage != NULL) {
    *p = garbage;
    garbage = garbage->next;
  } else {
    *p = (gbase *)Malloc(sizeof(gbase));
    (*p)->base = (baseptr)Malloc(nsites*sizeof(long));
  }
  (*p)->next = NULL;
}  /* gnu */


void chuck(p)
     gbase *p;
{
  /* collect garbage on p -- put it on front of garbage list */
  p->next = garbage;
  garbage = p;
}  /* chuck */


double randum(seed)
long *seed;
{
  /* random number generator -- slow but machine independent */
  long i, j, k, sum;
  longer mult, newseed;
  double x;

  mult[0] = 13;            /* multiplier is 1664525 */
  mult[1] = 24;
  mult[2] = 22;
  mult[3] = 6;
  for (i = 0; i <= 5; i++)
    newseed[i] = 0;
  for (i = 0; i <= 5; i++) {  /* multiple modulo 2**32 */
    sum = newseed[i];
    k = i;
    if (i > 3)
      k = 3;
    for (j = 0; j <= k; j++)
      sum += mult[j] * seed[i - j];
    newseed[i] = sum;
    for (j = i; j <= 4; j++) {
      newseed[j + 1] += newseed[j] / 64;
      newseed[j] &= 63;
    }
  }
  memcpy(seed, newseed, sizeof(longer));
  seed[5] &= 3;
  x = 0.0;
  for (i = 0; i <= 5; i++)
    x = x / 64.0 + seed[i];
  x /= 4.0;
  return x;
}  /* randum */


void setuptree()
{
  /* initialize nodes and pointers for not-yet-connected tree */
  long i, j;
  node *p, *q;

  treenode = (node **)Malloc(nonodes*sizeof(node *));
  for (i = 0; i < (spp); i++) {
    treenode[i] = (node *)Malloc(sizeof(node));
    treenode[i]->base = (baseptr)Malloc(nsites*sizeof(long));
  }
  for (i = spp; i < (nonodes); i++) {
    q = NULL;
    q = NULL;
  }
  for (i = spp; i < (nonodes); i++) {
    q = NULL;
    for (j = 1; j <= 3; j++) {
      p = (node *)Malloc(sizeof(node));
      p->base = (baseptr)Malloc(nsites*sizeof(long));
      p->next = q;
      q = p;
    }
    p->next->next->next = p;
    treenode[i] = p;
  }
  for (i = 1; i <= nonodes; i++) {
    treenode[i - 1]->back = NULL;
    treenode[i - 1]->tip = (i <= spp);
    treenode[i - 1]->index = i;
    if (i > spp) {
      p = treenode[i - 1]->next;
      while (p != treenode[i - 1]) {
        p->back = NULL;
        p->tip = false;
        p->index = i;
        p = p->next;
      }
    }
  }
} /* setuptree */


void configure()
{
  /* configure to machine -- set up special characters */
  chartype a;

  for (a = horiz; (long)a <= (long)quest; a = (chartype)((long)a + 1))
    reversed[(long)a] = false;
  for (a = horiz; (long)a <= (long)quest; a = (chartype)((long)a + 1))
    graphic[(long)a] = false;
  if (ibmpc) {
    chh[(long)horiz] = 205;
    graphic[(long)horiz] = true;
    chh[(long)vert] = 186;
    graphic[(long)vert] = true;
    chh[(long)up] = 186;
    graphic[(long)up] = true;
    chh[(long)over] = 205;
    graphic[(long)over] = true;
    chh[(long)upcorner] = 200;
    graphic[(long)upcorner] = true;
    chh[(long)downcorner] = 201;
    graphic[(long)downcorner] = true;
    chh[(long)aa] = 176;
    chh[(long)cc] = 178;
    chh[(long)gg] = 177;
    chh[(long)tt] = 219;
    chh[(long)quest] = '\001';
    return;
  }
  if (ansi) {
    chh[(long)horiz] = ' ';
    reversed[(long)horiz] = true;
    chh[(long)vert] = chh[(long)horiz];
    reversed[(long)vert] = true;
    chh[(long)up] = 'x';
    graphic[(long)up] = true;
    chh[(long)over] = 'q';
    graphic[(long)over] = true;
    chh[(long)upcorner] = 'm';
    graphic[(long)upcorner] = true;
    chh[(long)downcorner] = 'l';
    graphic[(long)downcorner] = true;
    chh[(long)aa] = 'a';
    reversed[(long)aa] = true;
    chh[(long)cc] = 'c';
    reversed[(long)cc] = true;
    chh[(long)gg] = 'g';
    reversed[(long)gg] = true;
    chh[(long)tt] = 't';
    reversed[(long)tt] = true;
    chh[(long)quest] = '?';
    reversed[(long)quest] = true;
    return;
  }
  chh[(long)horiz] = '=';
  chh[(long)vert] = ' ';
  chh[(long)up] = '!';
  chh[(long)upcorner] = '`';
  chh[(long)downcorner] = ',';
  chh[(long)over] = '-';
  chh[(long)aa] = 'a';
  chh[(long)cc] = 'c';
  chh[(long)gg] = 'g';
  chh[(long)tt] = 't';
  chh[(long)quest] = '.';
}  /* configure */


void inputtreedata()
{
  /* read parameters of run */
  long i;
  double x;
  Char c;

  treesim = true;
  ntarget = 10;
  nsites = 200;
  ttratio = 2.0;
  do {
    printf("\nDNATREE Molecular evolution and phylogeny version %s\n\n",
               VERSION);
    printf("Tree menu:\n");
    printf("   S                         Simulate a tree?  %s\n",
        treesim ? "Yes, simulate" : "No, read one in");
    if (treesim)
      printf("   N                       Number of species?  %ld\n", ntarget);
    printf("   C                         Number of sites?  %ld\n", nsites);
    printf("   T           Transition/transversion ratio?  %f\n", ttratio);
    printf("   0      Graphics type (IBM PC, ANSI, none)?  %s\n",
           ibmpc ? "IBM PC" :
           ansi  ? "ANSI"     : "(none)");
    printf("Type character or type Y to accept these settings\n");
    scanf("%c%*[^\n]", &c);
    uppercase(&c);
    getchar();
    if (c == 'N') {
      printf("number of total species?\n");
      scanf("%ld%*[^\n]", &ntarget);
      getchar();
    }
    if (c == 'S')
      treesim = !treesim;
    if (c == 'C') {
      printf("number of sites?\n");
      scanf("%ld%*[^\n]", &nsites);
      getchar();
    }
    if (c == 'T') {
      printf("transition/transversion ratio?\n");
      scanf("%lg%*[^\n]", &ttratio);
      getchar();
    }
    if (c == '0') {
      if (ibmpc) {
        ibmpc = false;
        ansi = false;
      } else {
        if (ansi) {
          ansi = false;
          ibmpc = true;
        } else {
          ibmpc = false;
          ansi = true;
        }
      }
    } 
  }
  while (!(c == 'Y'));
  printf("\nrandom number seed? (any whole number less than 1,000,000,000)\n");
  scanf("%ld%*[^\n]", &theseed);
  theseed = 4*theseed+1;
  seed[0] = theseed % 64;
  theseed = theseed / 64;
  seed[1] = theseed % 64;
  theseed = theseed / 64;
  seed[2] = theseed % 64;
  theseed = theseed / 64;
  seed[3] = theseed % 64;
  theseed = theseed / 64;
  seed[4] = theseed % 64;
  theseed = theseed / 64;
  seed[5] = theseed % 64;
  getchar();
  for (i = 1; i <= 100; i++)   /* to make the starting point more random */
    x = randum(seed);
  printf("prob of change per unit time?\n");
  scanf("%lf%*[^\n]", &pchange);
  getchar();
  spp = ntarget;
  nonodes = 2*spp - 1;
  configure();
}  /* inputtreedata */


void probtrav(r)
node *r;
{
  /* traverses the tree, filling in rates of change of a branch
     in the node at the top of the branch */
  double x, y;

  if (r != root)
    x = r->tyme - r->back->tyme;
  else
    x = r->tyme;
  y = pchange * x;
  if (r != root)
    r->cumprob = treenode[r->back->index-1]->cumprob + y;
  else
    r->cumprob = 0.0;
  if (!r->tip) {
    probtrav(r->next->back);
    probtrav(r->next->next->back);
  }
}  /* probtrav */


void simtree()
{
  long species, who, i, j;
  double when, increment;
  node *p, *q, *r;
  Char c;

  species = 1;
  root = treenode[0];        /* root points to first species */
  root->back = NULL;
  root->index = 1;
  when = 0.0;
  for (i = 1; i <= ntarget; i++) {  /* this loop does the branching */
    who = (long)(randum(seed) * species) + 1;   /* who will split */
      increment = -(log(randum(seed)) / species)*pchange;
                                            /* at an exponential time */
    when += increment;   /* when is time of next split */
    for (j = 0; j < i; j++)   /* update times of tips */
      treenode[j]->tyme = when;
    if (species < ntarget)   /* don't branch that last time */
      {  /* add a new tip */
      species++;
      p = treenode[ntarget+species-2];   /* new interior node */
      p->next->back = treenode[who-1];         /* connect it to  "who" */
      p->tyme = treenode[who-1]->tyme;   /* copy tip's time to new node */
      q = treenode[species-1];           /* new tip species */
      q->tyme = p->tyme;
      p->next->next->back = q;                 /* connect them all up */
      q->back = p->next->next;
      r = treenode[who-1]->back;
      treenode[who-1]->back = p;
      if (r != NULL) {
        r = treenode[r->index-1];   /* point to bottom of node below */
        if (r->next->back == treenode[who-1]) {
          r->next->back = p;
          p->back = r->next;
        }
        else {
          r->next->next->back = p;
          p->back = r->next->next;
        }
      }
      if (root == treenode[who-1]) {
        root = p;
        root->back = NULL;
      }
    }
  }
  probtrav(root);   /* assign rates of change to branches */
  for (i = 0; i < ntarget; i++) {
    if (i <= 26)
      nayme[i][1] = ' ';
    else
      nayme[i][1] = (Char)(64+(i/26));
    nayme[i][0] = (Char)(65+(i%26));
    for (j = 2; j < nmlngth; j++)
      nayme[i][j] = ' ';
  }
  for (i = 0; i < ntarget; i++) {  /* place tip names in random order */
    j = (long)(randum(seed) * ntarget) + 1;
    c = nayme[i][0];
    nayme[i][0] = nayme[j-1][0];
    nayme[j-1][0] = c;
    c = nayme[i][1];
    nayme[i][1] = nayme[j-1][1];
    nayme[j-1][1] = c;
  }
}  /* simtree */


void findch(c)
     Char c;
{
  /* scan forward until find character c */
  boolean done;

  done = false;
  while (!(done)) {
    if (c == ',') {
      if (ch == '(' || ch == ')' || ch == ';') {
   printf("\nERROR IN USER TREE: UNMATCHED PARENTHESIS OR MISSING COMMA\n");
        exit(-1);
      } else if (ch == ',')
        done = true;
    } else if (c == ')') {
      if (ch == '(' || ch == ',' || ch == ';') {
printf("\nERROR IN USER TREE: UNMATCHED PARENTHESIS OR NOT BIFURCATED NODE\n");
        exit(-1);
      } else {
        if (ch == ')')
          done = true;
      }
    } else if (c == ';') {
      if (ch != ';') {
  printf("\nERROR IN USER TREE: UNMATCHED PARENTHESIS OR MISSING SEMICOLON\n");
        exit(-1);
      } else
        done = true;
    }
    if ((done && ch == ')') || !(done)) {
      if (eoln(treefile)) {
        fscanf(treefile, "%*[^\n]");
        getc(treefile);
      }
      ch = getc(treefile);
    }
  }
}  /* findch */


void addelement(p, nextnode,lparens)
     node **p;
     long *nextnode,*lparens;
{
  /* recursive procedure adds nodes to user-defined tree */
  node *q;
  long i, n;
  boolean found;
  Char str[nmlngth];

  do {
    if (eoln(treefile)) {
      fscanf(treefile, "%*[^\n]");
      getc(treefile);
    }
    ch = getc(treefile);
  } while (ch == ' ' );
  if (ch == '(' ) {
    if (*lparens>= spp - 1) {
      printf("\nERROR IN USER TREE: TOO MANY LEFT PARENTHESES\n");
      exit(-1);
    }
    (*nextnode)++;
    (*lparens)++;
    q = treenode[*nextnode - 1];
    addelement(&q->next->back, nextnode,lparens);
    q->next->back->back = q->next;
    findch(',');
    addelement(&q->next->next->back, nextnode,lparens);
    q->next->next->back->back = q->next->next;
    findch(')');
    *p = q;
    return;
  }
  for (i = 0; i < nmlngth; i++)
    str[i] = ' ';
  n = 1;
  do {
    if (ch == '_')
      ch = ' ';
    str[n - 1] = ch;
    if (eoln(treefile)) {
      fscanf(treefile, "%*[^\n]");
      getc(treefile);
    }
    ch = getc(treefile);
    n++;
  } while (ch != ',' && ch != ')' && ch != ':' && n <= nmlngth);
  n = 1;
  do {
    found = true;
    for (i = 0; i < nmlngth; i++)
      found = (found && str[i] == nayme[n - 1][i]);
    if (found) {
      if (intree[n - 1] == false) {
        *p = treenode[n - 1];
        intree[n - 1] = true;
      } else {
        printf("\nERROR IN USER TREE: DUPLICATE NAME FOUND -- ");
        for (i = 0; i < nmlngth; i++)
          putchar(nayme[n - 1][i]);
        putchar('\n');
        exit(-1);
      }
    } else
      n++;
  } while (!(n > spp || found));
  if (n <= spp)
    return;
  printf("CANNOT FIND SPECIES: ");
  for (i = 0; i < nmlngth; i++)
    putchar(str[i]);
  putchar('\n');
}  /* addelement */


void treeread()
{
  /* read in user-defined tree and set it up */
  long nextnode, lparens;
  long i;

  openfile(&treefile, trfilename, "tree file", "r");
  root = treenode[spp];
  nextnode = spp;
  root->back = NULL;
  for (i = 0; i < (spp); i++)
    intree[i] = false;
  lparens = 0;
  addelement(&root, &nextnode,&lparens);
  findch(';');
  printf("\n\n");
  FClose(treefile);
}  /* treeread */


Void probmake(r)
node *r;
{
  /* traverses the tree, filling in rates of change of a branch
     in the node at the top of the branch */
  double x;

  if (r != root)
    x = r->tyme - treenode[r->back->index-1]->tyme;
  else
    x = 0.0;
  r->probxi = 0.25 + 0.25 * exp(-aaa * x) - 0.5 * exp(-(aaa + bbb) * x);
  r->probxv = 0.5 * (1.0 - exp(-aaa * x));
  if (!r->tip) {
    probmake(r->next->back);
    probmake(r->next->next->back);
  }
}  /* probmake */


Void evtrav(p, oldstate)
node *p;
long oldstate;
{
  /* preorder traverse the tree, doing the evolving on each branch */
  long state;
  double y;

  y = randum(seed);
  if (y < p->probxi) {
    if (oldstate < 2)
      state = 1 - oldstate;
    else
      state = 5 - oldstate;
  } else if (y < p->probxi + p->probxv) {
    if (oldstate < 2)
      state = (long)(2 * randum(seed)) + 2;
    else
      state = (long)(2 * randum(seed));
  } else
    state = oldstate;
  if (p->tip)
    data[p->index - 1][j - 1] = state;
  else {
    evtrav(p->next->back, state);
    evtrav(p->next->next->back, state);
  }
}  /* evtrav */


Void evolve()
{
  /* evolve the characters on that tree */
  evtrav(root, (long)(4 * randum(seed)));   /* starting from a random base */
}  /* evolve */


Void sort()
{
  /* Shell sort keeping species names in order */
  long gap, i, j, itemp;

  for (i = 1; i <= spp; i++)
    sporder[i - 1] = i;
  gap = spp / 2;
  while (gap > 0) {
    for (i = gap + 1; i <= spp; i++) {
      j = i - gap;
      while (j > 0) {
        if (nayme[sporder[j - 1] - 1][0] >
            nayme[sporder[j + gap - 1] - 1][0] ||
            (nayme[sporder[j - 1] - 1][0] ==
             nayme[sporder[j + gap - 1] - 1][0] &&
             nayme[sporder[j - 1] - 1][1] >
             nayme[sporder[j + gap - 1] - 1][1])) {
          itemp = sporder[j - 1];
          sporder[j - 1] = sporder[j + gap - 1];
          sporder[j + gap - 1] = itemp;
        }
        j -= gap;
      }
    }
    gap /= 2;
  }
}  /* sort */


void displaydata()
{ /* display part of data on screen */
   long i, j, startsite;
   Char ch;
   boolean dots;

  printf("\nHere are the sequences that were simulated along that tree.\n");
printf(" A dot (.) means that the base is the same as in the first species.\n");
   dots = true;
   startsite = 1;
   do {
     printf("\nHere are sites %ld through %ld out of a total of %ld.\n",
          startsite, (startsite+59 > nsites) ? nsites : startsite+59, nsites);
     printf("------------------------------------");
     printf("-----------------------------------------\n");
     for (i = 1; i <= spp; i++) {
       for (j = startsite; (j <= nsites) && (j < startsite+60); j++) {
         if (j == startsite)
           printf("%c%c          ", nayme[sporder[i-1]-1][0],
                                    nayme[sporder[i-1]-1][1]);
         if (dots && (i != 1)
                  && (data[sporder[i-1]-1][j-1] == data[sporder[0]-1][j-1]))
             printf(".");
         else
           switch (data[sporder[i-1] - 1][j-1]) {
  
           case 0:
             printf("A");
             break;
    
           case 1:
             printf("G");
             break;
    
           case 2:
             printf("C");
             break;
  
           case 3:
             printf("T");
             break;
  
         }
         if ((j > startsite) && ((j%10) == 0))
           printf(" ");
       }
       printf("\n");
     }
     printf("------------------------------------");
     printf("-----------------------------------------\n");
     printf("\nChoose what to do next:\n");
     if (startsite+60 <= nsites)
       printf("  N    See the next group of sites\n");
     if (startsite>1)
       printf("  P    See the previous group of sites\n");
     if (dots)
       printf("  D    Turn off using dot-differencing\n");
     else
       printf("  D    Turn on using dot-differencing\n");
     printf("  Y    Finish looking at sequences\n");
     printf("Type one of these letters, then press on Enter key\n");
     scanf("%c%*[^\n]", &ch);
     uppercase(&ch);
     getchar();
     if ((startsite+60 <= nsites) && (ch == 'N'))
       startsite = startsite + 60;
     if ((startsite > 1) && (ch == 'P'))
       startsite = startsite-60;
     if (ch == 'D')
       dots = !dots;
   } while (ch != 'Y');
} /* displaydata */


void outdata()
{
  /* output the data */
  long i, j, k;

  fprintf(outfile, "%5ld%5ld\n", spp, nsites);   /* output first line */
  for (k = 1; 60*(k-1) < nsites; k++) {
    for (i = 0; i < spp; i++) {  /* output the species one by one */
      if (k == 1) {
        fprintf(outfile, "%c%c", nayme[sporder[i] - 1][0],
                                 nayme[sporder[i] - 1][1]);
        fprintf(outfile, "        ");
      } else
        fprintf(outfile, "          ");
      for (j = 0; j < 60; j++) {
        if (60*(k-1)+j < nsites) {
          switch (data[sporder[i] - 1][60*(k-1)+j]) {
  
          case 0:
            putc('A', outfile);
            break;
  
          case 1:
            putc('G', outfile);
            break;
  
          case 2:
            putc('C', outfile);
            break;

          case 3:
            putc('T', outfile);
            break;
          }
        }
      }
      putc('\n', outfile);
    }
    putc('\n', outfile);
  }
  if (outfile != NULL)
    fclose(outfile);
}  /* outdata */


void treednamain()
{  /* main program -- pretty self explanatory */
  Char ch;

  outfile = NULL;
  garbage = NULL;
  aaa = 2.0 / (1.0 + ttratio);
  bbb = (2.0 * ttratio - 1.0) / (ttratio + 1.0);
/*  treeread();                       debug may use this later  */
  probmake(root);
  for (j = 1; j <= nsites; j++)
    evolve();
  sort();
  displaydata();
  do {
    printf("Do you want to write the data out onto a file? (Y or N) ");
    scanf("%c%*[^\n]", &ch);
    uppercase(&ch);
    getchar();
    if (ch == 'Y') {
      printf("please type name of that file: ");
      fgets(outfilename, 100, INPUT);
      openfile(&outfile, outfilename, "data file", "w");
      outdata();
    }
  } while ((ch  != 'Y') && (ch != 'N'));
} /* treednamain */


void inpnum(n, success)
long *n;
boolean *success;
{
  int fields;
  char line[100];
  fgets(line, 100, INPUT);
  *n = atof(line);
  fields = sscanf(line,"%ld",n);
  *success = (fields == 1);

}  /* inpnum */


void getoptions()
{
  /* interactively set options */
  Char ch;
  boolean done, done1, gotopt;

  how = ran;
  outgrno = 1;
  outgropt = false;
  thresh = false;
  weights = false;
  do {
    printf((ansi) ?  "\033[2J\033[H" : "\n");
    printf("\nInteractive search for most parsimonious tree\n\n");
    printf("Settings for this run:\n");
    printf("  O                           Outgroup root?");
    if (outgropt)
      printf("  Yes, at sequence number%3ld\n", outgrno);
    else
      printf("  No, use as outgroup species%3ld\n", outgrno);
    printf("  T                 Use Threshold parsimony?");
    if (thresh)
      printf("  Yes, count up to %4.1f per site\n", threshold);
    else
      printf("  No, use ordinary parsimony\n");
    printf("  L               Number of lines on screen?%4ld\n",screenlines);
    do {
      printf("\nAre these settings correct? ");
      printf("(type Y or the letter for one to change)\n");
      scanf("%c%*[^\n]", &ch);
      getchar();
      uppercase(&ch);
      done = (ch == 'Y');
      gotopt = (strchr("OTL",ch) != NULL) ? true : false;
      if (gotopt) {
        switch (ch) {
        
        case 'O':
          outgropt = !outgropt;
          if (outgropt) {
            done1 = true;
            do {
              printf("Type number of the outgroup:\n");
              scanf("%ld%*[^\n]", &outgrno);
              getchar();
              done1 = (outgrno >= 1 && outgrno <= spp);
              if (!done1) {
                printf("BAD OUTGROUP NUMBER: %4ld\n", outgrno);
                printf("  Must be in range 1 -%2ld\n", spp);
              }
            } while (done1 != true);
          }
          break;
        
        case 'T':
          thresh = !thresh;
          if (thresh) {
            done1 = false;
            do {
              printf("What will be the threshold value?\n");
              scanf("%lf%*[^\n]", &threshold);
              getchar();
              done1 = (threshold >= 1.0);
              if (!done1)
                printf("BAD THRESHOLD VALUE:  it must be greater than 1\n");
              else
                threshold = (long)(threshold * 10.0 + 0.5) / 10.0;
            } while (done1 != true);
          }
          break;
        
        case 'L':
          do {
            printf("Number of lines on screen?\n");
            scanf("%ld%*[^\n]", &screenlines);
            getchar();
          } while (screenlines <= 12);
          break;
        }
      }
      if (!(gotopt || done))
        printf("Not a possible option!\n");
    } while (!(gotopt || done));
  } while (!done);
}  /* getoptions */


void inputoptions()
{
  /* input the information on the options */
  long i;

  if (!thresh)
    threshold = spp;
  for (i = 0; i < (nsites); i++)
    threshwt[i] = threshold * weight[i];
}  /* inputoptions */


void doinput()
{ 
  /* reads the input data */
  long i;

  getoptions();
  intree = (boolean *)Malloc(nonodes*sizeof(boolean));
  weight = (steptr)Malloc(nsites*sizeof(long));
  numsteps = (steptr)Malloc(nsites*sizeof(long));
  necsteps = (steptr)Malloc(nsites*sizeof(long));
  threshwt = (double *)Malloc(nsites*sizeof(double));
  printf("\nReading input file ...\n\n");
  for (i=0; i<nsites; i++)
    weight[i] = 1;
  inputoptions();
}  /* doinput */


void pregraph()
{
  /* turn on graphic characters */
  if (ansi) {
    printf("\033(0");
    printf("\033[10m");
  }
}  /* pregraph */


void prereverse()
{
  /* turn on reverse video */
  if (ansi)
    printf("\033[7m");
}  /* prereverse */


void prefix(a)
     chartype a;
{
  /* give prefix appropriate for this character */
  if (reversed[(long)a])
    prereverse();
  if (graphic[(long)a])
    pregraph();
}  /* prefix */


void postgraph()
{
  /* turn off graphic characters */
  if (ansi) {
    printf("\033[11m");
    printf("\033(B");
  }
}  /* postgraph */


void postreverse()
{
  /* turn off reverse video */
  if (ansi)
    printf("\033[0m");
}  /* postreverse */


void postfix(a)
     chartype a;
{
  /* give postfix appropriate for this character */
  if (reversed[(long)a])
    postreverse();
  if (graphic[(long)a])
    postgraph();
}  /* postfix */


void makechar(a)
     chartype a;
{
  /* print out a character with appropriate prefix or postfix */
  prefix(a);
  putchar(chh[(long)a]);
  postfix(a);
}  /* makechar */


void add(below, newtip, newfork)
     node *below, *newtip, *newfork;
{
  /* inserts the nodes newfork and its left descendant, newtip,
     to the tree.  below becomes newfork's right descendant */
  boolean putleft;
  node *leftdesc, *rtdesc;

  if (below != treenode[below->index - 1])
    below = treenode[below->index - 1];
  if (below->back != NULL)
    below->back->back = newfork;
  newfork->back = below->back;
  putleft = true;
  if (restoring)
    putleft = wasleft;
  if (putleft) {
    leftdesc = newtip;
    rtdesc = below;
  } else {
    leftdesc = below;
    rtdesc = newtip;
  }
  rtdesc->back = newfork->next->next;
  newfork->next->next->back = rtdesc;
  newfork->next->back = leftdesc;
  leftdesc->back = newfork->next;
  if (root == below)
    root = newfork;
  root->back = NULL;
}  /* add */


void re_move(item, fork)
     node **item, **fork;
{
  /* removes nodes item and its ancestor, fork, from the tree.
     the new descendant of fork's ancestor is made to be
     fork's second descendant (other than item).  Also
     returns pointers to the deleted nodes, item and fork */
  node *p, *q;

  if ((*item)->back == NULL) {
    *fork = NULL;
    return;
  }
  *fork = treenode[(*item)->back->index - 1];
  if (*item == (*fork)->next->back) {
    if (root == *fork)
      root = (*fork)->next->next->back;
    wasleft = true;
  } else {
    if (root == *fork)
      root = (*fork)->next->back;
    wasleft = false;
  }
  p = (*item)->back->next->back;
  q = (*item)->back->next->next->back;
  if (p != NULL)
    p->back = q;
  if (q != NULL)
    q->back = p;
  (*fork)->back = NULL;
  p = (*fork)->next;
  while (p != *fork) {
    p->back = NULL;
    p = p->next;
  }
  (*item)->back = NULL;
}  /* remove */


void fillin(p)
     node *p;
{
  /* sets up for each node in the tree the base sequence
     at that point and counts the changes.  The program
     spends much of its time in this PROCEDURE */
  long i;
  long ns, rs, ls;

  for (i = 0; i < (nsites); i++) {
    ls = p->next->back->base[i];
    rs = p->next->next->back->base[i];
    ns = ls & rs;
    if (ns == 0) {
      ns = ls | rs;
      numsteps[i] += weight[i];
    }
    p->base[i] = ns;
  }
}  /* fillin */


void postorder(p)
     node *p;
{
  /* traverses a binary tree, calling PROCEDURE fillin at a
     node's descendants before calling fillin at the node */

  if (p->tip)
    return;
  postorder(p->next->back);
  postorder(p->next->next->back);
  fillin(p);
}  /* postorder */


void evaluate(r)
     node *r;
{
  /* determines the number of steps needed for a tree. this is
     the minimum number of steps needed to evolve sequences on
     this tree */
  long i, steps;
  double sum;

  compatible = 0;
  sum = 0.0;
  for (i = 0; i < (nsites); i++)
    numsteps[i] = 0;
  postorder(r);
  for (i = 0; i < (nsites); i++) {
    steps = numsteps[i];
    if (steps <= threshwt[i])
      sum += steps;
    else
      sum += threshwt[i];
    if (steps <= necsteps[i] && !earlytree)
      compatible += weight[i];
  }
  like = -sum;
  /*printf("like: %f\n",like);*/
}  /* evaluate */


void reroot(outgroup)
     node *outgroup;
{
  /* reorients tree, putting outgroup in desired position. */
  node *p, *q, *newbottom, *oldbottom;
  boolean onleft;

  if (outgroup->back->index == root->index)
    return;
  newbottom = outgroup->back;
  p = treenode[newbottom->index - 1]->back;
  while (p->index != root->index) {
    oldbottom = treenode[p->index - 1];
    treenode[p->index - 1] = p;
    p = oldbottom->back;
  }
  onleft = (p == root->next);
  if (restoring) {
    if (!onleft && wasleft){
      p = root->next->next;
      q = root->next;
    } else {
      p = root->next;
      q = root->next->next;
    }
  }
  else {
    if (onleft)
      oldoutgrno = root->next->next->back->index;
    else
      oldoutgrno = root->next->back->index;
    wasleft = onleft;
    p = root->next;
    q = root->next->next;
  }
  p->back->back = q->back;
  q->back->back = p->back;
  p->back = outgroup;
  q->back = outgroup->back;
  if (restoring) {
    if (!onleft && wasleft) {
      outgroup->back->back = root->next;
      outgroup->back = root->next->next;
    } else {
      outgroup->back->back = root->next->next;
      outgroup->back = root->next;
    }
  } else {
    outgroup->back->back = root->next->next;
    outgroup->back = root->next;
  }
  treenode[newbottom->index - 1] = newbottom;
}  /* reroot */


void ancestset(a, b, c)
     long a, b, *c;
{
  /* make the set of ancestral states below nodes
     whose base sets are a and b */
  *c = a & b;
  if (*c == 0)
    *c = a | b;
}  /* ancestset */


void firstrav(r,i)
     node *r;
     long i;
{
  /* initial traverse for hypothetical states */
  if (r->tip)
    return;
  firstrav(r->next->back, i);
  firstrav(r->next->next->back, i);
  ancestset(r->next->back->base[i - 1],
            r->next->next->back->base[i - 1], &r->base[i - 1]);
}  /* firstrav */


void hyptrav(r, hypset, i,bottom)
     node *r;
     long *hypset;
     long i;
     boolean *bottom;
{
  /*  compute, print out state at one interior node */
  long tempset, left, rt, anc;
  gbase *temparray, *ancset;

  gnu(&ancset);
  gnu(&temparray);
  anc = hypset[i - 1];
  if (!r->tip) {
    left = r->next->back->base[i - 1];
    rt = r->next->next->back->base[i - 1];
    tempset = left & rt & anc;
    if (tempset == 0) {
      tempset = (left & rt) | (left & anc) | (rt & anc);
      if (tempset == 0)
        tempset = left | rt | anc;
    }
    r->base[i - 1] = tempset;
  }
  if (!(*bottom))
    anc = treenode[r->back->index - 1]->base[i - 1];
  r->state = '?';
  if (r->base[dispchar - 1] == 1L << ((long)A))
    r->state = 'A';
  if (r->base[dispchar - 1] == 1L << ((long)C))
    r->state = 'C';
  if (r->base[dispchar - 1] == 1L << ((long)G))
    r->state = 'G';
  if (r->base[dispchar - 1] == 1L << ((long)U))
    r->state = 'T';
  /*printf("%c",r->state);*/
  *bottom = false;
  if (!r->tip) {
    memcpy(temparray->base, r->next->back->base, nsites*sizeof(long));
    ancestset(hypset[i - 1], r->next->next->back->base[i - 1],
              &ancset->base[i - 1]);
    hyptrav(r->next->back, ancset->base, i,bottom);
    ancestset(hypset[i - 1], temparray->base[i - 1],
              &ancset->base[i - 1]);
    hyptrav(r->next->next->back, ancset->base, i,bottom);
  }
  chuck(temparray);
  chuck(ancset);
}  /* hyptrav */


void hypstates()
{
  /* fill in and describe states at interior nodes */
  /* Local variables for hypstates: */
  long i;
  boolean bottom;
  baseptr nothing;

  i = dispchar;
  nothing = (baseptr)Malloc(nsites*sizeof(long));
  nothing[i - 1] = 0;
  bottom = true;
  firstrav(root, i);
  hyptrav(root, nothing, i,&bottom);
  free(nothing);
}  /* hypstates */


void coordinates(simulated, p, tipy)
     boolean simulated;
     node *p;
     long *tipy;
{
  /* establishes coordinates of nodes */
  node *q, *first, *last;

  if (p->tip) {
    p->xcoord = 0;
    p->ycoord = *tipy;
    p->ymin = *tipy;
    p->ymax = *tipy;
    (*tipy) += downn;
    return;
  }
  q = p->next;
  do {
    coordinates(simulated, q->back, tipy);
    q = q->next;
  } while (p != q);
  first = p->next->back;
  q = p->next;
  while (q->next != p)
    q = q->next;
  last = q->back;
  if (!simulated)
    p->xcoord = (last->ymax - first->ymin) * 3 / 2;
  else
    p->xcoord =
       ((treenode[0]->tyme-p->tyme)/(treenode[0]->tyme-root->tyme))*55;
  p->ycoord = (first->ycoord + last->ycoord) / 2;
  p->ymin = first->ymin;
  p->ymax = last->ymax;
}  /* coordinates */


void drawline(simulated, i, lastline)
     boolean simulated;
     long i, lastline;
{
  /* draws one row of the tree diagram by moving up tree */
  node *p, *q, *r, *first, *last;
  long n, j;
  boolean extra, done;
  Char st;
  chartype c, d;

  p = nuroot;
  q = nuroot;
  extra = false;
  if (i == p->ycoord && (p == root || subtree)) {
    c = over;
    if (display) {
      switch (p->state) {
        
      case 'A':
        c = aa;
        break;
        
      case 'C':
        c = cc;
        break;
        
      case 'G':
        c = gg;
        break;
        
      case 'T':
        c = tt;
        break;
        
      case '?':
        c = quest;
        break;
      }
    }
    if (p->index >= 10) {
      makechar(c);
      printf("%2ld", p->index);
    } else {
      makechar(c);
      makechar(c);
      printf("%ld", p->index);
    }
    extra = true;
  } else
    printf("  ");
  do {
    if (!p->tip) {
      r = p->next;
      done = false;
      do {
        if (i >= r->back->ymin && i <= r->back->ymax) {
          q = r->back;
          done = true;
        }
        r = r->next;
      } while (!(done || r == p));
      first = p->next->back;
      r = p->next;
      while (r->next != p)
        r = r->next;
      last = r->back;
    }
    done = (p == q);
    n = p->xcoord - q->xcoord;
    if (n < 3 && !q->tip)
      n = 3;
    if (extra) {
      n--;
      extra = false;
    }
    if (q->ycoord == i && !done) {
      if (q->ycoord > p->ycoord)
        d = upcorner;
      else
        d = downcorner;
      c = over;
      if (display) {
        switch (q->state) {
        
        case 'A':
          c = aa;
          break;
        
        case 'C':
          c = cc;
          break;
        
        case 'G':
          c = gg;
          break;
        
        case 'T':
          c = tt;
          break;
        
        case '?':
          c = quest;
          break;
        }
        d = c;
      }
      if (n > 1) {
        makechar(d);
        prefix(c);
        for (j = 1; j <= n - 2; j++)
          putchar(chh[(long)c]);
        postfix(c);
      }
      if (q->index >= 10)
        printf("%2ld", q->index);
      else {
        makechar(c);
        printf("%ld", q->index);
      }
      extra = true;
    } else if (!q->tip) {
      if (last->ycoord > i && first->ycoord < i && i != p->ycoord) {
        c = up;
        if (i < p->ycoord)
          st = p->next->back->state;
        else
          st = p->next->next->back->state;
        if (display) {
          switch (st) {
        
          case 'A':
            c = aa;
            break;
        
          case 'C':
            c = cc;
            break;
        
          case 'G':
            c = gg;
            break;
        
          case 'T':
            c = tt;
            break;
        
          case '?':
            c = quest;
            break;
          }
        }
        makechar(c);
        for (j = 1; j < n; j++)
          putchar(' ');
      } else {
        for (j = 1; j <= n; j++)
          putchar(' ');
      }
    } else {
      for (j = 1; j <= n; j++)
        putchar(' ');
    }
    if (p != q)
      p = q;
  } while (!done);
  if (p->ycoord == i && p->tip) {
    putchar(':');
    j = 0;
    while ((nayme[p->index-1][j] == ' ') && (j<2))
      j++;
    putchar(nayme[p->index - 1][j]);
    if ((j == 0) && (nayme[p->index-1][j+1] != ' '))
      putchar(nayme[p->index - 1][1]);
  }
  if (i == 1) {
    if (!simulated) {
      if (display)
        printf("       SITE%4ld", dispchar);
      else
        printf("           ");
      if (subtree)
        printf(" Subtree ");
      else
        printf("         ");
    }
  }
  if (i == 3 && display) {
    printf("   ");
    makechar(aa);
    printf(":A, ");
    makechar(cc);
    printf(":C, ");
    makechar(gg);
    printf(":G, ");
    makechar(tt);
    printf(":T, ");
    makechar(quest);
    printf(":?");
  }
  if ((i == 5 || (lastline < 5 && i == lastline)) && !earlytree) {
    if ((lastline < 5) && (!simulated))
      printf("\n                   ");
    if (!simulated)
      printf("   Steps:%10.5f", - like);
    gotlike = -like;
  }
  if ((i == 7 || (lastline < 7 && i == lastline)) && changed && !earlytree) {
    if ((lastline < 7) && (!simulated))
      printf("\n                     ");
    if (-like <bestyet) {
      printf("     BEST YET!");
      bestyet = -like;
    } else if (fabs(-like - bestyet) < 0.000001)
      printf("     (as good as best)");
    else {
      if (-like < gotlike)
        printf("     better");
      else if (-like > gotlike)
        printf("     worse!");
    }
  }
  if (!simulated && (i == 9 || (lastline < 9 && i == lastline))
       && !earlytree) {
    if ((lastline < 9) && (!simulated))
      printf("\n                       ");
    printf("     %3ld sites compatible", compatible);
  }
  putchar('\n');
}  /* drawline */


void printree(simulated)
boolean simulated;
{
  /* prints out diagram of the tree */
  long tipy;
  long i, rover, dow;

  if (!subtree)
    nuroot = root;
  if (changed || newtree)
    evaluate(root);
  if (display)
    hypstates();
  printf((ansi) ? "\033[2J\033[H" : "\n");
  tipy = 1;
  rover = 100 / spp;
  if (rover > overr)
    rover = overr;
  dow = downn;
  if (spp * dow > screenlines && !subtree) {
    dow--;
    rover--;
  }
  coordinates(simulated, nuroot, &tipy);
  for (i = 1; i <= (tipy - dow); i++)
    drawline(simulated, i, tipy-dow);
  if (spp <= screenlines / 2 || subtree)
    putchar('\n');
  changed = false;
}  /* printree */


void randomtree()
{ /* make up a random tree */
  long i, n;

  root = treenode[0];
  add(treenode[0], treenode[1], treenode[spp]);
  for (i = 3; i <= (spp); i++) {
    n = (long)((2*i-3)*randum(seed));
    if (n > (i-2))
      n += spp-i+1;
    add(treenode[n], treenode[i - 1], treenode[spp + i - 2]);
  }
  for (i = 0; i < nonodes;  i++)
    intree[i] = true;
}  /* randomtree */


void buildtree()
{ long i, j;

  changed = false;
  newtree = false;
  randomtree();
  if (!outgropt)
    outgrno = root->next->back->index;
  if (outgropt && intree[outgrno - 1])
    reroot(treenode[outgrno - 1]);
  for (i=0; i < spp; i++) {
    for (j=0; j < nsites; j++) {
       switch (data[i][j]) {

       case 0:
         treenode[i]->base[j] = 1L << (long)A;
         break;

       case 1:
         treenode[i]->base[j] = 1L << (long)G;
         break;

       case 2:
         treenode[i]->base[j] = 1L << (long)C;
         break;

       case 3:
         treenode[i]->base[j] = 1L << (long)U;
         break;

       }
    }
  }
}  /* buildtree */


void setorder()
{
  /* sets in order of number of members */
  sett[0] = 1L << ((long)A);
  sett[1] = 1L << ((long)C);
  sett[2] = 1L << ((long)G);
  sett[3] = 1L << ((long)U);
  sett[4] = 1L << ((long)O);
  sett[5] = (1L << ((long)A)) | (1L << ((long)C));
  sett[6] = (1L << ((long)A)) | (1L << ((long)G));
  sett[7] = (1L << ((long)A)) | (1L << ((long)U));
  sett[8] = (1L << ((long)A)) | (1L << ((long)O));
  sett[9] = (1L << ((long)C)) | (1L << ((long)G));
  sett[10] = (1L << ((long)C)) | (1L << ((long)U));
  sett[11] = (1L << ((long)C)) | (1L << ((long)O));
  sett[12] = (1L << ((long)G)) | (1L << ((long)U));
  sett[13] = (1L << ((long)G)) | (1L << ((long)O));
  sett[14] = (1L << ((long)U)) | (1L << ((long)O));
  sett[15] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)G));
  sett[16] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)U));
  sett[17] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)O));
  sett[18] = (1L << ((long)A)) | (1L << ((long)G)) | (1L << ((long)U));
  sett[19] = (1L << ((long)A)) | (1L << ((long)G)) | (1L << ((long)O));
  sett[20] = (1L << ((long)A)) | (1L << ((long)U)) | (1L << ((long)O));
  sett[21] = (1L << ((long)C)) | (1L << ((long)G)) | (1L << ((long)U));
  sett[22] = (1L << ((long)C)) | (1L << ((long)G)) | (1L << ((long)O));
  sett[23] = (1L << ((long)C)) | (1L << ((long)U)) | (1L << ((long)O));
  sett[24] = (1L << ((long)G)) | (1L << ((long)U)) | (1L << ((long)O));
  sett[25] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)G)) |
    (1L << ((long)U));
  sett[26] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)G)) |
    (1L << ((long)O));
  sett[27] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)U)) |
    (1L << ((long)O));
  sett[28] = (1L << ((long)A)) | (1L << ((long)G)) | (1L << ((long)U)) |
    (1L << ((long)O));
  sett[29] = (1L << ((long)C)) | (1L << ((long)G)) | (1L << ((long)U)) |
    (1L << ((long)O));
  sett[30] = (1L << ((long)A)) | (1L << ((long)C)) | (1L << ((long)G)) |
    (1L << ((long)U)) | (1L << ((long)O));
}  /* setorder */


void mincomp()
{
  /* computes for each site the minimum number of steps
     necessary to accomodate those species already
     in the analysis */
  long i, j, k;
  boolean done;

  for (i = 0; i < (nsites); i++) {
    done = false;
    j = 0;
    while (!done) {
      j++;
      done = true;
      k = 1;
      do {
        if (intree[k - 1])
          done = (done && (treenode[k - 1]->base[i] & sett[j - 1]) != 0);
        k++;
      } while (k <= spp && done);
    }
    if (j == 31)
      necsteps[i] = 4;
    if (j <= 30)
      necsteps[i] = 3;
    if (j <= 25)
      necsteps[i] = 2;
    if (j <= 15)
      necsteps[i] = 1;
    if (j <= 5)
      necsteps[i] = 0;
    necsteps[i] *= weight[i];
  }
}  /* mincomp */


void help()
{
  /* display help information */
  printf("\n\nR  Rearrange a tree by moving a node or group\n");
  printf("#  Show the states of the next site that doesn't fit tree\n");
  printf("+              ... of the next site\n");
  printf("-              ... of the previous site\n");
  printf("S  Show the states of a given site\n");
  printf(".  redisplay the same tree again\n");
  printf("T  Try all possible positions of a node or group\n");
  printf("U  Undo the most recent rearrangement\n");
  printf("W  Write tree to a file\n");
  printf("O  select an Outgroup for the tree\n");
  printf("F  Flip (rotate) branches at a node\n");
  printf("H  Help (this screen)\n");
  printf("?   \"     \"      \"   \n");
  printf("Q  (Quit) Exit from program\n");
  printf("X  Exit from program\n\n\n");
  printf("TO CONTINUE, PRESS ON THE Enter OR Return KEY");
  scanf("%c%*[^\n]", &ch);
  getchar();
  printree(false);
}  /* help */


void rearrange()
{
  long i, j;
  boolean ok1, ok2;
  node *p, *q;

  printf("Remove everything to the right of which node? ");
  inpnum(&i, &ok1);
  ok1 = (ok1 && i >= 1 && i < spp * 2 && i != root->index);
  if (ok1) {
    printf("Add before which node? ");
    inpnum(&j, &ok2);
    ok2 = (ok2 && j >= 1 && j < spp * 2);
    if (ok2) {
      ok2 = (treenode[j - 1] != treenode[treenode[i - 1]->back->index - 1]);
      p = treenode[j - 1];
      while (p != root) {
        ok2 = (ok2 && p != treenode[i - 1]);
        p = treenode[p->back->index - 1];
      }
      if (ok1 && ok2) {
        what = i;
        q = treenode[treenode[i - 1]->back->index - 1];
        if (q->next->back->index == i)
          fromwhere = q->next->next->back->index;
        else
          fromwhere = q->next->back->index;
        towhere = j;
        re_move(&treenode[i - 1], &q);
        add(treenode[j - 1], treenode[i - 1], q);
      }
      lastop = rearr;
    }
  }
  changed = (ok1 && ok2);
  printree(false);
  if (!(ok1 && ok2))
    printf("Not a possible rearrangement.  Try again: \n");
  else {
    oldwritten = written;
    written = false;
  }
}  /* rearrange */


Local Void nextinc()
{
  /* show next incompatible site */
  long disp0;
  boolean done;

  display = true;
  disp0 = dispchar;
  done = false;
  do {
    dispchar++;
    if (dispchar > nsites) {
      dispchar = 1;
      done = (disp0 == 0);
    }
  } while (!(necsteps[dispchar - 1] != numsteps[dispchar - 1] ||
             dispchar == disp0 || done));
  printree(false);
}  /* nextinc */


void nextchar()
{
  /* show next site */
  display = true;
  dispchar++;
  if (dispchar > nsites)
    dispchar = 1;
  printree(false);
}  /* nextchar */


void prevchar()
{
  /* show previous site */
  display = true;
  dispchar--;
  if (dispchar < 1)
    dispchar = nsites;
  printree(false);
}  /* prevchar */


void show()
{
  long i;
  boolean ok;

  do {
    printf("SHOW: (Character number or 0 to see none)? ");
    inpnum(&i, &ok);
    ok = (ok && (i == 0 || (i >= 1 && i <= nsites)));
    if (ok && i != 0) {
      display = true;
      dispchar = i;
    }
    if (ok && i == 0)
      display = false;
  } while (!ok);
  printree(false);
}  /* show */

/* Local variables for addpreorder: */


void tryadd(p,item,nufork,place)
     node *p,**item,**nufork;
     double *place;
{
  /* temporarily adds one fork and one tip to the tree.
     Records scores in ARRAY place */
  add(p, *item, *nufork);
  evaluate(root);
  place[p->index - 1] = -like;
  re_move(item, nufork);
}  /* tryadd */


void addpreorder(p, item_, nufork_,place)
     node *p, *item_, *nufork_;
     double *place;
{
  /* traverses a binary tree, calling PROCEDURE tryadd
     at a node before calling tryadd at its descendants */
  node *item, *nufork;

  item = item_;
  nufork = nufork_;
  if (p == NULL)
    return;
  tryadd(p,&item,&nufork,place);
  if (!p->tip) {
    addpreorder(p->next->back, item,nufork,place);
    addpreorder(p->next->next->back,item,nufork,place);
  }
}  /* addpreorder */


void try()
{
  /* Remove node, try it in all possible places */
  /* Local variables for try: */
  double *place;
  long i, j, oldcompat;
  double current;
  node *q, *dummy, *rute;
  boolean tied, better, ok;

  printf("Try other positions for which node? ");
  inpnum(&i, &ok);
  if (!(ok && i >= 1 && i <= nonodes && i != root->index)) {
    printf("Not a possible choice! ");
    return;
  }
  printf("WAIT ...\n");
  place = (double *)Malloc(nonodes*sizeof(double));
  for (j = 0; j < (nonodes); j++)
    place[j] = -1.0;
  evaluate(root);
  current = -like;
  oldcompat = compatible;
  what = i;
  q = treenode[treenode[i - 1]->back->index - 1];
  if (q->next->back->index == i)
    fromwhere = q->next->next->back->index;
  else
    fromwhere = q->next->back->index;
  rute = root;
  if (root == treenode[treenode[i - 1]->back->index - 1]) {
    if (treenode[treenode[i - 1]->back->index - 1]->next->back == treenode[i - 1])
      rute = treenode[treenode[i - 1]->back->index - 1]->next->next->back;
    else
      rute = treenode[treenode[i - 1]->back->index - 1]->next->back;
  }
  re_move(&treenode[i - 1], &dummy);
  oldleft = wasleft;
  root = rute;
  addpreorder(root, treenode[i - 1], dummy, place);
  wasleft = oldleft;
  restoring = true;
  add(treenode[fromwhere - 1], treenode[what - 1], dummy);
  like = -current;
  compatible = oldcompat;
  restoring = false;
  better = false;
  printf("       BETTER: ");
  for (j = 1; j <= (nonodes); j++) {
    if (place[j - 1] < current && place[j - 1] >= 0.0) {
      printf("%3ld:%6.2f", j, place[j - 1]);
      better = true;
    }
  }
  if (!better)
    printf(" NONE");
  printf("\n       TIED:    ");
  tied = false;
  for (j = 1; j <= (nonodes); j++) {
    if (fabs(place[j - 1] - current) < 1.0e-6 && j != fromwhere) {
      if (j < 10)
        printf("%2ld", j);
      else
        printf("%3ld", j);
      tied = true;
    }
  }
  if (tied)
    printf(":%6.2f\n", current);
  else
    printf("NONE\n");
  changed = true;
  free(place);
}  /* try */


void undo()
{
  /* restore to tree before last rearrangement */
  long temp;
  boolean btemp;
  node *q;

  switch (lastop) {

  case rearr:
    restoring = true;
    oldleft = wasleft;
    re_move(&treenode[what - 1], &q);
    btemp = wasleft;
    wasleft = oldleft;
    add(treenode[fromwhere - 1], treenode[what - 1],q);
    wasleft = btemp;
    restoring = false;
    temp = fromwhere;
    fromwhere = towhere;
    towhere = temp;
    changed = true;
    break;

  case flipp:
    q = treenode[atwhat - 1]->next->back;
    treenode[atwhat - 1]->next->back =treenode[atwhat - 1]->next->next->back;
    treenode[atwhat - 1]->next->next->back = q;
    treenode[atwhat - 1]->next->back->back = treenode[atwhat - 1]->next;
    treenode[atwhat - 1]->next->next->back->back =
      treenode[atwhat - 1]->next->next;
    break;

  case reroott:
    restoring = true;
    temp =oldoutgrno;
    oldoutgrno = outgrno;
    outgrno = temp;
    reroot(treenode[outgrno - 1]);
    restoring = false;
    break;

  case none:
    /* blank case */
    break;
  }
  printree(false);
  if (lastop == none) {
    printf("No operation to undo! ");
    return;
  }
  btemp = oldwritten;
  oldwritten = written;
  written = btemp;
}  /* undo */


void treeout(p, lengths)
     node *p;
     boolean lengths;
{
  /* write out file with representation of final tree */
  long i, n;
  Char c;

  if (p->tip) {
    n = 0;
    for (i = 1; i <= nmlngth; i++) {
      if (nayme[p->index - 1][i - 1] != ' ')
        n = i;
    }
    for (i = 0; i < n; i++) {
      c = nayme[p->index - 1][i];
      if (c == ' ')
        c = '_';
      putc(c, treefile);
    }
    col += n;
  } else {
    putc('(', treefile);
    col++;
    treeout(p->next->back, lengths);
    putc(',', treefile);
    col++;
    if (col > 65) {
      putc('\n', treefile);
      col = 0;
    }
    treeout(p->next->next->back, lengths);
    putc(')', treefile);
    col++;
  }
  if (p == root)
    fprintf(treefile, ";\n");
  else if (lengths)
    fprintf(treefile, ":%6.5f", p->tyme-treenode[p->back->index-1]->tyme);
}  /* treeout */


void treewrite(done, lengths)
     boolean done;
     boolean lengths;
{
  /* write out tree to a file */

  if (waswritten) {
    printf("\nTree file already was open.\n");
    printf("   A   Add to this tree to tree file\n");
    printf("   R   Replace tree file contents by this tree\n");
    printf("   F   Write out tree to a different tree file\n");
    printf("   N   Do Not write out this tree\n");
    do {
      printf("Which should we do? ");
      scanf("%c%*[^\n]", &ch);
      getchar();
      uppercase(&ch);
    } while (ch != 'A' && ch != 'R' && ch != 'N'  && ch != 'F');
  }
  if (!waswritten || (ch == 'F')) {
    trfilename[0] = '\0';
    while (trfilename[0] =='\0'){
      printf("Please enter a tree file name>");
      fgets(trfilename, 100, INPUT);
    }
  }
  if (ch == 'R' || ch == 'A' || !waswritten){
    openfile(&treefile, trfilename, "tree file",
             (waswritten && (ch == 'A')) ? "a" : "w");
  }
  if (!done)
    printree(false);
  if (waswritten && ch != 'A' && ch != 'R')
    return;
  col = 0;
  treeout(root, lengths);
  printf("\nTree written to file %s\n\n", trfilename);
  waswritten = true;
  written = true;
  FClose(treefile);
#ifdef MAC
  fixmacfile(trfilename);
#endif
}
/* treewrite */


void flip()
{
  /* flip at a node left-right */
  long i;
  boolean ok;
  node *p;

  printf("Flip branches at which node? ");
  inpnum(&i, &ok);
  ok = (ok && i > spp && i <= nonodes);
  if (ok) {
    p = treenode[i - 1]->next->back;
    treenode[i - 1]->next->back = treenode[i - 1]->next->next->back;
    treenode[i - 1]->next->next->back = p;
    treenode[i - 1]->next->back->back = treenode[i - 1]->next;
    treenode[i - 1]->next->next->back->back = treenode[i - 1]->next->next;
    atwhat = i;
    lastop = flipp;
  }
  printree(false);
  if (ok) {
    oldwritten = written;
    written = false;
    return;
  }
  if (i >= 1 && i <= spp)
    printf("Can't flip there. ");
  else
    printf("No such node. ");
}  /* flip */


void changeoutgroup()
{
  long i;
  boolean ok;

  oldoutgrno = outgrno;
  do {
    printf("Which node should be the new outgroup? ");
    inpnum(&i, &ok);
    ok = (ok && intree[i - 1] && i >= 1 && i <= nonodes &&
          i != root->index);
    if (ok)
      outgrno = i;
  } while (!ok);
  if (intree[outgrno - 1])
    reroot(treenode[outgrno - 1]);
  changed = true;
  lastop = reroott;
  printree(false);
  oldwritten = written;
  written = false;
}  /* changeoutgroup */


void redisplay()
{
  /* Local variables for redisplay: */
  boolean done = false;
  do {
    printf("NEXT? (Options: R # + - S . T U W O F C H ? X Q) ");
    printf("(H or ? for Help) ");
    scanf("%c%*[^\n]", &ch);
    getchar();
    uppercase(&ch);
    if (strchr("FORSTUXQ+#-.WH?",ch)){
      switch (ch) {
        
      case 'R':
        rearrange();
        break;
        
      case '#':
        nextinc();
        break;
        
      case '+':
        nextchar();
        break;
        
      case '-':
        prevchar();
        break;
        
      case 'S':
        show();
        break;
        
      case '.':
        printree(false);
        break;
        
      case 'T':
        try();
        break;
        
      case 'U':
        undo();
        break;
        
      case 'W':
        treewrite(done, false);
        break;
        
      case 'O':
        changeoutgroup();
        break;
        
      case 'F':
        flip();
        break;
        
      case 'H':
        help();
        break;
        
      case '?':
        help();
        break;
        
      case 'X':
        done = true;
        break;
        
      case 'Q':
        done = true;
        break;
      }
    }
  } while (!done);
  if (written)
    return;
  do {
    printf("Do you want to write out the tree to a file? (Y or N) ");
    scanf("%c%*[^\n]", &ch);
    getchar();
    if (ch == 'Y' || ch == 'y')
      treewrite(done, false);
  } while (ch != 'Y' && ch != 'y' && ch != 'N' && ch != 'n');
}  /* redisplay */


void treeconstruct()
{
  long i;
  node *p;

  /* constructs a binary tree from the pointers in treenode. */

  restoring = false;
  subtree = false;
  display = false;
  dispchar = 0;
  earlytree = true;
  waswritten = false;
  for (i = 0; i < spp; i++)
    treenode[i]->back = NULL;
  for (i = spp; i < nonodes; i++) {
    p = treenode[i];
    do {
      p->back = NULL;
      p = p->next;
    } while (p != treenode[i]);
  }
  buildtree();
  printf("\nComputing steps needed for compatibility in sites ...\n\n");
  setorder();
  mincomp();
  newtree = true;
  earlytree = false;
  printree(false);
  bestyet = -like;
  gotlike = -like;
  lastop = none;
  newtree = false;
  written = false;
  redisplay();
}  /* treeconstruct */


void simmain()
{  /* main program -- pretty self explanatory */
  Char ch;

  outfile = NULL;
  garbage = NULL;
  inputtreedata();
  setuptree();
  nayme = (naym *)Malloc(spp*sizeof(naym));
  simtree();
  printree(true);
  waswritten = false;
  do {
    printf("\nWrite this tree out to a tree file (Y or N)\n");
    scanf("%c%*[^\n]", &ch);
    getchar();
    uppercase(&ch);
    if (ch == 'Y') {
/* debug
      printf("please type name of that tree file: ");
      fgets(trfilename, 100, INPUT);
      openfile(&treefile, trfilename, "tree file", "w");
debug */
      treewrite(true, true);
    }
  } while ((ch  != 'Y') && (ch != 'N'));
} /* simmain */


int main(argc, argv)
     int argc;
     Char *argv[];

{ /* Interactive simulation of trees and data, and DNA parsimony */
  /* simulates a tree, or reads one in
     simulates data, or reads it in
     query the user who rearranges the tree */
#ifdef MAC
/*  macsetup("Dnatree","");*/
  argv[0] = "Dnatree";
#endif
  strcpy(progname,argv[0]);
  strcpy(outfilename,OUTFILE);
  strcpy(trfilename,TREEFILE);
  garbage = NULL;
  simmain();
  treednamain();
  doinput();
  treeconstruct();
  if (waswritten) {
    FClose(treefile);
#ifdef MAC
    fixmacfile(trfilename);
#endif
  }
#ifdef MAC
  fixmacfile(outfilename);
#endif
  printf("\nDone.\n\n");
  return(0);
  exit(0);
}  /* Interactive DNA evolution */


int eof(f)
     FILE *f;
{
  register int ch;

  if (feof(f))
    return 1;
  if (f == stdin)
    return 0;
  ch = getc(f);
  if (ch == EOF)
    return 1;
  ungetc(ch, f);
  return 0;
}


int eoln(f)
     FILE *f;
{
  register int ch;

  ch = getc(f);
  if (ch == EOF)
    return 1;
  ungetc(ch, f);
  return (ch == '\n');
}


void memerror()
{
printf("Error allocating memory\n");
exit(-1);
}


MALLOCRETURN *mymalloc(x)
long x;
{
MALLOCRETURN *mem;

mem = (MALLOCRETURN *)calloc(1, x);
if (!mem) {
  memerror();
  return((MALLOCRETURN *)(1));
}
else
  return (MALLOCRETURN *)mem;
}