getinput.c

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

#include  <stdlib.h>
#include  <stdio.h>
#include  <string.h>
#include  "fastDNAml_types.h"
#include  "fastDNAml_funcs.h"


/*==============================================================================
 *                      READ A DNA SEQUENCE DATA FILE                          *
 *============================================================================*/


#if 0
  void hang(char *msg)
   { printf("Hanging around: %s\n", msg);
     while(1);
     }
#endif


boolean getnums (FILE *seqfp, rawdata *rdta)
    /* input number of species, number of sites */
  { /* getnums */
    printf("\n%s, version %s, %s,\nCopyright (C) 1998, 1999, 2000 by Gary J. Olsen\n\n",
            programName,
            programVersion,
            programDate);
    printf("Based in part on Joseph Felsenstein's\n\n");
    printf("   Nucleic acid sequence Maximum Likelihood method, version 3.3\n\n\n");

    if (fscanf(seqfp, "%d %d", & rdta->numsp, & rdta->sites) != 2) {
      printf("ERROR: Problem reading number of species and sites\n");
      return FALSE;
      }
    printf("%d Species, %d Sites\n\n", rdta->numsp, rdta->sites);

    if (rdta->numsp < 4) {
      printf("TOO FEW SPECIES\n");
      return FALSE;
      }

    if (rdta->sites < 1) {
      printf("TOO FEW SITES\n");
      return FALSE;
      }

    return TRUE;
  } /* getnums */


boolean inputboot(FILE *fp, analdef *adef)
    /* read the bootstrap auxilliary info */
  { /* inputboot */
    if (! adef->boot) {
      printf("ERROR: Unexpected Bootstrap auxiliary data line\n");
      return FALSE;
      }
    else if (fscanf(fp, "%ld", & adef->boot) != 1 ||
             findch(fp,'\n') == EOF) {
      printf("ERROR: Problem reading boostrap random seed value\n");
      return FALSE;
      }

    return TRUE;
  } /* inputboot */


boolean inputcategories (FILE *infile, rawdata *rdta)
    /* read the category rates and the categories for each site */
  { /* inputcategories */
    int  i, j, ch, ci;

    if (rdta->categs >= 0) {
      printf("ERROR: Unexpected Categories auxiliary data line\n");
      return FALSE;
      }
    if (fscanf(infile, "%d", & rdta->categs) != 1) {
      printf("ERROR: Problem reading number of rate categories\n");
      return FALSE;
      }
    if (rdta->categs < 1 || rdta->categs > maxcategories) {
      printf("ERROR: Bad number of categories: %d\n", rdta->categs);
      printf("Must be in range 1 - %d\n", maxcategories);
      return FALSE;
      }

    for (j = 1; j <= rdta->categs
           && fscanf(infile, "%lf", &(rdta->catrat[j])) == 1; j++) ;

    if ((j <= rdta->categs) || (findch(infile,'\n') == EOF)) {
      printf("ERROR: Problem reading rate values\n");
      return FALSE;
      }

    for (i = 1; i <= nmlngth; i++)  (void) getc(infile);

    i = 1;
    while (i <= rdta->sites) {
      ch = getc(infile);
      ci = base36(ch);
      if (ci >= 0 && ci <= rdta->categs)
        rdta->sitecat[i++] = ci;
      else if (! whitechar(ch)) {
        printf("ERROR: Bad category character (%c) at site %d\n", ch, i);
        return FALSE;
        }
      }

    if (findch(infile,'\n') == EOF) {      /* skip to end of line */
      printf("ERROR: Missing newline at end of category data\n");
      return FALSE;
      }

    return TRUE;
  } /* inputcategories */


boolean inputextra (FILE *infile, analdef *adef)
  { /* inputextra */
    if (fscanf(infile,"%d", & adef->extra) != 1 ||
        findch(infile,'\n') == EOF) {
      printf("ERROR: Problem reading extra info value\n");
      return FALSE;
      }

    return TRUE;
  } /* inputextra */


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

    rdta->freqread = TRUE;
    return TRUE;
  } /* inputfreqs */


boolean inputglobal (FILE *infile, tree *tr)
    /* input the global option information */
  { /* inputglobal */
    int  ch;

    if (tr->global != -2) {
      printf("ERROR: Unexpected Global auxiliary data line\n");
      return FALSE;
      }
    if (fscanf(infile, "%d", &(tr->global)) != 1) {
      printf("ERROR: Problem reading rearrangement region size\n");
      return FALSE;
      }
    if (tr->global < 0) {
      printf("WARNING: Global region size too small;\n");
      printf("         value reset to local\n\n");
      tr->global = 1;
      }
    else if (tr->global == 0)  tr->partswap = 0;
    else if (tr->global > tr->mxtips - 3) {
      tr->global = tr->mxtips - 3;
      }

    while ((ch = getc(infile)) != '\n') {  /* Scan for second value */
      if (! whitechar(ch)) {
        if (ch != EOF)  (void) ungetc(ch, infile);
        if (ch == EOF || fscanf(infile, "%d", &(tr->partswap)) != 1
                      || findch(infile,'\n') == EOF) {
          printf("ERROR: Problem reading insert swap region size\n");
          return FALSE;
          }
        else if (tr->partswap < 0)  tr->partswap = 1;
        else if (tr->partswap > tr->mxtips - 3) {
          tr->partswap = tr->mxtips - 3;
          }

        if (tr->partswap > tr->global)  tr->global = tr->partswap;
        break;   /*  Break while loop */
        }
      }

    return TRUE;
  } /* inputglobal */


boolean inputjumble (FILE *infile, analdef *adef)
  { /* inputjumble */
    if (! adef->jumble) {
      printf("ERROR: Unexpected Jumble auxiliary data line\n");
      return FALSE;
      }
    else if (fscanf(infile, "%ld", & adef->jumble) != 1 ||
             findch(infile,'\n') == EOF) {
      printf("ERROR: Problem reading jumble random seed value\n");
      return FALSE;
      }
    else if (adef->jumble == 0) {
      printf("WARNING: Jumble random number seed is zero\n\n");
      }

    return TRUE;
  } /* inputjumble */


boolean inputkeep (FILE *infile, analdef *adef)
  { /* inputkeep */
    if (fscanf(infile, "%d", & adef->nkeep) != 1 ||
        findch(infile,'\n') == EOF || adef->nkeep < 1) {
      printf("ERROR: Problem reading number of kept trees\n");
      return FALSE;
      }

    return TRUE;
  } /* inputkeep */


boolean inputoutgroup (FILE *infile, analdef *adef, tree *tr)
  { /* inputoutgroup */
    if (! adef->root || tr->outgr > 0) {
      printf("ERROR: Unexpected Outgroup auxiliary data line\n");
      return FALSE;
      }
    else if (fscanf(infile, "%d", &(tr->outgr)) != 1 ||
             findch(infile,'\n') == EOF) {
      printf("ERROR: Problem reading outgroup number\n");
      return FALSE;
      }
    else if ((tr->outgr < 1) || (tr->outgr > tr->mxtips)) {
      printf("ERROR: Bad outgroup: '%d'\n", tr->outgr);
      return FALSE;
      }

    return TRUE;
  } /* inputoutgroup */


boolean inputratio (FILE *infile, rawdata *rdta)
  { /* inputratio */
    if (rdta->ttratio >= 0.0) {
      printf("ERROR: Unexpected Transition/transversion auxiliary data\n");
      return FALSE;
      }
    else if (fscanf(infile,"%lf", & rdta->ttratio)!=1 ||
             findch(infile,'\n') == EOF) {
      printf("ERROR: Problem reading transition/transversion ratio\n");
      return FALSE;
      }

    return TRUE;
  } /* inputratio */


/*  Y 0 is treeNone   (no tree)
    Y 1 is treeNewick
    Y 2 is treeProlog
    Y 3 is treePHYLIP
 */

boolean inputtreeopt (FILE *infile, analdef *adef)
  { /* inputtreeopt */
    if (! adef->trout) {
      printf("ERROR: Unexpected Treefile auxiliary data\n");
      return FALSE;
      }
    else if (fscanf(infile,"%d", & adef->trout) != 1 ||
             findch(infile,'\n') == EOF) {
      printf("ERROR: Problem reading output tree-type number\n");
      return FALSE;
      }
    else if ((adef->trout < 0) || (adef->trout > treeMaxType)) {
      printf("ERROR: Bad output tree-type number: '%d'\n", adef->trout);
      return FALSE;
      }

    return TRUE;
  } /* inputtreeopt */


boolean inputweights (FILE *infile, analdef *adef, rawdata *rdta, cruncheddata *cdta)
    /* input the character weights 0, 1, 2 ... 9, A, B, ... Y, Z */
  { /* inputweights */
    int i, ch, wi;

    if (! adef->userwgt || cdta->wgtsum > 0) {
      printf("ERROR: Unexpected Weights auxiliary data\n");
      return FALSE;
      }

    for (i = 2; i <= nmlngth; i++)  (void) getc(infile);
    cdta->wgtsum = 0;
    i = 1;
    while (i <= rdta->sites) {
      ch = getc(infile);
      wi = base36(ch);
      if (wi >= 0)
        cdta->wgtsum += rdta->wgt[i++] = wi;
      else if (! whitechar(ch)) {
        printf("ERROR: Bad weight character: '%c'", ch);
        printf("       Weights in dnaml must be a digit or a letter.\n");
        return FALSE;
        }
      }

    if (findch(infile,'\n') == EOF) {      /* skip to end of line */
      printf("ERROR: Missing newline at end of weight data\n");
      return FALSE;
      }

    return TRUE;
  } /* inputweights */


boolean getoptions (FILE *seqfp, analdef *adef, rawdata *rdta,
                    cruncheddata *cdta, tree *tr)
  { /* getoptions */
    int     ch, i, extranum;

    adef->boot    =           0;  /* Don't bootstrap column weights */
    adef->empf    =        TRUE;  /* Use empirical base frequencies */
    adef->extra   =           0;  /* No extra runtime info unless requested */
    adef->interleaved =    TRUE;  /* By default, data format is interleaved */
    adef->jumble  =       FALSE;  /* Use random addition sequence */
    adef->nkeep   =           0;  /* Keep only the one best tree */
    adef->prdata  =       FALSE;  /* Don't echo data to output stream */
    adef->qadd    =        TRUE;  /* Smooth branches globally in add */
    adef->restart =       FALSE;  /* Restart from user tree */
    adef->root    =       FALSE;  /* User-defined outgroup rooting */
    adef->trout   = treeDefType;  /* Output tree file */
    adef->trprint =        TRUE;  /* Print tree to output stream */
    rdta->categs  =           0;  /* No rate categories */
    rdta->catrat[1] =       1.0;  /* Rate values */
    rdta->freqread =      FALSE;  /* User-defined frequencies not read yet */
    rdta->ttratio =         2.0;  /* Transition/transversion rate ratio */
    tr->global    =          -1;  /* Default search locale for optimum */
    tr->mxtips    = rdta->numsp;
    tr->outgr     =           1;  /* Outgroup number */
    tr->partswap  =           1;  /* Default to swap locally after insert */
    tr->userlen   =       FALSE;  /* User-supplied branch lengths */
    adef->usertree =      FALSE;  /* User-defined tree topologies */
    adef->userwgt =       FALSE;  /* User-defined position weights */
    extranum      =           0;

    while ((ch = getc(seqfp)) != '\n' && ch != EOF) {
      uppercase(& ch);
      switch (ch) {
          case '1' : adef->prdata  = ! adef->prdata; break;
          case '3' : adef->trprint = ! adef->trprint; break;
          case '4' : adef->trout   = treeDefType - adef->trout; break;
          case 'B' : adef->boot    =  1; extranum++; break;
          case 'C' : rdta->categs  = -1; extranum++; break;
          case 'E' : adef->extra   = -1; break;
          case 'F' : adef->empf    = ! adef->empf; break;
          case 'G' : tr->global    = -2; break;
          case 'I' : adef->interleaved = ! adef->interleaved; break;
          case 'J' : adef->jumble  = 1; extranum++; break;
          case 'K' : extranum++; break;
          case 'L' : tr->userlen   = TRUE; break;
          case 'O' : adef->root    = TRUE; tr->outgr = 0; extranum++; break;
          case 'Q' : adef->qadd    = FALSE; break;
          case 'R' : adef->restart = TRUE; break;
          case 'T' : rdta->ttratio = -1.0; extranum++; break;
          case 'U' : adef->usertree = TRUE; break;
          case 'W' : adef->userwgt = TRUE; cdta->wgtsum = 0; extranum++; break;
          case 'Y' : adef->trout   = treeDefType - adef->trout; break;
          case ' ' : break;
          case '\t': break;
          default  :
              printf("ERROR: Bad option character: '%c'\n", ch);
              return FALSE;
          }
      }

    if (ch == EOF) {
      printf("ERROR: End-of-file in options list\n");
      return FALSE;
      }

    if (adef->usertree && adef->restart) {
      printf("ERROR:  The restart and user-tree options conflict:\n");
      printf("        Restart adds rest of taxa to a starting tree;\n");
      printf("        User-tree does not add any taxa.\n\n");
      return FALSE;
      }

    if (adef->usertree && adef->jumble) {
      printf("WARNING:  The jumble and user-tree options conflict:\n");
      printf("          Jumble adds taxa to a tree in random order;\n");
      printf("          User-tree does not use taxa addition.\n");
      printf("          Jumble option cancelled for this run.\n\n");
      adef->jumble = FALSE;
      }

    if (tr->userlen && tr->global != -1) {
      printf("ERROR:  The global and user-lengths options conflict:\n");
      printf("        Global optimizes a starting tree;\n");
      printf("        User-lengths constrain the starting tree.\n\n");
      return FALSE;
      }

    if (tr->userlen && ! adef->usertree) {
      printf("WARNING:  User lengths required user tree option.\n");
      printf("          User-tree option set for this run.\n\n");
      adef->usertree = TRUE;
      }

    rdta->wgt      = (int *)    Malloc((rdta->sites + 1) * sizeof(int));
    rdta->wgt2     = (int *)    Malloc((rdta->sites + 1) * sizeof(int));
    rdta->sitecat  = (int *)    Malloc((rdta->sites + 1) * sizeof(int));
    cdta->alias    = (int *)    Malloc((rdta->sites + 1) * sizeof(int));
    cdta->aliaswgt = (int *)    Malloc((rdta->sites + 1) * sizeof(int));
    cdta->patcat   = (int *)    Malloc((rdta->sites + 1) * sizeof(int));
    cdta->patrat   = (double *) Malloc((rdta->sites + 1) * sizeof(double));
    cdta->wr       = (double *) Malloc((rdta->sites + 1) * sizeof(double));
    cdta->wr2      = (double *) Malloc((rdta->sites + 1) * sizeof(double));
    if ( ! rdta->wgt || ! rdta->wgt2     || ! rdta->sitecat || ! cdta->alias
                     || ! cdta->aliaswgt || ! cdta->patcat  || ! cdta->patrat
                     || ! cdta->wr       || ! cdta->wr2) {
      fprintf(stderr, "getoptions: Malloc failure\n");
      return 0;
      }

    /*  process lines with auxiliary data */

    while (extranum--) {
      ch = getc(seqfp);
      uppercase(& ch);
      switch (ch) {
        case 'B':  if (! inputboot(seqfp,adef)) return FALSE; break;
        case 'C':  if (! inputcategories(seqfp,rdta)) return FALSE; break;
        case 'E':  if (! inputextra(seqfp,adef)) return FALSE; extranum++; break;
        case 'F':  if (! inputfreqs(seqfp,rdta)) return FALSE; break;
        case 'G':  if (! inputglobal(seqfp,tr)) return FALSE; extranum++; break;
        case 'J':  if (! inputjumble(seqfp,adef)) return FALSE; break;
        case 'K':  if (! inputkeep(seqfp,adef)) return FALSE; break;
        case 'O':  if (! inputoutgroup(seqfp,adef, tr)) return FALSE; break;
        case 'T':  if (! inputratio(seqfp,rdta)) return FALSE; break;
        case 'W':  if (! inputweights(seqfp,adef, rdta, cdta)) return FALSE; break;
        case 'Y':  if (! inputtreeopt(seqfp,adef)) return FALSE; extranum++; break;