From: David van der Spoel Date: Tue, 14 Jan 2014 18:17:13 +0000 (+0100) Subject: This patch moves the statistics library to its own dir. X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=630187a2143ed6a19e98ba4c7cdfd74496191e96;p=alexxy%2Fgromacs.git This patch moves the statistics library to its own dir. In the move to depopulate the legacyheaders directory http://redmine.gromacs.org/issues/1415 this is a small contribution. Histogram files removed. Added doxygen style comments as well. Change-Id: I9da9fefac0d4dd1c5de13862f348adc5c812ad4d --- diff --git a/src/gromacs/CMakeLists.txt b/src/gromacs/CMakeLists.txt index 378af42d21..64a948161a 100644 --- a/src/gromacs/CMakeLists.txt +++ b/src/gromacs/CMakeLists.txt @@ -62,6 +62,7 @@ add_subdirectory(pulling) if (NOT GMX_BUILD_MDRUN_ONLY) add_subdirectory(legacyheaders) add_subdirectory(gmxana) + add_subdirectory(statistics) add_subdirectory(analysisdata) add_subdirectory(selection) add_subdirectory(trajectoryanalysis) diff --git a/src/gromacs/gmxana/gmx_anadock.c b/src/gromacs/gmxana/gmx_anadock.c index 361853dd3d..58331c274b 100644 --- a/src/gromacs/gmxana/gmx_anadock.c +++ b/src/gromacs/gmxana/gmx_anadock.c @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2012,2013, by the GROMACS development team, led by + * Copyright (c) 2012,2013,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. @@ -45,7 +45,7 @@ #include "smalloc.h" #include "string2.h" #include "vec.h" -#include "gmx_statistics.h" +#include "gromacs/statistics/statistics.h" #include "gromacs/commandline/pargs.h" #include "typedefs.h" #include "xvgr.h" diff --git a/src/gromacs/gmxana/gmx_analyze.c b/src/gromacs/gmxana/gmx_analyze.c index 503a7d3bae..6551ef98c4 100644 --- a/src/gromacs/gmxana/gmx_analyze.c +++ b/src/gromacs/gmxana/gmx_analyze.c @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013, by the GROMACS development team, led by + * Copyright (c) 2013,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. @@ -51,7 +51,7 @@ #include "readinp.h" #include "txtdump.h" #include "gstat.h" -#include "gmx_statistics.h" +#include "gromacs/statistics/statistics.h" #include "xvgr.h" #include "gmx_ana.h" #include "geminate.h" diff --git a/src/gromacs/gmxana/gmx_current.c b/src/gromacs/gmxana/gmx_current.c index cf84d9d0b2..66825b4a3d 100644 --- a/src/gromacs/gmxana/gmx_current.c +++ b/src/gromacs/gmxana/gmx_current.c @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2008,2009,2010,2011,2012,2013, by the GROMACS development team, led by + * Copyright (c) 2008,2009,2010,2011,2012,2013,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. @@ -47,7 +47,7 @@ #include "pbc.h" #include "physics.h" #include "index.h" -#include "gmx_statistics.h" +#include "gromacs/statistics/statistics.h" #include "gmx_ana.h" #include "macros.h" diff --git a/src/gromacs/gmxana/gmx_dipoles.cpp b/src/gromacs/gmxana/gmx_dipoles.cpp index 61269ffbf6..2d1094ff20 100644 --- a/src/gromacs/gmxana/gmx_dipoles.cpp +++ b/src/gromacs/gmxana/gmx_dipoles.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013, by the GROMACS development team, led by + * Copyright (c) 2013,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. @@ -52,7 +52,7 @@ #include "gromacs/fileio/futil.h" #include "xvgr.h" #include "txtdump.h" -#include "gmx_statistics.h" +#include "gromacs/statistics/statistics.h" #include "gstat.h" #include "index.h" #include "random.h" diff --git a/src/gromacs/gmxana/gmx_msd.c b/src/gromacs/gmxana/gmx_msd.c index 94075a1134..bbb6c6c36a 100644 --- a/src/gromacs/gmxana/gmx_msd.c +++ b/src/gromacs/gmxana/gmx_msd.c @@ -52,7 +52,7 @@ #include "typedefs.h" #include "xvgr.h" #include "gstat.h" -#include "gmx_statistics.h" +#include "gromacs/statistics/statistics.h" #include "gromacs/fileio/tpxio.h" #include "gromacs/fileio/trxio.h" #include "pbc.h" diff --git a/src/gromacs/gmxlib/CMakeLists.txt b/src/gromacs/gmxlib/CMakeLists.txt index 5b6a510af4..b07fa65a9a 100644 --- a/src/gromacs/gmxlib/CMakeLists.txt +++ b/src/gromacs/gmxlib/CMakeLists.txt @@ -1,7 +1,7 @@ # # This file is part of the GROMACS molecular simulation package. # -# Copyright (c) 2009,2010,2012, by the GROMACS development team, led by +# Copyright (c) 2009,2010,2012,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. @@ -38,7 +38,7 @@ add_subdirectory(nonbonded) # The nonbonded directory contains subdirectories that are only # conditionally built, so we cannot use a GLOB_RECURSE here. -file(GLOB GMXLIB_SOURCES *.c *.cpp statistics/*.c) +file(GLOB GMXLIB_SOURCES *.c *.cpp) # This would be the standard way to include thread_mpi, but we want libgmx # to link the functions directly diff --git a/src/gromacs/gmxlib/statistics/histogram.c b/src/gromacs/gmxlib/statistics/histogram.c deleted file mode 100644 index 6fdfad79c5..0000000000 --- a/src/gromacs/gmxlib/statistics/histogram.c +++ /dev/null @@ -1,1179 +0,0 @@ -/* - * This file is part of the GROMACS molecular simulation package. - * - * Copyright (c) 2009,2010,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. - */ -/*! \page histograms Histogram calculation - * - * Functions to calculate histograms are defined in histogram.h. - * The principle is simple: histogram calculation is set up with - * gmx_histogram_create() and gmx_histogram_set_*(), the histogram is sampled - * using the provided functions, and gmx_histogram_finish() is called to - * finish the sampling. - * Various post-processing functions can then be used to normalize the - * histogram in various ways and to write it into a file. - * gmx_histogram_get_*() can be used to access data from the histogram. - * gmx_histogram_free() can be used to free the memory allocated for a - * histogram. - * - * - * \section histogram_sampling Initialization and sampling - * - * A histogram calculation is initialized by calling gmx_histogram_create() - * or gmx_histogram_create_range(). - * The first function takes the type of the histogram (see \ref e_histogram_t) - * and the number of bins, and allocates the necessary memory. - * The bin locations can be initialized with gmx_histogram_set_integerbins() - * and one of gmx_histogram_set_binwidth() and gmx_histogram_set_range(). - * Only uniformly spaced bins are currently supported. - * The second initialization function takes the beginning and end of the range - * and the bin width. - * The treatment of values that fall outside the histogram range can be - * set with gmx_histogram_set_all(). - * gmx_histogram_set_blocksize() can be used to specify a block size for error - * estimates. If not called, no error estimates are calculated. - * Otherwise, a histogram is calculated for each block of given number of - * frames, and the error estimated as the error of the mean of these block - * histograms. - * If the number of frames is not divisible by the block size, - * the last block is discarded. - * gmx_histogram_set_block_output() can be used to write the individual - * block histograms to a file (it is currently not very flexible). - * - * After initialization, three sets of sampling functions are provided - * to sample the differen types of histograms: - * - Use gmx_histogram_increment() or gmx_histogram_increment_bin() to - * calculate simple histograms (\ref HIST_SIMPLE histograms). - * - Use gmx_histogram_add() or gmx_histogram_add_to_bin() to calculate - * histograms where different points can have different weight - * (\ref HIST_WEIGHT histograms). - * - Use gmx_histogram_add_item() or gmx_histogram_add_item_to_bin() to - * calculate averages within each bin - * (\ref HIST_BINAVER histograms). - * - * The functions aboge that take bin indices directly do no range checking, - * i.e., one must ensure that the bin number is between 0 and the number of - * bins in the histogram. In contrast, the functions that use real values to - * locate the correct bin use range checking, and silently discard values - * outside the histogram (unless gmx_histogram_set_all() has been called). - * gmx_histogram_find_bin() can be used to find the bin number corresponding - * to a value, adhering to the gmx_histogram_set_all() setting. - * - * After a frame is sampled, gmx_histogram_finish_frame() needs to be - * called. - * After all data is sampled, gmx_histogram_finish() should be called. - * After these calls, the histogram will be normalized such that the - * sum over all bins equals the average number of samples per frame, and - * the error field contains the error of the mean at each bin. - * - * - * \section histogram_processing Post-processing - * - * The following functions can be used to process the histograms: - * - gmx_histogram_clone() makes a deep copy. - * - gmx_histogram_resample_dblbw() also makes a deep copy, - * but uses a binwidth that is double of the input one - * (see below for usage). - * - gmx_histogram_normalize_prob() normalizes a histogram such that its - * integral becomes one. - * - gmx_histogram_scale() and gmx_histogram_scale_vec() can be used - * to scale a histogram uniformly or non-uniformly, respectively. - * - * There are also functions to write histograms to a file: - * - gmx_histogram_write() writes a single histogram, optionally with - * errors. - * - gmx_histogram_write_array() writes multiple histograms with identical - * bins. Values and/or errors can be written. - * - gmx_histogram_write_cum_array() writes cumulative histograms for - * multiple histograms with identical bins. - * - * gmx_histogram_resample_dblbw() is useful if a histogram and its - * cumulative sum are required at same points. - * One can then sample the histogram with half the desired binwidth, and - * then use gmx_histogram_resample_dblbw() to construct two different - * versions of the histogram (by changing \p bIntegerBins), one suitable for - * writing out the histogram and the other for writing the cumulative - * distribution. - * - * Access to the histogram data is also provided through the following - * functions: - * - gmx_histogram_get_nbins() - * - gmx_histogram_get_binwidth() - * - gmx_histogram_get_value() - * - gmx_histogram_get_bin_value() - * - gmx_histogram_get_values() - * - gmx_histogram_get_errors() - * - * gmx_histogram_free() can be used to free memory associated with a - * histogram when it is no longer needed. - */ -/*! \internal \file - * \brief Implementation of functions in histogram.h. - */ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include - -#include "smalloc.h" -#include "vec.h" - -#include "histogram.h" - -/*! \internal \brief - * Stores data for a histogram. - */ -struct gmx_histogram_t -{ - /** The left edge of the first bin. */ - real start; - /** The right edge of the last bin. */ - real end; - /** Bin width. */ - real binwidth; - /** Number of bins. */ - int nbins; - /** The final histogram. */ - double *hist; - /** Standard deviation of each bin (calculated from block histograms). */ - double *histerr; - - /** Histogram type. */ - e_histogram_t type; - /** Histogram flags; */ - int flags; - /** Block size for averaging. */ - int bsize; - /** Output file for block histograms (can be NULL). */ - FILE *blockfp; - - /** Inverse bin width. */ - real invbw; - /** Current histogram value. */ - double *chist; - /** Current histogram count. */ - int *cn; - /** Number of frames read for the current block so far. */ - int nframes; - /** Number of blocks averaged so far. */ - int nblocks; -}; - -/*! - * \param[out] hp Histogram to create. - * \param[in] type Histogram type to create; - * \param[in] nbins Number of bins. - * \returns 0 on success, a non-zero error code on error. - * - * Initialized the histogram structure \p h with the provided values, - * allocating the necessary memory. - */ -int -gmx_histogram_create(gmx_histogram_t **hp, e_histogram_t type, int nbins) -{ - gmx_histogram_t *h; - - if (nbins <= 0) - { - *hp = NULL; - gmx_incons("number of histogram bins <= 0"); - return EINVAL; - } - - snew(h, 1); - h->start = 0; - h->end = 0; - h->binwidth = 0; - h->nbins = nbins; - h->hist = NULL; - h->histerr = NULL; - - h->type = type; - h->flags = 0; - - h->bsize = 0; - h->blockfp = NULL; - h->invbw = 0; - h->chist = NULL; - h->cn = NULL; - h->nframes = 0; - h->nblocks = 0; - - snew(h->hist, nbins + 1); - snew(h->histerr, nbins + 1); - snew(h->cn, nbins + 1); - if (type != HIST_SIMPLE) - { - snew(h->chist, nbins + 1); - } - gmx_histogram_clear(h); - - *hp = h; - return 0; -} - -/*! - * \param[out] hp Histogram to create. - * \param[in] type Histogram type to create; - * \param[in] start Left edge/center of the first bin - * (depending on \p bIntegerBins). - * \param[in] end Right edge/center of the last bin - * (depending on \p bIntegerBins). - * \param[in] binw Bin width. - * \param[in] bIntegerBins If true, histogram bins are centered at - * \c start+n*binwidth, otherwise the centers are at - * \c start+(n+0.5)*binwidth. - * \returns 0 on success, a non-zero error code on error. - * - * Initialized the histogram structure \p h with the provided values, - * allocating the necessary memory. - */ -int -gmx_histogram_create_range(gmx_histogram_t **hp, e_histogram_t type, - real start, real end, real binw, gmx_bool bIntegerBins) -{ - gmx_histogram_t *h; - int nbins; - int rc; - - *hp = NULL; - if (start >= end || binw <= 0) - { - gmx_incons("histogram left edge larger than right edge or bin width <= 0"); - return EINVAL; - } - - /* Adjust the end points and calculate the number of bins */ - if (bIntegerBins) - { - nbins = ceil((end - start) / binw) + 1; - end = start + (nbins - 1) * binw; - } - else - { - start = binw * floor(start / binw); - end = binw * ceil(end / binw); - if (start != 0) - { - start -= binw; - } - end += binw; - nbins = (int)((end - start) / binw + 0.5); - } - /* Create the histogram */ - rc = gmx_histogram_create(&h, type, nbins); - if (rc != 0) - { - return rc; - } - /* Set it up */ - gmx_histogram_set_integerbins(h, bIntegerBins); - gmx_histogram_set_binwidth(h, start, binw); - - *hp = h; - return 0; -} - -/*! - * \param[in,out] h Histogram to clear. - * - * Histograms are automatically cleared when initialized; you only need to - * call this function if you want to reuse a histogram structure that has - * already been used for sampling. - */ -void -gmx_histogram_clear(gmx_histogram_t *h) -{ - int i; - - if (h->nbins <= 0) - { - return; - } - for (i = 0; i <= h->nbins; ++i) - { - h->hist[i] = 0; - h->histerr[i] = 0; - if (h->chist) - { - h->chist[i] = 0; - } - h->cn[i] = 0; - } - h->nframes = 0; - h->nblocks = 0; -} - -/*! - * \param[in] h Histogram to free. - * - * The pointer \p h is invalid after the call. - */ -void -gmx_histogram_free(gmx_histogram_t *h) -{ - sfree(h->chist); - sfree(h->cn); - sfree(h->hist); - sfree(h->histerr); - sfree(h); -} - -/*! - * \param[in,out] h Histogram data structure. - * \param[in] start Left edge/center of the first bin - * (depending on gmx_histogram_set_integerbins()). - * \param[in] binw Bin width. - * \returns 0 on success, a non-zero error code on error. - */ -int -gmx_histogram_set_binwidth(gmx_histogram_t *h, real start, real binw) -{ - if (binw <= 0) - { - gmx_incons("histogram binwidth <= 0"); - return EINVAL; - } - if (h->flags & HIST_INTEGERBINS) - { - start -= 0.5*binw; - } - h->start = start; - h->binwidth = binw; - h->end = start + h->nbins * binw; - h->invbw = 1.0/binw; - h->flags |= HIST_INITBW; - return 0; -} - -/*! - * \param[in,out] h Histogram data structure. - * \param[in] start Left edge/center of the first bin - * (depending on gmx_histogram_set_integerbins()). - * \param[in] end Right edge/center of the last bin - * (depending on gmx_histogram_set_integerbins()). - * \returns 0 on success, a non-zero error code on error. - */ -int -gmx_histogram_set_range(gmx_histogram_t *h, real start, real end) -{ - if (start >= end) - { - gmx_incons("histogram left edge larger than right edge"); - return EINVAL; - } - h->start = start; - h->end = end; - if (h->flags & HIST_INTEGERBINS) - { - h->binwidth = (end - start) / (h->nbins - 1); - start -= 0.5*h->binwidth; - end += 0.5*h->binwidth; - } - else - { - h->binwidth = (end - start) / h->nbins; - } - h->invbw = 1.0/h->binwidth; - h->flags &= ~HIST_INITBW; - return 0; -} - -/*! - * \param[in,out] h Histogram data structure. - * \param[in] bIntegerBins If true, histogram bins are centered at - * \c start+n*binwidth, otherwise the centers are at - * \c start+(n+0.5)*binwidth. - */ -void -gmx_histogram_set_integerbins(gmx_histogram_t *h, gmx_bool bIntegerBins) -{ - /* Adjust the ranges if they have been initialized */ - if (h->start < h->end) - { - if (h->flags & HIST_INTEGERBINS) - { - h->start += 0.5*h->binwidth; - if (h->flags & HIST_INITBW) - { - h->end += 0.5*h->binwidth; - } - else - { - h->end -= 0.5*h->binwidth; - } - } - if (bIntegerBins) - { - h->start -= 0.5*h->binwidth; - if (h->flags & HIST_INITBW) - { - h->end -= 0.5*h->binwidth; - } - else - { - h->end += 0.5*h->binwidth; - } - } - } - if (bIntegerBins) - { - h->flags |= HIST_INTEGERBINS; - } - else - { - h->flags &= ~HIST_INTEGERBINS; - } -} - -/*! - * \param[in,out] h Histogram data structure. - * \param[in] bAll If true, values outside the histogram edges are added - * to the bins at the edges. - * - * \p bAll can be used to avoid rounding errors in cases where the histogram - * spans the full range of possible values. If set, the values that are at - * the exact maximum are still correctly included. - */ -void -gmx_histogram_set_all(gmx_histogram_t *h, gmx_bool bAll) -{ - if (bAll) - { - h->flags |= HIST_ALL; - } - else - { - h->flags &= ~HIST_ALL; - } -} - -/*! - * \param[in,out] h Histogram data structure. - * \param[in] bsize Block size for error estimates. - * \returns 0 on success, a non-zero error code on error. - * - * If \p bsize is zero, no block averaging or error estimates are done. - * This is also the case if this function is not called. - */ -int -gmx_histogram_set_blocksize(gmx_histogram_t *h, int bsize) -{ - if (bsize < 0) - { - gmx_incons("histogram block size < 0"); - return EINVAL; - } - h->bsize = bsize; - return 0; -} - -/*! - * \param[in,out] h Histogram data structure. - * \param[in] fp File for block output. - * \returns 0 on success, a non-zero error code on error. - * - * Sets a file into which each block histogram (the histogram after each - * \c gmx_histogram_t::bsize frames) is written. - * All histograms are written to the same file, separated by empty lines. - * If this function is not called, the block histograms are only used for - * error estimation, not written to a file. - */ -int -gmx_histogram_set_block_output(gmx_histogram_t *h, FILE *fp) -{ - if (h->bsize <= 0) - { - gmx_incons("histogram block size not set but output initialized"); - return EINVAL; - } - h->blockfp = fp; - return 0; -} - -/*! - * \param[in] h Histogram data. - * \param[in] pos Position. - * \returns Bin index in \p h corresponding to \p pos, or -1 if there is no - * such bin. - * - * If gmx_histogram_set_all() has been called with TRUE, values outside the - * histogram are mapped to the bins at the edges. - */ -int -gmx_histogram_find_bin(gmx_histogram_t *h, real pos) -{ - if (pos < h->start) - { - if (h->flags & HIST_ALL) - { - return 0; - } - else - { - return -1; - } - } - else if (pos >= h->end) - { - if (h->flags & HIST_ALL) - { - return h->nbins - 1; - } - else - { - return -1; - } - } - /* Calculate the bin index if we are inside the box */ - return (int)((pos - h->start)*h->invbw); -} - -/*! - * \param[in] h Histogram data. - * \returns Number of bins in \p h. - */ -int -gmx_histogram_get_nbins(gmx_histogram_t *h) -{ - return h->nbins; -} - -/*! - * \param[in] h Histogram data. - * \returns Bin width of \p h, 0 if no binwidth has been set. - */ -real -gmx_histogram_get_binwidth(gmx_histogram_t *h) -{ - return h->binwidth; -} - -/*! - * \param[in] h Histogram data. - * \param[in] pos Position. - * \param[out] value Pointer to receive the value (can be NULL). - * \param[out] err Pointer to receive the value (can be NULL). - * - * If \p pos is outside the range of the histogram, zeros are returned. - */ -void -gmx_histogram_get_value(gmx_histogram_t *h, real pos, double *value, double *err) -{ - int bin; - double v, e; - - if (pos < h->start || pos > h->end) - { - v = e = 0; - } - else - { - bin = gmx_histogram_find_bin(h, pos); - if (bin < 0) - { - v = e = 0; - } - else - { - v = h->hist[bin]; - e = h->histerr[bin]; - } - } - if (value) - { - *value = v; - } - if (err) - { - *err = e; - } -} - -/*! - * \param[in] h Histogram data. - * \param[in] bin Bin number. - * \param[out] value Pointer to receive the value (can be NULL). - * \param[out] err Pointer to receive the value (can be NULL). - * - * If \p bin is outside the valid range, zeros are returned. - */ -void -gmx_histogram_get_bin_value(gmx_histogram_t *h, int bin, double *value, double *err) -{ - double v, e; - - if (bin < 0 || bin >= h->nbins) - { - v = e = 0; - } - else - { - v = h->hist[bin]; - e = h->histerr[bin]; - } - if (value) - { - *value = v; - } - if (err) - { - *err = e; - } -} - -/*! - * \param[in] h Histogram data. - * \returns Pointer to an array of values for \p h. - * - * The returned array has one element for each bin of \p h. - * The returned pointer should not be freed. - */ -double * -gmx_histogram_get_values(gmx_histogram_t *h) -{ - return h->hist; -} - -/*! - * \param[in] h Histogram data. - * \returns Pointer to an array of errors for \p h. - * - * The returned array has one element for each bin of \p h. - * The returned pointer should not be freed. - */ -double * -gmx_histogram_get_errors(gmx_histogram_t *h) -{ - return h->histerr; -} - -/*! - * \param[in,out] h Histogram data. - * \param[in] pos Position. - */ -void -gmx_histogram_increment(gmx_histogram_t *h, real pos) -{ - int bin = gmx_histogram_find_bin(h, pos); - if (bin < 0) - { - return; - } - h->cn[bin]++; -} - -/*! - * \param[in,out] h Histogram data. - * \param[in] pos Position. - * \param[in] value Value to add. - */ -void -gmx_histogram_add(gmx_histogram_t *h, real pos, double value) -{ - int bin = gmx_histogram_find_bin(h, pos); - if (bin < 0) - { - return; - } - h->chist[bin] += value; -} - -/*! - * \param[in,out] h Histogram data. - * \param[in] pos Position. - * \param[in] value Value to add. - */ -void -gmx_histogram_add_item(gmx_histogram_t *h, real pos, double value) -{ - int bin = gmx_histogram_find_bin(h, pos); - if (bin < 0) - { - return; - } - h->chist[bin] += value; - h->cn[bin]++; -} - -/*! - * \param[in,out] h Histogram data. - * \param[in] bin Bin number. - * - * No checks for out-of-bound errors are done: \p bin should be >=0 and - * < \p h->nbins. - */ -void -gmx_histogram_increment_bin(gmx_histogram_t *h, int bin) -{ - h->cn[bin]++; -} - -/*! - * \param[in,out] h Histogram data. - * \param[in] bin Bin number. - * \param[in] value Value to add. - * - * No checks for out-of-bound errors are done: \p bin should be >=0 and - * < \p h->nbins. - */ -void -gmx_histogram_add_to_bin(gmx_histogram_t *h, int bin, double value) -{ - h->chist[bin] += value; -} - -/*! - * \param[in,out] h Histogram data. - * \param[in] bin Bin number. - * \param[in] value Value to add. - * - * No checks for out-of-bound errors are done: \p bin should be >=0 and - * < \p h->nbins. - */ -void -gmx_histogram_add_item_to_bin(gmx_histogram_t *h, int bin, double value) -{ - h->chist[bin] += value; - h->cn[bin]++; -} - -/*! \brief - * Processes a sampled block histogram. - * - * \param[in,out] h Histogram data structure. - */ -static void -finish_histogram_block(gmx_histogram_t *h) -{ - int i; - real v; - - if (h->nframes == 0) - { - return; - } - - if (h->flags & HIST_ALL) - { - if (h->chist) - { - h->chist[h->nbins-1] += h->chist[h->nbins]; - } - h->cn[h->nbins-1] += h->cn[h->nbins]; - } - - for (i = 0; i <= h->nbins; ++i) - { - if (h->chist) - { - v = h->chist[i] / (h->cn[i] > 0 ? h->cn[i] : h->nframes); - } - else - { - v = ((real)h->cn[i]) / h->nframes; - } - if (h->blockfp) - { - fprintf(h->blockfp, "%10g %10g\n", h->start + (i+0.5)*h->binwidth, v); - } - h->hist[i] += v; - h->histerr[i] += sqr(v); - if (h->chist) - { - h->chist[i] = 0; - } - h->cn[i] = 0; - } - if (h->blockfp) - { - fprintf(h->blockfp, "\n"); - } - h->nblocks++; - h->nframes = 0; -} - -/*! - * \param[in,out] h Histogram data structure. - * - * Should be called after each frame of data. - * - * \see gmx_histogram_finish() - */ -void -gmx_histogram_finish_frame(gmx_histogram_t *h) -{ - h->nframes++; - if (h->nframes == h->bsize) - { - finish_histogram_block(h); - } -} - -/*! - * \param[in,out] h Histogram data structure. - * - * Normalizes the histogram. - * Should be called after all frames have been sampled and before any - * post-processing of the histogram. - * If block size has not been set with gmx_histogram_set_blocksize(), - * this function does no normalization, but it should still be called, - * otherwise the \c gmx_histogram_t::hist array will contain only zeros. - * - * \see gmx_histogram_finish_frame() - */ -void -gmx_histogram_finish(gmx_histogram_t *h) -{ - int i; - - if (h->nframes > 0 || h->bsize == 0) - { - if (h->nframes < h->bsize) - { - fprintf(stderr, "Last block smaller (%d frames) than the target size (%d frames) skipped \n", - h->nframes, h->bsize); - } - else - { - finish_histogram_block(h); - } - } - if (h->nblocks == 0) - { - return; - } - - for (i = 0; i <= h->nbins; ++i) - { - h->hist[i] /= h->nblocks; - h->histerr[i] /= h->nblocks; - h->histerr[i] = sqrt((h->histerr[i] - sqr(h->hist[i])) / h->nblocks); - } -} - -/*! - * \param[out] destp Destination histogram. - * \param[in] src Source histogram. - * \param[in] bIntegerBins Control bin center position. - * - * Uses \p src to create a new histogram that has double the binwidth. - * If \p bIntegerBins is TRUE, the first new bin will be centered at - * \c src->start, otherwise the left edge will be at \c src->start. - * - * This function is mostly useful if one needs to sample both the - * histogram and the cumulative histogram at same points. - * To achieve this, first sample a histogram with half the desired - * binwidth, and then use this function to create two different verions - * of it with different values of \p bIntegerBins. - * gmx_histogram_write_array() and gmx_histogram_write_cum_array() can - * now used to write out the values correctly at identical values. - * - * Should be called only after gmx_histogram_finish() has been called. - * \p src should have been sampled without gmx_histogram_set_integerbins(). - */ -void -gmx_histogram_resample_dblbw(gmx_histogram_t **destp, gmx_histogram_t *src, - gmx_bool bIntegerBins) -{ - gmx_histogram_t *dest; - int i, j; - real v, ve; - - gmx_histogram_create(destp, src->type, src->nbins / 2); - dest = *destp; - sfree(dest->chist); dest->chist = NULL; - sfree(dest->cn); dest->cn = NULL; - gmx_histogram_set_binwidth(dest, src->start, src->binwidth * 2); - gmx_histogram_set_integerbins(dest, bIntegerBins); - - for (i = j = 0; i < dest->nbins; ++i) - { - if (bIntegerBins && i == 0) - { - v = src->hist[0]; - ve = sqr(src->histerr[0]); - ++j; - } - else - { - v = src->hist[j] + src->hist[j+1]; - ve = sqr(src->histerr[j]) + sqr(src->histerr[j+1]); - j += 2; - } - ve = sqrt(ve); - dest->hist[i] = v; - dest->histerr[i] = ve; - } - dest->hist[dest->nbins] = 0; - dest->histerr[dest->nbins] = 0; -} - -/*! - * \param[out] destp Destination histogram. - * \param[in] src Source histogram. - * - * Makes a clone of \p src for post-processing. - */ -void -gmx_histogram_clone(gmx_histogram_t **destp, gmx_histogram_t *src) -{ - gmx_histogram_t *dest; - - snew(dest, 1); - memcpy(dest, src, sizeof(*dest)); - - /* These are not needed in post-processing */ - dest->blockfp = NULL; - dest->chist = NULL; - dest->cn = NULL; - - /* Make a deep copy of the actual histograms */ - snew(dest->hist, src->nbins+1); - snew(dest->histerr, src->nbins+1); - memcpy(dest->hist, src->hist, (src->nbins+1)*sizeof(real)); - memcpy(dest->histerr, src->histerr, (src->nbins+1)*sizeof(real)); - - *destp = dest; -} - -/*! - * \param[in,out] h Histogram to normalize. - * - * Normalizes the histogram such that its integral equals one. - */ -void -gmx_histogram_normalize_prob(gmx_histogram_t *h) -{ - int i; - real sum; - real normfac; - - sum = 0; - for (i = 0; i <= h->nbins; ++i) - { - sum += h->hist[i]; - } - - normfac = h->invbw / sum; - gmx_histogram_scale(h, normfac); -} - -/*! - * \param[in,out] h Histogram to normalize. - * \param[in] norm Scaling factor. - * - * All bin values are multiplied by \p norm. - */ -void -gmx_histogram_scale(gmx_histogram_t *h, real norm) -{ - int i; - - for (i = 0; i <= h->nbins; ++i) - { - h->hist[i] *= norm; - h->histerr[i] *= norm; - } -} - -/*! - * \param[in,out] h Histogram to normalize. - * \param[in] norm Scaling vector. - * - * The i'th bin is multiplied by \p norm[i]. - */ -void -gmx_histogram_scale_vec(gmx_histogram_t *h, real norm[]) -{ - int i; - - for (i = 0; i < h->nbins; ++i) - { - h->hist[i] *= norm[i]; - h->histerr[i] *= norm[i]; - } - h->hist[h->nbins] *= norm[h->nbins-1]; - h->histerr[h->nbins] *= norm[h->nbins-1]; -} - -/*! \brief - * Makes some checks on output histograms and finds the maximum number of bins. - * - * \param[in] n Number of histograms in \p h. - * \param[in] h Array of histograms. - * \param[out] nbins Pointer to a value that will receive the maximum number - * of bins in \p h. - */ -static void -prepare_output(int n, gmx_histogram_t *h[], int *nbins) -{ - int j; - - *nbins = 0; - for (j = 0; j < n; ++j) - { - if (!gmx_within_tol(h[j]->start, h[0]->start, GMX_REAL_EPS)) - { - fprintf(stderr, "gmx_ana_histogram_write: histogram start values not identical\n"); - } - if (!gmx_within_tol(h[j]->binwidth, h[0]->binwidth, GMX_REAL_EPS)) - { - fprintf(stderr, "gmx_ana_histogram_write: bin widths not identical\n"); - } - if (*nbins < h[j]->nbins) - { - *nbins = h[j]->nbins; - } - } -} - -/*! - * \param[in] fp Output file. - * \param[in] h Array of histograms to write. - * \param[in] bErrors If TRUE, histogram errors are written. - * - * Convenience wrapper for gmx_histogram_write_array() for writing a single - * histogram. - * - * \see gmx_histogram_write_array() - */ -void -gmx_histogram_write(FILE *fp, gmx_histogram_t *h, gmx_bool bErrors) -{ - gmx_histogram_write_array(fp, 1, &h, TRUE, bErrors); -} - -/*! - * \param[in] fp Output file. - * \param[in] n Number of histograms in \p h. - * \param[in] h Array of histograms to write. - * \param[in] bValue If TRUE, histogram values are written. - * \param[in] bErrors If TRUE, histogram errors are written. - * - * All the histograms in the array \p h should have the same bin widths and - * left edges, otherwise the behavior is undefined. - * The output format is one bin per line. The first number on the line is the - * center of the bin, followed by one or two numbers for each histogram - * (depending on \p bValue and \p bErrors). - * If both \p bValue and \p bErrors are both TRUE, the values are written - * before the errors. - */ -void -gmx_histogram_write_array(FILE *fp, int n, gmx_histogram_t *h[], - gmx_bool bValue, gmx_bool bErrors) -{ - int i, j, nbins; - - prepare_output(n, h, &nbins); - for (i = 0; i < nbins; ++i) - { - fprintf(fp, "%10g", h[0]->start + (i+0.5)*h[0]->binwidth); - for (j = 0; j < n; ++j) - { - if (bValue) - { - fprintf(fp, " %10g", i < h[j]->nbins ? h[j]->hist[i] : 0.0); - } - if (bErrors) - { - fprintf(fp, " %10g", i < h[j]->nbins ? h[j]->histerr[i] : 0.0); - } - } - fprintf(fp, "\n"); - } - fprintf(fp, "\n"); -} - -/*! - * \param[in] fp Output file. - * \param[in] n Number of histograms in \p h. - * \param[in] h Array of histograms to write. - * \param[in] bValue If TRUE, histogram values are written. - * \param[in] bErrors If TRUE, histogram errors are written. - * - * Works as gmx_histogram_write_array(), but writes the cumulative - * histograms. - * The first column in output will be the right edges of the bins. - * - * \note - * Error output is not currently supported (zeros will be written if asked). - * - * \see gmx_histogram_write_array() - */ -void -gmx_histogram_write_cum_array(FILE *fp, int n, gmx_histogram_t *h[], - gmx_bool bValue, gmx_bool bErrors) -{ - int i, j, nbins; - double *sum; - - prepare_output(n, h, &nbins); - snew(sum, n); - - fprintf(fp, "%10g", h[0]->start); - for (j = 0; j < n; ++j) - { - if (bValue) - { - fprintf(fp, " %10g", 0.0); - } - if (bErrors) - { - fprintf(fp, " %10g", 0.0); - } - } - fprintf(fp, "\n"); - for (i = 0; i < nbins; ++i) - { - fprintf(fp, "%10g", h[0]->start + (i+1)*h[0]->binwidth); - for (j = 0; j < n; ++j) - { - sum[j] += i < h[j]->nbins ? h[j]->hist[i] : 0.0; - if (bValue) - { - fprintf(fp, " %10g", sum[j]); - } - /* TODO: Error output not implemented */ - if (bErrors) - { - fprintf(fp, " %10g", 0.0); - } - } - fprintf(fp, "\n"); - } - fprintf(fp, "\n"); - - sfree(sum); -} diff --git a/src/gromacs/legacyheaders/gmx_statistics.h b/src/gromacs/legacyheaders/gmx_statistics.h deleted file mode 100644 index 22379369c8..0000000000 --- a/src/gromacs/legacyheaders/gmx_statistics.h +++ /dev/null @@ -1,178 +0,0 @@ -/* - * 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) 2010, 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. - */ - -#ifndef _GMX_STATS_H -#define _GMX_STATS_H - -#ifdef __cplusplus -extern "C" { -#endif - -#include "typedefs.h" - -typedef struct gmx_stats *gmx_stats_t; - -/* Error codes returned by the routines */ -enum { - estatsOK, estatsNO_POINTS, estatsNO_MEMORY, estatsERROR, - estatsINVALID_INPUT, estatsNOT_IMPLEMENTED, estatsNR -}; - -enum { - elsqWEIGHT_NONE, elsqWEIGHT_X, elsqWEIGHT_Y, - elsqWEIGHT_XY, elsqWEIGHT_NR -}; - -enum { - ehistoX, ehistoY, ehistoNR -}; - -gmx_stats_t gmx_stats_init(); - -int gmx_stats_done(gmx_stats_t stats); - -/* Remove outliers from a straight line, where level in units of - sigma. Level needs to be larger than one obviously. */ -int gmx_stats_remove_outliers(gmx_stats_t stats, double level); - - - -int gmx_stats_add_point(gmx_stats_t stats, double x, double y, - double dx, double dy); - -/* The arrays dx and dy may be NULL if no uncertainties are available, - in that case zero uncertainties will be assumed. */ -int gmx_stats_add_points(gmx_stats_t stats, int n, real *x, real *y, - real *dx, real *dy); - -/* Return the data points one by one. Return estatsOK while there are - more points, and returns estatsNOPOINTS when the last point has - been returned. Should be used in a while loop. Variables for either - pointer may be NULL, in which case the routine can be used as an - expensive point counter. - If level > 0 then the outliers outside level*sigma are reported - only. */ -int gmx_stats_get_point(gmx_stats_t stats, real *x, real *y, - real *dx, real *dy, real level); - -/* Fit the data to y = ax + b, possibly weighted, if uncertainties - have been input. Returns slope in *a and intercept in b, *return - sigmas in *da and *db respectively. Returns normalized *quality of - fit in *chi2 and correlation of fit with data in Rfit. chi2, Rfit, - da and db may be NULL. */ -int gmx_stats_get_ab(gmx_stats_t stats, int weight, - real *a, real *b, - real *da, real *db, real *chi2, real *Rfit); - -/* Fit the data to y = ax, possibly weighted, if uncertainties have - been input. Returns slope in *a, sigma in a in *da, and normalized - quality of fit in *chi2 and correlation of fit with data in - Rfit. chi2, Rfit and da may be NULL. */ -int gmx_stats_get_a(gmx_stats_t stats, int weight, - real *a, real *da, real *chi2, real *Rfit); - -/* Return the correlation coefficient between the data (x and y) as - input to the structure. */ -int gmx_stats_get_corr_coeff(gmx_stats_t stats, real *R); - -/* Returns the root mean square deviation between x and y values. */ -int gmx_stats_get_rmsd(gmx_stats_t gstats, real *rmsd); - -int gmx_stats_get_npoints(gmx_stats_t stats, int *N); - -int gmx_stats_get_average(gmx_stats_t stats, real *aver); - -int gmx_stats_get_sigma(gmx_stats_t stats, real *sigma); - -int gmx_stats_get_error(gmx_stats_t stats, real *error); - -/* Get all three of the above. Pointers may be null, in which case no - assignment will be done. */ -int gmx_stats_get_ase(gmx_stats_t gstats, real *aver, real *sigma, real *error); - -/* Dump the x, y, dx, dy data to a text file */ -int gmx_stats_dump_xy(gmx_stats_t gstats, FILE *fp); - -/* Make a histogram of the data present. Uses either bindwith to - determine the number of bins, or nbins to determine the binwidth, - therefore one of these should be zero, but not the other. If *nbins = 0 - the number of bins will be returned in this variable. ehisto should be one of - ehistoX or ehistoY. If - normalized not equal to zero, the integral of the histogram will be - normalized to one. The output is in two arrays, *x and *y, to which - you should pass a pointer. Memory for the arrays will be allocated - as needed. Function returns one of the estats codes. */ -int gmx_stats_make_histogram(gmx_stats_t gstats, real binwidth, int *nbins, - int ehisto, - int normalized, real **x, real **y); - -/* Return message belonging to error code */ -const char *gmx_stats_message(int estats); - -/**************************************************** - * Some statistics utilities for convenience: useful when a complete data - * set is available already from another source, e.g. an xvg file. - ****************************************************/ -int lsq_y_ax(int n, real x[], real y[], real *a); -/* Fit a straight line y=ax thru the n data points x, y, return the - slope in *a. Return value can be estatsOK, or something else. */ - -int lsq_y_ax_b(int n, real x[], real y[], real *a, real *b, real *r, - real *chi2); -/* Fit a straight line y=ax+b thru the n data points x,y. - * Returns the "fit quality" sigma = sqrt(chi^2/(n-2)). - * The correlation coefficient is returned in r. - */ - -int lsq_y_ax_b_xdouble(int n, double x[], real y[], - real *a, real *b, real *r, real *chi2); -/* As lsq_y_ax_b, but with x in double precision. - */ - -int lsq_y_ax_b_error(int n, real x[], real y[], real dy[], - real *a, real *b, real *da, real *db, - real *r, real *chi2); -/* Fit a straight line y=ax+b thru the n data points x,y, with sigma dy - * Returns the "fit quality" sigma = sqrt(chi^2/(n-2)). - * The correlation coefficient is returned in r. - */ - -#ifdef __cplusplus -} -#endif - -#endif diff --git a/src/gromacs/legacyheaders/histogram.h b/src/gromacs/legacyheaders/histogram.h deleted file mode 100644 index bb71486ff0..0000000000 --- a/src/gromacs/legacyheaders/histogram.h +++ /dev/null @@ -1,220 +0,0 @@ -/* - * This file is part of the GROMACS molecular simulation package. - * - * Copyright (c) 2009,2010, 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. - */ -/*! \file - * \brief API for calculation of histograms with error estimates. - * - * The API is documented in more detail on a separate page: - * \ref histograms - * - * The functions within this file can be used and developed independently of - * the other parts of the library. - * Other parts of the library do not reference these functions. - */ -#ifndef HISTOGRAM_H -#define HISTOGRAM_H - -#include "typedefs.h" - -#ifdef __cplusplus -extern "C" -{ -#endif - -/** Type of histogram. */ -typedef enum -{ - /*! \brief - * Simple histogram. - * - * Use gmx_histogram_increment() or gmx_histogram_increment_bin() - * to sample. - */ - HIST_SIMPLE, - /*! \brief - * Weighted histogram where different points contribute different amounts. - * - * Use gmx_histogram_add() or gmx_histogram_add_to_bin() to sample. - */ - HIST_WEIGHT, - /*! \brief - * Calculate averages within each bin. - * - * Use gmx_histogram_add_item() or gmx_histogram_add_item_to_bin() to sample. - */ - HIST_BINAVER -} e_histogram_t; - -/** Whether bins are centered at integer values. */ -#define HIST_INTEGERBINS 1 -/** Whether the values outside the range should be included in the histogram. */ -#define HIST_ALL 2 -/** Whether the initialization used binwidths. */ -#define HIST_INITBW 128 - -/** Stores data for a histogram. */ -typedef struct gmx_histogram_t gmx_histogram_t; - -/*! \name Initialization functions - */ -/*@{*/ -/** Initialize calculation of a histogram. */ -int -gmx_histogram_create(gmx_histogram_t **h, e_histogram_t type, int nbins); -/** Initialize calculation of a histogram for a range. */ -int -gmx_histogram_create_range(gmx_histogram_t **h, e_histogram_t type, - real start, real end, real binw, gmx_bool bIntegerBins); -/** Clears the bins in the histogram. */ -void -gmx_histogram_clear(gmx_histogram_t *h); -/** Frees the memory allocated for a histogram. */ -void -gmx_histogram_free(gmx_histogram_t *h); -/** Sets histogram range using a starting point and a bin width. */ -int -gmx_histogram_set_binwidth(gmx_histogram_t *h, real start, real binw); -/** Sets histogram range using endpoint values. */ -int -gmx_histogram_set_range(gmx_histogram_t *h, real start, real end); -/** Sets histogram bins to center at integer values. */ -void -gmx_histogram_set_integerbins(gmx_histogram_t *h, gmx_bool bIntegerBins); -/** Sets histogram to include outlying values in the bins at the edges. */ -void -gmx_histogram_set_all(gmx_histogram_t *h, gmx_bool bAll); -/** Sets block size for histogram averaging. */ -int -gmx_histogram_set_blocksize(gmx_histogram_t *h, int bsize); -/** Sets output file for block histograms. */ -int -gmx_histogram_set_block_output(gmx_histogram_t *h, FILE *fp); -/*@}*/ - -/*! \name Access functions - */ -/*@{*/ -/** Finds the histogram bin corresponding to a value. */ -int -gmx_histogram_find_bin(gmx_histogram_t *h, real pos); -/** Returns the number of bins in a histogram. */ -int -gmx_histogram_get_nbins(gmx_histogram_t *h); -/** Returns the bin width of a histogram. */ -real -gmx_histogram_get_binwidth(gmx_histogram_t *h); -/** Returns the value of the histogram at a certain position. */ -void -gmx_histogram_get_value(gmx_histogram_t *h, real pos, double *val, double *err); -/** Returns the value of the histogram in a certain bin. */ -void -gmx_histogram_get_bin_value(gmx_histogram_t *h, int bin, double *val, double *err); -/** Returns an array of values for the histogram. */ -double * -gmx_histogram_get_values(gmx_histogram_t *h); -/** Returns an array of error values for the histogram. */ -double * -gmx_histogram_get_errors(gmx_histogram_t *h); -/*@}*/ - -/*! \name Sampling functions - */ -/*@{*/ -/** Increments the count in a histogram bin corresponding to \p pos. */ -void -gmx_histogram_increment(gmx_histogram_t *h, real pos); -/** Increments the count in a histogram bin. */ -void -gmx_histogram_increment_bin(gmx_histogram_t *h, int bin); - -/** Adds a value to a histogram bin corresponding to \p pos. */ -void -gmx_histogram_add(gmx_histogram_t *h, real pos, double value); -/** Adds a value to a histogram bin. */ -void -gmx_histogram_add_to_bin(gmx_histogram_t *h, int bin, double value); - -/** Adds a value to a histogram bin corresponding to \p pos. */ -void -gmx_histogram_add_item(gmx_histogram_t *h, real pos, double value); -/** Adds a value to a histogram bin. */ -void -gmx_histogram_add_item_to_bin(gmx_histogram_t *h, int bin, double value); - -/** Finishes histogram sampling for a frame. */ -void -gmx_histogram_finish_frame(gmx_histogram_t *h); -/** Normalizes a histogram. */ -void -gmx_histogram_finish(gmx_histogram_t *h); -/*@}*/ - - -/*! \name Post-processing functions - */ -/*@{*/ -/** Creates a new histogram with double the binwidth. */ -void -gmx_histogram_resample_dblbw(gmx_histogram_t **dest, gmx_histogram_t *src, - gmx_bool bIntegerBins); -/** Makes a clone of a histogram. */ -void -gmx_histogram_clone(gmx_histogram_t **dest, gmx_histogram_t *src); -/** Normalizes a histogram to a probability distribution. */ -void -gmx_histogram_normalize_prob(gmx_histogram_t *h); -/** Scales a histogram with a custom normalization factor. */ -void -gmx_histogram_scale(gmx_histogram_t *h, real norm); -/** Scales a histogram with a custom non-uniform normalization factor. */ -void -gmx_histogram_scale_vec(gmx_histogram_t *h, real norm[]); -/** Writes a single histogram to a file. */ -void -gmx_histogram_write(FILE *fp, gmx_histogram_t *h, gmx_bool bErrors); -/** Writes a set of histograms to a file. */ -void -gmx_histogram_write_array(FILE *fp, int n, gmx_histogram_t *h[], - gmx_bool bValue, gmx_bool bErrors); -/** Writes a set of cumulative histograms to a file. */ -void -gmx_histogram_write_cum_array(FILE *fp, int n, gmx_histogram_t *h[], - gmx_bool bValue, gmx_bool bErrors); -/*@}*/ - -#ifdef __cplusplus -} -#endif - -#endif diff --git a/src/gromacs/statistics/CMakeLists.txt b/src/gromacs/statistics/CMakeLists.txt new file mode 100644 index 0000000000..2f7f9f9c65 --- /dev/null +++ b/src/gromacs/statistics/CMakeLists.txt @@ -0,0 +1,39 @@ +# +# This file is part of the GROMACS molecular simulation package. +# +# Copyright (c) 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. + +file(GLOB STATISTICS_SOURCES statistics.c) + +set(LIBGROMACS_SOURCES + ${LIBGROMACS_SOURCES} ${STATISTICS_SOURCES} PARENT_SCOPE) + diff --git a/src/gromacs/gmxlib/statistics/gmx_statistics.c b/src/gromacs/statistics/statistics.c similarity index 99% rename from src/gromacs/gmxlib/statistics/gmx_statistics.c rename to src/gromacs/statistics/statistics.c index fb8726ffdc..0a984ae462 100644 --- a/src/gromacs/gmxlib/statistics/gmx_statistics.c +++ b/src/gromacs/statistics/statistics.c @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2012, by the GROMACS development team, led by + * Copyright (c) 2012,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. @@ -34,15 +34,14 @@ * 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 "statistics.h" #ifdef HAVE_CONFIG_H #include #endif - #include #include "typedefs.h" #include "smalloc.h" #include "vec.h" -#include "gmx_statistics.h" static int gmx_dnint(double x) { diff --git a/src/gromacs/statistics/statistics.h b/src/gromacs/statistics/statistics.h new file mode 100644 index 0000000000..4e0faa93e0 --- /dev/null +++ b/src/gromacs/statistics/statistics.h @@ -0,0 +1,331 @@ +/* + * 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) 2010,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. + */ + +#ifndef GMX_STATISTICS_H +#define GMX_STATISTICS_H + +/*! \libinternal \file + * + * \brief + * Declares simple statistics toolbox + * + * \authors David van der Spoel + * + * \inlibraryapi + */ +#ifdef __cplusplus +extern "C" { +#endif +#include +#include "types/simple.h" + +//! Abstract container type +typedef struct gmx_stats *gmx_stats_t; + +//! Error codes returned by the routines +enum { + estatsOK, estatsNO_POINTS, estatsNO_MEMORY, estatsERROR, + estatsINVALID_INPUT, estatsNOT_IMPLEMENTED, estatsNR +}; + +//! Enum for statistical weights +enum { + elsqWEIGHT_NONE, elsqWEIGHT_X, elsqWEIGHT_Y, + elsqWEIGHT_XY, elsqWEIGHT_NR +}; + +//! Enum determining which coordinate to histogram +enum { + ehistoX, ehistoY, ehistoNR +}; + +/*! \brief + * Initiate a data structure + * \return the data structure + */ +gmx_stats_t gmx_stats_init(); + +/*! \brief + * Destroy a data structure + * \param stats The data structure + * \return error code + */ +int gmx_stats_done(gmx_stats_t stats); + +/*! \brief + * Remove outliers from a straight line, where level in units of + * sigma. Level needs to be larger than one obviously. + * \param[in] stats The data structure + * \param[in] level The sigma level + * \return error code + */ +int gmx_stats_remove_outliers(gmx_stats_t stats, double level); + +/*! \brief + * Add a point to the data set + * \param[in] stats The data structure + * \param[in] x The x value + * \param[in] y The y value + * \param[in] dx The error in the x value + * \param[in] dy The error in the y value + * \return error code + */ +int gmx_stats_add_point(gmx_stats_t stats, double x, double y, + double dx, double dy); + +/*! \brief + * Add a series of datapoints at once. The arrays dx and dy may + * be NULL in that case zero uncertainties will be assumed. + * + * \param[in] stats The data structure + * \param[in] n Number of points + * \param[in] x The array of x values + * \param[in] y The array of y values + * \param[in] dx The error in the x value + * \param[in] dy The error in the y value + * \return error code + */ +int gmx_stats_add_points(gmx_stats_t stats, int n, real *x, real *y, + real *dx, real *dy); + +/*! \brief + * Delivers data points from the statistics. + * + * Should be used in a while loop. Variables for either + * pointer may be NULL, in which case the routine can be used as an + * expensive point counter. + * Return the data points one by one. Return estatsOK while there are + * more points, and returns estatsNOPOINTS when the last point has + * been returned. + * If level > 0 then the outliers outside level*sigma are reported + * only. + * \param[in] stats The data structure + * \param[out] x The array of x values + * \param[out] y The array of y values + * \param[out] dx The error in the x value + * \param[out] dy The error in the y value + * \param[in] level sigma level (see above) + * \return error code + */ +int gmx_stats_get_point(gmx_stats_t stats, real *x, real *y, + real *dx, real *dy, real level); + +/*! \brief + * Fit the data to y = ax + b, possibly weighted, if uncertainties + * have been input. da and db may be NULL. + * \param[in] stats The data structure + * \param[in] weight type of weighting + * \param[out] a slope + * \param[out] b intercept + * \param[out] da sigma in a + * \param[out] db sigma in b + * \param[out] chi2 normalized quality of fit + * \param[out] Rfit correlation coefficient + * \return error code + */ +int gmx_stats_get_ab(gmx_stats_t stats, int weight, + real *a, real *b, + real *da, real *db, real *chi2, real *Rfit); + +/*! \brief + * Fit the data to y = ax, possibly weighted, if uncertainties have + * have been input. da and db may be NULL. + * \param[in] stats The data structure + * \param[in] weight type of weighting + * \param[out] a slope + * \param[out] da sigma in a + * \param[out] chi2 normalized quality of fit + * \param[out] Rfit correlation coefficient + * \return error code + */ +int gmx_stats_get_a(gmx_stats_t stats, int weight, + real *a, real *da, real *chi2, real *Rfit); + +/*! \brief + * Get the correlation coefficient. + * \param[in] stats The data structure + * \param[out] R the correlation coefficient between the data (x and y) as input to the structure. + * \return error code + */ +int gmx_stats_get_corr_coeff(gmx_stats_t stats, real *R); + +/*! \brief + * Get the root mean square deviation. + * \param[in] stats The data structure + * \param[out] rmsd the root mean square deviation between x and y values. + * \return error code + */ +int gmx_stats_get_rmsd(gmx_stats_t stats, real *rmsd); + +/*! \brief + * Get the number of points. + * \param[in] stats The data structure + * \param[out] N number of data points + * \return error code + */ +int gmx_stats_get_npoints(gmx_stats_t stats, int *N); + +/*! \brief + * Computes and returns the average value. + * \param[in] stats The data structure + * \param[out] aver Average value + * \return error code + */ +int gmx_stats_get_average(gmx_stats_t stats, real *aver); + +/*! \brief + * Computes and returns the standard deviation. + * \param[in] stats The data structure + * \param[out] sigma Standard deviation + * \return error code + */ +int gmx_stats_get_sigma(gmx_stats_t stats, real *sigma); + +/*! \brief + * Computes and returns the standard error. + * \param[in] stats The data structure + * \param[out] error Standard error + * \return error code + */ +int gmx_stats_get_error(gmx_stats_t stats, real *error); + +/*! \brief + * Pointers may be null, in which case no assignment will be done. + * \param[in] stats The data structure + * \param[out] aver Average value + * \param[out] sigma Standard deviation + * \param[out] error Standard error + * \return error code + */ +int gmx_stats_get_ase(gmx_stats_t stats, real *aver, real *sigma, real *error); + +/*! \brief + * Dump the x, y, dx, dy data to a text file + * \param[in] stats The data structure + * \param[in] fp File pointer + * \return error code + */ +int gmx_stats_dump_xy(gmx_stats_t stats, FILE *fp); + +/*! \brief + * Make a histogram of the data present. + * + * Uses either binwidth to + * determine the number of bins, or nbins to determine the binwidth, + * therefore one of these should be zero, but not the other. If *nbins = 0 + * the number of bins will be returned in this variable. ehisto should be one of + * ehistoX or ehistoY. If + * normalized not equal to zero, the integral of the histogram will be + * normalized to one. The output is in two arrays, *x and *y, to which + * you should pass a pointer. Memory for the arrays will be allocated + * as needed. Function returns one of the estats codes. + * \param[in] stats The data structure + * \param[in] binwidth For the histogram + * \param[in] nbins Number of bins + * \param[in] ehisto Type (see enum above) + * \param[in] normalized see above + * \param[out] x see above + * \param[out] y see above + * \return error code + */ +int gmx_stats_make_histogram(gmx_stats_t stats, real binwidth, int *nbins, + int ehisto, + int normalized, real **x, real **y); + +/*! \brief + * Return message belonging to error code + * \param[in] estats error code + */ +const char *gmx_stats_message(int estats); + +/**************************************************** + * Some statistics utilities for convenience: useful when a complete data + * set is available already from another source, e.g. an xvg file. + ****************************************************/ +/*! \brief + * Fit a straight line y=ax thru the n data points x, y, return the + * slope in *a. + * \param[in] n number of points + * \param[in] x data points x + * \param[in] y data point y + * \param[out] a slope + * \return error code + */ +int lsq_y_ax(int n, real x[], real y[], real *a); + +/*! \brief + * Fit a straight line y=ax+b thru the n data points x, y. + * \param[in] n number of points + * \param[in] x data points x + * \param[in] y data point y + * \param[out] a slope + * \param[out] b intercept + * \param[out] r correlation coefficient + * \param[out] chi2 quality of fit + * \return error code + */ +int lsq_y_ax_b(int n, real x[], real y[], real *a, real *b, real *r, + real *chi2); + +/*! \copydoc lsq_y_ax_b + */ +int lsq_y_ax_b_xdouble(int n, double x[], real y[], + real *a, real *b, real *r, real *chi2); + +/*! \brief + * Fit a straight line y=ax+b thru the n data points x, y. + * \param[in] n number of points + * \param[in] x data points x + * \param[in] y data point y + * \param[in] dy uncertainty in data point y + * \param[out] a slope + * \param[out] b intercept + * \param[out] da error in slope + * \param[out] db error in intercept + * \param[out] r correlation coefficient + * \param[out] chi2 quality of fit + * \return error code + */ +int lsq_y_ax_b_error(int n, real x[], real y[], real dy[], + real *a, real *b, real *da, real *db, + real *r, real *chi2); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/src/gromacs/gmxlib/statistics/gmx_statistics_test.c b/src/gromacs/statistics/statistics_test.c similarity index 98% rename from src/gromacs/gmxlib/statistics/gmx_statistics_test.c rename to src/gromacs/statistics/statistics_test.c index 6742737aa2..8e00049661 100644 --- a/src/gromacs/gmxlib/statistics/gmx_statistics_test.c +++ b/src/gromacs/statistics/statistics_test.c @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2008, The GROMACS development team. - * Copyright (c) 2011, by the GROMACS development team, led by + * 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. @@ -43,7 +43,7 @@ #include "smalloc.h" #include "vec.h" #include "gmx_random.h" -#include "gmx_statistics.h" +#include "statistics.h" static void horizontal() {