Logo Search packages:      
Sourcecode: hmmer version File versions  Download package

viterbi_exercise.c

/* viterbi_exercise.c
 * SRE, Mon Mar  9 07:55:47 1998 [St. Louis]
 * 
 * Exercise the various Viterbi algorithms, big and small.
 * 
 * RCS $Id: viterbi_exercise.c,v 1.1 1998/03/11 17:18:26 eddy Exp $
 */


#include <stdio.h>
#include <time.h>
#include <math.h>
#include <string.h>

#include "structs.h"
#include "funcs.h"
#include "globals.h"
#include "squid.h"

#ifdef MEMDEBUG
#include "dbmalloc.h"
#endif

static char banner[] = "\
viterbi_exercise : testing of Plan7 Viterbi code";

static char usage[] = "\
Usage: testdriver [-options]\n\
  Available options are:\n\
  -h              : help; display this usage info\n\
  -v              : be verbose\n\
";

static char experts[] = "\
  --hmm <f>       : use HMM in file <f>\n\
\n";

static struct opt_s OPTIONS[] = {
  { "-h",       TRUE,  sqdARG_NONE },
  { "-v",       TRUE,  sqdARG_NONE },
  { "--hmm",    FALSE, sqdARG_STRING },
};
#define NOPTIONS (sizeof(OPTIONS) / sizeof(struct opt_s))

int
main(int argc, char **argv)
{
  char    *hmmfile;             /* file to read HMM(s) from                */
  HMMFILE *hmmfp;               /* opened hmmfile for reading              */
  struct plan7_s  *hmm;         /* the HMM to search with                  */ 
  char    *dsq;               /* digitized target sequence               */
  char    *seq;
  SQINFO   sqinfo;
  int      L;                 /* length of dsq                           */
  struct p7trace_s  *tr1;     /* traceback                               */
  struct p7trace_s  *tr2;     /* another traceback                       */
  int       nseq;
  float     sc1, sc2;         /* scores                                  */
  int       config;
  int       i;

  int be_verbose;

  char *optname;                /* name of option found by Getopt()         */
  char *optarg;                 /* argument found by Getopt()               */
  int   optind;                 /* index in argv[]                          */

#ifdef MEMDEBUG
  unsigned long histid1, histid2, orig_size, current_size;
  orig_size = malloc_inuse(&histid1);
  fprintf(stderr, "[... memory debugging is ON ...]\n");
#endif
  
  /*********************************************** 
   * Parse command line
   ***********************************************/

  be_verbose = FALSE;
  hmmfile    = "fn3.hmm";
  nseq       = 100;

  while (Getopt(argc, argv, OPTIONS, NOPTIONS, usage,
                &optind, &optname, &optarg))  {
    if      (strcmp(optname, "-v")       == 0) be_verbose = TRUE;
    else if (strcmp(optname, "--hmm")    == 0) hmmfile    = optarg;
    else if (strcmp(optname, "-h")       == 0) {
      Banner(stdout, banner);
      puts(usage);
      puts(experts);
      exit(0);
    }
  }
  if (argc - optind != 0)
    Die("Incorrect number of arguments.\n%s\n", usage);

  /*********************************************** 
   * Open HMM file 
   * Read a single HMM from it. 
   ***********************************************/

  if ((hmmfp = HMMFileOpen(hmmfile, NULL)) == NULL)
    Die("Failed to open HMM file %s\n%s", hmmfile, usage);
  if (!HMMFileRead(hmmfp, &hmm)) 
    Die("Failed to read any HMMs from %s\n", hmmfile);
  if (hmm == NULL) 
    Die("HMM file %s corrupt or in incorrect format? Parse failed", hmmfile);
  Plan7Renormalize(hmm);

  /*********************************************** 
   * We cycle through different model configurations.
   * For each configuration, we repeat 100 times:
   *    - generate a sequence
   *    - score it by Viterbi and by SmallViterbi  
   *    - make sure they give OK and identical results 
   ***********************************************/

  for (config = 1; config <= 5; config++)
    {
      switch (config) {
      case 1: Plan7NakedConfig(hmm);            break;
      case 2: Plan7GlobalConfig(hmm);           break;
      case 3: Plan7LSConfig(hmm);               break;
      case 4: Plan7FSConfig(hmm, 0.5, 0.5);     break;
      case 5: Plan7SWConfig(hmm, 0.5, 0.5);     break;
      default: Die("never happens");
      }
      P7Logoddsify(hmm, TRUE);

      for (i = 0; i < nseq; i++)
      {
        EmitSequence(hmm, &dsq, &L, NULL);
        sprintf(sqinfo.name, "seq%d", i+1);
        sqinfo.len   = L;
        sqinfo.flags = SQINFO_NAME | SQINFO_LEN;

        sc1 = P7Viterbi(dsq, L, hmm, &tr1);
        sc2 = P7SmallViterbi(dsq, L, hmm, &tr2);

        if (be_verbose)
          {
            printf("Viterbi score: %.1f   SmallViterbi: %.1f\n", sc1, sc2);
            P7PrintTrace(stdout, tr1, hmm, dsq);
            P7PrintTrace(stdout, tr2, hmm, dsq);
            
            seq = DedigitizeSequence(dsq, L);
            WriteSeq(stdout, kPearson, seq, &sqinfo);
            free(seq);
          }

        if (sc1 != sc2) 
          Die("Different scores from normal/small Viterbi");

        if (fabs(sc1 - P7TraceScore(hmm, dsq, tr1)) > 0.1)
          Die("P7Viterbi score doesn't match its TraceScore");
        if (fabs(sc2 - P7TraceScore(hmm, dsq, tr2)) > 0.1)
          Die("P7SmallViterbi score doesn't match its TraceScore");

        if (! TraceVerify(tr1, hmm->M, L))
          Die("TraceVerify() failed for a P7Viterbi trace");
        if (! TraceVerify(tr2, hmm->M, L))
          Die("TraceVerify() failed for a P7SmallViterbi trace");

        if (tr1->tlen != tr2->tlen)
          Die("Trace lengths differ for normal/small Viterbi");
        if (! TraceCompare(tr1, tr2))
          Die("Different traces from normal/small Viterbi");

        P7FreeTrace(tr1);
        P7FreeTrace(tr2);
        free(dsq);
      }
    }

#ifdef MEMDEBUG
  current_size = malloc_inuse(&histid2);
  if (current_size != orig_size) Die("viterbi_exercise failed memory test");
  else fprintf(stderr, "[No memory leaks.]\n");
#endif
  
  return EXIT_SUCCESS;
}

Generated by  Doxygen 1.6.0   Back to index