fastdnaml.c

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

      for (i = 0; i < tr->cdta->endsite; i++) {
        code = *yptr++;
        l->a =  code       & 1;
        l->c = (code >> 1) & 1;
        l->g = (code >> 2) & 1;
        l->t = (code >> 3) & 1;
        l->exp = 0;
        l++;
      }
      return TRUE;
    }

/* Otherwise, an internal node needs update. */
    q = p->next->back;
    r = p->next->next->back;

    while((!p->x) || (!q->x) || (!r->x)) {
      if(!q->x) if(!newview(tr,q))  return FALSE;
      if(!r->x) if(!newview(tr,r))  return FALSE;
      if(!p->x) if(!getxnode(p))    return FALSE;
    }

    lp = p->x->lv;

    lq = q->x->lv;
    zq = q->z;
    lzq = (zq > zmin) ? log(zq) : log(zmin);
    xvlzq = tr->rdta->xv * lzq;

    lr = r->x->lv;
    zr = r->z;
    lzr = (zr > zmin) ? log(zr) : log(zmin);
    xvlzr = tr->rdta->xv * lzr;

    { double  zzqtable[maxcategories+1], zvqtable[maxcategories+1],
              zzrtable[maxcategories+1], zvrtable[maxcategories+1],
             *zzqptr, *zvqptr, *zzrptr, *zvrptr, *rptr;
      double  fxqr, fxqy, fxqn, sumaq, sumgq, sumcq, sumtq,
              fxrr, fxry, fxrn, ki, tempi, tempj;
      int  *cptr;

#   ifdef Vectorize
      double zzq[maxpatterns], zvq[maxpatterns],
             zzr[maxpatterns], zvr[maxpatterns];
      int  cat;
#   else
      double zzq, zvq, zzr, zvr;
      int cat;
#   endif

      rptr   = &(tr->rdta->catrat[1]);
      zzqptr = &(zzqtable[1]);
      zvqptr = &(zvqtable[1]);
      zzrptr = &(zzrtable[1]);
      zvrptr = &(zvrtable[1]);

#   ifdef Vectorize
#     pragma IVDEP
#   endif

      for (i = 1; i <= tr->rdta->categs; i++) {   /* exps for each category */
        ki = *rptr++;
        *zzqptr++ = exp(ki *   lzq);
        *zvqptr++ = exp(ki * xvlzq);
        *zzrptr++ = exp(ki *   lzr);
        *zvrptr++ = exp(ki * xvlzr);
      }

      cptr = &(tr->cdta->patcat[0]);

#   ifdef Vectorize
#     pragma IVDEP
      for (i = 0; i < tr->cdta->endsite; i++) {
        cat = *cptr++;
        zzq[i] = zzqtable[cat];
        zvq[i] = zvqtable[cat];
        zzr[i] = zzrtable[cat];
        zvr[i] = zvrtable[cat];
      }

#     pragma IVDEP
      for (i = 0; i < tr->cdta->endsite; i++) {
        fxqr = tr->rdta->freqa * lq->a + tr->rdta->freqg * lq->g;
        fxqy = tr->rdta->freqc * lq->c + tr->rdta->freqt * lq->t;
        fxqn = fxqr + fxqy;
        tempi = fxqr * tr->rdta->invfreqr;
        tempj = zvq[i] * (tempi-fxqn) + fxqn;
        sumaq = zzq[i] * (lq->a - tempi) + tempj;
        sumgq = zzq[i] * (lq->g - tempi) + tempj;
        tempi = fxqy * tr->rdta->invfreqy;
        tempj = zvq[i] * (tempi-fxqn) + fxqn;
        sumcq = zzq[i] * (lq->c - tempi) + tempj;
        sumtq = zzq[i] * (lq->t - tempi) + tempj;

        fxrr = tr->rdta->freqa * lr->a + tr->rdta->freqg * lr->g;
        fxry = tr->rdta->freqc * lr->c + tr->rdta->freqt * lr->t;
        fxrn = fxrr + fxry;
        tempi = fxrr * tr->rdta->invfreqr;
        tempj = zvr[i] * (tempi-fxrn) + fxrn;
        lp->a = sumaq * (zzr[i] * (lr->a - tempi) + tempj);
        lp->g = sumgq * (zzr[i] * (lr->g - tempi) + tempj);
        tempi = fxry * tr->rdta->invfreqy;
        tempj = zvr[i] * (tempi-fxrn) + fxrn;
        lp->c = sumcq * (zzr[i] * (lr->c - tempi) + tempj);
        lp->t = sumtq * (zzr[i] * (lr->t - tempi) + tempj);
        lp->exp = lq->exp + lr->exp;

        if (lp->a < minlikelihood && lp->g < minlikelihood
         && lp->c < minlikelihood && lp->t < minlikelihood) {
          lp->a   *= twotothe256;
          lp->g   *= twotothe256;
          lp->c   *= twotothe256;
          lp->t   *= twotothe256;
          lp->exp += 1;
        }
        lp++;
        lq++;
        lr++;
      }

#   else  /*  Not Vectorize  */
      for (i = 0; i < tr->cdta->endsite; i++) {
        cat = *cptr++;
        zzq = zzqtable[cat];
        zvq = zvqtable[cat];
        fxqr = tr->rdta->freqa * lq->a + tr->rdta->freqg * lq->g;
        fxqy = tr->rdta->freqc * lq->c + tr->rdta->freqt * lq->t;
        fxqn = fxqr + fxqy;
        tempi = fxqr * tr->rdta->invfreqr;
        tempj = zvq * (tempi-fxqn) + fxqn;
        sumaq = zzq * (lq->a - tempi) + tempj;
        sumgq = zzq * (lq->g - tempi) + tempj;
        tempi = fxqy * tr->rdta->invfreqy;
        tempj = zvq * (tempi-fxqn) + fxqn;
        sumcq = zzq * (lq->c - tempi) + tempj;
        sumtq = zzq * (lq->t - tempi) + tempj;

        zzr = zzrtable[cat];
        zvr = zvrtable[cat];
        fxrr = tr->rdta->freqa * lr->a + tr->rdta->freqg * lr->g;
        fxry = tr->rdta->freqc * lr->c + tr->rdta->freqt * lr->t;
        fxrn = fxrr + fxry;
        tempi = fxrr * tr->rdta->invfreqr;
        tempj = zvr * (tempi-fxrn) + fxrn;
        lp->a = sumaq * (zzr * (lr->a - tempi) + tempj);
        lp->g = sumgq * (zzr * (lr->g - tempi) + tempj);
        tempi = fxry * tr->rdta->invfreqy;
        tempj = zvr * (tempi-fxrn) + fxrn;
        lp->c = sumcq * (zzr * (lr->c - tempi) + tempj);
        lp->t = sumtq * (zzr * (lr->t - tempi) + tempj);
        lp->exp = lq->exp + lr->exp;

        if (lp->a < minlikelihood && lp->g < minlikelihood
         && lp->c < minlikelihood && lp->t < minlikelihood) {
          lp->a   *= twotothe256;
          lp->g   *= twotothe256;
          lp->c   *= twotothe256;
          lp->t   *= twotothe256;
          lp->exp += 1;
        }
        lp++;
        lq++;
        lr++;
      }
#   endif  /*  Vectorize or not  */

      return TRUE;
    }
  } /* newview */


double evaluate (tree *tr, nodeptr p)
  { /* evaluate */
    double   sum, z, lz, xvlz,
             ki, fxpa, fxpc, fxpg, fxpt, fxpr, fxpy, fxqr, fxqy,
             suma, sumb, sumc, term;

# ifdef Vectorize
    double   zz[maxpatterns], zv[maxpatterns];
# else
    double   zz, zv;
# endif

    double        zztable[maxcategories+1], zvtable[maxcategories+1],
                 *zzptr, *zvptr;
    double       *log_f, *rptr;
    likelivector *lp, *lq;
    nodeptr       q;
    int           cat, *cptr, i, *wptr;

    q = p->back;
    while ((! p->x) || (! q->x)) {
      if (! (p->x)) if (! newview(tr, p)) return badEval;
      if (! (q->x)) if (! newview(tr, q)) return badEval;
    }

    lp = p->x->lv;
    lq = q->x->lv;

    z = p->z;
    if (z < zmin) z = zmin;
    lz = log(z);
    xvlz = tr->rdta->xv * lz;

    rptr  = &(tr->rdta->catrat[1]);
    zzptr = &(zztable[1]);
    zvptr = &(zvtable[1]);

# ifdef Vectorize
#   pragma IVDEP
# endif

    for (i = 1; i <= tr->rdta->categs; i++) {
      ki = *rptr++;
      *zzptr++ = exp(ki *   lz);
      *zvptr++ = exp(ki * xvlz);
    }

    wptr = &(tr->cdta->aliaswgt[0]);
    cptr = &(tr->cdta->patcat[0]);
    log_f = tr->log_f;
    tr->log_f_valid = tr->cdta->endsite;
    sum = 0.0;

# ifdef Vectorize
#   pragma IVDEP
    for (i = 0; i < tr->cdta->endsite; i++) {
      cat   = *cptr++;
      zz[i] = zztable[cat];
      zv[i] = zvtable[cat];
    }

#   pragma IVDEP
    for (i = 0; i < tr->cdta->endsite; i++) {
      fxpa = tr->rdta->freqa * lp->a;
      fxpg = tr->rdta->freqg * lp->g;
      fxpc = tr->rdta->freqc * lp->c;
      fxpt = tr->rdta->freqt * lp->t;
      suma = fxpa * lq->a + fxpc * lq->c + fxpg * lq->g + fxpt * lq->t;
      fxqr = tr->rdta->freqa * lq->a + tr->rdta->freqg * lq->g;
      fxqy = tr->rdta->freqc * lq->c + tr->rdta->freqt * lq->t;
      fxpr = fxpa + fxpg;
      fxpy = fxpc + fxpt;
      sumc = (fxpr + fxpy) * (fxqr + fxqy);
      sumb = fxpr * fxqr * tr->rdta->invfreqr + fxpy * fxqy * tr->rdta->invfreqy;
      suma -= sumb;
      sumb -= sumc;
      term = log(zz[i] * suma + zv[i] * sumb + sumc) + (lp->exp + lq->exp)*log(minlikelihood);
      sum += *wptr++ * term;
      *log_f++ = term;
      lp++;
      lq++;
    }

# else  /*  Not Vectorize  */
    for (i = 0; i < tr->cdta->endsite; i++) {
      cat  = *cptr++;
      zz   = zztable[cat];
      zv   = zvtable[cat];
      fxpa = tr->rdta->freqa * lp->a;
      fxpg = tr->rdta->freqg * lp->g;
      fxpc = tr->rdta->freqc * lp->c;
      fxpt = tr->rdta->freqt * lp->t;
      suma = fxpa * lq->a + fxpc * lq->c + fxpg * lq->g + fxpt * lq->t;
      fxqr = tr->rdta->freqa * lq->a + tr->rdta->freqg * lq->g;
      fxqy = tr->rdta->freqc * lq->c + tr->rdta->freqt * lq->t;
      fxpr = fxpa + fxpg;
      fxpy = fxpc + fxpt;
      sumc = (fxpr + fxpy) * (fxqr + fxqy);
      sumb = fxpr * fxqr * tr->rdta->invfreqr + fxpy * fxqy * tr->rdta->invfreqy;
      suma -= sumb;
      sumb -= sumc;
      term = log(zz * suma + zv * sumb + sumc) + (lp->exp + lq->exp)*log(minlikelihood);
/* printf("evaluate: %le\n", term); */
      sum += *wptr++ * term;
      *log_f++ = term;
      lp++;
      lq++;
    }
# endif  /*  Vectorize or not  */

    tr->likelihood = sum;
    return  sum;
  } /* evaluate */


double makenewz (tree *tr, nodeptr p, nodeptr q, double z0, int maxiter)
  { /* makenewz */
    likelivector *lp, *lq;
    double  *abi, *bci, *sumci, *abptr, *bcptr, *sumcptr;
    double   dlnLidlz, dlnLdlz, d2lnLdlz2, z, zprev, zstep, lz, xvlz,
             ki, suma, sumb, sumc, ab, bc, inv_Li, t1, t2,
             fx1a, fx1c, fx1g, fx1t, fx1r, fx1y, fx2r, fx2y;
    double   zztable[maxcategories+1], zvtable[maxcategories+1],
            *zzptr, *zvptr;
    double  *rptr, *wrptr, *wr2ptr;
    int      cat, *cptr, i, curvatOK;


    while ((! p->x) || (! q->x)) {
      if (! (p->x)) if (! newview(tr, p)) return badZ;
      if (! (q->x)) if (! newview(tr, q)) return badZ;
    }

    lp = p->x->lv;
    lq = q->x->lv;

    { unsigned scratch_size;
      scratch_size = sizeof(double) * tr->cdta->endsite;
      if ((abi   = (double *) Malloc(scratch_size)) &&
          (bci   = (double *) Malloc(scratch_size)) &&
          (sumci = (double *) Malloc(scratch_size))) ;
      else {
        printf("ERROR: makenewz unable to obtain space for arrays\n");
        return badZ;
      }
    }
    abptr   = abi;
    bcptr   = bci;
    sumcptr = sumci;

#   ifdef Vectorize
#     pragma IVDEP
#   endif

    for (i = 0; i < tr->cdta->endsite; i++) {
      fx1a = tr->rdta->freqa * lp->a;
      fx1g = tr->rdta->freqg * lp->g;
      fx1c = tr->rdta->freqc * lp->c;
      fx1t = tr->rdta->freqt * lp->t;
      suma = fx1a * lq->a + fx1c * lq->c + fx1g * lq->g + fx1t * lq->t;
      fx2r = tr->rdta->freqa * lq->a + tr->rdta->freqg * lq->g;
      fx2y = tr->rdta->freqc * lq->c + tr->rdta->freqt * lq->t;
      fx1r = fx1a + fx1g;
      fx1y = fx1c + fx1t;
      *sumcptr++ = sumc = (fx1r + fx1y) * (fx2r + fx2y);
      sumb       = fx1r * fx2r * tr->rdta->invfreqr + fx1y * fx2y * tr->rdta->invfreqy;
      *abptr++   = suma - sumb;
      *bcptr++   = sumb - sumc;
      lp++;
      lq++;
    }

    z = z0;
    do {
      zprev = z;
      zstep = (1.0 - zmax) * z + zmin;
      curvatOK = FALSE;

      do {
        if (z < zmin) z = zmin;
        else if (z > zmax) z = zmax;

        lz    = log(z);
        xvlz  = tr->rdta->xv * lz;
        rptr  = &(tr->rdta->catrat[1]);
        zzptr = &(zztable[1]);
        zvptr = &(zvtable[1]);

#       ifdef Vectorize
#         pragma IVDEP
#       endif

        for (i = 1; i <= tr->rdta->categs; i++) {
          ki = *rptr++;
          *zzptr++ = exp(ki *   lz);
          *zvptr++ = exp(ki * xvlz);
        }

        abptr   = abi;
        bcptr   = bci;
        sumcptr = sumci;
        cptr    = &(tr->cdta->patcat[0]);
        wrptr   = &(tr->cdta->wr[0]);
        wr2ptr  = &(tr->cdta->wr2[0]);
        dlnLdlz = 0.0;                 /*  = d(ln(likelihood))/d(lz) */
        d2lnLdlz2 = 0.0;               /*  = d2(ln(likelihood))/d(lz)2 */

#       ifdef Vectorize
#         pragma IVDEP
#       endif
/* Candidate loop for OpenMP */

        for (i = 0; i < tr->cdta->endsite; i++) {
          cat    = *cptr++;              /*  ratecategory(i) */
          ab     = *abptr++ * zztable[cat];
          bc     = *bcptr++ * zvtable[cat];
          sumc   = *sumcptr++;
          inv_Li = 1.0/(ab + bc + sumc);
          t1     = ab * inv_Li;
          t2     = tr->rdta->xv * bc * inv_Li;
          dlnLidlz   = t1 + t2;
          dlnLdlz   += *wrptr++  * dlnLidlz;
          d2lnLdlz2 += *wr2ptr++ * (t1 + tr->rdta->xv * t2 - dlnLidlz * dlnLidlz);
        }

        if ((d2lnLdlz2 >= 0.0) && (z < zmax))
          zprev = z = 0.37 * z + 0.63;  /*  Bad curvature, shorten branch */
        else
          curvatOK = TRUE;
      } while(!curvatOK);

      if (d2lnLdlz2 < 0.0) {
        z *= exp(-dlnLdlz / d2lnLdlz2);
        if (z < zmin) z = zmin;
        if (z > 0.25 * zprev + 0.75)    /*  Limit steps toward z = 1.0 */
          z = 0.25 * zprev + 0.75;
      }
      if(z>zmax) z=zmax;
    } while((--maxiter>0) && (ABS(z-zprev)>zstep));

    Free(abi);
    Free(bci);
    Free(sumci);

/* printf("makenewz: %le\n", z); */

    return  z;
  } /* makenewz */


boolean update (tree *tr, nodeptr p)
  { /* update */
    nodeptr  q;
    double   z0, z;