include/gmx_statistics.h

   1 /*
   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  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  10  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  11  * Copyright (c) 2001-2008, The GROMACS development team,
  12  * check out http://www.gromacs.org for more information.
  13
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * If you want to redistribute modifications, please consider that
  20  * scientific software is very special. Version control is crucial -
  21  * bugs must be traceable. We will be happy to consider code for
  22  * inclusion in the official distribution, but derived work must not
  23  * be called official GROMACS. Details are found in the README & COPYING
  24  * files - if they are missing, get the official version at www.gromacs.org.
  25  *
  26  * To help us fund GROMACS development, we humbly ask that you cite
  27  * the papers on the package - you can find them in the top README file.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  *
  31  * And Hey:
  32  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  33  */
  34
  35 #ifndef _GMX_STATS_H
  36 #define _GMX_STATS_H
  37
  38 #ifdef __cplusplus
  39 extern "C" {
  40 #endif
  41
  42 #include "typedefs.h"
  43
  44 typedef struct gmx_stats *gmx_stats_t;
  45
  46 /* Error codes returned by the routines */
  47 enum { estatsOK, estatsNO_POINTS, estatsNO_MEMORY, estatsERROR,
  48        estatsINVALID_INPUT, estatsNOT_IMPLEMENTED, estatsNR };
  49
  50 enum { elsqWEIGHT_NONE, elsqWEIGHT_X, elsqWEIGHT_Y,
  51        elsqWEIGHT_XY, elsqWEIGHT_NR  };
  52
  53 enum { ehistoX, ehistoY, ehistoNR };
  54
  55 gmx_stats_t gmx_stats_init();
  56
  57 int gmx_stats_done(gmx_stats_t stats);
  58
  59 /* Remove outliers from a straight line, where level in units of
  60    sigma. Level needs to be larger than one obviously. */
  61 int gmx_stats_remove_outliers(gmx_stats_t stats,double level);
  62
  63
  64
  65 int gmx_stats_add_point(gmx_stats_t stats,double x,double y,
  66                                double dx,double dy);
  67
  68 /* The arrays dx and dy may be NULL if no uncertainties are available,
  69    in that case zero uncertainties will be assumed. */
  70 int gmx_stats_add_points(gmx_stats_t stats,int n,real *x,real *y,
  71                                 real *dx,real *dy);
  72
  73 /* Return the data points one by one. Return estatsOK while there are
  74    more points, and returns estatsNOPOINTS when the last point has
  75    been returned. Should be used in a while loop. Variables for either
  76    pointer may be NULL, in which case the routine can be used as an
  77    expensive point counter.
  78    If level > 0 then the outliers outside level*sigma are reported
  79    only. */
  80 int gmx_stats_get_point(gmx_stats_t stats,real *x,real *y,
  81                         real *dx,real *dy,real level);
  82
  83 /* Fit the data to y = ax + b, possibly weighted, if uncertainties
  84    have been input. Returns slope in *a and intercept in b, *return
  85    sigmas in *da and *db respectively. Returns normalized *quality of
  86    fit in *chi2 and correlation of fit with data in Rfit. chi2, Rfit,
  87    da and db may be NULL. */
  88 int gmx_stats_get_ab(gmx_stats_t stats,int weight,
  89                             real *a,real *b,
  90                             real *da,real *db,real *chi2,real *Rfit);
  91
  92 /* Fit the data to y = ax, possibly weighted, if uncertainties have
  93    been input. Returns slope in *a, sigma in a in *da, and normalized
  94    quality of fit in *chi2 and correlation of fit with data in
  95    Rfit. chi2, Rfit and da may be NULL. */
  96 int gmx_stats_get_a(gmx_stats_t stats,int weight,
  97                            real *a,real *da,real *chi2,real *Rfit);
  98
  99 /* Return the correlation coefficient between the data (x and y) as
 100    input to the structure. */
 101 int gmx_stats_get_corr_coeff(gmx_stats_t stats,real *R);
 102
 103 /* Returns the root mean square deviation between x and y values. */
 104 int gmx_stats_get_rmsd(gmx_stats_t gstats,real *rmsd);
 105
 106 int gmx_stats_get_npoints(gmx_stats_t stats,int *N);
 107
 108 int gmx_stats_get_average(gmx_stats_t stats,real *aver);
 109
 110 int gmx_stats_get_sigma(gmx_stats_t stats,real *sigma);
 111
 112 int gmx_stats_get_error(gmx_stats_t stats,real *error);
 113
 114 /* Get all three of the above. Pointers may be null, in which case no
 115    assignment will be done. */
 116 int gmx_stats_get_ase(gmx_stats_t gstats,real *aver,real *sigma,real *error);
 117
 118 /* Dump the x, y, dx, dy data to a text file */
 119 int gmx_stats_dump_xy(gmx_stats_t gstats,FILE *fp);
 120
 121 /* Make a histogram of the data present. Uses either bindwith to
 122    determine the number of bins, or nbins to determine the binwidth,
 123    therefore one of these should be zero, but not the other. If *nbins = 0
 124    the number of bins will be returned in this variable. ehisto should be one of
 125    ehistoX or ehistoY. If
 126    normalized not equal to zero, the integral of the histogram will be
 127    normalized to one. The output is in two arrays, *x and *y, to which
 128    you should pass a pointer. Memory for the arrays will be allocated
 129    as needed. Function returns one of the estats codes. */
 130 int gmx_stats_make_histogram(gmx_stats_t gstats,real binwidth,int *nbins,
 131                                     int ehisto,
 132                                     int normalized,real **x,real **y);
 133
 134 /* Return message belonging to error code */
 135 const char *gmx_stats_message(int estats);
 136
 137 /****************************************************
 138  * Some statistics utilities for convenience: useful when a complete data
 139  * set is available already from another source, e.g. an xvg file.
 140  ****************************************************/
 141 int lsq_y_ax(int n, real x[], real y[], real *a);
 142 /* Fit a straight line y=ax thru the n data points x, y, return the
 143    slope in *a. Return value can be estatsOK, or something else. */
 144
 145 int lsq_y_ax_b(int n, real x[], real y[], real *a, real *b,real *r,
 146                       real *chi2);
 147 /* Fit a straight line y=ax+b thru the n data points x,y.
 148  * Returns the "fit quality" sigma = sqrt(chi^2/(n-2)).
 149  * The correlation coefficient is returned in r.
 150  */
 151
 152 int lsq_y_ax_b_xdouble(int n, double x[], real y[],
 153                                real *a, real *b,real *r,real *chi2);
 154 /* As lsq_y_ax_b, but with x in double precision.
 155  */
 156
 157 int lsq_y_ax_b_error(int n, real x[], real y[], real dy[],
 158                              real *a, real *b, real *da, real *db,
 159                              real *r,real *chi2);
 160 /* Fit a straight line y=ax+b thru the n data points x,y, with sigma dy
 161  * Returns the "fit quality" sigma = sqrt(chi^2/(n-2)).
 162  * The correlation coefficient is returned in r.
 163  */
 164
 165 #ifdef __cplusplus
 166 }
 167 #endif
 168
 169 #endif