*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2012,2014, by the GROMACS development team, led by
+ * Copyright (c) 2012,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "statistics.h"
-#include <math.h>
+#include <cmath>
#include "gromacs/math/vec.h"
#include "gromacs/utility/real.h"
static int gmx_dnint(double x)
{
- return (int) (x+0.5);
+ return static_cast<int>(x+0.5);
}
typedef struct gmx_stats {
int gmx_stats_add_point(gmx_stats_t gstats, double x, double y,
double dx, double dy)
{
- gmx_stats *stats = (gmx_stats *) gstats;
- int i;
+ gmx_stats *stats = gstats;
if (stats->np+1 >= stats->nalloc)
{
srenew(stats->y, stats->nalloc);
srenew(stats->dx, stats->nalloc);
srenew(stats->dy, stats->nalloc);
- for (i = stats->np; (i < stats->nalloc); i++)
+ for (int i = stats->np; (i < stats->nalloc); i++)
{
stats->x[i] = 0;
stats->y[i] = 0;
int gmx_stats_get_point(gmx_stats_t gstats, real *x, real *y,
real *dx, real *dy, real level)
{
- gmx_stats *stats = (gmx_stats *) gstats;
+ gmx_stats *stats = gstats;
int ok, outlier;
real rmsd, r;
outlier = 0;
while ((outlier == 0) && (stats->np_c < stats->np))
{
- r = fabs(stats->x[stats->np_c] - stats->y[stats->np_c]);
+ r = std::abs(stats->x[stats->np_c] - stats->y[stats->np_c]);
outlier = (r > rmsd*level);
if (outlier)
{
int gmx_stats_add_points(gmx_stats_t gstats, int n, real *x, real *y,
real *dx, real *dy)
{
- int i, ok;
-
- for (i = 0; (i < n); i++)
+ for (int i = 0; (i < n); i++)
{
+ int ok;
if ((ok = gmx_stats_add_point(gstats, x[i], y[i],
(NULL != dx) ? dx[i] : 0,
(NULL != dy) ? dy[i] : 0)) != estatsOK)
double yy, yx, xx, sx, sy, dy, chi2, chi2aa, d2;
double ssxx, ssyy, ssxy;
double w, wtot, yx_nw, sy_nw, sx_nw, yy_nw, xx_nw, dx2, dy2;
- int i, N;
- N = stats->np;
+ int N = stats->np;
+
if (stats->computed == 0)
{
if (N < 1)
sy = sy_nw = 0;
wtot = 0;
d2 = 0;
- for (i = 0; (i < N); i++)
+ for (int i = 0; (i < N); i++)
{
d2 += dsqr(stats->x[i]-stats->y[i]);
if ((stats->dy[i]) && (weight == elsqWEIGHT_Y))
/* Compute average, sigma and error */
stats->aver = sy_nw/N;
- stats->sigma_aver = sqrt(yy_nw/N - dsqr(sy_nw/N));
- stats->error = stats->sigma_aver/sqrt(N);
+ stats->sigma_aver = std::sqrt(yy_nw/N - dsqr(sy_nw/N));
+ stats->error = stats->sigma_aver/std::sqrt(static_cast<double>(N));
/* Compute RMSD between x and y */
- stats->rmsd = sqrt(d2/N);
+ stats->rmsd = std::sqrt(d2/N);
/* Correlation coefficient for data */
yx_nw /= N;
ssxx = N*(xx_nw - dsqr(sx_nw));
ssyy = N*(yy_nw - dsqr(sy_nw));
ssxy = N*(yx_nw - (sx_nw*sy_nw));
- stats->Rdata = sqrt(dsqr(ssxy)/(ssxx*ssyy));
+ stats->Rdata = std::sqrt(dsqr(ssxy)/(ssxx*ssyy));
/* Compute straight line through datapoints, either with intercept
zero (result in aa) or with intercept variable (results in a
chi2aa which returns the deviation from a line y = ax. */
chi2 = 0;
chi2aa = 0;
- for (i = 0; (i < N); i++)
+ for (int i = 0; (i < N); i++)
{
if (stats->dy[i] > 0)
{
}
if (N > 2)
{
- stats->chi2 = sqrt(chi2/(N-2));
- stats->chi2aa = sqrt(chi2aa/(N-2));
+ stats->chi2 = std::sqrt(chi2/(N-2));
+ stats->chi2aa = std::sqrt(chi2aa/(N-2));
/* Look up equations! */
dx2 = (xx-sx*sx);
dy2 = (yy-sy*sy);
- stats->sigma_a = sqrt(stats->chi2/((N-2)*dx2));
- stats->sigma_b = stats->sigma_a*sqrt(xx);
- stats->Rfit = fabs(ssxy)/sqrt(ssxx*ssyy);
- /*stats->a*sqrt(dx2/dy2);*/
- stats->Rfitaa = stats->aa*sqrt(dx2/dy2);
+ stats->sigma_a = std::sqrt(stats->chi2/((N-2)*dx2));
+ stats->sigma_b = stats->sigma_a*std::sqrt(xx);
+ stats->Rfit = std::abs(ssxy)/std::sqrt(ssxx*ssyy);
+ stats->Rfitaa = stats->aa*std::sqrt(dx2/dy2);
}
else
{
real *a, real *b, real *da, real *db,
real *chi2, real *Rfit)
{
- gmx_stats *stats = (gmx_stats *) gstats;
+ gmx_stats *stats = gstats;
int ok;
if ((ok = gmx_stats_compute(stats, weight)) != estatsOK)
int gmx_stats_get_a(gmx_stats_t gstats, int weight, real *a, real *da,
real *chi2, real *Rfit)
{
- gmx_stats *stats = (gmx_stats *) gstats;
+ gmx_stats *stats = gstats;
int ok;
if ((ok = gmx_stats_compute(stats, weight)) != estatsOK)
int gmx_stats_get_average(gmx_stats_t gstats, real *aver)
{
- gmx_stats *stats = (gmx_stats *) gstats;
+ gmx_stats *stats = gstats;
int ok;
if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != estatsOK)
int gmx_stats_get_ase(gmx_stats_t gstats, real *aver, real *sigma, real *error)
{
- gmx_stats *stats = (gmx_stats *) gstats;
+ gmx_stats *stats = gstats;
int ok;
if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != estatsOK)
int gmx_stats_get_sigma(gmx_stats_t gstats, real *sigma)
{
- gmx_stats *stats = (gmx_stats *) gstats;
+ gmx_stats *stats = gstats;
int ok;
if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != estatsOK)
int gmx_stats_get_error(gmx_stats_t gstats, real *error)
{
- gmx_stats *stats = (gmx_stats *) gstats;
+ gmx_stats *stats = gstats;
int ok;
if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != estatsOK)
int gmx_stats_get_corr_coeff(gmx_stats_t gstats, real *R)
{
- gmx_stats *stats = (gmx_stats *) gstats;
+ gmx_stats *stats = gstats;
int ok;
if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != estatsOK)
int gmx_stats_get_rmsd(gmx_stats_t gstats, real *rmsd)
{
- gmx_stats *stats = (gmx_stats *) gstats;
+ gmx_stats *stats = gstats;
int ok;
if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != estatsOK)
int gmx_stats_dump_xy(gmx_stats_t gstats, FILE *fp)
{
- gmx_stats *stats = (gmx_stats *) gstats;
- int i, ok;
+ gmx_stats *stats = gstats;
- for (i = 0; (i < stats->np); i++)
+ for (int i = 0; (i < stats->np); i++)
{
fprintf(fp, "%12g %12g %12g %12g\n", stats->x[i], stats->y[i],
stats->dx[i], stats->dy[i]);
int gmx_stats_remove_outliers(gmx_stats_t gstats, double level)
{
- gmx_stats *stats = (gmx_stats *) gstats;
- int i, iter = 1, done = 0, ok;
+ gmx_stats *stats = gstats;
+ int iter = 1, done = 0, ok;
real rmsd, r;
while ((stats->np >= 10) && !done)
return ok;
}
done = 1;
- for (i = 0; (i < stats->np); )
+ for (int i = 0; (i < stats->np); )
{
- r = fabs(stats->x[i]-stats->y[i]);
+ r = std::abs(stats->x[i]-stats->y[i]);
if (r > level*rmsd)
{
fprintf(stderr, "Removing outlier, iter = %d, rmsd = %g, x = %g, y = %g\n",
int gmx_stats_make_histogram(gmx_stats_t gstats, real binwidth, int *nb,
int ehisto, int normalized, real **x, real **y)
{
- gmx_stats *stats = (gmx_stats *) gstats;
- int i, ok, index = 0, nbins = *nb, *nindex;
+ gmx_stats *stats = gstats;
+ int index = 0, nbins = *nb, *nindex;
double minx, maxx, maxy, miny, delta, dd, minh;
if (((binwidth <= 0) && (nbins <= 0)) ||
}
minx = maxx = stats->x[0];
miny = maxy = stats->y[0];
- for (i = 1; (i < stats->np); i++)
+ for (int i = 1; (i < stats->np); i++)
{
miny = (stats->y[i] < miny) ? stats->y[i] : miny;
maxy = (stats->y[i] > maxy) ? stats->y[i] : maxy;
}
snew(*x, nbins);
snew(nindex, nbins);
- for (i = 0; (i < nbins); i++)
+ for (int i = 0; (i < nbins); i++)
{
(*x)[i] = minh + binwidth*(i+0.5);
}
}
snew(*y, nbins);
- for (i = 0; (i < stats->np); i++)
+ for (int i = 0; (i < stats->np); i++)
{
if (ehisto == ehistoY)
{
- index = (stats->y[i]-miny)/binwidth;
+ index = static_cast<int>((stats->y[i]-miny)/binwidth);
}
else if (ehisto == ehistoX)
{
- index = (stats->x[i]-minx)/binwidth;
+ index = static_cast<int>((stats->x[i]-minx)/binwidth);
}
if (index < 0)
{
{
*nb = nbins;
}
- for (i = 0; (i < nbins); i++)
+ for (int i = 0; (i < nbins); i++)
{
if (nindex[i] > 0)
{
return ok;
}
- /* int i;
- double xx,yx;
-
- yx=xx=0.0;
- for (i=0; i<n; i++) {
- yx+=y[i]*x[i];
- xx+=x[i]*x[i];
- }
- * a=yx/xx;
- */
return estatsOK;
}
static int low_lsq_y_ax_b(int n, real *xr, double *xd, real yr[],
real *a, real *b, real *r, real *chi2)
{
- int i, ok;
- gmx_stats_t lsq;
+ gmx_stats_t lsq = gmx_stats_init();
+ int ok;
- lsq = gmx_stats_init();
- for (i = 0; (i < n); i++)
+ for (int i = 0; (i < n); i++)
{
- if ((ok = gmx_stats_add_point(lsq, (NULL != xd) ? xd[i] : xr[i], yr[i], 0, 0))
- != estatsOK)
+ double pt;
+
+ if (xd != NULL)
+ {
+ pt = xd[i];
+ }
+ else if (xr != NULL)
+ {
+ pt = xr[i];
+ }
+ else
+ {
+ gmx_incons("Either xd or xr has to be non-NULL in low_lsq_y_ax_b()");
+ }
+
+ if ((ok = gmx_stats_add_point(lsq, pt, yr[i], 0, 0)) != estatsOK)
{
return ok;
}
}
return estatsOK;
- /*
- double x,y,yx,xx,yy,sx,sy,chi2;
-
- yx=xx=yy=sx=sy=0.0;
- for (i=0; i<n; i++) {
- if (xd != NULL) {
- x = xd[i];
- } else {
- x = xr[i];
- }
- y = yr[i];
-
- yx += y*x;
- xx += x*x;
- yy += y*y;
- sx += x;
- sy += y;
- }
- * a = (n*yx-sy*sx)/(n*xx-sx*sx);
- * b = (sy-(*a)*sx)/n;
- * r = sqrt((xx-sx*sx)/(yy-sy*sy));
-
- chi2 = 0;
- if (xd != NULL) {
- for(i=0; i<n; i++)
- chi2 += dsqr(yr[i] - ((*a)*xd[i] + (*b)));
- } else {
- for(i=0; i<n; i++)
- chi2 += dsqr(yr[i] - ((*a)*xr[i] + (*b)));
- }
-
- if (n > 2)
- return sqrt(chi2/(n-2));
- else
- return 0;
- */
}
int lsq_y_ax_b(int n, real x[], real y[], real *a, real *b, real *r, real *chi2)
real *a, real *b, real *da, real *db,
real *r, real *chi2)
{
- gmx_stats_t lsq;
- int i, ok;
+ gmx_stats_t lsq = gmx_stats_init();
+ int ok;
- lsq = gmx_stats_init();
- for (i = 0; (i < n); i++)
+ for (int i = 0; (i < n); i++)
{
if ((ok = gmx_stats_add_point(lsq, x[i], y[i], 0, dy[i])) != estatsOK)
{
sfree(lsq);
return estatsOK;
- /*
- double sxy,sxx,syy,sx,sy,w,s_2,dx2,dy2,mins;
-
- sxy=sxx=syy=sx=sy=w=0.0;
- mins = dy[0];
- for(i=1; (i<n); i++)
- mins = min(mins,dy[i]);
- if (mins <= 0)
- gmx_fatal(FARGS,"Zero or negative weigths in linear regression analysis");
-
- for (i=0; i<n; i++) {
- s_2 = dsqr(1.0/dy[i]);
- sxx += s_2*dsqr(x[i]);
- sxy += s_2*y[i]*x[i];
- syy += s_2*dsqr(y[i]);
- sx += s_2*x[i];
- sy += s_2*y[i];
- w += s_2;
- }
- sxx = sxx/w;
- sxy = sxy/w;
- syy = syy/w;
- sx = sx/w;
- sy = sy/w;
- dx2 = (sxx-sx*sx);
- dy2 = (syy-sy*sy);
- * a=(sxy-sy*sx)/dx2;
- * b=(sy-(*a)*sx);
-
- * chi2=0;
- for(i=0; i<n; i++)
- * chi2+=dsqr((y[i]-((*a)*x[i]+(*b)))/dy[i]);
- * chi2 = *chi2/w;
-
- * da = sqrt(*chi2/((n-2)*dx2));
- * db = *da*sqrt(sxx);
- * r = *a*sqrt(dx2/dy2);
-
- if (debug)
- fprintf(debug,"sx = %g, sy = %g, sxy = %g, sxx = %g, w = %g\n"
- "chi2 = %g, dx2 = %g\n",
- sx,sy,sxy,sxx,w,*chi2,dx2);
-
- if (n > 2)
- * chi2 = sqrt(*chi2/(n-2));
- else
- * chi2 = 0;
- */
}
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2008, The GROMACS development team.
- * Copyright (c) 2011,2014, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#include "gmxpre.h"
-
-#include <stdio.h>
-
-#include "gromacs/math/vec.h"
-#include "gromacs/random/random.h"
-#include "gromacs/statistics/statistics.h"
-#include "gromacs/utility/real.h"
-#include "gromacs/utility/smalloc.h"
-
-static void horizontal()
-{
- gmx_rng_t rng;
- gmx_stats_t straight;
- int i, ok, n = 1000;
- real y, a, b, da, db, aver, sigma, error, chi2, R, *xh, *yh;
- FILE *fp;
-
- rng = gmx_rng_init(13);
- straight = gmx_stats_init();
- for (i = 0; (i < n); i++)
- {
- y = gmx_rng_uniform_real(rng);
- if ((ok = gmx_stats_add_point(straight, i, y, 0, 0)) != estatsOK)
- {
- fprintf(stderr, "%s\n", gmx_stats_message(ok));
- }
- }
- /* Horizontal test */
- if ((ok = gmx_stats_get_ase(straight, &aver, &sigma, &error)) != estatsOK)
- {
- fprintf(stderr, "%s\n", gmx_stats_message(ok));
- }
- fp = fopen("straight.xvg", "w");
- if ((ok = gmx_stats_dump_xy(straight, fp)) != estatsOK)
- {
- fprintf(stderr, "%s\n", gmx_stats_message(ok));
- }
- fclose(fp);
- printf("Horizontal line: average %g, sigma %g, error %g\n", aver, sigma, error);
- if ((ok = gmx_stats_done(straight)) != estatsOK)
- {
- fprintf(stderr, "%s\n", gmx_stats_message(ok));
- }
-}
-
-static void line()
-{
- gmx_rng_t rng;
- gmx_stats_t line;
- int i, dy, ok, n = 1000;
- real y, a, b, da, db, aver, sigma, error, chi2, R, rfit;
- const real a0 = 0.23, b0 = 2.7;
- FILE *fp;
-
- for (dy = 0; (dy < 2); dy++)
- {
- rng = gmx_rng_init(13);
- line = gmx_stats_init();
- for (i = 0; (i < n); i++)
- {
- y = a0*i+b0+50*(gmx_rng_uniform_real(rng)-0.5);
- if ((ok = gmx_stats_add_point(line, i, y, 0, dy*0.1)) != estatsOK)
- {
- fprintf(stderr, "%s\n", gmx_stats_message(ok));
- }
- }
- /* Line with slope test */
- if ((ok = gmx_stats_get_ab(line, elsqWEIGHT_NONE, &a, &b, &da, &db, &chi2, &rfit)) != estatsOK)
- {
- fprintf(stderr, "%s\n", gmx_stats_message(ok));
- }
- if ((ok = gmx_stats_get_corr_coeff(line, &R)) != estatsOK)
- {
- fprintf(stderr, "%s\n", gmx_stats_message(ok));
- }
- if (dy == 0)
- {
- fp = fopen("line0.xvg", "w");
- }
- else
- {
- fp = fopen("line1.xvg", "w");
- }
- if ((ok = gmx_stats_dump_xy(line, fp)) != estatsOK)
- {
- fprintf(stderr, "%s\n", gmx_stats_message(ok));
- }
- fclose(fp);
- printf("Line with eqn. y = %gx + %g with noise%s\n", a0, b0,
- (dy == 0) ? "" : " and uncertainties");
- printf("Found: a = %g +/- %g, b = %g +/- %g\n", a, da, b, db);
- if ((ok = gmx_stats_done(line)) != estatsOK)
- {
- fprintf(stderr, "%s\n", gmx_stats_message(ok));
- }
- gmx_rng_destroy(rng);
- }
-}
-
-static void histogram()
-{
- gmx_rng_t rng;
- gmx_stats_t camel;
- int i, ok, n = 1000, norm;
- real y, a, b, da, db, aver, sigma, error, chi2, R, *xh, *yh;
- const real a0 = 0.23, b0 = 2.7;
- FILE *fp;
- char fn[256];
-
- for (norm = 0; (norm < 2); norm++)
- {
- rng = gmx_rng_init(13);
- camel = gmx_stats_init();
- for (i = 0; (i < n); i++)
- {
- y = sqr(gmx_rng_uniform_real(rng));
- if ((ok = gmx_stats_add_point(camel, i, y+1, 0, 0)) != estatsOK)
- {
- fprintf(stderr, "%s\n", gmx_stats_message(ok));
- }
- y = sqr(gmx_rng_uniform_real(rng));
- if ((ok = gmx_stats_add_point(camel, i+0.5, y+2, 0, 0)) != estatsOK)
- {
- fprintf(stderr, "%s\n", gmx_stats_message(ok));
- }
- }
- /* Histogram test */
- if ((ok = gmx_stats_make_histogram(camel, 0, 101, norm, &xh, &yh)) != estatsOK)
- {
- fprintf(stderr, "%s\n", gmx_stats_message(ok));
- }
- sprintf(fn, "histo%d-data.xvg", norm);
- fp = fopen(fn, "w");
- gmx_stats_dump_xy(camel, fp);
- fclose(fp);
- sprintf(fn, "histo%d.xvg", norm);
- fp = fopen(fn, "w");
- for (i = 0; (i < 101); i++)
- {
- fprintf(fp, "%12g %12g\n", xh[i], yh[i]);
- }
- fclose(fp);
- sfree(xh);
- sfree(yh);
- }
-}
-
-int main(int argc, char *argv[])
-{
- line();
- horizontal();
- histogram();
-
- return 0;
-}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff