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

   1 /* weeviterbi_test.c
   2  * Wed Mar  4 17:30:39 1998
   3  *
   4  * Test driver for Myers/Miller/Hirschberg linear memory Viterbi tracebacks.
   5  *
   6  * RCS $Id: weeviterbi_test.c,v 1.1.1.1 2005/03/22 08:34:47 cmzmasek Exp $
   7  */
   8
   9 #include <stdio.h>
  10 #include <time.h>
  11 #include <math.h>
  12
  13 #include "structs.h"
  14 #include "funcs.h"
  15 #include "globals.h"
  16 #include "squid.h"
  17
  18 static char banner[] = "\
  19 weeviterbi_test : testing of Plan7 Myers/Miller/Hirschberg Viterbi traceback code";
  20
  21 static char usage[] = "\
  22 Usage: testdriver [-options]\n\
  23   Available options are:\n\
  24   -h              : help; display this usage info\n\
  25   -v              : be verbose\n\
  26 ";
  27
  28 static char experts[] = "\
  29   --hmm <f>       : use HMM in file <f>\n\
  30   --seq <f>       : use seq(s) in file <f>\n\
  31 \n";
  32
  33 static struct opt_s OPTIONS[] = {
  34   { "-h",       TRUE,  sqdARG_NONE },
  35   { "-v",       TRUE,  sqdARG_NONE },
  36   { "--hmm",    FALSE, sqdARG_STRING },
  37   { "--seq",    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   char    *seqfile;             /* file to read target sequence(s) from    */
  47   SQFILE   *sqfp;               /* opened seqfile for reading              */
  48   char     *seq;                /* target sequence                         */
  49   SQINFO    sqinfo;             /* optional info for seq                   */
  50   char     *dsq;                /* digitized target sequence               */
  51   struct plan7_s  *hmm;         /* HMM to search with                      */
  52   struct p7trace_s  *t1;        /* standard Viterbi traceback              */
  53   struct p7trace_s  *t2;        /* WeeViterbi traceback                    */
  54   int       nseq;
  55   float     sc1,sc2;            /* scores from Viterbi, WeeViterbi         */
  56
  57   int be_verbose;
  58
  59   char *optname;                /* name of option found by Getopt()         */
  60   char *optarg;                 /* argument found by Getopt()               */
  61   int   optind;                 /* index in argv[]                          */
  62
  63   /***********************************************
  64    * Parse command line
  65    ***********************************************/
  66
  67   be_verbose = FALSE;
  68   hmmfile    = "weeviterbi_test.hmm";
  69   seqfile    = "weeviterbi_test.seq";
  70
  71   while (Getopt(argc, argv, OPTIONS, NOPTIONS, usage,
  72                 &optind, &optname, &optarg))  {
  73     if      (strcmp(optname, "-v")       == 0) be_verbose = TRUE;
  74     else if (strcmp(optname, "--hmm")    == 0) hmmfile    = optarg;
  75     else if (strcmp(optname, "--seq")    == 0) seqfile    = optarg;
  76     else if (strcmp(optname, "-h")       == 0) {
  77       Banner(stdout, banner);
  78       puts(usage);
  79       puts(experts);
  80       exit(0);
  81     }
  82   }
  83   if (argc - optind != 0)
  84     Die("Incorrect number of arguments.\n%s\n", usage);
  85
  86   /***********************************************
  87    * Open test sequence file
  88    ***********************************************/
  89
  90   if ((sqfp = SeqfileOpen(seqfile, SQFILE_UNKNOWN, "BLASTDB")) == NULL)
  91     Die("Failed to open sequence database file %s\n%s\n", seqfile, usage);
  92
  93   /***********************************************
  94    * Open HMM file
  95    * Read a single HMM from it. (Config HMM, if necessary).
  96    ***********************************************/
  97
  98   if ((hmmfp = HMMFileOpen(hmmfile, NULL)) == NULL)
  99     Die("Failed to open HMM file %s\n%s", hmmfile, usage);
 100   if (!HMMFileRead(hmmfp, &hmm))
 101     Die("Failed to read any HMMs from %s\n", hmmfile);
 102   if (hmm == NULL)
 103     Die("HMM file %s corrupt or in incorrect format? Parse failed", hmmfile);
 104   P7Logoddsify(hmm, TRUE);
 105
 106   /***********************************************
 107    * Search HMM against each sequence
 108    ***********************************************/
 109
 110   nseq = 0;
 111   while (ReadSeq(sqfp, sqfp->format, &seq, &sqinfo))
 112     {
 113       nseq++;
 114       dsq = DigitizeSequence(seq, sqinfo.len);
 115
 116       sc1 = P7Viterbi(dsq, sqinfo.len, hmm, &t1);
 117       sc2 = P7WeeViterbi(dsq, sqinfo.len, hmm, &t2);
 118
 119       if (be_verbose)
 120         {
 121           printf("test sequence %d: %s %s\n",
 122                  nseq, sqinfo.name,
 123                  sqinfo.flags & SQINFO_DESC ? sqinfo.desc : "");
 124           printf("** P7Viterbi trace:\n");
 125           P7PrintTrace(stdout, t1, hmm, dsq);
 126           printf("** P7WeeViterbi trace:\n");
 127           P7PrintTrace(stdout, t2, hmm, dsq);
 128         }
 129
 130       if (! TraceVerify(t1, hmm->M, sqinfo.len))
 131         Die("Trace verify failed on Viterbi for seq #%d, %s\n", nseq, sqinfo.name);
 132       if (! TraceVerify(t2, hmm->M, sqinfo.len))
 133         Die("Trace verify failed on WeeViterbi for seq #%d, %s\n", nseq, sqinfo.name);
 134       if (sc1 != sc2)
 135         Die("Scores for the two Viterbi implementations are unequal (%.1f,%.1f)", sc1, sc2);
 136       if (! TraceCompare(t1, t2))
 137         Die("WeeViterbi() trace is not identical to Viterbi() trace");
 138
 139       FreeSequence(seq, &sqinfo);
 140       P7FreeTrace(t1);
 141       P7FreeTrace(t2);
 142       free(dsq);
 143     }
 144
 145   FreePlan7(hmm);
 146   SeqfileClose(sqfp);
 147   HMMFileClose(hmmfp);
 148
 149   return EXIT_SUCCESS;
 150 }