dnapars.c

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

        printweights(outfile, 0, chars, weight, "Sites");
    }
  }
  makeweights();
  makevalues(treenode, zeros, usertree);
  if (!usertree) {
    allocnode(&temp, zeros, endsite);
    allocnode(&temp1, zeros, endsite);
    allocnode(&temp2, zeros, endsite);
    allocnode(&tempsum, zeros, endsite);
    allocnode(&temprm, zeros, endsite);
    allocnode(&tempadd, zeros, endsite);
    allocnode(&tempf, zeros, endsite);
    allocnode(&tmp, zeros, endsite);
    allocnode(&tmp1, zeros, endsite);
    allocnode(&tmp2, zeros, endsite);
    allocnode(&tmp3, zeros, endsite);
    allocnode(&tmprm, zeros, endsite);
    allocnode(&tmpadd, zeros, endsite);
  }
}  /* doinput */


void initdnaparsnode(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:
    gnutreenode(grbg, p, nodei, endsite, zeros);
    treenode[nodei - 1] = *p;
    break;
  case nonbottom:
    gnutreenode(grbg, p, nodei, endsite, zeros);
    break;
  case tip:
    match_names_to_data (str, treenode, p, spp);
    break;
  case length:         /* if there is a length, read it and discard value */
    processlength(&valyew, &divisor, ch, &minusread, intree, parens);
    break;
  default:            /*cases hslength,hsnolength,treewt,unittrwt,iter,*/
    break;
  }
} /* initdnaparsnode */


void evaluate(node *r)
{
  /* determines the number of steps needed for a tree. this is
     the minimum number of steps needed to evolve sequences on
     this tree */
  long i, steps;
  long term;
  double sum;

  sum = 0.0;
  for (i = 0; i < endsite; i++) {
    steps = r->numsteps[i];
    if ((long)steps <= threshwt[i])
      term = steps;
    else
      term = threshwt[i];
    sum += (double)term;
    if (usertree && which <= maxuser)
      fsteps[which - 1][i] = term;
  }
  if (usertree && which <= maxuser) {
    nsteps[which - 1] = sum;
    if (which == 1) {
      minwhich = 1;
      minsteps = sum;
    } else if (sum < minsteps) {
      minwhich = which;
      minsteps = sum;
    }
  }
  like = -sum;
}  /* evaluate */


void tryadd(node *p, node *item, node *nufork)
{
  /* temporarily adds one fork and one tip to the tree.
     if the location where they are added yields greater
     "likelihood" than other locations tested up to that
     time, then keeps that location as there */
  long pos;
  double belowsum, parentsum;
  boolean found, collapse, changethere, trysave;

  if (!p->tip) {
    memcpy(temp->base, p->base, endsite*sizeof(long));
    memcpy(temp->numsteps, p->numsteps, endsite*sizeof(long));
    memcpy(temp->numnuc, p->numnuc, endsite*sizeof(nucarray));
    temp->numdesc = p->numdesc + 1;
    if (p->back) {
      multifillin(temp, tempadd, 1);
      sumnsteps2(tempsum, temp, p->back, 0, endsite, threshwt);
    } else {
      multisumnsteps(temp, tempadd, 0, endsite, threshwt);
      tempsum->sumsteps = temp->sumsteps;
    }
    if (tempsum->sumsteps <= -bestyet) {
      if (p->back)
        sumnsteps2(tempsum, temp, p->back, endsite+1, endsite, threshwt);
      else {
        multisumnsteps(temp, temp1, endsite+1, endsite, threshwt);
        tempsum->sumsteps = temp->sumsteps;
      }
    }
    p->sumsteps = tempsum->sumsteps;
  }
  if (p == root)
    sumnsteps2(temp, item, p, 0, endsite, threshwt);
  else {
    sumnsteps(temp1, item, p, 0, endsite);
    sumnsteps2(temp, temp1, p->back, 0, endsite, threshwt);
  }
  if (temp->sumsteps <= -bestyet) {
    if (p == root)
      sumnsteps2(temp, item, p, endsite+1, endsite, threshwt);
    else {
      sumnsteps(temp1, item, p, endsite+1, endsite);
      sumnsteps2(temp, temp1, p->back, endsite+1, endsite, threshwt);
    }
  }
  belowsum = temp->sumsteps;
  multf = false;
  like = -belowsum;
  if (!p->tip && belowsum >= p->sumsteps) {
    multf = true;
    like = -p->sumsteps;
  }
  trysave = true;
  if (!multf && p != root) {
    parentsum = treenode[p->back->index - 1]->sumsteps;
    if (belowsum >= parentsum)
      trysave = false;
  }
  if (lastrearr) {
    changethere = true;
    if (like >= bstlike2 && trysave) {
      if (like > bstlike2)
        found = false;
      else {
        addnsave(p, item, nufork, &root, &grbg, multf, treenode, place, zeros);
        pos = 0;
        findtree(&found, &pos, nextree, place, bestrees);
      }
      if (!found) {
        collapse = collapsible(item, p, temp, temp1, temp2, tempsum, temprm,
                     tmpadd, multf, root, zeros, treenode);
        if (!thorough)
          changethere = !collapse;
        if (thorough || !collapse || like > bstlike2) {
          if (like > bstlike2) {
            addnsave(p, item, nufork, &root, &grbg, multf, treenode, 
                       place, zeros);
            bestlike = bstlike2 = like;
            addbestever(&pos, &nextree, maxtrees, collapse, place, bestrees);
          } else
            addtiedtree(pos, &nextree, maxtrees, collapse, place, bestrees);
        }
      }
    }
    if (like >= bestyet) {
      if (like > bstlike2)
        bstlike2 = like;
      if (changethere && trysave) {
        bestyet = like;
        there = p;
        mulf = multf;
      }
    }
  } else if ((like > bestyet) || (like >= bestyet && trysave)) {
    bestyet = like;
    there = p;
    mulf = multf;
  }
}  /* tryadd */


void addpreorder(node *p, node *item, node *nufork)
{
  /* traverses a n-ary tree, calling function tryadd
     at a node before calling tryadd at its descendants */
  node *q;

  if (p == NULL)
    return;
  tryadd(p, item, nufork);
  if (!p->tip) {
    q = p->next;
    while (q != p) {
      addpreorder(q->back, item, nufork);
      q = q->next;
    }
  }
}  /* addpreorder */


void trydescendants(node *item, node *forknode, node *parent,
                        node *parentback, boolean trybelow)
{
  /* tries rearrangements at parent and below parent's descendants */
  node *q, *tempblw;
  boolean bestever=0, belowbetter, multf=0, saved, trysave;
  double parentsum=0, belowsum;

  memcpy(temp->base, parent->base, endsite*sizeof(long));
  memcpy(temp->numsteps, parent->numsteps, endsite*sizeof(long));
  memcpy(temp->numnuc, parent->numnuc, endsite*sizeof(nucarray));
  temp->numdesc = parent->numdesc + 1;
  multifillin(temp, tempadd, 1);
  sumnsteps2(tempsum, parentback, temp, 0, endsite, threshwt);
  belowbetter = true;
  if (lastrearr) {
    parentsum = tempsum->sumsteps;
    if (-tempsum->sumsteps >= bstlike2) {
      belowbetter = false;
      bestever = false;
      multf = true;
      if (-tempsum->sumsteps > bstlike2)
        bestever = true;
      savelocrearr(item, forknode, parent, tmp, tmp1, tmp2, tmp3, tmprm,
                    tmpadd, &root, maxtrees, &nextree, multf, bestever,
                    &saved, place, bestrees, treenode, &grbg, zeros);
      if (saved) {
        like = bstlike2 = -tempsum->sumsteps;
        there = parent;
        mulf = true;
      }
    }
  } else if (-tempsum->sumsteps >= like) {
    there = parent;
    mulf = true;
    like = -tempsum->sumsteps;
  }
  if (trybelow) {
    sumnsteps(temp, parent, tempadd, 0, endsite);
    sumnsteps2(tempsum, temp, parentback, 0, endsite, threshwt);
    if (lastrearr) {
      belowsum = tempsum->sumsteps;
      if (-tempsum->sumsteps >= bstlike2 && belowbetter && 
            (forknode->numdesc > 2 ||
               (forknode->numdesc == 2 && 
                 parent->back->index != forknode->index))) {
        trysave = false;
        memcpy(temp->base, parentback->base, endsite*sizeof(long));
        memcpy(temp->numsteps, parentback->numsteps, endsite*sizeof(long));
        memcpy(temp->numnuc, parentback->numnuc, endsite*sizeof(nucarray));
        temp->numdesc = parentback->numdesc + 1;
        multifillin(temp, tempadd, 1);
        sumnsteps2(tempsum, parent, temp, 0, endsite, threshwt);
        if (-tempsum->sumsteps < bstlike2) {
          multf = false;
          bestever = false;
          trysave = true;
        }
        if (-belowsum > bstlike2) {
          multf = false;
          bestever = true;
          trysave = true;
        }
        if (trysave) {
          if (treenode[parent->index - 1] != parent)
            tempblw = parent->back;
          else
            tempblw = parent;
          savelocrearr(item, forknode, tempblw, tmp, tmp1, tmp2, tmp3, tmprm,
                         tmpadd, &root, maxtrees, &nextree, multf, bestever,
                         &saved, place, bestrees, treenode, &grbg, zeros);
          if (saved) {
            like = bstlike2 = -belowsum;
            there = tempblw;
            mulf = false;
          }
        }
      }
    } else if (-tempsum->sumsteps > like) {
      like = -tempsum->sumsteps;
      if (-tempsum->sumsteps > bestyet) {
        if (treenode[parent->index - 1] != parent)
          tempblw = parent->back;
        else
          tempblw = parent;
        there = tempblw;
        mulf = false;
      }
    }
  }
  q = parent->next;
  while (q != parent) {
    if (q->back && q->back != item) {
      memcpy(temp1->base, q->base, endsite*sizeof(long));
      memcpy(temp1->numsteps, q->numsteps, endsite*sizeof(long));
      memcpy(temp1->numnuc, q->numnuc, endsite*sizeof(nucarray));
      temp1->numdesc = q->numdesc;
      multifillin(temp1, parentback, 0);
      if (lastrearr)
        belowbetter = (-parentsum < bstlike2);
      if (!q->back->tip) {
        memcpy(temp->base, q->back->base, endsite*sizeof(long));
        memcpy(temp->numsteps, q->back->numsteps, endsite*sizeof(long));
        memcpy(temp->numnuc, q->back->numnuc, endsite*sizeof(nucarray));
        temp->numdesc = q->back->numdesc + 1;
        multifillin(temp, tempadd, 1);
        sumnsteps2(tempsum, temp1, temp, 0, endsite, threshwt);
        if (lastrearr) {
          if (-tempsum->sumsteps >= bstlike2) {
            belowbetter = false;
            bestever = false;
            multf = true;
            if (-tempsum->sumsteps > bstlike2)
              bestever = true;
            savelocrearr(item, forknode, q->back, tmp, tmp1, tmp2, tmp3, tmprm,
                          tmpadd, &root, maxtrees, &nextree, multf, bestever,
                          &saved, place, bestrees, treenode, &grbg, zeros);
            if (saved) {
              like = bstlike2 = -tempsum->sumsteps;
              there = q->back;
              mulf = true;
            }
          }
        } else if (-tempsum->sumsteps >= like) {
          like = -tempsum->sumsteps;
          there = q->back;
          mulf = true;
        }
      }
      sumnsteps(temp, q->back, tempadd, 0, endsite);
      sumnsteps2(tempsum, temp, temp1, 0, endsite, threshwt);
      if (lastrearr) {
        if (-tempsum->sumsteps >= bstlike2) {
          trysave = false;
          multf = false;
          if (belowbetter) {
            bestever = false;
            trysave = true;
          }
          if (-tempsum->sumsteps > bstlike2) {
            bestever = true;
            trysave = true;
          }
          if (trysave) {
            if (treenode[q->back->index - 1] != q->back)
              tempblw = q;
            else
              tempblw = q->back;
            savelocrearr(item, forknode, tempblw, tmp, tmp1, tmp2, tmp3, tmprm,
                        tmpadd, &root, maxtrees, &nextree, multf, bestever,
                        &saved, place, bestrees, treenode, &grbg, zeros);
            if (saved) {
              like = bstlike2 = -tempsum->sumsteps;
              there = tempblw;
              mulf = false;
            }
          }
        }
      } else if (-tempsum->sumsteps > like) {
        like = -tempsum->sumsteps;
        if (-tempsum->sumsteps > bestyet) {
          if (treenode[q->back->index - 1] != q->back)
            tempblw = q;
          else
            tempblw = q->back;
          there = tempblw;
          mulf = false;
        }
      }
    }
    q = q->next;
  }
} /* trydescendants */


void trylocal(node *item, node *forknode)
{
  /* rearranges below forknode, below descendants of forknode when
     there are more than 2 descendants, then unroots the back of
     forknode and rearranges on its descendants */
  node *q;
  boolean bestever, multf, saved;

  memcpy(temprm->base, zeros, endsite*sizeof(long));
  memcpy(temprm->numsteps, zeros, endsite*sizeof(long));
  memcpy(temprm->oldbase, item->base, endsite*sizeof(long));
  memcpy(temprm->oldnumsteps, item->numsteps, endsite*sizeof(long));
  memcpy(tempf->base, forknode->base, endsite*sizeof(long));
  memcpy(tempf->numsteps, forknode->numsteps, endsite*sizeof(long));
  memcpy(tempf->numnuc, forknode->numnuc, endsite*sizeof(nucarray));
  tempf->numdesc = forknode->numdesc - 1;
  multifillin(tempf, temprm, -1);
  if (!forknode->back) {
    sumnsteps2(tempsum, tempf, tempadd, 0, endsite, threshwt);
    if (lastrearr) {
      if (-tempsum->sumsteps > bstlike2) {
        bestever = true;
        multf = false;
        savelocrearr(item, forknode, forknode, tmp, tmp1, tmp2, tmp3, tmprm,
                       tmpadd, &root, maxtrees, &nextree, multf, bestever,
                       &saved, place, bestrees, treenode, &grbg, zeros);
        if (saved) {
          like = bstlike2 = -tempsum->sumsteps;
          there = forknode;
          mulf = false;
        }
      }
    } else if (-tempsum->sumsteps > like) {
      like = -tempsum->sumsteps;
      if (-tempsum->sumsteps > bestyet) {
        there = forknode;
        mulf = false;
      }
    }
  } else {
    sumnsteps(temp, tempf, tempadd, 0, endsite);
    sumnsteps2(tempsum, temp, forknode->back, 0, endsite, threshwt);
    if (lastrearr) {
      if (-tempsum->sumsteps > bstlike2) {
        bestever = true;
        multf = false;
        savelocrearr(item, forknode, forknode, tmp, tmp1, tmp2, tmp3, tmprm,
                       tmpadd, &root, maxtrees, &nextree, multf, bestever,
                       &saved, place, bestrees, treenode, &grbg, zeros);
        if (saved) {