fastdnaml.c

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

  { /* unlink_pid */
    char scr[512];
    (void) sprintf(scr, "%s.%d", filenm, getpid());
    unlink(scr);
  } /* unlink_pid */
#endif


void  writeCheckpoint (tree *tr)
  { /* writeCheckpoint */
    char   filename[128] = "";
    FILE  *checkpointf;
    static int  i=1;
  /*make_filename(filename,workdir,run_id,cp_ext,i++);*/  /* keep all cpt files */
    make_filename(filename,workdir,run_id,cp_ext,0);      /* keep only last one */
  /*checkpointf = fopen(filename,"a");*/
    checkpointf = fopen(filename,"w");
    if (checkpointf) {
      treeOut(checkpointf, tr, treeNewick);
      (void) fclose(checkpointf);
      }
  } /* writeCheckpoint */


node * findAnyTip(nodeptr p)
  { /* findAnyTip */
    return  p->tip ? p : findAnyTip(p->next->back);
  } /* findAnyTip */


boolean  optimize (tree *tr, int maxtrav, bestlist *bt)
  { /* optimize */
    nodeptr  p;
    int    mintrav, tested;

    if (tr->ntips < 4)  return  TRUE;

/* Checkpoint the starting tree. */
    writeCheckpoint(tr);

    if (maxtrav > tr->ntips - 3)  maxtrav = tr->ntips - 3;
    if (maxtrav <= tr->opt_level)  return  TRUE;

    printf("      Doing %s rearrangements\n",
             (maxtrav == 1)            ? "local" :
             (maxtrav < tr->ntips - 3) ? "regional" : "global");

/* Loop while tree gets better.  */
    do {
      (void) startOpt(bt, tr);
      mintrav = tr->opt_level + 1;

      /* rearrange must start from a tip or it will miss some trees */
      p = findAnyTip(tr->start);
      tested = rearrange(tr, p->back, mintrav, maxtrav, bt);
      if (tested == badRear)  return FALSE;

#     if Master
        if (! getReturnedTrees(tr, bt, tested)) return FALSE;
#     endif

      bt->numtrees += tested;
      (void) setOptLevel(bt, maxtrav);
      if (! recallBestTree(bt, 1, tr)) return FALSE;   /* recover best tree */

      printf("      Tested %d alternative trees\n", tested);
      if (bt->improved) {
        printf("      Ln Likelihood =%14.5f\n", tr->likelihood);
      }

      /* Checkpoint the new tree. */
      writeCheckpoint(tr);
    } while (maxtrav > tr->opt_level);

    return TRUE;
  } /* optimize */


nodeptr uprootTree (tree *tr, nodeptr p)
  { /* uprootTree */
    nodeptr  q, r, s, start;
    int      n;

    if (p->tip || p->back) {
      printf("ERROR: Unable to uproot tree.\n");
      printf("       Inappropriate node marked for removal.\n");
      return (nodeptr) NULL;
      }

    n = --(tr->nextnode);               /* last internal node added */
    if (n != tr->mxtips + tr->ntips - 1) {
      printf("ERROR: Unable to uproot tree.  Inconsistent\n");
      printf("       number of tips and nodes for rooted tree.\n");
      return (nodeptr) NULL;
      }

    q = p->next->back;                  /* remove p from tree */
    r = p->next->next->back;
    hookup(q, r, tr->userlen ? (q->z * r->z) : defaultz);

    start = (r->tip || (! q->tip)) ? r : r->next->next->back;

    if (tr->ntips > 2 && p->number != n) {
      q = tr->nodep[n];            /* transfer last node's conections to p */
      r = q->next;
      s = q->next->next;
      hookup(p,             q->back, q->z);   /* move connections to p */
      hookup(p->next,       r->back, r->z);
      hookup(p->next->next, s->back, s->z);
      if (start->number == q->number) start = start->back->back;
      q->back = r->back = s->back = (nodeptr) NULL;
      }
    else {
      p->back = p->next->back = p->next->next->back = (nodeptr) NULL;
      }

    tr->rooted = FALSE;
    return  start;
  } /* uprootTree */

/*=======================================================================*/
/*                 end of section E                                      */      
/*=======================================================================*/



/*=======================================================================*/
/*                       section G                                       */      
/*=======================================================================*/
/* Evaluate a user tree.
 */
boolean treeEvaluate (tree *tr, bestlist *bt)
  { /* treeEvaluate */

    if (Slave || ! tr->userlen) {
      if (! smoothTree(tr, 4 * smoothings)) return FALSE;
    /*fputc('S',dbgfp);fputc('\n',dbgfp);*/ /*DKB-debug*/
    }

    if (evaluate(tr, tr->start) == badEval)  return FALSE;

#   if ! Slave
      (void) saveBestTree(bt, tr);
#   endif
    return TRUE;
  } /* treeEvaluate */


boolean  showBestTrees (bestlist *bt, tree *tr, analdef *adef, FILE *treefile)
  { /* showBestTrees */
    int     rank;

    for (rank = 1; rank <= bt->nvalid; rank++) {
      if (rank > 1) {
        if (rank != recallBestTree(bt, rank, tr))  break;
        }
      if (evaluate(tr, tr->start) == badEval) return FALSE;
      if (tr->outgrnode->back)  tr->start = tr->outgrnode;
      printTree(tr, adef);
      summarize(tr);
      if (treefile)  treeOut(treefile, tr, adef->trout);
      }

    return TRUE;
  } /* showBestTrees */


boolean cmpBestTrees (bestlist *bt, tree *tr)
  { /* cmpBestTrees */
    double  sum, sum2, sd, temp, wtemp, bestscore;
    double *log_f0, *log_f0_ptr;      /* Save a copy of best log_f */
    double *log_f_ptr;
    int     i, j, num, besttips;

    num = bt->nvalid;
    if ((num <= 1) || (tr->cdta->wgtsum <= 1)) return TRUE;

    if (! (log_f0 = (double *) Malloc(sizeof(double) * tr->cdta->endsite))) {
      printf("ERROR: cmpBestTrees unable to obtain space for log_f0\n");
      return FALSE;
      }

    printf("Tree      Ln L        Diff Ln L       Its S.D.");
    printf("   Significantly worse?\n\n");

    for (i = 1; i <= num; i++) {
      if (i != recallBestTree(bt, i, tr))  break;
      if (! (tr->log_f_valid))  {
        if (evaluate(tr, tr->start) == badEval) return FALSE;
        }

      printf("%3d%14.5f", i, tr->likelihood);
      if (i == 1) {
        printf("  <------ best\n");
        besttips = tr->ntips;
        bestscore = tr->likelihood;
        log_f0_ptr = log_f0;
        log_f_ptr  = tr->log_f;
        for (j = 0; j < tr->cdta->endsite; j++)  *log_f0_ptr++ = *log_f_ptr++;
        }
      else if (tr->ntips != besttips)
        printf("  (different number of species)\n");
      else {
        sum = sum2 = 0.0;
        log_f0_ptr = log_f0;
        log_f_ptr  = tr->log_f;
        for (j = 0; j < tr->cdta->endsite; j++) {
          temp  = *log_f0_ptr++ - *log_f_ptr++;
          wtemp = tr->cdta->aliaswgt[j] * temp;
          sum  += wtemp;
          sum2 += wtemp * temp;
          }
        sd = sqrt( tr->cdta->wgtsum * (sum2 - sum*sum / tr->cdta->wgtsum)
                                    / (tr->cdta->wgtsum - 1) );
        printf("%14.5f%14.4f", tr->likelihood - bestscore, sd);
        printf("           %s\n", (sum > 1.95996 * sd) ? "Yes" : " No");
        }
      }

    Free(log_f0);
    printf("\n\n");

    return TRUE;
  } /* cmpBestTrees */


boolean  makeUserTree (FILE *seqfp, tree *tr, bestlist *bt, analdef *adef)
  { /* makeUserTree */
    char   filename[128];
    FILE  *treefile;
    int    nusertrees, which;

    nusertrees = adef->numutrees;

    printf("User-defined %s:\n\n", (nusertrees == 1) ? "tree" : "trees");

    if(adef->trout) {
      make_filename(filename,workdir,run_id,tr_ext,0);
      treefile = fopen(filename,"w");
    }
    else {
      treefile = (FILE*)NULL;
    }

    for (which = 1; which <= nusertrees; which++) {
      if (! treeReadLen(seqfp, tr)) return FALSE;
      if (! treeEvaluate(tr, bt))    return FALSE;
      if (tr->global <= 0) {
        if (tr->outgrnode->back)  tr->start = tr->outgrnode;
        printTree(tr, adef);
        summarize(tr);
        if (treefile)  treeOut(treefile, tr, adef->trout);
        }
      else {
        printf("%6d:  Ln Likelihood =%14.5f\n", which, tr->likelihood);
        }
      }

    if (tr->global > 0) {
      putchar('\n');
      if (! recallBestTree(bt, 1, tr))  return FALSE;
      printf("      Ln Likelihood =%14.5f\n", tr->likelihood);
      if (! optimize(tr, tr->global, bt))  return FALSE;
      if (tr->outgrnode->back)  tr->start = tr->outgrnode;
      printTree(tr, adef);
      summarize(tr);
      if (treefile)  treeOut(treefile, tr, adef->trout);
      }

    if (treefile) {
      (void) fclose(treefile);
      printf("Tree also written to %s\n", filename);
      }

    putchar('\n');

    (void) cmpBestTrees(bt, tr);
    return TRUE;
  } /* makeUserTree */


boolean makeDenovoTree (FILE *seqfp, tree *tr, bestlist *bt, analdef *adef)
  { /* makeDenovoTree */
    char     c;
    char     filename[128];
    FILE     *treefile;
    nodeptr  p;
    int      *enterorder;      /*  random entry order */
    int      i, j, k, nextsp, newsp, maxtrav, tested;

    double randum();

    enterorder = (int*)Malloc(sizeof(int)*(tr->mxtips+1));
    if (! enterorder) {
       fprintf(stderr, "makeDenovoTree: Malloc failure for enterorder\n");
       return 0;
     }

/* Either restart from a user-defined tree ... */
    if (adef->restart) {
      printf("Restarting from tree with the following sequences:\n");
      tr->userlen = TRUE;
      if(!treeReadLen(seqfp,tr))             return FALSE;
      if(!smoothTree(tr,smoothings))         return FALSE;
      if(evaluate(tr,tr->start) == badEval)  return FALSE;
      if(saveBestTree(bt,tr) < 1)            return FALSE;

      for(i=1, j=tr->ntips; i<=tr->mxtips; i++) { /* find loose tips */
        if(! tr->nodep[i]->back) {
          enterorder[++j] = i;
        }
        else {
          printf("   %s\n", tr->nodep[i]->name);
          if(infol>2 && monitor_id!=INVALID_ID) {
            if(i>3) send_msg(&c,0,monitor_id,DNAML_ADD_SPECS);
          }
        }
      }
      putchar('\n');
    }

/* ... or else start from scratch */
    else {
      tr->ntips = 0;
      for (i = 1; i <= tr->mxtips; i++) enterorder[i] = i;
    }

/* Set the order for adding taxa */
    if(adef->jumble)  for(i=tr->ntips+1; i<=tr->mxtips; i++) {
      j = randum(&(adef->jumble))*(tr->mxtips - tr->ntips) + tr->ntips + 1;
      k = enterorder[j];
      enterorder[j] = enterorder[i];
      enterorder[i] = k;
    }

    bt->numtrees = 1;
    if (tr->ntips < tr->mxtips)  printf("Adding species:\n");

    if (tr->ntips == 0) {
      for(i=1; i<=3; i++) {
        printf("   %s\n", tr->nodep[enterorder[i]]->name);
      }
      tr->nextnode = tr->mxtips + 1;
      if(!buildSimpleTree(tr,enterorder[1],enterorder[2],enterorder[3]))
        return FALSE;
    }

    while (tr->ntips < tr->mxtips || tr->opt_level < tr->global) {
      maxtrav = (tr->ntips == tr->mxtips) ? tr->global : tr->partswap;
      if (maxtrav > tr->ntips - 3)  maxtrav = tr->ntips - 3;

      if (tr->opt_level >= maxtrav) {
        nextsp = ++(tr->ntips);
        newsp = enterorder[nextsp];
        p = tr->nodep[newsp];
        printf("   %s\n", p->name);

        if(infol>1 && monitor_id!=INVALID_ID) {
          if(nextsp % DNAML_STEP_TIME_COUNT == 1) {
            record_times(&myproc.stats);
            send_msg(&myproc.stats,1,monitor_id,DNAML_STEP_TIME);
          }
          if(infol>2 && monitor_id!=INVALID_ID) {
            send_msg(&nextsp,1,monitor_id,DNAML_ADD_SPECS);
          }
        }

        (void) buildNewTip(tr, p);

        resetBestTree(bt);
        cacheZ(tr);
        tested = addTraverse(tr, p->back, findAnyTip(tr->start)->back,
                             1, tr->ntips - 2, bt, adef->qadd);
        if (tested == badRear) return FALSE;
        bt->numtrees += tested;

#       if Master
          getReturnedTrees(tr, bt, tested);
#       endif
        printf("      Tested %d alternative trees\n", tested);
        (void) recallBestTree(bt, 1, tr);
        if (! tr->smoothed) {
          if (! smoothTree(tr, smoothings))        return FALSE;
          if (evaluate(tr, tr->start) == badEval)  return FALSE;
          (void) saveBestTree(bt, tr);
        }

        if (tr->ntips == 4)  tr->opt_level = 1;  /* All 4 taxon trees done */
        maxtrav = (tr->ntips == tr->mxtips) ? tr->global : tr->partswap;
        if (maxtrav > tr->ntips - 3)  maxtrav = tr->ntips - 3;
      }

      printf("      Ln Likelihood =%14.5f\n", tr->likelihood);
      if (! optimize(tr, maxtrav, bt)) return FALSE;
    }

    printf("\nExamined %d %s\n", bt->numtrees,
                                 bt->numtrees != 1 ? "trees" : "tree");

    if(adef->trout) {
      make_filename(filename,workdir,run_id,tr_ext,0);
      treefile = fopen(filename,"w");
    }
    else {
      treefile = (FILE*)NULL;
    }

    (void) showBestTrees(bt, tr, adef, treefile);
    if (treefile) {
      (void) fclose(treefile);
      printf("Tree also written to %s\n\n", filename);
    }

    (void) cmpBestTrees(bt, tr);

#   if DeleteCheckpointFile
      unlink_pid(checkpointname);
#   endif

    Free(enterorder);

    return TRUE;
  } /* makeDenovoTree */