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