#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include <math.h>

#define VERSION "0.8"
#define DATAFILE_NAME "outfile"
#define NUMCHARS 1005
#define NMLTH 50

#define MAX_GAIN 10.0
#define MAX_LOSS 10.0
#define DEF_INC 0.10
#define DIST 800  /* amount of storage to allocate for distribution array */
#define NSTATES 2

/* Node structure */
struct node {
  struct node *left, *right, *back;
  int br_flag;  /* Indicates how many branches are attached */
  char name[NMLTH];
  char charstate[NUMCHARS];
  double downsteps1[NUMCHARS];
  double downsteps0[NUMCHARS];
  int node_id;
  double decay;
  int vi; /* visit index for writetree */
};

/* Tree storage structure */
struct treestruct {
  struct treestruct *next;
  struct node *root;
  char name[NMLTH];
  int id;
};

/* Floating point linked list structure */
struct fplist {
  double value;
  struct fplist *next, *back;
  int id;
};

/* Char linked list structure */
struct charlist {
	char state;
	struct charlist *next, *back;
	int id;
};

/********************************
 * Global variable declarations *
 ********************************/
FILE *infile;	/* input file */
FILE *outfile;	/* output file for plot data */
/* FILE *probfile;	  */
FILE *plotfile;	 /* plot file */
FILE *treefile;

char *infilename;  /* pointer to the infile name string */
char *outfilename;  /* pointer to the outfile name string */
char *tfname;  /* pointer to treefile name string */
char *pfname;  /* pointer to gnuplot commandfile */

double maxgain = MAX_GAIN;
double maxloss = MAX_LOSS;
double definc = DEF_INC;

int numtrees = 0;
int numchars = 0;
int treenum;
struct treestruct *tree_list;  /* pointer to the start of the tree list */
struct treestruct *tree;  /* pointer to the current tree */
int character;
struct node *root;  /* pointer to the root of the tree */
int num_taxa = 0;  /* number of tips on current tree */

double g_steps, l_steps;  /* values for step matrix */
int gainloss[5];
char tempstring[100];  /* crap command line string */
char tempname[15];  /* called in main for file opening */
char kar;

int distrib[DIST];
struct fplist *inflectlist;

double pi_gain, pi_loss;

/***********************************
 * Function prototype declarations *
 ***********************************/
void graph_char(int c, int t);
void sil_tul(void);
int countchange(struct node *startnode, int indch, char indstate, int depch, char startstate, char endstate);

struct node * newnode(struct node *input_node, char dir);

int count_trees(void);
struct treestruct * buildtree(void);
struct node * makenode(struct node *input_node);
struct node * maketree(void);
void getindex(struct node *start_node);
int count_taxa(struct node *start_node);
struct node * addstates(struct node *rootnode);
struct node * putdata(struct node *input_node, char name[NMLTH], char data[NUMCHARS]);

struct node * downpass(struct node *start_node, int c);
struct node * uppass(struct node *start_node, char state, int c);
double d_numsteps(struct node *input_node, char state, int c);
struct node * clearpass(struct node *start_node, int c);
int * countgl(struct node *start_node, int c);
void getstates(struct node * start_node, int c, struct charlist **listptr);

void banner(void);
void options(void);
void exec_file(void);
void list_trees(void);
struct treestruct * set_tree(int num);
void setvalues(void);
int getvalue(void);
double reconstruct(int c, struct treestruct *t, char flag);
void combine(int *a, int position, int nstates);

void putdecay(struct node **node, int nodeid, double value);
void getdecay(struct node *node);

void writetree(struct node **node);

/******************************
 *      the main program      *
 ******************************/
main(int argc, char *argv[])
{
  char *ctemp, *ttemp;
  char tempstring1[100], tempstring2[100], tempstring3[100];
  int c = 1, t = 1;

  banner();  /* print opening header */

  /* parse arguments */
  if (argc > 1) {
    infilename = &argv[(--argc)][0];

    while ((argc > 1) && (argv[1][0] == '-')) {
      switch (argv[1][1]) {
      case 'c':
	ctemp = &argv[1][2];
	if (*ctemp == 'a')
	  c = -2;
	else {
	  c = atoi(ctemp);
	  c--;
	}
	break;

      case 't':
	ttemp = &argv[1][2];
	if (*ttemp == 'a')
	  t = -1;
	else
	  t = atoi(ttemp);
	break;

      case 'm':
	if (argv[1][2] == 'g')
	  maxgain = atof(&argv[1][3]);
	else if (argv[1][2] == 'l')
	  maxloss = atof(&argv[1][3]);
	else
	  fprintf(stderr, "Bad option %s\n", argv[1]);
	break;

      case 'i':
	definc = atof(&argv[1][2]);
	break;

      default:
	fprintf(stderr, "Bad option %s\n", argv[1]);
      }
      argv++;
      argc--;
    }

    strcpy(tempstring1, infilename);
    strcat(tempstring1, ".c");
    strcat(tempstring1, ctemp);
    strcat(tempstring1, ".t");
    strcat(tempstring1, ttemp);

    strcpy(tempstring2, tempstring1);
    strcpy(tempstring3, tempstring1);

    outfilename = tempstring1;
    tfname = tempstring2;
    pfname = tempstring3;

    strcat(outfilename, ".data");
    strcat(tfname, ".tree");
    strcat(pfname, ".plot");

    exec_file();

    graph_char(c, t);
  }

  else {
    printf("No files opened.\n");
  }

  /* close the tree file */
  fclose(infile);

  return(0);
}

/**************************************************
 * newnode -- create and link a new node in the   *
 *            tree	                          *
 * Parameters -input_node - the current node      *
 *            -direction of intended branch       *
 * Returns: pointer to the current node           *
 **************************************************/
struct node * newnode(struct node *input_node, char dir)
{
  /* Create temporary node with branches pointing to nothing */
  struct node *temp_node;
  temp_node = (struct node *)malloc(sizeof(struct node));
  temp_node->br_flag = 0;

  /* attach new node to current node */
  if (dir == 'l') {
    input_node->left = temp_node;
    input_node->br_flag++;
  }
  else if (dir == 'r') {
    input_node->right = temp_node;
    input_node->br_flag++;
  }
  else {
    printf("Error in function call: newnode\n");
    exit(1);
  }
  temp_node->back = input_node;

  /*	input_node = temp_node;*/
  return(temp_node);
}

/****************************************************************
 * graph_char -- read the input file containing trees and data	*
 *   and graph changes over the step matrix range	        *
 * parameters: index int to a character and a tree              *
 * returns: void                                                *
 * subs called: buildtree, among others                         *
 ****************************************************************/
void graph_char(int c, int t)
{
  int i = 0;  /* index into the distribution array */
  char memflag = 'n';
  struct fplist *ip_list;
  int greater = 0, less = 0, equal = 0, total = 0;
  double testval;

  for (i = 0; i <= 5; i++) {
    gainloss[i] = 0;
  }

  g_steps = 1.0;
  l_steps = 1.0;

  inflectlist = (struct fplist *)malloc(sizeof(struct fplist));
  ip_list = inflectlist;

  /* set all values of distribution array to zero */
  for (i = 0; i < DIST; i++) {
    distrib[i] = 0;
  }

  /* open the data file to be written */
  outfile = fopen(outfilename, "w");
  fprintf(outfile, "#Datafile generated by crap v.%s\n", VERSION);
  fprintf(outfile, "#C(g) - C(l), gains, losses\n\n");
  /*	probfile = fopen(probfilename, "w"); */

  if (t == -1) {
    if (c == -2 | c >= numchars) {
      printf("Character specification is invalid\n");
      return;
    }
    else {
      memflag = 'y';
      tree = tree_list;
      while (1) {
	ip_list->id = (tree->id - 1);
	ip_list->value = reconstruct(c, tree, 'a');
	printf("%d %2.2f\n", ip_list->id, ip_list->value);
	if (! tree->next)
	  break;
	else {
	  tree = tree->next;
	  ip_list->next = (struct fplist *)malloc(sizeof(struct fplist));
	  ip_list = ip_list->next;
	}
      }
    }
  }

  else if (c == -2) {
    if (t == -1 | t > numtrees) {
      printf("Tree specification is invalid\n");
      return;
    }
    else {
      memflag = 'y';
      c = 0;
      tree = set_tree(t);
      while (1) {
	ip_list->id = c;
	ip_list->value = reconstruct(c, tree, 'a');
	if (c == (numchars - 1))
	  break;
	else {
	  c++;
	  ip_list->next = (struct fplist *)malloc(sizeof(struct fplist));
	  ip_list = ip_list->next;
	  tree = set_tree(t);
	}
      }
    }
  }

  else {
    tree = set_tree(t);
    ip_list->value = reconstruct(c, tree, 's');
  }

  fclose(outfile);
  /*	fclose(probfile); */

  plotfile = fopen(pfname, "w");

  fprintf(plotfile, "# Gnuplot commandfile generated by CRAP-%s\n", VERSION);
  fprintf(plotfile, "set xlabel '(cost of gain)-(cost of loss)'\n");
  fprintf(plotfile, "set ylabel 'number of changes'\n");
  fprintf(plotfile, "set title '%s'\n", infilename);
  fprintf(plotfile, "plot '%s' title '0->1' with linespoints, '%s' using '%cf%c*f%cf' title '1->0' with linespoints\npause -1\n", outfilename, outfilename, '%', '%', '%');

  fclose(plotfile);
  /*	system("gnuplot crap_plotfile"); */

  printf("Data written to %s.\n", outfilename);

  /* print the distribution of inflection points */
  printf("Horizontal distribution of inflection points:\n\n");
  i = 0;
  while (i < DIST) {
    if (distrib[i] > 0) {
      double j = ((i - (DIST/2 - 1))*0.1);
      int x = 0;

      printf("%2.1f: ", j);

      while (1) {
	if (x >= distrib[i]) {
	  printf("\n");
	  break;
	}
	else
	  printf("*");

	x++;
      }
    }
    else
      ;
    i++;
  }

  /* print out inflection point values */
  printf("Inflection point values:\n\n");
  ip_list = inflectlist;
  printf("%d %2.2f\n", ip_list->id, ip_list->value);
  while (memflag == 'y') {
    printf("%2.2f ", ip_list->value);
    if (ip_list->id == 0)
      testval = ip_list->value;
    else {
      if (ip_list->value > testval)
	greater++;
      else if (ip_list->value < testval)
	less++;
      else if (ip_list->value == testval)
	equal++;
      else
	;
    }
    total++;
    if (ip_list->id == (numchars - 1))
      break;
    else
      ip_list = ip_list->next;
  }
  printf("\n\n");
  printf("greater: %d, less: %d, equal: %d, total: %d\n", greater, less, equal, total);

  return;
}

/****************************************************
 * buildtree -- read in a parenthetical tree from   *
 *   an input file and create it using pointers     *
 * parameters - none                                *
 * subs called: makenode                            *
 * returns - pointer to root node                   *
 ****************************************************/
struct treestruct * buildtree(void)
{

  /* Local variable declarations */
  struct treestruct *first_tree, *current_tree, *temp_tree;
  char string[25];
  char nextline[200];
  char ch;
  int i;

  int index = 0;

  rewind(infile);

  /*
    find the TREES block
    */
  while (1) {
    fgets(nextline, sizeof(nextline), infile);
    ch = fgetc(infile);
    if (ch == EOF) {
      printf("No trees block found\n");
      exit(1);
    }
    else if (strstr(nextline, "GIN TREES;")) {
      printf("reading trees block...\n");
      break;
    }
    else {
      ;
    }
  }

  /*
    find the tree
    */
  while (1) {
    ch = fgetc(infile);
    if (ch == EOF) break;
    if (ch == 'T') ch = fgetc(infile);
    if (ch == 'R') ch = fgetc(infile);
    if (ch == 'E') ch = fgetc(infile);
    if (ch == 'E') {

      if (index == 0) {
	first_tree = (struct treestruct *)malloc(sizeof(struct treestruct));
	current_tree = first_tree;
	index++;
	current_tree->id = index;
      }
      else {
	temp_tree = (struct treestruct *)malloc(sizeof(struct treestruct));
	current_tree->next = temp_tree;
	current_tree = current_tree->next;
	index++;
	current_tree->id = index;
      }

      /* read along to the name of the tree */
      ch = fgetc(infile);
      while (1) {
	ch = fgetc(infile);
	if (isalnum(ch)) {
	  break;
	}
      }

      /* get the tree name */
      i = 0;
      string[i] = ch;
      while (1) {
	ch = fgetc(infile);
	i++;
	if (!isalnum(ch) && ch != '.') {
	  string[i] = '\0';
	  break;
	}
	string[i] = ch;
      }

      current_tree->root = maketree();

      strcpy(current_tree->name, string);
      printf("Tree name is %s\n", current_tree->name);
    }
  }

  printf("%d trees now in memory\n", index);

  return(first_tree);
}

/**************************************************
 * makenode -- create and link the next node in   *
 *   the tree                                     *
 * Parameters: input_node - the current node      *
 * Returns: pointer to the current node           *
 **************************************************/
struct node * makenode(struct node *input_node)
{
  /* Create temporary node with branches pointing to nothing */
  struct node *temp_node;
  temp_node = (struct node *)malloc(sizeof(struct node));
  temp_node->br_flag = 0;

  /* attach new node to current node */
  if (input_node->br_flag == 0) {
    input_node->left = temp_node;
    input_node->br_flag++;
  }
  else if (input_node->br_flag == 1) {
    input_node->right = temp_node;
    input_node->br_flag++;
  }
  else {
    printf("F: Error: trifurcation not allowed\n");
    exit(1);
  }
  temp_node->back = input_node;

  /*	input_node = temp_node;*/
  return(temp_node);
}

/****************************************************
 * getindex -- recurse through the tree and print   *
 *   the numbers of the nodes                       *
 * parameters - pointer to node to start at         *
 * returns - none                                   *
 ****************************************************/
void getindex(struct node *start_node)
{
  if (strcmp(start_node->name, "-") == 0) {
    printf("%d ", start_node->node_id);
    getindex(start_node->left);
    getindex(start_node->right);
  }
  else {
    printf("%d* ", start_node->node_id);
  }
  return;
}

/***************************************************
 * count_taxa -- count the descendant tips of a    *
 *   given node                                    *
 * parameters - pointer to node to start at        *
 * returns - number of taxa above the node         *
 ***************************************************/
int count_taxa(struct node *start_node) 
{
  int count = 0;

  if (strcmp(start_node->name, "-") == 0) {
    count_taxa(start_node->left);
    count_taxa(start_node->right);
  }
  else {
    num_taxa++;
  }

  count = num_taxa;
  return(count);
}

/*********************************************************
 * addstates -- put character state data into tip nodes  *
 * parameters: pointer to root of tree                   *
 * returns: same                                         *
 *********************************************************/
struct node * addstates(struct node *rootnode)
{
  char nextline[1500];
  char ch;
  char *tempname;
  char *tempchars;
  char *tmp_ptr;

  rewind(infile);
  while (1) {
    fgets(nextline, sizeof(nextline), infile);
    ch = fgetc(infile);
    if (ch == EOF) {
      printf("No data block found\n");
      exit(1);
    }
    if (strstr(nextline, "ATRIX")) {
      /*			printf("reading data block...");*/

      while (1) {
	fgets(nextline, sizeof(nextline), infile);
	if (strchr(nextline, ';')) {
	  /*					printf(" done\n");*/
	  break;
	}

	else if (strcmp(nextline, "\n") == 0) {
	  ;
	}

	else if (nextline[0] == '[') {
	  ;
	}

	else {
	  nextline[strlen(nextline) - 1] = '\0';
	  tempname = nextline;

	  tmp_ptr = strchr(nextline, ' ');

	  *tmp_ptr = '\0';
	  tmp_ptr++;

	  tempchars = strrchr(tmp_ptr, ' ');
	  *tempchars = '\0';
	  tempchars++;

	  rootnode = putdata(rootnode, tempname, tempchars);
	}
      }
    }

    else if (strstr(nextline, "ENDBLOCK;")) {
      break;
    }

    else {
      ;
    }
  }

  return (rootnode);
}


struct node * putdata(struct node *input_node, char name[NMLTH], char data[NUMCHARS])
{
  if (input_node->br_flag == 2) {
    input_node->left = putdata(input_node->left, name, data);
    input_node->right = putdata(input_node->right, name, data);
  }

  else if (strcmp(input_node->name, name) == 0) {
    strcpy(input_node->charstate, data);
  }

  else {
    ;
  }

  return (input_node);
}

/*****************************************************************
 * downpass -- recurse down through the tree, calling d_numsteps *
 * parameters - pointer to start node                            *
 *            - character of interest                            *
 * returns - pointer to same node                                *
 *****************************************************************/
struct node * downpass(struct node *start_node, int c)
{

  if (start_node->br_flag == 2) {
    start_node->left = downpass(start_node->left, c);
    start_node->right = downpass(start_node->right, c);

    start_node->downsteps1[c] = d_numsteps(start_node, '1', c);
    start_node->downsteps0[c] = d_numsteps(start_node, '0', c);

    return (start_node);
  }

  else {
    return (start_node);
  }
}

/****************************************************
 * uppass -- recurse up	through the tree            *
 * parameters - pointer to root of clade to pass    *
 *            - state of ancestral node             *
 *            - character of interest               *
 * returns - pointer to same node                   *
 ****************************************************/
struct node * uppass(struct node *start_node, char state, int c)
{
  char ifanc1, ifanc0;

  if (start_node->br_flag != 2) {
    return (start_node);
  }

  else if (state == 'r') {
    if (start_node->downsteps1[c] < start_node->downsteps0[c]) {
      start_node->charstate[c] = '1';
    }
    else if (start_node->downsteps1[c] > start_node->downsteps0[c]) {
      start_node->charstate[c] = '0';
    }
    else {
      start_node->charstate[c] = 'e';
    }

    start_node->left = uppass(start_node->left, start_node->charstate[c], c);
    start_node->right = uppass(start_node->right, start_node->charstate[c], c);

    return(start_node);
  }

  else if (state == '1') {

    if (start_node->downsteps1[c] < (start_node->downsteps0[c] + l_steps)) {
      start_node->charstate[c] = '1';
    }
    else if (start_node->downsteps1[c] > (start_node->downsteps0[c] + l_steps)) {
      start_node->charstate[c] = '0';
    }
    else {
      start_node->charstate[c] = 'e';
    }

    start_node->left = uppass(start_node->left, start_node->charstate[c], c);
    start_node->right = uppass(start_node->right, start_node->charstate[c], c);

    return(start_node);
  }

  else if (state == '0') {

    if ((start_node->downsteps1[c] + g_steps) < start_node->downsteps0[c]) {
      start_node->charstate[c] = '1';
    }
    else if ((start_node->downsteps1[c] + g_steps) > start_node->downsteps0[c]) {
      start_node->charstate[c] = '0';
    }
    else {
      start_node->charstate[c] = 'e';
    }

    start_node->left = uppass(start_node->left, start_node->charstate[c], c);
    start_node->right = uppass(start_node->right, start_node->charstate[c], c);

    return(start_node);
  }

  else if (state == 'e') {

    if (start_node->downsteps1[c] < (start_node->downsteps0[c] + l_steps)) {
      ifanc1 = '1';
    }
    else if (start_node->downsteps1[c] > (start_node->downsteps0[c] + l_steps)) {
      ifanc1 = '0';
    }
    else {
      ifanc1 = 'e';
    }

    if ((start_node->downsteps1[c] + g_steps) < start_node->downsteps0[c]) {
      ifanc0 = '1';
    }
    else if ((start_node->downsteps1[c] + g_steps) > start_node->downsteps0[c]) {
      ifanc0 = '0';
    }
    else {
      ifanc0 = 'e';
    }

    if (ifanc1 == ifanc0) {
      start_node->charstate[c] = ifanc1;
    }
    else if (ifanc1 == 'e') {
      start_node->charstate[c] = ifanc0;
    }
    else if (ifanc0 == 'e') {
      start_node->charstate[c] = ifanc1;
    }
    else {
      start_node->charstate[c] = 'e';
    }

    start_node->left = uppass(start_node->left, start_node->charstate[c], c);
    start_node->right = uppass(start_node->right, start_node->charstate[c], c);
			
    return(start_node);
  }

  else {
    printf("Error in function call: uppass (state = %c, gsteps=%2.1f, lsteps=%2.1f)\n", state, g_steps, l_steps);
  }
}

/*************************************************************
 * d_numsteps -- count the number of downsteps in a clade    *
 * parameters - root node of clade of interest               *
 *            - state assigned to input node                 *
 *            - character of interest                        *
 * returns - float number of steps                           *
 *************************************************************/
double d_numsteps(struct node *input_node, char state, int c)
{
  double left, right, total;
  double l_ds1, l_ds0, r_ds1, r_ds0;
	
  /*-----------------------------*
   |get info from left descendant|
   *-----------------------------*/
  if (input_node->left->br_flag == 2) {

    l_ds1 = input_node->left->downsteps1[c];
    l_ds0 = input_node->left->downsteps0[c];

    if (state == '1') {
      if (l_ds1 <= (l_ds0 + l_steps)) {
	left = l_ds1;
      }
      else {
	left = (l_ds0 + l_steps);
      }
    }
	
    else if (state == '0') {
      if (l_ds0 <= (l_ds1 + g_steps)) {
	left = l_ds0;
      }
      else {
	left = (l_ds1 + g_steps);
      }
    }

    else {
      printf("Error in function call: d_numsteps\n");
      exit(1);
    }
  }

  else if (input_node->left->charstate[c] == '1') {
    if (state == '1') {
      left = 0.0;
    }
    else if (state == '0') {
      left = g_steps;
    }
    else {
      printf("Error in function call: d_numsteps\n");
      exit(1);
    }
  }

  else if (input_node->left->charstate[c] == '0') {
    if (state == '1') {
      left = l_steps;
    }
    else if (state == '0') {
      left = 0.0;
    }
    else {
      printf("Error in function call: d_numsteps\n");
      exit(1);
    }
  }

  else if (input_node->left->charstate[c] == '?' | input_node->left->charstate[c] == '-' | input_node->left->charstate[c] == '2') {
    left = 0.0;
  }

  else {
    printf("Error in downpass: state is %c, %d\n", input_node->left->charstate[c], input_node->left->charstate[c]);
    exit(1);
  }
	
   /*------------------------------*
    |get info from right descendant|
    *------------------------------*/
  if (input_node->right->br_flag == 2) {

    r_ds1 = input_node->right->downsteps1[c];
    r_ds0 = input_node->right->downsteps0[c];

    if (state == '1') {
      if (r_ds1 <= (r_ds0 + l_steps)) {
	right = r_ds1;
      }
      else {
	right = (r_ds0 + l_steps);
      }
    }

    else if (state == '0') {
      if (r_ds0 <= (r_ds1 + g_steps)) {
	right = r_ds0;
      }
      else {
	right = (r_ds1 + g_steps);
      }
    }
	
    else {
      printf("Error in function call: d_numsteps\n");
      exit(1);
    }
  }

  else if (input_node->right->charstate[c] == '1') {
    if (state == '1') {
      right = 0.0;
    }
    else if (state == '0') {
      right = g_steps;
    }
    else {
      printf("Error in function call: d_numsteps\n");
    }
  }

  else if (input_node->right->charstate[c] == '0') {
    if (state == '1') {
      right = l_steps;
    }
    else if (state == '0') {
      right = 0.0;
    }
    else {
      printf("Error in function call: d_numsteps\n");
      exit(1);
    }
  }

  else if (input_node->right->charstate[c] == '?' | input_node->right->charstate[c] == '-' | input_node->right->charstate[c] == '2') {
    right = 0.0;
  }

  else {
    printf("Error in downpass: state is %c\n", input_node->right->charstate[c]);
    exit(1);
  }

  total = (left + right);
  return (total);
}

/*********************************************************
 * clearpass -- clear reconstructed states from internal *
 *  nodes between iterations                             *
 *********************************************************/
struct node * clearpass(struct node *start_node, int c)
{
  if (start_node->br_flag == 2) {
    start_node->left = clearpass(start_node->left, c);
    start_node->right = clearpass(start_node->right, c);

    start_node->charstate[c] = '-';
    return(start_node);
  }

  else {
    ;
  }
}

/*****************************************************
 * countgl -- count the number of changes from 0->1  *
 *  and 1->0, and the number of branches starting in *
 *  0, 1, or e                                       *
 *****************************************************/
int *countgl(struct node *start_node, int c)
{
  char left, right, current;
  int *tempcount;

  if (start_node->br_flag == 2) {
    left = start_node->left->charstate[c];
    right = start_node->right->charstate[c];
    current = start_node->charstate[c];

    if (current == '0') {
      gainloss[1]++;
      if (left == '1' | left == '2') {
	gainloss[0]++;
      }
      else {
	;
      }
      if (right == '1' | right == '2') {
	gainloss[0]++;
      }
      else {
	;
      }
    }

    else if (current == '1') {
      gainloss[3]++;
      if (left == '0' | left == '2') {
	gainloss[2]++;
      }
      else {
	;
      }
      if (right == '0' | right == '2') {
	gainloss[2]++;
      }
      else {
	;
      }
    }

    else if (current == 'e') {
      gainloss[4]++;
      if (left == '2' | right == '2') {
	printf("2 polymorphic sister taxa?\n");
	/*exit(1);*/
      }
    }

    else {
      ;
    }

    countgl(start_node->left, c);
    countgl(start_node->right, c);

    tempcount = &gainloss[0];
    return(tempcount);
  }

  else {
    ;
  }
}

/*
getstates --
*/
void getstates(struct node *start_node, int c, struct charlist **listptr)
{
  struct charlist *temp_ptr;

  if (start_node->br_flag == 2) { /* internal node */
    if ((*listptr) == NULL) {
      (*listptr) = (struct charlist *)malloc(sizeof(struct charlist));
    }
    else {
      (*listptr)->next = (struct charlist *)malloc(sizeof(struct charlist));
      temp_ptr = *listptr;
      *listptr = (*listptr)->next;
      (*listptr)->back = temp_ptr;
    }
    (*listptr)->state = start_node->charstate[c];
    (*listptr)->id = start_node->node_id;
    getstates(start_node->left, c, &(*listptr));
    getstates(start_node->right, c, &(*listptr));
  }
  else { /* is a tip */
    /*
      printf("%c* ", start_node->charstate[c]);
      */
  }

  return;
}

void banner(void)
{
  system("clear");

  printf("\t************************************************************\n");
  printf("\t  CRAP - a Character-state Reconstruction Analysis Program\n");
  printf("\t                         Version %s\n", VERSION);
  printf("\t               by Rick Ree <rree@oeb.harvard.edu>\n");
  printf("\t                  -- developed using LINUX --\n");
  printf("\t************************************************************\n");

  printf("\n");
}

int count_trees(void)
{
  int count = 0;
  char line[100];
  char ch;

  rewind(infile);

  while (1) {
    fgets(line, sizeof(line), infile);
    ch = fgetc(infile);
    if (ch == EOF) {
      break;
    }
    if (strstr(line, "TREE ")) {
      count++;
    }
    if (strstr(line, "IMENSIONS")) {
      char *tmpstring;
      int *nchars;

      tmpstring = (strrchr(line, '='));
      tmpstring++;
      sscanf(tmpstring, "%d", &nchars);
      numchars = (int) nchars;
    }
  }

  return (count);
}


struct node * maketree(void)
{
  struct node *current_node, *root_node;
  char ch; /* used to read the input file */
  int maderoot = 0;  /* Equal to 1 if root node has been made */
  int nodecount = 0;  /* number of nodes created, num assigned to each node */
  int i = 0; /* index into the name string */
  int eon = 0; /* end of name signal */
  char tempname[NMLTH];

  /* find the beginning of the tree */
  while (1) {
    ch = fgetc(infile);
    if (ch == '(') {
      eon = 1;
      break;
    }
    else {
      ;
    }
  }

  /* read the tree (in parenthetical notation) */
  while (1) {
    if (eon != 1) {
      ch = fgetc(infile);
    }

    eon = 0;

    if (ch == ';') {
      break;
    }

    else if (ch == '(') {  /* Make a node */

      /* Check to see if this is the root node */
      if (maderoot == 0) {
	root_node = (struct node *)malloc(sizeof(struct node));
	maderoot = 1;
	current_node = root_node;
	current_node->br_flag = 0;
	current_node->node_id = nodecount;
	current_node->vi = 0;
	current_node->decay = 1000;
	strcpy(current_node->name, "-");
      }

      /* If not, make a new internal node */
      else {
	nodecount++;
	current_node = makenode(current_node);
	current_node->node_id = nodecount;
	current_node->vi = 0;
	current_node->decay = 1000;
	strcpy(current_node->name, "-");
      }
    }

    /* Add a tip */
    else if (isalnum(ch) != 0) { /* if ch is a letter or digit */
      i = 0;
      nodecount++;
      current_node = makenode(current_node);
      current_node->node_id = nodecount;
      current_node->vi = 0;

      /* Read in chars and copy them into tempname */
      for (i = 0; isalnum(ch) != 0; i++) {
	tempname[i] = ch;
	ch = fgetc(infile);
      }

      tempname[i] = '\0';
      eon = 1;
      strcpy(current_node->name, tempname);
    }

    /* Move down a node */
    else if (ch == ')' | ch == ',') { 
      current_node = current_node->back;
    }

    else {
      printf("Error in tree file!\n");
      printf("%c\t%d\n", ch, ch);
      exit(1);
    }
  }

  return(root_node);
}

void exec_file(void)
{

  /* open the tree file */
  infile = fopen(infilename, "r");
  if (infile == NULL)	{
    (void) printf("Can not open %s\n", infilename);
    exit(1);
  }
  else {
    int i;

    printf("Executing %s\n", infilename);

    numtrees = count_trees();
    printf("Number of trees = %d\n", numtrees);
    printf("Number of characters = %d\n", numchars);

    /* build the tree */
    tree_list = buildtree();
    tree = tree_list;

    /* count the taxa added */
    /*
      num_taxa = count_taxa(tree->root);
      printf("Tree %s contains %d taxa\n", tree->name, num_taxa);
      */

    /* add character state data to the tree */
    printf("Reading data...");
    while (1) {
      tree->root = addstates(tree->root);
      if (tree->id == numtrees)
	break;
      else
	tree = tree->next;
    }
    printf(" done\n");
  }

  return;
}

void list_trees(void)
{
  int i;
  tree = tree_list;

  printf("\nTrees in memory are:\n");
  while (1) {
    printf("%d.  %s\n", tree->id, tree->name);
    if (tree->id < numtrees)
      tree = tree->next;
    else
      break;
  }

  return;
}

struct treestruct * set_tree(int num)
{
  struct treestruct *temp = tree_list;

  while (1) {
    if (temp->id == num)
      break;
    else
      temp = temp->next;
  }

  return(temp);
}

/*********************************************************************
 * reconstruct -- called by graph_char for each character, and loops *
 *    through  the step matrix range, reconstructing states and      *
 *    counting gains and losses                                      *
 * parameters: number of the character interest, tree of interest    *
 * returns: inflection point of the character wrt the x-axis         *
 *    (difference between gains cost and losses cost)                *
 *********************************************************************/
double reconstruct(int c, struct treestruct *t, char gr)
{
  double inf_point;
  double g_prob[2], l_prob[2];
  int flag = 0, ifexact = 0;
  int gindex, lindex, g, l, stop_val;
  int *gl_ptr;
  int i;

  struct charlist *curr = NULL;
  struct charlist *prev = NULL;

  gindex = ((maxgain + 0.0001)/definc);
  lindex = ((maxloss + 0.0001)/definc);
  stop_val = (1/definc);

  gl_ptr = gainloss;

  /* reconstruct ancestral states */
  l_steps = maxloss;
  for (l = lindex; l > stop_val; l--) {
    g_steps = 1.0;
    t->root = downpass(t->root, c);
    t->root = uppass(t->root, 'r', c);

    if (gr == 's') {
      getstates(t->root, c, &curr);

      if (prev != NULL) {
	struct charlist *ctemp = curr, *ptemp = prev;

	while (ctemp != NULL) {
	  if (ctemp->state == ptemp->state)
	    ;
	  else {
	    printf("node %d changes from %c to %c at %2.1f:%2.1f\n", ctemp->id, ptemp->state, ctemp->state, g_steps, l_steps);
	    putdecay(&(t->root), ctemp->id, -(l_steps - definc));
	  }
	  ctemp = ctemp->back;
	  ptemp = ptemp->back;
	}
      }
      prev = curr;
      curr = NULL;
    }

    /* count changes on the tree */
    gl_ptr = countgl(t->root, c);
    fprintf(outfile, "%2.2f     %d     ", (g_steps - l_steps), gl_ptr[0]);
    fprintf(outfile, "%d\n", gl_ptr[2]);
    g_prob[0] = gl_ptr[0];
    g_prob[1] = gl_ptr[1];
    l_prob[0] = gl_ptr[2];
    l_prob[1] = gl_ptr[3];
    /*		fprintf(probfile, "%2.2f %2.2f ", (g_steps - l_steps), (g_prob[0]/g_prob[1]));
		fprintf(probfile, "%2.2f ", (l_prob[0]/l_prob[1]));
		fprintf(probfile, "%d %d %d\n", gl_ptr[1], gl_ptr[3], gl_ptr[4]);
		*/
    t->root = clearpass(t->root, c);

    /* check for inflection point */
    if (gl_ptr[0] == gl_ptr[2]) {
      double temp = (g_steps - l_steps);
      int d = (DIST/2 - 1) + (temp*10 + 0.0001);
      printf("%d: gains = losses at exactly %2.1f\n", c, temp);
      distrib[d]++;
      inf_point = temp;
      ifexact = 1;
    }

    if ((gl_ptr[0] < gl_ptr[2]) && ifexact == 0) {
      double temp = (g_steps - l_steps);
      int d = (DIST/2 - 1) + (temp*10 + 0.0001);
      flag++;
      if (flag == 1) {
	printf("%d: inflection point prior to %2.1f\n", c, temp);
	distrib[d - 1]++;
	inf_point = temp;
      }
    }

    for (i = 0; i < 5; i++) {
      gainloss[i] = 0;
    }

    l_steps = (l_steps - definc);
  }

  g_steps = 1.0;
  for (g = stop_val; g < gindex; g++) {
    l_steps = 1.0;
    t->root = downpass(t->root, c);
    t->root = uppass(t->root, 'r', c);

    if (gr == 's') {
      getstates(t->root, c, &curr);

      if (prev != NULL) {
	struct charlist *ctemp = curr, *ptemp = prev;

	while (ctemp != NULL) {
	  if (ctemp->state == ptemp->state)
	    ;
	  else {
	    printf("node %d changes from %c to %c at %2.1f:%2.1f\n", ctemp->id, ptemp->state, ctemp->state, g_steps, l_steps);
	    putdecay(&(t->root), ctemp->id, (g_steps - definc));
	  }
	  ctemp = ctemp->back;
	  ptemp = ptemp->back;
	}
      }
      prev = curr;
      curr = NULL;
    }

    /* count changes on the tree */
    gl_ptr = countgl(t->root, c);
    fprintf(outfile, "%2.2f     %d     ", (g_steps - l_steps), gl_ptr[0]);
    fprintf(outfile, "%d\n", gl_ptr[2]);
    g_prob[0] = gl_ptr[0];
    g_prob[1] = gl_ptr[1];
    l_prob[0] = gl_ptr[2];
    l_prob[1] = gl_ptr[3];
    /*		fprintf(probfile, "%2.2f %2.2f ", (g_steps - l_steps), (g_prob[0]/g_prob[1]));
		fprintf(probfile, "%2.2f ", (l_prob[0]/l_prob[1]));
		fprintf(probfile, "%d %d %d\n", gl_ptr[1], gl_ptr[3], gl_ptr[4]);
		*/
    t->root = clearpass(t->root, c);

    /* check for inflection point */
    if (gl_ptr[0] == gl_ptr[2]) {
      double temp = (g_steps - l_steps);
      int d = (DIST/2 - 1) + (temp*10 + 0.0001);
      printf("%d: gains = losses at exactly %2.1f\n", c, temp);
      distrib[d]++;
      inf_point = temp;
      ifexact = 1;
    }

    if ((gl_ptr[0] < gl_ptr[2]) && ifexact == 0) {
      double temp = (g_steps - l_steps);
      int d = (DIST/2 - 1) + (temp*10 + 0.0001);
      flag++;
      if (flag == 1) {
	printf("%d: inflection point prior to %2.1f\n", c, temp);
	distrib[d - 1]++;
	inf_point = temp;
      }
    }

    for (i = 0; i < 5; i++) {
      gainloss[i] = 0;
    }

    g_steps = g_steps + definc;
  }
  fprintf(outfile, "\n");
  /*	fprintf(probfile, "\n");*/

  if (gr == 's') {
    treefile = fopen(tfname, "w");
    fprintf(treefile, "#NEXUS\n\nBEGIN TREES;\n\n");
    fprintf(treefile, "     TREE  * %s = ", t->name);
    writetree(&(t->root));
    fprintf(treefile, "ENDBLOCK;\n\n");
    fclose(treefile);
  }

  return(inf_point);
}

/********************************************************
 * combine - find all combinations of ancestral states	*
 ********************************************************/
void combine(int *a, int position, int nstates)
{
  int index;

  for (index = 0; index < nstates; index++) {
    int i = position;
    a[i] = index;
    if (i == 0) {
      /* one unique combination of states found; do whatever */;
    }
    else {
      i--;
      combine(a, i, NSTATES);
    }
  }

  return;
}


void putdecay(struct node **node, int nodeid, double value)
{
  if ((*node)->node_id == nodeid) {
    (*node)->decay = value;
    return;
  }

  else {
    if ((*node)->left != NULL)
      putdecay(&(*node)->left, nodeid, value);
    if ((*node)->right != NULL)
      putdecay(&(*node)->right, nodeid, value);
  }
}


void getdecay(struct node *node)
{
  printf("node %d has decay value %2.1f\n", node->node_id, node->decay);
  if (node->left->br_flag == 2)
    getdecay(node->left);
  if (node->right->br_flag == 2)
    getdecay(node->right);
  return;
}


void writetree(struct node **node)
{
  if ((*node)->br_flag == 2 && (*node)->vi == 0) {
    (*node)->vi++;
    fprintf(treefile, "(");
    writetree(&(*node)->left);
  }
  else if ((*node)->br_flag == 2 && (*node)->vi == 1) {
    (*node)->vi++;
    fprintf(treefile, ",");
    writetree(&(*node)->right);
  }
  else if ((*node)->br_flag == 2 && (*node)->vi == 2) {
    fprintf(treefile, ")");
    if ((*node)->decay >= 0.0 && (*node)->decay < 999) {
      fprintf(treefile, "%2.1f", (*node)->decay);
    }
    else if ((*node)->decay < 0.0) {
      fprintf(treefile, "999%2.1f", -((*node)->decay));
    }
    if ((*node)->node_id == 0) {
      fprintf(treefile, ";\n\n");
      return;
    }
    else {
      writetree(&(*node)->back);
    }
  }
  else if ((*node)->br_flag == 0) {
    /* write name */
    fprintf(treefile, "%s", (*node)->name);
    writetree(&(*node)->back);
  }
  else {
    printf("Error in writetree\n");
    exit(1);
  }
}