src/gromacs/gmxlib/statistics/gmx_statistics_test.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 4.5
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2008, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Groningen Machine for Chemical Simulation
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <stdio.h>
  40 #include "typedefs.h"
  41 #include "smalloc.h"
  42 #include "vec.h"
  43 #include "gmx_random.h"
  44 #include "gmx_statistics.h"
  45
  46 static void horizontal()
  47 {
  48     gmx_rng_t   rng;
  49     gmx_stats_t straight;
  50     int         i, ok, n = 1000;
  51     real        y, a, b, da, db, aver, sigma, error, chi2, R, *xh, *yh;
  52     FILE       *fp;
  53
  54     rng      = gmx_rng_init(13);
  55     straight = gmx_stats_init();
  56     for (i = 0; (i < n); i++)
  57     {
  58         y = gmx_rng_uniform_real(rng);
  59         if ((ok = gmx_stats_add_point(straight, i, y, 0, 0)) != estatsOK)
  60         {
  61             fprintf(stderr, "%s\n", gmx_stats_message(ok));
  62         }
  63     }
  64     /* Horizontal test */
  65     if ((ok = gmx_stats_get_ase(straight, &aver, &sigma, &error)) != estatsOK)
  66     {
  67         fprintf(stderr, "%s\n", gmx_stats_message(ok));
  68     }
  69     fp = fopen("straight.xvg", "w");
  70     if ((ok = gmx_stats_dump_xy(straight, fp)) != estatsOK)
  71     {
  72         fprintf(stderr, "%s\n", gmx_stats_message(ok));
  73     }
  74     fclose(fp);
  75     printf("Horizontal line: average %g, sigma %g, error %g\n", aver, sigma, error);
  76     if ((ok = gmx_stats_done(straight)) != estatsOK)
  77     {
  78         fprintf(stderr, "%s\n", gmx_stats_message(ok));
  79     }
  80 }
  81
  82 static void line()
  83 {
  84     gmx_rng_t   rng;
  85     gmx_stats_t line;
  86     int         i, dy, ok, n = 1000;
  87     real        y, a, b, da, db, aver, sigma, error, chi2, R, rfit;
  88     const real  a0 = 0.23, b0 = 2.7;
  89     FILE       *fp;
  90
  91     for (dy = 0; (dy < 2); dy++)
  92     {
  93         rng      = gmx_rng_init(13);
  94         line     = gmx_stats_init();
  95         for (i = 0; (i < n); i++)
  96         {
  97             y = a0*i+b0+50*(gmx_rng_uniform_real(rng)-0.5);
  98             if ((ok = gmx_stats_add_point(line, i, y, 0, dy*0.1)) != estatsOK)
  99             {
 100                 fprintf(stderr, "%s\n", gmx_stats_message(ok));
 101             }
 102         }
 103         /* Line with slope test */
 104         if ((ok = gmx_stats_get_ab(line, elsqWEIGHT_NONE, &a, &b, &da, &db, &chi2, &rfit)) != estatsOK)
 105         {
 106             fprintf(stderr, "%s\n", gmx_stats_message(ok));
 107         }
 108         if ((ok = gmx_stats_get_corr_coeff(line, &R)) != estatsOK)
 109         {
 110             fprintf(stderr, "%s\n", gmx_stats_message(ok));
 111         }
 112         if (dy == 0)
 113         {
 114             fp = fopen("line0.xvg", "w");
 115         }
 116         else
 117         {
 118             fp = fopen("line1.xvg", "w");
 119         }
 120         if ((ok = gmx_stats_dump_xy(line, fp)) != estatsOK)
 121         {
 122             fprintf(stderr, "%s\n", gmx_stats_message(ok));
 123         }
 124         fclose(fp);
 125         printf("Line with eqn. y = %gx + %g with noise%s\n", a0, b0,
 126                (dy == 0) ? "" : " and uncertainties");
 127         printf("Found: a = %g +/- %g, b = %g +/- %g\n", a, da, b, db);
 128         if ((ok = gmx_stats_done(line)) != estatsOK)
 129         {
 130             fprintf(stderr, "%s\n", gmx_stats_message(ok));
 131         }
 132         gmx_rng_destroy(rng);
 133     }
 134 }
 135
 136 static void histogram()
 137 {
 138     gmx_rng_t   rng;
 139     gmx_stats_t camel;
 140     int         i, ok, n = 1000, norm;
 141     real        y, a, b, da, db, aver, sigma, error, chi2, R, *xh, *yh;
 142     const real  a0 = 0.23, b0 = 2.7;
 143     FILE       *fp;
 144     char        fn[256];
 145
 146     for (norm = 0; (norm < 2); norm++)
 147     {
 148         rng      = gmx_rng_init(13);
 149         camel    = gmx_stats_init();
 150         for (i = 0; (i < n); i++)
 151         {
 152             y = sqr(gmx_rng_uniform_real(rng));
 153             if ((ok = gmx_stats_add_point(camel, i, y+1, 0, 0)) != estatsOK)
 154             {
 155                 fprintf(stderr, "%s\n", gmx_stats_message(ok));
 156             }
 157             y = sqr(gmx_rng_uniform_real(rng));
 158             if ((ok = gmx_stats_add_point(camel, i+0.5, y+2, 0, 0)) != estatsOK)
 159             {
 160                 fprintf(stderr, "%s\n", gmx_stats_message(ok));
 161             }
 162         }
 163         /* Histogram test */
 164         if ((ok = gmx_stats_make_histogram(camel, 0, 101, norm, &xh, &yh)) != estatsOK)
 165         {
 166             fprintf(stderr, "%s\n", gmx_stats_message(ok));
 167         }
 168         sprintf(fn, "histo%d-data.xvg", norm);
 169         fp = fopen(fn, "w");
 170         gmx_stats_dump_xy(camel, fp);
 171         fclose(fp);
 172         sprintf(fn, "histo%d.xvg", norm);
 173         fp = fopen(fn, "w");
 174         for (i = 0; (i < 101); i++)
 175         {
 176             fprintf(fp, "%12g  %12g\n", xh[i], yh[i]);
 177         }
 178         fclose(fp);
 179         sfree(xh);
 180         sfree(yh);
 181     }
 182 }
 183
 184 int main(int argc, char *argv[])
 185 {
 186     line();
 187     horizontal();
 188     histogram();
 189
 190     return 0;
 191 }