forester/archive/RIO/others/hmmer/testsuite/viterbi_exercise.c

   1 /* viterbi_exercise.c
   2  * SRE, Mon Mar  9 07:55:47 1998 [St. Louis]
   3  *
   4  * Exercise the various Viterbi algorithms, big and small.
   5  *
   6  * RCS $Id: viterbi_exercise.c,v 1.1.1.1 2005/03/22 08:34:50 cmzmasek Exp $
   7  */
   8
   9
  10 #include <stdio.h>
  11 #include <time.h>
  12 #include <math.h>
  13 #include <string.h>
  14
  15 #include "structs.h"
  16 #include "funcs.h"
  17 #include "globals.h"
  18 #include "squid.h"
  19
  20 static char banner[] = "\
  21 viterbi_exercise : testing of Plan7 Viterbi code";
  22
  23 static char usage[] = "\
  24 Usage: testdriver [-options]\n\
  25   Available options are:\n\
  26   -h              : help; display this usage info\n\
  27   -v              : be verbose\n\
  28 ";
  29
  30 static char experts[] = "\
  31   --hmm <f>       : use HMM in file <f>\n\
  32 \n";
  33
  34 static struct opt_s OPTIONS[] = {
  35   { "-h",       TRUE,  sqdARG_NONE },
  36   { "-v",       TRUE,  sqdARG_NONE },
  37   { "--hmm",    FALSE, sqdARG_STRING },
  38 };
  39 #define NOPTIONS (sizeof(OPTIONS) / sizeof(struct opt_s))
  40
  41 int
  42 main(int argc, char **argv)
  43 {
  44   char    *hmmfile;             /* file to read HMM(s) from                */
  45   HMMFILE *hmmfp;               /* opened hmmfile for reading              */
  46   struct plan7_s  *hmm;         /* the HMM to search with                  */
  47   char    *dsq;                 /* digitized target sequence               */
  48   char    *seq;
  49   SQINFO   sqinfo;
  50   int      L;                   /* length of dsq                           */
  51   struct p7trace_s  *tr1;       /* traceback                               */
  52   struct p7trace_s  *tr2;       /* another traceback                       */
  53   int       nseq;
  54   float     sc1, sc2;           /* scores                                  */
  55   int       config;
  56   int       i;
  57
  58   int be_verbose;
  59
  60   char *optname;                /* name of option found by Getopt()         */
  61   char *optarg;                 /* argument found by Getopt()               */
  62   int   optind;                 /* index in argv[]                          */
  63
  64
  65   /***********************************************
  66    * Parse command line
  67    ***********************************************/
  68
  69   be_verbose = FALSE;
  70   hmmfile    = "fn3.hmm";
  71   nseq       = 100;
  72
  73   while (Getopt(argc, argv, OPTIONS, NOPTIONS, usage,
  74                 &optind, &optname, &optarg))  {
  75     if      (strcmp(optname, "-v")       == 0) be_verbose = TRUE;
  76     else if (strcmp(optname, "--hmm")    == 0) hmmfile    = optarg;
  77     else if (strcmp(optname, "-h")       == 0) {
  78       Banner(stdout, banner);
  79       puts(usage);
  80       puts(experts);
  81       exit(0);
  82     }
  83   }
  84   if (argc - optind != 0)
  85     Die("Incorrect number of arguments.\n%s\n", usage);
  86
  87   /***********************************************
  88    * Open HMM file
  89    * Read a single HMM from it.
  90    ***********************************************/
  91
  92   if ((hmmfp = HMMFileOpen(hmmfile, NULL)) == NULL)
  93     Die("Failed to open HMM file %s\n%s", hmmfile, usage);
  94   if (!HMMFileRead(hmmfp, &hmm))
  95     Die("Failed to read any HMMs from %s\n", hmmfile);
  96   if (hmm == NULL)
  97     Die("HMM file %s corrupt or in incorrect format? Parse failed", hmmfile);
  98   Plan7Renormalize(hmm);
  99
 100   /***********************************************
 101    * We cycle through different model configurations.
 102    * For each configuration, we repeat 100 times:
 103    *    - generate a sequence
 104    *    - score it by Viterbi and by SmallViterbi
 105    *    - make sure they give OK and identical results
 106    ***********************************************/
 107
 108   for (config = 1; config <= 5; config++)
 109     {
 110       switch (config) {
 111       case 1: Plan7NakedConfig(hmm);            break;
 112       case 2: Plan7GlobalConfig(hmm);           break;
 113       case 3: Plan7LSConfig(hmm);               break;
 114       case 4: Plan7FSConfig(hmm, 0.5, 0.5);     break;
 115       case 5: Plan7SWConfig(hmm, 0.5, 0.5);     break;
 116       default: Die("never happens");
 117       }
 118       P7Logoddsify(hmm, TRUE);
 119
 120       for (i = 0; i < nseq; i++)
 121         {
 122           EmitSequence(hmm, &dsq, &L, NULL);
 123           sprintf(sqinfo.name, "seq%d", i+1);
 124           sqinfo.len   = L;
 125           sqinfo.flags = SQINFO_NAME | SQINFO_LEN;
 126
 127           sc1 = P7Viterbi(dsq, L, hmm, &tr1);
 128           sc2 = P7SmallViterbi(dsq, L, hmm, &tr2);
 129
 130           if (be_verbose)
 131             {
 132               printf("Viterbi score: %.1f   SmallViterbi: %.1f\n", sc1, sc2);
 133               P7PrintTrace(stdout, tr1, hmm, dsq);
 134               P7PrintTrace(stdout, tr2, hmm, dsq);
 135
 136               seq = DedigitizeSequence(dsq, L);
 137               WriteSeq(stdout, SQFILE_FASTA, seq, &sqinfo);
 138               free(seq);
 139             }
 140
 141           if (sc1 != sc2)
 142             Die("Different scores from normal/small Viterbi");
 143
 144           if (fabs(sc1 - P7TraceScore(hmm, dsq, tr1)) > 0.1)
 145             Die("P7Viterbi score doesn't match its TraceScore");
 146           if (fabs(sc2 - P7TraceScore(hmm, dsq, tr2)) > 0.1)
 147             Die("P7SmallViterbi score doesn't match its TraceScore");
 148
 149           if (! TraceVerify(tr1, hmm->M, L))
 150             Die("TraceVerify() failed for a P7Viterbi trace");
 151           if (! TraceVerify(tr2, hmm->M, L))
 152             Die("TraceVerify() failed for a P7SmallViterbi trace");
 153
 154           if (tr1->tlen != tr2->tlen)
 155             Die("Trace lengths differ for normal/small Viterbi");
 156           if (! TraceCompare(tr1, tr2))
 157             Die("Different traces from normal/small Viterbi");
 158
 159           P7FreeTrace(tr1);
 160           P7FreeTrace(tr2);
 161           free(dsq);
 162         }
 163     }
 164
 165   return EXIT_SUCCESS;
 166 }