xref: /dports/biology/recombine/recombine1.41/recombine.c

#include "recombine.h"
#include "getdata.h"

#ifdef DMEMDEBUG
#define MEMDEBUG
#include "memdebug.h"
#endif

#ifdef DMALLOC_FUNC_CHECK
#include "/usr/local/include/dmalloc.h"
#endif

/***************************************************************************
 *  Version 0.40. (c) Copyright 1986, 1991, 1992 by the University of      *
 *  Washington and Joseph Felsenstein.  Written by Joseph Felsenstein,     *
 *  Mary K. Kuhner and Jon A. Yamato, with some additional grunt work by   *
 *  Sean T. Lamont.  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.                                                 *
 *                                                                         *
 ***************************************************************************/

/* This program implements a Hastings-Metropolis Monte Carlo
  Maximum Likelihood sampling method for phylogenetic trees
  with recombination. */

/* Time, 'tyme', in this tree is measured from the tips to the root.
   I.e. the tips are at tyme '0', and the root node has the largest
   value for 'tyme'. */

/* Modified to accept multiple loci. */

/* Modified to "plump" input non-recombinant tree */

/* Modified to track final coalescence and begin to handle microsats */

/* Modified to handle panel snps, and real FC differentiation */

/* Being modified to handle basic haplotyping */
/* WARNING MARY temporary variables */
data_fmt truehaps;
data_fmt starthaps;

/* these definitions will be extern in other files*/
extern double **reci, **theti;
extern treerec ***sum;
extern unsigned int baseA, baseC, baseG, baseT;

extern long countcoalbranches(tree *tr);
extern void model_alloc(option_struct *op, data_fmt *data);
extern void scoretree(option_struct *op, data_fmt *data, long chain);
extern void rec_estimate(option_struct *op, data_fmt *data, long lowcus,
   long chain, boolean locusend);
extern boolean twiddle(option_struct *op, data_fmt *data);
extern boolean fliphap(option_struct *op, data_fmt *data, tree *tr);
extern boolean flipdrop(option_struct *op, data_fmt *data, tree *tr,
   long numdrop);
extern void add_one_recomb(option_struct *op, data_fmt *data, tree *tr,
   long *siteptr, long chain);
extern void addfractrecomb(option_struct *op, data_fmt *data,
   long chain, treerec ***treesum);
extern boolean makedrop(option_struct *op, data_fmt *data);
extern void setupdnadata(dnadata **dna, long *numsites, long numseq,
   long numloci, long numpop);
extern void freednadata(dnadata *dna);
extern void empiricaldnafreqs(option_struct *op, data_fmt *data,
   tree *curtree);
extern void inputdnaweights(long numchars, dnadata *dna, option_struct *op);
extern void remove_branch(tlist *start, node *target);
extern void rename_branch(tlist *start, node *source, node *target);
extern void rec_scoreread(option_struct *op, data_fmt *data);
extern void unflag(node *p);
extern void traverse_unflag(tree *tr, node *p);
extern void read_spacefile(FILE *file, option_struct *op, data_fmt *data);
extern void printdnadata(option_struct *op, data_fmt *data, FILE *out);

extern long countsites(option_struct *op, data_fmt *data);
extern long getdata_sitecount(option_struct *op, data_fmt *data,
   long marker);
extern double getnumlinks(option_struct *op, data_fmt *data, long start,
   long finish);

/* functions for haplotyping */
extern void read_flipfile(option_struct *op, data_fmt *data, tree *tr,
   FILE *file);
extern void pruneflips(option_struct *op, data_fmt *data, tree *tr);
extern void copyhaps(option_struct *op, data_fmt *oldhap, data_fmt *newhap,
  long numpop, long numloci, long numind, long *numsites);

/* functions for heating */
extern double coalprob(option_struct *op, data_fmt *data, tree *tr, double
  theta, double r);
extern int longcmp(const void *v1, const void *v2);

/* functions from traitlike.c for trait likelihoods */
extern void copytraits(node *source, node *target);
extern void traitlike(option_struct *op, data_fmt *data, tree *tr,
   long numsites, double mutrait, double traitratio, double pd,
   double *traitarray);
extern void traitread(tree *tr, long numseq);
extern void traitprint(long numsites, double *traitarray, FILE *outfile, long numout);
extern void traitsiteplot(option_struct *op, data_fmt *data, double *llike);
extern void traitresult(option_struct *op, data_fmt *data,
   double *traitarray, FILE *outfile, long numout);

tree **temptrees;/* the current tree for each temperature-chain */
tree *curtree;   /* the current tree */
long locus;      /* which locus are we currently working on? */
long population; /* which population are we currently working on? */
longer seed;     /* array used to hold randum seed */
long numdropped = 0;               /* number of trees dropped from
                                      consideration */

/* next are a largish set of temporary variables used to hold various
   values before their structures are setup.  Mostly data variables. */
double freqa, freqc, freqg, freqt; /* dna freqs */
double locus_ttratio;              /* dna transition/transversion ratio */

FILE *infile, *outfile, *treefile, *bestree, *seedfile, *simlog,
     *parmfile, *intree, *weightfile, *spacefile;
long rootnum, apps, col, numtrees, totchains;
boolean  **sametree;
double clearseed, theta0, rec0, branch0;
long		*category;
contribarr	*contribution;
node		*freenodes;
rlrec		*alogf;
char            gch;
long            nodectr;
double	 	sumweightrat, watttheta;
double	 	*weightrat;
valrec   	*tbl;
long		slid, slacc, hap, hacc, swap, swacc, indecks;


/* the following are for reading in parameters (readparmfile),
   and also writing them back out again (rec_parmfilewrite) */
char *booltokens[NUMBOOL] = {"interleaved","printdata","progress",
                          "print-trees","freqs-from-data","categories",
                          "watterson", "usertree", "autocorrelation",
                          "newdata","same-ne","interactive","mhmcsave",
                          "panel","map","final-coalescence","full-snp",
                          "haplotyping","profile-like","norecsnp"},
     *numbertokens[NUMNUMBER] = {"ttratio","short-chains",
                          "short-steps","short-inc","long-chains",
                          "long-inc","long-steps","rec-rate",
                          "holding","mu-ratio","trait-ratio",
                          "dis-freq","hapdrop","heating"};

/* the following are for managing x array recycling */
long numx;
datalike_fmt *sparex[XARRAYSIZE];

void openfile(FILE **fp, char *filename, char *mode, char *application,
   char *perm)
{
  FILE *of;
  char file[100];

  strcpy(file,filename);
  while (1){
    of = fopen(file,mode);
    if (of)
      break;
    else {
      switch (*mode){
      case 'r':
        fprintf(stdout,"%s:  can't read %s\n",application,file);
        file[0] = '\0';
        while (file[0] =='\0'){
          fprintf(stdout,"Please enter a new filename>");
          fgets(file,100,stdin);
          }
        break;
      case 'w':
        fprintf(stdout,"%s: can't write %s\n",application,file);
        file[0] = '\0';
        while (file[0] =='\0'){
          fprintf(stdout,"Please enter a new filename>");
          fgets(file,100,stdin);
          }
        break;
      }
    }
  }
  *fp=of;
  if (perm != NULL)
    strcpy(perm,file);
} /* openfile */


char lowercase(char c)
{
return((char)tolower((int) c));
} /* lowercase */


double randum(void)
/* Mary's version--faster but needs 32 bits.  Loops have been unrolled
   for speed. */
{
  long newseed0, newseed1, newseed2;

  newseed0 = 1549*seed[0];
    newseed1 = newseed0/2048;
    newseed0 &= 2047;
    newseed2 = newseed1/2048;
    newseed1 &= 2047;
  newseed1 += 1549*seed[1]
              + 812*seed[0];
    newseed2 += newseed1/2048;
    newseed1 &= 2047;
  newseed2 += 1549*seed[2]
              + 812*seed[1];

  seed[0] = newseed0;
  seed[1] = newseed1;
  seed[2] = newseed2 & 1023;
  return (((seed[0]/2048.0 + seed[1])/2048.0 + seed[2])/1024.0);
}  /* randum */


/****************************************************
 * setupoption_struct() creates a new option_struct */
void setupoption_struct(option_struct **op)
{

(*op) = (option_struct *)calloc(1,sizeof(option_struct));

(*op)->rate = (double *)calloc(1,sizeof(double));
(*op)->probcat = (double *)calloc(1,sizeof(double));
(*op)->temperature = NULL;
(*op)->numpanel = NULL;

} /* setupoption_struct */


/********************************************************************
 * istip() returns TRUE if the nodelet is a tip and FALSE otherwise */
boolean istip (node *p)
{
return(p->type == 't');
} /* istip */


/***********************************************************************
 * iscoal returns TRUE if the nodelet is a coalescent, FALSE otherwise */
boolean iscoal (node *p)
{
return(p->type == 'c');
} /* iscoal */


/************************************************************************
 * isrecomb returns TRUE if the nodelet is recombinant, FALSE otherwise */
boolean isrecomb (node *p)
{
return(p->type == 'r');
} /* isrecomb */


/******************************************************************
 * branchsub() substitutes one branch for another in one interval *
 * in a tymelist.                                                 */
boolean branchsub(tlist *t, node *oldbranch, node *newbranch)
{
long i;
boolean found;

found = FALSE;
for(i = 0; i < t->numbranch; i++)
   if (t->branchlist[i] == oldbranch) {
      t->branchlist[i] = newbranch;
      found = TRUE;
   }

return(found);
} /* branchsub */


/************************************************************************
 * insertaftertymelist adds the node "p" to a tymelist after, inserting *
 *    within the tymeslice pointed to by "t" using bisection.           */
void insertaftertymelist(tlist *t, node *p)
{
long i, j, temp;
tlist *s;
node *q;
boolean b1, b2;

newtymenode(&s);
s->eventnode = findtop(p);
s->age = t->age;
t->age = p->tyme;
s->prev = t;
s->succ = t->succ;
t->succ = s;
if (s->succ != NULL) s->succ->prev = s;

/* now deal with the branchlist for the new tymelist entry;
   the "if" setting "temp" is a workaround for not being able to
   allocate zero size chunks of memory by malloc-debug/gcc */
if (t->numbranch-1 > 0) temp = t->numbranch - 1;
else temp = t->numbranch;
s->branchlist = (node **)calloc(1,temp*sizeof(node *));
s->numbranch = 0;
q = findunique(p);
if (isrecomb(q)) { /* if we've added a recombinant node */
   for (i = 0; i < t->numbranch; i++) {
      if (t->branchlist[i] == q->back) continue;
      s->branchlist[s->numbranch] = t->branchlist[i];
      s->numbranch++;
   }
   s->branchlist =
      (node **)realloc(s->branchlist, (s->numbranch+2)*sizeof(node *));
   s->branchlist[s->numbranch] = q->next;
   s->branchlist[s->numbranch+1] = q->next->next;
   s->numbranch += 2;
} else { /* else we've added a coalescent node */
   for (i = 0; i < t->numbranch; i++)
     if(t->branchlist[i] != q->next->back &&
        t->branchlist[i] != q->next->next->back) {
        s->branchlist[s->numbranch] = t->branchlist[i];
        s->numbranch++;
     } else if (t->branchlist[i] == q->next->back) {
        s->branchlist[s->numbranch] = q;
        s->numbranch++;
     }
}

/* now remove the old branches from the rest of the tymelist */
if (isrecomb(q)) {
   for (s = s->succ; s != NULL; s = s->succ) {
      if (!branchsub(s,q->back,q->next)) break;
      s->branchlist =
         (node **)realloc(s->branchlist, (s->numbranch+1)*sizeof(node *));
      s->branchlist[s->numbranch] = q->next->next;
      s->numbranch++;
   }
} else {
   for (s = s->succ; s != NULL; s = s->succ) {
      b1 = branchsub(s,q->next->back,q);
      b2 = branchsub(s,q->next->next->back,q);
      if (b1 && b2) {
         for (i = 0; i < t->numbranch; i++)
            if (s->branchlist[i] == q) break;
         for (j = i; j < t->numbranch-1; j++)
            s->branchlist[j] = s->branchlist[j+1];
         s->numbranch--;
      }
      if (!b1 && !b2) break;
   }
}

} /* insertaftertymelist */


/******************************************************************
 * subtymelist() completely removes the tymelist entry passed in. *
 * the passed in branch should be the exact branch that wants to  *
 * be removed in the case of a recombinant node, or the branch    *
 * that will "replace" the missing branch in a coalescent node.   */
void subtymelist(tlist *t, node *branch, boolean all)
{

if (isrecomb(t->eventnode)) {
   remove_branch(t,branch);
   if (all) remove_branch(t,otherparent(branch));
   else rename_branch(t,findunique(t->eventnode)->back,otherparent(branch));
} else
   if (all) remove_branch(t,t->eventnode);
   else rename_branch(t,branch,t->eventnode);

if (t->prev) {t->prev->age = t->age; t->prev->succ = t->succ;}
else {fprintf(ERRFILE,"ERROR:tried to remove first tymelist entry!\n");}
if (t->succ) t->succ->prev = t->prev;

freetymenode(t);

} /* subtymelist */


/****************************************************************
 * printtymelist() prints the entire tymelist: a debug function */
void printtymelist(tlist *t)
{
  long i;

  fprintf(ERRFILE,"TYMELIST BEGINS\n");
  while (t != NULL) {
    fprintf(ERRFILE,"%3ld age%8.6f--",
           t->eventnode->number, t->age);
    for (i = 0; i < t->numbranch; i++) {
      fprintf(ERRFILE," %3ld", t->branchlist[i]->number);
      if (t->branchlist[i]->top) fprintf(ERRFILE,"t ");
    }
    fprintf(ERRFILE,"\n");
    t = t->succ;
  }
  fprintf(ERRFILE,"TYMELIST ENDS\n");
}  /* printtymelist */


/*****************************************************************
 * lengthof() returns the length of the branch "rootward" of the *
 * passed node                                                   */
double lengthof(node *p)
{
  return fabs(p->tyme - p->back->tyme);
}  /* lengthof */


/******************************************************************
 * fixlength() sets the length and v values for a mother-daughter *
 * pair of nodelets.                                              */
void fixlength(option_struct *op, data_fmt *data, node *p)
{
double newlength;

newlength = lengthof(p);
p->length = newlength;
ltov(op,data,p);
p->back->length = newlength;
ltov(op,data,p->back);

} /* fixlength */


/************************************************
 * findtop() finds a "top" nodelet in the node. */
node *findtop(node *p)
{

if (!istip(p))
   while (!p->top) p = p->next;

return(p);

} /* findtop */


/********************************************************************
 * findunique() returns the unique nodelet of a node (i.e. "bottom" *
 * of a recombinant node or "top" of a coalescent node).            */
node *findunique(node *p)
{

if (istip(p)) return(p);

if (isrecomb(p))
   while((p)->top) p = p->next;
else
   while (!p->top) p = p->next;

return(p);

}  /* findunique */


/*********************************************************************
 * otherdtr() finds the other daughter nodelet of a node, it assumes *
 * that a daughter nodelet was passed to it.                         */
node *otherdtr(node *p)
{
node *q;

if (isrecomb(p)) {
   fprintf(ERRFILE,"ERROR:otherdtr() passed a recombinant node!\n");
   return(NULL);
}

if (p->top) {
   fprintf(ERRFILE,"ERROR:otherdtr() passed a top nodelet!\n");
   return(NULL);
}

for(q = p->next; q->top; q = q->next)
   ;

return(q);

} /* otherdtr */


/*********************************************************************
 * findlink() takes a recombinant node and returns the # of the site *
 * to the left of the cut link; or FLAGLONG in error states.         */
long findlink(node *p)
{
node *q;

if (!isrecomb(p)) {
   fprintf(ERRFILE,"ERROR:findlink passed non-recombinant node\n");
   return(FLAGLONG);
}

q = findunique(p);
if (q->next->recstart) return(q->next->next->recend);
else return(q->next->recend);

} /* findlink */


/**********************************************************************
 * otherparent() finds the other parent nodelet of a node, it assumes *
 * that a parent nodelet was passed to it.                            */
node *otherparent(node *p)
{
node *q;

if (!isrecomb(p)) {
   fprintf(ERRFILE,"ERROR:otherparent() passed a non-recombinant node!\n");
   return(NULL);
}

if (!p->top) {
   fprintf(ERRFILE,"ERROR:otherparent() passed a non-top nodelet!\n");
   return(NULL);
}

for(q = p->next; !q->top; q = q->next)
   ;

return(q);

} /* otherdtr */


/*********************************************
 * free_z() frees the "z" field of a nodelet */
void free_z(option_struct *op, node *p)
{
if (p->z == NULL) return;

free(p->z);
p->z = NULL;

} /* free_z */


/***************************************************************
 * allocate_z() allocates space for the "z" field of a nodelet */
void allocate_z(option_struct *op, data_fmt *data, node *p)
{
if (p->z) free(p->z);

p->z = (double *)calloc(NUMTRAIT,sizeof(double));

} /* allocate_z */


/*********************************************
 * free_x() frees the "x" field of a nodelet */
void free_x(option_struct *op, node *p)
{

if (p->x == NULL) return;

if (numx < XARRAYSIZE) {
   sparex[numx] = p->x;
   numx++;
} else {
   switch(op->datatype) {
      case 'a':
         break;
      case 'b':
      case 'm':
         free(p->x->a);
         break;
      case 'n':
      case 's':
         free(p->x->s[0][0][0]);
         free(p->x->s[0][0]);
         free(p->x->s[0]);
         free(p->x->s);
         break;
      default:
         fprintf(ERRFILE,"ERROR:free_x: can't get here!\n");
         exit(-1);
   }
   free(p->x);
}

p->x = NULL;

} /* free_x */


/************************************************************
 * allocate_x() allocates space for the "x" field of a nodelet */
void allocate_x(option_struct *op, data_fmt *data, node *p)
{
long i, j, k, nummarkers = 0, numloci, numslice;

if (p->x != NULL) free_x(op,p);

if(numx > 0) {
   p->x = sparex[numx-1];
   sparex[numx-1] = NULL;
   numx--;
} else {
   p->x = (datalike_fmt *)calloc(1,sizeof(datalike_fmt));
   switch(op->datatype) {
      case 'a':
         p->x->a = NULL;
         p->x->s = NULL;
         break;
      case 'b':
      case 'm':
         numloci = getdata_numloci(op,data);
         p->x->s = NULL;
         p->x->a =
            (double **)calloc(numloci,sizeof(double *));
         p->x->a[0] = (double *)calloc(numloci*MICRO_ALLELEMAX,
            sizeof(double));
         for(i = 1; i < numloci; i++)
            p->x->a[i] = p->x->a[0] + i*MICRO_ALLELEMAX;
         break;
      case 'n':
         nummarkers += NUMINVAR;
      case 's':
         numslice = (op->panel) ? NUMSLICE : 1L;
         nummarkers += getdata_nummarkers(op,data);
         p->x->a = NULL;
         p->x->s = (double ****)calloc(nummarkers,sizeof(double ***));
         p->x->s[0] = (double ***)
            calloc(nummarkers*op->categs,sizeof(double **));
         for(i = 1; i < nummarkers; i++)
            p->x->s[i] = p->x->s[0] + i*op->categs;
         p->x->s[0][0] = (double **)
            calloc(nummarkers*op->categs*numslice,sizeof(double *));
         for(i = 0; i < nummarkers; i++)
            for(j = 0; j < op->categs; j++)
               p->x->s[i][j] = p->x->s[0][0] + i*op->categs*numslice +
                               j*numslice;
         p->x->s[0][0][0] = (double *)
            calloc(nummarkers*op->categs*numslice*4L,sizeof(double));
         for(i = 0; i < nummarkers; i++)
            for(j = 0; j < op->categs; j++)
               for(k = 0; k < numslice; k++)
                  p->x->s[i][j][k] = p->x->s[0][0][0] +
                                  i*op->categs*numslice*4L +
                                  j*numslice*4L + k*4L;
         break;
      default:
         fprintf(ERRFILE,"ERROR:alloc_x: can't get here!\n");
         exit(-1);
   }
}

} /* allocate_x */


/************************************************************
 * init_coal_alloc() callocs and initializes a range field. */
void init_coal_alloc(long **coal, long numcoalpairs)
{

if (*coal != NULL) free(*coal);

(*coal) = (long *)calloc(numcoalpairs * 2 + 2,sizeof(long));
(*coal)[0] = numcoalpairs;
(*coal)[numcoalpairs*2+1] = FLAGLONG;

} /* init_coal_alloc */


/****************************************************************
 * coal_Malloc() callocs or frees the "coal" field of a nodelet */
void coal_Malloc(node *p, boolean allokate, long numcoalpairs)
{

if (allokate && p->top) {
   init_coal_alloc(&(p->coal),numcoalpairs);
} else {
   if (p->coal != NULL) {
      free(p->coal);
      p->coal = NULL;
   }
}

} /* coal_Malloc */


/**************************************************************
 * init_ranges_alloc() callocs and initializes a range field. */
void init_ranges_alloc(long **ranges, long numrangepairs)
{

if (*ranges != NULL) free(*ranges);

(*ranges) = (long *)calloc(numrangepairs * 2 + 2,sizeof(long));
(*ranges)[0] = numrangepairs;
(*ranges)[numrangepairs*2+1] = FLAGLONG;

} /* init_ranges_alloc */


/********************************************************************
 * ranges_Malloc() callocs or frees the "ranges" field of a nodelet */
void ranges_Malloc(node *p, boolean allokate, long numrangepairs)
{

if (allokate && p->top) {
   init_ranges_alloc(&(p->ranges),numrangepairs);
} else {
   if (p->ranges != NULL) {
      free(p->ranges);
      p->ranges = NULL;
   }
}

} /* ranges_Malloc */


/****************************************************************
 * meld_adjacent_ranges() concatenates elements adjacent to the *
 * passed position if it is appropiate to do so.                */
void meld_adjacent_ranges(long *cranges, long newelem)
{
long numpairs, nummove, prevelem, nextelem, lastelem;

prevelem = newelem - 1;
nextelem = newelem + 2;
lastelem = 2*cranges[0]+1;

numpairs = 0;
if (cranges[newelem] == cranges[prevelem]+1) numpairs++;
if (cranges[newelem+1]+1 == cranges[nextelem]) {
   if (numpairs) numpairs++;
   else numpairs = -1L;
   }
if (numpairs == 2 && prevelem == 0) numpairs = -1L;
if (numpairs == 2 && nextelem == lastelem) numpairs--;

if (numpairs == 1 && prevelem != 0) {
   nummove = 2*cranges[0] - newelem;
   cranges[prevelem] = cranges[newelem+1];
   memmove(&cranges[newelem],&cranges[nextelem],nummove*sizeof(long));
   cranges[0]--;
}
if (numpairs == -1 && nextelem != lastelem) {
   nummove = 2*cranges[0] - newelem - 2;
   cranges[newelem+1] = cranges[nextelem+1];
   memmove(&cranges[nextelem],&cranges[nextelem+2],nummove*sizeof(long));
   cranges[0]--;
}
if (numpairs == 2) {
   nummove = 2*cranges[0] - newelem - 2;
   cranges[newelem-1] = cranges[nextelem+1];
   memmove(&cranges[newelem],&cranges[nextelem+2],nummove*sizeof(long));
   cranges[0] -= 2;
}

} /* meld_adjacent_ranges */


/****************************************************
 * addrange() adds a new element to the range field */
void addrange(long **nnewranges, long newstart, long newend)
{
long i, first, numpairs, nummove, *newranges;

if (newstart > newend)
   fprintf(ERRFILE,"ERROR:addrange--start > end\n");

newranges = (*nnewranges);
numpairs = newranges[0]*2;

/* was there nothing in the initial range list? */
if(!numpairs) {
   newranges[0]++;
   newranges = (long *)realloc(newranges,4*sizeof(long));
   (*nnewranges) = newranges;
   newranges[1] = newstart;
   newranges[2] = newend;
   newranges[3] = FLAGLONG;
   return;
}

/* do I start after everything else ends? */
if (newstart > newranges[numpairs]) {
   if (newstart == newranges[numpairs]+1) {
      newranges[numpairs] = newend;
   } else {
      newranges[0]++; numpairs += 2;
      newranges = (long *)realloc(newranges,(numpairs+2)*sizeof(long));
      (*nnewranges) = newranges;
      newranges[numpairs-1] = newstart;
      newranges[numpairs] = newend;
      newranges[numpairs+1] = FLAGLONG;
   }
   return;
}

/* do I end before anything begins? */
if (newend < newranges[1]) {
   if (newend == newranges[1]-1) {
      newranges[1] = newstart;
   } else {
      newranges[0]++;
      newranges =
         (long *)realloc(newranges,(newranges[0]*2+2)*sizeof(long));
      (*nnewranges) = newranges;
      memmove(&newranges[3],&newranges[1],(numpairs+1)*sizeof(long));
      newranges[1] = newstart;
      newranges[2] = newend;
   }
   return;
}

/* am I contained in another interval altogether? OR
   do I overlap one or more intervals?  OR
   am I between two already existing intervals? */
for(i = 1, numpairs = 0, first = 0; newranges[i] != FLAGLONG; i+=2) {
   if(newstart >= newranges[i] && newend <= newranges[i+1]) return;
   if (i != 1) {
      if (newstart > newranges[i-1] && newend < newranges[i]) {
         first = i;
         break;
      }
   }
   if((newranges[i] >= newstart && newranges[i] <= newend) ||
      (newranges[i+1] >= newstart && newranges[i+1] <= newend)) {
      numpairs++;
      if (!first) first = i;
   }
}

/* do I exist between two intervals? */
if (!numpairs) {
   nummove = 2*newranges[0]+2 - first;
   newranges[0]++;
   newranges =
      (long *)realloc(newranges,(newranges[0]*2+2)*sizeof(long));
   (*nnewranges) = newranges;
   memmove(&newranges[first+2],&newranges[first],nummove*sizeof(long));
   newranges[first] = newstart;
   newranges[first+1] = newend;
   meld_adjacent_ranges(newranges,first);
   return;
}

/* I must overlap one or more intervals */
if (numpairs == 1) {
   newranges[first] =
      (newranges[first] < newstart) ? newranges[first] : newstart;
   newranges[first+1] =
      (newranges[first+1] > newend) ? newranges[first+1] : newend;
} else {
   nummove = 2*newranges[0]+2 - first - 2*numpairs;
   newranges[0] -= numpairs-1;
   newranges[first] =
      (newranges[first] < newstart) ? newranges[first] : newstart;
   newranges[first+1] =
      (newranges[first+1+2*(numpairs-1)] > newend) ?
         newranges[first+1+2*(numpairs-1)] : newend;
   memmove(&newranges[first+2],&newranges[first+2*numpairs],
           nummove*sizeof(long));
}
meld_adjacent_ranges(newranges,first);

} /* addrange */


/******************************************************************
 * inrange() checks to see if a site is active on a given branch */
boolean inrange(long *ranges, long site)
{
long i;

for(i = 1; ranges[i] != FLAGLONG; i+=2) {
   if (ranges[i] > site) return(FALSE);
   if (ranges[i+1] >= site) return(TRUE);
}

return(FALSE);

} /* inrange */


/************************************************************
 * subrangefc() causes the passed range to be set to "dead" */
void subrangefc(long **newranges, long substart, long subend)
{
long i, *nranges, numpairs, first, nummove;

nranges = (*newranges);

/* was there nothing in the initial range list? */
/* do I start after everything else ends? */
/* do I end before anything begins? */
if (!nranges[0] || substart > nranges[2*nranges[0]] ||
    subend < nranges[1]) return;

/* am I contained in another interval altogether? OR
   do I overlap one or more intervals?  OR
   am I between two already existing intervals? */
for(i = 1, numpairs = 0, first = 0; nranges[i] != FLAGLONG; i+=2) {
   if(nranges[i] > subend) break;
   if (i != 1) {
      if (substart > nranges[i-1] && subend < nranges[i])
         return;
   }

/*    if the deletion starts in this range OR
            the deletion ends in this range OR
            the deletion overlaps this range THEN */
   if((substart >= nranges[i] && substart <= nranges[i+1]) ||
      (subend >= nranges[i] && subend <= nranges[i+1]) ||
      (substart <= nranges[i] && subend >= nranges[i+1])) {
      numpairs++;
      if (!first) first = i;
   }
}

/* am I contained within a single interval? */
if (numpairs == 1) {
/*    should I just remove the whole entry? */
   if ((substart <= nranges[first]) && (subend >= nranges[first+1])) {
      nummove = 2*nranges[0] - first;
      memmove(&nranges[first],&nranges[first+2],nummove*sizeof(long));
      nranges[0]--;
      return;
   }
/*    or chop off the beginning? */
   if (substart <= nranges[first]) {
      nranges[first] = subend+1;
      return;
   }
/*    or chop off the end? */
   if (subend >= nranges[first+1]) {
      nranges[first+1] = substart-1;
      return;
   }
/*    then split it in two! */
   nranges[0]++;
   nranges =
      (long *)realloc(nranges,(2*nranges[0]+2)*sizeof(long));
   (*newranges) = nranges;
   nummove = 2*nranges[0] - (first+1);
   memmove(&nranges[first+3],&nranges[first+1],nummove*sizeof(long));
   nranges[first+1] = substart-1;
   nranges[first+2] = subend+1;
   return;
}

/* I must overlap one or more intervals */
if (substart > nranges[first]) nranges[first+1] = substart-1;
if (subend < nranges[first+1+2*(numpairs-1)])
   nranges[first+2*(numpairs-1)] = subend + 1;
else {
/*   remove the last affected rangepair */
   nummove = 2*nranges[0] - (first+2*(numpairs-1));
   memmove(&nranges[first+2*(numpairs-1)],&nranges[first+2*(numpairs)],
      nummove*sizeof(long));
   nranges[0]--;
}
/*  remove the middle affected rangepairs, if any */
nummove = 2*nranges[0]+2 - (first+2*(numpairs-1));
memmove(&nranges[first+2],&nranges[first+2*(numpairs-1)],
   nummove*sizeof(long));
nranges[0] -= numpairs-2;
/*  remove the first affected rangepair, if necessary */
if (substart <= nranges[first]) {
   nummove = 2*nranges[0] - first;
   memmove(&nranges[first],&nranges[first+2],nummove*sizeof(long));
   nranges[0]--;
}

} /* subrangefc */


void printrange(long *ranges)
{
long i;

printf("\n%ld range pairs:  ",ranges[0]);
for(i = 1; ranges[i] != FLAGLONG; i+=2)
   printf("%ld to %ld ",ranges[i],ranges[i+1]);
printf("\n");

} /* printrange */


/* "newnode" & "freenode" are paired memory managers for tree nodes.
   They maintain a linked list of unused nodes which should be faster
   to use than "calloc" & "free" (at least for recombination).
   They can only be used for internal nodes, not tips. */
void newnode(node **p)
{
  long i;
  node *q;

  if (freenodes == NULL) { /* need a new node */
    (*p) = allocate_nodelet(3L,'i');
    (*p)->number = nodectr;
    for(q = (*p)->next; q != (*p); q = q->next) q->number = nodectr;
    nodectr++;
    return;
  } else { /* recycle an old node */
    q=freenodes;
    freenodes=freenodes->back;
    for(i = 0; i < 3; i++) {
      ranges_Malloc(q,FALSE,0L);
      coal_Malloc(q,FALSE,0L);
      q->updated = FALSE;
      q->futileflag = 0L;
      q->recstart = q->recend = -1L;
      q->back = NULL;
      q = q->next;
    }
    *p = q;
    return;
  }
}  /* newnode */

void freenode(node *p)
{
   p->back = freenodes;
   freenodes = p;
}  /* freenode */
/* END of treenode allocation functions */

void newtymenode(tlist **t)
{
  *t = (tlist *)calloc(1,sizeof(tlist));
  (*t)->prev = NULL;
  (*t)->succ = NULL;
}  /* newtymenode */

void freetymenode(tlist *t)
{
  if(t->branchlist!=NULL) free(t->branchlist);
  free(t);
}  /* freetymenode*/

void freetymelist(tlist *t)
{
   if (t->succ != NULL) freetymelist(t->succ);

   freetymenode(t);

} /* freetymelist */

void hookup(node *p, node *q)
{
  p->back = q;
  q->back = p;
}  /* hookup */

void atr(node *p)
/* "atr" prints out a text representation of a tree.  Pass
   curtree->root->back for normal results.  This is a debugging
   function which is not normally called. */
{
long i;
node *q;

  if (p->back == curtree->root) {
     fprintf(ERRFILE,"next node is root\n");
     fprintf(ERRFILE,"Node %4ld length %12.6f tyme %10.8f",
             p->back->number, lengthof(p), p->back->tyme);
     fprintf(ERRFILE," --> %4ld\n",p->number);
  }
  fprintf(ERRFILE,"Node %4ld update %ld length %10.8f tyme %10.8f -->",
         p->number, (long)p->updated, lengthof(p), p->tyme);
  if (!istip(p))
     for(i = 0, q = p; i < 3; i++, q = q->next)
        if (!q->top) fprintf(ERRFILE,"%4ld\n",q->back->number);
  fprintf(ERRFILE,"\n");
  if (p->top && p->back->top) fprintf(ERRFILE,"TWO TOPS HERE!!!!");
  if (!istip(p)) {
     if (!p->next->top) atr(p->next->back);
     if (!p->next->next->top) atr(p->next->next->back);
  }
  else fprintf(ERRFILE,"\n");
} /* atr */

void probatr(node *p)
/* "probatr" checks the internal representation of the tree.  Pass
   curtree->root->back for normal results.  This is a debugging
   function which is not normally called. */
{
long i;

if (p->back == curtree->root) {
}

if(p->top) {
   if(p->ranges[0] < 0)
      fprintf(ERRFILE,"node %ld has less then 0 ranges %ld %ld\n",
      p->number,indecks,apps);
   if(p->ranges[0] > 2)
      fprintf(ERRFILE,"node %ld has %ld ranges %ld %ld\n",
      p->number,p->ranges[0],indecks,apps);
   if(p->ranges[2*p->ranges[0]+1] != FLAGLONG)
      fprintf(ERRFILE,"node %ld has misformed end of ranges %ld %ld\n",
      p->number,indecks,apps);
   for(i = 1; p->ranges[i] != FLAGLONG; i+=2) {
      if(p->ranges[i] > p->ranges[i+1])
         fprintf(ERRFILE,"node %ld has begins > ends %ld %ld\n",
         p->number,indecks,apps);
      if(p->ranges[i] < 0 || p->ranges[i+1] < 0)
         fprintf(ERRFILE,"node %ld has negative starts or ends %ld %ld\n",
         p->number,indecks,apps);
   }
}

if (p->top && p->back->top) fprintf(ERRFILE,"TWO TOPS HERE!!!!");
if (!istip(p)) {
   if (!p->next->top) probatr(p->next->back);
   if (!p->next->next->top) probatr(p->next->next->back);
}


} /* probatr */


/***************************************************************
 * gettymenode() returns a pointer to the tymelist entry whose *
 * 'eventnode' has the number of 'target'.                     */
tlist *gettymenode(tree *tr, long target)
{
boolean found;
tlist *t;

if (target == tr->root->number) return (NULL);
if (tr->nodep[target]->type == 't') return(tr->tymelist);

for(t = tr->tymelist, found = FALSE; t != NULL; t = t->succ)
   if (t->eventnode->number == target) {found = TRUE; break;}

if (!found) {
   fprintf(ERRFILE,"ERROR:In gettymenode, failed to find node %12ld ", target);
   fprintf(ERRFILE,"CATASTROPHIC ERROR\n");
   exit(-1);
}

return(t);

}  /* gettymenode */


/***************************************************************
 * vtol() takes a "v" value for a branchlength and returns the *
 * branchlength.                                               */
double vtol(option_struct *op, data_fmt *data, double v)
{

switch(op->datatype) {
   case 'a':
      break;
   case 'b':
   case 'm':
      return(v);
      break;
   case 'n':
   case 's':
      return(data->dnaptr->fracchange * -1.0 * v);
      break;
   default:
      fprintf(ERRFILE,"ERROR:vtol, can't get here!\n");
      exit(-1);
}

return(-1.0);

} /* vtol */


/******************************************************************
 * ltov() recalculates the proper "v" value of a branch, from the *
 * tymes at either end of the branch.                             */
void ltov(option_struct *op, data_fmt *data, node *p)
{

switch(op->datatype) {
   case 'a':
      break;
   case 'b':
   case 'm':
      p->v = lengthof(p);
      break;
   case 'n':
   case 's':
      p->v = -1.0 * lengthof(p) / data->dnaptr->fracchange;
      break;
   default:
      fprintf(ERRFILE,"ERROR:ltov, can't get here!\n");
      exit(-1);
}
p->back->v = p->v;

}  /* ltov */


/***************************************************************
 * findcoal_ltov() is a duplicate ltov specially for findcoal. */
double findcoal_ltov(option_struct *op, data_fmt *data, double value)
{

switch(op->datatype) {
   case 'a':
      break;
   case 'b':
   case 'm':
      return(value);
      break;
   case 'n':
   case 's':
      return(-1.0 * value / data->dnaptr->fracchange);
      break;
   default:
      fprintf(ERRFILE,"ERROR:ltov, can't get here!\n");
      exit(-1);
}

return(-1.0);

} /* findcoal_ltov */


/* boolcheck(), booleancheck(), numbercheck(), and readparmfile()
   are used in reading the parameter file "parmfile" */
long boolcheck(char ch)
{
   ch = toupper((int)ch);
   if (ch == 'F' || ch == 'N') return 0;
   if (ch == 'T' || ch == 'Y') return 1;
   return -1;
} /* boolcheck */

boolean booleancheck(option_struct *op, char *var, char *value)
{
   long i, j, check;

   check = boolcheck(value[0]);
   if(check == -1) return FALSE;

   for(i = 0; i < NUMBOOL; i++) {
      if(!strcmp(var,booltokens[i])) {
         if(i == 0) op->interleaved = (boolean)(check);
         if(i == 1) op->printdata = (boolean)(check);
         if(i == 2) op->progress = (boolean)(check);
         if(i == 3) op->treeprint = (boolean)(check);
         if(i == 4) {
            op->freqsfrom = (boolean)(check);
            if(!op->freqsfrom) {
               strtok(value,":");
               freqa = (double)atof(strtok(NULL,";"));
               freqc = (double)atof(strtok(NULL,";"));
               freqg = (double)atof(strtok(NULL,";"));
               freqt = (double)atof(strtok(NULL,";"));
            }
         }
         if(i == 5) {
            op->ctgry = (boolean)(check);
            if(op->ctgry) {
               strtok(value,":");
               op->categs = (long)atof(strtok(NULL,";"));
               op->rate =
                  (double *)realloc(op->rate,op->categs*sizeof(double));
               op->probcat =
                  (double *)realloc(op->probcat,op->categs*sizeof(double));
               for(j = 0; j < op->categs; j++) {
                  op->rate[j] = (double)atof(strtok(NULL,";"));
                  op->probcat[j] = (double)atof(strtok(NULL,";"));
               }
            }
         }
         if(i == 6) {
            op->watt = (boolean)(check);
            if (!op->watt) {
               strtok(value,":");
               theta0 = (double)atof(strtok(NULL,";"));
            }
         }
         if(i == 7) op->usertree = (boolean)(check);
         if(i == 8) {
            op->autocorr = (boolean)(check);
            if (op->autocorr) {
               strtok(value,":");
               op->lambda = 1.0 / (double)atof(strtok(NULL,";"));
            }
         }
         if(i == 9) op->newdata = (boolean)(check);
         if(i == 10) {
            op->same_ne = (boolean)(check);
            if (!op->same_ne) {
               strtok(value,":");
               check = atol(strtok(NULL,";"));
               op->ne_ratio = (double *)calloc(check,sizeof(double));
               for(j = 0; j < check; j++)
                  op->ne_ratio[j] = (double)atof(strtok(NULL,";"));
            }
         }
         if(i == 11) op->interactive = (boolean)(check);
         if(i == 12) op->mhmcsave = (boolean)(check);
         if(i == 13) {
             op->panel = (boolean)(check);
             if (op->panel) {
                strtok(value,":");
                check = atol(strtok(NULL,";"));
                op->numpanel = (long *)calloc(check,sizeof(long));
                for (j = 0; j < check; j++)
                   op->numpanel[j] = atol(strtok(NULL,";"));
             }
         }
         if(i == 14) op->map = (boolean)(check);
         if(i == 15) op->fc = (boolean)(check);
         if(i == 16) {
            op->full = (boolean)(check);
            if (op->full) {
               strtok(value,":");
               op->chance_seen = (double)atof(strtok(NULL,";"));
            }
         }
         if(i == 17) {
            op->haplotyping = (boolean)(check);
            if (op->haplotyping) {
               strtok(value,":");
               op->happrob = (double)atof(strtok(NULL,";"));
            }
         }
         if(i == 18) op->profile = (boolean)(check);
         if(i == 19) op->norecsnp = (boolean)(check);
         return TRUE;
      }
   }
   return FALSE;
} /* booleancheck */

boolean numbercheck(option_struct *op, char *var, char *value)
{
long i, j;

for(i = 0; i < NUMNUMBER; i++) {
   if(!strcmp(var,numbertokens[i])) {
      if(i == 0) locus_ttratio = (double)atof(value);
      if(i == 1) op->numchains[0] = atol(value);
      if(i == 2) op->steps[0] = atol(value);
      if(i == 3) op->increm[0] = atol(value);
      if(i == 4) op->numchains[1] = atol(value);
      if(i == 5) op->increm[1] = atol(value);
      if(i == 6) op->steps[1] = atol(value);
      if(i == 7) rec0 = (double)atof(value);
      if(i == 8) {
         value[0] = toupper((int)value[0]);
         switch(value[0]) {
             case 'N':
                op->holding = 0;
                break;
             case 'T':
                op->holding = 1;
                break;
             case 'R':
                op->holding = 2;
                break;
             default:
                fprintf(ERRFILE,"unknown setting for holding: %s,",
                   value);
                fprintf(ERRFILE,"  no holding assumed\n");
                op->holding = 0;
                break;
         }
      }
      if(i == 9) op->mutrait = (double)atof(value);
      if(i == 10) op->traitratio = (double)atof(value);
      if(i == 11) op->pd = (double)atof(value);
      if(i == 12) op->hapdrop = atol(value);
      if(i == 13) {
         op->numtempchains = atol(strtok(value,";"));
         if (!op->temperature) free(op->temperature);
         op->temperature = (long *)calloc(op->numtempchains,sizeof(long));
         for(j = 0; j < op->numtempchains; j++)
            op->temperature[j] = atof(strtok(NULL,";"));
      qsort((void *)(op->temperature),op->numtempchains,sizeof(long),longcmp);
      }
      return TRUE;
   }
}
return FALSE;

} /* numbercheck */

void readparmfile(option_struct *op)
{
char fileline[LINESIZE],parmvar[LINESIZE],varvalue[LINESIZE];

#ifdef MAC
char parmfilename[100];

strcpy(parmfilename,"parmfile");
#endif

parmfile = fopen("parmfile","r");

if(parmfile) {
   while(fgets(fileline, LINESIZE, parmfile) != NULL) {
      if(fileline[0] == '#') continue;
      if(!strncmp(fileline,"end",3)) break;
      strcpy(parmvar,strtok(fileline,"="));
      strcpy(varvalue,strtok(NULL,"\n"));
      /* now to process... */
      if(!booleancheck(op,parmvar,varvalue))
         if(!numbercheck(op,parmvar,varvalue)) {
            fprintf(ERRFILE,
               "Inappropiate entry in parmfile: %s\n", fileline);
         }
   }
} else
   if (!MENU) {
      fprintf(simlog,"Parameter file (parmfile) missing\n");
      exit(-1);
   }

FClose(parmfile);

#ifdef MAC
fixmacfile(parmfilename);
#endif

} /* readparmfile */
/* END parameter file read */


boolean whichopbool(option_struct *op, long i)
{
if (i == 0) return(op->interleaved);
if (i == 1) return(op->printdata);
if (i == 2) return(op->progress);
if (i == 3) return(op->treeprint);
if (i == 4) return(op->freqsfrom);
if (i == 5) return(op->ctgry);
if (i == 6) return(op->watt);
if (i == 7) return(op->usertree);
if (i == 8) return(op->autocorr);
if (i == 9) return(op->newdata);
if (i == 10) return(op->same_ne);
if (i == 11) return(op->interactive);
if (i == 12) return(op->mhmcsave);
if (i == 13) return(op->panel);
if (i == 14) return(op->map);
if (i == 15) return(op->fc);
if (i == 16) return(op->full);
if (i == 17) return(op->haplotyping);
if (i == 18) return(op->profile);
if (i == 19) return(op->norecsnp);

fprintf(ERRFILE,"Warning:  whichopbool failed to identify option\n");
return(FALSE);
} /* whichopbool */


void rec_parmfilewrite(option_struct *op, long numloci, long numpop)
{
long i, j;
boolean b;
char parmfilename[100];


openfile(&parmfile,"parmfile","w+",NULL,parmfilename);

for(i = 0; i < NUMBOOL; i++) {
   b = whichopbool(op,i);
   fprintf(parmfile,"%s=%s",booltokens[i],BOOLPRINT(b));
   if (i == 4 && !b) { /* freqsfrom */
      fprintf(parmfile,":%f;%f;%f;%f;",
         freqa, freqc, freqg, freqt);
   }
   if (i == 5 && b) { /* category */
      fprintf(parmfile,":%ld;",op->categs);
      for(j = 0; j < op->categs; j++)
         fprintf(parmfile,"%f;%f;",op->rate[j],op->probcat[j]);
   }
   if (i == 6 && !b) { /* starting theta */
      fprintf(parmfile,":%f",theta0);
   }
   if (i == 8 && !b) { /* autocorrelation */
      fprintf(parmfile,":%f",1.0/op->lambda);
   }
   if (i == 10 && !b) { /* same_ne */
      fprintf(parmfile,":%ld;",numloci);
      for(j = 0; j < numloci; j++)
         fprintf(parmfile,"%f;",op->ne_ratio[j]);
   }
   if (i == 13 && b) { /* numpanels */
      fprintf(parmfile,":%ld;",numpop);
      for(j = 0; j < numpop; j++)
         fprintf(parmfile,"%ld;",op->numpanel[j]);
   }
   if (i == 16 && b) { /* full snp */
      fprintf(parmfile,":%f;",op->chance_seen);
   }
   fprintf(parmfile,"\n");
}

for(i = 0; i < NUMNUMBER; i++) {
   fprintf(parmfile,"%s=",numbertokens[i]);
   if (i == 0) fprintf(parmfile,"%f",locus_ttratio);
   if (i == 1) fprintf(parmfile,"%ld",op->numchains[0]);
   if (i == 2) fprintf(parmfile,"%ld",op->steps[0]);
   if (i == 3) fprintf(parmfile,"%ld",op->increm[0]);
   if (i == 4) fprintf(parmfile,"%ld",op->numchains[1]);
   if (i == 5) fprintf(parmfile,"%ld",op->increm[1]);
   if (i == 6) fprintf(parmfile,"%ld",op->steps[1]);
   if (i == 7) fprintf(parmfile,"%f",rec0);
   if (i == 8) {
      if (!op->holding) fprintf(parmfile,"none");
      if (op->holding == 1) fprintf(parmfile,"theta");
      if (op->holding == 2) fprintf(parmfile,"recombination");
   }
   if (i == 9) fprintf(parmfile,"%f",op->mutrait);
   if (i == 10) fprintf(parmfile,"%f",op->traitratio);
   if (i == 11) fprintf(parmfile,"%f",op->pd);
   if (i == 12) fprintf(parmfile,"%ld",op->hapdrop);
   if (i == 13) {
      fprintf(parmfile,"%ld;",op->numtempchains);
      for(j = 0; j < op->numtempchains; j++)
         fprintf(parmfile,"%ld;",op->temperature[j]);
   }
   fprintf(parmfile,"\n");
}

fprintf(parmfile,"end\n");

FClose(parmfile);

#ifdef MAC
   fixmacfile(parmfilename);
#endif

} /* rec_parmfilewrite */


void readseedfile(void)
{
long inseed;

seedfile = fopen("seedfile","r");

if (!MENU) {
   if (!seedfile) {
      fprintf(ERRFILE,"\nseedfile not present!\n");
      exit(-1);
   }
   fscanf(seedfile, "%ld%*[^\n]", &inseed);
   getc(seedfile);
} else {
   if (seedfile) {
      fscanf(seedfile, "%ld%*[^\n]", &inseed);
      getc(seedfile);
   } else {
      printf("Random number seed (must be odd)?\n");
      scanf("%ld%*[^\n]", &inseed);
      getchar();
   }
}

seed[0] = inseed & 2047;
inseed /= 2048;
seed[1] = inseed & 2047;
inseed /= 2048;
seed[2] = inseed & 1023;

} /* readseedfile */


void print_menuheader(option_struct *op)
{
printf("\n%s", op->ansi ? "\033[2J\033[H" : "\n");
printf("Metropolis-Hastings Markov Chain Monte Carlo");
printf(" method, version %3.2f\n\n",VERSION_NUM);
} /* print_menuheader */


void print_startmenu(option_struct *op, boolean writeout)
{
printf("STARTUP MENU\n");
printf("  #               Goto Data/Search Menus\n");
printf("  O         Save current options to file?  %s\n",
   writeout ? "Yes" : "No");
printf("  N          Use trees from previous run?  %s\n",
   op->newdata ? "No" : "Yes");
if (op->newdata)
   printf("  E        Echo the data at start of run?  %s\n",
      op->printdata ? "Yes" : "No");
else op->printdata = FALSE;
printf("  S            Save MHMC output to files?  %s\n",
   op->mhmcsave ? "Yes" : "No");
printf("  P Print indications of progress of run?  %s\n",
   op->progress ? "Yes" : "No");
#if 0 /* DEBUG debug WARNING warning--no tree writer yet! */
printf("  G                Print out genealogies?  %s\n",
   op->treeprint ? "Yes" : "No");
#endif
printf("  U      Use user tree in file \"intree\" ?  %s\n",
   op->usertree ? "Yes" : "No");
printf("  V          Number of temperatures used?  %ld\n",
   op->numtempchains);
printf("  H                     Infer haplotypes?  %s\n",
   op->haplotyping ? "Yes" : "No");
printf("  L       Calculate confidence intervals?  %s\n",
   op->profile ? "Yes" : "No");
if (op->haplotyping) {
   printf("  A Strategy for haplotype rearrangement?  ");
   if(op->hapdrop==0)printf("No resimulation\n");
   else if (op->hapdrop==1)printf("Single resimulation\n");
   else printf("Double resimulation\n");
   printf("  F Frequency of haplotype rearrangement?  %f\n",
      op->happrob);
}
printf("  M     Map trait information onto trees?  %s\n",
   op->map ? "Yes" : "No");
if (op->map) {
   printf("  R        Relative mutation rate of trait?  %f\n",
      op->mutrait);
   printf("  T     Forward versus back trait mutation?  %f\n",
      op->traitratio);
   printf("  D                     Frequency of trait?  %f\n",
      op->pd);
}

} /* print_startmenu */


void print_datamenu(option_struct *op)
{
printf("DATA MENU\n");
printf("  #                     Goto Startup Menu\n");
printf("  D                             Datatype:");
switch(op->datatype) {
   case 'a':
      printf("  Allelelic Markers\n");
      break;
   case 'b':
      printf("  Brownian-Motion Microsats\n");
   case 'm':
      if (op->datatype != 'b') printf("  Microsats\n");
      break;
   case 'n':
      printf("  SNPs\n");
      printf("  P              SNPs derived from panel?  %s\n",
         op->panel ? "Yes" : "No");
#if 0
      printf("  A   SNPs are all of the variable sites?  %s\n",
         op->full ? "No" : "Yes");
      if (op->full) {
         printf("  B           Probability of observation?  %6.4f\n",
            op->chance_seen);
      }
#endif
      printf("  N              SNPs with recombination?  %s\n",
         op->norecsnp ? "No" : "Yes");
   case 's':
      if (op->datatype != 'n') printf("  Sequence\n");
      printf("  I          Input sequences interleaved?  %s\n",
         op->interleaved ? "Yes" : "No, sequential");
      printf("  T        Transition/transversion ratio:  %8.4f\n",
         locus_ttratio);
      printf("  F       Use empirical base frequencies?  %s\n",
	 op->freqsfrom ? "Yes" : "No");
      printf("  C   One category of substitution rates?");
      if (!op->ctgry || op->categs == 1)
         printf("  Yes\n");
      else {
         printf("  %ld categories\n", op->categs);
         if (op->datatype != 'n') {
            printf("  R   Rates at adjacent sites correlated?");
            if (!op->autocorr) printf("  No, they are independent\n");
            else printf("  Yes, mean block length =%6.1f\n", 1.0 / op->lambda);
         }
      }
      printf("  V    Poplulation size equal among loci?  %s\n",
         op->same_ne ? "Yes" : "No");
      break;
   default:
      printf("  UNKNOWN\n");
      break;
}

} /* print_datamenu */


void print_searchmenu(option_struct *op)
{
printf("SEARCH MENU\n");
printf("  Q      Lots of recombinations expected?  %s\n",
   op->fc ? "Yes" : "No");
printf("  H                Hold parameters fixed?");
   if (op->holding == 0) printf("  No\n");
   else if (op->holding == 1) printf("  Yes, theta fixed\n");
        else if (op->holding == 2) printf("  Yes, rec-rate fixed\n");
             else printf("  Unknown option!!!\n");
printf("  W      Starting theta equals Watterson?  %s",
   op->watt ? "Yes\n" : "No");
if (!op->watt) printf(", initial theta = %6.4f\n", theta0);
printf("  Z          Starting recombination rate:  %e\n", rec0);
printf("  S               Number of short chains?  %6ld\n",
   op->numchains[0]);
if (op->numchains[0] > 0) {
   printf("  1             Short sampling increment?  %6ld\n",
      op->increm[0]);
   printf("  2   Number of steps along short chains?  %6ld\n",
      op->steps[0]);
}
printf("  L                Number of long chains?  %6ld\n",
   op->numchains[1]);
if (op->numchains[1] > 0) {
   printf("  3              Long sampling increment?  %6ld\n",
      op->increm[1]);
   printf("  4    Number of steps along long chains?  %6ld\n",
      op->steps[1]);
}

} /* print_searchmenu */


void print_menuend(void)
{
printf("\nAre these settings correct? (type Y or the letter for");
printf(" one to change)\n");
} /* print_menuend */


void initoptions(option_struct *op)
{
long i;

/* first some multiple rate-categories code stuff */
op->ctgry = FALSE;
op->rate[0] = 1.0;
op->probcat[0] = 1.0;
op->categs = 1;
op->lambda = 1.0;
op->autocorr = FALSE;  /* FALSE if categs == 1 */
/* end categories code stuff */

op->interactive = TRUE;
op->newdata = TRUE;
op->mhmcsave = FALSE;
op->interleaved = TRUE;
op->printdata = FALSE;
op->progress = TRUE;
op->treeprint = FALSE;
op->profile = TRUE;

weightfile = fopen("weightfile","r");
op->weights = (weightfile) ? TRUE : FALSE;

spacefile = fopen("spacefile","r");
op->spacing = (spacefile) ? TRUE : FALSE;

op->map = FALSE;
op->mutrait = 1.0;
op->traitratio = 1.0;
op->pd = 0.1;

op->panel = FALSE;
op->numpanel = NULL;

op->haplotyping = FALSE;
op->hapdrop = 0;
op->happrob = 0.2;

op->full = FALSE;
op->chance_seen = 1.0;
op->norecsnp = FALSE;

op->numtempchains = 1;
op->temperature = (long *)calloc(op->numtempchains,sizeof(long));
for(i = 0; i < op->numtempchains; i++) op->temperature[i] = 1;
op->ctemp = 1;

op->fc = FALSE;
op->freqsfrom = TRUE;
op->watt = FALSE;
op->usertree = FALSE;
op->plump = TRUE;
op->holding = 0;
op->same_ne = TRUE;
op->ne_ratio = NULL;
op->numchains[0] = 5;
op->increm[0] = 20;
op->steps[0] = 2000;
op->numchains[1] = 1;
op->increm[1] = 20;
op->steps[1] = 50000;
op->datatype = 's';
op->datatypeset = FALSE;

op->print_recbythmaxcurve = TRUE;
op->thlb = FLAGLONG;
op->thub = FLAGLONG;
op->reclb = FLAGLONG;
op->recub = FLAGLONG;

} /* initoptions */


void workmenu1(option_struct *op, char ch, boolean *menu1, boolean *writeout)
{
long i;
char input[LINESIZE];

if (strchr("#FLAHONESPGUMVRTD",ch) != NULL) {
   switch(ch) {
      case '#':
         *menu1 = !(*menu1);
         break;
      case 'V':
         printf("Number of different temperatures?\n");
         scanf("%ld",&(op->numtempchains));
         scanf("%*[^\n]");
         do {
            printf("Temperature of coldest chain? (must be positive)\n");
            scanf("%ld",&(op->temperature[0]));
            scanf("%*[^\n]");
         } while (op->temperature[0] <= 0.0);
         for (i = 1; i < op->numtempchains; i++) {
            do {
               printf("Temperature of temperature-chain %ld?\n",i+1);
               scanf("%ld",&(op->temperature[i]));
               scanf("%*[^\n]");
            } while(op->temperature[i] <= op->temperature[0]);
         }
         break;
      case 'F':
         printf("Probability of proposing a change in haplotype?\n");
         scanf("%lf",&(op->happrob));
         scanf("%*[^\n]");
         break;
      case 'L':
         op->profile = !op->profile;
         break;
      case 'H':
         op->haplotyping = !op->haplotyping;
         break;
      case 'O':
         *writeout = !(*writeout);
         break;
      case 'N':
         op->newdata = !op->newdata;
         break;
      case 'E':
         op->printdata = !op->printdata;
	 break;
      case 'P':
         op->progress = !op->progress;
         break;
      case 'G':
         op->treeprint = !op->treeprint;
         break;
      case 'U':
         op->usertree = !op->usertree;
         break;
      case 'S':
         op->mhmcsave = !op->mhmcsave;
         break;
      case 'M':
         op->map = !op->map;
         break;
      case 'R':
         do {
            printf("Relative mutation rate of trait?\n");
            fgets(input,LINESIZE,stdin);
            op->mutrait = atof(input);
         } while (op->mutrait <= 0.0);
         break;
      case 'T':
         do {
            printf("Ratio of forward to back trait mutation?\n");
            fgets(input,LINESIZE,stdin);
            op->traitratio = atof(input);
         } while (op->traitratio <= 0.0);
         break;
      case 'D':
         do {
            printf("Frequency of trait?\n");
            fgets(input,LINESIZE,stdin);
            op->pd = atof(input);
         } while (op->pd <= 0.0 || op->pd >= 1.0);
         break;
      case 'A':
         do {
            printf("Number of drops while resimulating (0-2)?\n");
            fgets(input,LINESIZE,stdin);
            op->hapdrop = atol(input);
         } while (op->hapdrop != 0 && op->hapdrop != 1 && op->hapdrop != 2);
      default:
         fprintf(stderr,"ERROR:Impossible option %c detected!\n",ch);
         break;
   }
} else fprintf(stderr,"%c is not an option here.\n",ch);

} /* workmenu1 */


void workmenu2(option_struct *op, char ch, boolean *menu1, long *numloci,
   long *numpop)
{
long i;
double probsum;
char input[LINESIZE];
boolean done;

if(strchr("#NQPDITFCRVHWZS12L34AB",ch) != NULL) {
   switch(ch) {
      case '#':
         *menu1 = !(*menu1);
         break;
      case 'N':
         op->norecsnp = !(op->norecsnp);
         if (op->norecsnp) {
            if (rec0 || op->holding != 2) {
               printf("\nThis SNP model requires recombination");
               printf(" rate to be fixed at zero.\n");
               printf("Fix starting rec-rate to zero first.\n");
               printf("----SNP model unchanged.----\n");
               op->norecsnp = FALSE;
            }
         }
         printf("press <enter> or <return> to continue\n");
         fgets(input,LINESIZE,stdin);
         break;
      case 'Q':
         op->fc = !op->fc;
         break;
      case 'D':
/* WARNING--most datatypes disabled!!! */
#if 0
         if (op->datatype == 'a') {op->datatype = 'b'; break;}
         if (op->datatype == 'b') {op->datatype = 'm'; break;}
         if (op->datatype == 'm') {op->datatype = 'n'; break;}
         if (op->datatype == 'n') {op->datatype = 's'; break;}
         if (op->datatype == 's') {op->datatype = 'a'; break;}
#endif
         op->datatypeset = TRUE;
         if (op->datatype == 'n') {op->datatype = 's'; break;}
         if (op->datatype == 's') {
            op->datatype = 'n';
            if (op->autocorr) {
               printf("Can't use autocorrelation with SNPs\n");
               printf("----autocorrelation disabled-----\n");
               printf("press <enter> or <return> to continue\n");
               fgets(input,LINESIZE,stdin);
               op->autocorr = FALSE;
            }
            break;
         }
         printf("ERROR:Impossible Datatype %c present\n",op->datatype);
         break;
#if 0
      case 'A': /* deliberate fall through to case 'B' */
         op->full = !op->full;
      case 'B':
         if (!op->full) break;
         do {
            printf("Percent chance that a variable site was observed");
            printf(" during sequencing?\n");
            scanf("%lf",&(op->chance_seen));
            scanf("%*[^\n]");
            if (op->chance_seen <= 0.0 || op->chance_seen > 1.0)
               printf("   the percentage must be between 0 and 1\n");
         } while (op->chance_seen <= 0.0 || op->chance_seen > 1.0);
         break;
#endif
      case 'I':
         op->interleaved = !op->interleaved;
         break;
      case 'T':
         do {
            printf("Transition/transversion ratio?\n");
            fgets(input,LINESIZE,stdin);
            locus_ttratio = atof(input);
            if (locus_ttratio < 0.5)
               printf("TTratio cannot be less than 0.5\n");
         } while (locus_ttratio < 0.5);
         break;
      case 'F':
         op->freqsfrom = !op->freqsfrom;
         if (!op->freqsfrom) {
            printf("Base frequencies for A, C, G, T/U (use blanks");
            printf(" to separate)?\n");
            scanf("%lf%lf%lf%lf", &freqa, &freqc, &freqg, &freqt);
            scanf("%*[^\n]");
	 }
         break;
      case 'P':
         op->panel = !op->panel;
         if (op->panel) {
            printf("Number of populations?\n");
            ;
            *numpop = atol(input);
            op->numpanel = (long *)calloc(*numpop,sizeof(long));
            for(i = 0; i < *numpop; i++) {
               printf("Number of panel haplotypes for population");
               printf(" %ld?\n",i+1);
               fgets(input,LINESIZE,stdin);
               op->numpanel[i] = atol(input);
            }
         } else
            if (op->numpanel) {
               free(op->numpanel);
               op->numpanel = NULL;
            }
         break;
      case 'C':
         op->ctgry = !op->ctgry;
         if (!op->ctgry) op->autocorr = FALSE;
         if (op->ctgry) {
            do {
               printf("Number of categories ?");
               fgets(input,LINESIZE,stdin);
               op->categs = atoi(input);
            } while (op->categs < 1);
            free(op->rate);
            free(op->probcat);
            printf("Rate for each category? (use a space to");
	    printf(" separate)\n");
            op->rate = (double *)calloc(op->categs,sizeof(double));
            op->probcat = (double *)calloc(op->categs,sizeof(double));
            for (i = 0; i < op->categs; i++)
               scanf("%lf*[^\n]", &(op->rate[i]));
            getchar();
            do {
               printf("Probability for each category?");
               printf(" (use a space to separate)\n");
               for (i = 0; i < op->categs; i++)
                  scanf("%lf", &(op->probcat[i]));
               scanf("%*[^\n]");
               getchar();
               done = TRUE;
               probsum = 0.0;
               for (i = 0; i < op->categs; i++)
                  probsum += op->probcat[i];
               if (fabs(1.0 - probsum) > 0.001) {
                  done = FALSE;
                  printf("Probabilities must add up to");
                  printf(" 1.0, plus or minus 0.001.\n");
               }
            } while (!done);
	 }
         break;
      case 'R':
         if (op->datatype == 'n') {
            printf("Can't use autocorrelation with SNPs\n");
            printf("----autocorrelation disabled-----\n");
            printf("press <enter> or <return> to continue\n");
            fgets(input,LINESIZE,stdin);
            op->autocorr = FALSE;
            break;
         }
         op->autocorr = !op->autocorr;
         if (op->autocorr) {
            do {
               printf("Mean block length of sites having the same ");
               printf("rate (greater than 1)?\n");
               scanf("%lf%*[^\n]", &(op->lambda));
               getchar();
            } while (op->lambda <= 1.0);
            op->lambda = 1.0 / op->lambda;
         }
         break;
      case 'V':
         op->same_ne = !op->same_ne;
         if (!op->same_ne) {
            printf("Enter number of loci: ");
            scanf("%ld",numloci);
            if (*numloci == 1) op->same_ne = !op->same_ne;
            if (op->ne_ratio) free(op->ne_ratio);
            op->ne_ratio = (double *)calloc(*numloci,sizeof(double));
            printf("\nEnter relative population size/mutation rate");
            printf(" for each locus in input order:\n");
            for(i = 0; i < *numloci; i++) {
               scanf("%lf",&(op->ne_ratio[i]));
               if (op->ne_ratio[i] <= 0) {
                  printf("\nratios must be positive, please reenter\n");
                  i--;
               }
            }
         }
         break;
      case 'H':
         op->holding += 1;
         if (op->holding > 2) op->holding = 0;
         break;
      case 'W':
         op->watt = !op->watt;
         if (!op->watt) {
            do {
               printf("Initial theta estimate?\n");
               fgets(input,LINESIZE,stdin);
               theta0 = atof(input);
            } while (theta0 <= 0.0);
         }
         break;
      case 'Z':
         printf("What recombination rate?\n");
         do {
            fgets(input,LINESIZE,stdin);
            rec0 = atof(input);
            if (rec0 < 0.0)
               printf("recombination rate must be non-negative\n");
         } while (rec0 < 0.0);
         break;
      case 'S':
         do {
            printf("How many Short Chains?\n");
            fgets(input,LINESIZE,stdin);
            op->numchains[0] = atoi(input);
            if (op->numchains[0] < 0)
            printf("Must be non-negative\n");
         } while (op->numchains[0] < 0);
         break;
      case '1':
         done = FALSE;
         while (!done) {
            printf("How often to sample trees?\n");
            fgets(input,LINESIZE,stdin);
            op->increm[0] = atoi(input);
            if (op->increm[0] > 0) done = TRUE;
            else printf("Must be a positive integer\n");
         }
         break;
      case '2':
         done = FALSE;
         while (!done) {
            printf("How many short steps?\n");
            fgets(input,LINESIZE,stdin);
            op->steps[0] = atoi(input);
            if (op->steps[0] > 0) done = TRUE;
            else printf("Must be a positive integer\n");
         }
         break;
      case 'L':
         do {
            printf("How many Long Chains?\n");
            fgets(input,LINESIZE,stdin);
            op->numchains[1] = atoi(input);
            if (op->numchains[1] < 1)
            printf("Must be a positive integer\n");
         } while (op->numchains[1] < 1);
         break;
      case '3':
         done = FALSE;
         while (!done) {
            printf("How often to sample trees?\n");
            fgets(input,LINESIZE,stdin);
            op->increm[1] = atoi(input);
            if (op->increm[1] > 0) done = TRUE;
            else printf("Must be a positive integer\n");
         }
         break;
      case '4':
         done = FALSE;
         while (!done) {
            printf("How many long steps?\n");
            fgets(input,LINESIZE,stdin);
            op->steps[1] = atoi(input);
            if (op->steps[1] > 0) done = TRUE;
            else printf("Must be a positive integer\n");
         }
         break;
      default:
         fprintf(stderr,"ERROR:Impossible option %c detected!\n",ch);
         break;
   }
} else fprintf(stderr,"%c is not an option here.\n",ch);

} /* workmenu2 */


/****************************************************************
 * checkparmfile() checks the parmfile for internal consistency */
void checkparmfile(option_struct *op)
{
boolean error = FALSE;

if (op->norecsnp) {
   if (rec0) {
      fprintf(ERRFILE,"ERROR--this snp model is incompatible with a");
      fprintf(ERRFILE," non-zero recombination rate.\n");
      fprintf(ERRFILE,"   Please change starting rec-rate.\n");
      fprintf(outfile,"ERROR--this snp model is incompatible with a");
      fprintf(outfile," non-zero recombination rate.\n");
      fprintf(outfile,"   Please change starting rec-rate.\n");
      error = TRUE;
   }
   if (op->holding != 2) {
      fprintf(ERRFILE,"ERROR--this snp model is incompatible with a");
      fprintf(ERRFILE," non-zero recombination rate.\n");
      fprintf(ERRFILE,"   Please set rec-rate constant.\n");
      fprintf(outfile,"ERROR--this snp model is incompatible with a");
      fprintf(outfile," non-zero recombination rate.\n");
      fprintf(outfile,"   Please set rec-rate constant.\n");
      error = TRUE;
   }
}

if (error) {
   fprintf(ERRFILE,"program finished\n");
   exit(-1);
}

} /* checkparmfile */


void getoptions(option_struct *op)
/* interactively set options using a very basic menu */
{
long numloci, numpop;
boolean writeout, done, menu1;
char ch, input[LINESIZE];

/* default initializations */
initoptions(op);
writeout = FALSE;
locus_ttratio = 2.0;
numloci = 1;
numpop = 1;
theta0 = 1.0;
rec0 = 0.0001;
/* end defaults */

readparmfile(op);
checkparmfile(op);
readseedfile();

fprintf(outfile, "\nRECOMBINE:  Recombinant HMMC ML method,");
fprintf(outfile, " version %3.2f\n\n",VERSION_NUM);

if (!MENU) return;

menu1 = TRUE;
do {
   print_menuheader(op);
   if (menu1) print_startmenu(op,writeout);
   else {print_datamenu(op); print_searchmenu(op);}
   print_menuend();
   fgets(input,LINESIZE,stdin);
   ch = toupper((int)input[0]);
   done = (ch == 'Y');
   if (!done) {
      if (menu1) workmenu1(op,ch,&menu1,&writeout);
      else workmenu2(op,ch,&menu1,&numloci,&numpop);
   }
} while (!done);

if (writeout) rec_parmfilewrite(op,numloci,numpop);

}  /* getoptions */


/*********************************************************************
 * firstinit() handles initialization for things that are recorded   *
 * over multiple loci/populations, and therefore are allocated once. */
void firstinit(option_struct *op, data_fmt *data)
{
long i, numloci;

numloci = getdata_numloci(op,data);

totchains = op->numchains[0] + op->numchains[1];

numtrees = MAX(op->steps[0]/op->increm[0],op->steps[1]/op->increm[1]);

if (op->same_ne) {
   op->ne_ratio = (double *)calloc(1,numloci*sizeof(double));
   for(i = 0; i < numloci; i++) op->ne_ratio[i] = 1.0;
}

sametree = (boolean **)calloc(1,numloci * sizeof(boolean *));
sametree[0] = (boolean *)calloc(1,numloci*numtrees * sizeof(boolean));
for(i = 1; i < numloci; i++)
   sametree[i] = sametree[0] + i*numtrees;

model_alloc(op,data);

freenodes = NULL;

}  /* firstinit */


/********************************************************************
 * popinit() initializes things that are specific to one population */
void popinit(option_struct *op, data_fmt *data)
{

switch(op->datatype) {
   case 'a':
      break;
   case 'b':
   case 'm':
      break;
   case 'n':
   case 's':
      rootnum = 0;
      break;
   default:
      fprintf(ERRFILE,"ERROR:popinit, can't get here!\n");
      exit(-1);
}

} /* popinit */


/***********************************************************************
 * makeinvarvalues() sets up the tip likelikelihoods for the invariant *
 * sites tacked on the end of the sequence for SNP data.               */
void makeinvarvalues(dnadata *dna, long numcategs, tree *tr)
{
long site, ind, categ, slice;

site = dna->sites[locus];
slice = 0L; /* we only set slice 0 here; the remainder, if any,
  are set elsewhere */

for(ind = 1; ind <= dna->numseq[population]; ind++)
   for(categ = 0; categ < numcategs; categ++) {
      tr->nodep[ind]->x->s[site][categ][slice][baseA] = 1.0;
      tr->nodep[ind]->x->s[site][categ][slice][baseC] = 0.0;
      tr->nodep[ind]->x->s[site][categ][slice][baseG] = 0.0;
      tr->nodep[ind]->x->s[site][categ][slice][baseT] = 0.0;

      tr->nodep[ind]->x->s[site+1][categ][slice][baseA] = 0.0;
      tr->nodep[ind]->x->s[site+1][categ][slice][baseC] = 1.0;
      tr->nodep[ind]->x->s[site+1][categ][slice][baseG] = 0.0;
      tr->nodep[ind]->x->s[site+1][categ][slice][baseT] = 0.0;

      tr->nodep[ind]->x->s[site+2][categ][slice][baseA] = 0.0;
      tr->nodep[ind]->x->s[site+2][categ][slice][baseC] = 0.0;
      tr->nodep[ind]->x->s[site+2][categ][slice][baseG] = 1.0;
      tr->nodep[ind]->x->s[site+2][categ][slice][baseT] = 0.0;

      tr->nodep[ind]->x->s[site+3][categ][slice][baseA] = 0.0;
      tr->nodep[ind]->x->s[site+3][categ][slice][baseC] = 0.0;
      tr->nodep[ind]->x->s[site+3][categ][slice][baseG] = 0.0;
      tr->nodep[ind]->x->s[site+3][categ][slice][baseT] = 1.0;
   }
} /* makeinvarvalues */


/*****************************************************************
 * locusinit() initializes things that are specific to one locus */
void locusinit(option_struct *op, data_fmt *data)
{
long i, nummarkers, numseq, numsites;
dnadata *dna;
FILE *flipfile;

dna = data->dnaptr;

if (op->spacing) read_spacefile(spacefile,op,data);
nummarkers = getdata_nummarkers(op,data);
/* hack done for speed-up */
dna->sitecount[locus][0] = getdata_space(op,data,0L);
dna->markersite[locus][0] = dna->sspace[0][0][getdata_nummarkers(op,data)];
for(i = 1; i < nummarkers; i++) {
   dna->sitecount[locus][i] = dna->sitecount[locus][i-1] +
                              getdata_space(op,data,i);
   dna->markersite[locus][i] = dna->markersite[locus][i-1] +
                               dna->sspace[0][0][i-1];
}
numsites = countsites(op,data);
dna->segranges0 = (int *)calloc(numsites,sizeof(int));
dna->segranges1 = (int *)calloc(numsites,sizeof(int));
dna->segranges2 = (int *)calloc(numsites,sizeof(int));
dna->segranges3 = (int *)calloc(numsites,sizeof(int));
for(i = 0; i < numsites; i++) {
   dna->segranges0[i] = 0;
   dna->segranges1[i] = 1;
   dna->segranges2[i] = 2;
   dna->segranges3[i] = -1;
}

treesetup(op,data);

numseq = getdata_numseq(op,data);

switch(op->datatype) {
   case 'a':
      break;
   case 'b':
   case 'm':
      break;
   case 'n':
      initinvartips(op,data,op->categs,curtree);
      if (!op->panel) invardatacheck(op,data);
   case 's':
      fprintf(outfile, "%4ld Sequences, %4ld Sites\n",numseq,nummarkers);
      alogf = (rlrec *)calloc(1,sizeof(rlrec));
      alogf->val = (double *)calloc(nummarkers+1,sizeof(double));
      contribution =
         (contribarr *)calloc(nummarkers+1,sizeof(contribarr));
      contribution[0] =
         (double *)calloc((nummarkers+1)*op->categs,sizeof(double));
      for(i=1;i<nummarkers+1;i++)
         contribution[i] = contribution[0] + i*op->categs;
      /* now read the categories from the infile, if applicable */
      inputcategories(op,data);
      /* setup fractional likelihoods at tree tips */
      makednavalues(op,data,op->categs,curtree);
      /* frequencies must be calculated after makednavalues() called */
      if (!op->freqsfrom) {
         data->dnaptr->freqa = freqa;
         data->dnaptr->freqc = freqc;
         data->dnaptr->freqg = freqg;
         data->dnaptr->freqt = freqt;
      } else empiricaldnafreqs(op,data,curtree);
      getbasednafreqs(data->dnaptr,op,locus_ttratio,outfile);
      /* table can only be initialized after frequencies are set */
      inittable(op,dna);
      makesiteptr(op,data);
      initweightrat(op,data);
      if (op->haplotyping) {
         flipfile = fopen("flipfile","r");
         read_flipfile(op,data,curtree,flipfile);
         fclose(flipfile);
         pruneflips(op,data,curtree);
      }
      break;
   default:
      fprintf(ERRFILE,"ERROR:locusinit, can't get here!\n");
      exit(-1);
}


if ((op->increm[0] < 0) || (op->increm[1] < 0)) {
   fprintf(ERRFILE,"Error in input sampling increment,");
   fprintf(ERRFILE," increment set to 10\n");
   if (op->increm[0] < 0)
      op->increm[0] = 10;
   if (op->increm[1] < 0)
      op->increm[1] = 10;
}
/* stuff for recycling x arrays */
numx = 0;
for (i=0;i<XARRAYSIZE;i++)
  sparex[i]=NULL;

} /* locusinit */


void inputcategories(option_struct *op, data_fmt *data)
{
/*  char ch; */
  long i, extranum, nummarkers;

  nummarkers = getdata_nummarkers(op,data);

  category = (long *)calloc(nummarkers,sizeof(long));

  for (i = 0; i < nummarkers; i++) category[i] = 1;
  extranum = 0;
/* DEBUG debug WARNING warning--categories input code commented out!
  while (!(eoln(infile))) {
    ch = getc(infile);
    if (ch == '\n')
      ch = ' ';
    ch = isupper(ch) ? ch : toupper((int)ch);
    if (ch == 'C')
      extranum++;
    else if (ch != ' ') {
      printf("BAD OPTION CHARACTER: %c\n", ch);
      exit(-1);
    }
  }
  fscanf(infile, "%*[^\n]");
  getc(infile);
  for (i = 1; i <= extranum; i++) {
    ch = getc(infile);
    if (ch == '\n')
      ch = ' ';
    ch = isupper(ch) ? ch : toupper((int)ch);
    if (ch != 'W'){
      printf("ERROR: INCORRECT AUXILIARY OPTIONS LINE WHICH STARTS WITH %c\n",
	     ch);
      exit(-1);
    }
  }
*/
  if (op->categs <= 1)
    return;
  fprintf(outfile, "\nSite category   Rate of change  Probability\n");
  for (i = 1; i <= op->categs; i++)
    fprintf(outfile, "%12ld%13.3f%13.3f\n",i,op->rate[i-1],op->probcat[i-1]);
  putc('\n', outfile);
}  /* inputcategories */


/*******************************************************************
 * allocate_nodelet() allocates and initializes a set of nodelets, *
 * which are set to be non-tops of the passed "type".  "x" and all *
 * other datatype specific fields are NULLed and must be allocated *
 * elsewhere.                                                      */
node *allocate_nodelet(long num, char type)
{
static long unique_id=0;
boolean isfirst=TRUE;
long j;
node *p, *q = NULL, *pfirst=NULL;

for (j = 0; j < num; j++) {
   p = (node *)calloc(1,sizeof(node));

   p->top = FALSE;
   p->updated = FALSE;
   p->number = FLAGLONG;
   p->type = type;
   p->id = unique_id++;
   p->next = q;
   p->back = NULL;
   p->nayme = NULL;
   p->x = NULL;
   p->v = p->tyme = p->length = 0.0;
   p->z = NULL;

   p->recstart = p->recend = -1L;
   p->ranges = NULL;
   p->coal = NULL;
   p->members = NULL;
   p->memberstrue = FALSE;
   p->futileflag = 0L;

   p->pop = p->actualpop = -1;
   p->lxmax = 0.0;

   if(isfirst){
      isfirst=FALSE;
      pfirst = p;
   }

   q = p;
}

pfirst->next = q;

return q;

} /* allocate_nodelet */


/*************************************************************
 * allocate_root() allocates and initializes a root for "tr" */
void allocate_root(option_struct *op, data_fmt *data, tree *tr)
{
long i, j, k, numslice, nummarkers;
dnadata *dna;
msatdata *ms;

tr->nodep[ROOTNUM] = allocate_nodelet(1,'t');
tr->nodep[ROOTNUM]->number = ROOTNUM;
allocate_x(op,data,tr->nodep[ROOTNUM]);
if (op->map) {
   tr->nodep[ROOTNUM]->z = NULL;
   allocate_z(op,data,tr->nodep[ROOTNUM]);
}
ranges_Malloc(tr->nodep[ROOTNUM],FALSE,0L);
coal_Malloc(tr->nodep[ROOTNUM],FALSE,0L);
tr->nodep[ROOTNUM]->nayme = (char *)calloc(NMLNGTH,sizeof(char));
strncpy(tr->nodep[ROOTNUM]->nayme,"ROOT",4);

/* guarantee that the root node contributes nothing to the likelihood
   of a tree (since its supposed to be at the end of a theoretically
   infinite root branch) */
switch(op->datatype) {
   case 'a':
      break;
   case 'b':
   case 'm':
      ms = data->msptr;
      for(i = 0; i < ms->numloci; i++)
         for(j = 0; j < MICRO_ALLELEMAX; j++)
            tr->nodep[ROOTNUM]->x->a[i][j] = 1.0;
      break;
   case 'n':
   case 's':
      dna = data->dnaptr;
      numslice = (op->panel) ? NUMSLICE : 1L;
      nummarkers = getdata_nummarkers(op,data);
      for (i = 0; i < nummarkers; i++) {
        for (j = 0; j < op->categs; j++) {
          for (k = 0; k < numslice; k++) {
            tr->nodep[ROOTNUM]->x->s[i][j][k][baseA] = 1.0;
            tr->nodep[ROOTNUM]->x->s[i][j][k][baseC] = 1.0;
            tr->nodep[ROOTNUM]->x->s[i][j][k][baseG] = 1.0;
            tr->nodep[ROOTNUM]->x->s[i][j][k][baseT] = 1.0;
          }
        }
      }
      break;
   default:
      fprintf(ERRFILE,"ERROR:allocate_root: can't get here!\n");
      exit(-1);
}

} /* allocate_root */


/***********************************************************
 * allocate_tip() allocates and initializes a tip nodelet. */
void allocate_tip(option_struct *op, data_fmt *data, tree *tr,
   long num)
{

tr->nodep[num] = allocate_nodelet(1,'t');
tr->nodep[num]->number = num;
tr->nodep[num]->top = TRUE;
tr->nodep[num]->tyme = 0.0;
tr->nodep[num]->nayme = (char *)calloc((NMLNGTH+1),sizeof(char));
allocate_x(op,data,tr->nodep[num]);
if (op->map) {
   tr->nodep[num]->z = NULL;
   allocate_z(op,data,tr->nodep[num]);
}
ranges_Malloc(tr->nodep[num],TRUE,1L);
if (op->fc) coal_Malloc(tr->nodep[num],TRUE,1L);

switch(op->datatype) {
   case 'a':
      break;
   case 'b':
   case 'm':
      break;
   case 'n':
   case 's':
      tr->nodep[num]->ranges[1] = 0L;
      tr->nodep[num]->ranges[2] = countsites(op,data)-1;
      if (op->fc) {
         tr->nodep[num]->coal[1] = 0L;
         tr->nodep[num]->coal[2] = countsites(op,data)-1;
      }
      break;
   default:
      fprintf(ERRFILE,"ERROR:allocate_tip: can't get here!\n");
      exit(-1);
}

} /* allocate_tip */


/*******************************************************************
 * allocate_interior() allocates and initializes an interior node. */
void allocate_interior(option_struct *op, data_fmt *data, tree *tr, long num)
{

if (freenodes == NULL) tr->nodep[num] = allocate_nodelet(3,'c');
else {tr->nodep[num] = freenodes; freenodes = freenodes->back;}
tr->nodep[num]->number = tr->nodep[num]->next->number =
   tr->nodep[num]->next->next->number = num;

} /* allocate_interior */


/***********************************************************************
 * init_creature() initializes an already allocated creature structure */
void init_creature(creature *cr, long numhaplotypes)
{

cr->numhaplotypes = numhaplotypes;
cr->numflipsites = 0;
cr->haplotypes = (node **)calloc(numhaplotypes,sizeof(node *));
cr->flipsites = NULL;

} /* init_creature */


/***************************************************************
 * treesetup() allocates and initializes the global "curtree". *
 * Curtree's tymelist will be handled later.                   */
void treesetup(option_struct *op, data_fmt *data)
{
long i, j, whichtip, numtips, numnodes, numcreatures;
node *p;

numtips = getdata_numtips(op,data);
numnodes = 2 * numtips;
numcreatures = numtips/NUMHAPLOTYPES;
curtree = (tree *)calloc(1,sizeof(tree));
curtree->nodep = (node **)calloc(numnodes,sizeof(node *));

allocate_root(op,data,curtree);

for(i = 1; i <= numtips; i++) allocate_tip(op,data,curtree,i);
for(i = numtips+1; i < numnodes; i++) allocate_interior(op,data,curtree,i);
for(p = freenodes; p != NULL; p = p->back, i++)
   p->number = p->next->number = p->next->next->number = i;
nodectr = i;
curtree->likelihood = NEGMAX;
curtree->coalprob = NEGMAX;
curtree->numcoals = numtips-1;
curtree->numrecombs = 0;

if (op->haplotyping) {
   curtree->creatures =
      (creature *)calloc(numcreatures,sizeof(creature));
   for(i = 0, whichtip = 1; i < numcreatures; i++) {
      init_creature(&curtree->creatures[i],NUMHAPLOTYPES);
      for(j = 0; j < curtree->creatures[i].numhaplotypes; j++)
         curtree->creatures[i].haplotypes[j] = curtree->nodep[whichtip++];
   }
   if (whichtip-1 > numtips)
      fprintf(ERRFILE,"\ntreesetup--whichtip too big!\n");
} else curtree->creatures = NULL;

newtymenode(&curtree->tymelist);
curtree->tymelist->branchlist = (node **)calloc(numtips,sizeof(node *));
curtree->tymelist->eventnode = curtree->nodep[1];

} /* treesetup */


/***************************************************************
 * end_of_population_free frees the working tree, curtree; and *
 * the freenode list plus the sparex array.                    */
void end_of_population_free(option_struct *op, data_fmt *data)
{
dnadata *dna;

dna = data->dnaptr;

free(op->ne_ratio);

switch(op->datatype) {
    case 'a':
       break;
    case 'b':
    case 'm':
       break;
    case 'n':
    case 's':
       break;
    default:
       fprintf(ERRFILE,"ERROR:end_of_population_free, can't be here!\n");
       exit(-1);
}

} /* end_of_population_free */


/***************************************************
 * end_of_locus_free() frees locus specific stuff. */
void end_of_locus_free(option_struct *op, data_fmt *data)
{
long i;
node *p, *q;

free(data->siteptr);

switch(op->datatype) {
    case 'a':
       break;
    case 'b':
    case 'm':
       break;
    case 'n':
    case 's':
       free(category);
       free(alogf->val);
       free(alogf);
       /* free the working arrays */
       free(weightrat);
       free(tbl);
       free(contribution[0]);
       free(contribution);
#if 0 /* pre-heating */
       freetree(op,data,curtree);
#endif
       for(i = 0; i < op->numtempchains; i++)
          freetree(op,data,temptrees[i]);
       for(p = freenodes; p != NULL; p = q) {
          q = p->back;
          freenodelet(op,data,p->next->next);
          freenodelet(op,data,p->next);
          freenodelet(op,data,p);
       }
       freenodes = NULL;
       /* the sparex arrays must be done after the tree and freenodes
          have been freed (or anything else that uses free_x to
          free it's xarrays). */
       for(i = 0; i < numx; i++) {
          free(sparex[i]->s[0][0][0]);
          free(sparex[i]->s[0][0]);
          free(sparex[i]->s[0]);
          free(sparex[i]->s);
          free(sparex[i]);
          sparex[i] = NULL;
       }
       numx = 0;
       break;
   default:
       fprintf(ERRFILE,"ERROR:end_of_locus_free, can't be here!\n");
       exit(-1);
}

} /* end_of_locus_free */


boolean sitecompare(dnadata *dna, long site1, long site2)
{
long i;

for(i = 0; i < dna->numseq[population]; i++)
   if(dna->seqs[population][locus][i][site1] !=
      dna->seqs[population][locus][i][site2])
      return FALSE;

return TRUE;
} /* sitecompare */


/*****************************************************************
 * makesiteptr creates the siteptr array:                        *
 * if (!siteptr[site]) there is no alias, otherwise potentially  *
 * use siteptr[site]-1 for the alias site;                       *
 * if (siteptr[site]-1 < 0) then the alias is currently unusable *
 * due to recombination.                                         */
void makesiteptr(option_struct *op, data_fmt *data)
{
long whichsite, lastsite, i, *ptrarray;
boolean found;
dnadata *dna;

dna = data->dnaptr;

if (op->datatype == 'n') lastsite = dna->sites[locus]+NUMINVAR;
else lastsite = dna->sites[locus];

ptrarray = (long *)calloc(lastsite,sizeof(long));
ptrarray[0] = 0;

#if ALIASING
for(whichsite = 1; whichsite < lastsite; whichsite++) {
   if (op->datatype == 'n')
      if (whichsite >= dna->sites[locus]) {
         ptrarray[whichsite] = 0;
         continue;
      }
   found = FALSE;
   for(i = whichsite - 1; i >= 0; i--) {
      if(sitecompare(dna,i,whichsite)) {
         ptrarray[whichsite] = i+1;
         found = TRUE;
         break;
      }
   }
   if (found) continue;
   ptrarray[whichsite] = 0;
}
#endif

#if !ALIASING
for(i = 0; i < dna->sites[locus]; i++) ptrarray[i] = 0;
#endif

data->siteptr = ptrarray;

} /* makesiteptr */


void read_line(FILE *source, char *line)
{
char *p;

fgets(line,LINESIZE,source);
while(line[0] == '\n' || line[0] == '#') fgets(line,LINESIZE,source);
if((p = (char *)strchr(line,'\n')) != NULL) *p = '\0';

while(isspace(*line)) *line++;
if(*line == '\0') read_line(source,line);

} /* read_line */


void read_header(option_struct *op, long *numpop, long *numloci,
   long **numsites, char *title, char *dlm)
{
long i;
char *input, *tok, temp;

input = (char *)calloc(LINESIZE,sizeof(char));

read_line(infile,input);

switch(lowercase(input[0])) {
   case 'a':
      break;
   case 'b':
   case 'm':
      sscanf(input,"%c%ld%ld%c%[^\n]",&temp,numpop,numloci,dlm,title);
      if (op->datatypeset)
         if (op->datatype != temp) {
            fprintf(ERRFILE,"\nYou chose a datatype of %c,",op->datatype);
            fprintf(ERRFILE,"but your infile shows a datatype of %c\n",temp);
            exit(-1);
         }
      op->datatype = temp;
      numsites = NULL;
      break;
   case 'n':
      if (op->autocorr) {
         printf("\nSNP data type is incompatible with autocorrelation");
         printf("\n----autocorrelation disabled-----\n");
         printf("press <enter> or <return> to continue\n");
         fgets(input,LINESIZE,stdin);
         fprintf(outfile,"\nSNP data type is incompatible with");
         fprintf(outfile," autocorrelation");
         fprintf(outfile,"\n----autocorrelation disabled-----\n");
         op->autocorr = FALSE;
      }
   case 's':
      sscanf(input,"%c%ld%ld%[^\n]",&temp,numpop,numloci,title);
      if (op->datatypeset)
         if (op->datatype != temp) {
            fprintf(ERRFILE,"\nYou chose a datatype of \"%c\",",op->datatype);
            fprintf(ERRFILE," but your infile shows a datatype");
            fprintf(ERRFILE," of \"%c\".\n",temp);
            exit(-1);
         }
      op->datatype = temp;
      read_line(infile,input);
      (*numsites) = (long *)calloc(*numloci,sizeof(long));
      for(i = 0; i < *numloci; i++) {
          while(isspace((int)*input)) (*input)++;
          if(i == 0) tok = strtok(input," ");
          else tok = strtok(NULL," ");
          (*numsites)[i] = atol(tok);
      }
      dlm = NULL;
      break;
   default:
      fprintf(ERRFILE,"WARNING--datatype not specified in infile\n");
      switch(op->datatype) {
         case 'a':
            break;
         case 'b':
            break;
         case 'm':
            sscanf(input,"%ld%ld%c%[^\n]",numpop,numloci,dlm,title);
            break;
         case 'n':
         case 's':
            sscanf(input,"%ld%ld%[^\n]",numpop,numloci,title);
            read_line(infile,input);
            (*numsites) = (long *)calloc(*numloci,sizeof(long));
            for(i = 0; i < *numloci; i++) {
                while(isspace((int)*input)) (*input)++;
                if(locus == 0) tok = strtok(input," ");
                else tok = strtok(NULL," ");
                (*numsites)[i] = atol(tok);
            }
            break;
         default:
            fprintf(ERRFILE,"ERROR:unknown datatype %c\n",op->datatype);
            fprintf(ERRFILE,"only the types a, m, s");
            fprintf(ERRFILE,"(electrophoretic alleles,\n");
            fprintf(ERRFILE,"microsatellite data, sequence data)");
            fprintf(ERRFILE,"are allowed.\n");
            exit(-1L);
            break;
      }
      break;
}

free(input);

} /* read_header */


void read_popheader(long *numind, char *poptitle)
{
char *input;

input = (char *)calloc(LINESIZE,sizeof(char));

read_line(infile,input);
sscanf(input,"%ld%[^\n]",numind,poptitle);

free(input);

} /* read_popheader */


void getinput(option_struct *op, data_fmt *data)
{
long i, numpop, numloci, numind, *numsites;
char *title, *poptitle, dlm;

title = (char *)calloc(LINESIZE,sizeof(char));
poptitle = (char *)calloc(LINESIZE,sizeof(char));

/* setup data structures and read in the first population's
   worth of data */
read_header(op,&numpop,&numloci,&numsites,title,&dlm);
read_popheader(&numind,poptitle);
setupdata(data,op->datatype,numpop,numloci,numind,numsites,title);
setpopstuff(data,0L,numind,poptitle,op->datatype);
read_popdata(infile,data,0L,op);
/* "true haplotype" code */
if (TRUEHAP) {
  readtruehaps(op,&truehaps,numpop,numloci,numind,numsites);
  copyhaps(op,data,&starthaps,numpop,numloci,numind,numsites);
}

/* now read in all other populations' data */
for(i = 1; i < numpop; i++) {
   read_popheader(&numind,poptitle);
   setpopstuff(data,i,numind,poptitle,op->datatype);
   read_popdata(infile,data,i,op);
}


if (op->datatype == 's' || op->datatype == 'n') {
   if (op->weights)
      inputdnaweights(data->dnaptr->sites[locus],data->dnaptr,op);
   free(numsites);
}

free(poptitle);
free(title);

}  /* getinput */


void getnodata(option_struct *op, data_fmt *data)
{
long i, numpop, numloci, numchains, numind, *numsites;
char *title, *poptitle, dlm;
FILE *kks;

/* debug DEBUG warning WARNING--bogus default settings */
op->datatype = 's';
numpop = 1;
numind = 1;
title = "no data run";
dlm = ' ';
poptitle = "first pop";

kks = fopen("kks","r");
fscanf(kks,"%ld %ld",&numloci,&numchains);
numsites = (long *)calloc(numloci,sizeof(long));

/* debug DEBUG warning WARNING--bogus default settings */
for(i = 0; i < numloci; i++) numsites[i] = 1;

setupdata(data,op->datatype,numpop,numloci,numind,numsites,title);
setpopstuff(data,0L,numind,poptitle,op->datatype);

fseek(kks,0,SEEK_SET);
fclose(kks);

} /* getnodata */


double watterson(option_struct *op, data_fmt *data)
{
  /* estimate theta using method of Watterson */
  long i, j, k, kn, numprivate, numsites;
  boolean varies, private;
  double watter;
  dnadata *dna;

  dna = data->dnaptr;
  numprivate = kn = 0;

  for (i = 0; i < dna->sites[locus]; i++) {
    varies = FALSE;
    private = TRUE;
    for (j = 1; j < dna->numseq[population]; j++) {
      if (dna->seqs[population][locus][j][i] !=
          dna->seqs[population][locus][0][i]) {
        if (varies) {
           for(k = 1; k < dna->numseq[population]; k++) {
              if (dna->seqs[population][locus][k][i] !=
                  dna->seqs[population][locus][j][i])
                 private = FALSE;
           }
        }
	varies = TRUE;
      }
    }
    if (varies) kn++;
    if (private && varies) numprivate++;
  }
  watter = 0.0;
  if (kn > 0) {
    for (i = 1; i < dna->numseq[population]; i++)
      watter += 1.0 / i;
    numsites = countsites(op,data);
    watter = kn / (numsites * watter);
    fprintf(outfile,"\nThere are %ld variable sites in the data\n",kn);

    if (kn == 1) {
       op->holding = 2;
       fprintf(outfile,"\nWARNING:  There is only 1 variable site in");
       fprintf(outfile," this data set.\nRecombination-rate cannot");
       fprintf(outfile," be accurately estimated in this case.\n");
       fprintf(outfile,"The \"estimate\" of rec-rate has been fixed");
       fprintf(outfile," to its starting value.\n\n");

       printf("\nWARNING:  There is only 1 variable site in");
       printf(" this data set.\nRecombination-rate cannot");
       printf(" be accurately estimated in this case.\n");
       printf("The \"estimate\" of rec-rate has been fixed");
       printf(" to its starting value.\n\n");
    } else {
       if (kn-numprivate < 2) {
          op->holding = 2;
          fprintf(outfile,"\nWARNING:  Too many of the variable sites ");
          fprintf(outfile,"are uninformative.\nRecombination-rate ");
          fprintf(outfile,"cannot be accurately estimated in this ");
          fprintf(outfile,"case.\nThe \"estimate\" of rec-rate has ");
          fprintf(outfile,"been fixed to its starting value.\n\n");

          printf("\nWARNING:  Too many of the variable sites ");
          printf("are uninformative.\nRecombination-rate ");
          printf("cannot be accurately estimated in this ");
          printf("case.\nThe \"estimate\" of rec-rate has ");
          printf("been fixed to its starting value.\n\n");
       }
    }

    return watter;
  }
  fprintf(outfile, "Warning:  There are no variable sites");
  fprintf(outfile, " in this data set.\n\n");
  if (MENU) printf("Warning:  There are no variable sites in this data set.\n");
  else {
     fprintf(simlog, "Warning:  There are no variable sites");
     fprintf(simlog, " in this data set.\n\n");
  }
  exit(-1);
}  /* watterson */


/**********************************************************************
 * invardatacheck() is called to check non-panel SNP data for         *
 * invariant sites, which are illegal under that model.  If it finds  *
 * any, it changes that site to a datatype of "unknown"; "x" is used. */
void invardatacheck(option_struct *op, data_fmt *data)
{
long marker, seq, nummarkers, numseq;
char **dna;
boolean varies;

nummarkers = getdata_nummarkers(op,data);
numseq = getdata_numseq(op,data);

dna = data->dnaptr->seqs[population][locus];

for(marker = 0; marker < nummarkers; marker++) {
   varies = FALSE;
   for(seq = 1; seq < numseq; seq++) {
      if(dna[0][marker] != dna[seq][marker]) {
         varies = TRUE;
         break;
      }
   }
   if (!varies)
      for(seq = 0; seq < numseq; seq++) dna[seq][marker] = 'X';
}

} /* invardatacheck */


/*********************************************************************
 * buildtymelist() allocates space and initializes all the fields of *
 * a tymelist except for failure to initialize the branchlists.      *
 * Those are handled by finishsetup().                               *
 *                                                                   *
 * WARNING buildtymelist() currently assumes a non-recombinant tree! */
void buildtymelist(tree *tr, option_struct *op, data_fmt *data, node *p)
{
long numentries;
tlist *t, *u;

t = tr->tymelist;

newtymenode(&u);

numentries = getdata_numtips(op,data);
u->branchlist = (node **)calloc(numentries,sizeof(node *));

u->eventnode = p;

while (t != NULL) {
   if (u->eventnode->tyme < t->eventnode->tyme) {
      u->prev = t->prev;
      t->prev = u;
      u->succ = t;
      u->prev->succ = u;
      break;
   }
   if (t->succ != NULL)
      t = t->succ;
   else {
      t->succ = u;
      u->prev = t;
      u->succ = NULL;
      break;
   }
}

} /* buildtymelist */


/* WARNING warning orient assumes a non-recombinant tree */
void orient(tree *tr, option_struct *op, data_fmt *data, node *p)
{
long lastsite;

lastsite = countsites(op,data)-1;

  if (istip(p)) {
    return;
  }

  curtree->nodep[p->number] = p;  /* insure that curtree->nodep points
                                    to nodes with info */

 /* since p is a top nodelet, it needs to actually store
    likelihood information, x is a NULL pointer
    in all other non-tip nodelets */
  p->top = TRUE;
  allocate_x(op,data,p);
  if (op->map) {
     p->z = NULL;
     allocate_z(op,data,p);
  }
  ranges_Malloc(p,TRUE,1L);
if (op->fc) {
  coal_Malloc(p,TRUE,1L);
  p->ranges[1] = p->coal[1] = 0;
  p->ranges[2] = p->coal[2] = lastsite;
} else {
  p->ranges[1] = 0;
  p->ranges[2] = lastsite;
}
  p->next->top = FALSE;
  free_x(op,p->next);
  ranges_Malloc(p->next,FALSE,0L);
  coal_Malloc(p->next,FALSE,0L);
  p->next->next->top = FALSE;
  free_x(op,p->next->next);
  ranges_Malloc(p->next->next,FALSE,0L);
  coal_Malloc(p->next->next,FALSE,0L);

  orient(tr,op,data,p->next->back);
  orient(tr,op,data,p->next->next->back);

  p->tyme = p->next->length + p->next->back->tyme;
  p->next->tyme = p->tyme;
  p->next->next->tyme = p->tyme;
  if (p->number == curtree->root->back->number) {
    p->back->top = FALSE;
    p->back->tyme = ROOTLENGTH;
  }

  buildtymelist(tr,op,data,p);

}  /* orient */

void plumptree(option_struct *op, data_fmt *data, double th0)
/* change an input tree into a perfect coalescent tree for theta0 */
/* WARNING--doesn't work on a tree with recombinations!!!!! */
{
  tlist *t;
  long k;
  double tyme;

  /* we assume that tips are tyme 0 already */
  t = curtree->tymelist->succ;

  k = getdata_numtips(op,data);
  tyme = 0.0;
  do {
    tyme += th0/(k*(k-1));
    if (!istip(t->eventnode)) {
      t->eventnode->tyme = tyme;
      t->eventnode->next->tyme = tyme;
      t->eventnode->next->next->tyme = tyme;
    }
    k--;
    t = t->succ;
  } while (t != NULL);

  /* now you must call finishsetup and initbranchlist or a disaster
  will happen! */
} /* plumptree */

void finishsetup(option_struct *op, data_fmt *data, node *p)
{
  if (istip(p)) {
    ltov(op,data,p);
    return;
  }
  ltov(op,data,p);
  finishsetup(op,data,p->next->back);
  finishsetup(op,data,p->next->next->back);
  return;
} /* finishsetup */

void initbranchlist(option_struct *op, data_fmt *data)
{
  tlist *t;
  node *p, *q;
  long i, j, k, n, numtips;

  t = curtree->tymelist;
  numtips = n = getdata_numtips(op,data);
  t->numbranch = numtips;
  for(i = 0; i < numtips; i++) t->branchlist[i] = curtree->nodep[i+1];
  t->age = t->succ->eventnode->tyme;
  t = t->succ;
  /* WARNING assumes initial tree is not recombinant! */
  for (i = 0; i < numtips-1; i++) {
    /* for each interior node, do... */
    n--;
    t->numbranch = n;
    if (n == 1)
      t->age = t->eventnode->tyme + ROOTLENGTH;
    else
      t->age = t->succ->eventnode->tyme;
    p = t->eventnode->next->back;
    q = t->eventnode->next->next->back;
    k = 0;
    for (j = 0; j < t->prev->numbranch ; j++) {
      /* for the number of branches above the coalescent node, do...*/
      if (t->prev->branchlist[j] != p && t->prev->branchlist[j] != q) {
	t->branchlist[k] = t->prev->branchlist[j];
	k++;
      }
    }
    t->branchlist[t->numbranch - 1] = t->eventnode;
    t = t->succ;
  }
}  /* initbranchlist */

void inittable(option_struct *op, dnadata *dna)
{
  long i;
  tbl = (valrec *)calloc(1,op->categs*sizeof(valrec));
  /* Define a lookup table. Precompute values and store them in a table */
  for (i = 0; i < op->categs; i++) {
    tbl[i].rat_xi = op->rate[i] * dna->xi;
    tbl[i].rat_xv = op->rate[i] * dna->xv;
  }
}  /* inittable */

void initweightrat(option_struct *op, data_fmt *data)
{
/* WARNING non-DNA data is not yet able to do this! */
  long i, nummarkers;

  nummarkers = getdata_nummarkers(op,data);
  weightrat = (double *)calloc(nummarkers,sizeof(double));
  sumweightrat = 0.0;
  for (i = 0; i < nummarkers; i++) {
    weightrat[i] = data->dnaptr->dnaweight[i] * op->rate[category[i]-1];
    sumweightrat += weightrat[i];
  }
}  /* initweightrat */

void rec_outtree(node *p, boolean first, FILE **usefile)
/*  write out a recombinant tree in KYB format; call with
    curtree.root->back.  Calling program should append a
    semicolon.  "first" is used to avoid an extra comma
    in the tree list and should be TRUE when this is called. */
{
  long i;

  if(istip(p)) {
    if(!first)fprintf(*usefile,",");
    else first=FALSE;
    fprintf(*usefile,"\"");
    for(i=0;i<NMLNGTH;i++) {
      if(p->nayme[i]==' ') break;
      else putc(p->nayme[i],*usefile);
    }
    fprintf(*usefile,"\":%f",p->length);
    return;
  }
  if(isrecomb(p)) {
    if(!first)fprintf(*usefile,",");
    else first=FALSE;
    fprintf(*usefile,"%ld",p->number);
    fprintf(*usefile,"{%ld-%ld}",p->recstart,p->recend);
    fprintf(*usefile,":%f",p->length);
    p->updated = !p->updated;
    if (p->updated) {
      if (p->next->top) rec_outtree(p->next->back, first, usefile);
      else rec_outtree(p->next->next->back, first, usefile);
    }
    return;
  } else { /* p is coalescent */
    if(!first)fprintf(*usefile,",");
    else first=FALSE;
    fprintf(*usefile,"*%ld",p->number);
    fprintf(*usefile,":%f",p->length);
    rec_outtree(p->next->back, first, usefile);
    rec_outtree(p->next->next->back, first, usefile);
    return;
  }
} /* rec_outtree */

void treeout(node *p, long s, FILE **usefile)
{
  /* write out file with representation of final tree */
  long i, n, w;
  char c;
  double x;

  if (istip(p)) {
    n = 0;
    for (i = 1; i <= NMLNGTH; i++) {
      if (p->nayme[i - 1] != ' ')
	n = i;
    }
    for (i = 0; i < n; i++) {
      c = p->nayme[i];
      if (c == ' ')
	c = '_';
      putc(c, *usefile);
    }
    col += n;
  } else {
    putc('(', *usefile);
    col++;
    treeout(p->next->back, s, usefile);
    putc(',', *usefile);
    col++;
    if (col > 45) {
      putc('\n', *usefile);
      col = 0;
    }
    treeout(p->next->next->back, s, usefile);
    putc(')', *usefile);
    col++;
  }
  if (p->v >= 1.0)
    x = -1.0;
  else
    x = lengthof(p);
  if (x > 0.0)
    w = (long)(0.4343 * log(x));
  else if (x == 0.0)
    w = 0;
  else
    w = (long)(0.4343 * log(-x)) + 1;
  if (w < 0)
    w = 0;
  if (p == curtree->root->back)
    putc(';', *usefile);
  else {
    fprintf(*usefile, ":%*.10f", (int)(w + 7), x);
    col += w + 8;
  }
}  /* treeout */


double snp_eval_calcrange(option_struct *op, data_fmt *data, tree *tr,
   boolean first, long start, long finish, long numcategs)
{
double sumterm, sumterm2, lterm, lterm2, temp;
long i, j, k, slice, snpstart, snpfinish, nummarkers, numinvar, *siteptr;
contribarr term, term2, clai;
node *p;
double *x1;
dnadata *dna;

dna = data->dnaptr;

term = (double *)calloc(numcategs,sizeof(double));
term2 = (double *)calloc(numcategs,sizeof(double));
clai = (double *)calloc(numcategs,sizeof(double));

slice = 0;
temp = 0.0;
p = tr->root->back;
siteptr = data->siteptr;

nummarkers = getdata_nummarkers(op, data);
findsubtreemarkers(op,data,start,finish,&snpstart,&snpfinish);

if (snpstart == FLAGLONG) numinvar = finish-start+1;
else numinvar = (finish-start+1) - (snpfinish-snpstart+1);

sumterm2 = 0.0;
for (j = 0; j < numcategs; j++) { /* compute P(Invar) */
   term2[j] = 0.0;
   for(k = nummarkers; k < nummarkers+NUMINVAR; k++) {
      x1 = p->x->s[k][j][slice];
      term2[j] += dna->freqa * x1[baseA] +
          dna->freqc * x1[baseC] + dna->freqg * x1[baseG] +
          dna->freqt * x1[baseT];
   }
   sumterm2 += op->probcat[j] * term2[j];
}
lterm2 = log(sumterm2);
for (j = 0; j < numcategs; j++) {
   clai[j] = term2[j] / sumterm2;
}
memcpy(contribution[nummarkers], clai, numcategs*sizeof(double));
if (!op->autocorr)
   alogf->val[nummarkers] = lterm2;
if (!op->norecsnp) temp += numinvar * lterm2;

if (snpstart != FLAGLONG) { /* that is, if there are SNPs */
   for (i = snpstart; i <= snpfinish; i++) {
      for (j = 0; j < numcategs; j++) {
         x1 = p->x->s[i][j][slice];
         term[j] = dna->freqa * x1[baseA] + dna->freqc * x1[baseC] +
            dna->freqg * x1[baseG] + dna->freqt * x1[baseT];
         if(term[j] == 0) {
            fprintf(ERRFILE,"Encountered tree incompatible with data\n");
            if(first) {
               fprintf(ERRFILE,"starting tree needs to be legal\n");
               exit(-1);
            }
            curtree->likelihood = NEGMAX;
            return(-1);
         }
      }
      sumterm = 0.0;
      for (j = 0; j < numcategs; j++) {
         sumterm += op->probcat[j] * term[j];
      }
      lterm = log(sumterm);
      if (op->norecsnp && !op->full) lterm -= log(1.0-sumterm2);
      for (j = 0; j < numcategs; j++) {
            clai[j] = ((op->norecsnp) ?
               (term[j] / (1.0-term2[j])) / (sumterm / (1.0-sumterm2))
               : (term[j] / sumterm));
      }
      memcpy(contribution[i], clai, numcategs*sizeof(double));
      if (!op->autocorr)
         alogf->val[i] = lterm;
      temp += dna->dnaweight[i] * lterm;
   }

#if 0 /* This is probably dead DEAD dead--warning WARNING DEBUG debug */
   if (op->datatype == 'n' && op->full) {
      temp += log(op->chance_seen);
      numunobs = getnumlinks(op,data,start,finish)-(finish-start);
      if (numunobs < 0 || (long)numunobs != numunobs) {
         fprintf(ERRFILE,"non integral SNP distances with Full model\n");
         exit(-1);
      }
      temp += numunobs * log(1.0 - (op->chance_seen*(1.0-sumterm2)));
   }
#endif
}

free(term);
free(term2);
free(clai);

return(temp);

} /* snp_eval_calcrange */


double panel_eval_calcrange(option_struct *op, data_fmt *data, tree *tr,
   boolean first, long start, long finish, long numcategs)
{
double sumterm, sumterm2, lterm, lterm2, temp, invarterm;
long i, j, k, numunobs, nummarkers, numinvar, snpstart, snpfinish,
   *siteptr;
contribarr term, term2, clai;
node *p;
double **x1;
dnadata *dna;

dna = data->dnaptr;

term = (double *)calloc(numcategs,sizeof(double));
term2 = (double *)calloc(numcategs,sizeof(double));
clai = (double *)calloc(numcategs,sizeof(double));

temp = 0.0;
p = tr->root->back;
siteptr = data->siteptr;

/* compute P(INVAR) for each category, used as P(SNP) = 1.0 - P(INVAR) */
nummarkers = getdata_nummarkers(op,data);
findsubtreemarkers(op,data,start,finish,&snpstart,&snpfinish);

if (snpstart == FLAGLONG) numinvar = finish-start+1;
else numinvar = (finish-start+1) - (snpfinish-snpstart+1);

sumterm2 = 0.0;
for (j = 0; j < numcategs; j++) {
   x1 = p->x->s[nummarkers][j];
   term2[j] = 0.0;
    for (k = 1; k < NUMSLICE; k++) {
      term2[j] += dna->freqa * x1[k][baseA] +
          dna->freqc * x1[k][baseC] + dna->freqg * x1[k][baseG] +
          dna->freqt * x1[k][baseT];
   }
   sumterm2 += op->probcat[j] * term2[j];
}
lterm2 = log(sumterm2);
for (j = 0; j < numcategs; j++) {
   clai[j] = term2[j] / sumterm2;
}
memcpy(contribution[nummarkers], clai, numcategs*sizeof(double));
if (!op->autocorr)
   alogf->val[nummarkers] = lterm2;
if (!op->norecsnp) temp += numinvar * lterm2;

if (snpstart != FLAGLONG) {
   for (i = snpstart; i <= snpfinish; i++) {
      for (j = 0; j < numcategs; j++) {
      /* compute P(data) for this category */
         x1 = p->x->s[i][j];
         term[j] = dna->freqa * x1[0][baseA] + dna->freqc * x1[0][baseC] +
           dna->freqg * x1[0][baseG] + dna->freqt * x1[0][baseT];
         invarterm = 0.0;
      /* compute P(data && panel invariant) for this category */
         for (k = 1; k < NUMSLICE; k++) {
            invarterm += dna->freqa * x1[k][baseA] + dna->freqc *
              x1[k][baseC] + dna->freqg * x1[k][baseG] + dna->freqt *
              x1[k][baseT];
         }
         term[j] -= invarterm;
         if(term[j] == 0) {
            fprintf(ERRFILE,"Encountered tree incompatible with data\n");
            if(first) {
               fprintf(ERRFILE,"starting tree needs to be legal\n");
               exit(-1);
            }
            curtree->likelihood = NEGMAX;
            return(-1);
         }
      }
      sumterm = 0.0;
      for (j = 0; j < numcategs; j++) {
         sumterm += op->probcat[j] * term[j];
      }
      lterm = log(sumterm);
      if (op->norecsnp && !op->full) lterm -= log(1.0-sumterm2);
      for (j = 0; j < numcategs; j++) {
         clai[j] = ((op->norecsnp) ?
            (term[j] / (1.0 - term2[j])) / (sumterm / (1.0 - sumterm2))
            : (term[j] / sumterm));
      }
      memcpy(contribution[i], clai, numcategs*sizeof(double));
      if (!op->autocorr)
         alogf->val[i] = lterm;
      temp += dna->dnaweight[i] * lterm;
   }

   if (op->full) {
      temp += log(op->chance_seen);
      numunobs = (long)(getnumlinks(op,data,start,finish)-(finish-start));
      if (numunobs < 0 || (long)numunobs != numunobs) {
         fprintf(ERRFILE,"non integral SNP distances with Full model\n");
         exit(-1);
      }
      temp += numunobs * log(1.0 - (op->chance_seen*(1.0-sumterm2)));
   }
}

free(term);
free(term2);
free(clai);

return(temp);

} /* panel_eval_calcrange */


double eval_calcrange(option_struct *op, data_fmt *data, tree *tr,
   boolean first, long start, long finish, long numcategs)
{
double sumterm, lterm, temp;
long i, j, slice, nummarkers, *siteptr;
contribarr term, term2, clai;
node *p;
double *x1;
dnadata *dna;

dna = data->dnaptr;

term = (double *)calloc(numcategs,sizeof(double));
term2 = (double *)calloc(numcategs,sizeof(double));
clai = (double *)calloc(numcategs,sizeof(double));

slice = 0;
temp = 0.0;
p = tr->root->back;
siteptr = data->siteptr;
nummarkers = getdata_nummarkers(op,data);

for (i = start; i <= finish; i++) {
   for (j = 0; j < numcategs; j++) {
      x1 = p->x->s[i][j][slice];
      term[j] = dna->freqa * x1[baseA] + dna->freqc * x1[baseC] +
         dna->freqg * x1[baseG] + dna->freqt * x1[baseT];
      if(term[j] == 0) {
         fprintf(ERRFILE,"Encountered tree incompatible with data\n");
         if(first) {
            fprintf(ERRFILE,"starting tree needs to be legal\n");
            exit(-1);
         }
         curtree->likelihood = NEGMAX;
         return(-1);
      }
   }
   sumterm = 0.0;
   for (j = 0; j < numcategs; j++) {
      sumterm += op->probcat[j] * term[j];
   }
   lterm = log(sumterm);
   for (j = 0; j < numcategs; j++) {
        clai[j] = term[j] / sumterm;
   }
   memcpy(contribution[i], clai, numcategs*sizeof(double));
   if (!op->autocorr)
      alogf->val[i] = lterm;
   temp += dna->dnaweight[i] * lterm;
}

free(term);
free(term2);
free(clai);

return(temp);

} /* eval_calcrange */


double evaluate(option_struct *op, data_fmt *data, tree *tr,
   double llike)
{
double sum2, sumc, one_minus_lambda;
contribarr like, nulike, clai;
long i, j, k, nummarkers;

if (op->categs > 1) {
   nummarkers = getdata_nummarkers(op,data);

   like = (double *)calloc(op->categs,sizeof(double));
   nulike = (double *)calloc(op->categs,sizeof(double));

   one_minus_lambda = 1.0 - op->lambda;

   for (j = 0; j < op->categs; j++) like[j] = 1.0;
   for (i = 0; i < nummarkers; i++) {
      sumc = 0.0;
      for (k = 1; k <= op->categs; k++)
         sumc += op->probcat[k - 1] * like[k - 1];
      sumc *= op->lambda;
      clai = contribution[i];
      for (j = 0; j < op->categs; j++)
         nulike[j] = (one_minus_lambda * like[j] + sumc) * clai[j];
      memcpy(like, nulike, op->categs*sizeof(double));
   }
   sum2 = 0.0;
   for (i = 0; i < op->categs; i++)
      sum2 += op->probcat[i] * like[i];
   llike += log(sum2);
   tr->likelihood = llike;

   free(like);
   free(nulike);

   return(llike);
} else {
   tr->likelihood = llike;
   return(llike);
}

}  /* evaluate */


double snp_evaluate(option_struct *op, data_fmt *data, tree *tr,
   double llike, long start, long finish)
{
double sum2, sumc;
contribarr like, nulike, clai;
long i, j, snpstart, snpfinish, nummarkers, numinvar;

if (op->categs > 1) {
   nummarkers = getdata_nummarkers(op,data);
   findsubtreemarkers(op,data,start,finish,&snpstart,&snpfinish);

   like = (double *)calloc(op->categs,sizeof(double));
   nulike = (double *)calloc(op->categs,sizeof(double));

   for (j = 0; j < op->categs; j++) like[j] = 1.0;

   if (snpstart == FLAGLONG) numinvar = finish-start+1;
   else numinvar = (finish-start+1) - (snpfinish-snpstart+1);

   if (snpstart != FLAGLONG) {
      for (i = snpstart; i <= snpfinish; i++) {
         sumc = 0.0;
         for (j = 0; j < op->categs; j++)
            sumc += op->probcat[j] * like[j];
         clai = contribution[i];
         for (j = 0; j < op->categs; j++)
            nulike[j] = sumc * clai[j];
         memcpy(like, nulike, op->categs*sizeof(double));
      }
   }

   for (i = 0; i < numinvar; i++) { /* invariant sites */
      sumc = 0.0;
      for (j = 0; j < op->categs; j++)
         sumc += op->probcat[j] * like[j];
      clai = contribution[nummarkers];
      for (j = 0; j < op->categs; j++)
         nulike[j] = sumc * clai[j];
      memcpy(like, nulike, op->categs*sizeof(double));
   }

   sum2 = 0.0;
   for (j = 0; j < op->categs; j++)
      sum2 += op->probcat[j] * like[j];
   llike += log(sum2);

   free(like);
   free(nulike);

   tr->likelihood = llike;
   return(llike);
} else {
   tr->likelihood = llike;
   return(llike);
}

}  /* snp_evaluate */


boolean nuview_usebranch(data_fmt *data, node *p, long site)
{

if(!p->top) return(FALSE);

return(inrange(p->ranges,site));

} /* nuview_usebranch */


double prob_micro(msatdata *ms, double t, long diff)
{
double **steps = ms->steps;
long stepnum = MICRO_MAXCHANGE;
long k;
double newsum = 0.0, oldsum=0.0;
double logt = log(t);

if(diff>=stepnum)
     return newsum;
for (k = diff; k < diff + stepnum; k+=2){
  newsum += exp(-t + logt * k - steps[diff][k-diff]);
      if(oldsum-newsum < DBL_EPSILON)
             break;
       oldsum = newsum;
}
return newsum;

} /* prob_micro */


/********************************************************
 * findcoal() takes a non-top coalescent nodelet and    *
 * returns the closest tipwards top coalescent nodelet. *
 * It also returns the "v" value between the two nodes. *
 * In the case of failure it returns a NULL pointer     *
 * with v set to FLAGLONG.                              *
 *                                                      *
 * findcoal assumes that branchlengths are additive.    */
node *findcoal(option_struct *op, data_fmt *data, node *p, double *v)
{
double total;
node *q;

if(!iscoal(p) || p->top) {*v = FLAGLONG; return(NULL);}

for (q = p->back, total = 0; !iscoal(q) && !istip(q); q = findunique(q)->back)
   total += lengthof(q);

total += lengthof(q);
*v = findcoal_ltov(op,data,total);

return(q);

} /* findcoal */


void nuview_micro(option_struct *op, data_fmt *data, node *p,
   long indexsite, long startsite, long endsite)
{
long i, a, s, smax, margin, diff;
double *xx1, *xx2, *xx3, vv1 = 0.0, vv2 = 0.0, pija1s, pija2s,
   lx1, lx2, x3m;
node *q, *r;
boolean qactive, ractive;

smax = MICRO_ALLELEMAX;
margin = MICRO_MAXCHANGE;
x3m = NEGMAX;

xx1 = NULL;
xx2 = NULL;
vv1 = 0.0;
vv2 = 0.0;
xx3 = (double *)calloc(smax,sizeof(double));

q = findcoal(op,data,p->next,&vv1);
r = findcoal(op,data,p->next->next,&vv2);

/* are these sites active and on an upwards branch?  */
qactive = nuview_usebranch(data,p->next,indexsite);
ractive = nuview_usebranch(data,p->next->next->back,indexsite);

for(i = startsite; i <= endsite; i++) {
   if (qactive) {
      xx1 = q->x->a[i];
      lx1 = q->lxmax;
   } else {
      lx1 = 0.0;
   }

   if (ractive) {
      xx2 = r->x->a[i];
      lx2 = r->lxmax;
   } else {
      lx2 = 0.0;
   }

   for(s = 0; s < smax; s++) {
      if (qactive) pija1s = 0.0;
      else pija1s = 1.0;
      if (ractive) pija2s = 0.0;
      else pija2s = 1.0;
      for(a = MAX(0,s-margin); a < s + margin && a < smax; a++) {
         diff = labs(s-a);
         if(xx1[a] > 0 && qactive) {
            pija1s += prob_micro(data->msptr,vv1,diff) * xx1[a];
         }
         if(xx2[a] > 0 && ractive) {
            pija2s += prob_micro(data->msptr,vv2,diff) * xx2[a];
         }
      }
      xx3[s] = pija1s * pija2s;
      if(xx3[s] > x3m) x3m = xx3[s];
   }
   if(x3m == 0.0) p->lxmax = NEGMAX;
   else {
      for(s = 0; s < smax; s++) xx3[s] /= x3m;
      p->lxmax = log(x3m) + lx1 + lx2;
   }
   memcpy(p->x->a[i],xx3,smax*sizeof(double));
}

free(xx3);

} /* nuview_micro */


void nuview(option_struct *op, data_fmt *data, node *p, long indexsite,
   long startsite, long endsite)
{
long i;
double lw1 = 0.0, lw2 = 0.0;
node *q, *r;
boolean toobig;

/* set "q" and "r" for use in calcrange(), either to a valid nodelet with
   datalikelihoods, or NULL if there is no such nodelet. */
if (nuview_usebranch(data,p->next->back,indexsite)) {
   q = findcoal(op,data,p->next,&lw1);
   toobig = (exp(lw1) <= 0.0);
} else {q = NULL; toobig = TRUE;}

if (toobig) {
   for (i = 0; i < op->categs; i++) {
      tbl[i].ww1 = tbl[i].zz1 = 0.0;
      tbl[i].ww1zz1 = tbl[i].vv1zz1 = 0.0;
   }
} else {
   for (i = 0; i < op->categs; i++) {
      tbl[i].ww1 = exp(tbl[i].rat_xi * lw1);
      tbl[i].zz1 = exp(tbl[i].rat_xv * lw1);
      tbl[i].ww1zz1 = tbl[i].ww1 * tbl[i].zz1;
      tbl[i].vv1zz1 = (1.0 - tbl[i].ww1) * tbl[i].zz1;
   }
}

if (nuview_usebranch(data,p->next->next->back,indexsite)) {
   r = findcoal(op,data,p->next->next,&lw2);
   toobig = (exp(lw2) <= 0.0);
} else {r = NULL; toobig = TRUE;}

if (toobig) {
   for (i = 0; i < op->categs; i++) {
      tbl[i].ww2 = tbl[i].zz2 = 0.0;
      tbl[i].ww2zz2 = tbl[i].vv2zz2 = 0.0;
   }
} else {
   for (i = 0; i < op->categs; i++) {
      tbl[i].ww2 = exp(tbl[i].rat_xi * lw2);
      tbl[i].zz2 = exp(tbl[i].rat_xv * lw2);
      tbl[i].ww2zz2 = tbl[i].ww2 * tbl[i].zz2;
      tbl[i].vv2zz2 = (1.0 - tbl[i].ww2) * tbl[i].zz2;
   }
}

calcrange(op,data,p,q,r,indexsite,startsite,endsite,op->categs);

}  /* nuview */


void calcrange(option_struct *op, data_fmt *data, node *p, node * q,
   node *r, long whichtree, long start, long finish, long numcategs)
/* whichtree indicates which tree should be used to calculate the
   likelihood of this set of sites:  pass the number of a site which
   has the desired tree.  (This is helpful in calculations for the
   SNP invariant sites.) */
{
long i, j, k, numslice, *siteptr;
double yy1, yy2, ww1zz1, vv1zz1, ww2zz2, vv2zz2, vv1zz1_sumr1,
   vv2zz2_sumr2, vv1zz1_sumy1, vv2zz2_sumy2, sum1, sum2,
   sumr1, sumr2, sumy1, sumy2;
boolean qactive, ractive;
dnadata *dna;
double **allones, **xx1, **xx2, **xx3;
long istart, ifinish;

// DEBUG
double temptotal;

dna = data->dnaptr;
numslice = (op->panel) ? NUMSLICE : 1L;
siteptr = data->siteptr;

if (op->datatype == 'n') {
  if (start == FLAGLONG) { /* extra sites */
     istart = getdata_nummarkers(op,data);
     ifinish = istart+NUMINVAR-1;
  } else {
     findsubtreemarkers(op,data,start,finish,&istart,&ifinish);
     if(istart==FLAGLONG) return;
  }
} else {
  istart = start;
  ifinish = finish;
}

/* allocate some working space--this should probably be moved
upstream at some point to save time! */

allones = (double **)calloc(numslice,sizeof(double *));
xx3 = (double **)calloc(numslice,sizeof(double *));
allones[0] = (double *)calloc(numslice*4,sizeof(double));
xx3[0] = (double *)calloc(numslice*4,sizeof(double));
for (i = 0; i < numslice; i++) {
  allones[i] = allones[0] + i * 4;
  xx3[i] = xx3[0] + i * 4;
  for (j = 0; j < 4; j++)
    allones[i][j] = 1.0;
}

/* are these sites active and on an upwards branch?  */
qactive = (q) ? TRUE : FALSE;
ractive = (r) ? TRUE : FALSE;


for(i = istart; i <= ifinish; i++) {
   if(siteptr[i] > 0) { /* is an alias site available? */
      memcpy(p->x->s[i][0][0], p->x->s[siteptr[i]-1][0][0],
        op->categs*numslice*4L*sizeof(double));
// DEBUG
      temptotal = p->x->s[i][0][0][baseA] + p->x->s[i][0][0][baseC] + p->x->s[i][0][0][baseG] + p->x->s[i][0][0][baseT];
      if (temptotal == 0.0)
        printf("bad likelihood in memcpy!");
//
      continue;  /* go to next site */
   }

   for (j = 0; j < numcategs; j++) {
      memcpy(xx3[0], allones[0], numslice * 4L * sizeof(double));
      if (qactive) {
          ww1zz1 = tbl[j].ww1zz1;
          vv1zz1 = tbl[j].vv1zz1;
          yy1 = 1.0 - tbl[j].zz1;
          xx1 = q->x->s[i][j];

          for (k=0; k<numslice; k++) {
            sum1 = yy1 * (dna->freqa * xx1[k][baseA] +
              dna->freqc * xx1[k][baseC] +
    	      dna->freqg * xx1[k][baseG] + dna->freqt * xx1[k][baseT]);
            sumr1 = dna->freqar * xx1[k][baseA] + dna->freqgr * xx1[k][baseG];
            sumy1 = dna->freqcy * xx1[k][baseC] + dna->freqty * xx1[k][baseT];
            vv1zz1_sumr1 = vv1zz1 * sumr1;
            vv1zz1_sumy1 = vv1zz1 * sumy1;
            xx3[k][baseA] *=
              (sum1 + ww1zz1 * xx1[k][baseA] + vv1zz1_sumr1);
            xx3[k][baseC] *=
      	      (sum1 + ww1zz1 * xx1[k][baseC] + vv1zz1_sumy1);
            xx3[k][baseG] *=
      	      (sum1 + ww1zz1 * xx1[k][baseG] + vv1zz1_sumr1);
            xx3[k][baseT] *=
      	      (sum1 + ww1zz1 * xx1[k][baseT] + vv1zz1_sumy1);
// DEBUG
            temptotal = xx3[k][baseA] + xx3[k][baseC] + xx3[k][baseG] + xx3[k][baseT];
            if (temptotal == 0.0)
              printf("bad likelihood in q!");
//
           }
        }

        if (ractive) {
          ww2zz2 = tbl[j].ww2zz2;
          vv2zz2 = tbl[j].vv2zz2;
          yy2 = 1.0 - tbl[j].zz2;
          xx2 = r->x->s[i][j];

          for (k=0; k<numslice; k++) {
            sum2 = yy2 * (dna->freqa * xx2[k][baseA] +
              dna->freqc * xx2[k][baseC] +
      	      dna->freqg * xx2[k][baseG] + dna->freqt * xx2[k][baseT]);
            sumr2 = dna->freqar * xx2[k][baseA] + dna->freqgr * xx2[k][baseG];
            sumy2 = dna->freqcy * xx2[k][baseC] + dna->freqty * xx2[k][baseT];
            vv2zz2_sumr2 = vv2zz2 * sumr2;
            vv2zz2_sumy2 = vv2zz2 * sumy2;
            xx3[k][baseA] *=
              (sum2 + ww2zz2 * xx2[k][baseA] + vv2zz2_sumr2);
            xx3[k][baseC] *=
              (sum2 + ww2zz2 * xx2[k][baseC] + vv2zz2_sumy2);
            xx3[k][baseG] *=
              (sum2 + ww2zz2 * xx2[k][baseG] + vv2zz2_sumr2);
            xx3[k][baseT] *=
              (sum2 + ww2zz2 * xx2[k][baseT] + vv2zz2_sumy2);
// DEBUG
            temptotal = xx3[k][baseA] + xx3[k][baseC] + xx3[k][baseG] + xx3[k][baseT];
            if (temptotal == 0.0)
              printf("bad likelihood in r!");
//
          }
        }

      memcpy(p->x->s[i][j][0], xx3[0], numslice * 4L * sizeof(double));
    }
}

free(allones[0]);
free(allones);
free(xx3[0]);
free(xx3);

} /* calcrange */


void localsmooth(option_struct *op, data_fmt *data, node *p,
   long indexsite, long startsite, long endsite)
{
/* never mind if it's a tip */
if (istip(p)) return;

/* recurse left */
if (!p->next->top)
  if (!p->next->back->updated)
     if (inrange(p->next->back->ranges,indexsite))
       localsmooth(op,data,p->next->back,indexsite,startsite,endsite);

/* recurse right */
if (!p->next->next->top)
  if (!p->next->next->back->updated)
     if (inrange(p->next->next->back->ranges,indexsite))
       localsmooth(op,data,p->next->next->back,indexsite,startsite,
         endsite);

/* update likelihoods */
if (!p->updated && iscoal(p))
  switch(op->datatype) {
  case 'a':
     break;
  case 'b':
  case 'm':
     nuview_micro(op,data,p,indexsite,startsite,endsite);
     break;
  case 'n':
  case 's':
     nuview(op,data,p,indexsite,startsite,endsite);
     break;
  default:
     fprintf(ERRFILE,"localsmooth:unknown datatype %c\n",op->datatype);
     fprintf(ERRFILE,"only the types a, m, s");
     fprintf(ERRFILE,"(electrophoretic alleles,\n");
     fprintf(ERRFILE,"microsatellite data, sequence data)");
     fprintf(ERRFILE,"are allowed.\n");
     exit(-1);
  }

} /* localsmooth */


void snpsmooth (option_struct *op, data_fmt *data, tree *tr,
   long indexsite, long startsite, long endsite)
{
tlist *t;

for(t = tr->tymelist->succ; t != NULL; t = t->succ)
   if (iscoal(t->eventnode))
      nuview(op,data,t->eventnode,indexsite,startsite,endsite);

} /* snpsmooth */


void localeval(option_struct *op, data_fmt *data, node *p, boolean first)
{
long subtree, *subtree_ranges, substart, subend, lastmarker, *siteptr;
double dummy, llike;
dnadata *dna;
msatdata *ms;

dna = data->dnaptr;
ms = data->msptr;

subtree_ranges = (long *)calloc(2*curtree->numrecombs+4,sizeof(long));
subtree_ranges[0] = 1;

siteptr = data->siteptr;

findsubtrees(op,data,curtree->tymelist,subtree_ranges);

llike = 0.0;
for(subtree = 0; subtree < subtree_ranges[0]; subtree++) {
   substart = subtree_ranges[2*subtree+1];
   subend = subtree_ranges[2*subtree+2];
   /* do the regular sites */
   localsmooth(op,data,p,substart,substart,subend);
   switch(op->datatype) {
      case 'a':
         break;
      case 'b':
      case 'm':
/* WARNING this throws away llike, which must be wrong! */
         llike += micro_evaluate(op,ms,curtree,substart,subend);
         break;
      case 'n':
         lastmarker = getdata_nummarkers(op,data)-1;
         snpsmooth(op,data,curtree,substart,FLAGLONG,FLAGLONG);
         if (op->panel)
            llike += panel_eval_calcrange(op,data,curtree,first,substart,
               subend,op->categs);
         else
            llike += snp_eval_calcrange(op,data,curtree,first,substart,
               subend,op->categs);
         dummy = snp_evaluate(op,data,curtree,llike,substart,subend);
         break;
      case 's':
         llike += eval_calcrange(op,data,curtree,first,substart,subend,
            op->categs);
         dummy = evaluate(op,data,curtree,llike);
         break;
      default:
         fprintf(ERRFILE,"\nlocaleval--impossible datatype %c\n",
            op->datatype);
         break;
   }
}

traverse_unflag(curtree,curtree->nodep[1]);
free(subtree_ranges);

} /* localeval */


/********************************************************************
 * micro_evaluate() calculates the tree data likelihood at the root *
 * assuming that the tree has been properly "smoothed" first.       */
double micro_evaluate(option_struct *op, msatdata *ms, tree *tr,
  long start, long end)
{
long site, a;
double term;
node *nn;

nn = tr->root->back;
term = 0.0;

for(site = start; site <= end; site++)
   for(a = 0; a < MICRO_ALLELEMAX-1; a++)
      term += nn->x->a[site][a];

if(term == 0.0) return(NEGMAX);
else return(log(term) + nn->lxmax);

} /* micro_evaluate */


/******************************************************************
 * countbranches returns the total number of branches in the tree *
 *    not including the root.                                     */
long countbranches(option_struct *op, data_fmt *data, tree *tr)
{

return(getdata_numtips(op,data) + 2 * tr->numrecombs + tr->numcoals - 1);

} /* countbranches */

/******************************************************************
 * testratio() returns TRUE to accept a tree, FALSE to reject it. *
 * ratiotype = "d" if called to decide a "drop"                   *
 *             "t" if called to decide a "twiddle"                *
 *             "f" if called to decide a "flip"                   */
boolean testratio(option_struct *op, data_fmt *data, tree *oldtree,
   tree *newtree, char ratiotype)
{
double test, x, numoldbranches, numnewbranches;

if(newtree->likelihood == NEGMAX)
   return FALSE;
#if 0 /* first heating */
test = (newtree->likelihood - oldtree->likelihood) * (1.0/op->temperature[0]);
#endif
test = (newtree->likelihood - oldtree->likelihood) * (1.0/op->ctemp);
switch((int)ratiotype) {
case 'd' :
   numoldbranches = (double)countbranches(op,data,oldtree);
   numnewbranches = (double)countbranches(op,data,newtree);
   test += log(numoldbranches/numnewbranches);
   break;
case 't' :
   test += log((double)oldtree->numrecombs/(double)newtree->numrecombs);
   break;
case 'f' :
   break;
default :
   fprintf(ERRFILE,"ERROR:testratio encounters unknown type--no Hastings\n");
  break;
}

#if 0 /* first heating */
if (op->temperature[0] != 1) { /* heating code */
  test += (newtree->coalprob - oldtree->coalprob) *
             (1.0/op->temperature[0] - 1.0);
}
#endif

#if 0 /* don't heat this */
test += (newtree->coalprob - oldtree->coalprob) *
             (1.0/op->ctemp - 1.0);
#endif

if (test >= 0.0)
   return TRUE;
else {
   x = log(randum());
   if (test >= x)
      return TRUE;
   else
      return FALSE;
}
} /* testratio */


void seekch(char c) /* use only in reading file intree! */
{
  if (gch == c)
    return;
  do {
    if (eoln(intree)) {
      fscanf(intree, "%*[^\n]");
      getc(intree);
    }
    gch = getc(intree);
    if (gch == '\n')
      gch = ' ';
  } while (gch != c);
}  /* seekch */

void getch(char *c) /* use only in reading file intree! */
{
  /* get next nonblank character */
  do {
    if (eoln(intree)) {
      fscanf(intree, "%*[^\n]");
      getc(intree);
    }
    *c = getc(intree);
    if (*c == '\n')
      *c = ' ';
  } while (*c == ' ');
}  /* getch */

void processlength(node *p)
{
  long digit;
  double valyew, divisor;
  boolean pointread;

  pointread = FALSE;
  valyew = 0.0;
  divisor = 1.0;
  getch(&gch);
  digit = gch - '0';
  while (((unsigned long)digit <= 9) || gch == '.'){
    if (gch == '.')
      pointread = TRUE;
    else {
      valyew = valyew * 10.0 + digit;
      if (pointread)
	divisor *= 10.0;
    }
    getch(&gch);
    digit = gch - '0';
  }
  p->length = valyew / divisor;
  p->back->length = p->length;
}  /* processlength */

void addelement(option_struct *op, data_fmt *data, node *p, long *nextnode)
{
  node *q;
  long i, n;
  boolean found;
  char str[NMLNGTH];

  getch(&gch);
  if (gch == '(') {
    (*nextnode)++;
    newnode(&q);
    q = curtree->nodep[(*nextnode)];
    hookup(p, q);
    addelement(op,data,q->next,nextnode);
    seekch(',');
    addelement(op,data,q->next->next, nextnode);
    seekch(')');
    getch(&gch);
  } else {
    for (i = 0; i < NMLNGTH; i++)
      str[i] = ' ';
    n = 1;
    do {
      if (gch == '_')
	gch = ' ';
      str[n - 1] = gch;
      if (eoln(intree)) {
	fscanf(intree, "%*[^\n]");
	getc(intree);
      }
      gch = getc(intree);
      if (gch == '\n')
	gch = ' ';
      n++;
    } while (gch != ':' && gch != ',' && gch != ')' && n <= NMLNGTH);
    n = 1;
    do {
      found = TRUE;
      for (i = 0; i < NMLNGTH; i++)
	found = (found && str[i] == curtree->nodep[n]->nayme[i]);
      if (!found)
	n++;
    } while (!(n > getdata_numseq(op,data) || found));
    if (n > getdata_numseq(op,data)) {
      printf("Cannot find sequence: ");
      for (i = 0; i < NMLNGTH; i++)
	putchar(str[i]);
      putchar('\n');
    }
    hookup(curtree->nodep[n], p);
  }
  if (gch == ':')
    processlength(p);
}  /* addelement */

void treeread(option_struct *op, data_fmt *data)
/* WARNING:  this only works with sequence data, not
  microsatellites, not panel SNPs, not allozymes! */
{
  long nextnode;
  node *p;

  curtree->root = curtree->nodep[rootnum];
  getch(&gch);
  if (gch == '(') {
    nextnode = getdata_numtips(op,data) + 1;
    p = curtree->nodep[nextnode];
    addelement(op,data,p, &nextnode);
    seekch(',');
    addelement(op,data,p->next, &nextnode);
    hookup(p->next->next, curtree->nodep[rootnum]);
    p->next->next->length = ROOTLENGTH;
    curtree->nodep[rootnum]->length = p->next->next->length;
    ltov(op,data,curtree->nodep[rootnum]);
  }
  fscanf(intree, "%*[^\n]");
  getc(intree);
}  /* treeread */


/* DEBUG debug WARNING warning --
   actually still used in traitlike.c:traitlike() */
void finddnasubtrees(option_struct *op, data_fmt *data,
   tlist *tstart, long *sranges)
{
long i, j, numsites, *temp;
tlist *t;

numsites = countsites(op,data);

temp = (long *)calloc(numsites,sizeof(long));

for(t = tstart, i = 1; t != NULL; t = t->succ)
   if(isrecomb(t->eventnode))
      temp[findlink(t->eventnode)] = 1;

for(i = 0, j = 2, sranges[1] = 0; i < numsites; i++) {
   if(temp[i] == 0) continue;
   sranges[0]++;
   sranges[j] = i;
   sranges[j+1] = i+1;
   j+=2;
}
sranges[j] = numsites-1;
sranges[j+1] = FLAGLONG;

if (j+1 > 2*sranges[0]+3)
   fprintf(ERRFILE,"ERROR:finddnasubtree: j calculated wrong\n");

free(temp);

} /* finddnasubtrees */


void findsubtrees_node(option_struct *op, data_fmt *data,
   tlist *tlast, tree *tr, long **coal)
{
long i, j, numsites, *temp;
tlist *t;

numsites = countsites(op,data);

temp = (long *)calloc(numsites,sizeof(long));
for(t = tr->tymelist, i = 1; t != tlast->succ; t = t->succ)
   if(isrecomb(t->eventnode)) {
      j = findlink(t->eventnode);
      if (temp[j]) continue;
      temp[j] = 1;
      i++;
   }

init_coal_alloc(coal,i);

for(i = 0, j = 2; i < numsites; i++) {
   if (temp[i] == 0) continue;
   (*coal)[j] = i;
   (*coal)[j+1] = i+1;
   j += 2;
}
(*coal)[j] = numsites-1;

free(temp);

} /* findsubtrees_node */


void drop_findsubtrees(option_struct *op, data_fmt *data, tree *tr,
   node *p)
{
node *q;
tlist *t;

t = gettymenode(tr,p->number);

findsubtrees_node(op,data,t,tr,&(p->coal));

if (isrecomb(p)) {
   q = otherdtr(p);
   findsubtrees_node(op,data,t,tr,&(q->coal));
}

} /* drop_findsubtrees */


void findsubtrees(option_struct *op, data_fmt *data, tlist *tstart,
   long *sranges)
{
long i, j, numsites, *temp;
tlist *t;

numsites = countsites(op,data);

temp = (long *)calloc(numsites,sizeof(long));

for(t = tstart; t != NULL; t = t->succ)
   if(isrecomb(t->eventnode))
      temp[findlink(t->eventnode)] = 1;

for(i = 0, j = 2, sranges[1] = 0; i < numsites; i++) {
   if(temp[i] == 0) continue;
   sranges[0]++;
   sranges[j] = i;
   sranges[j+1] = i+1;
   j+=2;
}
sranges[j] = numsites-1;
sranges[j+1] = FLAGLONG;

if (j+1 > 2*sranges[0]+3)
   fprintf(ERRFILE,"ERROR:findsubtree: j calculated wrong\n");

free(temp);

} /* findsubtrees */


/***************************************************************
 * markertosite() returns the exact "site" that corresponds to *
 * the passed marker position.                                 */
long markertosite(option_struct *op, data_fmt *data, long marker)
{

return(getdata_markersite(op,data,marker));

} /* markertosite */


/**********************************************************************
 * sitetorightmarker() returns the closest marker to the right of the *
 * given site, if that marker exists.  If that marker doesn't exist,  *
 * FLAGLONG is returned.  If the site exactly corresponds to a marker *
 * then the corresponding marker will be returned.                    */
long sitetorightmarker(option_struct *op, data_fmt *data, long site)
{
long marker, nummarkers;
double msite;

nummarkers = getdata_nummarkers(op,data);
for(marker = 0; marker < nummarkers; marker++) {
   msite = markertosite(op,data,marker);
   if (msite >= site) return(marker);
}

return(FLAGLONG);

} /* sitetorightmarker */


/**********************************************************************
 * sitetomarker() returns the closest marker to the left of the given *
 * site, if that marker exists.  If that marker doesn't exist, the    *
 * closest marker on the right is returned.                           */
long sitetomarker(option_struct *op, data_fmt *data, long site)
{
long marker, nummarkers;
double sitecounts;

nummarkers = getdata_nummarkers(op,data);
for(marker = 0, sitecounts = 0.0; marker < nummarkers; marker++) {
   sitecounts = data->dnaptr->sitecount[locus][marker];
   if (site < sitecounts) return(marker);
}

fprintf(ERRFILE,"ERROR--failure to find marker for site %ld\n\n",
   site);
return(FLAGLONG);


} /* sitetomarker */


/**********************************************************
 * findsubtreemarkers() gets passed the beginning and end *
 * of a subtree (in psuedosites) and returns the position *
 * of the first and last markers found in that subtree.   *
 * Return FLAGLONG in both marker positions if there are  *
 * no markers present in the subtree.                     */
void findsubtreemarkers(option_struct *op, data_fmt *data, long pstart,
  long pend, long *mstart, long *mend)
{
long nummarkers, marker, site;

nummarkers = getdata_nummarkers(op,data);

for(marker = 0; marker < nummarkers; marker++) {
   site = markertosite(op,data,marker);
   if (site >=  pstart) break;
}

if (site > pend) {
    (*mstart) = (*mend) = FLAGLONG;
    return;
}

(*mstart) = marker;

for(;marker < nummarkers; marker++) {
   site = markertosite(op,data,marker);
   if (site > pend) break;
}

(*mend) = marker-1;

} /* findsubtreemarkers */


/**********************************************************
 * sameranges() returns TRUE if the ranges are identical, *
 * FALSE otherwise.                                       */
boolean sameranges(long *range1, long *range2)
{
long i;

for(i = 0; range1[i] != FLAGLONG && range2[i] != FLAGLONG; i++)
   if(range1[i] != range2[i]) return(FALSE);

if (range1[i] != range2[i]) {
   printf("ERROR:sameranges bailed on FLAGLONG but not done yet!!! %ld %ld\n",
      indecks,apps);
   return(FALSE);
}

return(TRUE);

} /* sameranges */


/*************************************************
 * copycoal() copies the coal array of a nodelet */
void copycoal(node *source, node *target)
{

if (source->coal != NULL) {
   coal_Malloc(target,TRUE,source->coal[0]);
   memcpy(target->coal,source->coal,
      (source->coal[0]*2+2)*sizeof(long));
} else coal_Malloc(target,FALSE,0L);

} /* copycoal */


/*****************************************************
 * copyranges() copies the ranges array of a nodelet */
void copyranges(node *source, node *target)
{

if (source->ranges != NULL) {
   ranges_Malloc(target,TRUE,source->ranges[0]);
   memcpy(target->ranges,source->ranges,
      (source->ranges[0]*2+2)*sizeof(long));
} else ranges_Malloc(target,FALSE,0L);

} /* copyranges */


/***************************************************************
 * This function copies the likelihood, "x" and "z", arrays of *
 * a nodelet.                                                  */
void copylikes(option_struct *op, data_fmt *data, node *source,
   node *target)
{
long nummarkers = 0, numslice;

numslice = (op->panel) ? NUMSLICE : 1L;

if (op->map) {
   if (source->z != NULL) {
      allocate_z(op,data,target);
      memcpy(target->z,source->z,NUMTRAIT*sizeof(double));
   } else free_z(op,target);
}

if (source->x != NULL) {
   allocate_x(op,data,target);
   switch (op->datatype) {
      case 'a':
         break;
      case 'b':
      case 'm':
         memcpy(target->x->a[0],source->x->a[0],
         getdata_numloci(op,data)*MICRO_ALLELEMAX*sizeof(double));
         break;
      case 'n':
         nummarkers += NUMINVAR;
      case 's':
         nummarkers += getdata_nummarkers(op,data);
         memcpy(target->x->s[0][0][0],source->x->s[0][0][0],
            nummarkers*op->categs*numslice*4L*sizeof(double));
         break;
      default:
         fprintf(ERRFILE,"\ncopylikes--impossible datatype %c\n",
            op->datatype);
         exit(-1);
         break;
   }
} else free_x(op,target);

} /* copylikes */


/************************************************************
 * copynode copies the "source" node onto the "target" node */
void copynode(option_struct *op, data_fmt *data, node *source,
   node *target)
{
  long i, nodecount;

  if (istip(source)) nodecount = 1;
  else nodecount = 3;

  for (i = 0; i < nodecount; i++) {
    target->type = source->type;
    /* NEVER! target->next := source->next; */
    target->back = source->back;
    target->top = source->top;
    /* but NOT target->number := source->number; */
    copylikes(op,data,source,target);
    copyranges(source,target);
    if (op->fc) copycoal(source,target);
    if (istip(source)) memcpy(target->nayme,source->nayme,sizeof(naym));
    target->v = source->v;
    target->lxmax = source->lxmax;
    target->tyme = source->tyme;
    target->length = source->length;
    target->updated = source->updated;
    target->recstart = source->recstart;
    target->recend = source->recend;
    source = source->next;
    target = target->next;
  }
}  /* copynode */

/**********************************************************************
 * addnode hooks nodelet "p" to node "q" by free back ptrs in p and q */
void addnode(option_struct *op, data_fmt *data, node *p, node *q)
{

if (p->top) while(q->top && q->back) q = q->next;
else while(!q->top && q->back) q = q->next;

hookup(p,q);
fixlength(op,data,p);

} /* addnode */


/****************************************************************
 * getrecnodelet finds which nodelet in node "copy" corresponds *
 * to the exact nodelet pointed to by "orig".                   */
node *getrecnodelet (node *copy, node *orig)
{
long i;
node *p;

p = copy;
for(i = 0; i < 3; i++) {
   if ((p->recstart == orig->recstart) && (p->recend == orig->recend))
      return(p);
   p = p->next;
}

fprintf(ERRFILE,"ERROR:getrecnodelet failed to match\n");
return(NULL);

} /* getrecnodelet */


/**********************************************************
 * make_tree_copy copies tree "source" into tree "target" */
void make_tree_copy(option_struct *op, data_fmt *data, tree *source,
   tree *target)
{
long i, nodenumber, *tnum, numno, numtips;
tlist *t, *addhere;
node *p, *q, *r;
dnadata *dna;

dna = data->dnaptr;
numtips = getdata_numtips(op,data);
numno = numtips + source->numcoals + source->numrecombs + 1;
tnum = (long *)calloc(1,numno*sizeof(long));

/* make the tips */
target->nodep = (node **)calloc(1,numno*sizeof(node *));
for (i = 0; i < numtips+1; i++) {
   allocate_tip(op,data,target,source->nodep[i]->number);
   strcpy(target->nodep[i]->nayme,source->nodep[i]->nayme);
   if (i) {
      target->nodep[i]->top = TRUE;
      target->nodep[i]->tyme = 0.0;
   } else {
      target->nodep[i]->top = FALSE;
      target->nodep[i]->tyme = ROOTLENGTH;
   }
   target->nodep[i]->updated = TRUE;
   copyranges(source->nodep[i],target->nodep[i]);
   if(op->map) copytraits(source->nodep[i],target->nodep[i]);
   if(op->fc) copycoal(source->nodep[i],target->nodep[i]);
   copylikes(op,data,source->nodep[i],target->nodep[i]);
   tnum[i] = i;
}

/* now construct the first tymelist entry */
newtymenode(&target->tymelist);
target->tymelist->numbranch = numtips;
target->tymelist->branchlist =
   (node **)calloc(numtips,sizeof(node *));
for(i = 0; i < numtips; i++)
   target->tymelist->branchlist[i] = target->nodep[i+1];
target->tymelist->eventnode = target->nodep[1];
target->tymelist->age = source->root->tyme;

/* now do the rest of the tree and tymelist */
for (t = source->tymelist->succ; t != NULL; t = t->succ) {
   nodenumber = t->eventnode->number;
   newnode(&p);
   p->number = nodenumber;
   p->next->number = nodenumber;
   p->next->next->number = nodenumber;
   target->nodep[nodenumber] = p;
   copynode(op,data,t->eventnode,p);
   p->back = NULL;
   p->next->back = NULL;
   p->next->next->back = NULL;
   tnum[t->eventnode->number] = nodenumber;
}

addhere = target->tymelist;
for (t = source->tymelist->succ; t != NULL; t = t->succ) {
   p = findunique(target->nodep[tnum[t->eventnode->number]]);
   q = findunique(t->eventnode);
   if (isrecomb(p)) {
      r = target->nodep[tnum[q->back->number]];
      if (isrecomb(r)) r = getrecnodelet(r,q->back);
      addnode(op,data,p,r);
   } else {
      r = target->nodep[tnum[q->next->back->number]];
      if (isrecomb(r)) r = getrecnodelet(r,q->next->back);
      addnode(op,data,p->next,r);
      r = target->nodep[tnum[q->next->next->back->number]];
      if (isrecomb(r)) r = getrecnodelet(r,q->next->next->back);
      addnode(op,data,p->next->next,r);
   }

   insertaftertymelist(addhere,p);
   addhere = addhere->succ;
}

addnode(op,data,target->nodep[0],
   target->nodep[tnum[source->root->back->number]]);
target->root = target->nodep[0];

free(tnum);

} /* make_tree_copy */


/********************************************************************
 * copycreature() copies a single creature to the "target" pointer. *
 * The passed tree should be the tree associated with the "target"  *
 * creature.                                                        */
void copycreature(option_struct *op, data_fmt *data, creature *source,
   creature *target, tree *tr)
{
long i;

target->numflipsites = source->numflipsites;
target->flipsites = (long *)calloc(target->numflipsites,sizeof(long));
memcpy(target->flipsites,source->flipsites,
   (target->numflipsites)*sizeof(long));

target->numhaplotypes = source->numhaplotypes;
target->haplotypes = (node **)calloc(target->numhaplotypes,sizeof(node *));
for(i = 0; i < target->numhaplotypes; i++) {
   target->haplotypes[i] = tr->nodep[source->haplotypes[i]->number];
}

} /* copycreature */


/********************************************************
 * copycreatures() copies the creatures array of a tree */
void copycreatures(option_struct *op, data_fmt *data, tree *source,
   tree *target)
{
long cr, numcreatures;

if (!op->haplotyping) return;

numcreatures = getdata_numtips(op,data)/NUMHAPLOTYPES;
target->creatures = (creature *)calloc(numcreatures,sizeof(creature));

for(cr = 0; cr < numcreatures; cr++) {
   copycreature(op,data,&(source->creatures[cr]),
      &(target->creatures[cr]),target);
}

} /* copycreatures */


/**********************************************************************
 * copytree makes a copy from scratch of the "source" tree, returning *
 * a pointer to the copy that it makes.                               */
tree *copytree(option_struct *op, data_fmt *data, tree *source)
{
tree *target;

target = (tree *)calloc(1,sizeof(tree));

target->likelihood = source->likelihood;
target->coalprob = source->coalprob;
target->numcoals = source->numcoals;
target->numrecombs = source->numrecombs;
make_tree_copy(op,data,source,target);
copycreatures(op,data,source,target);

return(target);
} /* copytree */


/***************************************************
 * freetree frees the "target" tree from memory */
void freetree(option_struct *op, data_fmt *data, tree *target)
{
long i;
tlist *t;

if (op->haplotyping) {
   for(i = 0; i < getdata_numtips(op,data)/NUMHAPLOTYPES; i++) {
      free(target->creatures[i].haplotypes);
      free(target->creatures[i].flipsites);
   }
   free(target->creatures);
}

t = target->tymelist->succ;
while(t != NULL) {
   freenode(t->eventnode);
   t = t->succ;
}
freetymelist(target->tymelist);

for (i = 0; i < getdata_numtips(op,data)+1; i++) {
   freenodelet(op,data,target->nodep[i]);
}
free(target->nodep);

free(target);

} /* freetree */

/* joinnode and constructtree are used for constructing a rather bad
   starting tree if the user doesn't provide one */
void joinnode(option_struct *op, data_fmt *data, double length,
   node *p, node *q)
{
   hookup(p,q);
   p->length = length;
   q->length = length;
   ltov(op,data,p);
} /* joinnode */

void constructtree(option_struct *op, data_fmt *data, double branchlength)
{
long i, j, numtips, nextnode;
double height;
node *p, *q;

numtips = getdata_numtips(op,data);

curtree->root = curtree->nodep[rootnum];
nextnode = numtips+1;
 p = curtree->root;
 q = curtree->nodep[nextnode];

p->back = q;
q->back = p;
p->length = ROOTLENGTH;
q->length = ROOTLENGTH;
ltov(op,data,p);

height = (numtips - 1) * branchlength;
p->tyme = ROOTLENGTH + height;
for (i = 1; i < numtips; i++) {
   p = curtree->nodep[i];
   q = curtree->nodep[nextnode]->next;
   joinnode(op,data,height,p,q);
   q = q->next;
   if (i != numtips-1) {
      nextnode++;
      p = curtree->nodep[nextnode];
      joinnode(op,data,branchlength,p,q);
      height -= branchlength;
   } else {
      p = curtree->nodep[numtips];
      joinnode(op,data,height,p,q);
   }
   for (j = 0; j < 3; j++)
      q->tyme = height;
}

} /* constructtree */
/* End bad starting tree construction */


/******************************************************************
 * min_theta_calc() calculates the minimum value which we wish to *
 * let the mutation-rate parameter, theta, attain before final    *
 * estimation.                                                    */
double min_theta_calc(option_struct *op, data_fmt *data, double th)
{
double value;

value = THETAMIN;

return(value);

} /* min_theta_calc */


/***********************************************************
 * hasrec_chain() returns TRUE if the passed chain has any *
 * recombinations currently sampled, FALSE otherwise.      */
boolean hasrec_chain(option_struct *op, long chain)
{
long i, chaintype, refchain;
treerec *trii;

refchain = REF_CHAIN(chain);
chaintype = TYPE_CHAIN(chain);

for(i = 0; i < op->numout[chaintype]; i++) {
   trii = &sum[locus][refchain][i];
   if (trii->numrecombs > 0) return(TRUE);
}

return(FALSE);

} /* hasrec_chain */


/***********************************************************************
 * chainendcheck() runs at the end of a chain and makes sure the chain *
 * had the following desired properties:                               *
 *                                                                     *
 *    --at least 1 sampled tree contains 1+ recombinations             *
 *      FAILURE = call addfractrecomb()                                */
void chainendcheck(option_struct *op, data_fmt *data, tree *tr,
   long chain, boolean locusend)
{

if(!locusend && !hasrec_chain(op,chain) && op->holding != 2) {
   if (apps == totchains - 1) {
      fprintf(outfile,"WARNING--forced a recombination");
      fprintf(outfile," which may have affected final result\n");
   }
   addfractrecomb(op,data,chain,sum);
}

} /* chainendcheck */


/******************************************************************
 * rearrange() is the driver for basic rearrangement of a tree,   *
 * returning TRUE if it successfully rearranged a tree, and FALSE *
 * otherwise.                                                     */
boolean rearrange(option_struct *op, data_fmt *data, long whichtchain)
{
double chance;
boolean accepted;

curtree = temptrees[whichtchain];
chance = randum();
accepted = FALSE;
op->ctemp = op->temperature[whichtchain];

if (chance < TWIDDLE_PROB) accepted = twiddle(op,data);
else {
   chance -= TWIDDLE_PROB;
   if (chance < op->happrob && op->haplotyping) {
      if (op->hapdrop == 0)
         accepted = fliphap(op,data,curtree);
      else if (op->hapdrop == 1)
         accepted = flipdrop(op,data,curtree,1L);
      else
         accepted = flipdrop(op,data,curtree,2L);
      hap++;
      if (accepted && whichtchain == 0) hacc++;
   } else {
      accepted = makedrop(op,data);
      slid++;
      if (accepted && whichtchain == 0) slacc++;
   }
}
#ifdef MAC
   eventloop();
#endif

temptrees[whichtchain] = curtree;

return(accepted);

} /* rearrange */


/**********************************************************************
 * temptreeswap() checks to see if a swap of trees between different  *
 * temperatures is called for, executing the swap if so.              *
 *                                                                    *
 * it accepts a proposed swap using                                   *
 *    [ P(Gh|Tl) * P(D|Gh) ] ^ 1/l * [ P(Gl|Th) * P(D|Gl) ] ^ 1/h     *
 *    -----------------------------------------------------------     *
 *    [ P(Gh|Th) * P(D|Gh) ] ^ 1/h * [ P(Gl|Tl) * P(D|Gl) ] ^ 1/l     *
 *                                                                    *
 * where Gh = the tree generated at high temp                         *
 *       Gl = the tree generated at low temp                          *
 *       Tl = the parameter values at the low temp                    *
 *       Th = the parameter values at the high temp                   */
void temptreeswap(option_struct *op, data_fmt *data, boolean *changed,
   long chain)
{
long chl, chh;
double chance, num, denom, templ, temph, recl, rech, thetal, thetah;
boolean picked;
tree *trh, *trl;

if (op->numtempchains < 2) return;

/* first pick which two chains to swap */

picked = FALSE;
while(!picked) {
   chl = (long)(randum()*op->numtempchains);
   chh = (long)(randum()*op->numtempchains);
   if (chl == op->numtempchains || chh == op->numtempchains) continue;
   if (chh - chl > 0) {
      if (chh - chl == 1) picked = TRUE;
   } else {
      if (chl - chh == 1) picked = TRUE;
   }
}

trh = temptrees[chh];
temph = op->temperature[chh];
thetah = theti[locus][chain];
rech = reci[locus][chain];

trl = temptrees[chl];
templ = op->temperature[chl];
thetal = theti[locus][chain];
recl = reci[locus][chain];

num = (coalprob(op,data,trh,thetal,recl)+trh->likelihood)*(1.0/templ) +
   (coalprob(op,data,trl,thetah,rech)+trl->likelihood)*(1.0/temph);

denom = (trh->coalprob+trh->likelihood)*(1.0/temph) +
     (trl->coalprob+trl->likelihood)*(1.0/templ);

chance = num - denom;

swap++;
if (chance >= log(randum())) {
   temptrees[chh] = trl;
   temptrees[chl] = trh;
   changed[chh] = TRUE;
   changed[chl] = TRUE;
   swacc++;
}

} /* temptreeswap */


#define PRINTNUM     8 /* the number of digits used for printing
                          the number of steps done */

void maketree(option_struct *op, data_fmt *data)
{
long tempchain, incrprog, metout, progout, chaintype, i;
double bestlike;
boolean runend, *changetree, first;
char chainlit[2][6] = {"Short","Long"};
dnadata *dna;
double *traitarray;
long hapscore;

traitarray = NULL;

changetree = (boolean *)calloc(op->numtempchains,sizeof(boolean));

if (op->map)
  traitarray = (double *)calloc(countsites(op,data),sizeof(double));

chainlit[0][5] = chainlit[1][5] = '\0';

/* WARNING DEBUG assumes DNA data (and shouldn't)! */
dna = data->dnaptr;

  getc(infile);
if(op->haplotyping) {
  fprintf(outfile,"Haplotypes being inferred");
  if (op->hapdrop == 0) fprintf(outfile," without resimulation\n");
  if (op->hapdrop == 1) fprintf(outfile," with single resimulation\n");
  if (op->hapdrop == 2) fprintf(outfile," with double resimulation\n");
}
  fprintf(outfile,"Watterson estimate of theta is %12.8f\n", watttheta);
#if !DISTRIBUTION
  fprintf(outfile,"\n\nIntermediate chains (only)\n");
  fprintf(outfile,"chain#     theta     rec-rate\n");
  fprintf(outfile,"-----------------------------\n");
#endif
  if (op->usertree)
     treeread(op,data);
  else {
#if !ALWAYS_REJECT
     branch0 = watttheta/getdata_numtips(op,data);
     constructtree(op,data,branch0);
#else
     constructtree(op,data,0.0);
#endif
  }
  orient(curtree,op,data,curtree->root->back);
#if !ALWAYS_REJECT
  if (op->plump) plumptree(op,data,watttheta);
#endif
  finishsetup(op,data,curtree->root->back);
  initbranchlist(op,data);
#if !ALWAYS_REJECT
  localeval(op,data,curtree->root->back,TRUE);
  curtree->coalprob = coalprob(op,data,curtree,theta0,rec0);
#endif
  bestlike = NEGMAX;
  theti[locus][0] = theta0;
  reci[locus][0] = rec0;
  runend = FALSE;

/* WARNING DEBUG MARY */
    if (TRUEHAP) {
      hapscore = hapdist(op,&truehaps,data,curtree->creatures);
      printf("Distance from truth%ld\n",hapscore);
      fprintf(outfile,"Distance from truth%ld\n",hapscore);
      hapscore = hapdist(op,&starthaps,data,curtree->creatures);
      fprintf(outfile,"Distance from start%ld\n",hapscore);
    }

  if(op->map) traitread(curtree, getdata_numseq(op,data));
  /**********************************/
  /* Begin Hastings-Metropolis loop */
  /**********************************/
  for (apps = 0; apps < totchains; apps++) {
    if (apps >= op->numchains[0]) chaintype = 1;
    else chaintype = 0;
    if (op->progress) {
      printf("%s chain %ld ",chainlit[chaintype],
      ((chaintype) ? (apps + 1 - op->numchains[0]) : apps + 1));
      fflush(stdout);
    }
    numdropped = 0;
    metout = op->increm[chaintype] - 1 + NUMBURNIN;
  /* print a "." to stdout every 10% of the chain(s) finished */
    incrprog = (long)(op->steps[chaintype] / 10.0);
    progout = incrprog - 1 + NUMBURNIN;
    op->numout[chaintype] = 0;
    slacc = hacc = swacc = 0;
    slid = hap = swap = 0;

/* if apps == 0, start all temperature chains with the same
   initial tree */
    if (apps == 0) {
       for(tempchain = 0; tempchain < op->numtempchains; tempchain++) {
          temptrees[tempchain] = copytree(op,data,curtree);
       }
       /* curtree is no longer needed as a base tree!, it's just a ptr */
       freetree(op,data,curtree);
    }
    for(tempchain = 0; tempchain < op->numtempchains; tempchain++)
       changetree[tempchain] = FALSE;
    first = TRUE;
    for (indecks=0; indecks < op->steps[chaintype]+NUMBURNIN; indecks++) {
#if 0 /* first heating */
/* heating code */
      if (max_heat != 1.0) {
         if (heating_cnt == heating_interval) {
            heating_cnt = 0;
            if (op->temperature[0] >= max_heat) direction = -1;
            if (op->temperature[0] <= 1.0) direction = 1;
            op->temperature[0] += direction * heating_by;
         }
         heating_cnt++;
      }
#endif
      col = 0; /* column number, used in treeout */

/* only count acceptances after burn-in */
      if (indecks == NUMBURNIN)
         hacc = slacc = 0;

      for(tempchain = 0; tempchain < op->numtempchains; tempchain++) {
         if (rearrange(op,data,tempchain))
            changetree[tempchain] = TRUE;
      }

      temptreeswap(op,data,changetree,apps);


      if (apps == totchains - 1) /* end of run? */
        runend = TRUE;
      if (temptrees[0]->likelihood > bestlike) {
        if (ONEBESTREE) {
           FClose(bestree);
           bestree = fopen("bestree","w+");
	   fprintf(bestree, "Chain #%2ld (%s) Step:%8ld\n",apps+1,
             chainlit[chaintype], indecks+1);
/* debug DEBUG warning WARNING--no rectreeout routine */
           /*rec_outtree(temptrees[0]->root->back,TRUE, &bestree);*/
	   fprintf(bestree, "; [%12.10f]\n", temptrees[0]->likelihood);
	   bestlike = temptrees[0]->likelihood;
        }
        else {
	   fprintf(bestree, "Chain #%2ld (%s) Step:%8ld\n",apps+1,
             chainlit[chaintype], indecks+1);
/* debug DEBUG warning WARNING--no rectreeout routine */
           /*rec_outtree(temptrees[0]->root->back,TRUE, &bestree); */
	   fprintf(bestree, "; [%12.10f]\n", temptrees[0]->likelihood);
	   bestlike = temptrees[0]->likelihood;
        }
      }
      if (indecks == metout) {
        if (op->numout[chaintype] == 0)
          sametree[locus][0] = FALSE;
        else
          sametree[locus][op->numout[chaintype]] = !changetree[0];
        changetree[0] = FALSE;
#if !ALWAYS_REJECT
#if 0 /* old heating */
        if (op->temperature[0] == 1.0) { /* not too heated */
	  op->numout[chaintype]++;
       	  scoretree(op,data,apps);
          if(op->map)
            traitlike(op,data,curtree,countsites(op,data),
              op->mutrait,op->traitratio, op->pd, traitarray);
        }
#endif
	op->numout[chaintype]++;
/* set curtree to the coldest tree for scoretree() scoring */
        curtree = temptrees[0];
       	scoretree(op,data,apps);
        if(op->map)
          traitlike(op,data,curtree,countsites(op,data),
            op->mutrait,op->traitratio, op->pd, traitarray);
#endif
	metout += op->increm[chaintype];
        if (op->treeprint) {
           fprintf(treefile,"\nlocus = %ld, chain = %ld, tree = %ld\n",
                   locus,apps,indecks);
/* debug DEBUG warning WARNING--no rectreeout routine */
           /*rec_outtree(temptrees[0]->root->back,TRUE, &treefile);*/
           fprintf(treefile,";\n");
        }
      }
#if 0
      if (op->progress && indecks == progout) {
	printf(".");
        fflush(stdout);
	progout += incrprog;
      }
#endif
      if (op->progress) {
         if (!first) for(i = 0; i < PRINTNUM; i++) printf("\b");
         else printf(" examined: ");
         printf("%*ld",(int)PRINTNUM,indecks-NUMBURNIN+1);
         fflush(stdout);
         first = FALSE;
      }
    }
    if(runend) {
       if(op->map) {
          traitprint(countsites(op,data),traitarray,
            outfile,op->numout[chaintype]);
#if TRUTHKNOWN
          traitresult(op,data,traitarray,outfile,op->numout[chaintype]);
#endif
       }
    }
    printf(" trees");
    chainendcheck(op,data,curtree,apps,runend);
    rec_estimate(op,data,locus,apps,runend);
    if (TRUEHAP) {
      hapscore = hapdist(op,&truehaps,data,curtree->creatures);
      printf("Distance from truth: %ld\n",hapscore);
      fprintf(outfile,"Distance from truth: %ld\n",hapscore);
      hapscore = hapdist(op,&starthaps,data,curtree->creatures);
      printf("Distance from start: %ld\n",hapscore);
      fprintf(outfile,"Distance from start: %ld\n",hapscore);
    }

    if (op->map && runend) traitsiteplot(op,data,traitarray);
    theta0 = theti[locus][apps+1];
    rec0 = reci[locus][apps+1];
    if (min_theta_calc(op,data,theta0) > theta0 &&
       op->holding != 1) {
       if (apps == totchains - 2) {
          fprintf(outfile,"WARNING--lower bound on estimate of");
          fprintf(outfile," theta may have affected final result\n");
       }
       theta0 = min_theta_calc(op,data,theta0);
       theti[locus][apps+1] = theta0;
    }
    if (op->progress) {
       printf("accepted %ld/%ld trees",slacc,slid-NUMBURNIN);
       if (numdropped) printf(", dropped %ld",numdropped);
       if (swap) printf(", swapped %ld/%ld",swacc,swap);
       if (op->haplotyping)
          printf("and %ld/%ld haplotype changes\n",hacc,hap);
       else printf("\n");
    }
    if ((apps == totchains-1) && numdropped) {
       fprintf(outfile,"%ld trees were dropped from",numdropped);
       fprintf(outfile," the final chain\n");
    }
  }
  if(slacc == 0) {
     fprintf(outfile,"WARNING--no proposed trees ever accepted\n");
     fprintf(ERRFILE,"WARNING--no proposed trees ever accepted\n");
  }
  free(changetree);
}  /* maketree */

/************************************************************
 * freenodelet frees all fields of a nodelet and the actual *
 * nodelet too.                                             */
void freenodelet(option_struct *op, data_fmt *data, node *p)
{
free_x(op,p);
if(op->map)free_z(op,p);
if (p->nayme) free(p->nayme);
ranges_Malloc(p,FALSE,0L);
coal_Malloc(p,FALSE,0L);
free(p);
} /* freenodelet */

/*******************************************************************
 * finalfree() frees                sum, reci, sametree,           *
 * theti,           the options structure, and the data structure. */
void finalfree(option_struct *op, data_fmt *data)
{
long i, j, k, chtype, numloci;

numloci = 0;

numloci = getdata_numloci(op,data);

free(reci[0]);
free(reci);
free(theti[0]);
free(theti);
free(sametree[0]);
free(sametree);

for(i = 0; i < numloci; i++)
   for(j = 0; j < 1+op->numchains[1]; j++) {
      if (j) chtype = 1;
      else chtype = 0;
      for(k = 0; k < op->numout[chtype]; k++) {
#if GROWTHUSED
         free(sum[i][j][k].kk);
         free(sum[i][j][k].kend);
         free(sum[i][j][k].actives);
#endif
         free(sum[i][j][k].eventtype);
         free(sum[i][j][k].sitescore);
      }
}
free(sum[0][0]);
free(sum[0]);
free(sum);

freedata(data);

free(op->probcat);
free(op->rate);
free(op->temperature);
if (!op->same_ne) free(op->ne_ratio);
if (op->panel) free(op->numpanel);
free(op);

} /* finalfree */


int main(int argc, char *argv[])
{  /* Recombine */
long i, numloci, numpop;
data_fmt data;
option_struct *options;
char infilename[100], outfilename[100], logfilename[100],
   bestrfilename[100], intrfilename[100], trfilename[100];

#ifdef MAC
char spafilename[100], seedfilename[100];

   argv[0] = "Recombine";
#endif

/* Open various filenames. */

openfile(&infile,INFILE,"r",argv[0],infilename);
openfile(&simlog,"simlog","w",argv[0],logfilename);
openfile(&outfile,OUTFILE,"w",argv[0],outfilename);

/* warning WARNING debug DEBUG--initialization handled wrong */
population = 0;

setupoption_struct(&options);
options->ibmpc = IBMCRT;
options->ansi = ANSICRT;
getoptions(options);
for (i = 1; i <= 1000; i++)
  clearseed = randum();
if (options->usertree)
  openfile(&intree,INTREE,"r",argv[0],intrfilename);
if (options->treeprint)
  openfile(&treefile,TREEFILE,"w",argv[0],trfilename);
openfile(&bestree,"bestree","w",argv[0],bestrfilename);

if (options->newdata) getinput(options,&data);
else getnodata(options,&data);

if (options->printdata) printdnadata(options,&data,outfile);

firstinit(options,&data);

if (ALWAYS_ACCEPT) fprintf(outfile,"***ALWAYS ACCEPT MODE!***");
if (ALWAYS_REJECT) fprintf(outfile,"***ALWAYS REJECT MODE!***");

if (options->newdata) {
   temptrees = (tree **)calloc(options->numtempchains,sizeof(tree *));
   numpop = getdata_numpop(options,&data);
   for(population = 0; population < numpop; population++) {
      popinit(options,&data);
      numloci = getdata_numloci(options,&data);
      for (locus = 0; locus < numloci; locus++) {
         if (options->progress) printf("Locus %ld\n",locus+1);
         fprintf(outfile,"\n---------------------------------------\n");
         fprintf(outfile, "      Locus %ld\n",locus+1);
         fprintf(outfile,"---------------------------------------\n");
         if (!options->same_ne)
            fprintf(outfile,"Assumed relative Ne = %f\n",
            options->ne_ratio[locus]);
         if (options->holding) {
            if (options->holding == 1)
               fprintf(outfile,"Theta was held contant.\n");
            else
               fprintf(outfile,"Recombination was held constant.\n");
         }
         locusinit(options,&data);
#if !ALWAYS_REJECT
         watttheta = watterson(options,&data);
#else
         watttheta = ARBITRARY_THETA;
#endif
         if (options->watt) theta0 = watttheta;
         maketree(options,&data);
         end_of_locus_free(options,&data);
      }
      locus--;
      if (numloci > 1) rec_estimate(options,&data,-1L,totchains-1,TRUE);
      end_of_population_free(options,&data);
   }
   free(temptrees);
} else {
   rec_scoreread(options,&data);
   numloci = getdata_numloci(options,&data);
   if (numloci > 1) rec_estimate(options,&data,-1L,totchains-1,TRUE);
   else rec_estimate(options,&data,0L,totchains-1,TRUE);
}

finalfree(options,&data);

FClose(infile);
FClose(outfile);
FClose(treefile);
FClose(bestree);
FClose(seedfile);
FClose(simlog);
FClose(spacefile);

#ifdef MAC
   strcpy(seedfilename,"seedfile");
   strcpy(spafilename,"spacefile");
   fixmacfile(outfilename);
   fixmacfile(infilename);
   fixmacfile(bestrfilename);
   fixmacfile(logfilename);
   fixmacfile(intrfilename);
   fixmacfile(seedfilename);
   fixmacfile(trfilename);
   fixmacfile(spafilename);
#endif

printf("PROGRAM DONE\n");
exit(0);

}  /* recombine */

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

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

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

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