
/* version 3.6. (c) Copyright 1993-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 "phylip.h"
#include "seq.h"

#define initialv        0.1     /* starting value of branch length          */

#define over            60   /* maximum width of a tree on screen */


#ifndef OLDC
/* function prototypes */
void   restml_inputnumbers(void);
void   getoptions(void);
void   allocrest(void);
void   setuppie(void);
void   doinit(void);
void   inputoptions(void);
void   restml_inputdata(void);
void   restml_sitesort(void);
void   restml_sitecombine(void);
void   makeweights(void);

void   restml_makevalues(void);
void   getinput(void);
void   copymatrix(transmatrix, transmatrix);
void   maketrans(double, boolean);
void   branchtrans(long, double);
double evaluate(tree *, node *);
void   nuview(node *);
void   makenewv(node *);
void   update(node *);
void   smooth(node *);

void   insert_(node *p, node *);
void   restml_re_move(node **, node **);
void   restml_copynode(node *, node *);
void   restml_copy_(tree *, tree *);
void   buildnewtip(long , tree *);
void   buildsimpletree(tree *);
void   addtraverse(node *, node *, boolean);
void   rearrange(node *, node *);
void   restml_coordinates(node *, double, long *,double *, double *);
void   restml_printree(void);

double sigma(node *, double *);
void   describe(node *);
void   summarize(void);
void   restml_treeout(node *);
void   inittravtree(tree* t,node *);
void   treevaluate(void);
void   maketree(void);
void   globrearrange(void); 
void   adjust_lengths(tree *);
double adjusted_v(double v);
sitelike2 init_sitelike(long sitelength);
void   free_sitelike(sitelike2 sl);
void   copy_sitelike(sitelike2 dest, sitelike2 src,long sitelength);
void   reallocsites(void);
void   alloctrans(tree *t, long nonodes, long sitelength);
long   get_trans(tree* t);
void   free_trans(tree* t, long trans);
void   free_all_trans(tree* t);
void   alloclrsaves(void);
void   freelrsaves(void);
void   resetlrsaves(void);
/* function prototypes */
#endif

Char infilename[FNMLNGTH], outfilename[FNMLNGTH], intreename[FNMLNGTH], 
     outtreename[FNMLNGTH];
long nonodes2, sites, enzymes, weightsum, sitelength, datasets, 
        ith, njumble, jumb=0;
long inseed, inseed0;
boolean  global, jumble, lengths, weights, trout, trunc8, usertree,
           progress, mulsets, firstset, improve, smoothit,inserting = false;
double bestyet;
tree curtree, priortree, bestree, bestree2;
longer seed;
long *enterorder;
steptr aliasweight;
char *progname;
node *qwhere,*addwhere;
/* local rearrangements need to save views.  
   created globally so that reallocation of the same variable is unnecessary */
node **lrsaves;

/* Local variables for maketree, propagated globally for C version: */
long       nextsp, numtrees, maxwhich, col, shimotrees;
double      maxlogl;
boolean     succeeded, smoothed;
#define NTEMPMATS 7
transmatrix *tempmatrix, tempslope, tempcurve;
sitelike2    pie;
double      *l0gl;
double     **l0gf;
Char ch;

/* variables added to keep treeread2() happy */
boolean goteof;
double trweight;
boolean haslengths;


void alloctrans(tree *t, long nonodes, long sitelength)
{
  /* used by restml */
  long i, j;

  t->trans = (transptr)Malloc(nonodes*sizeof(transmatrix));
  for (i = 0; i < nonodes; ++i){
    t->trans[i] = (transmatrix)Malloc((sitelength + 1) * sizeof(double *));
    for (j = 0;j < sitelength + 1; ++j)
      t->trans[i][j] = (double *)Malloc((sitelength + 1) * sizeof(double));
  }
 
  t->freetrans = Malloc(nonodes* sizeof(long));
  for ( i = 0; i < nonodes; i++ )
    t->freetrans[i] = i+1;
  t->transindex = nonodes - 1;
}  /* alloctrans */


long get_trans(tree* t)
{
  long ret;
  ret = t->freetrans[t->transindex];
  t->transindex--;
  return ret;
}


void free_trans(tree* t, long trans)
{
  t->transindex++;
  t->freetrans[t->transindex] = trans;
}


void free_all_trans(tree* t) 
{
  long i;

  for ( i = 0; i < nonodes2; i++ )
    t->freetrans[i] = i;
  t->transindex = nonodes2 - 1;
}


sitelike2 init_sitelike(long sitelength) 
{
  return Malloc((sitelength+1) * sizeof(double));
}


void free_sitelike(sitelike2 sl) 
{
  free(sl);
}


void copy_sitelike(sitelike2 dest, sitelike2 src,long sitelength)
{
  memcpy(dest,src,(sitelength+1)*sizeof(double));
}


void restml_inputnumbers()
{
  /* read and print out numbers of species and sites */
  fscanf(infile, "%ld%ld%ld", &spp, &sites, &enzymes);
  nonodes2 = spp * 2 - 1;
}  /* restml_inputnumbers */


void getoptions()
{
  /* interactively set options */
  long loopcount, loopcount2;
  Char ch;

  fprintf(outfile, "\nRestriction site Maximum Likelihood");
  fprintf(outfile, " method, version %s\n\n",VERSION);
  putchar('\n');
  sitelength = 6;
  trunc8 = true;
  global = false;
  improve = false;
  jumble = false;
  njumble = 1;
  lengths = false;
  outgrno = 1;
  outgropt = false;
  trout = true;
  usertree = false;
  weights = false;
  printdata = false;
  progress = true;
  treeprint = true;
  interleaved = true;
  loopcount = 0;
  for (;;) {
    cleerhome();
    printf("\nRestriction site Maximum Likelihood");
    printf(" method, version %s\n\n",VERSION);
    printf("Settings for this run:\n");
    printf("  U                 Search for best tree?  %s\n",
           (usertree ? "No, use user trees in input file" : "Yes"));
    if (usertree) {
      printf("  N          Use lengths from user trees?  %s\n",
              (lengths ? "Yes" : "No"));
    }
    printf("  A               Are all sites detected?  %s\n",
           (trunc8 ? "No" : "Yes"));
    if (!usertree) {
      printf("  S        Speedier but rougher analysis?  %s\n",
             (improve ? "No, not rough" : "Yes"));
      printf("  G                Global rearrangements?  %s\n",
             (global ? "Yes" : "No"));
      printf("  J   Randomize input order of sequences?  ");
      if (jumble)
        printf("Yes (seed =%8ld,%3ld times)\n", inseed0, njumble);
      else
        printf("No. Use input order\n");
    }
    printf("  L                          Site length?%3ld\n",sitelength);
    printf("  O                        Outgroup root?  %s%3ld\n",
           (outgropt ? "Yes, at sequence number" :
                       "No, use as outgroup species"),outgrno);

    printf("  M           Analyze multiple data sets?");
    if (mulsets)
      printf("  Yes, %2ld sets\n", datasets);
    else
      printf("  No\n");
    printf("  I          Input sequences interleaved?  %s\n",
           (interleaved ? "Yes" : "No, sequential"));
    printf("  0   Terminal type (IBM PC, ANSI, none)?  %s\n",
           ibmpc ? "IBM PC" : ansi  ? "ANSI" : "(none)");
    printf("  1    Print out the data at start of run  %s\n",
           (printdata ? "Yes" : "No"));
    printf("  2  Print indications of progress of run  %s\n",
           (progress ? "Yes" : "No"));
    printf("  3                        Print out tree  %s\n",
           (treeprint ? "Yes" : "No"));
    printf("  4       Write out trees onto tree file?  %s\n",
           (trout ? "Yes" : "No"));
    printf("\n  Y to accept these or type the letter for one to change\n");
    scanf("%c%*[^\n]", &ch);
    getchar();
    if (ch == '\n')
      ch = ' ';
    uppercase(&ch);
    if (ch == 'Y')
      break;
    if (strchr("UNASGJLOTMI01234",ch) != NULL){
      switch (ch) {

      case 'A':
        trunc8 = !trunc8;
        break;
        
      case 'S':
        improve = !improve;
        break;

      case 'G':
        global = !global;
        break;
        
      case 'J':
        jumble = !jumble;
        if (jumble)
          initjumble(&inseed, &inseed0, seed, &njumble);
        else njumble = 1;
        break;
        
      case 'L':
        loopcount2 = 0;
        do {
          printf("New Sitelength?\n");
          scanf("%ld%*[^\n]", &sitelength);
          getchar();
          if (sitelength < 1)
            printf("BAD RESTRICTION SITE LENGTH: %ld\n", sitelength);
          countup(&loopcount2, 10);
        } while (sitelength < 1);
        break;
        
      case 'N':
        lengths = !lengths;
        break;

      case 'O':
        outgropt = !outgropt;
        if (outgropt)
          initoutgroup(&outgrno, spp);
        else outgrno = 1;
        break;
        
      case 'U':
        usertree = !usertree;
        break;
        
      case 'M':
        mulsets = !mulsets;
        if (mulsets)
          initdatasets(&datasets);
        break;
        
      case 'I':
        interleaved = !interleaved;
        break;
        
      case '0':
        initterminal(&ibmpc, &ansi);
        break;
        
      case '1':
        printdata = !printdata;
        break;
        
      case '2':
        progress = !progress;
        break;
        
      case '3':
        treeprint = !treeprint;
        break;
        
      case '4':
        trout = !trout;
        break;
      }
    } else
      printf("Not a possible option!\n");
    countup(&loopcount, 100);
  }
}  /* getoptions */

void reallocsites() 
{
  long i;
  for (i = 0; i < spp; i++) {
    free(y[i]);
    y[i] = (Char *)Malloc(sites*sizeof(Char));
  }
  free(weight);
  free(alias);
  free(aliasweight);
  weight = (steptr)Malloc((sites+1)*sizeof(long));
  alias = (steptr)Malloc((sites+1)*sizeof(long));
  aliasweight = (steptr)Malloc((sites+1)*sizeof(long));
}


void allocrest()
{
  long i;

  y = (Char **)Malloc(spp*sizeof(Char *));
  for (i = 0; i < spp; i++)
    y[i] = (Char *)Malloc(sites*sizeof(Char));
  nayme = (naym *)Malloc(spp*sizeof(naym));
  enterorder = (long *)Malloc(spp*sizeof(long));
  weight = (steptr)Malloc((sites+1)*sizeof(long));
  alias = (steptr)Malloc((sites+1)*sizeof(long));
  aliasweight = (steptr)Malloc((sites+1)*sizeof(long));
}  /* allocrest */


void freelrsaves()
{
  long i,j;
  for ( i = 0 ; i < NLRSAVES ; i++ ) {
    for (j = 0; j < endsite; j++)
      free(lrsaves[i]->x2[j]);
    free(lrsaves[i]->x2);
    free(lrsaves[i]->underflows);
    free(lrsaves[i]);
  }
  free(lrsaves);
}


void resetlrsaves() 
{
  freelrsaves();
  alloclrsaves();
}


void alloclrsaves()
{
  long i,j;
  lrsaves = Malloc(NLRSAVES * sizeof(node*));
  for ( i = 0 ; i < NLRSAVES ; i++ ) {
    lrsaves[i] = Malloc(sizeof(node));
    lrsaves[i]->x2 = Malloc((endsite + 1)*sizeof(sitelike2));
    for ( j = 0 ; j <= endsite ; j++ ) {
      lrsaves[i]->x2[j] = Malloc((sitelength + 1) * sizeof(double));
    }
  }
} /* alloclrsaves */


void setuppie()
{
  /* set up equilibrium probabilities of being a given
     number of bases away from a restriction site */
  long i;
  double sum;

  pie = init_sitelike(sitelength);
  pie[0] = 1.0;
  sum = pie[0];
  for (i = 1; i <= sitelength; i++) {
    pie[i] = 3 * pie[i - 1] * (sitelength - i + 1) / i;
    sum += pie[i];
  }
  for (i = 0; i <= sitelength; i++)
    pie[i] /= sum;
}  /* setuppie */


void doinit()
{
  /* initializes variables */
  long i,j;

  restml_inputnumbers();
  getoptions();
  if (!usertree)
    nonodes2--;
  if (printdata)
    fprintf(outfile, "%4ld Species, %4ld Sites,%4ld Enzymes\n",
            spp, sites, enzymes);
  tempmatrix = Malloc(NTEMPMATS * sizeof(transmatrix));
  for ( i = 0 ; i < NTEMPMATS ; i++ ) {
    tempmatrix[i] = Malloc((sitelength+1) * sizeof(double *));
    for ( j = 0 ; j <= sitelength ; j++)
      tempmatrix[i][j] = (double *)Malloc((sitelength+1) * sizeof(double));
  }
  tempslope = (transmatrix)Malloc((sitelength+1) * sizeof(double *));
  for (i=0; i<=sitelength; i++)
    tempslope[i] = (double *)Malloc((sitelength+1) * sizeof(double));
  tempcurve = (transmatrix)Malloc((sitelength+1) * sizeof(double *));
  for (i=0; i<=sitelength; i++)
    tempcurve[i] = (double *)Malloc((sitelength+1) * sizeof(double));
  setuppie();
  alloctrans(&curtree, nonodes2, sitelength);
  alloctree(&curtree.nodep, nonodes2, false);
  allocrest();
  if (usertree)
    return;
  alloctrans(&bestree, nonodes2, sitelength);
  alloctree(&bestree.nodep, nonodes2, 0);
  alloctrans(&priortree, nonodes2, sitelength);
  alloctree(&priortree.nodep, nonodes2, 0);
  if (njumble == 1) return;
  alloctrans(&bestree2, nonodes2, sitelength);
  alloctree(&bestree2.nodep, nonodes2, 0);
}  /* doinit */


void inputoptions()
{
  /* read the options information */
  Char ch;
  long i, extranum, cursp, curst, curenz;

  if (!firstset) {
    if (eoln(infile))
      scan_eoln(infile);
    fscanf(infile, "%ld%ld%ld", &cursp, &curst, &curenz);
    if (cursp != spp) {
      printf("\nERROR: INCONSISTENT NUMBER OF SPECIES IN DATA SET %4ld\n",
             ith);
      exxit(-1);
    }
    if (curenz != enzymes) {
      printf("\nERROR: INCONSISTENT NUMBER OF ENZYMES IN DATA SET %4ld\n",
             ith);
      exxit(-1);
    }
    sites = curst;
  }
  if ( !firstset )
    reallocsites();

  for (i = 1; i <= sites; i++)
    weight[i] = 1;
  weightsum = sites;
  extranum = 0;
  readoptions(&extranum, "W");
  for (i = 1; i <= extranum; i++) {
    matchoptions(&ch, "W");
    if (ch == 'W')
      inputweights2(1, sites+1, &weightsum, weight, &weights, "RESTML");
  }
  fprintf(outfile, "\n  Recognition sequences all%2ld bases long\n",
          sitelength);
  if (trunc8)
    fprintf(outfile,
      "\nSites absent from all species are assumed to have been omitted\n\n");
  if (weights)
    printweights(outfile, 1, sites, weight, "Sites");
}  /* inputoptions */


void restml_inputdata()
{
  /* read the species and sites data */
  long i, j, k, l, sitesread, sitesnew=0;
  Char ch;
  boolean allread, done;

  if (printdata)
    putc('\n', outfile);
  j = nmlngth + (sites + (sites - 1) / 10) / 2 - 5;
  if (j < nmlngth - 1)
    j = nmlngth - 1;
  if (j > 39)
    j = 39;
  if (printdata) {
    fprintf(outfile, "Name");
    for (i = 1; i <= j; i++)
      putc(' ', outfile);
    fprintf(outfile, "Sites\n");
    fprintf(outfile, "----");
    for (i = 1; i <= j; i++)
      putc(' ', outfile);
    fprintf(outfile, "-----\n\n");
  }
  sitesread = 0;
  allread = false;
  while (!(allread)) {
    /* eat white space -- if the separator line has spaces on it*/
    do {
      ch = gettc(infile);
    } while (ch == ' ' || ch == '\t');
    ungetc(ch, infile);
    if (eoln(infile))
      scan_eoln(infile);
    i = 1;
    while (i <= spp ) {
      if ((interleaved && sitesread == 0) || !interleaved)
        initname(i - 1);
      if (interleaved)
        j = sitesread;
      else
        j = 0;
      done = false;
      while (!done && !eoff(infile)) {
        if (interleaved)
          done = true;
        while (j < sites && !(eoln(infile) || eoff(infile))) {
          ch = gettc(infile);
          if (ch == '\n' || ch == '\t')
            ch = ' ';
          if (ch == ' ')
            continue;
          uppercase(&ch);
          if (ch != '1' && ch != '0' && ch != '+' && ch != '-' && ch != '?') {
            printf(" ERROR: Bad symbol %c", ch);
            printf(" at position %ld of species %ld\n", j+1, i);
            exxit(-1);
          }
          if (ch == '1')
            ch = '+';
          if (ch == '0')
            ch = '-';
          j++;
          y[i - 1][j - 1] = ch;
        }
        if (interleaved)
          continue;
        if (j < sites) 
          scan_eoln(infile);
        else if (j == sites)
          done = true;
      }
      if (interleaved && i == 1)
        sitesnew = j;
      scan_eoln(infile);
      if ((interleaved && j != sitesnew ) || ((!interleaved) && j != sites)){
        printf("ERROR: SEQUENCES OUT OF ALIGNMENT\n");
        exxit(-1);}
      i++;
    }
    if (interleaved) {
      sitesread = sitesnew;
      allread = (sitesread == sites);
    } else
      allread = (i > spp);
  }
  if (printdata) {
    for (i = 1; i <= ((sites - 1) / 60 + 1); i++) {
      for (j = 0; j < spp; j++) {
        for (k = 0; k < nmlngth; k++)
          putc(nayme[j][k], outfile);
        fprintf(outfile, "   ");
        l = i * 60;
        if (l > sites)
          l = sites;
        for (k = (i - 1) * 60 + 1; k <= l; k++) {
          putc(y[j][k - 1], outfile);
          if (k % 10 == 0 && k % 60 != 0)
            putc(' ', outfile);
        }
        putc('\n', outfile);
      }
      putc('\n', outfile);
    }
    putc('\n', outfile);
  }
  putc('\n', outfile);
}  /* restml_inputdata */


void restml_sitesort()
{
  /* Shell sort keeping alias, aliasweight in same order */
  long gap, i, j, jj, jg, k, itemp;
  boolean flip, tied;

  gap = sites / 2;
  while (gap > 0) {
    for (i = gap + 1; i <= sites; i++) {
      j = i - gap;
      flip = true;
      while (j > 0 && flip) {
        jj = alias[j];
        jg = alias[j + gap];
        flip = false;
        tied = true;
        k = 1;
        while (k <= spp && tied) {
          flip = (y[k - 1][jj - 1] > y[k - 1][jg - 1]);
          tied = (tied && y[k - 1][jj - 1] == y[k - 1][jg - 1]);
          k++;
        }
        if (tied) {
          aliasweight[j] += aliasweight[j + gap];
          aliasweight[j + gap] = 0;
        }
        if (!flip)
          break;
        itemp = alias[j];
        alias[j] = alias[j + gap];
        alias[j + gap] = itemp;
        itemp = aliasweight[j];
        aliasweight[j] = aliasweight[j + gap];
        aliasweight[j + gap] = itemp;
        j -= gap;
      }
    }
    gap /= 2;
  }
}  /* restml_sitesort */


void restml_sitecombine()
{
  /* combine sites that have identical patterns */
  long i, j, k;
  boolean tied;

  i = 1;
  while (i < sites) {
    j = i + 1;
    tied = true;
    while (j <= sites && tied) {
      k = 1;
      while (k <= spp && tied) {
        tied = (tied && y[k - 1][alias[i] - 1] == y[k - 1][alias[j] - 1]);
        k++;
      }
      if (tied && aliasweight[j] > 0) {
        aliasweight[i] += aliasweight[j];
        aliasweight[j] = 0;
        alias[j] = alias[i];
      }
      j++;
    }
    i = j - 1;
  }
}  /* restml_sitecombine */


void makeweights()
{
  /* make up weights vector to avoid duplicate computations */
  long i;

  for (i = 1; i <= sites; i++) {
    alias[i] = i;
    aliasweight[i] = weight[i];
  }
  restml_sitesort();
  restml_sitecombine();
  sitescrunch2(sites + 1, 2, 3, aliasweight);
  for (i = 1; i <= sites; i++) {
    weight[i] = aliasweight[i];
    if (weight[i] > 0)
      endsite = i;
  }
  weight[0] = 1;
}  /* makeweights */


void restml_makevalues()
{
  /* set up fractional likelihoods at tips */
  long i, j, k, l, m;
  boolean found;

  for (k = 1; k <= endsite; k++) {
    j = alias[k];
    for (i = 0; i < spp; i++) {
      for (l = 0; l <= sitelength; l++)
        curtree.nodep[i]->x2[k][l] = 1.0;
      switch (y[i][j - 1]) {

      case '+':
        for (m = 1; m <= sitelength; m++)
          curtree.nodep[i]->x2[k][m] = 0.0;
        break;

      case '-':
        curtree.nodep[i]->x2[k][0] = 0.0;
        break;

      case '?':
        /* blank case */
        break;
      }
    }
  }
  for (i = 0; i < spp; i++) {
    for (k = 1; k <= sitelength; k++)
      curtree.nodep[i]->x2[0][k] = 1.0;
    curtree.nodep[i]->x2[0][0] = 0.0;
  }
  if (trunc8)
    return;
  found = false;
  i = 1;
  while (!found && i <= endsite) {
    found = true;
    for (k = 0; k < spp; k++)
      found = (found && y[k][alias[i] - 1] == '-');
    if (!found)
      i++;
  }
  if (found) {
    weightsum += (enzymes - 1) * weight[i];
    weight[i] *= enzymes;
  }
}  /* restml_makevalues */


void getinput()
{
  /* reads the input data */
  inputoptions();
  restml_inputdata();
  if ( !firstset ) freelrsaves();
  makeweights();
  alloclrsaves();
  setuptree2(curtree);
  if (!usertree) {
    setuptree2(priortree);
    setuptree2(bestree);
    if (njumble > 1) 
      setuptree2(bestree2);
  }
  allocx2(nonodes2, endsite+1, sitelength, curtree.nodep, false);
  if (!usertree) {
    allocx2(nonodes2, endsite+1, sitelength, priortree.nodep, 0);
    allocx2(nonodes2, endsite+1, sitelength, bestree.nodep, 0);
    if (njumble > 1)
      allocx2(nonodes2, endsite+1, sitelength, bestree2.nodep, 0);
  }
  restml_makevalues();
}  /* getinput */


void copymatrix(transmatrix tomat, transmatrix frommat)
{
  /* copy a matrix the size of transition matrix */
  int i,j;

  for (i=0;i<=sitelength;++i){
    for (j=0;j<=sitelength;++j)
      tomat[i][j] = frommat[i][j];
  }
} /* copymatrix */


void maketrans(double p, boolean nr)
{
  /* make transition matrix, product matrix with change
     probability p.  Put the results in tempmatrix, tempslope, tempcurve */
  long i, j, k, m1, m2;
  double sump, sums=0, sumc=0, pover3, pijk, term;
  sitelike2 binom1, binom2;

  binom1 = init_sitelike(sitelength);
  binom2 = init_sitelike(sitelength);
  pover3 = p / 3;
  for (i = 0; i <= sitelength; i++) {
    if (p > 1.0 - epsilon)
      p = 1.0 - epsilon;
    if (p < epsilon)
      p = epsilon;
    binom1[0] = exp((sitelength - i) * log(1 - p));
    for (k = 1; k <= sitelength - i; k++)
      binom1[k] = binom1[k - 1] * (p / (1 - p)) * (sitelength - i - k + 1) / k;
    binom2[0] = exp(i * log(1 - pover3));
    for (k = 1; k <= i; k++)
      binom2[k] = binom2[k - 1] * (pover3 / (1 - pover3)) * (i - k + 1) / k;
    for (j = 0; j <= sitelength; ++j) {
      sump = 0.0;
      if (nr) {
        sums = 0.0;
        sumc = 0.0;
      }
      if (i - j > 0)
        m1 = i - j;
      else
        m1 = 0;
      if (sitelength - j < i)
        m2 = sitelength - j;
      else
        m2 = i;
      for (k = m1; k <= m2; k++) {
        pijk = binom1[j - i + k] * binom2[k];
        sump += pijk;
        if (nr) {
          term = (j-i+2*k)/p - (sitelength-j-k)/(1.0-p) - (i-k)/(3.0-p);
          sums += pijk * term;
          sumc += pijk * (term * term
                            - (j-i+2*k)/(p*p)
                            - (sitelength-j-k)/((1.0-p)*(1.0-p))
                            - (i-k)/((3.0-p)*(3.0-p)) );
        }
      }
      tempmatrix[0][i][j] = sump;
      if (nr) {
        tempslope[i][j] = sums;
        tempcurve[i][j] = sumc;
      }
    }
  }
  free_sitelike(binom1);
  free_sitelike(binom2);
}  /* maketrans */


void branchtrans(long i, double p)
{
  /* make branch transition matrix for branch i with probability of change p */
  boolean nr;

  nr = false;
  maketrans(p, nr);
  copymatrix(curtree.trans[i - 1], tempmatrix[0]);
}  /* branchtrans */


double evaluate(tree *tr, node *p)
{
  /* evaluates the likelihood, using info. at one branch */
  double sum, sum2, y, liketerm, like0, lnlike0=0, term;
  long i, j, k,branchnum;
  node *q;
  sitelike2 x1, x2;

  x1 = init_sitelike(sitelength);
  x2 = init_sitelike(sitelength);
  sum = 0.0;
  q = p->back;
  if ( p->initialized  == false && p->tip == false)  nuview(p);
  if ( q->initialized  == false && q->tip == false)  nuview(q);
  y = p->v;
  branchnum = p->branchnum;
  copy_sitelike(x1,p->x2[0],sitelength);
  copy_sitelike(x2,q->x2[0],sitelength);
  if (trunc8) {
    like0 = 0.0;
    for (j = 0; j <= sitelength; j++) {
      liketerm = pie[j] * x1[j];
      for (k = 0; k <= sitelength; k++)
        like0 += liketerm * tr->trans[branchnum-1][j][k] * x2[k];
    }
    lnlike0 = log(enzymes * (1.0 - like0));
  }
  for (i = 1; i <= endsite; i++) {
    copy_sitelike(x1,p->x2[i],sitelength);
    copy_sitelike(x2,q->x2[i],sitelength);
    sum2 = 0.0;
    for (j = 0; j <= sitelength; j++) {
      liketerm = pie[j] * x1[j];
      for (k = 0; k <= sitelength; k++)
        sum2 += liketerm * tr->trans[branchnum-1][j][k] * x2[k];
    }
    term = log(sum2);
    if (trunc8)
      term -= lnlike0;
    if (usertree && (which <= shimotrees))
      l0gf[which - 1][i - 1] = term;
    sum += weight[i] * term;
  }
/* *** debug  put a variable "saveit" in evaluate as third argument as to
   whether to save the KHT suff */
  if (usertree) {
    if(which <= shimotrees)
      l0gl[which - 1] = sum;
    if (which == 1) {
      maxwhich = 1;
      maxlogl = sum;
    } else if (sum > maxlogl) {
      maxwhich = which;
      maxlogl = sum;
    }
  }
  tr->likelihood = sum;
  free_sitelike(x1);
  free_sitelike(x2);
  return sum;
}  /* evaluate */


void nuview(node *p)
{
  /* recompute fractional likelihoods for one part of tree */
  long i, j, k, lowlim;
  double sumq, sumr;
  node *q, *r;
  sitelike2 xq, xr, xp;
  transmatrix tempq, tempr;
  double *tq, *tr;

  xq = init_sitelike(sitelength);
  xr = init_sitelike(sitelength);
  xp = init_sitelike(sitelength);
  if (!p->next->back->tip && !p->next->back->initialized)
    nuview (p->next->back);
  if (!p->next->next->back->tip && !p->next->next->back->initialized)
    nuview (p->next->next->back);
  tempq = (transmatrix)Malloc((sitelength+1) * sizeof(double *));
  tempr = (transmatrix)Malloc((sitelength+1) * sizeof(double *));
  for (i=0;i<=sitelength;++i){
    tempq[i] = (double *)Malloc((sitelength+1) * sizeof (double));
    tempr[i] = (double *)Malloc((sitelength+1) * sizeof (double));
  }
  if (trunc8)
    lowlim = 0;
  else
    lowlim = 1;
  q = p->next->back;
  r = p->next->next->back;
  copymatrix(tempq,curtree.trans[q->branchnum - 1]);
  copymatrix(tempr,curtree.trans[r->branchnum - 1]);
  for (i = lowlim; i <= endsite; i++) {
    copy_sitelike(xq,q->x2[i],sitelength);
    copy_sitelike(xr,r->x2[i],sitelength);
    for (j = 0; j <= sitelength; j++) {
      sumq = 0.0;
      sumr = 0.0;
      tq = tempq[j];
      tr = tempr[j];
      for (k = 0; k <= sitelength; k++) {
        sumq += tq[k] * xq[k];
        sumr += tr[k] * xr[k];
      }
      xp[j] = sumq * sumr;
    }
    copy_sitelike(p->x2[i],xp,sitelength);
  }
  for (i=0;i<=sitelength;++i){
    free(tempq[i]);
    free(tempr[i]);
  }
  free(tempq);
  free(tempr);
  p->initialized = true;
  free_sitelike(xq);
  free_sitelike(xr);
  free_sitelike(xp);
}  /* nuview */


void makenewv(node *p)
{
  /* Newton-Raphson algorithm improvement of a branch length */
  long i, j, k, lowlim, it, ite;
  double sum, sums, sumc, like, slope, curve, liketerm, liket,
         y, yold=0, yorig, like0=0, slope0=0, curve0=0, oldlike=0, temp;
  boolean done, nr, firsttime, better;
  node *q;
  sitelike2 xx1, xx2;
  double *tm, *ts, *tc;

  xx1 = init_sitelike(sitelength);
  xx2 = init_sitelike(sitelength);
  q = p->back;
  y = p->v;
  yorig = y;
  if (trunc8)
    lowlim = 0;
  else
    lowlim = 1;
  done = false;
  nr = true;
  firsttime = true;
  it = 1;
  ite = 0;
  while ((it < iterations) && (ite < 20) && (!done)) {
    like = 0.0;
    slope = 0.0;
    curve = 0.0;
    maketrans(y, nr);
    for (i = lowlim; i <= endsite; i++) {
      copy_sitelike(xx1,p->x2[i],sitelength);
      copy_sitelike(xx2,q->x2[i],sitelength);
      sum = 0.0;
      sums = 0.0;
      sumc = 0.0;
      for (j = 0; j <= sitelength; j++) {
        liket = xx1[j] * pie[j];
        tm = tempmatrix[0][j];
        ts = tempslope[j];
        tc = tempcurve[j];
        for (k = 0; k <= sitelength; k++) {
          liketerm = liket * xx2[k];
          sum += tm[k] * liketerm;
          sums += ts[k] * liketerm;
          sumc += tc[k] * liketerm;
        }
      }
      if (i == 0) {
        like0 = sum;
        slope0 = sums;
        curve0 = sumc;
      } else {
        like += weight[i] * log(sum);
        slope += weight[i] * sums/sum;
        temp = sums/sum;
        curve += weight[i] * (sumc/sum-temp*temp);
      }
    }
    if (trunc8 && fabs(like0 - 1.0) > 1.0e-10) {
      like -= weightsum * log(enzymes * (1.0 - like0));
      slope += weightsum * slope0 /(1.0 - like0);
      curve += weightsum * (curve0 /(1.0 - like0)
                            + slope0*slope0/((1.0 - like0)*(1.0 - like0)));
    }
    better = false;
    if (firsttime) {
      yold = y;
      oldlike = like;
      firsttime = false;
      better = true;
    } else {
      if (like > oldlike) {
        yold = y;
        oldlike = like;
        better = true;
        it++;
      }
    }
    if (better) {
      y = y + slope/fabs(curve);
      if (y < epsilon)
        y = 10.0 * epsilon;
      if (y > 0.75)
        y = 0.75;
    } else {
        if (fabs(y - yold) < epsilon)
          ite = 20;
        y = (y + yold) / 2.0;
    }
    ite++;
    done = fabs(y-yold) < epsilon;
  }
  smoothed = (fabs(yold-yorig) < epsilon) && (yorig > 1000.0*epsilon);
  p->v = yold;
  q->v = yold;
  branchtrans(p->branchnum, yold);
  curtree.likelihood = oldlike;
  free_sitelike(xx1);
  free_sitelike(xx2);

}  /* makenewv */


void update(node *p)
{
  /* improve branch length and views for one branch */
  if (!p->tip && !p->initialized)
    nuview(p);
  if (!p->back->tip && !p->back->initialized)
    nuview(p->back);
  if ( !(usertree && lengths) ) {
    makenewv(p);
    if (smoothit ) {
      inittrav(p);
      inittrav(p->back);
    }
    else {
     if (inserting && !p->tip) {
        p->next->initialized = false;
        p->next->next->initialized = false;
      }
    }
  }
}  /* update */


void smooth(node *p)
{
  /* update nodes throughout the tree, recursively */

  smoothed = false;
  update(p);
  if (!p->tip) {
    if (smoothit && !smoothed) {
      smooth(p->next->back);
    }
    if (smoothit && !smoothed) {
      smooth(p->next->next->back);
    }
  }
}  /* smooth */


void insert_(node *p, node *q)
{
  /* insert a subtree into a branch, improve lengths in tree */
  long i;
  node *r;

  r = p->next->next;
  hookup(r, q->back);
  hookup(p->next, q);
  if (q->v >= 0.75)
    q->v = 0.75;
  else
    q->v = 0.75 * (1 - sqrt(1 - 1.333333 * q->v));
  if ( q->v < epsilon)
    q->v = epsilon;

  q->back->v = q->v;
  r->v = q->v;
  r->back->v = r->v;
 
  q->back->branchnum = q->branchnum;
  r->branchnum = get_trans(&curtree);
  r->back->branchnum = r->branchnum;
  branchtrans(q->branchnum, q->v);
  branchtrans(r->branchnum, r->v);

  if ( smoothit ) {
    inittrav(p);
    inittrav(p->back);
  } 
  p->initialized = false;

  i = 1;
  inserting = true;
  while (i <= smoothings) {
    smooth(p);
    if (!p->tip) {
      smooth (p->next->back);
      smooth (p->next->next->back);
    }
    i++;
  }
  inserting = false;
}  /* insert_ */


void restml_re_move(node **p, node **q)
{
  /* remove p and record in q where it was */
  long i;

  *q = (*p)->next->back;
  hookup(*q, (*p)->next->next->back);
  free_trans(&curtree,(*q)->back->branchnum);
  (*q)->back->branchnum = (*q)->branchnum;
  (*q)->v = 0.75*(1 - (1 - 1.333333*(*q)->v) * (1 - 1.333333*(*p)->next->v));

  if ( (*q)->v > 1 - epsilon)
    (*q)->v = 1 - epsilon;
  else if ( (*q)->v < epsilon)
    (*q)->v = epsilon;
                  
  (*q)->back->v = (*q)->v;
  branchtrans((*q)->branchnum, (*q)->v);

  (*p)->next->back = NULL;
  (*p)->next->next->back = NULL;
 
  if ( smoothit ) {
    inittrav((*q)->back);
    inittrav(*q);
  }
  
  if ( smoothit ) {
    for ( i = 0 ; i < smoothings ; i++ ) {
      smooth(*q);
      smooth((*q)->back);
    }
  }
  else ( smooth(*q));
}  /* restml_re_move */


void restml_copynode(node *c, node *d)
{
  /* copy a node */
  long i;

  d->branchnum = c->branchnum;
  
  for ( i = 0 ; i <= endsite ; i++)
    copy_sitelike(d->x2[i],c->x2[i],sitelength);
  d->v = c->v;
  d->iter = c->iter;
  d->xcoord = c->xcoord;
  d->ycoord = c->ycoord;
  d->ymin = c->ymin;
  d->ymax = c->ymax;
  d->initialized = c->initialized;
}  /* restml_copynode */


void restml_copy_(tree *a, tree *b)
{
  /* copy tree a to tree b */
  long i,j;
  node *p, *q;

  for (i = 0; i < spp; i++) {
    restml_copynode(a->nodep[i], b->nodep[i]);
    if (a->nodep[i]->back) {
      if (a->nodep[i]->back == a->nodep[a->nodep[i]->back->index - 1])
        b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1];
      else if (a->nodep[i]->back
                == a->nodep[a->nodep[i]->back->index - 1]->next)
          b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1]->next;
        else
          b->nodep[i]->back
                         = b->nodep[a->nodep[i]->back->index - 1]->next->next;
    }
    else b->nodep[i]->back = NULL;
  }
  for (i = spp; i < nonodes2; i++) {
    p = a->nodep[i];
    q = b->nodep[i];
    for (j = 1; j <= 3; j++) {
      restml_copynode(p, q);
      if (p->back) {
        if (p->back == a->nodep[p->back->index - 1])
          q->back = b->nodep[p->back->index - 1];
        else if (p->back == a->nodep[p->back->index - 1]->next)
          q->back = b->nodep[p->back->index - 1]->next;
        else
          q->back = b->nodep[p->back->index - 1]->next->next;
      }
      else
        q->back = NULL;
      p = p->next;
      q = q->next;
    }
  }
  b->likelihood = a->likelihood;
  for (i=0;i<nonodes2;++i)
    copymatrix(b->trans[i],a->trans[i]);
  b->transindex = a->transindex;
  memcpy(b->freetrans,a->freetrans,nonodes*sizeof(long));

  b->start = a->start;
}  /* restml_copy */


void buildnewtip(long m,tree *tr)
{
  /* set up a new tip and interior node it is connected to */
  node *p;
  long i, j;

  p = tr->nodep[nextsp + spp - 3];
  for (i = 0; i <= endsite; i++) {  /* FIXME, ask joe about this */
    for (j = 0; j < sitelength; j++) { /* trunc8 */
      p->x2[i][j] = 1.0; 
      p->next->x2[i][j] = 1.0; 
      p->next->next->x2[i][j] = 1.0; 
    }
  }

  hookup(tr->nodep[m - 1], p);
  p->v = initialv;
  p->back->v = initialv;
  p->branchnum = get_trans(tr);
  p->back->branchnum = p->branchnum;
  branchtrans(p->branchnum, initialv);
}  /* buildnewtip */


void buildsimpletree(tree *tr)
{
  /* set up and adjust branch lengths of a three-species tree */
  long branch;

  hookup(tr->nodep[enterorder[0] - 1], tr->nodep[enterorder[1] - 1]);
  tr->nodep[enterorder[0] - 1]->v = initialv;
  tr->nodep[enterorder[1] - 1]->v = initialv;
  branchtrans(enterorder[1], initialv);
  branch = get_trans(tr);
  tr->nodep[enterorder[0] - 1]->branchnum = branch;
  tr->nodep[enterorder[1] - 1]->branchnum = branch;
  buildnewtip(enterorder[2], tr);
  insert_(tr->nodep[enterorder[2] - 1]->back, tr->nodep[enterorder[1] - 1]);
  tr->start = tr->nodep[enterorder[2]-1]->back;
}  /* buildsimpletree */


void addtraverse(node *p, node *q, boolean contin)
{
  /* try adding p at q, proceed recursively through tree */
  double like, vsave = 0;
  node *qback =NULL;

  if (!smoothit) {
    copymatrix (tempmatrix[1], curtree.trans[q->branchnum - 1]);
    vsave = q->v;
    qback = q->back;
  }
  insert_(p, q);
  like = evaluate(&curtree, p);
  if (like > bestyet) {
    bestyet = like;
    if (smoothit) {
      restml_copy_(&curtree, &bestree);
      addwhere = q;
    }
    else
      qwhere = q;
    succeeded = true;
  }
  if (smoothit)
    restml_copy_(&priortree, &curtree);
  else {
    hookup (q, qback);
    q->v = vsave;
    q->back->v = vsave;
    free_trans(&curtree,q->back->branchnum);
    q->back->branchnum = q->branchnum;
    copymatrix (curtree.trans[q->branchnum - 1], tempmatrix[1]);
    curtree.likelihood = bestyet;
  }
  if (!q->tip && contin) {
    addtraverse(p, q->next->back, contin);
    addtraverse(p, q->next->next->back, contin);
  }
  if ( contin && q == curtree.root ) {
    addtraverse(p,q->back,contin);
  }
}  /* addtraverse */


void globrearrange() 
{
  /* does global rearrangements */
  tree globtree;
  tree oldtree;
  int i,j,k,l,num_sibs,num_sibs2;
  node *where,*sib_ptr,*sib_ptr2;
  double oldbestyet = curtree.likelihood;
  int success = false;
  printf("\n   ");
  alloctree(&globtree.nodep,nonodes2,0);
  alloctree(&oldtree.nodep,nonodes2,0);
  alloctrans(&globtree, nonodes2, sitelength);
  alloctrans(&oldtree, nonodes2, sitelength);
  setuptree2(globtree);
  setuptree2(oldtree);
  allocx2(nonodes2, endsite + 1, sitelength,globtree.nodep, 0);
  allocx2(nonodes2, endsite + 1, sitelength,oldtree.nodep, 0);
  restml_copy_(&curtree,&globtree);
  restml_copy_(&curtree,&oldtree);
  bestyet = curtree.likelihood;
  for ( i = spp ; i < nonodes2 ; i++ ) {
    num_sibs = count_sibs(curtree.nodep[i]);
    sib_ptr  = curtree.nodep[i];
    if ( (i - spp) % (( nonodes2 / 72 ) + 1 ) == 0 )
      putchar('.');
    fflush(stdout);
    for ( j = 0 ; j <= num_sibs ; j++ ) {
      restml_re_move(&sib_ptr,&where);
      restml_copy_(&curtree,&priortree);
      qwhere = where;
      
      if (where->tip) {
        restml_copy_(&oldtree,&curtree);
        restml_copy_(&oldtree,&bestree);
        sib_ptr=sib_ptr->next;
        continue;
      }
      else num_sibs2 = count_sibs(where);
      sib_ptr2 = where;
      for ( k = 0 ; k < num_sibs2 ; k++ ) {
        addwhere = NULL;
        addtraverse(sib_ptr,sib_ptr2->back,true);
        if ( !smoothit ) {
          if (succeeded && qwhere != where && qwhere != where->back) {
            insert_(sib_ptr,qwhere);
            smoothit = true;
            for (l = 1; l<=smoothings; l++) {
              smooth (where);
              smooth (where->back);
            }
            smoothit = false;
            success = true;
            restml_copy_(&curtree,&globtree);
          }
          restml_copy_(&priortree,&curtree);
        }
        else if ( addwhere && where != addwhere && where->back != addwhere
              && bestyet > globtree.likelihood) {
            restml_copy_(&bestree,&globtree);
            success = true;
        }
        sib_ptr2 = sib_ptr2->next;
      } 
      restml_copy_(&oldtree,&curtree);
      restml_copy_(&oldtree,&bestree);
      sib_ptr = sib_ptr->next;
    }
  }
  restml_copy_(&globtree,&curtree);
  restml_copy_(&globtree,&bestree);
  if (success && globtree.likelihood > oldbestyet)  {
    succeeded = true;
  }
  else  {
    succeeded = false;
  }
  bestyet = globtree.likelihood;
  freetrans(&globtree.trans, nonodes2, sitelength);
  freetrans(&oldtree.trans, nonodes2, sitelength);
  freetree2(globtree.nodep,nonodes2);
  freetree2(oldtree.nodep,nonodes2);
}

void printnode(node* p);
void printnode(node* p)
{
  if (p->back)
  printf("p->index = %3ld, p->back->index = %3ld, p->branchnum = %3ld,evaluates"
      " to %lf\n",p->index,p->back->index,p->branchnum,evaluate(&curtree,p));
  
  else 
  printf("p->index = %3ld, p->back->index =none, p->branchnum = %3ld,evaluates"
      " to nothing\n",p->index,p->branchnum);
}

void printvals(void);
void printvals(void) 
{
  int i;
  node* p;

  for ( i = 0 ; i < nextsp ; i++ )  {
    p = curtree.nodep[i];
    printnode(p);
  }
  for ( i = spp ; i <= spp + nextsp - 3 ; i++ ) {
    p = curtree.nodep[i];
    printnode(p);
    printnode(p->next);
    printnode(p->next->next);
  }
}

void rearrange(node *p, node *pp)
{
  /* rearranges the tree locally */
  long i;
  node *q, *r;
  node *rnb, *rnnb;

  if (!p->tip && !p->back->tip) {
    bestyet = curtree.likelihood;
    if (p->back->next != pp)
      r = p->back->next;
    else
      r = p->back->next->next;
    if (!smoothit) {
      rnb = r->next->back;
      rnnb = r->next->next->back;
      restml_copynode(r,lrsaves[0]);
      restml_copynode(r->next,lrsaves[1]);
      restml_copynode(r->next->next,lrsaves[2]);
      restml_copynode(p->next,lrsaves[3]);
      restml_copynode(p->next->next,lrsaves[4]);
      copymatrix (tempmatrix[2], curtree.trans[r->branchnum - 1]);
      copymatrix (tempmatrix[3], curtree.trans[r->next->branchnum - 1]);
      copymatrix (tempmatrix[4], curtree.trans[r->next->next->branchnum-1]);
      copymatrix (tempmatrix[5], curtree.trans[p->next->branchnum-1]);
      copymatrix (tempmatrix[6], curtree.trans[p->next->next->branchnum-1]);
    }
    else
      restml_copy_(&curtree, &bestree);
    restml_re_move(&r, &q);
    nuview(p->next);
    nuview(p->next->next);
    if (smoothit)
      restml_copy_(&curtree, &priortree);
    else
      qwhere = q;
    addtraverse(r, p->next->back, false);
    addtraverse(r, p->next->next->back, false);
    if (smoothit)
      restml_copy_(&bestree, &curtree);
    else {
      if (qwhere == q) {
        hookup(rnb,r->next);
        hookup(rnnb,r->next->next);
        restml_copynode(lrsaves[0],r);
        restml_copynode(lrsaves[1],r->next);
        restml_copynode(lrsaves[2],r->next->next);
        restml_copynode(lrsaves[3],p->next);
        restml_copynode(lrsaves[4],p->next->next);
        r->back->v = r->v;
        r->next->back->v = r->next->v;
        r->next->next->back->v = r->next->next->v;
        p->next->back->v = p->next->v;
        p->next->next->back->v = p->next->next->v;
        r->back->branchnum = r->branchnum;
        r->next->back->branchnum = r->next->branchnum;
        r->next->next->back->branchnum = r->next->next->branchnum;
        p->next->back->branchnum = p->next->branchnum;
        p->next->next->back->branchnum = p->next->next->branchnum;
        copymatrix (curtree.trans[r->branchnum-1], tempmatrix[2]);
        copymatrix (curtree.trans[r->next->branchnum-1], tempmatrix[3]);
        copymatrix (curtree.trans[r->next->next->branchnum-1], tempmatrix[4]);
        copymatrix (curtree.trans[p->next->branchnum-1], tempmatrix[5]);
        copymatrix (curtree.trans[p->next->next->branchnum-1], tempmatrix[6]);
        curtree.likelihood = bestyet;
      }
      else {
        smoothit = true;
        insert_(r, qwhere);
        
        for (i = 1; i<=smoothings; i++) {
          smooth (r);
          smooth (r->back);
        }
        smoothit = false;
      }
    }
  }
  if (!p->tip) {
    rearrange(p->next->back, p);
    rearrange(p->next->next->back, p);
  }
}  /* rearrange */


void restml_coordinates(node *p, double lengthsum, long *tipy,
                double *tipmax, double *x)
{
  /* establishes coordinates of nodes */
  node *q, *first, *last;

  if (p->tip) {
    p->xcoord = (long)(over * lengthsum + 0.5);
    p->ycoord = (*tipy);
    p->ymin = (*tipy);
    p->ymax = (*tipy);
    (*tipy) += down;
    if (lengthsum > (*tipmax))
      (*tipmax) = lengthsum;
    return;
  }
  q = p->next;
  do {
    (*x) = -0.75 * log(1.0 - 1.333333 * q->v);
    restml_coordinates(q->back, lengthsum + (*x),tipy,tipmax,x);
    q = q->next;
  } while ((p == curtree.start || p != q) &&
           (p != curtree.start || p->next != q));
  first = p->next->back;
  q = p;
  while (q->next != p)
    q = q->next;
  last = q->back;
  p->xcoord = (long)(over * lengthsum + 0.5);
  if (p == curtree.start)
    p->ycoord = p->next->next->back->ycoord;
  else
    p->ycoord = (first->ycoord + last->ycoord) / 2;
  p->ymin = first->ymin;
  p->ymax = last->ymax;
}  /* restml_coordinates */


void restml_printree()
{
  /* prints out diagram of the tree */
  long tipy,i;
  double scale, tipmax, x;

  putc('\n', outfile);
  if (!treeprint)
    return;
  putc('\n', outfile);
  tipy = 1;
  tipmax = 0.0;
  restml_coordinates(curtree.start, 0.0, &tipy,&tipmax,&x);
  scale = 1.0 / (tipmax + 1.000);
  for (i = 1; i <= tipy - down; i++)
    drawline2(i, scale, curtree);
  putc('\n', outfile);
}  /* restml_printree */


double sigma(node *q, double *sumlr)
{
  /* get 1.95996 * approximate standard error of branch length */
  double sump, sumr, sums, sumc, p, pover3, pijk, Qjk, liketerm, f;
  double  slopef,curvef;
  long i, j, k, m1, m2;
  sitelike2 binom1, binom2;
  transmatrix Prob, slopeP, curveP;
  node *r;
  sitelike2 x1, x2;
  double term, TEMP;

  x1 = init_sitelike(sitelength);
  x2 = init_sitelike(sitelength);
  binom1 = init_sitelike(sitelength);
  binom2 = init_sitelike(sitelength);
  Prob   = (transmatrix)Malloc((sitelength+1) * sizeof(double *));
  slopeP = (transmatrix)Malloc((sitelength+1) * sizeof(double *));
  curveP = (transmatrix)Malloc((sitelength+1) * sizeof(double *));
  for (i=0; i<=sitelength; ++i) {
    Prob[i]   = (double *)Malloc((sitelength+1) * sizeof(double));
    slopeP[i] = (double *)Malloc((sitelength+1) * sizeof(double));
    curveP[i] = (double *)Malloc((sitelength+1) * sizeof(double));
  }
  p = q->v;
  pover3 = p / 3;
  for (i = 0; i <= sitelength; i++) {
    binom1[0] = exp((sitelength - i) * log(1 - p));
    for (k = 1; k <= (sitelength - i); k++)
      binom1[k] = binom1[k - 1] * (p / (1 - p)) * (sitelength - i - k + 1) / k;
    binom2[0] = exp(i * log(1 - pover3));
    for (k = 1; k <= i; k++)
      binom2[k] = binom2[k - 1] * (pover3 / (1 - pover3)) * (i - k + 1) / k;
    for (j = 0; j <= sitelength; j++) {
      sump = 0.0;
      sums = 0.0;
      sumc = 0.0;
      if (i - j > 0)
        m1 = i - j;
      else
        m1 = 0;
      if (sitelength - j < i)
        m2 = sitelength - j;
      else
        m2 = i;
      for (k = m1; k <= m2; k++) {
        pijk = binom1[j - i + k] * binom2[k];
        sump += pijk;
        term = (j-i+2*k)/p - (sitelength-j-k)/(1.0-p) - (i-k)/(3.0-p);
        sums += pijk * term;
        sumc += pijk * (term * term
                          - (j-i+2*k)/(p*p)
                          - (sitelength-j-k)/((1.0-p)*(1.0-p))
                          - (i-k)/((3.0-p)*(3.0-p)) );
      }
      Prob[i][j] = sump;
      slopeP[i][j] = sums;
      curveP[i][j] = sumc;
    }
  }
  (*sumlr) = 0.0;
  sumc = 0.0;
  sums = 0.0;
  r = q->back;
  for (i = 1; i <= endsite; i++) {
    f = 0.0;
    slopef = 0.0;
    curvef = 0.0;
    sumr = 0.0;
    copy_sitelike(x1,q->x2[i],sitelength);
    copy_sitelike(x2,r->x2[i],sitelength);
    for (j = 0; j <= sitelength; j++) {
      liketerm = pie[j] * x1[j];
      sumr += liketerm * x2[j];
      for (k = 0; k <= sitelength; k++) {
        Qjk = liketerm * x2[k];
        f += Qjk * Prob[j][k];
        slopef += Qjk * slopeP[j][k];
        curvef += Qjk * curveP[j][k];
      }
    }
    (*sumlr) += weight[i] * log(f / sumr);
    sums += weight[i] * slopef / f;
    TEMP = slopef / f;
    sumc += weight[i] * (curvef / f - TEMP * TEMP);
  }
  if (trunc8) {
    f = 0.0;
    slopef = 0.0;
    curvef = 0.0;
    sumr = 0.0;
    copy_sitelike(x1,q->x2[0],sitelength);
    copy_sitelike(x2,r->x2[0],sitelength);
    for (j = 0; j <= sitelength; j++) {
      liketerm = pie[j] * x1[j];
      sumr += liketerm * x2[j];
      for (k = 0; k <= sitelength; k++) {
        Qjk = liketerm * x2[k];
        f += Qjk * Prob[j][k];
        slopef += Qjk * slopeP[j][k];
        curvef += Qjk * curveP[j][k];
      }
    }
    (*sumlr) += weightsum * log((1.0 - sumr) / (1.0 - f));
    sums += weightsum * slopef / (1.0 - f);
    TEMP = slopef / (1.0 - f);
    sumc += weightsum * (curvef / (1.0 - f) + TEMP * TEMP);
  }
  for (i=0;i<=sitelength;++i){
    free(Prob[i]);
    free(slopeP[i]);
    free(curveP[i]);
  }
  free(Prob);
  free(slopeP);
  free(curveP);
  free_sitelike(x1);
  free_sitelike(x2);
  free_sitelike(binom1);
  free_sitelike(binom2);
  if (sumc < -1.0e-6)
    return ((-sums - sqrt(sums * sums - 3.841 * sumc)) / sumc);
  else
    return -1.0;
}  /* sigma */


void describe(node *p)
{
  /* print out information on one branch */
  double sumlr;
  long i;
  node *q;
  double s;

  q = p->back;
  fprintf(outfile, "%4ld      ", q->index - spp);
  fprintf(outfile, "    ");
  if (p->tip) {
    for (i = 0; i < nmlngth; i++)
      putc(nayme[p->index - 1][i], outfile);
  } else
    fprintf(outfile, "%4ld      ", p->index - spp);
  if (q->v >= 0.75)
    fprintf(outfile, "     infinity");
  else
    fprintf(outfile, "%13.5f", -0.75 * log(1 - 1.333333 * q->v));
  if (p->iter) {
    s = sigma(q, &sumlr);
    if (s < 0.0)
      fprintf(outfile, "     (     zero,    infinity)");
    else {
      fprintf(outfile, "     (");
      if (q->v - s <= 0.0)
        fprintf(outfile, "     zero");
      else
        fprintf(outfile, "%9.5f", -0.75 * log(1 - 1.333333 * (q->v - s)));
      putc(',', outfile);
      if (q->v + s >= 0.75)
        fprintf(outfile, "    infinity");
      else
        fprintf(outfile, "%12.5f", -0.75 * log(1 - 1.333333 * (q->v + s)));
      putc(')', outfile);
      }
    if (sumlr > 1.9205)
      fprintf(outfile, " *");
    if (sumlr > 2.995)
      putc('*', outfile);
    }
  else
    fprintf(outfile, "            (not varied)");
  putc('\n', outfile);
  if (!p->tip) {
    describe(p->next->back);
    describe(p->next->next->back);
  }
}  /* describe */


void summarize()
{
  /* print out information on branches of tree */

  fprintf(outfile, "\nremember: ");
  if (outgropt)
    fprintf(outfile, "(although rooted by outgroup) ");
  fprintf(outfile, "this is an unrooted tree!\n\n");
  fprintf(outfile, "Ln Likelihood = %11.5f\n\n", curtree.likelihood);
  fprintf(outfile, " \n");
  fprintf(outfile, " Between        And            Length");
  fprintf(outfile, "      Approx. Confidence Limits\n");
  fprintf(outfile, " -------        ---            ------");
  fprintf(outfile, "      ------- ---------- ------\n");
  describe(curtree.start->next->back);
  describe(curtree.start->next->next->back);
  describe(curtree.start->back);
  fprintf(outfile, "\n     *  = significantly positive, P < 0.05\n");
  fprintf(outfile, "     ** = significantly positive, P < 0.01\n\n\n");
}  /* summarize */


void restml_treeout(node *p)
{
  /* write out file with representation of final tree */
  long i, n, w;
  Char c;
  double x;

  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, outtree);
    }
    col += n;
  } else {
    putc('(', outtree);
    col++;
    restml_treeout(p->next->back);
    putc(',', outtree);
    col++;
    if (col > 45) {
      putc('\n', outtree);
      col = 0;
    }
    restml_treeout(p->next->next->back);
    if (p == curtree.start) {
      putc(',', outtree);
      col++;
      if (col > 45) {
        putc('\n', outtree);
        col = 0;
      }
      restml_treeout(p->back);
    }
    putc(')', outtree);
    col++;
  }
  if (p->v >= 0.75)
    x = -1.0;
  else
    x = -0.75 * log(1 - 1.333333 * p->v);
  if (x > 0.0)
    w = (long)(0.43429448222 * log(x));
  else if (x == 0.0)
    w = 0;
  else
    w = (long)(0.43429448222 * log(-x)) + 1;
  if (w < 0)
    w = 0;
  if (p == curtree.start)
    fprintf(outtree, ";\n");
  else {
    fprintf(outtree, ":%*.5f", (int)(w + 7), x);
    col += w + 8;
  }
}  /* restml_treeout */


void inittravtree(tree* t,node *p)
{ /* traverse tree to set initialized and v to initial values */
  node* q; 

  if ( p->branchnum == 0) {
    p->branchnum = get_trans(t);
    p->back->branchnum = p->branchnum;
  }
  p->initialized = false;
  p->back->initialized = false;
  if ((!lengths) || p->iter) {
    branchtrans(p->branchnum, initialv);
    p->v = initialv;
    p->back->v = initialv;
  }
  else branchtrans(p->branchnum, p->v);

  if ( !p->tip ) {
    q = p->next;
    while ( q != p ) {
      inittravtree(t,q->back);
      q = q->next;
    }
  }
} /* inittravtree */


double adjusted_v(double v) {
  return 3.0/4.0 * (1.0-exp(-4.0/3.0 * v));
}

void adjust_lengths(tree* t) {
  long i;
  for ( i = 0 ; i < spp ; i++) {
    t->nodep[i]->v = adjusted_v(t->nodep[i]->v);
  }
  for ( i = spp ; i < nonodes2 ; i++) {
    t->nodep[i]->v = adjusted_v(t->nodep[i]->v);
    t->nodep[i]->next->v = adjusted_v(t->nodep[i]->next->v);
    t->nodep[i]->next->next->v = adjusted_v(t->nodep[i]->next->next->v);
  }
}

void treevaluate()
{
  /* find maximum likelihood branch lengths of user tree */
  long i;
  unroot(&curtree,nonodes2);
  if ( lengths)  
    adjust_lengths(&curtree);
  nonodes2--;

  inittravtree(&curtree,curtree.start);
  inittravtree(&curtree,curtree.start->back);
  smoothit = true;
  for (i = 1; i <= smoothings * 4; i++)  {
    smooth (curtree.start);
    smooth (curtree.start->back);
  }
  evaluate(&curtree, curtree.start);
  
  nonodes2++;
}  /* treevaluate */


void maketree()
{
  /* construct and rearrange tree */
  long i,j;

  if (usertree) {
    openfile(&intree,INTREE,"input tree file","r",progname,intreename);
    numtrees = countsemic(&intree);
    if(numtrees > MAXSHIMOTREES)
      shimotrees = MAXSHIMOTREES;
    else
      shimotrees = numtrees;
    if (numtrees > 2)
      initseed(&inseed, &inseed0, seed);
    l0gl = (double *) Malloc(shimotrees * sizeof(double));
    l0gf = (double **) Malloc(shimotrees * sizeof(double *));
    for (i=0; i < shimotrees; ++i)
      l0gf[i] = (double *)Malloc(endsite * sizeof(double));
    if (treeprint) {
      fprintf(outfile, "User-defined tree");
      if (numtrees > 1)
        putc('s', outfile);
      fprintf(outfile, ":\n\n");
    }
    which = 1;
    while (which <= numtrees) {
      treeread2 (intree, &curtree.start, curtree.nodep, lengths, &trweight, 
                      &goteof, &haslengths, &spp,false,nonodes2);
      if ( outgropt )
        curtree.start = curtree.nodep[outgrno - 1]->back;
      treevaluate();
      restml_printree();
      summarize();
      if (trout) {
        col = 0;
        restml_treeout(curtree.start);
      }
      clear_connections(&curtree,nonodes2); 

      which++;
    }
    FClose(intree);
    if (numtrees > 1 && weightsum > 1 )
    standev2(numtrees, maxwhich, 0, endsite-1, maxlogl, l0gl, l0gf,
              aliasweight, seed);
  } else {
    free_all_trans(&curtree);
    for (i = 1; i <= spp; i++)
      enterorder[i - 1] = i;
    if (jumble)
      randumize(seed, enterorder);
    if (progress) {
      printf("\nAdding species:\n");
      writename(0, 3, enterorder);
    }
    nextsp = 3;
    buildsimpletree(&curtree);
    curtree.start = curtree.nodep[enterorder[0] - 1]->back;
    smoothit = improve;
    nextsp = 4;
    while (nextsp <= spp) {
      buildnewtip(enterorder[nextsp - 1], &curtree);
      bestyet = - nextsp*sites*sitelength*log(4.0);
      if (smoothit)
        restml_copy_(&curtree, &priortree);
      addtraverse(curtree.nodep[enterorder[nextsp - 1] - 1]->back,
                  curtree.nodep[enterorder[0]-1]->back, true);
      if (smoothit)
        restml_copy_(&bestree, &curtree);
      else {
        smoothit = true;
        insert_(curtree.nodep[enterorder[nextsp - 1] - 1]->back, qwhere);
        for (i = 1; i<=smoothings; i++) {
          smooth (curtree.start);
          smooth (curtree.start->back);
        }
        smoothit = false;
        bestyet = curtree.likelihood;
      }
      if (progress)
        writename(nextsp - 1, 1, enterorder);
      if (global && nextsp == spp) {
        if (progress) {
          printf("Doing global rearrangements\n");
          printf("  !");
          for (j = spp ; j < nonodes2 ; j++)
            if ( (j - spp) % (( nonodes2 / 72 ) + 1 ) == 0 )
              putchar('-');
          putchar('!');
        }
      }
      succeeded = true;
      while (succeeded) {
        succeeded = false;
        if (global && nextsp == spp)
          globrearrange();
        else 
          rearrange(curtree.start, curtree.start->back);
      }
      nextsp++;
    }
    if (global && progress) {
      putchar('\n');
      fflush(stdout);
    }
    restml_copy_(&curtree, &bestree);
    if (njumble > 1) {
      if (jumb == 1)
        restml_copy_(&bestree, &bestree2);
      else
        if (bestree2.likelihood < bestree.likelihood)
          restml_copy_(&bestree, &bestree2);
    }
    if (jumb == njumble) {
      if (njumble > 1)
        restml_copy_(&bestree2, &curtree);
      curtree.start = curtree.nodep[outgrno - 1]->back;
      restml_printree();
      summarize();
      if (trout) {
        col = 0;
        restml_treeout(curtree.start);
      }
    }
  }
  if ( jumb < njumble )
    return;
  freex2(nonodes2, curtree.nodep);
   if (!usertree) {
     freex2(nonodes2, priortree.nodep);
     freex2(nonodes2, bestree.nodep);
     if (njumble > 1)
       freex2(nonodes2, bestree2.nodep);
   } else {
     free(l0gl);
     for (i=0;i<shimotrees;++i)
       free(l0gf[i]);
     free(l0gf);
   }
  if (jumb == njumble) {
    if (progress) {
      printf("\nOutput written to file \"%s\"\n\n", outfilename);
      if (trout)
        printf("Tree also written onto file \"%s\"\n", outtreename);
      putchar('\n');
    }
  }
}  /* maketree */


int main(int argc, Char *argv[])
{  /* maximum likelihood phylogenies from restriction sites */
#ifdef MAC
  argc = 1;                /* macsetup("Restml","");        */
  argv[0] = "RestML";
#endif
  init(argc, argv);
  progname = argv[0];
  openfile(&infile,INFILE,"input file","r",argv[0],infilename);
  openfile(&outfile,OUTFILE,"output file","w",argv[0],outfilename);
  ibmpc = IBMCRT;
  ansi = ANSICRT;
  mulsets = false;
  datasets = 1;
  firstset = true;
  doinit();
  if (trout)
    openfile(&outtree,OUTTREE,"output tree file","w",argv[0],outtreename);
  for (ith = 1; ith <= datasets; ith++) {
    if (datasets > 1) {
      fprintf(outfile, "Data set # %ld:\n",ith);
      if (progress)
        printf("\nData set # %ld:\n",ith);
    }
    getinput();
    if (ith == 1)
      firstset = false;
    for (jumb = 1; jumb <= njumble; jumb++)
      maketree();
  }
  FClose(infile);
  FClose(outfile);
  FClose(outtree);
#ifdef MAC
  fixmacfile(outfilename);
  fixmacfile(outtreename);
#endif
  printf("Done.\n\n");
  return 0;
}  /* maximum likelihood phylogenies from restriction sites */
