protpars.c

一个神经网络工具箱

      case '6':
        trout = !trout;
        break;
      }
    } else
        printf("Not a possible option!\n");
    countup(&loopcount, 100);
  }
  */

}  /* getoptions */


void protalloctree()
{ /* allocate treenode dynamically */
  long i, j;
  node *p, *q;

  treenode = (pointarray)Malloc(nonodes*sizeof(node *));
  for (i = 0; i < (spp); i++) {
    treenode[i] = (node *)Malloc(sizeof(node));
    treenode[i]->numsteps = (steptr)Malloc(chars*sizeof(long));
    treenode[i]->siteset = (seqptr)Malloc(chars*sizeof(sitearray));
    treenode[i]->seq = (aas *)Malloc(chars*sizeof(aas));
  }
  for (i = spp; i < (nonodes); i++) {
    q = NULL;
    for (j = 1; j <= 3; j++) {
      p = (node *)Malloc(sizeof(node));
      p->numsteps = (steptr)Malloc(chars*sizeof(long));
      p->siteset = (seqptr)Malloc(chars*sizeof(sitearray));
      p->seq = (aas *)Malloc(chars*sizeof(aas));
      p->next = q;
      q = p;
    }
    p->next->next->next = p;
    treenode[i] = p;
  }
}  /* protalloctree */


void reallocnode(node* p) 
{
  free(p->numsteps);
  free(p->siteset);
  free(p->seq);
  p->numsteps = (steptr)Malloc(chars*sizeof(long));
  p->siteset = (seqptr)Malloc(chars*sizeof(sitearray));
  p->seq = (aas *)Malloc(chars*sizeof(aas));
}


void reallocchars(void) 
{ /* reallocates variables that are dependand on the number of chars
   * do we need to reallocate the garbage list too? */
  long i;
  node *p;

  if (usertree)
    for (i = 0; i < maxuser; i++) {
      free(fsteps[i]);
      fsteps[i] = (long *)Malloc(chars*sizeof(long));
    }
  
  for (i = 0; i < nonodes; i++) {
    reallocnode(treenode[i]);  
    if (i >= spp) {
      p=treenode[i]->next;
      while (p != treenode[i])  {
        reallocnode(p); 
        p = p->next;
      }
    }
  }

  free(weight);
  free(threshwt);
  free(temp->numsteps);
  free(temp->siteset);
  free(temp->seq);
  free(temp1->numsteps);
  free(temp1->siteset); 
  free(temp1->seq);

  weight = (steptr)Malloc(chars*sizeof(long));
  threshwt = (steptr)Malloc(chars*sizeof(long));
  temp->numsteps = (steptr)Malloc(chars*sizeof(long));
  temp->siteset = (seqptr)Malloc(chars*sizeof(sitearray));
  temp->seq = (aas *)Malloc(chars*sizeof(aas));
  temp1->numsteps = (steptr)Malloc(chars*sizeof(long));
  temp1->siteset = (seqptr)Malloc(chars*sizeof(sitearray));
  temp1->seq = (aas *)Malloc(chars*sizeof(aas));
}


void allocrest()
{ /* allocate remaining global arrays and variables dynamically */
  long i;

  if (usertree) {
    fsteps = (long **)Malloc(maxuser*sizeof(long *));
    for (i = 0; i < maxuser; i++)
      fsteps[i] = (long *)Malloc(chars*sizeof(long));
  }
  bestrees = (bestelm *)Malloc(maxtrees*sizeof(bestelm));
  for (i = 1; i <= maxtrees; i++)
    bestrees[i - 1].btree = (long *)Malloc(spp*sizeof(long));
  nayme = (naym *)Malloc(spp*sizeof(naym));
  enterorder = (long *)Malloc(spp*sizeof(long));
  place = (long *)Malloc(nonodes*sizeof(long));
  weight = (steptr)Malloc(chars*sizeof(long));
  threshwt = (steptr)Malloc(chars*sizeof(long));
  temp = (node *)Malloc(sizeof(node));
  temp->numsteps = (steptr)Malloc(chars*sizeof(long));
  temp->siteset = (seqptr)Malloc(chars*sizeof(sitearray));
  temp->seq = (aas *)Malloc(chars*sizeof(aas));
  temp1 = (node *)Malloc(sizeof(node));
  temp1->numsteps = (steptr)Malloc(chars*sizeof(long));
  temp1->siteset = (seqptr)Malloc(chars*sizeof(sitearray));
  temp1->seq = (aas *)Malloc(chars*sizeof(aas));
}  /* allocrest */


void doinit()
{
  /* initializes variables */

  inputnumbers(&spp, &chars, &nonodes, 1);
  getoptions();
  if (printdata)
    fprintf(outfile, "%2ld species, %3ld  sites\n\n", spp, chars);
  protalloctree();
  allocrest();
}  /* doinit*/

void protinputdata()
{
  /* input the names and sequences for each species */
  long i, j, k, l, aasread, aasnew = 0;
  Char charstate;
  boolean allread, done;
  aas aa;   /* temporary amino acid for input */

  if (printdata)
    headings(chars, "Sequences", "---------");
  aasread = 0;
  allread = false;
  while (!(allread)) {
    allread = true;
    if (eoln(infile)) {
      fscanf(infile, "%*[^\n]");
      gettc(infile);
    }
    i = 1;
    while (i <= spp) {
      if ((interleaved && aasread == 0) || !interleaved)
        initname(i - 1);
      j = interleaved ? aasread : 0;
      done = false;
      while (!done && !eoff(infile)) {
        if (interleaved)
          done = true;
        while (j < chars && !(eoln(infile) || eoff(infile))) {
          charstate = gettc(infile);
          if (charstate == ' ' || (charstate >= '0' && charstate <= '9'))
            continue;
          uppercase(&charstate);
          if ((!isalpha(charstate) && charstate != '?' &&
               charstate != '-' && charstate != '*') || charstate == 'J' ||
              charstate == 'O' || charstate == 'U') {
            printf("WARNING -- BAD AMINO ACID:%c",charstate);
            printf(" AT POSITION%5ld OF SPECIES %3ld\n",j,i);
            exxit(-1);
          }
          j++;
          aa = (charstate == 'A') ?  ala :
               (charstate == 'B') ?  asx :
               (charstate == 'C') ?  cys :
               (charstate == 'D') ?  asp :
               (charstate == 'E') ?  glu :
               (charstate == 'F') ?  phe :
               (charstate == 'G') ?  gly : aa;
          aa = (charstate == 'H') ?  his :
               (charstate == 'I') ? ileu :
               (charstate == 'K') ?  lys :
               (charstate == 'L') ?  leu :
               (charstate == 'M') ?  met :
               (charstate == 'N') ?  asn :
               (charstate == 'P') ?  pro :
               (charstate == 'Q') ?  gln :
               (charstate == 'R') ?  arg : aa;
          aa = (charstate == 'S') ?  ser :
               (charstate == 'T') ?  thr :
               (charstate == 'V') ?  val :
               (charstate == 'W') ?  trp :
               (charstate == 'X') ?  unk :
               (charstate == 'Y') ?  tyr :
               (charstate == 'Z') ?  glx :
               (charstate == '*') ? stop :
               (charstate == '?') ? quest:
               (charstate == '-') ? del  :  aa;

          treenode[i - 1]->seq[j - 1] = aa;
          memcpy(treenode[i - 1]->siteset[j - 1],
                 translate[(long)aa - (long)ala], sizeof(sitearray));
        }
        if (interleaved)
          continue;
        if (j < chars) 
          scan_eoln(infile);
        else if (j == chars)
          done = true;
      }
      if (interleaved && i == 1)
        aasnew = j;
      scan_eoln(infile);
      if ((interleaved && j != aasnew) || ((!interleaved) && j != chars)){
        printf("ERROR: SEQUENCES OUT OF ALIGNMENT\n");
        exxit(-1);}
      i++;
    }
    if (interleaved) {
      aasread = aasnew;
      allread = (aasread == chars);
    } else
      allread = (i > spp);
  }
  if (printdata) {
    for (i = 1; i <= ((chars - 1) / 60 + 1); i++) {
      for (j = 1; j <= (spp); j++) {
        for (k = 0; k < nmlngth; k++)
          putc(nayme[j - 1][k], outfile);
        fprintf(outfile, "   ");
        l = i * 60;
        if (l > chars)
          l = chars;
        for (k = (i - 1) * 60 + 1; k <= l; k++) {
          if (j > 1 && treenode[j - 1]->seq[k - 1] == treenode[0]->seq[k - 1])
            charstate = '.';
          else {
            tmpa = treenode[j-1]->seq[k-1];
              charstate =  (tmpa == ala) ? 'A' :
                           (tmpa == asx) ? 'B' :
                           (tmpa == cys) ? 'C' :
                           (tmpa == asp) ? 'D' :
                           (tmpa == glu) ? 'E' :
                           (tmpa == phe) ? 'F' :
                           (tmpa == gly) ? 'G' :
                           (tmpa == his) ? 'H' :
                           (tmpa ==ileu) ? 'I' :
                           (tmpa == lys) ? 'K' :
                           (tmpa == leu) ? 'L' : charstate;
              charstate =  (tmpa == met) ? 'M' :
                           (tmpa == asn) ? 'N' :
                           (tmpa == pro) ? 'P' :
                           (tmpa == gln) ? 'Q' :
                           (tmpa == arg) ? 'R' :
                           (tmpa == ser) ? 'S' :
                           (tmpa ==ser1) ? 'S' :
                           (tmpa ==ser2) ? 'S' : charstate;
              charstate =  (tmpa == thr) ? 'T' :
                           (tmpa == val) ? 'V' :
                           (tmpa == trp) ? 'W' :
                           (tmpa == unk) ? 'X' :
                           (tmpa == tyr) ? 'Y' :
                           (tmpa == glx) ? 'Z' :
                           (tmpa == del) ? '-' :
                           (tmpa ==stop) ? '*' :
                           (tmpa==quest) ? '?' : charstate;
        }
          putc(charstate, outfile);
          if (k % 10 == 0 && k % 60 != 0)
            putc(' ', outfile);
        }
        putc('\n', outfile);
      }
      putc('\n', outfile);
    }
    putc('\n', outfile);
  }
  putc('\n', outfile);
}  /* protinputdata */


void protmakevalues()
{
  /* set up fractional likelihoods at tips */
  long i, j;
  node *p;

  for (i = 1; i <= nonodes; i++) {
    treenode[i - 1]->back = NULL;
    treenode[i - 1]->tip = (i <= spp);
    treenode[i - 1]->index = i;
    for (j = 0; j < (chars); j++)
      treenode[i - 1]->numsteps[j] = 0;
    if (i > spp) {
      p = treenode[i - 1]->next;
      while (p != treenode[i - 1]) {
        p->back = NULL;
        p->tip = false;
        p->index = i;
        for (j = 0; j < (chars); j++)
          p->numsteps[j] = 0;
        p = p->next;
      }
    }
  }
}  /* protmakevalues */


void doinput()
{
  /* reads the input data */
  long i;

  if (justwts) {
    if (firstset)
      protinputdata();
    for (i = 0; i < chars; i++)
      weight[i] = 1;
    inputweights(chars, weight, &weights);
    if (justwts) {
      fprintf(outfile, "\n\nWeights set # %ld:\n\n", ith);
      if (progress)
        printf("\nWeights set # %ld:\n\n", ith);
    }
    if (printdata)
      printweights(outfile, 0, chars, weight, "Sites");
  } else {
    if (!firstset){
      samenumsp(&chars, ith);
      reallocchars();
    }
    for (i = 0; i < chars; i++)
      weight[i] = 1;
    if (weights) {
      inputweights(chars, weight, &weights);
    }
    if (weights)
      printweights(outfile, 0, chars, weight, "Sites");
    protinputdata();
  }
  if(!thresh)
    threshold = spp * 3.0;
  for(i = 0 ; i < (chars) ; i++){
    weight[i]*=10;
    threshwt[i] = (long)(threshold * weight[i] + 0.5);      
  }

  protmakevalues();
}  /* doinput */


void protfillin(node *p, node *left, node *rt)
{
  /* sets up for each node in the tree the aa set for site m
     at that point and counts the changes.  The program
     spends much of its time in this function */
  boolean counted;
  aas aa;
  long s = 0;
  sitearray ls, rs, qs;
  long i, m, n;

  for (m = 0; m < chars; m++) {
    if (left != NULL)
      memcpy(ls, left->siteset[m], sizeof(sitearray));
    if (rt != NULL)
      memcpy(rs, rt->siteset[m], sizeof(sitearray));
    if (left == NULL) {
      n = rt->numsteps[m];
      memcpy(qs, rs, sizeof(sitearray));
    }
    else if (rt == NULL) {
      n = left->numsteps[m];
      memcpy(qs, ls, sizeof(sitearray));
    }
    else {
      n = left->numsteps[m] + rt->numsteps[m];
      if (ls[0] == rs[0]) {
        qs[0] = ls[0];
        qs[1] = ls[1];
        qs[2] = ls[2];
      }
      else {
        counted = false;
        for (i = 0; (!counted) && (i <= 3); i++) {
          switch (i) {
 
            case 0:
              s = ls[0] & rs[0];
              break;

            case 1:
              s = (ls[0] & rs[1]) | (ls[1] & rs[0]);
              break;

            case 2:
              s = (ls[0] & rs[2]) | (ls[1] & rs[1]) | (ls[2] & rs[0]);
              break;

            case 3:
              s = ls[0] | (ls[1] & rs[2]) | (ls[2] & rs[1]) | rs[0];
              break;

          }
          if (s != 0) {
            qs[0] = s;
            counted = true;
          } else
              n += weight[m];
        }
        qs[1] = 0;
        qs[2] = 0;
        for (i = 0; i <= 1; i++) {
          for (aa = ala; (long)aa <= (long)stop; aa = (aas)((long)aa + 1)) {
            if (((1L << ((long)aa)) & qs[i]) != 0)
              qs[i+1] |= translate[(long)aa - (long)ala][1];
          }
        }
      }
    }
    p->numsteps[m] = n;
    memcpy(p->siteset[m], qs, sizeof(sitearray));
  }
}  /* protfillin */


void protpreorder(node *p)
{
  /* recompute number of steps in preorder taking both ancestoral and
     descendent steps into account */
  if (p != NULL && !p->tip) {
    protfillin (p->next, p->next->next->back, p->back);
    protfillin (p->next->next, p->back, p->next->back);
    protpreorder (p->next->back);
    protpreorder (p->next->next->back);
  }
} /* protpreorder */


void protadd(node *below, node *newtip, node *newfork)
{
  /* inserts the nodes newfork and its left descendant, newtip,
     to the tree.  below becomes newfork's right descendant */

  if (below != treenode[below->index - 1])
    below = treenode[below->index - 1];
  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;
  root->back = NULL;

  if (recompute) {
    protfillin (newfork, newfork->next->back, newfork->next->next->back);
    protpreorder(newfork);
    if (newfork != root)
      protpreorder(newfork->back);
  }
}  /* protadd */


void protre_move(node **item, node **fork)
{
  /* 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 */
  node *p, *q, *other;

  if ((*item)->back == NULL) {
    *fork = NULL;
    return;
  }
  *fork = treenode[(*item)->back->index - 1];
  if ((*item) == (*fork)->next->back)
    other = (*fork)->next->next->back;
  else other = (*fork)->next->back;
  if (root == *fork)
    root = other;
  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;