src/contrib/gmx_stats_test.c

   1 #include <stdio.h>
   2 #include "typedefs.h"
   3 #include "gromacs/utility/smalloc.h"
   4 #include "gromacs/math/vec.h"
   5 #include "gmx_random.h"
   6 #include "gmx_statistics.h"
   7
   8 static void horizontal()
   9 {
  10   gmx_rng_t   rng;
  11   gmx_stats_t straight;
  12   int i,ok,n=1000;
  13   real y,a,b,da,db,aver,sigma,error,chi2,R,*xh,*yh;
  14   FILE *fp;
  15
  16   rng      = gmx_rng_init(13);
  17   straight = gmx_stats_init();
  18   for(i=0; (i<n); i++) {
  19     y = gmx_rng_uniform_real(rng);
  20     if ((ok = gmx_stats_add_point(straight,i,y,0,0)) != estatsOK)
  21       fprintf(stderr,"%s\n",gmx_stats_message(ok));
  22   }
  23   /* Horizontal test */
  24   if ((ok = gmx_stats_get_ase(straight,&aver,&sigma,&error)) != estatsOK)
  25     fprintf(stderr,"%s\n",gmx_stats_message(ok));
  26   fp = fopen("straight.xvg","w");
  27   if ((ok = gmx_stats_dump_xy(straight,fp)) != estatsOK)
  28     fprintf(stderr,"%s\n",gmx_stats_message(ok));
  29   fclose(fp);
  30   printf("Horizontal line: average %g, sigma %g, error %g\n",aver,sigma,error);
  31   if ((ok = gmx_stats_done(straight)) != estatsOK)
  32     fprintf(stderr,"%s\n",gmx_stats_message(ok));
  33 }
  34
  35 static void line()
  36 {
  37   gmx_rng_t   rng;
  38   gmx_stats_t line;
  39   int  i,dy,ok,n=1000;
  40   real y,a,b,da,db,aver,sigma,error,chi2,R,rfit;
  41   const real a0=0.23,b0=2.7;
  42   FILE *fp;
  43
  44   for(dy=0; (dy<2); dy++) {
  45     rng      = gmx_rng_init(13);
  46     line     = gmx_stats_init();
  47     for(i=0; (i<n); i++) {
  48       y = a0*i+b0+50*(gmx_rng_uniform_real(rng)-0.5);
  49       if ((ok = gmx_stats_add_point(line,i,y,0,dy*0.1)) != estatsOK)
  50         fprintf(stderr,"%s\n",gmx_stats_message(ok));
  51     }
  52     /* Line with slope test */
  53     if ((ok = gmx_stats_get_ab(line,elsqWEIGHT_NONE,&a,&b,&da,&db,&chi2,&rfit)) != estatsOK)
  54       fprintf(stderr,"%s\n",gmx_stats_message(ok));
  55     if ((ok = gmx_stats_get_corr_coeff(line,&R)) != estatsOK)
  56       fprintf(stderr,"%s\n",gmx_stats_message(ok));
  57     if (dy == 0)
  58       fp = fopen("line0.xvg","w");
  59     else
  60       fp = fopen("line1.xvg","w");
  61     if ((ok = gmx_stats_dump_xy(line,fp)) != estatsOK)
  62       fprintf(stderr,"%s\n",gmx_stats_message(ok));
  63     fclose(fp);
  64     printf("Line with eqn. y = %gx + %g with noise%s\n",a0,b0,
  65            (dy == 0) ? "" : " and uncertainties");
  66     printf("Found: a = %g +/- %g, b = %g +/- %g\n",a,da,b,db);
  67     if ((ok = gmx_stats_done(line)) != estatsOK)
  68       fprintf(stderr,"%s\n",gmx_stats_message(ok));
  69     gmx_rng_destroy(rng);
  70   }
  71 }
  72
  73 static void histogram()
  74 {
  75   gmx_rng_t   rng;
  76   gmx_stats_t camel;
  77   int i,ok,n=1000,norm;
  78   real y,a,b,da,db,aver,sigma,error,chi2,R,*xh,*yh;
  79   const real a0=0.23,b0=2.7;
  80   FILE *fp;
  81   char fn[256];
  82
  83   for(norm=0; (norm<2); norm++) {
  84     rng      = gmx_rng_init(13);
  85     camel    = gmx_stats_init();
  86     for(i=0; (i<n); i++) {
  87       y = sqr(gmx_rng_uniform_real(rng));
  88       if ((ok = gmx_stats_add_point(camel,i,y+1,0,0)) != estatsOK)
  89         fprintf(stderr,"%s\n",gmx_stats_message(ok));
  90       y = sqr(gmx_rng_uniform_real(rng));
  91       if ((ok = gmx_stats_add_point(camel,i+0.5,y+2,0,0)) != estatsOK)
  92         fprintf(stderr,"%s\n",gmx_stats_message(ok));
  93     }
  94     /* Histogram test */
  95     if ((ok = gmx_stats_make_histogram(camel,0,101,norm,&xh,&yh)) != estatsOK)
  96       fprintf(stderr,"%s\n",gmx_stats_message(ok));
  97     sprintf(fn,"histo%d-data.xvg",norm);
  98     fp = fopen(fn,"w");
  99     gmx_stats_dump_xy(camel,fp);
 100     fclose(fp);
 101     sprintf(fn,"histo%d.xvg",norm);
 102     fp = fopen(fn,"w");
 103     for(i=0; (i<101); i++)
 104       fprintf(fp,"%12g  %12g\n",xh[i],yh[i]);
 105     fclose(fp);
 106     sfree(xh);
 107     sfree(yh);
 108   }
 109 }
 110
 111 int main(int argc,char *argv[])
 112 {
 113   line();
 114   horizontal();
 115   histogram();
 116
 117   return 0;
 118 }