seq.c

一个神经网络工具箱

        case '-':
          for (b = A; (long)b <= (long)T; b = (bases)((long)b + 1))
            treenode[i]->x[k][l][(long)b - (long)A] = 1.0;
          break;
        }
      }
    }
  }
}  /* makevalues2 */


void fillin(node *p, node *left, node *rt)
{
  /* sets up for each node in the tree the base sequence
     at that point and counts the changes.  */
  long i, j, k, n, purset, pyrset;
  node *q;

  purset = (1 << (long)A) + (1 << (long)G);
  pyrset = (1 << (long)C) + (1 << (long)T);
  if (!left) {
    memcpy(p->base, rt->base, endsite*sizeof(long));
    memcpy(p->numsteps, rt->numsteps, endsite*sizeof(long));
    q = rt;
  } else if (!rt) {
    memcpy(p->base, left->base, endsite*sizeof(long));
    memcpy(p->numsteps, left->numsteps, endsite*sizeof(long));
    q = left;
  } else {
    for (i = 0; i < endsite; i++) {
      p->base[i] = left->base[i] & rt->base[i];
      p->numsteps[i] = left->numsteps[i] + rt->numsteps[i];
      if (p->base[i] == 0) {
        p->base[i] = left->base[i] | rt->base[i];
        if (transvp) {
          if (!((p->base[i] == purset) || (p->base[i] == pyrset)))
            p->numsteps[i] += weight[i];
        }
        else p->numsteps[i] += weight[i];
      }
    }
    q = rt;
  }
  if (left && rt) n = 2;
  else n = 1;
  for (i = 0; i < endsite; i++)
    for (j = (long)A; j <= (long)O; j++)
      p->numnuc[i][j] = 0;
  for (k = 1; k <= n; k++) {
    if (k == 2) q = left;
    for (i = 0; i < endsite; i++) {
      for (j = (long)A; j <= (long)O; j++) {
        if (q->base[i] & (1 << j))
          p->numnuc[i][j]++;
      }
    }
  }
}  /* fillin */


long getlargest(long *numnuc)
{
  /* find the largest in array numnuc */
  long i, largest;

  largest = 0;
  for (i = (long)A; i <= (long)O; i++)
    if (numnuc[i] > largest)
      largest = numnuc[i];
  return largest;
} /* getlargest */


void multifillin(node *p, node *q, long dnumdesc)
{
  /* sets up for each node in the tree the base sequence
     at that point and counts the changes according to the
     changes in q's base */
  long i, j, b, largest, descsteps, purset, pyrset;

  memcpy(p->oldbase, p->base, endsite*sizeof(long));
  memcpy(p->oldnumsteps, p->numsteps, endsite*sizeof(long));
  purset = (1 << (long)A) + (1 << (long)G);
  pyrset = (1 << (long)C) + (1 << (long)T);
  for (i = 0; i < endsite; i++) {
    descsteps = 0;
    for (j = (long)A; j <= (long)O; j++) {
      b = 1 << j;
      if ((descsteps == 0) && (p->base[i] & b)) 
        descsteps = p->numsteps[i] 
                      - (p->numdesc - dnumdesc - p->numnuc[i][j]) * weight[i];
    }
    if (dnumdesc == -1)
      descsteps -= q->oldnumsteps[i];
    else if (dnumdesc == 0)
      descsteps += (q->numsteps[i] - q->oldnumsteps[i]);
    else
      descsteps += q->numsteps[i];
    if (q->oldbase[i] != q->base[i]) {
      for (j = (long)A; j <= (long)O; j++) {
        b = 1 << j;
        if (transvp) {
          if (b & purset) b = purset;
          if (b & pyrset) b = pyrset;
        }
        if ((q->oldbase[i] & b) && !(q->base[i] & b))
          p->numnuc[i][j]--;
        else if (!(q->oldbase[i] & b) && (q->base[i] & b))
          p->numnuc[i][j]++;
      }
    }
    largest = getlargest(p->numnuc[i]);
    if (q->oldbase[i] != q->base[i]) {
      p->base[i] = 0;
      for (j = (long)A; j <= (long)O; j++) {
        if (p->numnuc[i][j] == largest)
            p->base[i] |= (1 << j);
      }
    }
    p->numsteps[i] = (p->numdesc - largest) * weight[i] + descsteps;
  }
} /* multifillin */


void sumnsteps(node *p, node *left, node *rt, long a, long b)
{
  /* sets up for each node in the tree the base sequence
     at that point and counts the changes. */
  long i;
  long ns, rs, ls, purset, pyrset;

  if (!left) {
    memcpy(p->numsteps, rt->numsteps, endsite*sizeof(long));
    memcpy(p->base, rt->base, endsite*sizeof(long));
  } else if (!rt) {
    memcpy(p->numsteps, left->numsteps, endsite*sizeof(long));
    memcpy(p->base, left->base, endsite*sizeof(long));
  } else  {
    purset = (1 << (long)A) + (1 << (long)G);
    pyrset = (1 << (long)C) + (1 << (long)T);
    for (i = a; i < b; i++) {
      ls = left->base[i];
      rs = rt->base[i];
      ns = ls & rs;
      p->numsteps[i] = left->numsteps[i] + rt->numsteps[i];
      if (ns == 0) {
        ns = ls | rs;
        if (transvp) {
          if (!((ns == purset) || (ns == pyrset)))
            p->numsteps[i] += weight[i];
        }
        else p->numsteps[i] += weight[i];
      }
      p->base[i] = ns;
      }
    }
}  /* sumnsteps */


void sumnsteps2(node *p,node *left,node *rt,long a,long b,long *threshwt)
{
  /* counts the changes at each node.  */
  long i, steps;
  long ns, rs, ls, purset, pyrset;
  long term;

  if (a == 0) p->sumsteps = 0.0;
  if (!left)
    memcpy(p->numsteps, rt->numsteps, endsite*sizeof(long));
  else if (!rt)
    memcpy(p->numsteps, left->numsteps, endsite*sizeof(long));
  else {
    purset = (1 << (long)A) + (1 << (long)G);
    pyrset = (1 << (long)C) + (1 << (long)T);
    for (i = a; i < b; i++) {
      ls = left->base[i];
      rs = rt->base[i];
      ns = ls & rs;
      p->numsteps[i] = left->numsteps[i] + rt->numsteps[i];
      if (ns == 0) {
        ns = ls | rs;
        if (transvp) {
          if (!((ns == purset) || (ns == pyrset)))
            p->numsteps[i] += weight[i];
        }
        else p->numsteps[i] += weight[i];
      }
    }
  }
  for (i = a; i < b; i++) {
    steps = p->numsteps[i];
    if ((long)steps <= threshwt[i])
      term = steps;
    else
      term = threshwt[i];
    p->sumsteps += (double)term;
  }
}  /* sumnsteps2 */


void multisumnsteps(node *p, node *q, long a, long b, long *threshwt)
{
  /* computes the number of steps between p and q */
  long i, j, steps, largest, descsteps, purset, pyrset, b1;
  long term;

  if (a == 0) p->sumsteps = 0.0;
  purset = (1 << (long)A) + (1 << (long)G);
  pyrset = (1 << (long)C) + (1 << (long)T);
  for (i = a; i < b; i++) {
    descsteps = 0;
    for (j = (long)A; j <= (long)O; j++) {
      if ((descsteps == 0) && (p->base[i] & (1 << j))) 
        descsteps = p->numsteps[i] -
                        (p->numdesc - 1 - p->numnuc[i][j]) * weight[i];
    }
    descsteps += q->numsteps[i];
    largest = 0;
    for (j = (long)A; j <= (long)O; j++) {
      b1 = (1 << j);
      if (transvp) {
        if (b1 & purset) b1 = purset;
        if (b1 & pyrset) b1 = pyrset;
      }
      if (q->base[i] & b1)
        p->numnuc[i][j]++;
      if (p->numnuc[i][j] > largest)
        largest = p->numnuc[i][j];
    }
    steps = (p->numdesc - largest) * weight[i] + descsteps;
    if ((long)steps <= threshwt[i])
      term = steps;
    else
      term = threshwt[i];
    p->sumsteps += (double)term;
  }
} /* multisumnsteps */


void multisumnsteps2(node *p)
{
  /* counts the changes at each multi-way node. Sums up
     steps of all descendants */
  long i, j, largest, purset, pyrset, b1;
  node *q;
  baseptr b;

  purset = (1 << (long)A) + (1 << (long)G);
  pyrset = (1 << (long)C) + (1 << (long)T);
  for (i = 0; i < endsite; i++) {
    p->numsteps[i] = 0;
    q = p->next;
    while (q != p) {
      if (q->back) {
        p->numsteps[i] += q->back->numsteps[i];
        b = q->back->base;
        for (j = (long)A; j <= (long)O; j++) {
          b1 = (1 << j);   
          if (transvp) {
            if (b1 & purset) b1 = purset;
            if (b1 & pyrset) b1 = pyrset;
          }
          if (b[i] & b1) p->numnuc[i][j]++;
        }
      }
      q = q->next;
    }
    largest = getlargest(p->numnuc[i]);
    p->base[i] = 0;
    for (j = (long)A; j <= (long)O; j++) {
      if (p->numnuc[i][j] == largest)
        p->base[i] |= (1 << j);
    }
    p->numsteps[i] += ((p->numdesc - largest) * weight[i]);
  }
}  /* multisumnsteps2 */

boolean alltips(node *forknode, node *p)
{
  /* returns true if all descendants of forknode except p are tips; 
     false otherwise.  */
  node *q, *r;
  boolean tips;

  tips = true;
  r = forknode;
  q = forknode->next;
  do {
    if (q->back && q->back != p && !q->back->tip)
      tips = false;
    q = q->next;
  } while (tips && q != r);
  return tips;
} /* alltips */


void gdispose(node *p, node **grbg, pointarray treenode)
{
  /* go through tree throwing away nodes */
  node *q, *r;

  p->back = NULL;
  if (p->tip)
    return;
  treenode[p->index - 1] = NULL;
  q = p->next;
  while (q != p) {
    gdispose(q->back, grbg, treenode);
    q->back = NULL;
    r = q;
    q = q->next;
    chucktreenode(grbg, r);
  }
  chucktreenode(grbg, q);
}  /* gdispose */


void preorder(node *p, node *r, node *root, node *removing, node *adding,
                        node *changing, long dnumdesc)
{
  /* recompute number of steps in preorder taking both ancestoral and
     descendent steps into account. removing points to a node being 
     removed, if any */
  node *q, *p1, *p2;

  if (p && !p->tip && p != adding) {
    q = p;
    do {
      if (p->back != r) {
        if (p->numdesc > 2) {
          if (changing)
            multifillin (p, r, dnumdesc);
          else
            multifillin (p, r, 0);
        } else {
          p1 = p->next;
          if (!removing)
            while (!p1->back)
              p1 = p1->next;
          else
            while (!p1->back || p1->back == removing)
              p1 = p1->next;
          p2 = p1->next;
          if (!removing)
            while (!p2->back)
              p2 = p2->next;
          else
            while (!p2->back || p2->back == removing)
              p2 = p2->next;
          p1 = p1->back;
          p2 = p2->back;
          if (p->back == p1) p1 = NULL;
          else if (p->back == p2) p2 = NULL;
          memcpy(p->oldbase, p->base, endsite*sizeof(long));
          memcpy(p->oldnumsteps, p->numsteps, endsite*sizeof(long));
          fillin(p, p1, p2);
        }
      }
      p = p->next;
    } while (p != q);
    q = p;
    do {
      preorder(p->next->back, p->next, root, removing, adding, NULL, 0);
      p = p->next;
    } while (p->next != q);
  }
} /* preorder */

void updatenumdesc(node *p, node *root, long n)
{
  /* set p's numdesc to n.  If p is the root, numdesc of p's
  descendants are set to n-1. */
  node *q;

  q = p;
  if (p == root && n > 0) {
    p->numdesc = n;
    n--;
    q = q->next;
  }
  do {
    q->numdesc = n;
    q = q->next;
  } while (q != p);
}  /* updatenumdesc */


void add(node *below,node *newtip,node *newfork,node **root,
        boolean recompute,pointarray treenode,node **grbg,long *zeros)
{
  /* inserts the nodes newfork and its left descendant, newtip,
     to the tree.  below becomes newfork's right descendant.
     if newfork is NULL, newtip is added as below's sibling */
  /* used in dnacomp & dnapars */
  node *p;

  if (below != treenode[below->index - 1])
    below = treenode[below->index - 1];
  if (newfork) {
    if (below->back != NULL)
      below->back->back = newfork;
    newfork->back = below->back;
    below->back = newfork->next->next;
    newfork->next->next->back = below;
    newfork->next->back = newtip;
    newtip->back = newfork->next;
    if (*root == below)
      *root = newfork;
    updatenumdesc(newfork, *root, 2);
  } else {
    gnutreenode(grbg, &p, below->index, endsite, zeros);
    p->back = newtip;
    newtip->back = p;
    p->next = below->next;
    below->next = p;
    updatenumdesc(below, *root, below->numdesc + 1);
  }
  if (!newtip->tip)
    updatenumdesc(newtip, *root, newtip->numdesc);
  (*root)->back = NULL;
  if (!recompute)
    return;
  if (!newfork) {
    memcpy(newtip->back->base, below->base, endsite*sizeof(long));
    memcpy(newtip->back->numsteps, below->numsteps, endsite*sizeof(long));
    memcpy(newtip->back->numnuc, below->numnuc, endsite*sizeof(nucarray));
    if (below != *root) {
      memcpy(below->back->oldbase, zeros, endsite*sizeof(long));
      memcpy(below->back->oldnumsteps, zeros, endsite*sizeof(long));
      multifillin(newtip->back, below->back, 1);
    }
    if (!newtip->tip) {
      memcpy(newtip->back->oldbase, zeros, endsite*sizeof(long));
      memcpy(newtip->back->oldnumsteps, zeros, endsite*sizeof(long));
      preorder(newtip, newtip->back, *root, NULL, NULL, below, 1);
    }
    memcpy(newtip->oldbase, zeros, endsite*sizeof(long));
    memcpy(newtip->oldnumsteps, zeros, endsite*sizeof(long));
    preorder(below, newtip, *root, NULL, newtip, below, 1);
    if (below != *root)
      preorder(below->back, below, *root, NULL, NULL, NULL, 0);
  } else {
    fillin(newtip->back, newtip->back->next->back,
             newtip->back->next->next->back);
    if (!newtip->tip) {
      memcpy(newtip->back->oldbase, zeros, endsite*sizeof(long));
      memcpy(newtip->back->oldnumsteps, zeros, endsite*sizeof(long));
      preorder(newtip, newtip->back, *root, NULL, NULL, newfork, 1);
    }
    if (newfork != *root) {
      memcpy(below->back->base, newfork->back->base, endsite*sizeof(long));
      memcpy(below->back->numsteps, newfork->back->numsteps, endsite*sizeof(long));
      preorder(newfork, newtip, *root, NULL, newtip, NULL, 0);
    } else {
      fillin(below->back, newtip, NULL);
      fillin(newfork, newtip, below);
      memcpy(below->back->oldbase, zeros, endsite*sizeof(long));
      memcpy(below->back->oldnumsteps, zeros, endsite*sizeof(long));
      preorder(below, below->back, *root, NULL, NULL, newfork, 1);
    }
    if (newfork != *root) {
      memcpy(newfork->oldbase, below->base, endsite*sizeof(long));
      memcpy(newfork->oldnumsteps, below->numsteps, endsite*sizeof(long));
      preorder(newfork->back, newfork, *root, NULL, NULL, NULL, 0);
    }
  }
}  /* add */