include/gmx_statistics.h

   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,2013, 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
  39 #ifndef _GMX_STATS_H
  40 #define _GMX_STATS_H
  41 #include "visibility.h"
  42
  43 #ifdef __cplusplus
  44 extern "C" {
  45 #endif
  46
  47 #include "typedefs.h"
  48
  49 typedef struct gmx_stats *gmx_stats_t;
  50
  51 /* Error codes returned by the routines */
  52 enum {
  53     estatsOK, estatsNO_POINTS, estatsNO_MEMORY, estatsERROR,
  54     estatsINVALID_INPUT, estatsNOT_IMPLEMENTED, estatsNR
  55 };
  56
  57 enum {
  58     elsqWEIGHT_NONE, elsqWEIGHT_X, elsqWEIGHT_Y,
  59     elsqWEIGHT_XY, elsqWEIGHT_NR
  60 };
  61
  62 enum {
  63     ehistoX, ehistoY, ehistoNR
  64 };
  65
  66 GMX_LIBGMX_EXPORT
  67 gmx_stats_t gmx_stats_init();
  68
  69 GMX_LIBGMX_EXPORT
  70 int gmx_stats_done(gmx_stats_t stats);
  71
  72 /* Remove outliers from a straight line, where level in units of
  73    sigma. Level needs to be larger than one obviously. */
  74 int gmx_stats_remove_outliers(gmx_stats_t stats, double level);
  75
  76 GMX_LIBGMX_EXPORT
  77 int gmx_stats_add_point(gmx_stats_t stats, double x, double y,
  78                         double dx, double dy);
  79
  80 /* The arrays dx and dy may be NULL if no uncertainties are available,
  81    in that case zero uncertainties will be assumed. */
  82 int gmx_stats_add_points(gmx_stats_t stats, int n, real *x, real *y,
  83                          real *dx, real *dy);
  84
  85 /* Return the data points one by one. Return estatsOK while there are
  86    more points, and returns estatsNOPOINTS when the last point has
  87    been returned. Should be used in a while loop. Variables for either
  88    pointer may be NULL, in which case the routine can be used as an
  89    expensive point counter. */
  90 GMX_LIBGMX_EXPORT
  91 int gmx_stats_get_point(gmx_stats_t stats, real *x, real *y,
  92                         real *dx, real *dy);
  93
  94 /* Fit the data to y = ax + b, possibly weighted, if uncertainties
  95    have been input. Returns slope in *a and intercept in b, *return
  96    sigmas in *da and *db respectively. Returns normalized *quality of
  97    fit in *chi2 and correlation of fit with data in Rfit. chi2, Rfit,
  98    da and db may be NULL. */
  99 GMX_LIBGMX_EXPORT
 100 int gmx_stats_get_ab(gmx_stats_t stats, int weight,
 101                      real *a, real *b,
 102                      real *da, real *db, real *chi2, real *Rfit);
 103
 104 /* Fit the data to y = ax, possibly weighted, if uncertainties have
 105    been input. Returns slope in *a, sigma in a in *da, and normalized
 106    quality of fit in *chi2 and correlation of fit with data in
 107    Rfit. chi2, Rfit and da may be NULL. */
 108 int gmx_stats_get_a(gmx_stats_t stats, int weight,
 109                     real *a, real *da, real *chi2, real *Rfit);
 110
 111 /* Return the correlation coefficient between the data (x and y) as
 112    input to the structure. */
 113 int gmx_stats_get_corr_coeff(gmx_stats_t stats, real *R);
 114
 115 /* Returns the root mean square deviation between x and y values. */
 116 int gmx_stats_get_rmsd(gmx_stats_t gstats, real *rmsd);
 117
 118 GMX_LIBGMX_EXPORT
 119 int gmx_stats_get_npoints(gmx_stats_t stats, int *N);
 120
 121 GMX_LIBGMX_EXPORT
 122 int gmx_stats_get_average(gmx_stats_t stats, real *aver);
 123
 124 int gmx_stats_get_sigma(gmx_stats_t stats, real *sigma);
 125
 126 int gmx_stats_get_error(gmx_stats_t stats, real *error);
 127
 128 /* Get all three of the above. Pointers may be null, in which case no
 129    assignment will be done. */
 130 GMX_LIBGMX_EXPORT
 131 int gmx_stats_get_ase(gmx_stats_t gstats, real *aver, real *sigma, real *error);
 132
 133 /* Dump the x, y, dx, dy data to a text file */
 134 int gmx_stats_dump_xy(gmx_stats_t gstats, FILE *fp);
 135
 136 /* Make a histogram of the data present. Uses either bindwith to
 137    determine the number of bins, or nbins to determine the binwidth,
 138    therefore one of these should be zero, but not the other. If *nbins = 0
 139    the number of bins will be returned in this variable. ehisto should be one of
 140    ehistoX or ehistoY. If
 141    normalized not equal to zero, the integral of the histogram will be
 142    normalized to one. The output is in two arrays, *x and *y, to which
 143    you should pass a pointer. Memory for the arrays will be allocated
 144    as needed. Function returns one of the estats codes. */
 145 int gmx_stats_make_histogram(gmx_stats_t gstats, real binwidth, int *nbins,
 146                              int ehisto,
 147                              int normalized, real **x, real **y);
 148
 149 /* Return message belonging to error code */
 150 GMX_LIBGMX_EXPORT
 151 const char *gmx_stats_message(int estats);
 152
 153 /****************************************************
 154  * Some statistics utilities for convenience: useful when a complete data
 155  * set is available already from another source, e.g. an xvg file.
 156  ****************************************************/
 157 int lsq_y_ax(int n, real x[], real y[], real *a);
 158 /* Fit a straight line y=ax thru the n data points x, y, return the
 159    slope in *a. Return value can be estatsOK, or something else. */
 160
 161 GMX_LIBGMX_EXPORT
 162 int lsq_y_ax_b(int n, real x[], real y[], real *a, real *b, real *r,
 163                real *chi2);
 164 /* Fit a straight line y=ax+b thru the n data points x,y.
 165  * Returns the "fit quality" sigma = sqrt(chi^2/(n-2)).
 166  * The correlation coefficient is returned in r.
 167  */
 168
 169 int lsq_y_ax_b_xdouble(int n, double x[], real y[],
 170                        real *a, real *b, real *r, real *chi2);
 171 /* As lsq_y_ax_b, but with x in double precision.
 172  */
 173
 174 GMX_LIBGMX_EXPORT
 175 int lsq_y_ax_b_error(int n, real x[], real y[], real dy[],
 176                      real *a, real *b, real *da, real *db,
 177                      real *r, real *chi2);
 178 /* Fit a straight line y=ax+b thru the n data points x,y, with sigma dy
 179  * Returns the "fit quality" sigma = sqrt(chi^2/(n-2)).
 180  * The correlation coefficient is returned in r.
 181  */
 182
 183 #ifdef __cplusplus
 184 }
 185 #endif
 186
 187 #endif