proml.c

一个神经网络工具箱

  (*p)->next->back = NULL;
  (*p)->next->next->back = NULL;
  inittrav((*q));
  inittrav((*q)->back);
  i = 1;
  while (i <= smoothings) {
    smooth(*q);
    if (smoothit)
      smooth((*q)->back);
    i++;
  }
}  /* proml_re_move */


void insert_(node *p, node *q, boolean dooinit)
{
  /* Insert q near p */
  long i, j, num_sibs;
  node *r, *sib_ptr;

  r = p->next->next;
  hookup(r, q->back);
  hookup(p->next, q);
  q->v = 0.5 * q->v;
  q->back->v = q->v;
  r->v = q->v;
  r->back->v = r->v;
  p->initialized = false;
  p->next->initialized = false;
  p->next->next->initialized = false;
  if (dooinit) {
    inittrav(p);
    inittrav(q);
    inittrav(q->back);
  }
  i = 1;
  while (i <= smoothings) {
    smooth(p);
    if (!smoothit) {
      if (!p->tip) {
        num_sibs = count_sibs(p);
        sib_ptr  = p;
        for (j=0; j < num_sibs; j++) {
          smooth(sib_ptr->next->back);
          sib_ptr = sib_ptr->next;
        }
      }
    }
    else
      smooth(p->back);
    i++;
  }
}  /* insert_ */


void addtraverse(node *p, node *q, boolean contin)
{
  /* try adding p at q, proceed recursively through tree */
  long i, num_sibs;
  double like, vsave = 0;
  node *qback = NULL, *sib_ptr;

  if (!smoothit) {
    vsave = q->v;
    qback = q->back;
  }
  insert_(p, q, false);
  like = prot_evaluate(p, false);
  if (like > bestyet || bestyet == UNDEFINED) {
    bestyet = like;
    if (smoothit)
      promlcopy(&curtree, &bestree, nonodes2, rcategs);
    else
      qwhere = q;
    succeeded = true;
  }
  if (smoothit)
    promlcopy(&priortree, &curtree, nonodes2, rcategs);
  else {
    hookup (q, qback);
    q->v = vsave;
    q->back->v = vsave;
    curtree.likelihood = bestyet;
  }
  if (!q->tip && contin) {
    num_sibs = count_sibs(q);
    if (q == curtree.start)
      num_sibs++;
    sib_ptr  = q;
    for (i=0; i < num_sibs; i++) {
      addtraverse(p, sib_ptr->next->back, contin);
      sib_ptr = sib_ptr->next;
    }
  }
}  /* addtraverse */


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
      promlcopy(&curtree, &bestree, nonodes2, rcategs);
    proml_re_move(&r, &q);
    if (smoothit)
      promlcopy(&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)
      promlcopy(&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;
        promlcopy(&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 proml_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 = q->v;
    if (xx > 100.0)
      xx = 100.0;
    proml_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;
}  /* proml_coordinates */


void proml_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;
  proml_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);
}  /* proml_printree */


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

  prot_slopecurv(p, p->v, &like, &slope, &curv);
  tt = p->v;
  p->v = epsilon;
  p->back->v = epsilon;
  aa = prot_evaluate(p, false);
  p->v = tt;
  p->back->v = tt;
  (*sumlr) = prot_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)
    prot_nuview(p);
  if (!p->back->tip && !p->back->initialized)
    prot_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);
  if (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);
      fprintf(outfile, ",%12.5f", sigma1);
      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 prot_reconstr(node *p, long n)
{
  /* reconstruct and print out acid at site n+1 at node p */
  long i, j, k, first, num_sibs = 0;
  double f, sum, xx[20];
  node *q = NULL;

  if (p->tip)
    putc(y[p->index-1][n], outfile);
  else {
    num_sibs = count_sibs(p);
    if ((ally[n] == 0) || (location[ally[n]-1] == 0))
      putc('.', outfile);
    else {
      j = location[ally[n]-1] - 1;
      sum = 0;
      for (i = 0; i <= 19; i++) {
        f = p->protx[j][mx-1][i];
        if (!p->tip) {
          q = p;
          for (k = 0; k < num_sibs; k++) {
            q = q->next;
            f *= q->protx[j][mx-1][i];
          }
        }
        f = sqrt(f);
        xx[i] = f * freqaa[i];
        sum += xx[i];
      }
      for (i = 0; i <= 19; i++)
        xx[i] /= sum;
      first = 0;
      for (i = 0; i <= 19; i++)
        if (xx[i] > xx[first])
          first = i;
      if (xx[first] > 0.95)
        putc(aachar[first], outfile);
      else
        putc(tolower(aachar[first]), outfile);
      if (rctgry && rcategs > 1)
        mx = mp[n][mx - 1];
      else
        mx = 1;
    }
  }
} /* prot_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);
    prot_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))
    fprintf(outfile, "      Approx. Confidence Limits");
  fprintf(outfile, "\n");
  fprintf(outfile, " -------        ---            ------");
  if (!(usertree && lngths))
    fprintf(o