seq.c

来自「一个神经网络工具箱」· C语言 代码 · 共 2,339 行 · 第 1/5 页

void findbelow(node **below, node *item, node *fork)
{
  /* decide which of fork's binary children is below */

  if (fork->next->back == item)
    *below = fork->next->next->back;
  else
    *below = fork->next->back;
} /* findbelow */


void re_move(node *item, node **fork, node **root, boolean recompute,
                        pointarray treenode, node **grbg, long *zeros)
{
  /* removes nodes item and its ancestor, fork, from the tree.
     the new descendant of fork's ancestor is made to be
     fork's second descendant (other than item).  Also
     returns pointers to the deleted nodes, item and fork.
     If item belongs to a node with more than 2 descendants,
     fork will not be deleted */
  /* used in dnacomp & dnapars */
  node *p, *q, *other = NULL, *otherback = NULL;

  if (item->back == NULL) {
    *fork = NULL;
    return;
  }
  *fork = treenode[item->back->index - 1];
  if ((*fork)->numdesc == 2) {
    updatenumdesc(*fork, *root, 0);
    findbelow(&other, item, *fork);
    otherback = other->back;
    if (*root == *fork) {
      *root = other;
      if (!other->tip)
        updatenumdesc(other, *root, other->numdesc);
    }
    p = item->back->next->back;
    q = item->back->next->next->back;
    if (p != NULL)
      p->back = q;
    if (q != NULL)
      q->back = p;
    (*fork)->back = NULL;
    p = (*fork)->next;
    while (p != *fork) {
      p->back = NULL;
      p = p->next;
    }
  } else {
    updatenumdesc(*fork, *root, (*fork)->numdesc - 1);
    p = *fork;
    while (p->next != item->back)
      p = p->next;
    p->next = item->back->next;
  }
  if (!item->tip) {
    updatenumdesc(item, item, item->numdesc);
    if (recompute) {
      memcpy(item->back->oldbase, item->back->base, endsite*sizeof(long));
      memcpy(item->back->oldnumsteps, item->back->numsteps, endsite*sizeof(long));
      memcpy(item->back->base, zeros, endsite*sizeof(long));
      memcpy(item->back->numsteps, zeros, endsite*sizeof(long));
      preorder(item, item->back, *root, item->back, NULL, item, -1);
    }
  }
  if ((*fork)->numdesc >= 2)
    chucktreenode(grbg, item->back);
  item->back = NULL;
  if (!recompute)
    return;
  if ((*fork)->numdesc == 0) {
    memcpy(otherback->oldbase, otherback->base, endsite*sizeof(long));  
    memcpy(otherback->oldnumsteps, otherback->numsteps, endsite*sizeof(long));
    if (other == *root) {
      memcpy(otherback->base, zeros, endsite*sizeof(long));
      memcpy(otherback->numsteps, zeros, endsite*sizeof(long));
    } else {
      memcpy(otherback->base, other->back->base, endsite*sizeof(long));
      memcpy(otherback->numsteps, other->back->numsteps, endsite*sizeof(long));
    }
    p = other->back;
    other->back = otherback;
    if (other == *root)
      preorder(other, otherback, *root, otherback, NULL, other, -1);
    else
      preorder(other, otherback, *root, NULL, NULL, NULL, 0);
    other->back = p;
    if (other != *root) {
      memcpy(other->oldbase,(*fork)->base, endsite*sizeof(long));
      memcpy(other->oldnumsteps,(*fork)->numsteps, endsite*sizeof(long));
      preorder(other->back, other, *root, NULL, NULL, NULL, 0);
    }
  } else {
    memcpy(item->oldbase, item->base, endsite*sizeof(long));
    memcpy(item->oldnumsteps, item->numsteps, endsite*sizeof(long));
    memcpy(item->base, zeros, endsite*sizeof(long));
    memcpy(item->numsteps, zeros, endsite*sizeof(long));
    preorder(*fork, item, *root, NULL, NULL, *fork, -1);
    if (*fork != *root)
      preorder((*fork)->back, *fork, *root, NULL, NULL, NULL, 0);
    memcpy(item->base, item->oldbase, endsite*sizeof(long));
    memcpy(item->numsteps, item->oldnumsteps, endsite*sizeof(long));
  }
}  /* remove */


void postorder(node *p)
{
  /* traverses an n-ary tree, suming up steps at a node's descendants */
  /* used in dnacomp, dnapars, & dnapenny */
  node *q;

  if (p->tip)
    return;
  q = p->next;
  while (q != p) {
    postorder(q->back);
    q = q->next;
  }
  zeronumnuc(p, endsite);
  if (p->numdesc > 2)
    multisumnsteps2(p);
  else
    fillin(p, p->next->back, p->next->next->back);
}  /* postorder */


void getnufork(node **nufork,node **grbg,pointarray treenode,long *zeros)
{
  /* find a fork not used currently */
  long i;

  i = spp;
  while (treenode[i] && treenode[i]->numdesc > 0) i++;
  if (!treenode[i])
    gnutreenode(grbg, &treenode[i], i, endsite, zeros);
  *nufork = treenode[i];
} /* getnufork */


void reroot(node *outgroup, node *root)
{
  /* reorients tree, putting outgroup in desired position. used if
     the root is binary. */
  /* used in dnacomp & dnapars */
  node *p, *q;

  if (outgroup->back->index == root->index)
    return;
  p = root->next;
  q = root->next->next;
  p->back->back = q->back;
  q->back->back = p->back;
  p->back = outgroup;
  q->back = outgroup->back;
  outgroup->back->back = q;
  outgroup->back = p;
}  /* reroot */


void reroot2(node *outgroup, node *root)
{
  /* reorients tree, putting outgroup in desired position. */
  /* used in dnacomp & dnapars */
  node *p;

  p = outgroup->back->next;
  while (p->next != outgroup->back)
    p = p->next;
  root->next = outgroup->back;
  p->next = root;
}  /* reroot2 */


void reroot3(node *outgroup, node *root, node *root2, node *lastdesc,
                        node **grbg)
{
  /* reorients tree, putting back outgroup in original position. */
  /* used in dnacomp & dnapars */
  node *p;

  p = root->next;
  while (p->next != root)
    p = p->next;
  chucktreenode(grbg, root);
  p->next = outgroup->back;
  root2->next = lastdesc->next;
  lastdesc->next = root2;
}  /* reroot3 */


void savetraverse(node *p)
{
  /* sets BOOLEANs that indicate which way is down */
  node *q;

  p->bottom = true;
  if (p->tip)
    return;
  q = p->next;
  while (q != p) {
    q->bottom = false;
    savetraverse(q->back);
    q = q->next;
  }
}  /* savetraverse */


void newindex(long i, node *p)
{
  /* assigns index i to node p */

  while (p->index != i) {
    p->index = i;
    p = p->next;
  }
} /* newindex */


void flipindexes(long nextnode, pointarray treenode)
{
  /* flips index of nodes between nextnode and last node.  */
  long last;
  node *temp;

  last = nonodes;
  while (treenode[last - 1]->numdesc == 0)
    last--;
  if (last > nextnode) {
    temp = treenode[nextnode - 1];
    treenode[nextnode - 1] = treenode[last - 1];
    treenode[last - 1] = temp;
    newindex(nextnode, treenode[nextnode - 1]);
    newindex(last, treenode[last - 1]);
  }
} /* flipindexes */  


boolean parentinmulti(node *anode)
{
  /* sees if anode's parent has more than 2 children */
  node *p;

  while (!anode->bottom) anode = anode->next;
  p = anode->back;
  while (!p->bottom)
    p = p->next;
  return (p->numdesc > 2);
} /* parentinmulti */


long sibsvisited(node *anode, long *place)
{
  /* computes the number of nodes which are visited earlier than anode among 
  its siblings */
  node *p;
  long nvisited;

  while (!anode->bottom) anode = anode->next;
  p = anode->back->next;
  nvisited = 0;
  do {
    if (!p->bottom && place[p->back->index - 1] != 0)
      nvisited++;
    p = p->next;
  } while (p != anode->back);
  return nvisited;
}  /* sibsvisited */


long smallest(node *anode, long *place)
{
  /* finds the smallest index of sibling of anode */
  node *p;
  long min;

  while (!anode->bottom) anode = anode->next;
  p = anode->back->next;
  if (p->bottom) p = p->next;
  min = nonodes;
  do {
    if (p->back && place[p->back->index - 1] != 0) {
      if (p->back->index <= spp) {
        if (p->back->index < min)
          min = p->back->index;
      } else {
        if (place[p->back->index - 1] < min)
          min = place[p->back->index - 1];
      }
    }
    p = p->next;
    if (p->bottom) p = p->next;
  } while (p != anode->back);
  return min;
}  /* smallest */


void bintomulti(node **root, node **binroot, node **grbg, long *zeros)
{  /* attaches root's left child to its right child and makes
      the right child new root */
  node *left, *right, *newnode, *temp;

  right = (*root)->next->next->back;
  left = (*root)->next->back;
  if (right->tip) {
    (*root)->next = right->back;
    (*root)->next->next = left->back;
    temp = left;
    left = right;
    right = temp;
    right->back->next = *root;
  }
  gnutreenode(grbg, &newnode, right->index, endsite, zeros);
  newnode->next = right->next;
  newnode->back = left;
  left->back = newnode;
  right->next = newnode;
  (*root)->next->back = (*root)->next->next->back = NULL;
  *binroot = *root;
  (*binroot)->numdesc = 0;
  *root = right;
  (*root)->numdesc++;
  (*root)->back = NULL;
} /* bintomulti */


void backtobinary(node **root, node *binroot, node **grbg)
{ /* restores binary root */
  node *p;

  binroot->next->back = (*root)->next->back;
  (*root)->next->back->back = binroot->next;
  p = (*root)->next;
  (*root)->next = p->next;
  binroot->next->next->back = *root;
  (*root)->back = binroot->next->next;
  chucktreenode(grbg, p);
  (*root)->numdesc--;
  *root = binroot;
  (*root)->numdesc = 2;
} /* backtobinary */


boolean outgrin(node *root, node *outgrnode)
{ /* checks if outgroup node is a child of root */
  node *p;

  p = root->next;
  while (p != root) {
    if (p->back == outgrnode)
      return true;
    p = p->next;
  }
  return false;
} /* outgrin */


void flipnodes(node *nodea, node *nodeb)
{ /* flip nodes */
  node *backa, *backb;

  backa = nodea->back;
  backb = nodeb->back;
  backa->back = nodeb;
  backb->back = nodea;
  nodea->back = backb;
  nodeb->back = backa;
} /* flipnodes */


void moveleft(node *root, node *outgrnode, node **flipback)
{ /* makes outgroup node to leftmost child of root */
  node *p;
  boolean done;

  p = root->next;
  done = false;
  while (p != root && !done) {
    if (p->back == outgrnode) {
      *flipback = p;
      flipnodes(root->next->back, p->back);
      done = true;
    }
    p = p->next;
  }
} /* moveleft */


void savetree(node *root,  long *place, pointarray treenode,
                        node **grbg, long *zeros)
{ /* record in place where each species has to be
     added to reconstruct this tree */
  /* used by dnacomp & dnapars */
  long i, j, nextnode, nvisited;
  node *p, *q, *r = NULL, *root2, *lastdesc, 
                *outgrnode, *binroot, *flipback;
  boolean done, newfork;

  binroot = NULL;
  lastdesc = NULL;
  root2 = NULL;
  flipback = NULL;
  outgrnode = treenode[outgrno - 1];
  if (root->numdesc == 2)
    bintomulti(&root, &binroot, grbg, zeros);
  if (outgrin(root, outgrnode)) {
    if (outgrnode != root->next->back)
      moveleft(root, outgrnode, &flipback);
  } else {
    root2 = root;
    lastdesc = root->next;
    while (lastdesc->next != root)
      lastdesc = lastdesc->next;
    lastdesc->next = root->next;
    gnutreenode(grbg, &root, outgrnode->back->index, endsite, zeros);
    root->numdesc = root2->numdesc;
    reroot2(outgrnode, root);
  }
  savetraverse(root);
  nextnode = spp + 1;
  for (i = nextnode; i <= nonodes; i++)
    if (treenode[i - 1]->numdesc == 0)
      flipindexes(i, treenode);
  for (i = 0; i < nonodes; i++)
    place[i] = 0;
  place[root->index - 1] = 1;
  for (i = 1; i <= spp; i++) {
    p = treenode[i - 1];
    while (place[p->index - 1] == 0) {
      place[p->index - 1] = i;
      while (!p->bottom)
        p = p->next;
      r = p;
      p = p->back;
    }
    if (i > 1) {
      q = treenode[i - 1]; 
      newfork = true;
      nvisited = sibsvisited(q, place);
      if (nvisited == 0) {
        if (parentinmulti(r)) {
          nvisited = sibsvisited(r, place);
          if (nvisited == 0)
            place[i - 1] = place[p->index - 1];
          else if (nvisited == 1)
            place[i - 1] = smallest(r, place);
          else {
            place[i - 1] = -smallest(r, place);
            newfork = false;
          }
        } else
          place[i - 1] = place[p->index - 1];
      } else if (nvisited == 1) {
        place[i - 1] = place[p->index - 1];
      } else {
        place[i - 1] = -smallest(q, place);
        newfork = false;
      }
      if (newfork) {
        j = place[p->index - 1];
        done = false;
        while (!done) {
          place[p->index - 1] = nextnode;
          while (!p->bottom)
            p = p->next;
          p = p->back;