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

#ifndef DROP_INCLUDE
#include "jdrop.h"
#endif

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

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

extern tree *curtree;
extern double theta0, rec0;
extern long population, locus;
extern FILE *simlog;

extern boolean istip(node *p);
extern node *otherparent(node *p);
extern void fixlength(option_struct *op, data_fmt *data, node *p);
extern void ltov(option_struct *op, data_fmt *data, node *p);
extern boolean branchsub(tlist *t, node *oldbranch, node *newbranch);
extern void init_ranges_alloc(long **ranges, long numranges);
extern boolean testratio(option_struct *op, data_fmt *data, tree
   *oldtree, tree *newtree, char ratiotype);
extern void ranges_Malloc(node *p, boolean allokate, long numrangepairs);
extern void subtymelist(tlist *t, node *branch, boolean all);
extern double getdata_space(option_struct *op, data_fmt *data,
   long target);
extern long countrbanches(option_struct *op, data_fmt *data,
   tree *tr);
extern void localeval(option_struct *op, data_fmt *data, node *p,
   boolean first);
extern void scoretree(option_struct *op, data_fmt *data, long chain);
extern double coalprob(option_struct *op, data_fmt *data, tree *tr,
   double theta, double r);

extern long locus; /* which locus are we currently working on? */
extern long numdropped; /* how many trees are dropped due to excessive
                           amounts of recombination */
extern long indecks, apps; /* what step and chain are we on? */
/* This file contains all the functions involved in modifying the
tree */


long countsites(option_struct *op, data_fmt *data)
{
long nummarkers;

nummarkers = getdata_nummarkers(op,data);

switch(op->datatype) {
   case 'a':
      break;
   case 'b':
   case 'm':
      break;
   case 'n':
      return(data->dnaptr->sitecount[locus][nummarkers-1]);
      break;
   case 's':
      return(nummarkers);
      break;
   default:
      fprintf(ERRFILE,"countsites, found an unknown datatype %c\n",
         op->datatype);
      break;
}

return(FLAGLONG);

} /* countsites */


void rename_branch(tlist *start, node *newbranch, node *oldbranch)
{
tlist *t;
boolean found;

if (!newbranch) return;

found = TRUE;
for(t = start; found && t != NULL; t = t->succ)
   found = branchsub(t,oldbranch,newbranch);

} /* rename_branch */


void remove_branch(tlist *start, node *target)
{
long i, j;
node *p;
tlist *t;
boolean found;

if (!target) return;

if(!target->top) p = target->back;
else p = target;

for(t = start, found = TRUE; t != NULL && found; t = t->succ) {
   found = FALSE;
   for(i = 0; i < t->numbranch && !found; i++) {
      if(t->branchlist[i] == p) {
          t->numbranch--;
          for(j = i; j < t->numbranch; j++)
             t->branchlist[j] = t->branchlist[j+1];
          found = TRUE;
      }
   }
}

} /* remove_branch */


/* OBSOLETE */
/*********************************************************************
 * subtymenode removes the tymenode containing the node "p" from the *
 * tymelist.  The branch pointed to by the nodelet "p" will be       *
 * completely erased from the branchlist.                            */
void subtymenode(tree *tr, node *p)
{
tlist *t;
node *q, *r;

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

q = findunique(p);
for(r = q->next; r == p; r = r->next) ;
if (isrecomb(p)) rename_branch(t,q->back,r);
else rename_branch(t,r->back,q);
remove_branch(t,p);

} /* subtymenode */


/********************************************************************
 * fix_treenodep() needs to called everytime a node is removed from *
 * the tree.  "Number" is the number of the node to be removed.     *
 * fix_treenodep assumes that the field  containing the number of   *
 * coalescences is correct.                                         */
void fix_treenodep(tree *tr, long number)
{
long i, numnodes;

numnodes = 2 * tr->numcoals + 2; /* strange but true! */
for(i = number; i < numnodes-1; i++) tr->nodep[i] = tr->nodep[i + 1];

} /* fix_treenodep */


boolean isdead(option_struct *op, node *p)
{

if (op->fc) return(!p->coal[0]);
else return(!p->ranges[0]);

} /* isdead */


node *findpaired_coalnode(option_struct *op, node *p)
{

if (istip(p)) {
   fprintf(ERRFILE,"ERROR:findpaired_coalnode found a tip!\n");
   exit(-1);
}

if (!isrecomb(p)) return(p);

printf("ERROR:Shouldn't be here in findpaired_coalnode!\n");
if (isdead(op,p->next)) findpaired_coalnode(op,p->next->back);
else if (isdead(op,p->next->next))
        findpaired_coalnode(op,p->next->next->back);
     else {fprintf(ERRFILE,"ERROR:findpaired_coalnode failed!\n"); exit(-1);}

fprintf(ERRFILE,"ERROR:findpaired_coalnode failure.  %ld %ld\n",indecks,apps);
return(NULL);

} /* findpaired_coalnode */


/**********************************************************************
 * remove_pairedcoalnode removes node "p" from tree "tr", and assumes *
 * that the branch to be removed is the one which has it's bottom at  *
 * the specific nodelet pointed to by "p".                            */
void remove_pairedcoalnode(option_struct *op, data_fmt *data,
   tree *tr, node *p)
{
tlist *t;
node *q;

tr->numcoals--;

hookup(p->next->next->back,p->next->back);
fixlength(op,data,p->next->back);
fix_treenodep(tr,p->number);

t = gettymenode(tr,p->number);
q = (p->next->top) ? p->next->next->back : p->next->back;
subtymelist(t,q,FALSE);
freenode(p);
} /* remove_pairedcoalnode */


/****************************************************************
 * remove_pairedrecombnode removes node "p" from tree "tr", and *
 * assumes that "p" points to the unique nodelet of the         *
 * recombination node.                                          */
void remove_pairedrecombnode(option_struct *op, data_fmt *data,
   tree *tr, node *p, node *dead)
{
node *q;
tlist *t;

tr->numrecombs--;

for(q = p->next; q == dead; q = q->next) ;
hookup(p->back,q->back);
fixlength(op,data,p->back);
fix_treenodep(tr,p->number);

t = gettymenode(tr,p->number);
subtymelist(t,dead,FALSE);
freenode(p);
} /* remove_pairedrecombnode */


void remove_pairednodes(option_struct *op, data_fmt *data, tree *tr,
   node *p, node *remove_thisway)
{
node *q;

if (!remove_thisway->top) {
   fprintf(ERRFILE,"ERROR:remove_pairednodes passed a non-top!\n");
   exit(-1);
}
q = findpaired_coalnode(op,remove_thisway->back);
remove_pairedcoalnode(op,data,tr,q);
remove_pairedrecombnode(op,data,tr,p,remove_thisway);

} /* remove_pairednodes */


long setnodenumber(boolean **nodenumber, long *numnodes)
{
long i;

for (i = 0; i < (*numnodes); i++)
   if (!(*nodenumber)[i]) {
      (*nodenumber)[i] = TRUE;
      return(i);
   }

(*numnodes)++;
*nodenumber = (boolean *)realloc(*nodenumber,(*numnodes)*sizeof(boolean));
(*nodenumber)[(*numnodes) - 1] = TRUE;

curtree->nodep = (node **)
   realloc(curtree->nodep,((*numnodes)+1)*sizeof(node *));

return((*numnodes) - 1);

} /*setnodenumber */

boolean foundbranch(tlist *t, node *p)
{
long i;

for(i = 0; i < t->numbranch; i++)
   if (t->branchlist[i] == p) return(TRUE);

return(FALSE);

} /* foundbranch */

void readdactive(tlist *stop, node *p)
{
long i, j;
tlist *t;

t = gettymenode(curtree,p->number);
do {
   if (!foundbranch(t,p)) {
      t->numbranch++;
      t->branchlist =
      (node **)realloc(t->branchlist, t->numbranch*sizeof(node *));
      t->branchlist[t->numbranch-1] = p;
   }
   t = t->succ;
} while (t != stop && t != NULL);

/* remove the FALSE branchlist entry from the rest of the tymelist */
while (t != NULL) {
   if (!foundbranch(t,p)) break;
   for (i = 0; i < t->numbranch; i++)
      if (t->branchlist[i] == p) break;
   for (j = i; j < t->numbranch - 1; j++)
      t->branchlist[j] = t->branchlist[j+1];
   t->numbranch--;
   t = t->succ;
}

} /* readdactive */

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

void freelin(linlist *t)
{
  free(t);
} /* freelin */


void freelinlist(linlist *t)
{
linlist *u, *usucc;

for(u = t; u != NULL; u = usucc) {usucc = u->succ; freelin(u);}

} /* freelinlist */


void addlinlist(linlist **t, node *target, double activelinks)
{
  linlist *r, *s, *u;

  s = *t;
  if (s == NULL) {
    newlin(&s);
    s->branchtop = target;
    s->activesites = activelinks;
    s->succ = NULL;
    s->prev = NULL;
    (*t)=s;
  } else {
    newlin(&r);
    r->branchtop = target;
    r->activesites = activelinks;
    /* insert new entry right after s */
    u = s->succ;
    s->succ = r;
    r->prev = s;
    r->succ = u;
    if(u!=NULL) u->prev = r;
  }
} /* addlinlist */


void sublinlist(linlist **t, node *target)
{
  boolean found;
  linlist *u;

  u = *t;
  found = FALSE;
  while (!found) {
    if (u->branchtop==target) {
      found = TRUE;
      if (u->succ != NULL) u->succ->prev = u->prev;
      if (u->prev != NULL) u->prev->succ = u->succ;
      /* if this was the first entry, reset the lineage pointer */
      if (u==(*t)) *t = u->succ;
      freelin(u);
    } else if (u->succ != NULL) u = u->succ;
    else fprintf(ERRFILE, "ERROR:SUBLINLIST failed to find target\n");
  }
} /* sublinlist */


/******************************************************
 * printlinlist() prints a linlist:  a debug function */
void printlinlist(linlist *u)
{
if (u) {
   do {
     if(u->branchtop!=NULL)
       fprintf(ERRFILE,"Branchtop %ld activesites %g\n",u->branchtop->number,
       u->activesites);
     else fprintf(ERRFILE,"This entry has no branchtop!\n");
     u=u->succ;
   } while (u!=NULL);
}

}  /* printlinlist */


/***************************************************************
 * foundlin() finds whether a given nodelet is already present *
 * in the passed linlist.  TRUE = found, FALSE = not found.    */
boolean foundlin(linlist *u, node *p)
{
for(;u != NULL; u = u->succ)
   if(u->branchtop == p) return (TRUE);

return(FALSE);

} /* foundlin */


void findlin(linlist **u, long which)
{
  long i;

  i = 1;
  do {
    if (i==which) return;
    (*u) = (*u)->succ;
    i++;
  } while ((*u)!=NULL);
  printf("ERROR:Unable to find linlist entry!\n");
} /* findlin */


/************************************************************************
 * getnumlinks() returns the number of links between the passed "sites" */
double getnumlinks(option_struct *op, data_fmt *data, long start, long end)
{
long i;
double sum, *spaces;

spaces = NULL;

switch(op->datatype) {
   case 'a':
      break;
   case 'b':
   case 'm':
      spaces = data->msptr->mspace[population][locus];
      break;
   case 'n':
      return((double)(end - start));
      break;
   case 's':
      spaces = data->dnaptr->sspace[population][locus];
      for(i = start, sum = 0.0; i < end; i++)
         sum += spaces[i];
      return(sum);
      break;
   default:
      fprintf(ERRFILE,"unknown datatype in getnumlinks\n");
      break;
}

return((double)FLAGLONG);

} /* getnumlinks */


/*******************************************************************
 * count_active() counts the number of active links leaving a node */
double count_active(option_struct *op, data_fmt *data, node *p)
{
double value;
node *q;

q = p;
if (!q->top) q = q->back;

if (istip(q))
return(getnumlinks(op,data,0L,countsites(op,data)-1));

if (isrecomb(q)) value = count_rec_active(op,data,q);
else value = count_coal_active(op,data,q);

return(value);

} /* count_active */


/*********************************************************************
 * count_coal_active() returns the number of active links along the  *
 * branch "p", where the node containing "p" must be coalescent; and *
 * the branch "p" may not exist yet.                                 */
double count_coal_active(option_struct *op, data_fmt *data, node *p)
{
long *range1, *range2, newstart, newend;

range1 = p->next->back->ranges;
range2 = p->next->next->back->ranges;

if (!range1[0] && !range2[0])
   fprintf(ERRFILE, "ERROR:Dead branch in count_coal_active! %ld %ld\n",
      indecks,apps);

if (!range1[0]) return(range2[2*range2[0]] - range2[1]);
if (!range2[0]) return(range1[2*range1[0]] - range1[1]);

newstart = (range1[1] < range2[1]) ? range1[1] : range2[1];
newend = (range1[2*range1[0]] > range2[2*range2[0]]) ?
            range1[2*range1[0]] : range2[2*range2[0]];

if (newend - newstart < 0)
   fprintf(ERRFILE,"ERROR:negative link count in count_coal_active! %ld %ld\n",
      indecks,apps);

return(getnumlinks(op,data,newstart,newend));

} /* count_coal_active */


/*******************************************************************
 * count_rec_active() returns the number of active links along the *
 * branch "p", where the node containing "p" must be recombinant;  *
 * and the branch "p" may not exist yet.                           */
double count_rec_active(option_struct *op, data_fmt *data, node *p)
{
long i, newstart, newend;
node *q;

q = findunique(p)->back;

if (!q->ranges[0])
   fprintf(ERRFILE,"ERROR:count_rec_active counted dead, %ld %ld\n",
      indecks,apps);

for(i = 1, newstart = -1L; q->ranges[i] != FLAGLONG; i+=2) {
   if(p->recstart > q->ranges[i+1]) continue;
   newstart = (q->ranges[i] > p->recstart) ? q->ranges[i] : p->recstart;
   break;
}

if(newstart == -1L) return(0L);

for(i=2*q->ranges[0], newend = -1L; i != 0; i-=2) {
   if(p->recend < q->ranges[i-1]) continue;
   newend = (q->ranges[i] < p->recend) ? q->ranges[i] : p->recend;
   break;
}

if(newend == -1L) return(0L);

if (newend - newstart < 0)
   fprintf(ERRFILE,"ERROR:negative link count in count_rec_active! %ld %ld\n",
      indecks,apps);

return(getnumlinks(op,data,newstart,newend));

} /* count_rec_active */


/*********************************************************************
 * count_activefc() counts the number of active links leaving a node *
 * taking into account which sites have achieved final coalescence.  */
double count_activefc(option_struct *op, data_fmt *data, node *p)
{
double value;
node *q;

q = p;
if (!q->top) q = q->back;

if (istip(q))
   return(getnumlinks(op,data,0L,countsites(op,data)-1));

if (op->datatype == 'n')
   value = ((q->coal[0]) ? q->coal[2*q->coal[0]] - q->coal[1] : 0L);
else
   value = ((q->coal[0]) ?
   getnumlinks(op,data,q->coal[1],q->coal[2*q->coal[0]]) : 0L);

return(value);

} /* count_activefc */


/******************************************************************
 * is_fc() returns TRUE if site has achieved fc by treesec, FALSE *
 * otherwise.                                                     */
boolean is_fc(tlist *treesec, long site)
{
long i, count;

for(i = 0, count = 0; i < treesec->numbranch; i++) {
   if(inrange(treesec->branchlist[i]->ranges,site)) count++;
   if (count > 1) return(FALSE);
}

return(TRUE);

} /* is_fc */


/***********************************************************
 * fix_coal() upodates the coal array for the passed node. *
 * "tfc" is the tymeslice used for FC calculations and all *
 * tipwards tymeslices plus their eventnode's coal arrays  *
 * are assumed to be correct.                              */
void fix_coal(option_struct *op, data_fmt *data, tree *tr, tlist *tfc,
  node *p)
{
long *subtrees;
node *q, *r, *s;

if (isrecomb(p)) {
   q = findunique(p);
   r = q->next;
   s = q->next->next;
   q = q->back;
   copycoal(q,r);
   copycoal(q,s);
   subrangefc(&(r->coal),s->recstart,s->recend);
   subrangefc(&(s->coal),r->recstart,r->recend);
} else {
   subtrees = NULL;
   q = findunique(p);
   findsubtrees_node(op,data,tfc,tr,&subtrees);
   makecoal(op,data,tr,tfc,&(q->coal),q,subtrees);
   free(subtrees);
}

} /* fix_coal */


/******************************************************************
 * makecoal() creates the appropiate coal array for the passed in *
 * tymeslice.                                                     *
 * makecoal() is nodelet specific, and the passed in nodelet must *
 * have correct information contained in its ranges array.        *
 * makecoal() is nodelet specific iff "p" is passed in as non-NULL*/
void makecoal(option_struct *op, data_fmt *data, tree *tr,
   tlist *tfc, long **coal, node *p, long *subtrees)
{
long i;

init_coal_alloc(coal,1L);
(*coal)[1] = 0;
(*coal)[2] = countsites(op,data)-1;

for(i = 0; i < subtrees[0]; i++) {
   if (p)
      if (!inrange(p->ranges,subtrees[2*i+1])) {
         subrangefc(coal,subtrees[2*i+1],subtrees[2*i+2]);
         continue;
      }

   if (is_fc(tfc,subtrees[2*i+1]))
      subrangefc(coal,subtrees[2*i+1],subtrees[2*i+2]);

}

} /* makecoal */


/* three functions involving updating the alias array "siteptr" */

void edit_alias(option_struct *op, data_fmt *data, long *sp, long cutpoint)
{
  long i, alias_site, cutmarker, numsites;

  /* this looks awful because sp[i] does not actually contain
   the alias site, it contains the alias site + 1 (to avoid zero). */

  /* cutpoint is the first site *after* the recombination */

  numsites = getdata_nummarkers(op,data);
  cutmarker = sitetorightmarker(op,data,cutpoint);
  if (cutmarker == FLAGLONG) return;

  for(i=cutmarker;i<numsites;i++) {
    alias_site = sp[i]-1;
    if(alias_site < cutmarker && alias_site+1 > 0) sp[i]*= -1;
  }
} /* edit_alias */


void traverse_rebuild_alias(option_struct *op, data_fmt *data,
  tlist *tstart, long *sp)
{
tlist *t;
node *p;

for(t = tstart; t != NULL; t = t->succ) {
   p = t->eventnode;
   if (!isrecomb(p)) continue;
   if (p->recstart == 0) edit_alias(op,data,sp,p->recend+1);
   else edit_alias(op,data,sp,p->recstart);
}

} /* traverse_rebuild_alias */


void rebuild_alias(option_struct *op, data_fmt *data, long *sp)
{
  long i, numsites;

  numsites = getdata_nummarkers(op,data);

  for(i=0;i<numsites;i++) if (sp[i]<0) sp[i]*= -1;
  traverse_rebuild_alias(op,data,curtree->tymelist,sp);
} /* rebuild_alias */


/****************************************************************
 * contrib() sets the passed in "ranges" field, using "p".  "p" *
 * points to the nodelet at the top of the relevant branch.     */
void contrib(option_struct *op, data_fmt *data, node *p, long **newranges)
{
node *q, *r;
long i, numremove;


if(istip(p)) {
   (*newranges)[0] = 1L;
   (*newranges)[1] = 0L;
   (*newranges)[2] = countsites(op,data)-1;
   (*newranges)[3] = FLAGLONG;
   return;
}

if(isrecomb(p)) {
   q = findunique(p)->back;
   init_ranges_alloc(newranges,q->ranges[0]);
   memcpy((*newranges),q->ranges,(2*q->ranges[0]+2)*sizeof(long));
/* first remove the leading excess ranges */
   for(i = 1, numremove = 0; (*newranges)[i] != FLAGLONG; i+=2) {
      if(p->recstart > (*newranges)[i+1]) {numremove+=2; continue;}
      if(p->recstart > (*newranges)[i]) (*newranges)[i] = p->recstart;
      break;
   }
   memmove(&(*newranges)[1],&(*newranges)[1+numremove],
      (2*(*newranges)[0]-numremove+1)*sizeof(long));
   (*newranges)[0] -= numremove/2;
/* then removing the trailing excess ranges */
   for(i = 2*(*newranges)[0]-1, numremove = 0; i+1 != 0; i-=2) {
      if(p->recend < (*newranges)[i]) {numremove++; continue;}
      if(p->recend < (*newranges)[i+1]) (*newranges)[i+1] = p->recend;
      break;
   }
   (*newranges)[i+2] = FLAGLONG;
   (*newranges)[0] -= numremove;
} else {
   q = p->next->back;
   r = p->next->next->back;
   init_ranges_alloc(newranges,q->ranges[0]);
   memcpy((*newranges),q->ranges,(2*q->ranges[0]+2)*sizeof(long));
   for(i = 1; r->ranges[i] != FLAGLONG; i+=2)
      addrange(newranges,r->ranges[i],r->ranges[i+1]);
}

} /* contrib */


/***************************************************************
 * fix_coal_ranges() updates the "ranges" of a coalescent node */
void fix_coal_ranges(option_struct *op, data_fmt *data, node *p)
{
  p = findunique(p);
  contrib(op,data,p,&p->ranges);
} /* fix_coal_ranges */


/***************************************************************
 * fix_rec_ranges() updates the "ranges" of a recombinant node */
void fix_rec_ranges(option_struct *op, data_fmt *data, node *p)
{
  p = findunique(p);
  contrib(op,data,p->next,&p->next->ranges);
  contrib(op,data,p->next->next,&p->next->next->ranges);
} /* fix_rec_ranges */


boolean popbranch(tlist *t, long branch)
{
  long i;
  boolean found;

  found = FALSE;
  for(i=0;i<t->numbranch;i++) {
    if (t->branchlist[i]!=NULL) {
      if (t->branchlist[i]->number==branch) {
        t->branchlist[i]=NULL;
        found = TRUE;
      }
    }
  }
  return(found);
} /* popbranch */


boolean subbranch(tlist *t, long oldbranch, node *newbranch)
{
  long i;
  boolean found;

  found = FALSE;
  for(i=0;i<t->numbranch;i++) {
    if (t->branchlist[i]!=NULL) {
      if (t->branchlist[i]->number==oldbranch) {
        t->branchlist[i]=newbranch;
        found = TRUE;
      }
    }
  }
  return(found);
} /* subbranch */


void poptymelist (node *p)
{  /* remove futile tymelist entry */
   boolean done;
   node *q;
   tlist *s, *t, *u;

   t = gettymenode(curtree, p->number);
   u = t;
   if(isrecomb(p)) {
     s = u->succ;
     do {
       if(s==NULL) break;
       done = !popbranch(s,p->number);
       s = s->succ;
     } while (!done);
     u->prev->age=t->age;
     u->prev->succ=u->succ;
     if (t->succ != NULL) u->succ->prev=u->prev;
     freetymenode(u);
   } else {
     s = u->succ;
     if(p->next->top) q=p->next->next->back;
     else q=p->next->back;
     do {
       if(s==NULL) break;
       done = !subbranch(s,p->number,q);
       s = s->succ;
     } while (!done);
     u->prev->age=t->age;
     u->prev->succ=u->succ;
     if (u->succ != NULL) u->succ->prev=u->prev;
     freetymenode(u);
   }
} /* poptymelist */


void fixbranch(tlist *t)
{
  long i, j, nullcnt, newnum;
  node **newarray;

  nullcnt = 0;
  for(i=0;i<t->numbranch;i++) {
    if(t->branchlist[i]==NULL) nullcnt++;
  }
  if(nullcnt==0) return;
  /* workaround to avoid callocing 0 (or less) elements */
  newnum = (t->numbranch-nullcnt < 1) ? 1 : t->numbranch-nullcnt;
  newarray = (node **)calloc(1,newnum*sizeof(node *));
  j = 0;
  for(i=0;i<t->numbranch;i++) {
    if(t->branchlist[i]!=NULL) {
      newarray[j]=t->branchlist[i];
      j++;
    }
  }
  free(t->branchlist);
  t->numbranch = newnum;
  t->branchlist=newarray;
} /* fixbranch */


/***********************************************************
 * is_futile() checks to see if a node has been tagged for *
 * futile removal.                                         */
boolean is_futile(node *p)
{

return((p->futileflag == FLAGLONG));

} /* is_futile */


/**********************************************************
 * find_nonfutile() finds the first non-futile node along *
 * the path of nodes going in the direction indicated by  *
 * upwards (TRUE = up, FALSE = down).                     */
node *find_nonfutile(node *p, node *cutnode, boolean upwards)
{
node *q;

if (p == cutnode) return(NULL);
if (istip(p) || !is_futile(p)) return(p);

p = findunique(p);
if (isrecomb(p)) {
   if (upwards) q = find_nonfutile(p->back,cutnode,upwards);
   else {
      q = find_nonfutile(p->next->back,cutnode,upwards);
      if (!q) q = find_nonfutile(p->next->next->back,cutnode,upwards);
   }
} else {
   if (upwards) {
      q = find_nonfutile(p->next->back,cutnode,upwards);
      if (!q) q = find_nonfutile(p->next->next->back,cutnode,upwards);
   } else q = find_nonfutile(p->back,cutnode,upwards);
}

return(q);

} /* find_nonfutile */


/**********************************************************
 * tag_futile() sets a flag to indicate that it should be *
 * removed.  The flag set is for recstart & recend set to *
 * FLAGLONG;                                              */
void tag_futile(node *p)
{

if (isrecomb(p)) {
   p->futileflag = p->next->futileflag = FLAGLONG;
   p->next->next->futileflag = FLAGLONG;
   tag_futile(p->next->back);
   if (p->next->next->back != curtree->root)
      tag_futile(p->next->next->back);
   return;
} else {
   if (!is_futile(p)) {
      p->futileflag = p->next->futileflag = FLAGLONG;
      p->next->next->futileflag = FLAGLONG;
      return;
   } else {
      if (p->next->top) tag_futile(p->next->back);
      else tag_futile(p->next->next->back);
   }
}

} /* tag_futile */


/**********************************************************************
 * renamebrlist() goes through a brlist replacing old branchtops with *
 * the new branchtop.                                                 */
void renamebrlist(brlist **bigbr, brlist *start, node *oldtop, node *newtop)
{
brlist *br, *brsucc;

for(br = start; br != NULL; br = brsucc) {
   brsucc = br->succ;
   if (br->branchtop == oldtop) {
      br->branchtop = newtop;
      if (newtop->tyme >= br->endtyme) {
         br_remove(bigbr,br);
         continue;
      }
      if (br->starttyme < newtop->tyme) {
         br->starttyme = newtop->tyme;
         if (br->prev) br->prev->succ = br->succ;
         else (*bigbr) = br->succ;
         if (br->succ) br->succ->prev = br->prev;
         hookup_brlist(bigbr,br);
      }
   }
}

} /* renamebrlist */


/**********************************************************************
 * drop_renamebrlist() is a coalesce specific driver for renamebrlist */
void drop_renamebrlist(brlist **bigbr, node *oldtop, node *newtop)
{
brlist *br;

if (!bigbr) return;

br = (*bigbr);

renamebrlist(bigbr,br,oldtop,newtop);

} /* drop_renamebrlist */


/********************************************************************
 * drop_brfix() updates the ranges & coal information, and changes  *
 * brlist info, for all branches which start at or after the passed *
 * in tymelist's eventnode.                                         */
void drop_brfix(option_struct *op, data_fmt *data, tree *newtree,
   tlist *start, brlist **oldbr)
{
node *p, *q, *r;
tlist *t;
long *oldcoal1, *oldcoal2;
dnadata *dna;

oldcoal1 = oldcoal2 = NULL;
dna = data->dnaptr;

for (t = start; t != NULL; t = t->succ) {
   p = findunique(t->eventnode);
   if (istip(p)) continue;
   if (isrecomb(p)) {
      q = p->next;
      r = p->next->next;
      if (op->fc) {
         fix_rec_ranges(op,data,p);
         init_coal_alloc(&oldcoal1,q->coal[0]);
         init_coal_alloc(&oldcoal2,r->coal[0]);
         memcpy(oldcoal1,q->coal,(2*q->coal[0]+2)*sizeof(long));
         memcpy(oldcoal2,r->coal,(2*r->coal[0]+2)*sizeof(long));
         fix_coal(op,data,newtree,t,p);
         if (!sameranges(oldcoal1,q->coal))
            addbrlist(op,data,oldbr,q,oldcoal1,q->coal,NULL,q->tyme,
               q->back->tyme,FALSE);
         if (!sameranges(oldcoal2,r->coal))
            addbrlist(op,data,oldbr,r,oldcoal2,r->coal,NULL,r->tyme,
               r->back->tyme,FALSE);
      } else {
         init_ranges_alloc(&oldcoal1,q->ranges[0]);
         init_ranges_alloc(&oldcoal2,r->ranges[0]);
         memcpy(oldcoal1,q->ranges,(2*q->ranges[0]+2)*sizeof(long));
         memcpy(oldcoal2,r->ranges,(2*r->ranges[0]+2)*sizeof(long));
         fix_rec_ranges(op,data,p);
         if (!sameranges(oldcoal1,q->ranges))
            addbrlist(op,data,oldbr,q,oldcoal1,q->ranges,NULL,q->tyme,
               q->back->tyme,FALSE);
         if (!sameranges(oldcoal2,r->ranges))
            addbrlist(op,data,oldbr,r,oldcoal2,r->ranges,NULL,r->tyme,
               r->back->tyme,FALSE);
      }
   } else {
      if (op->fc) {
         fix_coal_ranges(op,data,p);
         init_coal_alloc(&oldcoal1,p->coal[0]);
         memcpy(oldcoal1,p->coal,(2*p->coal[0]+2)*sizeof(long));
         fix_coal(op,data,newtree,t,p);
         if (!sameranges(oldcoal1,p->coal))
            addbrlist(op,data,oldbr,p,oldcoal1,p->coal,NULL,
               p->tyme,p->back->tyme,FALSE);
      } else {
         init_ranges_alloc(&oldcoal1,p->ranges[0]);
         memcpy(oldcoal1,p->ranges,(2*p->ranges[0]+2)*sizeof(long));
         fix_coal_ranges(op,data,p);
         if (!sameranges(oldcoal1,p->ranges))
            addbrlist(op,data,oldbr,p,oldcoal1,p->ranges,NULL,
               p->tyme,p->back->tyme,FALSE);
      }
   }
}

if (oldcoal1) free(oldcoal1);
if (oldcoal2) free(oldcoal2);

} /* drop_brfix */


/************************************************************
 * samesites() checks to see if the segranges array and the *
 * ranges array could be describing the same set of sites.  */
boolean samesites(int *segranges, long *ranges, long numsites)
{
long i;

for(i = 0; i < numsites; i++) {
   if(segranges[i] == 1 || segranges[i] == 2) {
      if(!inrange(ranges,i)) return(FALSE);
   } else {
      if(inrange(ranges,i)) return(FALSE);
   }
}

return(TRUE);

} /* samesites */


/****************************************************************
 * treelegalbrlist() checks to make sure that a brlist entry is *
 * legal given the tree.                                        */
void treelegalbrlist(option_struct *op, data_fmt *data, tree *tr,
   brlist **bigbr, brlist *target)
{
long *comp;
node *top, *bottom, *treenode;

top = target->branchtop;
bottom = top->back;
treenode = tr->nodep[top->number];

if (isrecomb(top)) {
   if(treenode != top && treenode != otherparent(top)) {
      br_remove(bigbr,target);
      return;
   }
} else {
   if(treenode != top) {
      br_remove(bigbr,target);
      return;
   }
}

if (target->endtyme > bottom->tyme) {
   if (isrecomb(bottom)) {
      bottom = findunique(bottom);
      if (!bottom->next->back) top = bottom->next->next;
      else if (!bottom->next->next->back) top = bottom->next;
      else {
         comp = ((op->fc) ? bottom->next->coal : bottom->next->ranges);
         top =
         (samesites(target->segranges,comp,countsites(op,data))
         ? bottom->next : bottom->next->next);
      }
   } else top = findunique(bottom);

   if (target->endtyme < curtree->root->back->tyme)
      addbrlist(op,data,bigbr,top,NULL,NULL,target->segranges,bottom->tyme,
         target->endtyme,FALSE);
   target->endtyme = bottom->tyme;
}

} /* treelegalbrlist */


/******************************************************************
 * samebranch() checks to see if the two segments are on the same *
 * branch (i.e. whether they have the same tops).                 */
boolean samebranch(brlist *br1, brlist *br2)
{

return(br1->branchtop == br2->branchtop);

} /* samebranch */


/********************************************************************
 * checktreebranch() checks to see if a segment is properly part of *
 * the branch whose top is p.                                       *
 *                                                                  *
 * It returns: 0 if no match is found                               *
 *             1 if starts match                                    *
 *            -1 if ends match                                      *
 *             2 if both start and end match                        */
long checktreebranch(brlist *br, node *p, double starttyme,
   double endtyme)
{

if (br->branchtop != p) return(0L);

if (br->starttyme == starttyme && br->endtyme == endtyme) return(2L);
if (br->starttyme == starttyme) return(1L);
if (br->endtyme == endtyme) return(-1L);

return(0L);

} /* checktreebranch */


/********************************************************************
 * sametreebranch() checks to see if a segment was ever on the same *
 * branch as the branch possibly bisected by the node p.            *
 *                                                                  *
 * It returns: 0 if no match is found                               *
 *             1 if starts match                                    *
 *            -1 if ends match                                      *
 *             2 if both start and end match                        */
long sametreebranch(brlist *br, node *p, double starttyme,
   double endtyme)
{
node *target;
long match;

match = checktreebranch(br,p,starttyme,endtyme);

if (match /*|| istip(p)*/) return(match);

target = findunique(p);
if (isrecomb(target)) {
   match = checktreebranch(br,target->back,starttyme,endtyme);
} else {
   match = checktreebranch(br,target->next->back,starttyme,endtyme);
   if (match) return(match);
   match = checktreebranch(br,target->next->next->back,starttyme,endtyme);
}

return(match);

} /* sametreebranch */


/****************************************************************
 * isoverlap() checks to see if two brlist entries overlap each *
 * other on the tree.                                           */
boolean isoverlap(brlist *br1, brlist *br2)
{
double br1start, br2start;

if (!samebranch(br1,br2)) return(FALSE);

br1start = br1->starttyme;
br2start = br2->starttyme;

if (((br1start > br2start) && (br1start < br2->endtyme)) ||
    ((br2start > br1start) && (br2start < br1->endtyme))) return(TRUE);

return(FALSE);


} /* isoverlap */


/******************************************************************
 * iscontainedin() checks to see if one brlist entry is entirely  *
 * within another.                                                */
boolean iscontainedin(brlist *br1, brlist *br2)
{

if (!samebranch(br1,br2)) return(FALSE);

if(((br1->starttyme > br2->starttyme) && (br1->endtyme < br2->endtyme))
  || ((br2->starttyme > br1->starttyme) && (br2->endtyme < br1->endtyme)))
  return(TRUE);

return(FALSE);

} /* iscontainedin */


/****************************************************************
 * mash_brlist() searches a brlist for a match to its "target". *
 * If it finds one, it mashes the two entries together.         */
void mash_brlist(tree *tr, brlist *start, brlist *target)
{
brlist *br, *brupper, *brlower;
double temp;

for(br = start; br != NULL; br = br->succ) {
   if(!isoverlap(target,br) && !iscontainedin(target,br)) continue;
   if(target->starttyme < br->starttyme) {
      brupper = target;
      brlower = br;
   } else {
      brupper = br;
      brlower = target;
   }
   if (iscontainedin(target,br)) {
      brupper->endtyme = brlower->starttyme;
      continue;
   }
   temp = brupper->endtyme;
   brupper->endtyme = brlower->starttyme;
   brlower->starttyme = temp;
}

} /* mash_brlist */


/*****************************************************************
 * consolidate_brlist() attempts to remove and replace redundate *
 * entries from a brlist.  The tree must be correct.             */
void consolidate_brlist(option_struct *op, data_fmt *data, tree *tr,
   brlist **bigbr)
{
brlist *br, *brsucc;

if (!bigbr) return;

for(br = (*bigbr); br != NULL; br = brsucc) {
   brsucc = br->succ;
   treelegalbrlist(op,data,tr,bigbr,br);
}

for(br = (*bigbr); br != NULL; br = br->succ) {
   mash_brlist(tr,br->succ,br);
}

} /* consolidate_brlist */


/******************************************************************
 * br_remove() performs a very basic removal of the passed brlist *
 * element from whatever brlist it is attached to, then frees it. */
void br_remove(brlist **bigbr, brlist *target)
{

if (target->prev) target->prev->succ = target->succ;
else (*bigbr) = target->succ;
if (target->succ) target->succ->prev = target->prev;

freebrlist(target);

} /* br_remove */


/*************************************
 * printbrlist() prints out a brlist */
void printbrlist(brlist *start)
{
brlist *br;

fprintf(ERRFILE,"Brlist follows\n");
for(br = start; br != NULL; br = br->succ) {
   fprintf(ERRFILE,"   top %3ld at tyme %12.8f [%12.8f/%12.8f]",
      br->branchtop->number,br->branchtop->tyme,br->starttyme, br->endtyme);
   fprintf(ERRFILE," with newlinks %f\n",br->numnewactives);
}

} /* printbrlist */


/*********************************************************************
 * find_futilecoaldtr() returns the "traversed" daughter node from a *
 * the passed node "p".                                              */
node *find_futilecoaldtr(node *cutnode, node *p)
{
node *q;

q = findtop(p);

if (is_futile(q->next->back) || q->next->back == cutnode)
   return(q->next->next->back);
else return(q->next->back);

} /* find_futilecoaldtr */


/********************************************************************
 * remove_futile() removes the nodes from the tree that have been   *
 * tagged by tag_futile, while maintaing brlist entries for them.   *
 * remove_futile() assumes that the first entry in the tymelist     *
 * contains the tips and so will never be removed.  It also assumes *
 * that coalescent nodes require tree clean-up, whereas recombinant *
 * nodes may be simply removed.                                     *
 * Traverse the tymelist in reverse order so that path of removal   *
 * info is preserved.                                               */
void remove_futile(option_struct *op, data_fmt *data, tree *newtree,
   tree *oldtree, node *cutnode, boolean *nodenumber, brlist **br)
{
long i, numnodes, **oranges, *oldranges, topcount, *actives;
double btyme;
tlist *t, *tsucc, *firstfutile;
node *p, *q, *dtr, *par, **obnode, **otnode, *oldnode, *newtopnode,
  **tops;
dnadata *dna;

dna = data->dnaptr;

numnodes = 2 * oldtree->numcoals + 2;
obnode = (node **)calloc(numnodes,sizeof(node *));
otnode = (node **)calloc(numnodes,sizeof(node *));
tops = (node **)calloc(numnodes,sizeof(node *));
oranges = (long **)calloc(numnodes,sizeof(long *));

topcount = 0;
firstfutile = NULL;

/* scanning tymelist from the bottom up */
for(t = newtree->tymelist; t->succ != NULL; t = t->succ) ;

/* modify the basic underlying tree structure */
for(; t->prev != NULL; t = t->prev) {
   obnode[t->eventnode->number] = otnode[t->eventnode->number] = NULL;
   if (!is_futile(t->eventnode)) continue;
   firstfutile = t;
   p = t->eventnode;
   nodenumber[p->number] = FALSE;
   if (isrecomb(p)) {
      newtree->numrecombs--;
      q = findunique(p);
      otnode[p->number] = (is_futile(q->next->back)) ?
                             q->next : q->next->next;
   } else {
      obnode[p->number] = p->back;
      dtr = find_nonfutile(p,cutnode,TRUE);
      otnode[p->number] = dtr;
      if (dtr != NULL) {
        tops[topcount] = dtr;
        topcount++;
      }
      actives = (op->fc) ? p->coal : p->ranges;
      init_coal_alloc(&oranges[p->number],actives[0]);
      memcpy(oranges[p->number],actives,(2*actives[0]+2)*sizeof(long));
      if (!is_futile(p->back)) {
         par = p->back;
         hookup(dtr,par);
         fixlength(op,data,dtr);
      }
      newtree->numcoals--;
   }
}

/* remove futile branches from the branchlists */
for(t = firstfutile; t != NULL; t = tsucc) {
   tsucc = t->succ;
   if (!is_futile(t->eventnode)) continue;
   p = t->eventnode;
   if (isrecomb(p)) {
      if (p == cutnode || otherparent(p) == cutnode) {
         if (p == cutnode) {
            remove_branch(t,p);
            rename_branch(t,findunique(p)->back,otherparent(p));
         } else {
            remove_branch(t,otherparent(p));
            rename_branch(t,findunique(p)->back,p);
         }
      } else {
         remove_branch(t,p);
         remove_branch(t,otherparent(p));
      }
   } else {
      if (otnode[p->number]) rename_branch(t,otnode[p->number],p);
      else remove_branch(t,p);
   }
}

extbranch(curtree,gettymenode(newtree,cutnode->number),cutnode);

/* scanning tymelist from the top down, updating coal and ranges arrays */
for (t = firstfutile; t != NULL; t = t->succ) {
   p = findunique(t->eventnode);
   if (!is_futile(p)) {
      if (isrecomb(p)) {
         fix_rec_ranges(op,data,p);
      } else {
         fix_coal_ranges(op,data,p);
      }
      if (op->fc) fix_coal(op,data,newtree,t,p);
   }
}


/* scanning tymelist from the top down again, building the brlist,
   and finally removing the remnant tymelist entries */
for (t = firstfutile; t != NULL; t = tsucc) {
   p = findunique(t->eventnode);
   tsucc = t->succ;
   if (!is_futile(p)) {
      oldnode = findunique(oldtree->nodep[p->number]);
      if (isrecomb(p)) {
         if (op->fc) {
            oldranges = oldnode->next->coal;
            newtopnode = p->next;
            btyme = oldnode->next->back->tyme;
            addbrlist(op,data,br,newtopnode,oldranges,newtopnode->coal,
               NULL,newtopnode->tyme,btyme,TRUE);

            oldranges = oldnode->next->next->coal;
            newtopnode = p->next->next;
            btyme = oldnode->next->next->back->tyme;
            addbrlist(op,data,br,newtopnode,oldranges,newtopnode->coal,
               NULL,newtopnode->tyme,btyme,TRUE);
         } else {
            oldranges = oldnode->next->ranges;
            newtopnode = p->next;
            btyme = oldnode->next->back->tyme;
            addbrlist(op,data,br,newtopnode,oldranges,newtopnode->ranges,
               NULL,newtopnode->tyme,btyme,TRUE);

            oldranges = oldnode->next->next->ranges;
               newtopnode = p->next->next;
            btyme = oldnode->next->next->back->tyme;
            addbrlist(op,data,br,newtopnode,oldranges,newtopnode->ranges,
               NULL,newtopnode->tyme,btyme,TRUE);
         }

      } else {
         if (op->fc) {
            oldranges = oldnode->coal;
            newtopnode = p;
            btyme = oldnode->back->tyme;
            addbrlist(op,data,br,newtopnode,oldranges,newtopnode->coal,
               NULL,newtopnode->tyme,btyme,TRUE);
         } else {
            oldranges = oldnode->ranges;
            newtopnode = p;
            btyme = oldnode->back->tyme;
            addbrlist(op,data,br,newtopnode,oldranges,newtopnode->ranges,
               NULL,newtopnode->tyme,btyme,TRUE);
         }
      }
   } else {
      if (!isrecomb(p)) {
         if (op->fc) {
            addbrlist(op,data,br,otnode[p->number],oranges[p->number],
               obnode[p->number]->back->coal,NULL,p->tyme,
               obnode[p->number]->tyme,TRUE);
         } else {
            addbrlist(op,data,br,otnode[p->number],oranges[p->number],
               obnode[p->number]->back->ranges,NULL,p->tyme,
               obnode[p->number]->tyme,TRUE);
         }
      }
      if (isrecomb(p)) subtymelist(t,p->next,TRUE);
      else subtymelist(t,p,TRUE);
      freenode(p);
   }
}

/* The following loop deals with the case (and probably others as well)
   where the cutnode is the top of a recombinant loop and the non-cut
   branch needs to be extended down to replace the bottom nodes' branch
*/
for(i=0; i<topcount; i++) {
    t = gettymenode(newtree,tops[i]->number);
    readdbranch(newtree,t,tops[i]);
  }

cutnode->back = NULL;

free(obnode);
free(otnode);
free(tops);
for(i = 0; i < numnodes; i++) if (oranges[i]) free(oranges[i]);
free(oranges);

} /* remove_futile */


/**************************************************************
*  Re-introduces a branch into the tymelist after it has      *
*  been stripped out during remove_futile.  This code is      *
*  quite specific to remove_futile; don't use it generically. */
void readdbranch(tree *newtree, tlist *u, node *addbranch)
{
  long i;
  tlist *t;
  boolean found;

  for (t = u; t!=NULL; t=t->succ) {
    if(t->eventnode->number==addbranch->back->number) return;
    found = FALSE;
    for (i = 0; i < t->numbranch; i++) {
      if (t->branchlist[i]==addbranch) found = TRUE;
    }
    if (!found) {
      t->numbranch++;
      t->branchlist = (node **)realloc(t->branchlist,
        t->numbranch*sizeof(node *));
      t->branchlist[t->numbranch - 1] = addbranch;
    }
  }
} /* addbranch */


/***********************************************************************
 * rembranch() removes the passed branch from all tymelist branchlists *
 * from "start" to the bottom/end of the tree.                         */
void rembranch(tree *newtree, tlist *start, node *sbranch)
{
long i;
tlist *t;
boolean found;
node **newlist;

for (t = start; t!=NULL; t=t->succ) {
    newlist = (node **)calloc(t->numbranch,sizeof(node *));
    found = FALSE;
    for(i = 0; i < t->numbranch; i++) {
       if (!found) newlist[i] = t->branchlist[i];
       else newlist[i-1] = t->branchlist[i];
       if (t->branchlist[i] == sbranch) found = TRUE;
    }
    free(t->branchlist);
    t->branchlist = newlist;
    if (found) t->numbranch--;
}

} /* rembranch */



/******************************************************************
 * extbranch() adds the passed branch to all tymelist branchlists *
 * from "start" to the bottom/end of the tree.                    */
void extbranch(tree *newtree, tlist *start, node *addbranch)
{
  long i;
  tlist *t;
  boolean found;

  for (t = start; t!=NULL; t=t->succ) {
    found = FALSE;
    for (i = 0; i < t->numbranch; i++) {
      if (t->branchlist[i]==addbranch) found = TRUE;
    }
    if (!found) {
      t->numbranch++;
      t->branchlist = (node **)realloc(t->branchlist,
        t->numbranch*sizeof(node *));
      t->branchlist[t->numbranch - 1] = addbranch;
    }
  }

} /* extbranch */


boolean isactive(node *p, linlist *u)
{
linlist *t;

for(t = u; t != NULL; t = t->succ)
   if (p == t->branchtop) return(TRUE);

return(FALSE);

} /* isactive */


long countinactive(option_struct *op, tlist *t, linlist *u)
{
long i, numinactive, *actives;

numinactive = t->numbranch;
for (i = 0; i < t->numbranch; i++) {
   if (isactive(t->branchlist[i],u)) {numinactive--; continue;}
   actives = (op->fc) ? t->branchlist[i]->coal :
                        t->branchlist[i]->ranges;
   if (!actives[0]) {numinactive--; continue;}
}

return(numinactive);

} /* countinactive */


/*********************************************************************
 * pickbrtarget() returns a randum brlist segment from the passed in *
 * brlist using the passed in number of active sites as a guide.     */
brlist *pickbrtarget(brlist *start, double numactive)
{
double pick;
brlist *br;

pick = randum() * numactive - start->numnewactives;
for(br = start; pick > 0; br = br->succ) pick -= br->succ->numnewactives;

return(br);

} /* pickbrtarget */


/*********************************************************************
 * pickbrcross() picks a randum point of crossover within the passed *
 * in brsegment.                                                     */
long pickbrcross(option_struct *op, data_fmt *data, brlist *source)
{
long cross;
double numactive;
boolean found;

numactive = randum() * source->numnewactives;
found = FALSE;

if (source->nfsite < source->ofsite && source->nlsite > source->olsite) {
   for(cross = source->nfsite; cross < source->ofsite && !found; cross++) {
      if (op->datatype == 'n') numactive--;
      else numactive -= getdata_space(op,data,cross);
      if (numactive < 0) {found = TRUE; break;}
   }
   if (!found)
      for(cross = source->olsite; cross < source->nlsite && !found; cross++) {
         if (op->datatype == 'n') numactive--;
         else numactive -= getdata_space(op,data,cross);
         if (numactive < 0) {found = TRUE; break;}
      }
} else {
   if (source->nlsite > source->olsite && source->nfsite < source->olsite)
      for(cross = source->olsite; cross < source->nlsite && !found; cross++) {
         if (op->datatype == 'n') numactive--;
         else numactive -= getdata_space(op,data,cross);
         if (numactive < 0) {found = TRUE; break;}
      }
   else
      for(cross = source->nfsite; cross < source->nlsite && !found; cross++) {
         if (op->datatype == 'n') numactive--;
         else numactive -= getdata_space(op,data,cross);
         if (numactive < 0) {found = TRUE; break;}
      }
}

return(cross);


} /* pickbrcross */


/*********************************************************************
 * eventprobbr() is an auxiliary function to eventprob that helps it *
 * deal with brlist segments.                                        *
 *                                                                   *
 * Due to brsegment endpoint problems, (what endpoints to include?), *
 * eventprobbr() now also picks both the brlist segment that will    *
 * have an event, if one is desired; and what the crossover point    *
 * will be.  Segment and point info are passed in the bseg struct.   */
double eventprobbr(option_struct *op, data_fmt *data, brlist **bigbr,
   tlist *tint, double pc, double pr, double *pb, bseg *brs)
{
double ttyme, btyme, tyme, offset, brlength, stoptyme, bsitesum;
brlist *br, *brsucc, *finish;

bsitesum = 0.0;
offset = 0.0;
stoptyme = tint->age;
tyme = BOGUSTREETYME;
brlength = BOGUSTREETYME - 1.0;
btyme = tint->eventnode->tyme;
finish = NULL;

while (tyme > brlength) {
   if (!(*bigbr)) {
      ttyme = btyme;
      btyme = stoptyme;
      (*pb) = 0.0;
   } else {
      if ((*bigbr)->starttyme > btyme) {
         ttyme = btyme;
         btyme = (((*bigbr)->starttyme < stoptyme) ?
                    (*bigbr)->starttyme : stoptyme);
         (*pb) = 0.0;
         finish = (*bigbr);
      } else {
         ttyme = (*bigbr)->starttyme;
         btyme = DBL_MAX;
         bsitesum = 0;
         for(br = (*bigbr); br->starttyme == ttyme; br = br->succ) {
            if (br->endtyme < btyme) btyme = br->endtyme;
            bsitesum += br->numnewactives;
            if (br->succ == NULL) {br = NULL; break;}
         }
         if (btyme > stoptyme) btyme = stoptyme;
/* now also find anything that starts before the shortest entry ends;
   finish will end up pointing to the first entry we don't want to do
   anything with.  Assumption of contingous brlist entries contingent
   upon the list being sorted by starttymes, least to greatest! */
         finish = br;
         if (br) if (br->starttyme < btyme) btyme = br->starttyme;
         *pb = bsitesum*rec0;
      }
   }

   if (*pb) {
      brs->target = pickbrtarget((*bigbr),bsitesum);
      brs->cross = pickbrcross(op,data,brs->target);
   }

   /* when? */
   if (pc || pr || *pb) tyme = -log(randum())/(pc + pr + (*pb));
   else tyme = btyme + 0.1;

   brlength = btyme - ttyme;
   if (tyme > brlength) {
      if (*bigbr)
         for(br = (*bigbr); br != finish; br = brsucc) {
            brsucc = br->succ;
            br->starttyme = btyme;
            if (br->starttyme == br->endtyme) br_remove(bigbr,br);
         }
      offset += btyme - ttyme;
   }

   if (btyme == stoptyme) break;

}

return(tyme + offset);

} /* eventprobbr */


/**********************************************************************
 * eventprob() calculates both what kind of event is the "next" event *
 * and what the time is to that event.                                */
double eventprob(option_struct *op, data_fmt *data, linlist *alines,
   tlist *t, brlist **bigbr, char *event, bseg *brs)
{
long numilines, numalines;
double pc, pr, pb, tester, tyme, asitesum;
linlist *u;

if (!alines && (!*bigbr)) {(*event) = '0'; return(0.0);}

/* number of inactive lineages */
numilines = countinactive(op,t,alines);

/* count the active lineages and sum active sites */
numalines = 0;
asitesum = 0.0;
for (u = alines; u != NULL; u = u->succ) {
   numalines++;
   asitesum += u->activesites;
}

pc = (numalines*(numalines-1) + 2.0*numalines*numilines) / theta0;
pr = asitesum*rec0;

if (!(*bigbr)) {pb = 0.0; tyme = -log(randum())/(pc + pr);}
else tyme = eventprobbr(op,data,bigbr,t,pc,pr,&pb,brs);

/* what kind of event? */
if (!pc && !pr && !pb) {(*event) = '0'; return(tyme);}
tester = randum();
if (tester <= pr/(pc+pr+pb)) (*event) = 'r';
else if (tester <= (pc+pr)/(pc+pr+pb)) (*event) = 'c';
else (*event) = 'b';

return(tyme);

} /* eventprob */


/******************************************************************
 * traverse_flagbelow sets the "updated" field of the passed node *
 * and all nodes rootwards of it to FALSE.                        */
void traverse_flagbelow(tree *tr, node *p)
{
tlist *t;

for(t = gettymenode(tr,p->number); t != NULL; t = t->succ) {
   t->eventnode->updated = FALSE;
   if (!istip(t->eventnode)) {
      t->eventnode->next->updated = FALSE;
      t->eventnode->next->next->updated = FALSE;
   }
}

} /* traverse_flagbelow */


/*********************************************************************
 * flag_below marks the node "p" iff p is not a tip, and in any case *
 *    marks all nodes "rootwards" of p.                              */
void flag_below(node *p)
{

if (istip(p) && p != curtree->root) flag_below(p->back);

if(!istip(p)) {
   while (p->top) p = p->next;
   p->updated = FALSE;
   p->next->updated = FALSE;
   if(p->next->top) flag_below(p->next->back);
   p->next->next->updated = FALSE;
   if(p->next->next->top) flag_below(p->next->next->back);
}
} /* flag_below */


/*********************************************************************
 * traverse_unflag() sets the "updated" field of the passed node and *
 * all nodes rootwards of it to TRUE.                                */
void traverse_unflag(tree *tr, node *p)
{
tlist *t;

for(t = gettymenode(tr,p->number); t != NULL; t = t->succ) {
   t->eventnode->updated = TRUE;
   if (!istip(t->eventnode)) {
      t->eventnode->next->updated = TRUE;
      t->eventnode->next->next->updated = TRUE;
   }
}

} /* traverse_unflag */


void unflag(node *p)
/* remove flags.  call with curtree->root->back. */
{
   if(!istip(p)) {
     if(!p->next->top) unflag(p->next->back);
     if(!p->next->next->top) unflag(p->next->next->back);
     p->updated = TRUE;
     p->next->updated = TRUE;
     p->next->next->updated = TRUE;
   }
} /* unflag */

node *pickactive(long numalin, linlist **lineages, node *prohibited)
{
  long daughter;
  linlist *u;

  do {
     daughter = (long)(randum() * numalin) + 1;
     u = *lineages;
     findlin(&u,daughter);
  } while (u->branchtop == prohibited);

  return(u->branchtop);

} /* pickactive */

node *pickinactive(option_struct *op, tlist *t, node *prohibited,
   linlist *alines)
{
node *daughter;

do {
   daughter = t->branchlist[(long)(randum() * t->numbranch)];
} while (daughter == prohibited || isactive(daughter,alines) ||
         (isrecomb(daughter) && isdead(op,daughter)));

return(daughter);

} /* pickinactive */

void coalesce(option_struct *op, data_fmt *data, tlist *t,
   linlist **lineages, long *lines, double tyme, boolean **nodenumber,
   long *numnodes, boolean rootdropping, brlist **br)
{
  long numalin, numilin, *active1, *active2;
  double Pbactive, activelinks;
  node *p, *daughter1, *daughter2;
  boolean bothactive, rootpick;
  dnadata *dna;

  dna = data->dnaptr;

/* find the number of active (numalin) and inactive (numilin) lineages,
   when rootdropping, then the last inactive lineage (the root)
   becomes active; leaving no inactive lineages. */
  numalin = *lines;
  if (!rootdropping) numilin = countinactive(op,t,*lineages);
  else numilin = 0;


  /* pick the 2 lineages to coalesce, daughter1 and daughter2
     Pbactive = the chance that they are both active */
  Pbactive = (numalin)*(numalin-1) /
  ((numalin)*(numalin-1) + 2.0*numalin*numilin);

  if (Pbactive >= randum()) {
     daughter1 = pickactive(numalin,lineages,NULL);
     daughter2 = pickactive(numalin,lineages,daughter1);
     bothactive = TRUE;
  } else {
     daughter1 = pickactive(numalin,lineages,NULL);
     daughter2 = pickinactive(op,t,daughter1,*lineages);
     bothactive = FALSE;
  }

  rootpick = FALSE;
  if (daughter1 == curtree->root->back ||
      daughter2 == curtree->root->back)
     rootpick = TRUE;

  /* coalesce the two partners */
  newnode(&p);
  p->number = setnodenumber(nodenumber,numnodes);
  p->next->number = p->number;
  p->next->next->number = p->number;
  p->top = FALSE;
  free_x(op,p);
  p->next->top = FALSE;
  free_x(op,p->next);
  p->next->next->top = TRUE;
  allocate_x(op,data,p->next->next);
  if (op->map) allocate_z(op,data,p->next->next);
  ranges_Malloc(p->next->next,TRUE,0L);
  p->tyme = tyme;
  p->next->tyme = tyme;
  p->next->next->tyme = tyme;
  p->updated = FALSE;
  p->next->updated = FALSE;
  p->next->next->updated = FALSE;
  p->type = p->next->type = p->next->next->type = 'c';

  curtree->numcoals += 1;
  curtree->nodep[p->number] = p->next->next;

  if(!bothactive) {
    hookup(p->next->next,daughter2->back);
    fixlength(op,data,p->next->next); /* but only if not bothactive! */
  }
  hookup(p,daughter1);
  fixlength(op,data,p);
  hookup(p->next,daughter2);
  fixlength(op,data,p->next);
  if (rootpick) {
     hookup(p->next->next,curtree->root);
     fixlength(op,data,p->next->next);
  }

/* now set the ranges for the new node */
  contrib(op,data,p->next->next,&p->next->next->ranges);

/* update the tymelist */
  /*readdactive(t->succ,daughter1);*/
  /*if (bothactive) readdactive(t->succ,daughter2);*/
  rembranch(curtree,t->succ,daughter1);
  if (bothactive) {
     rembranch(curtree,t->succ,daughter2);
  }
  insertaftertymelist(t,p);
  if (bothactive) {
     extbranch(curtree,t->succ,p->next->next);
  }

/* now set the coal array for the new node */
  if (op->fc) fix_coal(op,data,curtree,t->succ,p);

  sublinlist(lineages,daughter1);
  if (bothactive) {
    if (op->fc) activelinks = count_activefc(op,data,p->next->next);
    else activelinks = count_active(op,data,p->next->next);
    addlinlist(lineages,p->next->next, activelinks);
    sublinlist(lineages,daughter2);
  }
  (*lines)--;

  if(!bothactive) {
  /* rename the changed branchtop segments */
     drop_renamebrlist(br,daughter2,p->next->next);
  /* then fix the rest of the ranges and all brlist segments */
     if (op->fc) {
        active1 = p->next->next->coal;
        active2 = daughter2->coal;
     } else {
        active1 = p->next->next->ranges;
        active2 = daughter2->ranges;
     }
     if(!sameranges(active1,active2))
        addbrlist(op,data,br,p->next->next,active2,
           active1,NULL,p->next->next->tyme,
           p->next->next->back->tyme,FALSE);
     drop_brfix(op,data,curtree,t->succ,br);
     consolidate_brlist(op,data,curtree,br);
  }

} /* coalesce */


/***********************************************************************
 * init_numnewactive() sets the numnewactives field for an entire      *
 * brlist, removing any brlist segments whose tyme has passed.         */
void init_numnewactive(option_struct *op, data_fmt *data, tree *tr,
   brlist **bigbr, node *branchtop)
{
brlist *br, *brsucc;

for(br = (*bigbr); br != NULL; br = brsucc) {
   brsucc = br->succ;
   if (branchtop->tyme >= br->endtyme) {br_remove(bigbr,br); continue;}
   if (branchtop->tyme > br->starttyme) {
      br->starttyme = branchtop->tyme;
      if (br->starttyme >= br->endtyme) {br_remove(bigbr,br); continue;}
   }
   br->numnewactives = count_activebr(op,data,br);
   }

} /* init_numnewactive */


/**********************************************************************
 * pickbr_recomb() picks the proper segment to use for re-droping the *
 * new lineage, given the tyme at which the drop must occur.  It then *
 * removes/truncates appropiately all non-chosen segments in the list.*
 * It also picks which site will be the point of crossover.           */
brlist *pickbr_recomb(dnadata *dna, brlist **bigbr, double droptyme, long *cross)
{
brlist *br, *brsucc, *start, *picked;
long numbr;
double pick, numactive;
boolean found;

/* initializations to make lint happy */
start = NULL;
numbr = 0;
numactive = 0.0;

for(br = (*bigbr); br != NULL; br = brsucc) {
   brsucc = br->succ;
   if(!isintyme(br,droptyme)) {br_remove(bigbr,br); continue;}
   numbr = 0;
   numactive = 0.0;
   for(start = br; isintyme(br,droptyme); br = br->succ) {
      numbr++;
      numactive += br->numnewactives;
      if (br->succ == NULL) break;
   }
   break;
}
if (!numbr)
   fprintf(ERRFILE,"ERROR:pickbr_recomb failed to do any picking!\n");

if (!numactive)
   fprintf(ERRFILE,"ERROR:pickbr_recomb has no active links!\n");


pick = randum() * numactive - start->numnewactives;
for(br = start; pick > 0; br = br->succ) pick -= br->succ->numnewactives;
picked = br;

numactive = pick+br->numnewactives;
found = FALSE;

if (br->nfsite < br->ofsite && br->nlsite > br->olsite) {
   for(*cross = br->nfsite; *cross < br->ofsite && !found; (*cross)++) {
      numactive -= dna->sspace[population][locus][*cross];
      if (numactive < 0) {found = TRUE; break;}
   }
   if (!found)
      for(*cross = br->olsite; *cross < br->nlsite && !found;
          (*cross)++) {
         numactive -= dna->sspace[population][locus][*cross];
         if (numactive < 0) {found = TRUE; break;}
      }
} else {
   if (br->nlsite > br->olsite && br->nfsite < br->olsite)
      for(*cross = br->olsite; *cross < br->nlsite && !found;
          (*cross)++) {
         numactive -= dna->sspace[population][locus][*cross];
         if (numactive < 0) {found = TRUE; break;}
      }
   else
      for(*cross = br->nfsite; *cross < br->nlsite && !found;
          (*cross)++) {
         numactive -= dna->sspace[population][locus][*cross];
         if (numactive < 0) {found = TRUE; break;}
      }
}

for(pick = 0, br = start; pick < numbr; pick++, br = brsucc) {
   brsucc = br->succ;
   if (br != picked) {
      if(droptyme < br->endtyme) br->starttyme = droptyme;
      else br_remove(bigbr,br);
   }
}

return(picked);

} /* pickbr_recomb */


/*********************************************************************
 * boolean brlistrecomb() inserts a new recombination into a partial *
 * tree, using a brlist entry as the basis for the recombination.    */
boolean brlistrecomb(option_struct *op, data_fmt *data, tlist *t,
   linlist **lineages, long *lines, brlist **bigbr, double tyme,
   long *siteptr, boolean **nodenumber, long *numnodes, bseg *brs)
{
long target, numsites;
brlist *br;
node *d, *p;
dnadata *dna;

dna = data->dnaptr;
numsites = countsites(op,data);

br = brs->target;
target = brs->cross;
d = br->branchtop;

/* create recombination node */
newnode(&p);
p->number = setnodenumber(nodenumber,numnodes);
p->next->number = p->number;
p->next->next->number = p->number;
p->top = FALSE;
free_x(op,p);
p->next->top = TRUE;
allocate_x(op,data,p->next);
if(op->map)allocate_z(op,data,p->next);
ranges_Malloc(p->next,TRUE,0L);
p->next->next->top = TRUE;
allocate_x(op,data,p->next->next);
if(op->map)allocate_z(op,data,p->next->next);
ranges_Malloc(p->next->next,TRUE,0L);
p->tyme = tyme;
p->next->tyme = tyme;
p->next->next->tyme = tyme;
p->updated = FALSE;
p->next->updated = FALSE;
p->next->next->updated = FALSE;
p->type = p->next->type = p->next->next->type = 'r';

curtree->numrecombs += 1;
curtree->nodep[p->number] = p->next;

hookup(p->next,d->back);
fixlength(op,data,p->next);
hookup(p,d);
fixlength(op,data,p);
fix_rec_ranges(op,data,p);

/* partition the sites */
if(randum()>0.5) {
   p->next->recstart = 0;
   p->next->recend = target;
   p->next->next->recstart = target + 1;
   p->next->next->recend = numsites-1;
} else {
   p->next->recstart = target + 1;
   p->next->recend = numsites-1;
   p->next->next->recstart = 0;
   p->next->next->recend = target;
}

/* now update the ranges for the new node */
  contrib(op,data,p->next,&p->next->ranges);
  contrib(op,data,p->next->next,&p->next->next->ranges);

/* update the tymelist */
insertaftertymelist(t,p);
extbranch(curtree,t->succ,p->next->next);

/* now set the coal arrays */
if (op->fc) fix_coal(op,data,curtree,t->succ,p);

/* now add brlist stuff for the new node's branch segment */
if (op->fc) addbrlist(op,data,bigbr,p->next,d->coal,p->next->coal,
               NULL,p->next->tyme,p->next->back->tyme,FALSE);
else addbrlist(op,data,bigbr,p->next,d->ranges,p->next->ranges,
        NULL,p->next->tyme,p->next->back->tyme,FALSE);


/* fix alias entries */
if (op->datatype == 'n' || op->datatype == 's')
   edit_alias(op,data,siteptr,target);

/* add new entry to linlist */
if (op->fc) addlinlist(lineages,p->next->next,
   count_activefc(op,data,p->next->next));
else addlinlist(lineages,p->next->next,count_active(op,data,p->next->next));
(*lines)++;

if (tyme < br->endtyme) br->starttyme = tyme;
drop_renamebrlist(bigbr,br->branchtop,p->next);
drop_brfix(op,data,curtree,t->succ,bigbr);
consolidate_brlist(op,data,curtree,bigbr);

if (curtree->numrecombs > RECOMB_MAX) return(FALSE);
return(TRUE);

} /* brlistrecomb */


/********************************************************************
 * choosepsite() chooses the sequence site that a new recombination *
 * will start at.                                                   */
long choosepsite(option_struct *op, data_fmt *data, node *p,
   double targetlink)
{
long target;


for(target = p->ranges[1]; target < p->ranges[2*p->ranges[0]]; target++) {
   if (op->datatype == 'n') targetlink--;
   else targetlink -= getdata_space(op,data,target);
   if (targetlink <= 0) break;
}

if (target == p->ranges[2*p->ranges[0]]) {
   fprintf(ERRFILE,"ERROR:choosepsite: failure to choose %ld %ld\n",indecks,
      apps);
}

return(target);

} /* choosepsite */


/**********************************************************************
 * choosepsitefc() chooses the sequence site that a new recombination *
 * will start at.                                                     */
long choosepsitefc(option_struct *op, data_fmt *data, node *p,
   double targetlink)
{
long target;

for(target = p->coal[1]; target < p->coal[2*p->coal[0]]; target++) {
   if (op->datatype == 'n') targetlink--;
   else targetlink -= getdata_space(op,data,target);
   if (targetlink <= 0) break;
}

if (target == p->coal[2*p->coal[0]]) {
   fprintf(ERRFILE,"ERROR:choosepsitefc: failure to choose %ld %ld\n",indecks,
      apps);
}

return(target);

} /* choosepsite */


boolean recomb(option_struct *op, data_fmt *data, tlist *t,
  linlist **lineages, long *lines, double tyme, long *siteptr,
  boolean **nodenumber, long *numnodes)
{
  long target, numsites;
  double ntarget, numactive;
  linlist *u;
  node *d, *p;
  boolean rootpick;
  dnadata *dna;

  dna = data->dnaptr;
  numsites = countsites(op,data);

  u = *lineages;
  /* find splitting lineage, which must be active */
  /* we first count up all active sites on active lineages -- */
  numactive = 0;
  do {
    numactive += u->activesites;
    u = u->succ;
  } while (u != NULL);
  /* -- then we choose a lineage proportional to the active sites */
  ntarget = randum()*numactive;
  numactive = 0;
  u = *lineages;
  do {
    numactive += u->activesites;
    if (u->activesites != 0 && numactive >= ntarget) break;
    u = u->succ;
  } while(1);
  d = u->branchtop;

  rootpick = FALSE;
  if (d == curtree->root->back) rootpick = TRUE;

/* create recombination node */
  newnode(&p);
  p->number = setnodenumber(nodenumber,numnodes);
  p->next->number = p->number;
  p->next->next->number = p->number;
  p->top = FALSE;
  free_x(op,p);
  p->next->top = TRUE;
  allocate_x(op,data,p->next);
  if(op->map)allocate_z(op,data,p->next);
  ranges_Malloc(p->next,TRUE,0L);
  p->next->next->top = TRUE;
  allocate_x(op,data,p->next->next);
  if(op->map)allocate_z(op,data,p->next->next);
  ranges_Malloc(p->next->next,TRUE,0L);
  p->tyme = tyme;
  p->next->tyme = tyme;
  p->next->next->tyme = tyme;
  p->updated = FALSE;
  p->next->updated = FALSE;
  p->next->next->updated = FALSE;
  p->type = p->next->type = p->next->next->type = 'r';

  curtree->numrecombs += 1;
  curtree->nodep[p->number] = p->next;

  hookup(p,d);
  fixlength(op,data,p);
  if (rootpick) {
     hookup(p->next,curtree->root);
     fixlength(op,data,p->next);
  }

/* partition the sites */
  numactive = randum()*u->activesites;
  if (op->fc) target = choosepsitefc(op,data,d,numactive);
  else target = choosepsite(op,data,d,numactive);
  if(randum()>0.5) {
    p->next->recstart = 0;
    p->next->recend = target;
    p->next->next->recstart = target+1;
    p->next->next->recend = numsites-1;
  } else {
    p->next->recstart = target+1;
    p->next->recend = numsites-1;
    p->next->next->recstart = 0;
    p->next->next->recend = target;
  }

/* set the ranges fields for the 2 rootwards branches */
  contrib(op,data,p->next,&p->next->ranges);
  contrib(op,data,p->next->next,&p->next->next->ranges);

  /* update the tymelist */
  /*readdactive(t->succ,d);*/
  rembranch(curtree,t->succ,d);
  insertaftertymelist(t,p);
  extbranch(curtree,t->succ,p->next->next);
  extbranch(curtree,t->succ,p->next);

/* set the coal fields for the 2 rootwards branches */
  if (op->fc) fix_coal(op,data,curtree,t->succ,p);

/* fix alias entries */
  if (op->datatype == 'n' || op->datatype == 's')
     edit_alias(op,data,siteptr,target+1);

/* remove old entry from linlist */
  sublinlist(lineages,d);

/* add two new entries to linlist */
  if (op->fc) {
     addlinlist(lineages,p->next,count_activefc(op,data,p->next));
     addlinlist(lineages,p->next->next,count_activefc(op,data,p->next->next));
  } else {
     addlinlist(lineages,p->next,count_active(op,data,p->next));
     addlinlist(lineages,p->next->next,count_active(op,data,p->next->next));
  }
  (*lines)++;

  if (curtree->numrecombs > RECOMB_MAX) return(FALSE);
  return(TRUE);
} /* recomb */


boolean find_rootloop(tlist *start, tlist **found)
{
tlist *t;

for (t = start; t != NULL; t = t->succ)
   if (t->numbranch == 1)
      if (t->eventnode != curtree->root->back) {
         (*found) = t;
         return(TRUE);
      }

(*found) = NULL;
return(FALSE);

} /* find_rootloop */


/******************************************************************
 * free_from_tymelist truncates the tymelist starting at argument *
 * "start", then cleans up the revenant time slices.              */
void free_from_tymelist(tlist *start)
{
tlist *t, *temp;

start->prev->succ = NULL;

for (t = start; t != NULL; t = temp) {
   temp = t->succ;
   if (isrecomb(t->eventnode)) curtree->numrecombs--;
   else curtree->numcoals--;
   fix_treenodep(curtree,t->eventnode->number);
   freenode(t->eventnode);
   freetymenode(t);
}

} /* free_from_tymelist */


/***********************************************************************
 * update_rangecoal() updates all ranges and coal arrays in the passed *
 * in tree.  It assumes that both the tree and tymelist are correct.   */
void update_rangecoal(option_struct *op, data_fmt *data, tree *tr)
{
tlist *t;
node *p;

for(t = tr->tymelist->succ; t != NULL; t = t->succ) {
   p = findunique(t->eventnode);
   if (isrecomb(p)) {
      fix_rec_ranges(op,data,p);
   } else {
      fix_coal_ranges(op,data,p);
   }
   if (op->fc) fix_coal(op,data,tr,t,p);
}

} /* update_rangecoal */


void remove_rootloop(option_struct *op, data_fmt *data, long *sp)
{
tlist *t;
dnadata *dna;

dna = data->dnaptr;

if (curtree->numrecombs == 0) return;

if (find_rootloop(curtree->tymelist,&t)) {
   hookup(t->eventnode,curtree->root);
   t->age = ROOTLENGTH + t->eventnode->tyme;
   curtree->root->tyme = t->age;
   curtree->root->back->length = ROOTLENGTH;
   ltov(op,data,curtree->root->back);
   free_from_tymelist(t->succ);
   if (op->datatype == 'n' || op->datatype == 's')
      rebuild_alias(op,data,sp);
}

} /* remove_rootloop */


void remove_deadrecombs(option_struct *op, data_fmt *data, tree *tr)
{
long deadcount;
node *p;
tlist *t;

if (tr->numrecombs < 2) return;

deadcount = 0;
for(t = tr->tymelist; t->succ != NULL; t = t->succ) ;
for(; t != NULL; t = t->prev) {
   p = findunique(t->eventnode);
   if (isrecomb(p)) {
      if (isdead(op,p->next)) {
         remove_pairednodes(op,data,tr,p,p->next);
         deadcount++;
         for(t = tr->tymelist; t->succ != NULL; t = t->succ) ;
      } else {
         if (isdead(op,p->next->next)) {
              remove_pairednodes(op,data,tr,p,p->next->next);
              deadcount++;
              for(t = tr->tymelist; t->succ != NULL; t = t->succ) ;
           }
      }
   }
}

if (deadcount) update_rangecoal(op,data,tr);

} /* remove_deadrecombs */


/*********************************************************************
 * remove_excess_tree() removes "dead" branches from the tree:       *
 *    all rootloops are removed, all nodes and branches after the    *
 *        total number of branches in an interval has reached 1,     *
 *    all branches that no longer contain any "live" sites, those    *
 *        sites that eventually reach a tip.                         *
 *                                                                   *
 * Note: these two procedures destroy the exact correspondence       *
 * the "nodep" array and the tree, renumber_nodes() should be called *
 * aftercalling this proceedure.                                     */
void remove_excess_tree(option_struct *op, data_fmt *data, tree *tr,
   long *siteptr)
{
remove_deadrecombs(op,data,tr);
remove_rootloop(op,data,siteptr);
} /* remove_excess_tree */


long counttymelist(tlist *t)
/* count number of entries in tymelist */
{
  long cntr;

  cntr = 0;
  do {
     cntr++;
     t = t->succ;
  } while (t!=NULL);
  return(cntr);
} /* counttymelist */


boolean rootdrop(option_struct *op, data_fmt *data, tlist *t, double starttyme, linlist **lineages,
   long *lines, long *siteptr, boolean **nodenumber, long *numnodes,
   brlist **br)
{
double basetyme, newtyme;
char eventtype;
boolean not_aborted;
node *p;

/* if necessary, add the root to the active lineages */
if (!foundlin(*lineages,curtree->root->back)) {
   if (op->fc)
      addlinlist(lineages,curtree->root->back,
         count_activefc(op,data,curtree->root->back));
   else
      addlinlist(lineages,curtree->root->back,
         count_active(op,data,curtree->root->back));
   (*lines)++;
}

basetyme = starttyme;
p = curtree->root->back;
not_aborted = TRUE;
delete_brlist(br);

while (*lines > 1) {
   newtyme = basetyme + eventprob(op,data,*lineages,t,br,&eventtype,NULL);
   switch ((int)eventtype) {
      case 'c' :
         coalesce(op,data,t,lineages,lines,newtyme,nodenumber,
            numnodes,TRUE,br);
         break;
      case 'r' :
         not_aborted = recomb(op,data,t,lineages,lines,
            newtyme,siteptr,nodenumber,numnodes);
         break;
      case 'b' :
         fprintf(ERRFILE,"ERROR:rootdrop case b: can't get here\n");
         break;
      case '0' :
         fprintf(ERRFILE,"ERROR:rootdrop case 0: can't get here\n");
         break;
      default :
         fprintf(ERRFILE,"ERROR:rootdrop case: can't get here\n");
         break;
      }
   if (!not_aborted) break;
   t = t->succ;
   basetyme = newtyme;

}


if (not_aborted) traverse_flagbelow(curtree,p);
return(not_aborted);

} /* rootdrop */


node *pickbranch(option_struct *op, data_fmt *data, tree *source,
   tlist **tymeslice, double *offset)
{
long branch, numbranches;
node *p;
tlist *t;

(*offset) = 0.0;
numbranches = countbranches(op,data,source);

branch = (long)(randum() * numbranches) + 1;

if (branch <= getdata_numtips(op,data)) {
   if (tymeslice) (*tymeslice) = source->tymelist;
   return(source->nodep[branch]);
}
branch -= getdata_numtips(op,data);

for (t = source->tymelist->succ; t != NULL; t = t->succ) {
   branch -= (isrecomb(t->eventnode) ? 2 : 1);
   if (branch < 1) {
      (*tymeslice) = t;
      p = findunique(t->eventnode);
      if (isrecomb(p))
         return((randum() < 0.5) ? p->next : p->next->next);
      else
         return(findtop(p));
   }
}

fprintf(ERRFILE,"ERROR:pickbranch--failure to find branch\n");
return(NULL);

} /* pickbranch */


void renumber_nodes(option_struct *op, data_fmt *data, tree *trii)
{
tlist *t;
long nodenumber;

/* since the root is "0", the first non-tip node will be "numseq+1" */
nodenumber = getdata_numtips(op,data) + 1;

for (t = trii->tymelist->succ; t != NULL; t = t->succ) {
   t->eventnode->number = nodenumber;
   t->eventnode->next->number = nodenumber;
   t->eventnode->next->next->number = nodenumber;
   trii->nodep[nodenumber] = t->eventnode;
   nodenumber++;
}

} /* renumber_nodes */

void finishbadtree(option_struct *op, data_fmt *data, tree *tr, linlist **lineages, long *numlins,
   boolean **nodenumber, long *numnodes, boolean rootdropped,
   brlist **brlines)
{
double ntyme;
tlist *t;

/* go to the end of the tymelist */
for(t = tr->tymelist; t->succ != NULL; t = t->succ) ;

while (*numlins > 1) {
   ntyme = t->eventnode->tyme + 0.01;
   coalesce(op,data,t,lineages,numlins,ntyme,nodenumber,numnodes,rootdropped,
      brlines);
   t = t->succ;
}
delete_brlist(brlines);

} /* finishbadtree */


/***********************************************************
 * initbrlist() allocates and initializes a brlist element */
brlist *initbrlist(node *branchtop, double starttyme, double endtyme,
   long numsites)
{
brlist *newbr;

newbr = (brlist *)calloc(1,sizeof(brlist));
newbr->branchtop = branchtop;
newbr->segranges = (int *)calloc(numsites,sizeof(int));
newbr->starttyme = starttyme;
newbr->endtyme = endtyme;
newbr->updated = FALSE;

return(newbr);

} /* initbrlist */


/**************************************************
 * freebrlist() just frees an element of a brlist */
void freebrlist(brlist *br)
{

free(br->segranges);
free(br);

} /* freebrlist */


/*********************************************************
 * delete_brlist() deletes an entire brlist, setting the *
 * base pointer to NULL.                                 */
void delete_brlist(brlist **bigbr)
{
brlist *br, *brsucc;

if (!bigbr) return;

for(br = (*bigbr); br != NULL; br = brsucc) {
   brsucc = br->succ;
   br_remove(bigbr,br);
}

(*bigbr) = NULL;

} /* delete_brlist */


/******************************************************************
 * hookup_brlist() mechanically adds a new element into a brlist. *
 * New elements are inserted so that the list stays time-sorted   *
 * (least to greatest).                                           */
void hookup_brlist(brlist **bigbr, brlist *newbr)
{
brlist *br;

if (!*bigbr) { /* new list being created */
   (*bigbr) = newbr;
   newbr->succ = newbr->prev = NULL;
} else {
   br = getbrlistbytyme((*bigbr),newbr->starttyme);
   if (!br) { /* adding to end of list */
      for (br = (*bigbr); br->succ != NULL; br = br->succ)
         ;
      br->succ = newbr;
      newbr->prev = br;
      newbr->succ = NULL;
   } else {
      newbr->succ = br;
      newbr->prev = br->prev;
      if (br->prev) br->prev->succ = newbr;
      br->prev = newbr;
      if (*bigbr == br) (*bigbr) = newbr;
   }
}

} /* hookup_brlist */


/*******************************************************************
 * set_sitesbr() sets the "sites" information for a brlist segment */
void set_sitesbr(option_struct *op, data_fmt *data, brlist *br)
{
long i, numsites;
boolean done1, done2;

numsites = countsites(op,data);

done1 = FALSE;
done2 = FALSE;

br->nfsite = numsites;
br->ofsite = numsites;

for (i = 0; i < numsites; i++) {
   if (!done1)
      if (br->segranges[i] >= 1) {
         br->nfsite = i;
         done1 = TRUE;
         }
   if (!done2)
      if (br->segranges[i] == -1 ||
          br->segranges[i] == 2) {
         br->ofsite = i;
         done2 = TRUE;
         }
   if (done1 && done2) break;
}

done1 = FALSE;
done2 = FALSE;

br->nlsite = br->nfsite - 1;
br->olsite = br->ofsite - 1;

for (i = numsites-1; i >= 0; i--) {
   if (!done1)
      if (br->segranges[i] >= 1) {
          br->nlsite = i;
          done1 = TRUE;
          }
   if (!done2)
      if (br->segranges[i] == -1 ||
          br->segranges[i] == 2) {
          br->olsite = i;
          done2 = TRUE;
          }
   if (done1 && done2) break;
}

} /* set_sitesbr */


/***************************************************************
 * count_activebr() returns the number of active links present *
 * in a brlist segment.                                        */
double count_activebr(option_struct *op, data_fmt *data, brlist *br)
{

if (!br->updated) {
   set_sitesbr(op,data,br);
   br->updated = TRUE;
   }

if (br->nfsite == countsites(op,data)) return(0L);

if (br->olsite <= br->nfsite || br->ofsite >= br->nlsite) {
   return(getnumlinks(op,data,br->nfsite,br->nlsite));
}

if (br->nfsite < br->ofsite && br->nlsite > br->olsite) {
   return(getnumlinks(op,data,br->nfsite,br->ofsite) +
       getnumlinks(op,data,br->olsite,br->nlsite));
}

if (br->nfsite < br->ofsite) {
   return(getnumlinks(op,data,br->nfsite,br->ofsite));
}

if (br->nlsite > br->olsite) {
   return(getnumlinks(op,data,br->olsite,br->nlsite));
}

return(0.0);


} /* count_activebr */


/***************************************************************
 * isintyme() checks to see if a given tyme is within the range *
 * covered by a particular brlist entry.                       */
boolean isintyme(brlist *br, double target)
{

return((target >= br->starttyme) && (target <= br->endtyme) ?
          TRUE : FALSE);

} /* isintyme */


/******************************************************************
 * getbrlistbystart() returns a pointer to the brlist entry which *
 * contains the nodelet; returning NULL if nothing is found.      */
brlist *getbrlistbystart(brlist *start, node *target, double starttyme)
{
brlist *br;
long match;

for(br = start; br != NULL; br = br->succ) {
   match = sametreebranch(br,target,starttyme,BOGUSTREETYME);
   if (match == 1 || match == 2) return(br);
}

return(NULL);

} /* getbrlistbystart */


/****************************************************************
 * getbrlistbyend() returns a pointer to the brlist entry which *
 * contains the nodelet; returning NULL if nothing is found.    */
brlist *getbrlistbyend(brlist *start, node *target, double endtyme)
{
brlist *br;
long match;

for(br = start; br != NULL; br = br->succ) {
   match = sametreebranch(br,target,BOGUSTREETYME,endtyme);
   if (match == -1 || match == 2) return(br);
}

return(NULL);

} /* getbrlistbyend */


/********************************************************************
 * getbrlistbytyme() finds the first element in a brlist that has a *
 * starttyme >= the searchtyme, returning NULL if nothing is found. */
brlist *getbrlistbytyme(brlist *start, double searchtyme)
{
brlist *br;

for(br = start; br != NULL; br = br->succ)
   if (br->starttyme >= searchtyme) return(br);

return(NULL);

} /* getbrlistbytyme */


/*******************************************************
 * updatesegranges() updates the segranges of a brlist *
 * given the passed arrays.                            *
 *                                                     *
 * segranges is 0 if site has remained dead            *
 *           is +1 for becoming live                   *
 *           is -1 for becoming dead                   *
 *           is +2 for remaining alive                 */
void updatesegranges(option_struct *op, data_fmt *data, brlist *br,
   int *oldsegranges, long *oldranges, long *newranges)
{
long i, newi, newstart, newend, numsites;

long oldind, newind, cstart, cend, csize, oval, nval;
boolean oldlive, newlive;
dnadata *dna;
int *bstart;

br->updated = FALSE;

numsites = countsites(op,data);

if (!oldranges && !newranges && oldsegranges) {
   memcpy(br->segranges,oldsegranges,numsites*sizeof(int));
   return;
}

dna = data->dnaptr;
newi = 1;
newstart = newranges[newi];
if (newstart == FLAGLONG) {
   newstart = numsites;
   newend = numsites;
} else newend = newranges[newi+1];
if (oldranges) {
if (oldranges[0]) {
   oldlive = FALSE;
   newlive = FALSE;
   oldind = newind = 1;
   oval = oldranges[oldind];
   nval = ((newstart != numsites) ? newranges[newind] : numsites);
   cstart = 0;
   cend = ((nval < oval) ? nval : oval);
   while (1) {
      csize = cend - cstart;
      bstart = &(br->segranges[cstart]);
      if (csize) {
         if (oldlive && newlive)
           memcpy(bstart,dna->segranges2,csize*sizeof(int));
         else if (!oldlive && newlive)
           memcpy(bstart,dna->segranges1,csize*sizeof(int));
         else if (!oldlive && !newlive)
           memcpy(bstart,dna->segranges0,csize*sizeof(int));
         else if (oldlive && !newlive)
           memcpy(bstart,dna->segranges3,csize*sizeof(int));
      }
      if (cend == numsites) break;
      if (oval == cend) {
         oldind++;
         oldlive = !oldlive;
         if (oldlive) oval = oldranges[oldind] + 1;
         else oval = oldranges[oldind];
         if (oval == FLAGLONG)
            oval = numsites;
      }
      if (nval == cend) {
         newind++;
         newlive = !newlive;
         if (newlive) nval = newranges[newind] + 1;
         else nval = newranges[newind];
         if (nval == FLAGLONG)
            nval = numsites;
      }
      cstart = cend;
      cend = (nval < oval) ? nval : oval;
   }
}
} else {
   for(i = 0; i < numsites; i++) {
      if (i > newend) {
         newi+=2;
         if (newranges[newi] == FLAGLONG) {
            newstart = numsites;
            newend = numsites;
         } else {
            newstart = newranges[newi];
            newend = newranges[newi+1];
         }
      }
      if (i >= newstart && i <= newend) {
         if (oldsegranges[i] == 0) br->segranges[i] = 1;
         else if (oldsegranges[i] == -1) br->segranges[i] = 2;
         else br->segranges[i] = oldsegranges[i];
     } else {
         if (oldsegranges[i] == 2) br->segranges[i] = -1;
         else if (oldsegranges[i] == 1) br->segranges[i] = 0;
         else br->segranges[i] = oldsegranges[i];
     }

   }
}

} /* updatesegranges */


/************************************************************
 * segranges_to_ranges() converts from segranges to ranges. *
 *                                                          *
 * This code may be wrong--it converts segranges to current *
 * range info, not the original?                            */
void segranges_to_ranges(int *segranges, long **ranges, long numsites)
{
long i, newstart, newend;
boolean live;

/* initialization for lint */
newend = 0;

for(i = 0, newstart = FLAGLONG, live = FALSE; i < numsites; i++) {
   if((!live && (segranges[i] == -1 || segranges == 0)) ||
     (live && (segranges[i] == 1 || segranges[i] == 2)))
      continue;
   live = !live;
   if(live) {newstart = i; newend = FLAGLONG;}
   if(!live) newend = i-1;
   if(newstart != FLAGLONG && newend != FLAGLONG) {
      addrange(ranges,newstart,newend);
      newstart = FLAGLONG;
   }
}

} /* segranges_to_ranges */


/*********************************************************
 * addbrlist() adds an element to an existing brlist and *
 * depends on both the tree and tymelist being correct   *
 * from the tips to the tyme of the nodelet newtop.      */
void addbrlist(option_struct *op, data_fmt *data, brlist **bigbr,
   node *newtop, long *oldranges, long *newranges, int *segranges,
   double starttyme, double endtyme, boolean building)
{
double oldendtyme;
int *tsegranges, *oldsegs;
brlist *tnewbr, *bnewbr, *newbr;
long numsites;

numsites = countsites(op,data);
if (starttyme == endtyme || !newtop) return;

tnewbr = getbrlistbystart(*bigbr,newtop,starttyme);
bnewbr = getbrlistbyend(*bigbr,newtop,endtyme);

if (!tnewbr && !bnewbr) {
   newbr = initbrlist(newtop,starttyme,endtyme,countsites(op,data));
   if (oldranges && newranges) {
      updatesegranges(op,data,newbr,NULL,oldranges,newranges);
   } else {
      updatesegranges(op,data,newbr,segranges,NULL,NULL);
   }
   hookup_brlist(bigbr,newbr);
   return;
}

/* we have the same start & end */
if (tnewbr && bnewbr) {
   if (tnewbr->endtyme == endtyme) { /* we have the same segment */
      if (newranges) {
         if (building) {
            updatesegranges(op,data,tnewbr,NULL,oldranges,newranges);
         } else {
            updatesegranges(op,data,tnewbr,tnewbr->segranges,NULL,newranges);
         }
      } else {
         updatesegranges(op,data,tnewbr,segranges,NULL,NULL);
      }
   } else {
      if (newranges) {
         for(; ; tnewbr = tnewbr->succ) {
            if (tnewbr->branchtop != bnewbr->branchtop) continue;
            if (building) {
               updatesegranges(op,data,tnewbr,NULL,oldranges,newranges);
            } else {
               updatesegranges(op,data,tnewbr,tnewbr->segranges,NULL,
                  newranges);
            }
         if (tnewbr == bnewbr) break;
         if (tnewbr->endtyme != bnewbr->starttyme)
            addbrlist(op,data,bigbr,newtop,oldranges,newranges,segranges,
               tnewbr->endtyme,bnewbr->starttyme,building);
         }
      }
   }
   return;
}

/* we only have the same start */
if (tnewbr) {
   oldendtyme = tnewbr->endtyme;
   if (oldendtyme > endtyme) {
      tnewbr->starttyme = endtyme;
      addbrlist(op,data,bigbr,newtop,oldranges,newranges,segranges,
         starttyme,endtyme,building);
   } else {
      updatesegranges(op,data,tnewbr,tnewbr->segranges,NULL,newranges);
      addbrlist(op,data,bigbr,newtop,oldranges,newranges,segranges,
         oldendtyme,endtyme,building);
      }
   return;
} else { /* we only have the same end */
   if (starttyme > bnewbr->starttyme) {
      bnewbr->endtyme = starttyme;
      oldsegs = (int *)calloc(numsites,sizeof(int));
      memcpy(oldsegs,bnewbr->segranges,numsites*sizeof(int));
      updatesegranges(op,data,bnewbr,bnewbr->segranges,NULL,newranges);
      tsegranges = (int *)calloc(numsites,sizeof(int));
      memcpy(tsegranges,bnewbr->segranges, numsites*sizeof(int));
      memcpy(bnewbr->segranges,oldsegs,numsites*sizeof(int));
      addbrlist(op,data,bigbr,newtop,NULL,NULL,tsegranges,starttyme,
         endtyme,building);
      free(oldsegs);
      free(tsegranges);
   } else {
      addbrlist(op,data,bigbr,newtop,oldranges,newranges,segranges,
         starttyme,bnewbr->starttyme,building);
      updatesegranges(op,data,bnewbr,bnewbr->segranges,NULL,newranges);
   }
   return;
}

} /* addbrlist */


boolean makeevent(option_struct *op, data_fmt *data, tree *growtree,
   tlist *tymeslice, brlist **brlines, linlist **alines, long *numalines,
   double newtyme, boolean **nodenumber, long *numnodes, long *siteptr,
   char eventtype, boolean rootdropped, bseg *brs)
{
boolean succeeded;

succeeded = TRUE;
switch((int)eventtype) {
   case 'c' :
      coalesce(op,data,tymeslice,alines,numalines,newtyme,nodenumber,
         numnodes,rootdropped,brlines);
      break;
   case 'r' :
      succeeded = recomb(op,data,tymeslice,alines,numalines,newtyme,
         siteptr,nodenumber,numnodes);
      break;
   case 'b' :
      succeeded = brlistrecomb(op,data,tymeslice,alines,numalines,
         brlines,newtyme,siteptr,nodenumber,numnodes,brs);
      break;
   case '0' :
      break;
   default :
      fprintf(ERRFILE,"ERROR:makeevent impossible case: can't get here\n");
      break;
}

return(succeeded);
} /* makeevent */


boolean growlineages(option_struct *op, data_fmt *data, tree *growtree,
   tree *oldtree, tlist *tstart, brlist **brlines, linlist **alines,
   long *numalines, double offsetstart, boolean **nodenumber,
   long *numnodes, long *siteptr)
{
tlist *tymeslice;
double btyme, obtyme, newtyme, offset;
boolean succeeded, rootdropped;
char eventtype;
bseg brs;

tymeslice = tstart;
offset = offsetstart;
rootdropped = FALSE;
succeeded = TRUE;

while(tymeslice != NULL) {

/* just drop the root if you're in the last interval but don't make the
   root active til you're below the oldtree->root->back */
   if (tymeslice->succ == NULL) {
      btyme = curtree->root->back->tyme;
      obtyme = oldtree->root->back->tyme;
      while (obtyme > btyme) {
         init_numnewactive(op,data,growtree,brlines,tymeslice->eventnode);
         newtyme = tymeslice->eventnode->tyme + offset +
            eventprob(op,data,(*alines),tymeslice,brlines,&eventtype,&brs);
         if (newtyme <= obtyme) {
            succeeded = makeevent(op,data,growtree,tymeslice,brlines,alines,
               numalines,newtyme,nodenumber,numnodes,siteptr,eventtype,
               rootdropped,&brs);
            if (!succeeded) break;
            if (succeeded && eventtype != '0') tymeslice = tymeslice->succ;
         }
         btyme = newtyme;
         offset = 0.0;
         if ((*numalines) == 0 && !(*brlines)) break;
      }
      if ((*numalines) == 0 && !(*brlines)) break;
      if (!succeeded) break;
      succeeded = rootdrop(op,data,tymeslice,obtyme,alines,numalines,
         siteptr,nodenumber,numnodes,brlines);
      rootdropped = TRUE;
      break;
   }

   init_numnewactive(op,data,growtree,brlines,tymeslice->eventnode);
   newtyme = tymeslice->eventnode->tyme + offset +
      eventprob(op,data,(*alines),tymeslice,brlines,&eventtype,&brs);
   if (newtyme <= tymeslice->age) {
      succeeded = makeevent(op,data,growtree,tymeslice,brlines,alines,
         numalines,newtyme,nodenumber,numnodes,siteptr,eventtype,
         rootdropped,&brs);
      if ((*numalines) == 0 && !(*brlines)) break;
   }
   if ((*numalines) == 0 && !(*brlines)) break;
   if (!succeeded) break;
   if (tymeslice->succ) tymeslice = tymeslice->succ;
   offset = 0.0;
}

if (!succeeded) {
   numdropped++;
   finishbadtree(op,data,growtree,alines,numalines,nodenumber,numnodes,
      rootdropped,brlines);
}

return(succeeded);

} /* growlineages */

boolean makedrop(option_struct *op, data_fmt *data)
{
long i, numnodes, numalines, nummarkers, *oldsiteptr, *siteptr;
double offset;
linlist *alines;
brlist *brlines;
tlist *tymeslice;
boolean accept_change, *nodenumber;
node *p;
tree *oldtree;

oldsiteptr = NULL;
siteptr = data->siteptr;

/* copy the original tree configuration */
oldtree = copytree(op,data,curtree);
numnodes = 2 * oldtree->numcoals + 2;
nodenumber = (boolean *)calloc(1,numnodes*sizeof(boolean));
nummarkers = getdata_nummarkers(op,data);
for(i = 0; i < numnodes; i++) nodenumber[i] = TRUE;
if(op->datatype == 'n' || op->datatype == 's') {
   oldsiteptr = (long *)calloc(nummarkers,sizeof(long));
   memcpy(oldsiteptr,siteptr,nummarkers*sizeof(long));
}

/* pick a spot to cut */
alines = NULL;
p = pickbranch(op,data,curtree,&tymeslice,&offset);
if (op->fc) addlinlist(&alines,p,count_activefc(op,data,p));
else addlinlist(&alines,p,count_active(op,data,p));
numalines = 1;

brlines = NULL;
tag_futile(p->back);
traverse_flagbelow(curtree,p->back);
remove_futile(op,data,curtree,oldtree,p,nodenumber,&brlines);
accept_change = TRUE;

accept_change = growlineages(op,data,curtree,oldtree,tymeslice,&brlines,
   &alines,&numalines,offset,&nodenumber,&numnodes,siteptr);

if (accept_change) {
   renumber_nodes(op,data,curtree);
   remove_excess_tree(op,data,curtree,siteptr);
   renumber_nodes(op,data,curtree);
#if ALWAYS_ACCEPT
   curtree->likelihood = oldtree->likelihood;
   accept_change = testratio(op,data,oldtree,curtree,'d');
#endif
#if !ALWAYS_ACCEPT
   traverse_flagbelow(curtree,p->back);
   if (op->datatype == 'n' || op->datatype == 's')
      rebuild_alias(op,data,siteptr);
   localeval(op,data,curtree->root->back,FALSE);
   curtree->coalprob = coalprob(op,data,curtree,theta0,rec0);
   accept_change = testratio(op,data,oldtree,curtree,'d');
#endif

}

#if ALWAYS_REJECT
   op->numout[1]++;
   scoretree(op,data,0L);
   accept_change = FALSE;
#endif

free(nodenumber);
freelinlist(alines);

if (accept_change) {
   freetree(op,data,oldtree);
   if (op->datatype == 'n' || op->datatype == 's') free(oldsiteptr);
   return(TRUE);
} else {
   freetree(op,data,curtree);
   curtree = oldtree;
   if (op->datatype == 'n' || op->datatype == 's') {
      memcpy(siteptr,oldsiteptr,nummarkers*sizeof(long));
      free(oldsiteptr);
   }
   return(FALSE);
}

} /* makedrop */


long countrec(tlist *start)
{
long accum;
tlist *t;

accum = 0;
for (t = start; t != NULL; t = t->succ)
   if (isrecomb(t->eventnode)) accum++;

return(accum);

} /* countrec */

node *findrec(tlist *start, long target)
{
long accum;
tlist *t;

accum = 0;
for (t = start; t != NULL; t = t->succ) {
   if (isrecomb(t->eventnode)) accum++;
   if (accum == target) return(findunique(t->eventnode));
}

fprintf(ERRFILE,"ERROR:findrec--unable to find recombination #%ld\n",target);
exit(-1);

} /* findrec */

boolean twiddle(option_struct *op, data_fmt *data)
{
long i, numrecs, target, *oldsiteptr, numnodes, numalines, *siteptr;
double activelinks;
node *p;
boolean accept_change, *nodenumber;
tlist *tstart;
brlist *brlines;
linlist *alines;
tree *oldtree;

numrecs = countrec(curtree->tymelist);
if (!numrecs) return(FALSE);

oldtree = copytree(op,data,curtree);

target = (long)(randum()*numrecs) + 1;
p = findrec(curtree->tymelist,target);
tstart = gettymenode(curtree,p->number);

siteptr = data->siteptr;
oldsiteptr = (long *)calloc(getdata_nummarkers(op,data),sizeof(long));
memcpy(oldsiteptr,siteptr,getdata_nummarkers(op,data)*sizeof(long));

/* now choose the new recombination point, adjusting for the fact that
   the active sites don't necessairly start at the beginning of the
   sequence */
if (op->fc) {
   activelinks = count_activefc(op,data,p);
   activelinks = randum()*activelinks;
   target = choosepsitefc(op,data,p->back,activelinks);
} else {
   activelinks = count_active(op,data,p);
   activelinks = randum()*activelinks;
   target = choosepsite(op,data,p->back,activelinks);
}

if (!p->back->ranges[0])
    fprintf(ERRFILE,"ERROR:twiddle chose a dead branch!\n");

/* 50% chance which branch represents the "first" part of the sites */
if (randum() > 0.5) {
   p->next->recstart = 0;
   p->next->recend = target;
   p->next->next->recstart = target + 1;
   p->next->next->recend = countsites(op,data) - 1;
} else {
   p->next->recstart = target + 1;
   p->next->recend = countsites(op,data) - 1;
   p->next->next->recstart = 0;
   p->next->next->recend = target;
}

/* fix up the ranges and siteptr info, construct the brlist;
   startting drop_brifx_ranges at tstart->prev because a recombination
   should never be the eventnode of the first entry of the tymelist */
brlines = NULL;
drop_brfix(op,data,curtree,tstart->prev,&brlines);
if (op->datatype == 'n' || op->datatype == 's')
   rebuild_alias(op,data,siteptr);

/* now reevaluate the tree for added recombinations */
alines = NULL;
numalines = 0;
numnodes = getdata_numtips(op,data) + counttymelist(curtree->tymelist);
nodenumber = (boolean *)calloc(numnodes,sizeof(boolean));
for(i = 0; i < numnodes; i++) nodenumber[i] = TRUE;

accept_change = growlineages(op,data,curtree,oldtree,tstart,&brlines,
   &alines,&numalines,0.0,&nodenumber,&numnodes,siteptr);

if (accept_change) {
   renumber_nodes(op,data,curtree);
   remove_excess_tree(op,data,curtree,siteptr);
   renumber_nodes(op,data,curtree);
/* root->back's range info may not be set; and must be hand set here. */
   if(curtree->root->back->ranges[0] != 1L)
      addrange(&curtree->root->back->ranges,0L,countsites(op,data)-1);
#if ALWAYS_ACCEPT
   curtree->likelihood = oldtree->likelihood;
   accept_change = testratio(op,data,oldtree,curtree,'t');
#endif
#if !ALWAYS_ACCEPT
   traverse_flagbelow(curtree,p->back);
   if (op->datatype == 'n' || op->datatype == 's')
      rebuild_alias(op,data,siteptr);
   localeval(op,data,curtree->root->back,FALSE);
   curtree->coalprob = coalprob(op,data,curtree,theta0,rec0);
   accept_change = testratio(op,data,oldtree,curtree,'t');
#endif
}

/* if the change was accepted then nothing further needs to be done
   except for cleanup */
if (!accept_change) {
   freetree(op,data,curtree);
   curtree = oldtree;
   memcpy(siteptr,oldsiteptr,getdata_nummarkers(op,data)*sizeof(long));
} else freetree(op,data,oldtree);

free(nodenumber);
free(oldsiteptr);

return(accept_change);

} /* twiddle */


/******************************************************************
 * isprohibited() returns TRUE if the "value" is contained in the *
 * passed "badstuff" and FALSE otherwise.                         */
boolean isprohibited(long value, long numbad, long *badstuff)
{
long bad;

for(bad = 0; bad < numbad; bad++)
   if (value == badstuff[bad]) return(TRUE);

return(FALSE);

} /* isprohibited */


/*******************************************************************
 * isinvariant() returns TRUE if the site is invariant for a given *
 * creature, otherwise FALSE                                       */
boolean isinvariant(option_struct *op, data_fmt *data, long site,
   creature *cr)
{
char **dna;

switch(op->datatype) {
case 'a':
    break;
case 'b':
    break;
case 'm':
    break;
case 'n':
case 's':
   dna = data->dnaptr->seqs[population][locus];
   if (dna[cr->haplotypes[0]->number-1][site] ==
      dna[cr->haplotypes[1]->number-1][site])
      return(TRUE);
   break;
}

return(FALSE);

} /* isinvariant */


/**********************************************
 * choose_flipsite() returns a valid flipsite */
long choose_flipsite(option_struct *op, data_fmt *data, creature *cr,
   long numprohibited, long *prohibited)
{
long choice;

do {
   choice = (long)(randum()*cr->numflipsites);
} while (isprohibited(choice,numprohibited,prohibited));

return(cr->flipsites[choice]);

} /* choose_flipsite */


/********************************************************************
 * datachangesite() flips a passed site between two elements of the *
 * actual data matrix.                                              */
void datachangesite(option_struct *op, data_fmt *data, node *p1, node *p2,
   long site)
{
char tempbase;

tempbase = data->dnaptr->seqs[population][locus][p1->number-1][site];
data->dnaptr->seqs[population][locus][p1->number-1][site] =
   data->dnaptr->seqs[population][locus][p2->number-1][site];
data->dnaptr->seqs[population][locus][p2->number-1][site] = tempbase;

} /* datachangesite */


/*****************************************************************
 * fliphap() flips NUMHAPFLIP sites among the haplotypes of a    *
 * randomly chosen creature of the passed in tree.  The sites to *
 * be flipped are chosen randomly amongst the eligible sites of  *
 * that creature.                                                */
boolean fliphap(option_struct *op, data_fmt *data, tree *tr)
{
long i, j, flipsite, numcreatures, numslice, numcategs, datanumsites,
   *oldsiteptr, *flippedsites, *siteptr;
double ***temp;
creature *cr;
tree *oldtree;
boolean accept_change;

datanumsites = getdata_nummarkers(op,data);
siteptr = data->siteptr;

oldtree = copytree(op,data,tr);
oldsiteptr = (long *)calloc(datanumsites,sizeof(long));
memcpy(oldsiteptr,siteptr,datanumsites*sizeof(long));

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

temp = (double ***)calloc(op->categs,sizeof(double **));
temp[0] = (double **)calloc(op->categs*numslice,sizeof(double *));
for(i = 1; i < op->categs; i++) temp[i] = temp[0] + i*numslice;
temp[0][0] = (double *)calloc(op->categs*numslice*4,sizeof(double));
for(i = 0; i < op->categs; i++)
   for(j = 0; j < numslice; j++)
      temp[i][j] = temp[0][0] + i*numslice*4 + j*4;


flippedsites = (long *)calloc(NUMHAPFLIP,sizeof(long));

numcreatures = getdata_numtips(op,data)/NUMHAPLOTYPES;
do {
   cr = &(tr->creatures[(long)(randum()*numcreatures)]);
} while(cr->numflipsites == 0);

/* debug DEBUG warning WARNING--this code is two haplotype specific */
/* and it is also dna/snp specific!!!!                              */
for(i = 0; i < NUMHAPFLIP; i++) {
   flippedsites[i] = flipsite = choose_flipsite(op,data,cr,i,flippedsites);
   /* first flip the x-arrays */
   memcpy(temp[0][0],cr->haplotypes[0]->x->s[flipsite][0][0],
      numcategs*numslice*4*sizeof(double));
   memcpy(cr->haplotypes[0]->x->s[flipsite][0][0],
      cr->haplotypes[1]->x->s[flipsite][0][0],
      numcategs*numslice*4*sizeof(double));
   memcpy(cr->haplotypes[1]->x->s[flipsite][0][0],temp[0][0],
      numcategs*numslice*4*sizeof(double));
   /* then flip the actual data */
   datachangesite(op,data,cr->haplotypes[0],cr->haplotypes[1],flipsite);
}

/* can just rebuild_alias() because flippable sites can no
   longer be aliased */
if (op->datatype == 'n' || op->datatype == 's')
   rebuild_alias(op,data,siteptr);

for(i = 0; i < NUMHAPLOTYPES; i++)
   traverse_flagbelow(tr,cr->haplotypes[i]->back);
localeval(op,data,tr->root->back,FALSE);
tr->coalprob = coalprob(op,data,tr,theta0,rec0);
accept_change = testratio(op,data,oldtree,tr,'f');
#if ALWAYS_ACCEPT
tr->likelihood = oldtree->likelihood;
accept_change = testratio(op,data,oldtree,tr,'f');
#endif

if (!accept_change) {
   for(i = 0; i < NUMHAPFLIP; i++)
      datachangesite(op,data,cr->haplotypes[1],cr->haplotypes[0],
         flippedsites[i]);
   freetree(op,data,curtree);
   curtree = oldtree;
   memcpy(siteptr,oldsiteptr,datanumsites*sizeof(long));
} else freetree(op,data,oldtree);

free(oldsiteptr);
free(flippedsites);
free(temp[0][0]);
free(temp[0]);
free(temp);

return(accept_change);

} /* fliphap */


/***********************************************************
 * isintymelist() returns TRUE if the nodelet "p" is among *
 * the nodelets in the tymelist; FALSE otherwise.          */
boolean isintymelist(option_struct *op, data_fmt *data,
   tlist *start, node *p)
{
tlist *t;

for(t = start; t != NULL; t = t->succ) {
   if (p == t->eventnode) return(TRUE);
   if (p->next) {
      if (p->next == t->eventnode) return(TRUE);
      if (p->next->next)
         if (p->next->next == t->eventnode) return(TRUE);
   }
}

return(FALSE);

} /* isintymelist */


/**************************************************************
 * flipdrop_prunebr() removes from the brlist the bad entries *
 * created by serial-calling remove_futile().                 */
void flipdrop_prunebr(option_struct *op, data_fmt *data, tree *tr,
   brlist **brlines)
{
brlist *br, *brsucc;

for(br = (*brlines); br != NULL; br = brsucc) {
   brsucc = br->succ;
   if (isintymelist(op,data,tr->tymelist,br->branchtop)) continue;
   br_remove(brlines,br);
}

} /* flipdrop_prunebr */


/*****************************************************************
 * flipdrop() flips NUMHAPFLIP sites among the haplotypes of a   *
 * randomly chosen creature of the passed in tree, and then      *
 * does a drop on one or both affected tips.  The sites to       *
 * be flipped are chosen randomly amongst the eligible sites of  *
 * that creature.                                                */

/* WARNING:  the logic for dropping multiple lines here works    *
 * only for two tips, not in the general case!                   */

boolean flipdrop(option_struct *op, data_fmt *data, tree *tr,
   long numdrop)
{
long i, j, flipsite, numcreatures, numslice, numcategs, datanumsites,
   *oldsiteptr, *flippedsites, *siteptr;
double ***temp;
creature *cr;
tree *oldtree;
linlist *alines;
brlist *brlines;
boolean accept_change, *nodenumber;
node *p;
long numnodes, numalines;

datanumsites = getdata_nummarkers(op,data);
siteptr = data->siteptr;

oldtree = copytree(op,data,tr);
oldsiteptr = (long *)calloc(datanumsites,sizeof(long));
memcpy(oldsiteptr,siteptr,datanumsites*sizeof(long));
numnodes = 2 * oldtree->numcoals + 2;
nodenumber = (boolean *)calloc(1,numnodes*sizeof(boolean));
for(i = 0; i < numnodes; i++) nodenumber[i] = TRUE;

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

temp = (double ***)calloc(op->categs,sizeof(double **));
temp[0] = (double **)calloc(op->categs*numslice,sizeof(double *));
for(i = 1; i < op->categs; i++) temp[i] = temp[0] + i*numslice;
temp[0][0] = (double *)calloc(op->categs*numslice*4,sizeof(double));
for(i = 0; i < op->categs; i++)
   for(j = 0; j < numslice; j++)
      temp[i][j] = temp[0][0] + i*numslice*4 + j*4;


flippedsites = (long *)calloc(NUMHAPFLIP,sizeof(long));

numcreatures = getdata_numtips(op,data)/NUMHAPLOTYPES;
do {
   cr = &(tr->creatures[(long)(randum()*numcreatures)]);
} while(cr->numflipsites == 0);

/* debug DEBUG warning WARNING--this code is two haplotype specific */
/* and it is also dna/snp specific!!!!                              */
for(i = 0; i < NUMHAPFLIP; i++) {
   flippedsites[i] = flipsite = choose_flipsite(op,data,cr,i,flippedsites);
   /* first flip the x-arrays */
   memcpy(temp[0][0],cr->haplotypes[0]->x->s[flipsite][0][0],
      numcategs*numslice*4*sizeof(double));
   memcpy(cr->haplotypes[0]->x->s[flipsite][0][0],
      cr->haplotypes[1]->x->s[flipsite][0][0],
      numcategs*numslice*4*sizeof(double));
   memcpy(cr->haplotypes[1]->x->s[flipsite][0][0],temp[0][0],
      numcategs*numslice*4*sizeof(double));
   /* then flip the actual data */
   datachangesite(op,data,cr->haplotypes[0],cr->haplotypes[1],flipsite);
}

/* now cut off the two tips */
alines = NULL;
brlines = NULL;
numalines = 0;
if (numdrop == 1) { /* drop one lineage at random */
  if (randum()>= 0.5) i=0;
  else i=1;
  p = cr->haplotypes[i];
  if (op->fc) addlinlist(&alines,p,count_activefc(op,data,p));
  else addlinlist(&alines,p,count_active(op,data,p));
  numalines++;
  tag_futile(p->back);
  traverse_flagbelow(curtree,p->back);
  remove_futile(op,data,curtree,oldtree,p,nodenumber,&brlines);
} else { /* drop both lineages */
  for (i = 0; i < NUMHAPLOTYPES ; i++) {
    if (randum()>= 0.5) i=0;
    else i=1;
    p = cr->haplotypes[i];
    if (op->fc) addlinlist(&alines,p,count_activefc(op,data,p));
    else addlinlist(&alines,p,count_active(op,data,p));
    numalines++;
    tag_futile(p->back);
    traverse_flagbelow(curtree,p->back);
    remove_futile(op,data,curtree,oldtree,p,nodenumber,&brlines);
  }
  flipdrop_prunebr(op,data,curtree,&brlines);
}

accept_change = growlineages(op,data,curtree,oldtree,
  curtree->tymelist,&brlines,&alines,&numalines,0L,&nodenumber,
  &numnodes,siteptr);
if (accept_change) {
   renumber_nodes(op,data,curtree);
   remove_excess_tree(op,data,curtree,siteptr);
   renumber_nodes(op,data,curtree);
   rebuild_alias(op,data,siteptr);
   for(i = 0; i < NUMHAPLOTYPES; i++)
      traverse_flagbelow(curtree,cr->haplotypes[i]);
#if !ALWAYS_ACCEPT
   localeval(op,data,tr->root->back,FALSE);
   tr->coalprob = coalprob(op,data,tr,theta0,rec0);
#endif
#if ALWAYS_ACCEPT
   tr->likelihood = oldtree->likelihood;
#endif
   accept_change = testratio(op,data,oldtree,tr,'f');
}

free(nodenumber);
freelinlist(alines);

if (!accept_change) {
   for(i = 0; i < NUMHAPFLIP; i++)
      datachangesite(op,data,cr->haplotypes[1],cr->haplotypes[0],
         flippedsites[i]);
   freetree(op,data,curtree);
   curtree = oldtree;
   memcpy(siteptr,oldsiteptr,datanumsites*sizeof(long));
} else freetree(op,data,oldtree);

free(oldsiteptr);
free(flippedsites);
free(temp[0][0]);
free(temp[0]);
free(temp);

return(accept_change);

} /* flipdrop */


/*******************************************************************
 * addfractrecomb() adds the constant FRACTRECOMB to the number of *
 * recombinations scored for the last tree.                        *
 * This function should only be called at the end of a chain, and  *
 * then only if the trees scored during that chain have no         *
 * recombinations in any of them!                                  */
void addfractrecomb(option_struct *op, data_fmt *data, long chain,
   treerec ***treesum)
{
long refchain, chaintype;
treerec *lasttree;

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

lasttree = &treesum[locus][refchain][op->numout[chaintype]-1];

lasttree->numrecombs += FRACTRECOMB;

} /* addfractrecomb */


/***********************************************************************
* add_one_recomb() runs a special rearrangement that always begins    *
* with a recombination event.  The resulting tree is always accepted, *
* and then is forcibly sampled (replacing the last tree of the chain).*/
void add_one_recomb(option_struct *op, data_fmt *data, tree *tr,
   long *siteptr, long chain)
{
long i, numalines, numnodes, *oldsiteptr, chaintype;
double offset, newtyme, pc, pr;
linlist *alines;
brlist *brlines;
tlist *tymeslice;
node *cutpt;
boolean *nodenumber, succeeded;
tree *oldtree;

oldtree = copytree(op,data,curtree);
numnodes = 2 * curtree->numcoals + 2;
nodenumber = (boolean *)calloc(numnodes,sizeof(boolean));
for(i = 0; i < numnodes; i++) nodenumber[i] = TRUE;
if(op->datatype == 'n' || op->datatype == 's') {
   oldsiteptr = (long *)calloc(getdata_nummarkers(op,data),sizeof(long));
   memcpy(oldsiteptr,siteptr,getdata_nummarkers(op,data)*sizeof(long));
}
alines = NULL;
brlines = NULL;

cutpt = tr->nodep[1];
tymeslice = tr->tymelist;
offset = 0.0;

if (op->fc) addlinlist(&alines,cutpt,count_activefc(op,data,cutpt));
else addlinlist(&alines,cutpt,count_active(op,data,cutpt));
numalines = 1;

tag_futile(cutpt->back);
traverse_flagbelow(curtree,cutpt->back);
remove_futile(op,data,curtree,oldtree,cutpt,nodenumber,&brlines);

/* startoff with 1 recombination */
/* generate the time the recombination will occur at */
pc = 2.0 * countinactive(op,tymeslice,alines)/theta0;
pr = rec0 * alines->activesites;
newtyme = tymeslice->eventnode->tyme + offset -
   log(randum())/(pc + pr);

/* find the tymeslice that corresponds with that time */
for(; tymeslice->age < newtyme; tymeslice = tymeslice->succ)
   ;

/* put in the recombination */
succeeded = recomb(op,data,tymeslice,&alines,&numalines,newtyme,
   siteptr,&nodenumber,&numnodes);
tymeslice = tymeslice->succ;

succeeded = growlineages(op,data,curtree,oldtree,tymeslice,&brlines,
   &alines,&numalines,offset,&nodenumber,&numnodes,siteptr);


free(nodenumber);
freelinlist(alines);

if (succeeded) {
   renumber_nodes(op,data,curtree);
   remove_excess_tree(op,data,curtree,siteptr);
   renumber_nodes(op,data,curtree);
   traverse_flagbelow(curtree,cutpt->back);
   if (op->datatype == 'n' || op->datatype == 's')
      rebuild_alias(op,data,siteptr);
   localeval(op,data,curtree->root->back,FALSE);
   curtree->coalprob = coalprob(op,data,curtree,theta0,rec0);

   /* now actually sample the generated tree */
   chaintype = TYPE_CHAIN(chain);
   op->numout[chaintype]--;

   scoretree(op,data,chain);
   op->numout[chaintype]++;
}

freetree(op,data,curtree);
curtree = oldtree;
if (op->datatype == 'n' || op->datatype == 's') {
   memcpy(siteptr,oldsiteptr,getdata_nummarkers(op,data)*sizeof(long));
   free(oldsiteptr);
}


} /* add_one_recomb */