getinput.c

fastDNAml is an attempt to solve the same problem as DNAML, but to do so faster and using less memo

        default:
            printf("ERROR: Auxiliary options line starts with '%c'\n", ch);
            return FALSE;
        }
      }

    if (! adef->userwgt) {
      for (i = 1; i <= rdta->sites; i++) rdta->wgt[i] = 1;
      cdta->wgtsum = rdta->sites;
      }

    if (adef->userwgt && cdta->wgtsum < 1) {
      printf("ERROR:  Missing or bad user-supplied weight data.\n");
      return FALSE;
      }

    if (adef->boot) {
      printf("Bootstrap random number seed = %ld\n\n", adef->boot);
      }

    if (adef->jumble) {
      printf("Jumble random number seed = %ld\n\n", adef->jumble);
      }

    if (adef->qadd) {
      printf("Quick add (only local branches initially optimized) in effect\n\n");
      }

    if (rdta->categs > 0) {
      printf("Site category   Rate of change\n\n");
      for (i = 1; i <= rdta->categs; i++)
        printf("           %c%13.3f\n", itobase36(i), rdta->catrat[i]);
      putchar('\n');
      for (i = 1; i <= rdta->sites; i++) {
        if ((rdta->wgt[i] > 0) && (rdta->sitecat[i] < 1)) {
          printf("ERROR: Bad category (%c) at site %d\n",
                  itobase36(rdta->sitecat[i]), i);
          return FALSE;
          }
        }
      }
    else if (rdta->categs < 0) {
      printf("ERROR: Category auxiliary data missing from input\n");
      return FALSE;
      }
    else {                                        /* rdta->categs == 0 */
      for (i = 1; i <= rdta->sites; i++) rdta->sitecat[i] = 1;
      rdta->categs = 1;
      }

    if (tr->outgr < 1) {
      printf("ERROR: Outgroup auxiliary data missing from input\n");
      return FALSE;
      }

    if (rdta->ttratio < 0.0) {
      printf("ERROR: Transition/transversion auxiliary data missing from input\n");
      return FALSE;
      }

    if (tr->global < 0) {
      if (tr->global == -2) tr->global = tr->mxtips - 3;  /* Default global */
      else                  tr->global = adef->usertree ? 0 : 1;/* No global */
      }

    if (adef->restart) {
      printf("Restart option in effect.  ");
      printf("Sequence addition will start from appended tree.\n\n");
      }

    if (adef->usertree && ! tr->global) {
      printf("User-supplied tree topology%swill be used.\n\n",
        tr->userlen ? " and branch lengths " : " ");
      }
    else {
      if (! adef->usertree) {
        printf("Rearrangements of partial trees may cross %d %s.\n",
               tr->partswap, tr->partswap == 1 ? "branch" : "branches");
        }
      printf("Rearrangements of full tree may cross %d %s.\n\n",
             tr->global, tr->global == 1 ? "branch" : "branches");
      }

    if (! adef->usertree && adef->nkeep == 0) adef->nkeep = 1;

    return TRUE;
  } /* getoptions */


boolean getbasefreqs (FILE *seqfp, rawdata *rdta)
  { /* getbasefreqs */
    int  ch;

    if (rdta->freqread) return TRUE;

    ch = getc(seqfp);
    if (! ((ch == 'F') || (ch == 'f')))  (void) ungetc(ch, seqfp);

    if (fscanf(seqfp, "%lf%lf%lf%lf",
               & rdta->freqa, & rdta->freqc,
               & rdta->freqg, & rdta->freqt) != 4 ||
        findch(seqfp,'\n') == EOF) {
      printf("ERROR: Problem reading user base frequencies\n");
      return  FALSE;
      }

    return TRUE;
  } /* getbasefreqs */


boolean getyspace (rawdata *rdta)
  { /* getyspace */
    long   size;
    int    i;
    yType *y0;

    if (! (rdta->y = (yType **) Malloc((rdta->numsp + 1) * sizeof(yType *)))) {
      printf("ERROR: Unable to obtain space for data array pointers\n");
      return  FALSE;
      }

    size = 4 * (rdta->sites / 4 + 1);
    if (! (y0 = (yType *) Malloc((rdta->numsp + 1) * size * sizeof(yType)))) {
      printf("ERROR: Unable to obtain space for data array\n");
      return  FALSE;
      }

    for (i = 0; i <= rdta->numsp; i++) {
      rdta->y[i] = y0;
      y0 += size;
      }

    return  TRUE;
  } /* getyspace */


boolean setupTree (tree *tr, int nsites)
  { /* setupTree */
    nodeptr  p0, p, q;
    int      i, j, tips, inter;

    tips  = tr->mxtips;
    inter = tr->mxtips - 1;

    if (!(p0 = (nodeptr) Malloc((tips + 3*inter) * sizeof(node)))) {
      printf("ERROR: Unable to obtain sufficient tree memory\n");
      return  FALSE;
      }

    if (!(tr->nodep = (nodeptr *) Malloc((2*tr->mxtips) * sizeof(nodeptr)))) {
      printf("ERROR: Unable to obtain sufficient tree memory, too\n");
      return  FALSE;
      }

    tr->nodep[0] = (node *) NULL;    /* Use as 1-based array */

    for (i = 1; i <= tips; i++) {    /* Set-up tips */
      p = p0++;
      p->x      = (xarray *) NULL;
      p->tip    = (yType *) NULL;
      p->number = i;
      p->next   = p;
      p->back   = (node *) NULL;

      tr->nodep[i] = p;
      }

    for (i = tips + 1; i <= tips + inter; i++) { /* Internal nodes */
      q = (node *) NULL;
      for (j = 1; j <= 3; j++) {
        p = p0++;
        p->x      = (xarray *) NULL;
        p->tip    = (yType *) NULL;
        p->number = i;
        p->next   = q;
        p->back   = (node *) NULL;
        q = p;
        }
      p->next->next->next = p;
      tr->nodep[i] = p;
      }

    tr->likelihood  = unlikely;
    tr->start       = (node *) NULL;
    tr->outgrnode   = tr->nodep[tr->outgr];
    tr->ntips       = 0;
    tr->nextnode    = 0;
    tr->opt_level   = 0;
    tr->prelabeled  = TRUE;
    tr->smoothed    = FALSE;
    tr->log_f_valid = 0;

    tr->log_f = (double *) Malloc(nsites * sizeof(double));
    if (! tr->log_f) {
      printf("ERROR: Unable to obtain sufficient tree memory, trey\n");
      return  FALSE;
      }

    return TRUE;
  } /* setupTree */


void freeTreeNode (nodeptr p)   /* Free tree node (sector) associated data */
  { /* freeTreeNode */
    if (p) {
      if (p->x) {
        if (p->x->lv) Free(p->x->lv);
        Free(p->x);
        }
      }
  } /* freeTreeNode */


void freeTree (tree *tr)
  { /* freeTree */
    nodeptr  p, q;
    int  i, tips, inter;

    tips  = tr->mxtips;
    inter = tr->mxtips - 1;

    for (i = 1; i <= tips; i++) freeTreeNode(tr->nodep[i]);

    for (i = tips + 1; i <= tips + inter; i++) {
      if (p = tr->nodep[i]) {
        if (q = p->next) {
          freeTreeNode(q->next);
          freeTreeNode(q);
          }
        freeTreeNode(p);
        }
      }

    Free(tr->nodep[1]);       /* Free the actual nodes */
  } /* freeTree */


boolean getdata (FILE *seqfp, analdef *adef, rawdata *rdta, tree *tr)
    /* read sequences */
  { /* getdata */
    int   i, j, k, l, basesread, basesnew, ch;
    int   meaning[256];          /*  meaning of input characters */ 
    char *nameptr;
    boolean  allread, firstpass;

    for (i = 0; i <= 255; i++) meaning[i] = 0;
    meaning['A'] =  1;
    meaning['B'] = 14;
    meaning['C'] =  2;
    meaning['D'] = 13;
    meaning['G'] =  4;
    meaning['H'] = 11;
    meaning['K'] = 12;
    meaning['M'] =  3;
    meaning['N'] = 15;
    meaning['O'] = 15;
    meaning['R'] =  5;
    meaning['S'] =  6;
    meaning['T'] =  8;
    meaning['U'] =  8;
    meaning['V'] =  7;
    meaning['W'] =  9;
    meaning['X'] = 15;
    meaning['Y'] = 10;
    meaning['?'] = 15;
    meaning['-'] = 15;

    basesread = basesnew = 0;

    allread = FALSE;
    firstpass = TRUE;
    ch = ' ';

    while (! allread) {
      for (i = 1; i <= tr->mxtips; i++) {     /*  Read data line */

        if (firstpass) {                      /*  Read species names */
          j = 1;
          while (whitechar(ch = getc(seqfp))) {  /*  Skip blank lines */
            if (ch == '\n')  j = 1;  else  j++;
            }

          if (j > nmlngth) {
            printf("ERROR: Blank name for species %d; ", i);
            printf("check number of species,\n");
            printf("       number of sites, and interleave option.\n");
            return  FALSE;
            }

          nameptr = tr->nodep[i]->name;
          for (k = 1; k < j; k++)  *nameptr++ = ' ';

          while (ch != '\n' && ch != EOF) {
            if (whitechar(ch))  ch = ' ';
            *nameptr++ = ch;
            if (++j > nmlngth) break;
            ch = getc(seqfp);
            }

          while (*(--nameptr) == ' ') ;          /*  remove trailing blanks */
          *(++nameptr) = '\0';                   /*  add null termination */

          if (ch == EOF) {
            printf("ERROR: End-of-file in name of species %d\n", i);
            return  FALSE;
            }
          }    /* if (firstpass) */

        j = basesread;
        while ((j < rdta->sites)
               && ((ch = getc(seqfp)) != EOF)
               && ((! adef->interleaved) || (ch != '\n'))) {
          uppercase(& ch);
          if (meaning[ch] || ch == '.') {
            j++;
            if (ch == '.') {
              if (i != 1) ch = rdta->y[1][j];
              else {
                printf("ERROR: Dot (.) found at site %d of sequence 1\n", j);
                return  FALSE;
                }
              }
            rdta->y[i][j] = ch;
            }
          else if (whitechar(ch) || digitchar(ch)) ;
          else {
            printf("ERROR: Bad base (%c) at site %d of sequence %d\n",
                    ch, j, i);
            return  FALSE;
            }
          }

        if (ch == EOF) {
          printf("ERROR: End-of-file at site %d of sequence %d\n", j, i);
          return  FALSE;
          }

        if (! firstpass && (j == basesread)) i--;        /* no data on line */
        else if (i == 1) basesnew = j;
        else if (j != basesnew) {
          printf("ERROR: Sequences out of alignment\n");
          printf("%d (instead of %d) residues read in sequence %d\n",
                  j - basesread, basesnew - basesread, i);
          return  FALSE;
          }

        while (ch != '\n' && ch != EOF) ch = getc(seqfp);  /* flush line */
        }                                                  /* next sequence */
      firstpass = FALSE;
      basesread = basesnew;
      allread = (basesread >= rdta->sites);
      }

    /*  Print listing of sequence alignment */

    if (adef->prdata) {
      j = nmlngth - 5 + ((rdta->sites + ((rdta->sites - 1)/10))/2);
      if (j < nmlngth - 1) j = nmlngth - 1;
      if (j > 37) j = 37;
      printf("Name"); for (i=1;i<=j;i++) putchar(' '); printf("Sequences\n");
      printf("----"); for (i=1;i<=j;i++) putchar(' '); printf("---------\n");
      putchar('\n');

      for (i = 1; i <= rdta->sites; i += 60) {
        l = i + 59;
        if (l > rdta->sites) l = rdta->sites;

        if (adef->userwgt) {
          printf("Weights   ");
          for (j = 11; j <= nmlngth+3; j++) putchar(' ');
          for (k = i; k <= l; k++) {
            putchar(itobase36(rdta->wgt[k]));
            if (((k % 10) == 0) && (k < l)) putchar(' ');
            }
          putchar('\n');
          }

        if (rdta->categs > 1) {
          printf("Categories");
          for (j = 11; j <= nmlngth+3; j++) putchar(' ');
          for (k = i; k <= l; k++) {
            putchar(itobase36(rdta->sitecat[k]));
            if (((k % 10) == 0) && (k < l)) putchar(' ');
            }
          putchar('\n');
          }

        for (j = 1; j <= tr->mxtips; j++) {
          nameptr = tr->nodep[j]->name;
          k = nmlngth+3;
          while (ch = *nameptr++) {putchar(ch); k--;}
          while (--k >= 0) putchar(' ');

          for (k = i; k <= l; k++) {
            ch = rdta->y[j][k];
            if ((j > 1) && (ch == rdta->y[1][k])) ch = '.';
            putchar(ch);
            if (((k % 10) == 0) && (k < l)) putchar(' ');
            }
          putchar('\n');
          }
        putchar('\n');
        }
      }

    for (j = 1; j <= tr->mxtips; j++)    /* Convert characters to meanings */
      for (i = 1; i <= rdta->sites; i++) {
        rdta->y[j][i] = meaning[rdta->y[j][i]];
        }

    return  TRUE;
  } /* getdata */


boolean  getntrees (FILE *seqfp, analdef *adef)
  { /* getntrees */

    if (fscanf(seqfp, "%d", &(adef->numutrees)) != 1 || findch(seqfp,'\n') == EOF) {
      printf("ERROR: Problem reading number of user trees\n");
      return  FALSE;
      }

    if (adef->nkeep == 0) adef->nkeep = adef->numutrees;

    return  TRUE;
  } /* getntrees */


boolean getinput (FILE *seqfp, analdef *adef, rawdata *rdta, cruncheddata *cdta,
                  tree *tr)
  { /* getinput */
    if (! getnums(seqfp,rdta))                       return FALSE;
    if (! getoptions(seqfp,adef,rdta,cdta,tr))       return FALSE;
    if (! adef->empf && ! getbasefreqs(seqfp,rdta))  return FALSE;
    if (! getyspace(rdta))                           return FALSE;
    if (! setupTree(tr, rdta->sites))                return FALSE;
    if (! getdata(seqfp,adef, rdta, tr))             return FALSE;
    if (adef->usertree && ! getntrees(seqfp,adef))   return FALSE;

    return TRUE;
  } /* getinput */