fastdnaml.c

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

    q = p->back;
    z0 = q->z;
    if ((z = makenewz(tr,p,q,z0,newzpercycle)) == badZ) return FALSE;
  /*fprintf(dbgfp," r%er ",z);*/ /*DKB-debug*/
    p->z = q->z = z;
    if(ABS(z-z0)>deltaz) tr->smoothed = FALSE;
    return TRUE;
  } /* update */


boolean smooth (tree *tr, nodeptr p)
  { /* smooth */
    nodeptr  q;
    if(!update(tr,p))               return FALSE; /*  Adjust branch */
    if(!p->tip) {                                 /*  Adjust descendants */
        q = p->next;
        while (q != p) {
          if(!smooth(tr,q->back))   return FALSE;
          q = q->next;
        }
#     if ReturnSmoothedView
        if(!newview(tr,p)) return FALSE;
#     endif
    }
    return TRUE;
  } /* smooth */


boolean smoothTree (tree *tr, int maxtimes)
  { /* smoothTree */
    nodeptr  p, q;

    p = tr->start;
  /*fputc('\n',dbgfp);*/ /*DKB-debug*/
    while (--maxtimes >= 0) {
      tr->smoothed = TRUE;
      if(!smooth(tr,p->back))          return FALSE;
      if(!p->tip) {
        q = p->next;
        while (q != p) {
          if(!smooth(tr,q->back))      return FALSE;
          q = q->next;
        }
      }
    /*fputc('s',dbgfp);*/  /*DKB-debug*/
      if (tr->smoothed)  break;
    }
    return TRUE;
  } /* smoothTree */


boolean localSmooth (tree *tr, nodeptr p, int maxtimes)
  { /* localSmooth -- Smooth branches around p */
    nodeptr  q;

    if (p->tip) return FALSE;            /* Should be an error */

    while (--maxtimes >= 0) {
      tr->smoothed = TRUE;
      q = p;
      do {
        if (! update(tr, q)) return FALSE;
        q = q->next;
      } while (q != p);
      if (tr->smoothed)  break;
    }

    tr->smoothed = FALSE;             /* Only smooth locally */
    return TRUE;
  } /* localSmooth */


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


/* Insert node p into branch q <-> q->back
 * glob -- Smooth tree globally?
 *
 *                q
 *               /.
 *           add/ .
 *             /  .
 *           pn   .
 *  s ---- p      .remove
 *           pnn  .
 *             \  .
 *           add\ .
 *               \.      pn  = p->next;
 *                r      pnn = p->next->next;
 */
boolean insert (tree *tr, nodeptr p, nodeptr q, boolean glob)
  { /* insert */
    nodeptr  r, s;

    r = q->back;
    s = p->back;

#if BestInsertAverage && ! Master
    { double  zqr, zqs, zrs, lzqr, lzqs, lzrs, lzsum, lzq, lzr, lzs, lzmax;

      if ((zqr = makenewz(tr, q, r, q->z,     iterations)) == badZ) return FALSE;
      if ((zqs = makenewz(tr, q, s, defaultz, iterations)) == badZ) return FALSE;
      if ((zrs = makenewz(tr, r, s, defaultz, iterations)) == badZ) return FALSE;

      lzqr = (zqr > zmin) ? log(zqr) : log(zmin);  /* long branches */
      lzqs = (zqs > zmin) ? log(zqs) : log(zmin);
      lzrs = (zrs > zmin) ? log(zrs) : log(zmin);
      lzsum = 0.5 * (lzqr + lzqs + lzrs);

      lzq = lzsum - lzrs;
      lzr = lzsum - lzqs;
      lzs = lzsum - lzqr;
      lzmax = log(zmax);

      if      (lzq > lzmax) {lzq = lzmax; lzr = lzqr; lzs = lzqs;} /* short */
      else if (lzr > lzmax) {lzr = lzmax; lzq = lzqr; lzs = lzrs;}
      else if (lzs > lzmax) {lzs = lzmax; lzq = lzqs; lzr = lzrs;}

      hookup(p->next,       q, exp(lzq));
      hookup(p->next->next, r, exp(lzr));
      hookup(p,             s, exp(lzs));
    }
#else
    { double  z;
      z = sqrt(q->z);
    /*fprintf(dbgfp," f%lff ",z);*/ /*DKB-debug*/
      hookup(p->next,       q, z);
      hookup(p->next->next, r, z);
    }
#endif

    if(!newview(tr,p)) return FALSE;   /*  So that p is valid at update */
    tr->opt_level = 0;

#if ! Master         /*  Smoothings are done by slave */
      if (glob) {                                    /* Smooth whole tree */
        if (! smoothTree(tr, smoothings)) return FALSE;
      }
      else {                                         /* Smooth locale of p */
        if (! localSmooth(tr, p, smoothings)) return FALSE;
      }
#else
      tr->likelihood = unlikely;
#endif
    return  TRUE;
  } /* insert */


nodeptr  removeNode (tree *tr, nodeptr p)

/*                q
                 .|
          remove. |
               .  |
             pn   |
    s ---- p      |add
             pnn  |
               .  |
          remove. |
                 .|      pn  = p->next;
                  r      pnn = p->next->next;
 */

    /* remove p and return where it was */
  { /* removeNode */
    double   zqr;
    nodeptr  q, r;

    q = p->next->back;
    r = p->next->next->back;
    zqr = q->z * r->z;
#   if ! Master
      if ((zqr = makenewz(tr, q, r, zqr, iterations)) == badZ) return (node *) NULL;
#   endif
    hookup(q, r, zqr);

    p->next->next->back = p->next->back = (node *) NULL;
    return  q;
  } /* removeNode */


boolean initrav (tree *tr, nodeptr p)
  { /* initrav */
    nodeptr  q;
    if (! p->tip) {
      q = p->next;
      do {
        if(!initrav(tr,q->back))  return FALSE;
        q = q->next;
      } while (q != p);
      if(!newview(tr,p))          return FALSE;
    }
    return TRUE;
  } /* initrav */


nodeptr buildNewTip (tree *tr, nodeptr p)
  { /* buildNewTip */
    nodeptr  q;
    q = tr->nodep[(tr->nextnode)++];
    hookup(p, q, defaultz);
    return  q;
  } /* buildNewTip */


boolean buildSimpleTree (tree *tr, int ip, int iq, int ir)
  { /* buildSimpleTree */
    /* p, q and r are tips meeting at s */
    nodeptr  p, s;
    int  i;
    i = MIN(ip, iq);
    if (ir < i)  i = ir; 
    tr->start = tr->nodep[i];
    tr->ntips = 3;
    p = tr->nodep[ip];
  /*fprintf(dbgfp,"@@@%lf@@@",defaultz);*/ /*DKB-debug*/
    hookup(p, tr->nodep[iq], defaultz);
    s = buildNewTip(tr, tr->nodep[ir]);
    return insert(tr, s, p, FALSE);  /* Smoothing is local to s */
  } /* buildSimpleTree */


void cacheZ (tree *tr)
  { /* cacheZ */
    nodeptr  p;
    int  nodes;
    nodes = tr->mxtips  +  3 * (tr->mxtips - 2);
    p = tr->nodep[1];
    while (nodes-- > 0) {p->z0 = p->z; p++;}
  } /* cacheZ */


void restoreZ (tree *tr)
  { /* restoreZ */
    nodeptr  p;
    int  nodes;
    nodes = tr->mxtips  +  3 * (tr->mxtips - 2);
    p = tr->nodep[1];
    while (nodes-- > 0) {p->z = p->z0; p++;}
  } /* restoreZ */


boolean testInsert (tree *tr, nodeptr p, nodeptr q, bestlist *bt, boolean  fast)
  { /* testInsert */
    double  qz;
    nodeptr  r;

#if ! Master
    send_msg(NULL,0,0,DNAML_SEND_TREE);
#endif
    r = q->back;             /* Save original connection */
    qz = q->z;
    if(!insert(tr,p,q,!fast)) return FALSE;

#if ! Master
    if(evaluate(tr, fast ? p->next->next : tr->start) == badEval) return FALSE;
    send_msg(NULL,0,0,DNAML_RECV_TREE);
    (void) saveBestTree(bt, tr);
#else  /* Master */
    /*fprintf(dbgfp," c%lfc \n",q->z);*/ /*DKB-debug*/
    tr->likelihood = unlikely;
    if (! sendTree(tr))  return FALSE;
#endif

/* Remove p from this branch */
    hookup(q, r, qz);
    p->next->next->back = p->next->back = (nodeptr) NULL;
    if (! fast) {            /* With fast add, other values are still OK */
      restoreZ(tr);          /*   Restore branch lengths */
#     if ! Master            /*   Regenerate x values */
        if(!initrav(tr,p->back))    return FALSE;
        if(!initrav(tr,q))          return FALSE;
        if(!initrav(tr,r))          return FALSE;
#     endif
    }
    return TRUE;
  } /* testInsert */


int addTraverse (tree *tr, nodeptr p, nodeptr q,
                 int mintrav, int maxtrav, bestlist *bt, boolean fast)
  { /* addTraverse */
    int  tested, newtested;

    tested = 0;
    if (--mintrav <= 0) {           /* Moved minimum distance? */
      if(!testInsert(tr,p,q,bt,fast))  return badRear;
      tested++;
    }

    if ((! q->tip) && (--maxtrav > 0)) {    /* Continue traverse? */
      newtested = addTraverse(tr,p,q->next->back,mintrav,maxtrav,bt,fast);
      if (newtested == badRear) return badRear;
      tested += newtested;
      newtested = addTraverse(tr,p,q->next->next->back,mintrav,maxtrav,bt,fast);
      if (newtested == badRear) return badRear;
      tested += newtested;
    }

    return tested;
  } /* addTraverse */


int  rearrange (tree *tr, nodeptr p, int mintrav, int maxtrav, bestlist *bt)
    /* rearranges the tree, globally or locally */
  { /* rearrange */
    double   p1z, p2z, q1z, q2z;
    nodeptr  p1, p2, q, q1, q2;
    int      tested, mintrav2, newtested;

    tested = 0;
    if (maxtrav < 1 || mintrav > maxtrav)  return tested;

/* Moving subtree forward in tree. */

    if (! p->tip) {
      p1 = p->next->back;
      p2 = p->next->next->back;
      if (! p1->tip || ! p2->tip) {
        p1z = p1->z;
        p2z = p2->z;
        if (! removeNode(tr, p)) return badRear;
        cacheZ(tr);
        if (! p1->tip) {
          newtested = addTraverse(tr, p, p1->next->back,
                                  mintrav, maxtrav, bt, FALSE);
          if (newtested == badRear) return badRear;
          tested += newtested;
          newtested = addTraverse(tr, p, p1->next->next->back,
                                  mintrav, maxtrav, bt, FALSE);
          if (newtested == badRear) return badRear;
          tested += newtested;
          }

        if (! p2->tip) {
          newtested = addTraverse(tr, p, p2->next->back,
                                  mintrav, maxtrav, bt, FALSE);
          if (newtested == badRear) return badRear;
          tested += newtested;
          newtested = addTraverse(tr, p, p2->next->next->back,
                                  mintrav, maxtrav, bt, FALSE);
          if (newtested == badRear) return badRear;
          tested += newtested;
          }

        hookup(p->next,       p1, p1z);  /*  Restore original tree */
        hookup(p->next->next, p2, p2z);
        if (! (initrav(tr, tr->start)
            && initrav(tr, tr->start->back))) return badRear;
        }
      }   /* if (! p->tip) */

/* Moving subtree backward in tree.  Minimum move is 2 to avoid duplicates */

    q = p->back;
    if (! q->tip && maxtrav > 1) {
      q1 = q->next->back;
      q2 = q->next->next->back;
      if (! q1->tip && (!q1->next->back->tip || !q1->next->next->back->tip) ||
          ! q2->tip && (!q2->next->back->tip || !q2->next->next->back->tip)) {
        q1z = q1->z;
        q2z = q2->z;
        if (! removeNode(tr, q)) return badRear;
        cacheZ(tr);
        mintrav2 = mintrav > 2 ? mintrav : 2;

        if (! q1->tip) {
          newtested = addTraverse(tr, q, q1->next->back,
                                  mintrav2 , maxtrav, bt, FALSE);
          if (newtested == badRear) return badRear;
          tested += newtested;
          newtested = addTraverse(tr, q, q1->next->next->back,
                                  mintrav2 , maxtrav, bt, FALSE);
          if (newtested == badRear) return badRear;
          tested += newtested;
          }

        if (! q2->tip) {
          newtested = addTraverse(tr, q, q2->next->back,
                                  mintrav2 , maxtrav, bt, FALSE);
          if (newtested == badRear) return badRear;
          tested += newtested;
          newtested = addTraverse(tr, q, q2->next->next->back,
                                  mintrav2 , maxtrav, bt, FALSE);
          if (newtested == badRear) return badRear;
          tested += newtested;
          }

        hookup(q->next,       q1, q1z);  /*  Restore original tree */
        hookup(q->next->next, q2, q2z);
        if (! (initrav(tr, tr->start)
            && initrav(tr, tr->start->back))) return badRear;
        }
      }   /* if (! q->tip && maxtrav > 1) */

/* Move other subtrees */

    if (! p->tip) {
      newtested = rearrange(tr, p->next->back,       mintrav, maxtrav, bt);
      if (newtested == badRear) return badRear;
      tested += newtested;
      newtested = rearrange(tr, p->next->next->back, mintrav, maxtrav, bt);
      if (newtested == badRear) return badRear;
      tested += newtested;
      }

    return  tested;
  } /* rearrange */
/*=======================================================================*/
/*                 end of section D                                      */
/*=======================================================================*/



/*=======================================================================*/
/*                 section E                                             */
/*=======================================================================*/

#if DeleteCheckpointFile
void  unlink_pid (char *filenm)