print_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  <math.h>
#include  "fastDNAml_types.h"
#include  "fastDNAml_funcs.h"



void coordinates (tree *tr, nodeptr p, double lengthsum, drawdata *tdptr)
  { /* coordinates */
    /* establishes coordinates of nodes */
    double  x, z;
    nodeptr  q, first, last;

    if (p->tip) {
      p->xcoord = NINT(over * lengthsum);
      p->ymax = p->ymin = p->ycoord = tdptr->tipy;
      tdptr->tipy += down;
      if (lengthsum > tdptr->tipmax) tdptr->tipmax = lengthsum;
      }

    else {
      q = p->next;
      do {
        z = q->z;
        if (z < zmin) z = zmin;
        x = lengthsum - tr->rdta->fracchange * log(z);
        coordinates(tr, q->back, x, tdptr);
        q = q->next;
        } while (p == tr->start->back ? q != p->next : q != p);

      first = p->next->back;
      q = p;
      while (q->next != p) q = q->next;
      last = q->back;
      p->xcoord = NINT(over * lengthsum);
      p->ycoord = (first->ycoord + last->ycoord)/2;
      p->ymin = first->ymin;
      p->ymax = last->ymax;
      }
  } /* coordinates */


void drawline (tree *tr, int i, double scale)
    /* draws one row of the tree diagram by moving up tree */
    /* Modified to handle 1000 taxa, October 16, 1991 */
  { /* drawline */
    nodeptr  p, q, r, first, last;
    int  n, j, k, l, extra;
    boolean  done;

    p = q = tr->start->back;
    extra = 0;

    if (i == p->ycoord) {
      k = q->number - tr->mxtips;
      for (j = k; j < 1000; j *= 10) putchar('-');
      printf("%d", k);
      extra = 1;
      }
    else printf("   ");

    do {
      if (! p->tip) {
        r = p->next;
        done = FALSE;
        do {
          if ((i >= r->back->ymin) && (i <= r->back->ymax)) {
            q = r->back;
            done = TRUE;
            }
          r = r->next;
          } while (! done && (p == tr->start->back ? r != p->next : r != p));

        first = p->next->back;
        r = p;
        while (r->next != p) r = r->next;
        last = r->back;
        if (p == tr->start->back) last = p->back;
        }

      done = (p->tip) || (p == q);
      n = NINT(scale*(q->xcoord - p->xcoord));
      if ((n < 3) && (! q->tip)) n = 3;
      n -= extra;
      extra = 0;

      if ((q->ycoord == i) && (! done)) {
        if (p->ycoord != q->ycoord) putchar('+');
        else                        putchar('-');

        if (! q->tip) {
          k = q->number - tr->mxtips;
          l = n - 3;
          for (j = k; j < 100; j *= 10)  l++;
          for (j = 1; j <= l; j++) putchar('-');
          printf("%d", k);
          extra = 1;
          }
        else for (j = 1; j <= n-1; j++) putchar('-');
        }

      else if (! p->tip) {
        if ((last->ycoord > i) && (first->ycoord < i) && (i != p->ycoord)) {
          putchar('!');
          for (j = 1; j <= n-1; j++) putchar(' ');
          }
        else for (j = 1; j <= n; j++) putchar(' ');
        }

      else
        for (j = 1; j <= n; j++) putchar(' ');

      p = q;
      } while (! done);

    if ((p->ycoord == i) && p->tip) {
      printf(" %s", p->name);
      }

    putchar('\n');
  } /* drawline */


void printTree (tree *tr, analdef *adef)
    /* prints out diagram of the tree */
  { /* printTree */
    drawdata  tipdata;
    double  scale;
    int  i, imax;

    if (adef->trprint) {
      putchar('\n');
      tipdata.tipy = 1;
      tipdata.tipmax = 0.0;
      coordinates(tr, tr->start->back, (double) 0.0, & tipdata);
      scale = 1.0 / tipdata.tipmax;
      imax = tipdata.tipy - down;
      for (i = 1; i <= imax; i++)  drawline(tr, i, scale);
      printf("\nRemember: ");
      if (adef->root) printf("(although rooted by outgroup) ");
      printf("this is an unrooted tree!\n\n");
      }
  } /* printTree */


double sigma (tree *tr, nodeptr p, double *sumlrptr)
    /* compute standard deviation */
  { /* sigma */
    likelivector *lp, *lq;
    double  slope, sum, sumlr, z, zv, zz, lz,
            rat, suma, sumb, sumc, d2, d, li, temp, abzz, bczv, t3,
            fxpa, fxpc, fxpg, fxpt, fxpr, fxpy, fxqr, fxqy, w;
    double  *rptr;
    nodeptr  q;
    int  i, *wptr;

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

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

    z = p->z;
    if (z < zmin) z = zmin;
    lz = log(z);

    wptr = &(tr->cdta->aliaswgt[0]);
    rptr = &(tr->cdta->patrat[0]);
    sum = sumlr = slope = 0.0;

#   ifdef Vectorize
#     pragma IVDEP
#   endif

    for (i = 0; i < tr->cdta->endsite; i++) {
      rat  = *rptr++;
      zz   = exp(rat                * lz);
      zv   = exp(rat * tr->rdta->xv * lz);

      fxpa = tr->rdta->freqa * lp->a;
      fxpg = tr->rdta->freqg * lp->g;
      fxpc = tr->rdta->freqc * lp->c;
      fxpt = tr->rdta->freqt * lp->t;
      fxpr = fxpa + fxpg;
      fxpy = fxpc + fxpt;
      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;
      sumc = (fxpr + fxpy) * (fxqr + fxqy);
      sumb = fxpr * fxqr * tr->rdta->invfreqr + fxpy * fxqy * tr->rdta->invfreqy;
      abzz = zz * (suma - sumb);
      bczv = zv * (sumb - sumc);
      li = sumc + abzz + bczv;
      t3 = tr->rdta->xv * bczv;
      d  = abzz + t3;
      d2 = rat * (abzz*(rat-1.0) + t3*(rat * tr->rdta->xv - 1.0));

      temp = rat * d / li;
      w = *wptr++;
      slope += w *  temp;
      sum   += w * (temp * temp - d2/li);
      sumlr += w * log(li / (suma + 1.0E-300));
      lp++;
      lq++;
      }

    *sumlrptr = sumlr;
    return (sum > 1.0E-300) ? z*(-slope + sqrt(slope*slope + 3.841*sum))/sum
                            : 1.0;
  } /* sigma */


void describe (tree *tr, nodeptr p)
    /* print out information for one branch */
  { /* describe */
    double   z, s, sumlr;
    nodeptr  q;
    char    *nameptr;
    int      k, ch;

    q = p->back;
    printf("%4d          ", q->number - tr->mxtips);
    if (p->tip) {
      nameptr = p->name;
      k = nmlngth;
      while (ch = *nameptr++) {putchar(ch); k--;}
      while (--k >= 0) putchar(' ');
      }
    else {
      printf("%4d", p->number - tr->mxtips);
      for (k = 4; k < nmlngth; k++) putchar(' ');
      }

    z = q->z;
    if (z <= zmin) printf("    infinity");
    else printf("%15.5f", -log(z) * tr->rdta->fracchange);

    s = sigma(tr, q, & sumlr);
    printf("     (");
    if (z + s >= zmax) printf("     zero");
    else printf("%9.5f", (double) -log(z + s) * tr->rdta->fracchange);
    putchar(',');
    if (z - s <= zmin) printf("    infinity");
    else printf("%12.5f", (double) -log(z - s) * tr->rdta->fracchange);
    putchar(')');

    if      (sumlr > 2.995 ) printf(" **");
    else if (sumlr > 1.9205) printf(" *");
    putchar('\n');

    if (! p->tip) {
      describe(tr, p->next->back);
      describe(tr, p->next->next->back);
      }
  } /* describe */


void summarize (tree *tr)
    /* print out branch length information and node numbers */
  { /* summarize */
    printf("Ln Likelihood =%14.5f\n", tr->likelihood);
    putchar('\n');
    printf(" Between        And             Length");
    printf("      Approx. Confidence Limits\n");
    printf(" -------        ---             ------");
    printf("      ------- ---------- ------\n");

    describe(tr, tr->start->back->next->back);
    describe(tr, tr->start->back->next->next->back);
    describe(tr, tr->start);
    putchar('\n');
    printf("     *  = significantly positive, P < 0.05\n");
    printf("     ** = significantly positive, P < 0.01\n\n\n");
  } /* summarize */