if (NOT GMX_BUILD_MDRUN_ONLY)
add_subdirectory(legacyheaders)
add_subdirectory(gmxana)
+ add_subdirectory(statistics)
add_subdirectory(analysisdata)
add_subdirectory(selection)
add_subdirectory(trajectoryanalysis)
*
* 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.
#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"
*
* 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.
#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"
/*
* 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.
#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"
*
* 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.
#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"
#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"
#
# 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.
# 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
+++ /dev/null
-/*
- * 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 <config.h>
-#endif
-
-#include <math.h>
-#include <string.h>
-
-#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);
-}
+++ /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) 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
+++ /dev/null
-/*
- * 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
--- /dev/null
+#
+# 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)
+
*
* 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.
* 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 <config.h>
#endif
-
#include <math.h>
#include "typedefs.h"
#include "smalloc.h"
#include "vec.h"
-#include "gmx_statistics.h"
static int gmx_dnint(double x)
{
--- /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) 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 <david.vanderspoel@icm.uu.se>
+ *
+ * \inlibraryapi
+ */
+#ifdef __cplusplus
+extern "C" {
+#endif
+#include <stdio.h>
+#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
*
* 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.
#include "smalloc.h"
#include "vec.h"
#include "gmx_random.h"
-#include "gmx_statistics.h"
+#include "statistics.h"
static void horizontal()
{
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff