src/gromacs/analysisdata/modules/frameaverager.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013, by the GROMACS development team, led by
   5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   6  * others, as listed in the AUTHORS file in the top-level source
   7  * directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Declares gmx::AnalysisDataFrameAverager.
  38  *
  39  * \author Teemu Murtola <teemu.murtola@gmail.com>
  40  * \ingroup module_analysisdata
  41  */
  42 #ifndef GMX_ANALYSISDATA_MODULES_FRAMEAVERAGER_H
  43 #define GMX_ANALYSISDATA_MODULES_FRAMEAVERAGER_H
  44
  45 #include <vector>
  46
  47 #include "../../legacyheaders/types/simple.h"
  48
  49 #include "../../utility/gmxassert.h"
  50
  51 namespace gmx
  52 {
  53
  54 class AnalysisDataPointSetRef;
  55
  56 /*! \internal \brief
  57  * Helper class for modules that average values over frames.
  58  *
  59  * This class implements common functionality for analysis data modules that
  60  * need to average a set of values over frames.  Currently, it is designed for
  61  * computing averages for each input column independently, but should be
  62  * relatively easy to make more general if required.
  63  *
  64  * This class takes care of accumulating the values and computing their
  65  * variance.  It allows different number of samples for each input column.
  66  * Accumulation is always in double precision and uses a formula that is
  67  * relatively stable numerically.  For now, does nothing fancy,
  68  * but provides ground for other implementation (e.g., related to
  69  * parallelization) that would benefit all such modules.
  70  *
  71  * Methods in this class do not throw unless otherwise indicated.
  72  *
  73  * \ingroup module_analysisdata
  74  */
  75 class AnalysisDataFrameAverager
  76 {
  77     public:
  78         AnalysisDataFrameAverager() : bFinished_(false) {}
  79
  80         /*! \brief
  81          * Sets the number of columns in the input data.
  82          *
  83          * \throws std::bad_alloc if out of memory.
  84          *
  85          * Typically called from AnalysisDataModuleInterface::dataStarted().
  86          *
  87          * Must be called exactly once, before setting calling any other method
  88          * in the class.
  89          */
  90         void setColumnCount(int columnCount);
  91         /*! \brief
  92          * Adds a single value to the average for a given column.
  93          *
  94          * \param[in] index  Index of the column to add the value to.
  95          * \param[in] value  Value to add to the sample.
  96          */
  97         void addValue(int index, real value);
  98         /*! \brief
  99          * Accumulates data from a given point set into the average.
 100          *
 101          * Typically called from AnalysisDataModuleInterface::pointsAdded().
 102          *
 103          * Each call accumulates the values for those columns that are present
 104          * in the point set.  Can be called multiple times for a frame, and
 105          * does not need to be called for every frame.
 106          */
 107         void addPoints(const AnalysisDataPointSetRef &points);
 108         /*! \brief
 109          * Finalizes the calculation of the averages and variances.
 110          *
 111          * Does any computation that is not done during the accumulation in
 112          * addPoints().  Currently, does nothing, but provided as a placeholder
 113          * for more complex implementation.
 114          *
 115          * Typically called from AnalysisDataModuleInterface::dataFinished().
 116          */
 117         void finish();
 118
 119         /*! \brief
 120          * Returns the computed average for a given column.
 121          *
 122          * If called before finish(), the results are undefined.
 123          */
 124         real average(int index) const
 125         {
 126             GMX_ASSERT(index >= 0 && index <= static_cast<int>(values_.size()),
 127                        "Invalid column index");
 128             GMX_ASSERT(bFinished_,
 129                        "Values available only after finished() has been called");
 130             return values_[index].average;
 131         }
 132         /*! \brief
 133          * Returns the computed (sample) variance for a given column.
 134          *
 135          * If called before finish(), the results are undefined.
 136          */
 137         real variance(int index) const
 138         {
 139             GMX_ASSERT(index >= 0 && index <= static_cast<int>(values_.size()),
 140                        "Invalid column index");
 141             GMX_ASSERT(bFinished_,
 142                        "Values available only after finished() has been called");
 143             const AverageItem &item = values_[index];
 144             return item.samples > 1 ? item.squaredSum / (item.samples - 1) : 0.0;
 145         }
 146         /*! \brief
 147          * Returns the number of samples for a given column.
 148          *
 149          * If called before finish(), the results are undefined.
 150          */
 151         int sampleCount(int index) const
 152         {
 153             GMX_ASSERT(index >= 0 && index <= static_cast<int>(values_.size()),
 154                        "Invalid column index");
 155             GMX_ASSERT(bFinished_,
 156                        "Values available only after finished() has been called");
 157             return values_[index].samples;
 158         }
 159
 160     private:
 161         struct AverageItem
 162         {
 163             AverageItem() : average(0.0), squaredSum(0.0), samples(0) {}
 164
 165             //! Average of the values so far.
 166             double               average;
 167             //! Sum of squared deviations from the average for values so far.
 168             double               squaredSum;
 169             //! Number of values so far.
 170             int                  samples;
 171         };
 172
 173         std::vector<AverageItem> values_;
 174         bool                     bFinished_;
 175 };
 176
 177 } // namespace gmx
 178
 179 #endif