seq.c

一个神经网络工具箱

        sum = (*freqa) * treenode[i]->x[j][0][0];
        sum += (*freqc) * treenode[i]->x[j][0][(long)C - (long)A];
        sum += (*freqg) * treenode[i]->x[j][0][(long)G - (long)A];
        sum += (*freqt) * treenode[i]->x[j][0][(long)T - (long)A];
        suma += w * (*freqa) * treenode[i]->x[j][0][0] / sum;
        sumc += w * (*freqc) * treenode[i]->x[j][0][(long)C - (long)A] / sum;
        sumg += w * (*freqg) * treenode[i]->x[j][0][(long)G - (long)A] / sum;
        sumt += w * (*freqt) * treenode[i]->x[j][0][(long)T - (long)A] / sum;
      }
    }
    sum = suma + sumc + sumg + sumt;
    *freqa = suma / sum;
    *freqc = sumc / sum;
    *freqg = sumg / sum;
    *freqt = sumt / sum;
  }
  if (*freqa <= 0.0)
    *freqa = 0.000001;
  if (*freqc <= 0.0)
    *freqc = 0.000001;
  if (*freqg <= 0.0)
    *freqg = 0.000001;
  if (*freqt <= 0.0)
    *freqt = 0.000001;
}  /* empiricalfreqs */


void sitesort(long chars, steptr weight)
{
  /* Shell sort keeping sites, weights in same order */
  /* used in dnainvar, dnapars, dnacomp & dnapenny */
  long gap, i, j, jj, jg, k, itemp;
  boolean flip, tied;

  gap = chars / 2;
  while (gap > 0) {
    for (i = gap + 1; i <= chars; i++) {
      j = i - gap;
      flip = true;
      while (j > 0 && flip) {
        jj = alias[j - 1];
        jg = alias[j + gap - 1];
        tied = true;
        k = 1;
        while (k <= spp && tied) {
          flip = (y[k - 1][jj - 1] > y[k - 1][jg - 1]);
          tied = (tied && y[k - 1][jj - 1] == y[k - 1][jg - 1]);
          k++;
        }
        if (!flip)
          break;
        itemp = alias[j - 1];
        alias[j - 1] = alias[j + gap - 1];
        alias[j + gap - 1] = itemp;
        itemp = weight[j - 1];
        weight[j - 1] = weight[j + gap - 1];
        weight[j + gap - 1] = itemp;
        j -= gap;
      }
    }
    gap /= 2;
  }
}  /* sitesort */


void sitecombine(long chars)
{
  /* combine sites that have identical patterns */
  /* used in dnapars, dnapenny, & dnacomp */
  long i, j, k;
  boolean tied;

  i = 1;
  while (i < chars) {
    j = i + 1;
    tied = true;
    while (j <= chars && tied) {
      k = 1;
      while (k <= spp && tied) {
        tied = (tied &&
            y[k - 1][alias[i - 1] - 1] == y[k - 1][alias[j - 1] - 1]);
        k++;
      }
      if (tied) {
        weight[i - 1] += weight[j - 1];
        weight[j - 1] = 0;
        ally[alias[j - 1] - 1] = alias[i - 1];
      }
      j++;
    }
    i = j - 1;
  }
}  /* sitecombine */


void sitescrunch(long chars)
{
  /* move so one representative of each pattern of
     sites comes first */
  /* used in dnapars & dnacomp */
  long i, j, itemp;
  boolean done, found;

  done = false;
  i = 1;
  j = 2;
  while (!done) {
    if (ally[alias[i - 1] - 1] != alias[i - 1]) {
      if (j <= i)
        j = i + 1;
      if (j <= chars) {
        do {
          found = (ally[alias[j - 1] - 1] == alias[j - 1]);
          j++;
        } while (!(found || j > chars));
        if (found) {
          j--;
          itemp = alias[i - 1];
          alias[i - 1] = alias[j - 1];
          alias[j - 1] = itemp;
          itemp = weight[i - 1];
          weight[i - 1] = weight[j - 1];
          weight[j - 1] = itemp;
        } else
          done = true;
      } else
        done = true;
    }
    i++;
    done = (done || i >= chars);
  }
}  /* sitescrunch */


void sitesort2(long sites, steptr aliasweight)
{
  /* Shell sort keeping sites, weights in same order */
  /* used in dnaml & dnamnlk */
  long gap, i, j, jj, jg, k, itemp;
  boolean flip, tied, samewt;

  gap = sites / 2;
  while (gap > 0) {
    for (i = gap + 1; i <= sites; i++) {
      j = i - gap;
      flip = true;
      while (j > 0 && flip) {
        jj = alias[j - 1];
        jg = alias[j + gap - 1];
        samewt = ((weight[jj - 1] != 0) && (weight[jg - 1] != 0))
                 || ((weight[jj - 1] == 0) && (weight[jg - 1] == 0));
        tied = samewt && (category[jj - 1] == category[jg - 1]);
        flip = ((!samewt) && (weight[jj - 1] == 0))
               || (samewt && (category[jj - 1] > category[jg - 1]));
        k = 1;
        while (k <= spp && tied) {
          flip = (y[k - 1][jj - 1] > y[k - 1][jg - 1]);
          tied = (tied && y[k - 1][jj - 1] == y[k - 1][jg - 1]);
          k++;
        }
        if (!flip)
          break;
        itemp = alias[j - 1];
        alias[j - 1] = alias[j + gap - 1];
        alias[j + gap - 1] = itemp;
        itemp = aliasweight[j - 1];
        aliasweight[j - 1] = aliasweight[j + gap - 1];
        aliasweight[j + gap - 1] = itemp;
        j -= gap;
      }
    }
    gap /= 2;
  }
}  /* sitesort2 */


void sitecombine2(long sites, steptr aliasweight)
{
  /* combine sites that have identical patterns */
  /* used in dnaml & dnamlk */
  long i, j, k;
  boolean tied, samewt;

  i = 1;
  while (i < sites) {
    j = i + 1;
    tied = true;
    while (j <= sites && tied) {
      samewt = ((aliasweight[i - 1] != 0) && (aliasweight[j - 1] != 0))
               || ((aliasweight[i - 1] == 0) && (aliasweight[j - 1] == 0));
      tied = samewt
             && (category[alias[i - 1] - 1] == category[alias[j - 1] - 1]);
      k = 1;
      while (k <= spp && tied) {
        tied = (tied &&
            y[k - 1][alias[i - 1] - 1] == y[k - 1][alias[j - 1] - 1]);
        k++;
      }
      if (!tied)
        break;
      aliasweight[i - 1] += aliasweight[j - 1];
      aliasweight[j - 1] = 0;
      ally[alias[j - 1] - 1] = alias[i - 1];
      j++;
    }
    i = j;
  }
}  /* sitecombine2 */


void sitescrunch2(long sites, long i, long j, steptr aliasweight)
{
  /* move so positively weighted sites come first */
  /* used by dnainvar, dnaml, dnamlk, & restml */
  long itemp;
  boolean done, found;

  done = false;
  while (!done) {
    if (aliasweight[i - 1] > 0)
      i++;
    else {
      if (j <= i)
        j = i + 1;
      if (j <= sites) {
        do {
          found = (aliasweight[j - 1] > 0);
          j++;
        } while (!(found || j > sites));
        if (found) {
          j--;
          itemp = alias[i - 1];
          alias[i - 1] = alias[j - 1];
          alias[j - 1] = itemp;
          itemp = aliasweight[i - 1];
          aliasweight[i - 1] = aliasweight[j - 1];
          aliasweight[j - 1] = itemp;
        } else
          done = true;
      } else
        done = true;
    }
    done = (done || i >= sites);
  }
}  /* sitescrunch2 */


void makevalues(pointarray treenode, long *zeros, boolean usertree)
{
  /* set up fractional likelihoods at tips */
  /* used by dnacomp, dnapars, & dnapenny */
  long i, j;
  char ns = 0;
  node *p;

  setuptree(treenode, nonodes, usertree);
  for (i = 0; i < spp; i++)
    alloctip(treenode[i], zeros);
  if (!usertree) {
    for (i = spp; i < nonodes; i++) {
      p = treenode[i];
      do {
        allocnontip(p, zeros, endsite);
        p = p->next;
      } while (p != treenode[i]);
    }
  }
  for (j = 0; j < endsite; j++) {
    for (i = 0; i < spp; i++) {
      switch (y[i][alias[j] - 1]) {

      case 'A':
        ns = 1 << A;
        break;

      case 'C':
        ns = 1 << C;
        break;

      case 'G':
        ns = 1 << G;
        break;

      case 'U':
        ns = 1 << T;
        break;

      case 'T':
        ns = 1 << T;
        break;

      case 'M':
        ns = (1 << A) | (1 << C);
        break;

      case 'R':
        ns = (1 << A) | (1 << G);
        break;

      case 'W':
        ns = (1 << A) | (1 << T);
        break;

      case 'S':
        ns = (1 << C) | (1 << G);
        break;

      case 'Y':
        ns = (1 << C) | (1 << T);
        break;

      case 'K':
        ns = (1 << G) | (1 << T);
        break;

      case 'B':
        ns = (1 << C) | (1 << G) | (1 << T);
        break;

      case 'D':
        ns = (1 << A) | (1 << G) | (1 << T);
        break;

      case 'H':
        ns = (1 << A) | (1 << C) | (1 << T);
        break;

      case 'V':
        ns = (1 << A) | (1 << C) | (1 << G);
        break;

      case 'N':
        ns = (1 << A) | (1 << C) | (1 << G) | (1 << T);
        break;

      case 'X':
        ns = (1 << A) | (1 << C) | (1 << G) | (1 << T);
        break;

      case '?':
        ns = (1 << A) | (1 << C) | (1 << G) | (1 << T) | (1 << O);
        break;

      case 'O':
        ns = 1 << O;
        break;

      case '-':
        ns = 1 << O;
        break;
      }
      treenode[i]->base[j] = ns;
      treenode[i]->numsteps[j] = 0;
    }
  }
}  /* makevalues */


void makevalues2(long categs, pointarray treenode, long endsite,
                        long spp, sequence y, steptr alias)
{
  /* set up fractional likelihoods at tips */
  /* used by dnaml & dnamlk */
  long i, j, k, l;
  bases b;

  for (k = 0; k < endsite; k++) {
    j = alias[k];
    for (i = 0; i < spp; i++) {
      for (l = 0; l < categs; l++) {
        for (b = A; (long)b <= (long)T; b = (bases)((long)b + 1))
          treenode[i]->x[k][l][(long)b - (long)A] = 0.0;
        switch (y[i][j - 1]) {

        case 'A':
          treenode[i]->x[k][l][0] = 1.0;
          break;

        case 'C':
          treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
          break;

        case 'G':
          treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
          break;

        case 'T':
          treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
          break;

        case 'U':
          treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
          break;

        case 'M':
          treenode[i]->x[k][l][0] = 1.0;
          treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
          break;

        case 'R':
          treenode[i]->x[k][l][0] = 1.0;
          treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
          break;

        case 'W':
          treenode[i]->x[k][l][0] = 1.0;
          treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
          break;

        case 'S':
          treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
          treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
          break;

        case 'Y':
          treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
          treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
          break;

        case 'K':
          treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
          treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
          break;

        case 'B':
          treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
          treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
          treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
          break;

        case 'D':
          treenode[i]->x[k][l][0] = 1.0;
          treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
          treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
          break;

        case 'H':
          treenode[i]->x[k][l][0] = 1.0;
          treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
          treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
          break;

        case 'V':
          treenode[i]->x[k][l][0] = 1.0;
          treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
          treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
          break;

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

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

        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;

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