dnaml.c

一个神经网络工具箱

void rearrange(node *p, node *pp)
{
  /* rearranges the tree, globally or locally moving pp around near p */
  long i, num_sibs;
  double v3 = 0, v4 = 0, v5 = 0;
  node *q, *r, *sib_ptr;

  if (!p->tip && !p->back->tip) {
    curtree.likelihood = bestyet;
    if (p->back->next != pp)
      r = p->back->next;
    else
      r = p->back->next->next;
    /* assumes bifurcations? */
    if (!smoothit) {
      v3 = r->v;
      v4 = r->next->v;
      v5 = r->next->next->v;
    }
    else
      dnamlcopy(&curtree, &bestree, nonodes2, rcategs);
    dnaml_re_move(&r, &q);
    if (smoothit)
      dnamlcopy(&curtree, &priortree, nonodes2, rcategs);
    else
      qwhere = q;
    num_sibs = count_sibs (p);
    sib_ptr  = p;
    for (i=0; i < num_sibs; i++) {
      sib_ptr = sib_ptr->next;
      addtraverse(r, sib_ptr->back, (boolean)(global && (nextsp == spp)));
    }
    if (global && nextsp == spp && !succeeded) {
      p = p->back;
      if (!p->tip) {
        num_sibs = count_sibs (p);
        sib_ptr  = p;
        for (i=0; i < num_sibs; i++) {
          sib_ptr = sib_ptr->next;
          addtraverse(r, sib_ptr->back,(boolean)(global && (nextsp == spp)));
        }
      }
      p = p->back;
    }
    if (smoothit)
      dnamlcopy(&bestree, &curtree, nonodes2, rcategs);
    else {
      insert_(r, qwhere, true);
      if (qwhere == q) {
        r->v = v3;
        r->back->v = v3;
        r->next->v = v4;
        r->next->back->v = v4;
        r->next->next->v = v5;
        r->next->next->back->v = v5;
        curtree.likelihood = bestyet;
      }
      else {
        smoothit = true;
        for (i = 1; i<=smoothings; i++) {
          smooth (r);
          smooth (r->back);
        }
        smoothit = false;
        dnamlcopy(&curtree, &bestree, nonodes2, rcategs);
      }
    }
    if (global && nextsp == spp && progress) {
      putchar('.');
      fflush(stdout);
    }
  }
  if (!p->tip) {
    num_sibs = count_sibs (p);
    if (p == curtree.start)
      num_sibs++;
    sib_ptr  = p;
    for (i=0; i < num_sibs; i++) {
      sib_ptr = sib_ptr->next;
      rearrange(sib_ptr->back, p);
    }
  }
}  /* rearrange */


void initdnamlnode(node **p, node **grbg, node *q, long len, long nodei,
                        long *ntips, long *parens, initops whichinit,
                        pointarray treenode, pointarray nodep, Char *str, Char *ch,
                        FILE *intree)
{
  /* initializes a node */
  boolean minusread;
  double valyew, divisor;

  switch (whichinit) {
  case bottom:
    gnu(grbg, p);
    (*p)->index = nodei;
    (*p)->tip = false;
    malloc_pheno((*p), endsite, rcategs);
    nodep[(*p)->index - 1] = (*p);
    break;
  case nonbottom:
    gnu(grbg, p);
    malloc_pheno(*p, endsite, rcategs);
    (*p)->index = nodei;
    break;
  case tip:
    match_names_to_data (str, nodep, p, spp);
    break;
  case iter:
    (*p)->initialized = false;
    (*p)->v = initialv;
    (*p)->iter = true;
    if ((*p)->back != NULL)
      (*p)->back->iter = true;
    break;
  case length:
    processlength(&valyew, &divisor, ch, &minusread, intree, parens);
    (*p)->v = valyew / divisor / fracchange;
    (*p)->iter = false;
    if ((*p)->back != NULL) {
      (*p)->back->v = (*p)->v;
      (*p)->back->iter = false;
    }
    break;
  case hsnolength:
    haslengths = false;
    break;
  default:        /* cases hslength, treewt, unittrwt */
    break;        /* should never occur                                */
  }
} /* initdnamlnode */


void dnaml_coordinates(node *p, double lengthsum, long *tipy,
                        double *tipmax)
{
  /* establishes coordinates of nodes */
  node *q, *first, *last;
  double xx;

  if (p->tip) {
    p->xcoord = (long)(over * lengthsum + 0.5);
    p->ycoord = (*tipy);
    p->ymin = (*tipy);
    p->ymax = (*tipy);
    (*tipy) += down;
    if (lengthsum > (*tipmax))
      (*tipmax) = lengthsum;
    return;
  }
  q = p->next;
  do {
    xx = fracchange * q->v;
    if (xx > 100.0)
      xx = 100.0;
    dnaml_coordinates(q->back, lengthsum + xx, tipy,tipmax);
    q = q->next;
  } while ((p == curtree.start || p != q) &&
           (p != curtree.start || p->next != q));
  first = p->next->back;
  q = p;
  while (q->next != p)
    q = q->next;
  last = q->back;
  p->xcoord = (long)(over * lengthsum + 0.5);
  if (p == curtree.start)
    p->ycoord = p->next->next->back->ycoord;
  else
    p->ycoord = (first->ycoord + last->ycoord) / 2;
  p->ymin = first->ymin;
  p->ymax = last->ymax;
}  /* dnaml_coordinates */


void dnaml_printree()
{
  /* prints out diagram of the tree2 */
  long tipy;
  double scale, tipmax;
  long i;

  if (!treeprint)
    return;
  putc('\n', outfile);
  tipy = 1;
  tipmax = 0.0;
  dnaml_coordinates(curtree.start, 0.0, &tipy, &tipmax);
  scale = 1.0 / (long)(tipmax + 1.000);
  for (i = 1; i <= (tipy - down); i++)
    drawline2(i, scale, curtree);
  putc('\n', outfile);
}  /* dnaml_printree */


void sigma(node *p, double *sumlr, double *s1, double *s2)
{
  /* compute standard deviation */
  double tt, aa, like, slope, curv;

  slopecurv (p, p->v, &like, &slope, &curv);
  tt = p->v;
  p->v = epsilon;
  p->back->v = epsilon;
  aa = evaluate(p, false);
  p->v = tt;
  p->back->v = tt;
  (*sumlr) = evaluate(p, false) - aa;
  if (curv < -epsilon) {
    (*s1) = p->v + (-slope - sqrt(slope * slope -  3.841 * curv)) / curv;
    (*s2) = p->v + (-slope + sqrt(slope * slope -  3.841 * curv)) / curv;
  }
  else {
    (*s1) = -1.0;
    (*s2) = -1.0;
  }
}  /* sigma */


void describe(node *p)
{
  /* print out information for one branch */
  long i, num_sibs;
  node *q, *sib_ptr;
  double sumlr, sigma1, sigma2;

  if (!p->tip && !p->initialized)
    nuview(p);
  if (!p->back->tip && !p->back->initialized)
    nuview(p->back);
  q = p->back;
  if (q->tip) {
    fprintf(outfile, " ");
    for (i = 0; i < nmlngth; i++)
      putc(nayme[q->index-1][i], outfile);
    fprintf(outfile, "    ");
  } else
    fprintf(outfile, "  %4ld          ", q->index - spp);
  if (p->tip) {
    for (i = 0; i < nmlngth; i++)
      putc(nayme[p->index-1][i], outfile);
  } else
    fprintf(outfile, "%4ld      ", p->index - spp);
  fprintf(outfile, "%15.5f", q->v * fracchange);
  if (!usertree || (usertree && !lngths) || p->iter) {
    sigma(q, &sumlr, &sigma1, &sigma2);
    if (sigma1 <= sigma2)
      fprintf(outfile, "     (     zero,    infinity)");
    else {
      fprintf(outfile, "     (");
      if (sigma2 <= 0.0)
        fprintf(outfile, "     zero");
      else
        fprintf(outfile, "%9.5f", sigma2 * fracchange);
      fprintf(outfile, ",%12.5f", sigma1 * fracchange);
      putc(')', outfile);
      }
    if (sumlr > 1.9205)
      fprintf(outfile, " *");
    if (sumlr > 2.995)
      putc('*', outfile);
    }
  putc('\n', outfile);
  if (!p->tip) {
    num_sibs = count_sibs (p);
    sib_ptr  = p;
    for (i=0; i < num_sibs; i++) {
      sib_ptr = sib_ptr->next;
      describe(sib_ptr->back);
    }
  }
}  /* describe */


void reconstr(node *p, long n)
{
  /* reconstruct and print out base at site n+1 at node p */
  long i, j, k, m, first, second, num_sibs;
  double f, sum, xx[4];
  node *q;

  if ((ally[n] == 0) || (location[ally[n]-1] == 0))
    putc('.', outfile);
  else {
    j = location[ally[n]-1] - 1;
    for (i = 0; i < 4; i++) {
      f = p->x[j][mx-1][i];
      num_sibs = count_sibs(p);
      q = p;
      for (k = 0; k < num_sibs; k++) {
        q = q->next;
        f *= q->x[j][mx-1][i];
      }
      f = sqrt(f);
      xx[i] = f;
    }
    xx[0] *= freqa;
    xx[1] *= freqc;
    xx[2] *= freqg;
    xx[3] *= freqt;
    sum = xx[0]+xx[1]+xx[2]+xx[3];
    for (i = 0; i < 4; i++)
      xx[i] /= sum;
    first = 0;
    for (i = 1; i < 4; i++)
      if (xx [i] > xx[first])
        first = i;
    if (first == 0)
      second = 1;
    else
      second = 0;
    for (i = 0; i < 4; i++)
      if ((i != first) && (xx[i] > xx[second]))
        second = i;
    m = 1 << first;
    if (xx[first] < 0.4999995)
      m = m + (1 << second);
    if (xx[first] > 0.95)
      putc(toupper(basechar[m - 1]), outfile);
    else
      putc(basechar[m - 1], outfile);
    if (rctgry && rcategs > 1)
      mx = mp[n][mx - 1];    
    else
      mx = 1;
  }
} /* reconstr */


void rectrav(node *p, long m, long n)
{
  /* print out segment of reconstructed sequence for one branch */
  long i;

  putc(' ', outfile);
  if (p->tip) {
    for (i = 0; i < nmlngth; i++)
      putc(nayme[p->index-1][i], outfile);
  } else
    fprintf(outfile, "%4ld      ", p->index - spp);
  fprintf(outfile, "  ");
  mx = mx0;
  for (i = m; i <= n; i++) {
    if ((i % 10 == 0) && (i != m))
      putc(' ', outfile);
    if (p->tip)
      putc(y[p->index-1][i], outfile);
    else
      reconstr(p, i);
  }
  putc('\n', outfile);
  if (!p->tip) {
    rectrav(p->next->back, m, n);
    rectrav(p->next->next->back, m, n);
  }
  mx1 = mx;
}  /* rectrav */


void summarize()
{
  /* print out branch length information and node numbers */
  long i, j, mm, num_sibs;
  double mode, sum;
  double like[maxcategs], nulike[maxcategs];
  double **marginal;
  node   *sib_ptr;

  if (!treeprint)
    return;
  fprintf(outfile, "\nremember: ");
  if (outgropt)
    fprintf(outfile, "(although rooted by outgroup) ");
  fprintf(outfile, "this is an unrooted tree!\n\n");
  fprintf(outfile, "Ln Likelihood = %11.5f\n", curtree.likelihood);
  fprintf(outfile, "\n Between        And            Length");
  if (!(usertree && lngths && haslengths))
    fprintf(outfile, "      Approx. Confidence Limits");
  fprintf(outfile, "\n");
  fprintf(outfile, " -------        ---            ------");
  if (!(usertree && lngths && haslengths))
    fprintf(outfile, "      ------- ---------- ------");
  fprintf(outfile, "\n\n");
  for (i = spp; i < nonodes2; i++) {
    /* So this works with arbitrary multifurcations */
    if (curtree.nodep[i]) {
      num_sibs = count_sibs (curtree.nodep[i]);
      sib_ptr  = curtree.nodep[i];
      for (j = 0; j < num_sibs; j++) {
        sib_ptr->initialized = false;
        sib_ptr              = sib_ptr->next;
      }
    }
  }

  describe(curtree.start->back);

  /* So this works with arbitrary multifurcations */
  num_sibs = count_sibs (curtree.start);
  sib_ptr  = curtree.start;
  for (i=0; i < num_sibs; i++) {
    sib_ptr = sib_ptr->next;
    describe(sib_ptr->back);
  }

  fprintf(outfile, "\n");
  if (!(usertree && lngths && haslengths)) {
    fprintf(outfile, "     *  = significantly positive, P < 0.05\n");
    fprintf(outfile, "     ** = significantly positive, P < 0.01\n\n");
  }
  dummy = evaluate(curtree.start, false);
  if (rctgry && rcategs > 1) {
    for (i = 0; i < rcategs; i++)
      like[i] = 1.0;
    for (i = sites - 1; i >= 0; i--) {
      sum = 0.0;
      for (j = 0; j < rcategs; j++) {
        nulike[j] = (1.0 - lambda + lambda * probcat[j]) * like[j];
        mp[i][j] = j + 1;
        for (k = 1; k <= rcategs; k++) {
          if (k != j + 1) {
            if (lambda * probcat[k - 1] * like[k - 1] > nulike[j]) {
              nulike[j] = lambda * probcat[k - 1] * like[k - 1];
              mp[i][j] = k;
            }
          }
        }
        if ((ally[i] > 0) && (location[ally[i]-1] > 0))
          nulike[j] *= contribution[location[ally[i] - 1] - 1][j];
        sum += nulike[j];
      }
      for (j = 0; j < rcategs; j++)
        nulike[j] /= sum;
      memcpy(like, nulike, rcategs * sizeof(double));
    }
    mode = 0.0;
    mx = 1;
    for (i = 1; i <= rcategs; i++) {
      if (probcat[i - 1] * like[i - 1] > mode) {
        mx = i;
        mode = probcat[i - 1] * like[i - 1];
      }
    }
    mx0 = mx;
    fprintf(outfile,
 "Combination of categories that contributes the most to the likelihood:\n\n");
    for (i = 1; i <= nmlngth + 3; i++)
      putc(' ', outfile);
    for (i = 1; i <= sites; i++) {
      fprintf(outfile, "%ld", mx);
      if (i % 10 == 0)
        putc(' ', outfile);
      if (i % 60 == 0 && i != sites) {
        putc('\n', outfile);
        for (j = 1; j <= nmlngth + 3; j++)
          putc(' ', outfile);
      }
      mx = mp[i - 1][mx - 1];
    }
    fprintf(outfile, "\n\n");
    marginal = (double **) Malloc( sites*sizeof(double *));
    for (i = 0; i < sites; i++)
      marginal[i] = (double *) Malloc( rcategs*sizeof(double));
    for (i = 0; i < rcategs; i++)
      like[i] = 1.0;
    for (i = sites - 1; i >= 0; i--) {