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

   1 /* parsingviterbi_test.c
   2  * Wed Mar  4 15:07:37 1998
   3  * cp trace_test.c ../src/testdriver.c; cd ../src; make testdriver
   4  *
   5  * Test driver for P7ParsingViterbi(); alignment in linear memory.
   6  *
   7  * CVS $Id: parsingviterbi_test.c,v 1.1.1.1 2005/03/22 08:34:47 cmzmasek Exp $
   8  */
   9
  10 #include <stdio.h>
  11 #include <time.h>
  12 #include <math.h>
  13
  14 #include "structs.h"
  15 #include "funcs.h"
  16 #include "globals.h"
  17 #include "squid.h"
  18
  19 static char banner[] = "\
  20 parsingviterbi_test : testing of Plan7 linear memory alignment code";
  21
  22 static char usage[] = "\
  23 Usage: parsingviterbi_test [-options]\n\
  24   Available options are:\n\
  25   -h              : help; display this usage info\n\
  26   -v              : be verbose\n\
  27 ";
  28
  29 static char experts[] = "\
  30 \n";
  31
  32 static struct opt_s OPTIONS[] = {
  33   { "-h",       TRUE,  sqdARG_NONE },
  34   { "-v",       TRUE,  sqdARG_NONE },
  35 };
  36 #define NOPTIONS (sizeof(OPTIONS) / sizeof(struct opt_s))
  37
  38 int
  39 main(int argc, char **argv)
  40 {
  41   char    *hmmfile;             /* file to read HMM(s) from                */
  42   HMMFILE *hmmfp;               /* opened hmmfile for reading              */
  43   char    *seqfile;             /* file to read target sequence(s) from    */
  44   SQFILE   *sqfp;               /* opened seqfile for reading              */
  45   char     *seq;                /* target sequence                         */
  46   SQINFO    sqinfo;             /* optional info for seq                   */
  47   char     *dsq;                /* digitized target sequence               */
  48   struct plan7_s  *hmm;         /* HMM to search with                      */
  49   struct p7trace_s  *tr1;       /* traceback from P7Viterbi()              */
  50   struct p7trace_s  *tr2;       /* traceback from P7ParsingViterbi()       */
  51   int       nseq;
  52   float     sc1, sc2;           /* scores from Viterbi, ParsingViterbi()   */
  53
  54   struct p7trace_s **tarr;      /* array of decomposed Viterbi traces      */
  55   int    ntr;                   /* number of traces                        */
  56   int    i1,i2,k1,k2;           /* starts, stops in seq, model for Viterbi */
  57   int    idx;                   /* index of a decomposed trace             */
  58
  59   int be_verbose;
  60
  61   char *optname;                /* name of option found by Getopt()         */
  62   char *optarg;                 /* argument found by Getopt()               */
  63   int   optind;                 /* index in argv[]                          */
  64
  65   /***********************************************
  66    * Parse command line
  67    ***********************************************/
  68
  69   be_verbose = FALSE;
  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, "-h")       == 0) {
  75       Banner(stdout, banner);
  76       puts(usage);
  77       puts(experts);
  78       exit(0);
  79     }
  80   }
  81   if (argc - optind != 0)
  82     Die("Incorrect number of arguments.\n%s\n", usage);
  83
  84   hmmfile = "fn3.hmm";
  85   seqfile = "titin.fa";
  86
  87   /***********************************************
  88    * Open test sequence file
  89    ***********************************************/
  90
  91   if ((sqfp = SeqfileOpen(seqfile, SQFILE_UNKNOWN, "BLASTDB")) == NULL)
  92     Die("Failed to open sequence database file %s\n%s\n", seqfile, usage);
  93
  94   /***********************************************
  95    * Open HMM file
  96    * Read a single HMM from it. (Config HMM, if necessary).
  97    ***********************************************/
  98
  99   if ((hmmfp = HMMFileOpen(hmmfile, NULL)) == NULL)
 100     Die("Failed to open HMM file %s\n%s", hmmfile, usage);
 101   if (!HMMFileRead(hmmfp, &hmm))
 102     Die("Failed to read any HMMs from %s\n", hmmfile);
 103   if (hmm == NULL)
 104     Die("HMM file %s corrupt or in incorrect format? Parse failed", hmmfile);
 105   P7Logoddsify(hmm, TRUE);
 106
 107   /***********************************************
 108    * Search HMM against each sequence, using both
 109    * normal Viterbi and P7ParsingViterbi.
 110    ***********************************************/
 111
 112   nseq = 0;
 113   while (ReadSeq(sqfp, sqfp->format, &seq, &sqinfo))
 114     {
 115       nseq++;
 116       dsq = DigitizeSequence(seq, sqinfo.len);
 117
 118       sc1  = P7Viterbi(dsq, sqinfo.len, hmm, &tr1);
 119       sc2  = P7ParsingViterbi(dsq, sqinfo.len, hmm, &tr2);
 120
 121       if (be_verbose)
 122         {
 123           printf("test sequence %d: %s %s\n",
 124                  nseq, sqinfo.name,
 125                  sqinfo.flags & SQINFO_DESC ? sqinfo.desc : "");
 126           for (idx = 0; idx < tr2->tlen; idx++)
 127             printf("%1s  %d\n", Statetype(tr2->statetype[idx]), tr2->pos[idx]);
 128         }
 129
 130       if (sc1 != sc2)
 131         Die("Scores for the two Viterbi implementations are unequal (%d,%d)", sc1, sc2);
 132
 133       TraceDecompose(tr1, &tarr, &ntr);
 134       if (ntr == 0)
 135         Die("ntr == 0 can't happen");
 136       if (ntr != (tr2->tlen/2) -1)
 137         Die("# of domains for the two Viterbi implementations are unequal (%d, %d)",
 138             ntr, (tr2->tlen/2) -1);
 139
 140       for (idx = 0; idx < ntr; idx++)
 141         {
 142           TraceSimpleBounds(tarr[idx], &i1, &i2, &k1, &k2);
 143
 144           if (i1 != tr2->pos[idx*2 + 1] + 1)
 145             Die("Start positions %d and %d disagree for domain %d\n",
 146                 i1, tr2->pos[idx*2 + 1] + 1, idx);
 147           if (i2 != tr2->pos[idx*2 + 2])
 148             Die("End positions %d and %d disagree for domain %d\n",
 149                 i2, tr2->pos[idx*2 + 2], idx);
 150         }
 151
 152
 153       for (idx = 0; idx < ntr; idx++)
 154         P7FreeTrace(tarr[idx]);
 155       free(tarr);
 156       FreeSequence(seq, &sqinfo);
 157       P7FreeTrace(tr1);
 158       P7FreeTrace(tr2);
 159       free(dsq);
 160     }
 161
 162   FreePlan7(hmm);
 163   HMMFileClose(hmmfp);
 164   SeqfileClose(sqfp);
 165
 166   return EXIT_SUCCESS;
 167 }