src/gromacs/analysisdata/modules/average.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2010,2011,2012,2013,2014,2019, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source 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  * Implements gmx::AnalysisDataAverageModule.
  38  *
  39  * \author Teemu Murtola <teemu.murtola@gmail.com>
  40  * \ingroup module_analysisdata
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "average.h"
  45
  46 #include <cmath>
  47
  48 #include <algorithm>
  49 #include <vector>
  50
  51 #include "gromacs/analysisdata/dataframe.h"
  52 #include "gromacs/analysisdata/datastorage.h"
  53
  54 #include "frameaverager.h"
  55
  56 namespace gmx
  57 {
  58
  59 /********************************************************************
  60  * AnalysisDataAverageModule
  61  */
  62
  63 class AnalysisDataAverageModule::Impl
  64 {
  65 public:
  66     Impl() : bDataSets_(false) {}
  67
  68     //! Averaging helper objects for each input data set.
  69     std::vector<AnalysisDataFrameAverager> averagers_;
  70     //! Whether to average all columns in a data set into a single value.
  71     bool bDataSets_;
  72 };
  73
  74 AnalysisDataAverageModule::AnalysisDataAverageModule() : impl_(new Impl()) {}
  75
  76 AnalysisDataAverageModule::~AnalysisDataAverageModule() {}
  77
  78 void AnalysisDataAverageModule::setAverageDataSets(bool bDataSets)
  79 {
  80     impl_->bDataSets_ = bDataSets;
  81 }
  82
  83 int AnalysisDataAverageModule::flags() const
  84 {
  85     return efAllowMultipoint | efAllowMulticolumn | efAllowMissing | efAllowMultipleDataSets;
  86 }
  87
  88 void AnalysisDataAverageModule::dataStarted(AbstractAnalysisData* data)
  89 {
  90     if (impl_->bDataSets_)
  91     {
  92         setColumnCount(1);
  93         setRowCount(data->dataSetCount());
  94         impl_->averagers_.resize(1);
  95         impl_->averagers_[0].setColumnCount(data->dataSetCount());
  96     }
  97     else
  98     {
  99         setColumnCount(data->dataSetCount());
 100         impl_->averagers_.resize(data->dataSetCount());
 101         int rowCount = 0;
 102         for (int i = 0; i < data->dataSetCount(); ++i)
 103         {
 104             impl_->averagers_[i].setColumnCount(data->columnCount(i));
 105             rowCount = std::max(rowCount, data->columnCount(i));
 106         }
 107         setRowCount(rowCount);
 108     }
 109 }
 110
 111 void AnalysisDataAverageModule::frameStarted(const AnalysisDataFrameHeader& /*header*/) {}
 112
 113 void AnalysisDataAverageModule::pointsAdded(const AnalysisDataPointSetRef& points)
 114 {
 115     if (impl_->bDataSets_)
 116     {
 117         const int dataSet = points.dataSetIndex();
 118         for (int i = 0; i < points.columnCount(); ++i)
 119         {
 120             if (points.present(i))
 121             {
 122                 impl_->averagers_[0].addValue(dataSet, points.y(i));
 123             }
 124         }
 125     }
 126     else
 127     {
 128         impl_->averagers_[points.dataSetIndex()].addPoints(points);
 129     }
 130 }
 131
 132 void AnalysisDataAverageModule::frameFinished(const AnalysisDataFrameHeader& /*header*/) {}
 133
 134 void AnalysisDataAverageModule::dataFinished()
 135 {
 136     allocateValues();
 137     for (int i = 0; i < columnCount(); ++i)
 138     {
 139         impl_->averagers_[i].finish();
 140         int j = 0;
 141         for (; j < impl_->averagers_[i].columnCount(); ++j)
 142         {
 143             value(j, i).setValue(impl_->averagers_[i].average(j),
 144                                  std::sqrt(impl_->averagers_[i].variance(j)));
 145         }
 146         for (; j < rowCount(); ++j)
 147         {
 148             value(j, i).setValue(0.0, 0.0, false);
 149         }
 150     }
 151     valuesReady();
 152 }
 153
 154 real AnalysisDataAverageModule::average(int dataSet, int column) const
 155 {
 156     if (impl_->bDataSets_)
 157     {
 158         GMX_ASSERT(column == 0, "Column should be zero with setAverageDataSets(true)");
 159         std::swap(dataSet, column);
 160     }
 161     return value(column, dataSet).value();
 162 }
 163
 164 real AnalysisDataAverageModule::standardDeviation(int dataSet, int column) const
 165 {
 166     if (impl_->bDataSets_)
 167     {
 168         GMX_ASSERT(column == 0, "Column should be zero with setAverageDataSets(true)");
 169         std::swap(dataSet, column);
 170     }
 171     return value(column, dataSet).error();
 172 }
 173
 174 int AnalysisDataAverageModule::sampleCount(int dataSet, int column) const
 175 {
 176     if (impl_->bDataSets_)
 177     {
 178         GMX_ASSERT(column == 0, "Column should be zero with setAverageDataSets(true)");
 179         std::swap(dataSet, column);
 180     }
 181     return impl_->averagers_[dataSet].sampleCount(column);
 182 }
 183
 184
 185 /********************************************************************
 186  * AnalysisDataFrameAverageModule
 187  */
 188
 189 class AnalysisDataFrameAverageModule::Impl
 190 {
 191 public:
 192     //! Storage implementation object.
 193     AnalysisDataStorage storage_;
 194     //! Number of samples in a frame for each data set.
 195     std::vector<int> sampleCount_;
 196 };
 197
 198 AnalysisDataFrameAverageModule::AnalysisDataFrameAverageModule() : impl_(new Impl()) {}
 199
 200 AnalysisDataFrameAverageModule::~AnalysisDataFrameAverageModule() {}
 201
 202 int AnalysisDataFrameAverageModule::frameCount() const
 203 {
 204     return impl_->storage_.frameCount();
 205 }
 206
 207 int AnalysisDataFrameAverageModule::flags() const
 208 {
 209     return efAllowMultipoint | efAllowMulticolumn | efAllowMissing | efAllowMultipleDataSets;
 210 }
 211
 212 void AnalysisDataFrameAverageModule::dataStarted(AbstractAnalysisData* data)
 213 {
 214     setColumnCount(0, data->dataSetCount());
 215     impl_->sampleCount_.resize(data->dataSetCount());
 216     impl_->storage_.startDataStorage(this, &moduleManager());
 217 }
 218
 219 void AnalysisDataFrameAverageModule::frameStarted(const AnalysisDataFrameHeader& header)
 220 {
 221     AnalysisDataStorageFrame& frame = impl_->storage_.startFrame(header);
 222     for (int i = 0; i < columnCount(); ++i)
 223     {
 224         impl_->sampleCount_[i] = 0;
 225         frame.setValue(i, 0.0);
 226     }
 227 }
 228
 229 void AnalysisDataFrameAverageModule::pointsAdded(const AnalysisDataPointSetRef& points)
 230 {
 231     const int                 dataSet = points.dataSetIndex();
 232     AnalysisDataStorageFrame& frame   = impl_->storage_.currentFrame(points.frameIndex());
 233     for (int i = 0; i < points.columnCount(); ++i)
 234     {
 235         if (points.present(i))
 236         {
 237             // TODO: Consider using AnalysisDataFrameAverager
 238             const real y     = points.y(i);
 239             const real delta = y - frame.value(dataSet);
 240             impl_->sampleCount_[dataSet] += 1;
 241             frame.value(dataSet) += delta / impl_->sampleCount_[dataSet];
 242         }
 243     }
 244 }
 245
 246 void AnalysisDataFrameAverageModule::frameFinished(const AnalysisDataFrameHeader& header)
 247 {
 248     impl_->storage_.finishFrame(header.index());
 249 }
 250
 251 void AnalysisDataFrameAverageModule::dataFinished()
 252 {
 253     impl_->storage_.finishDataStorage();
 254 }
 255
 256 AnalysisDataFrameRef AnalysisDataFrameAverageModule::tryGetDataFrameInternal(int index) const
 257 {
 258     return impl_->storage_.tryGetDataFrame(index);
 259 }
 260
 261 bool AnalysisDataFrameAverageModule::requestStorageInternal(int nframes)
 262 {
 263     return impl_->storage_.requestStorage(nframes);
 264 }
 265
 266 } // namespace gmx