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

   1 /* evd_test.c
   2  * SRE, Wed Nov 12 11:17:27 1997 [St. Louis]
   3  *
   4  * Test driver for EVD distribution support in histogram.c
   5  * Generates random EVD samples; fits them; checks fitted mu, lambda
   6  * against parametric mu, lambda. If they differ badly, calls Die().
   7  * If OK, returns EXIT_SUCCESS.
   8  *
   9  * RCS $Id: evd_test.c,v 1.1.1.1 2005/03/22 08:34:45 cmzmasek Exp $
  10  */
  11
  12
  13 #include <stdio.h>
  14 #include <time.h>
  15 #include <math.h>
  16
  17 #include "structs.h"
  18 #include "funcs.h"
  19 #include "globals.h"
  20 #include "squid.h"
  21
  22 #ifdef MEMDEBUG
  23 #include "dbmalloc.h"
  24 #endif
  25
  26 static char banner[] = "\
  27 evd_test : testing of EVD code in histogram.c";
  28
  29 static char usage[] = "\
  30 Usage: testdriver [-options]\n\
  31   Available options are:\n\
  32   -h              : help; display this usage info\n\
  33   -c <x>          : censor data below <x>\n\
  34   -e <n>          : sample <n> times from EVD\n\
  35   -g <n>          : add <n> Gaussian samples of \"noise\"\n\
  36   -n <n>          : set number of trials to <n>\n\
  37   -s <n>          : set random seed to <n>\n\
  38   -v              : be verbose (default is to simply exit with status 1 or 0)\n\
  39 ";
  40
  41 static char experts[] = "\
  42   --xmgr <file>   : save graphical data to <file>\n\
  43   --hist          : fit to histogram instead of raw samples\n\
  44   --loglog <file> : save log log regression line to <file>\n\
  45   --regress       : do old-style linear regression fit, not ML\n\
  46   --mu <x>        : set EVD mu to <x>\n\
  47   --lambda <x>    : set EVD lambda to <x>\n\
  48   --mean <x>      : set Gaussian mean to <x>\n\
  49   --sd   <x>      : set Gaussian std. dev. to <x>\n\
  50 \n";
  51
  52 static struct opt_s OPTIONS[] = {
  53   { "-h",       TRUE,  sqdARG_NONE  },
  54   { "-c",       TRUE,  sqdARG_FLOAT },
  55   { "-e",       TRUE,  sqdARG_INT },
  56   { "-g",       TRUE,  sqdARG_INT },
  57   { "-n",       TRUE,  sqdARG_INT },
  58   { "-s",       TRUE,  sqdARG_INT   },
  59   { "-v",       TRUE,  sqdARG_NONE  },
  60   { "--xmgr",   FALSE, sqdARG_STRING},
  61   { "--hist",   FALSE, sqdARG_NONE},
  62   { "--loglog", FALSE, sqdARG_STRING},
  63   { "--regress",FALSE, sqdARG_NONE},
  64   { "--mu",     FALSE, sqdARG_FLOAT},
  65   { "--lambda", FALSE, sqdARG_FLOAT},
  66   { "--mean",   FALSE, sqdARG_FLOAT},
  67   { "--sd",     FALSE, sqdARG_FLOAT},
  68 };
  69 #define NOPTIONS (sizeof(OPTIONS) / sizeof(struct opt_s))
  70
  71 int
  72 main(int argc, char **argv)
  73 {
  74   struct histogram_s *h;        /* histogram structure          */
  75   int ntrials;                  /* number of different fits     */
  76   int be_verbose;               /* option: TRUE to show output  */
  77   int seed;                     /* option: random number seed   */
  78   int   nevd;                   /* # of samples from EVD        */
  79   float mu;                     /* EVD mu parameter             */
  80   float lambda;                 /* EVD lambda parameter         */
  81   int   ngauss;                 /* # of samples from Gaussian   */
  82   float mean;                   /* Gaussian "noise" mean        */
  83   float sd;                     /* Gaussian "noise" std. dev.   */
  84   float x;                      /* a random sample              */
  85   int   i, idx;
  86   float *val;                   /* array of samples             */
  87   float mlmu;                   /* estimate of mu               */
  88   float mllambda;               /* estimate of lambda           */
  89
  90   char *xmgrfile;               /* output file for XMGR graph data */
  91   char *logfile;                /* output file for regression line */
  92   FILE *xmgrfp;                 /* open output file                */
  93   FILE *logfp;                  /* open log log file               */
  94   int   do_ml;                  /* TRUE to do a max likelihood fit */
  95   int   fit_hist;               /* TRUE to fit histogram instead of samples */
  96   int   censoring;              /* TRUE to left-censor the data    */
  97   float censorlevel;            /* value to censor at              */
  98
  99   char *optname;                /* name of option found by Getopt()         */
 100   char *optarg;                 /* argument found by Getopt()               */
 101   int   optind;                 /* index in argv[]                          */
 102
 103 #ifdef MEMDEBUG
 104   unsigned long histid1, histid2, orig_size, current_size;
 105   orig_size = malloc_inuse(&histid1);
 106   fprintf(stderr, "[... memory debugging is ON ...]\n");
 107 #endif
 108
 109   /***********************************************
 110    * Parse command line
 111    ***********************************************/
 112   be_verbose = FALSE;
 113   seed       = (int) time ((time_t *) NULL);
 114   ntrials    = 1;
 115   nevd       = 1000;
 116   mu         = -20.0;
 117   lambda     = 0.4;
 118   ngauss     = 0;
 119   mean       = 20.;
 120   sd         = 20.;
 121   xmgrfile   = NULL;
 122   logfile    = NULL;
 123   xmgrfp     = NULL;
 124   logfp      = NULL;
 125   do_ml      = TRUE;
 126   censoring  = FALSE;
 127   censorlevel= 0.;
 128   fit_hist   = FALSE;
 129
 130   while (Getopt(argc, argv, OPTIONS, NOPTIONS, usage,
 131                 &optind, &optname, &optarg))  {
 132     if      (strcmp(optname, "-e")       == 0) { nevd       = atoi(optarg); }
 133     else if (strcmp(optname, "-c")       == 0) { censoring  = TRUE;
 134                                                  censorlevel= atof(optarg); }
 135     else if (strcmp(optname, "-g")       == 0) { ngauss     = atoi(optarg); }
 136     else if (strcmp(optname, "-n")       == 0) { ntrials    = atoi(optarg); }
 137     else if (strcmp(optname, "-s")       == 0) { seed       = atoi(optarg); }
 138     else if (strcmp(optname, "-v")       == 0) { be_verbose = TRUE;         }
 139     else if (strcmp(optname, "--xmgr")   == 0) { xmgrfile   = optarg; }
 140     else if (strcmp(optname, "--hist")   == 0) { fit_hist   = TRUE; }
 141     else if (strcmp(optname, "--loglog") == 0) { logfile    = optarg; }
 142     else if (strcmp(optname, "--regress")== 0) { do_ml      = FALSE; }
 143     else if (strcmp(optname, "--mu")     == 0) { mu         = atof(optarg); }
 144     else if (strcmp(optname, "--lambda") == 0) { lambda     = atof(optarg); }
 145     else if (strcmp(optname, "--mean")   == 0) { mean       = atof(optarg); }
 146     else if (strcmp(optname, "--sd")     == 0) { sd         = atof(optarg); }
 147     else if (strcmp(optname, "-h")       == 0) {
 148       Banner(stdout, banner);
 149       puts(usage);
 150       puts(experts);
 151       exit(0);
 152     }
 153   }
 154   if (argc - optind != 0)
 155     Die("Incorrect number of arguments.\n%s\n", usage);
 156
 157   sre_srandom(seed);
 158
 159   /****************************************************************
 160    * Print options
 161    ****************************************************************/
 162
 163   if (be_verbose)
 164     {
 165       puts("--------------------------------------------------------");
 166       printf("EVD samples    = %d\n", nevd);
 167       printf("mu, lambda     = %f, %f\n", mu, lambda);
 168       if (ngauss > 0) {
 169         printf("Gaussian noise = %d\n", ngauss);
 170         printf("mean, sd       = %f, %f\n", mean, sd);
 171       }
 172       if (censoring) printf("pre-censoring  = ON, at %f\n", censorlevel);
 173       printf("total trials   = %d\n", ntrials);
 174       printf("random seed    = %d\n", seed);
 175       printf("fit method     = %s\n", do_ml ? "ML" : "linear regression");
 176       printf("fit is to      = %s\n", fit_hist ? "histogram" : "list");
 177       puts("--------------------------------------------------------");
 178     }
 179
 180   if (xmgrfile != NULL)
 181     if ((xmgrfp = fopen(xmgrfile, "w")) == NULL)
 182       Die("Failed to open output file %s", xmgrfile);
 183   if (logfile != NULL)
 184     if ((logfp = fopen(logfile, "w")) == NULL)
 185       Die("Failed to open output file %s", logfile);
 186
 187   /* Generate random EVD "signal" (and Gaussian "noise")
 188    * samples and put them in the histogram
 189    */
 190   while (ntrials--)
 191     {
 192       val = MallocOrDie(sizeof(double) * (nevd+ngauss));
 193       h   = AllocHistogram(-20, 20, 10);
 194
 195                                 /* EVD signal */
 196       idx = 0;
 197       for (i = 0; i < nevd; i++)
 198         {
 199           x = EVDrandom(mu, lambda);
 200           if (! censoring || x > censorlevel)
 201             {
 202               AddToHistogram(h, x);
 203               val[idx] = x;
 204               idx++;
 205             }
 206         }
 207                                 /* Gaussian noise */
 208       for (; i < nevd + ngauss; i++)
 209         {
 210           x = Gaussrandom(mean, sd);
 211           if (! censoring || x > censorlevel)
 212             {
 213               AddToHistogram(h, x);
 214               val[idx] = x;
 215               idx++;
 216             }
 217         }
 218
 219       if (do_ml)
 220         {
 221
 222           if (censoring)
 223             {
 224               if (be_verbose)
 225                 printf("I have censored the data at %f: %d observed, %d censored\n", censorlevel, idx, (nevd+ngauss)-idx);
 226
 227               EVDCensoredFit(val, NULL, idx,
 228                              (nevd+ngauss)-idx, censorlevel,
 229                              &mlmu, &mllambda);
 230               ExtremeValueSetHistogram(h, (float) mlmu, (float) mllambda,
 231                                        censorlevel, h->highscore, 1);
 232             }
 233           else
 234             {
 235               if (fit_hist)
 236                 {
 237                   ExtremeValueFitHistogram(h, TRUE, 20.);
 238                 }
 239               else
 240                 {
 241                   EVDMaxLikelyFit(val, NULL, idx, &mlmu, &mllambda);
 242                   ExtremeValueSetHistogram(h, (float) mlmu, (float) mllambda,
 243                                            h->lowscore, h->highscore, 2);
 244                 }
 245             }
 246         }
 247       else
 248         EVDBasicFit(h);
 249
 250       if (be_verbose) {
 251         printf("%f\tmu\n",     h->param[EVD_MU]);
 252         printf("%f\tlambda\n", h->param[EVD_LAMBDA]);
 253         printf("%f\t%% error on mu\n",
 254                fabs(100. * (h->param[EVD_MU] - mu) / mu));
 255         printf("%f\t%% error on lambda\n",
 256                fabs(100. * (h->param[EVD_LAMBDA] - lambda) / lambda));
 257         printf("%f\tchi-squared P value\n", h->chip);
 258       }
 259       if (xmgrfp != NULL) PrintXMGRHistogram(xmgrfp, h);
 260       /*      if (xmgrfp != NULL) PrintXMGRDistribution(xmgrfp, h); */
 261       if (logfp  != NULL) PrintXMGRRegressionLine(logfp, h);
 262
 263       /* Generate the expected lines: sets 5,7 of xmgrfile (manually delete 4,6)
 264        *                              set 3 of loglogfile  (manually delete 2)
 265        */
 266       ExtremeValueSetHistogram(h, mu, lambda, h->lowscore, h->highscore, 0);
 267       if (xmgrfp != NULL) PrintXMGRHistogram(xmgrfp, h);
 268       /*      if (xmgrfp != NULL) PrintXMGRDistribution(xmgrfp, h); */
 269       if (logfp  != NULL) PrintXMGRRegressionLine(logfp, h);
 270
 271       /* Do the internal test.
 272        * Criterion: on a 1000 sample EVD of u = -40 and lambda = 0.4,
 273        * estimate u to within +/- 2 and lambda to within +/- 0.05.
 274        */
 275       if (fabs(h->param[EVD_MU] - mu) > 2.)
 276         Die("evd_test: tolerance to mu exceeded (%f)",
 277             fabs(h->param[EVD_MU] - mu));
 278       if (fabs(h->param[EVD_LAMBDA] - lambda) > 0.05)
 279         Die("evd_test: tolerance to lambda exceeded (%f)",
 280             fabs(h->param[EVD_LAMBDA] - lambda));
 281
 282       FreeHistogram(h);
 283       free(val);
 284     }
 285
 286 #ifdef MEMDEBUG
 287   current_size = malloc_inuse(&histid2);
 288   if (current_size != orig_size) Die("evd_test failed memory test");
 289   else fprintf(stderr, "[No memory leaks.]\n");
 290 #endif
 291
 292   if (xmgrfp != NULL) fclose(xmgrfp);
 293   if (logfp != NULL)  fclose(logfp);
 294   return EXIT_SUCCESS;
 295 }