src/gmxlib/statistics/gmx_statistics_test.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2008, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <stdio.h>
  43 #include "typedefs.h"
  44 #include "smalloc.h"
  45 #include "vec.h"
  46 #include "gmx_random.h"
  47 #include "gmx_statistics.h"
  48
  49 static void horizontal()
  50 {
  51   gmx_rng_t   rng;
  52   gmx_stats_t straight;
  53   int i,ok,n=1000;
  54   real y,a,b,da,db,aver,sigma,error,chi2,R,*xh,*yh;
  55   FILE *fp;
  56
  57   rng      = gmx_rng_init(13);
  58   straight = gmx_stats_init();
  59   for(i=0; (i<n); i++) {
  60     y = gmx_rng_uniform_real(rng);
  61     if ((ok = gmx_stats_add_point(straight,i,y,0,0)) != estatsOK)
  62       fprintf(stderr,"%s\n",gmx_stats_message(ok));
  63   }
  64   /* Horizontal test */
  65   if ((ok = gmx_stats_get_ase(straight,&aver,&sigma,&error)) != estatsOK)
  66     fprintf(stderr,"%s\n",gmx_stats_message(ok));
  67   fp = fopen("straight.xvg","w");
  68   if ((ok = gmx_stats_dump_xy(straight,fp)) != estatsOK)
  69     fprintf(stderr,"%s\n",gmx_stats_message(ok));
  70   fclose(fp);
  71   printf("Horizontal line: average %g, sigma %g, error %g\n",aver,sigma,error);
  72   if ((ok = gmx_stats_done(straight)) != estatsOK)
  73     fprintf(stderr,"%s\n",gmx_stats_message(ok));
  74 }
  75
  76 static void line()
  77 {
  78   gmx_rng_t   rng;
  79   gmx_stats_t line;
  80   int  i,dy,ok,n=1000;
  81   real y,a,b,da,db,aver,sigma,error,chi2,R,rfit;
  82   const real a0=0.23,b0=2.7;
  83   FILE *fp;
  84
  85   for(dy=0; (dy<2); dy++) {
  86     rng      = gmx_rng_init(13);
  87     line     = gmx_stats_init();
  88     for(i=0; (i<n); i++) {
  89       y = a0*i+b0+50*(gmx_rng_uniform_real(rng)-0.5);
  90       if ((ok = gmx_stats_add_point(line,i,y,0,dy*0.1)) != estatsOK)
  91         fprintf(stderr,"%s\n",gmx_stats_message(ok));
  92     }
  93     /* Line with slope test */
  94     if ((ok = gmx_stats_get_ab(line,elsqWEIGHT_NONE,&a,&b,&da,&db,&chi2,&rfit)) != estatsOK)
  95       fprintf(stderr,"%s\n",gmx_stats_message(ok));
  96     if ((ok = gmx_stats_get_corr_coeff(line,&R)) != estatsOK)
  97       fprintf(stderr,"%s\n",gmx_stats_message(ok));
  98     if (dy == 0)
  99       fp = fopen("line0.xvg","w");
 100     else
 101       fp = fopen("line1.xvg","w");
 102     if ((ok = gmx_stats_dump_xy(line,fp)) != estatsOK)
 103       fprintf(stderr,"%s\n",gmx_stats_message(ok));
 104     fclose(fp);
 105     printf("Line with eqn. y = %gx + %g with noise%s\n",a0,b0,
 106            (dy == 0) ? "" : " and uncertainties");
 107     printf("Found: a = %g +/- %g, b = %g +/- %g\n",a,da,b,db);
 108     if ((ok = gmx_stats_done(line)) != estatsOK)
 109       fprintf(stderr,"%s\n",gmx_stats_message(ok));
 110     gmx_rng_destroy(rng);
 111   }
 112 }
 113
 114 static void histogram()
 115 {
 116   gmx_rng_t   rng;
 117   gmx_stats_t camel;
 118   int i,ok,n=1000,norm;
 119   real y,a,b,da,db,aver,sigma,error,chi2,R,*xh,*yh;
 120   const real a0=0.23,b0=2.7;
 121   FILE *fp;
 122   char fn[256];
 123
 124   for(norm=0; (norm<2); norm++) {
 125     rng      = gmx_rng_init(13);
 126     camel    = gmx_stats_init();
 127     for(i=0; (i<n); i++) {
 128       y = sqr(gmx_rng_uniform_real(rng));
 129       if ((ok = gmx_stats_add_point(camel,i,y+1,0,0)) != estatsOK)
 130         fprintf(stderr,"%s\n",gmx_stats_message(ok));
 131       y = sqr(gmx_rng_uniform_real(rng));
 132       if ((ok = gmx_stats_add_point(camel,i+0.5,y+2,0,0)) != estatsOK)
 133         fprintf(stderr,"%s\n",gmx_stats_message(ok));
 134     }
 135     /* Histogram test */
 136     if ((ok = gmx_stats_make_histogram(camel,0,101,norm,&xh,&yh)) != estatsOK)
 137       fprintf(stderr,"%s\n",gmx_stats_message(ok));
 138     sprintf(fn,"histo%d-data.xvg",norm);
 139     fp = fopen(fn,"w");
 140     gmx_stats_dump_xy(camel,fp);
 141     fclose(fp);
 142     sprintf(fn,"histo%d.xvg",norm);
 143     fp = fopen(fn,"w");
 144     for(i=0; (i<101); i++)
 145       fprintf(fp,"%12g  %12g\n",xh[i],yh[i]);
 146     fclose(fp);
 147     sfree(xh);
 148     sfree(yh);
 149   }
 150 }
 151
 152 int main(int argc,char *argv[])
 153 {
 154   line();
 155   horizontal();
 156   histogram();
 157
 158   return 0;
 159 }