best_tree.c

fastDNAml is an attempt to solve the same problem as DNAML, but to do so faster and using less memo

#include  <stdlib.h>
#include  <stdio.h>
#include  <string.h>
#include  "fastDNAml_types.h"
#include  "fastDNAml_funcs.h"


/*=======================================================================*/
/*                    Best tree handling for dnaml                       */
/*=======================================================================*/

/*  Tip value comparisons
 *
 *  Use void pointers to hide type from other routines.  Only tipValPtr and
 *  cmpTipVal need to be changed to alter the nature of the values compared
 *  (e.g., names instead of node numbers).
 *
 *    cmpTipVal(tipValPtr(nodeptr p), tipValPtr(nodeptr q)) == -1, 0 or 1.
 *
 *  This provides direct comparison of tip values (for example, for
 *  definition of tr->start).
 */


void  *tipValPtr (nodeptr p)
  { return  (void *) & p->number;
    }


int  cmpTipVal (void *v1, void *v2)
  { /* cmpTipVal */
    int  i1, i2;

    i1 = *((int *) v1);
    i2 = *((int *) v2);
    return  (i1 < i2) ? -1 : ((i1 == i2) ? 0 : 1);
  } /* cmpTipVal */


/*  These are the only routines that need to UNDERSTAND topologies */


topol  *setupTopol (int maxtips, int nsites)
  { /* setupTopol */
    topol   *tpl;

    if (! (tpl = (topol *) Malloc(sizeof(topol))) || 
        ! (tpl->links = (connptr) Malloc((2*maxtips-3) * sizeof(connect))) || 
        (nsites && ! (tpl->log_f
                = (double *) Malloc(nsites * sizeof(double))))) {
      printf("ERROR: Unable to get topology memory");
      tpl = (topol *) NULL;
      }

    else {
      if (nsites == 0)  tpl->log_f = (double *) NULL;
      tpl->likelihood  = unlikely;
      tpl->start       = (node *) NULL;
      tpl->nextlink    = 0;
      tpl->ntips       = 0;
      tpl->nextnode    = 0;
      tpl->opt_level   = 0;     /* degree of branch swapping explored */
      tpl->scrNum      = 0;     /* position in sorted list of scores */
      tpl->tplNum      = 0;     /* position in sorted list of trees */
      tpl->log_f_valid = 0;     /* log_f value sites */
      tpl->prelabeled  = TRUE;
      tpl->smoothed    = FALSE; /* branch optimization converged? */
      }

    return  tpl;
  } /* setupTopol */


void  freeTopol (topol *tpl)
  { /* freeTopol */
    Free(tpl->links);
    if (tpl->log_f)  Free(tpl->log_f);
    Free(tpl);
  } /* freeTopol */


int  saveSubtree (nodeptr p, topol *tpl)
    /*  Save a subtree in a standard order so that earlier branches
     *  from a node contain lower value tips than do second branches from
     *  the node.  This code works with arbitrary furcations in the tree.
     */
  { /* saveSubtree */
    connptr  r, r0;
    nodeptr  q, s;
    int      t, t0, t1;

    r0 = tpl->links;
    r = r0 + (tpl->nextlink)++;
    r->p = p;
    r->q = q = p->back;
    r->z = p->z;
    r->descend = 0;                     /* No children (yet) */

    if (q->tip) {
      r->valptr = tipValPtr(q);         /* Assign value */
      }

    else {                              /* Internal node, look at children */
      s = q->next;                      /* First child */
      do {
        t = saveSubtree(s, tpl);        /* Generate child's subtree */

        t0 = 0;                         /* Merge child into list */
        t1 = r->descend;
        while (t1 && (cmpTipVal(r0[t1].valptr, r0[t].valptr) < 0)) {
          t0 = t1;
          t1 = r0[t1].sibling;
          }
        if (t0) r0[t0].sibling = t;  else  r->descend = t;
        r0[t].sibling = t1;

        s = s->next;                    /* Next child */
        } while (s != q);

      r->valptr = r0[r->descend].valptr;   /* Inherit first child's value */
      }                                 /* End of internal node processing */

    return  r - r0;
  } /* saveSubtree */


nodeptr  minSubtreeTip (nodeptr  p0)
  { /* minTreeTip */
    nodeptr  minTip, p, testTip;

    if (p0->tip) return p0;

    p = p0->next;
    minTip = minSubtreeTip(p->back);
    while ((p = p->next) != p0) {
      testTip = minSubtreeTip(p->back);
      if (cmpTipVal(tipValPtr(testTip), tipValPtr(minTip)) < 0)
        minTip = testTip;
      }
    return minTip;
  } /* minTreeTip */


nodeptr  minTreeTip (nodeptr  p)
  { /* minTreeTip */
    nodeptr  minp, minpb;

    minp  = minSubtreeTip(p);
    minpb = minSubtreeTip(p->back);
    return cmpTipVal(tipValPtr(minp), tipValPtr(minpb)) < 0 ? minp : minpb;
  } /* minTreeTip */


void saveTree (tree *tr, topol *tpl)
    /*  Save a tree topology in a standard order so that first branches
     *  from a node contain lower value tips than do second branches from
     *  the node.  The root tip should have the lowest value of all.
     */
  { /* saveTree */
    connptr  r;
    double  *tr_log_f, *tpl_log_f;
    int  i;

    tpl->nextlink = 0;                             /* Reset link pointer */
    r = tpl->links + saveSubtree(minTreeTip(tr->start), tpl);  /* Save tree */
    r->sibling = 0;

    tpl->likelihood = tr->likelihood;
    tpl->start      = tr->start;
    tpl->ntips      = tr->ntips;
    tpl->nextnode   = tr->nextnode;
    tpl->opt_level  = tr->opt_level;
    tpl->prelabeled = tr->prelabeled;
    tpl->smoothed   = tr->smoothed;

    if (tpl_log_f = tpl->log_f) {
      tr_log_f  = tr->log_f;
      i = tpl->log_f_valid = tr->log_f_valid;
      while (--i >= 0)  *tpl_log_f++ = *tr_log_f++;
      }
    else {
      tpl->log_f_valid = 0;
      }
  } /* saveTree */


void copyTopol (topol *tpl1, topol *tpl2)
  { /* copyTopol */
    connptr  r1, r2, r10, r20;
    double  *tpl1_log_f, *tpl2_log_f;
    int  i;

    r10 = tpl1->links;
    r20 = tpl2->links;
    tpl2->nextlink = tpl1->nextlink; 

    r1 = r10;
    r2 = r20;
    i = 2 * tpl1->ntips - 3;
    while (--i >= 0) {
      r2->z = r1->z;
      r2->p = r1->p;
      r2->q = r1->q;
      r2->valptr = r1->valptr;
      r2->descend = r1->descend; 
      r2->sibling = r1->sibling; 
      r1++;
      r2++;
      }

    if (tpl1->log_f_valid && tpl2->log_f) {
      tpl1_log_f = tpl1->log_f;
      tpl2_log_f = tpl2->log_f;
      tpl2->log_f_valid = i = tpl1->log_f_valid;
      while (--i >= 0)  *tpl2_log_f++ = *tpl1_log_f++;
      }
    else {
      tpl2->log_f_valid = 0;
      }

    tpl2->likelihood = tpl1->likelihood;
    tpl2->start      = tpl1->start;
    tpl2->ntips      = tpl1->ntips;
    tpl2->nextnode   = tpl1->nextnode;
    tpl2->opt_level  = tpl1->opt_level;
    tpl2->prelabeled = tpl1->prelabeled;
    tpl2->scrNum     = tpl1->scrNum;
    tpl2->tplNum     = tpl1->tplNum;
    tpl2->smoothed   = tpl1->smoothed;
  } /* copyTopol */


boolean restoreTree (topol *tpl, tree *tr)
  { /* restoreTree */
    void  hookup();
    boolean  initrav();

    connptr  r;
    nodeptr  p, p0;
    double  *tr_log_f, *tpl_log_f;
    int  i;

/*  Clear existing connections */

    for (i = 1; i <= 2*(tr->mxtips) - 2; i++) {  /* Uses p = p->next at tip */
      p0 = p = tr->nodep[i];
      do {
        p->back = (nodeptr) NULL;
        p = p->next;
        } while (p != p0);
      }

/*  Copy connections from topology */

    for (r = tpl->links, i = 0; i < tpl->nextlink; r++, i++) {
      hookup(r->p, r->q, r->z);
      }

    tr->likelihood = tpl->likelihood;
    tr->start      = tpl->start;
    tr->ntips      = tpl->ntips;
    tr->nextnode   = tpl->nextnode;
    tr->opt_level  = tpl->opt_level;
    tr->prelabeled = tpl->prelabeled;
    tr->smoothed   = tpl->smoothed;

    if (tpl_log_f = tpl->log_f) {
      tr_log_f = tr->log_f;
      i = tr->log_f_valid = tpl->log_f_valid;
      while (--i >= 0)  *tr_log_f++ = *tpl_log_f++;
      }
    else {
      tr->log_f_valid = 0;
      }

    return (initrav(tr, tr->start) && initrav(tr, tr->start->back));
  } /* restoreTree */


int initBestTree (bestlist *bt, int newkeep, int numsp, int sites)
  { /* initBestTree */
    int  i, nlogf;