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 by the GROMACS development team.
   5  * Copyright (c) 2019,2020, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 /*! \internal \file
  37  * \brief
  38  * Implements gmx::AnalysisDataAverageModule.
  39  *
  40  * \author Teemu Murtola <teemu.murtola@gmail.com>
  41  * \ingroup module_analysisdata
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "average.h"
  46
  47 #include <cmath>
  48
  49 #include <algorithm>
  50 #include <vector>
  51
  52 #include "gromacs/analysisdata/dataframe.h"
  53 #include "gromacs/analysisdata/datastorage.h"
  54
  55 #include "frameaverager.h"
  56
  57 namespace gmx
  58 {
  59
  60 /********************************************************************
  61  * AnalysisDataAverageModule
  62  */
  63
  64 class AnalysisDataAverageModule::Impl
  65 {
  66 public:
  67     Impl() : bDataSets_(false) {}
  68
  69     //! Averaging helper objects for each input data set.
  70     std::vector<AnalysisDataFrameAverager> averagers_;
  71     //! Whether to average all columns in a data set into a single value.
  72     bool bDataSets_;
  73 };
  74
  75 AnalysisDataAverageModule::AnalysisDataAverageModule() : impl_(new Impl()) {}
  76
  77 AnalysisDataAverageModule::~AnalysisDataAverageModule() {}
  78
  79 void AnalysisDataAverageModule::setAverageDataSets(bool bDataSets)
  80 {
  81     impl_->bDataSets_ = bDataSets;
  82 }
  83
  84 int AnalysisDataAverageModule::flags() const
  85 {
  86     return efAllowMultipoint | efAllowMulticolumn | efAllowMissing | efAllowMultipleDataSets;
  87 }
  88
  89 void AnalysisDataAverageModule::dataStarted(AbstractAnalysisData* data)
  90 {
  91     if (impl_->bDataSets_)
  92     {
  93         setColumnCount(1);
  94         setRowCount(data->dataSetCount());
  95         impl_->averagers_.resize(1);
  96         impl_->averagers_[0].setColumnCount(data->dataSetCount());
  97     }
  98     else
  99     {
 100         setColumnCount(data->dataSetCount());
 101         impl_->averagers_.resize(data->dataSetCount());
 102         int rowCount = 0;
 103         for (int i = 0; i < data->dataSetCount(); ++i)
 104         {
 105             impl_->averagers_[i].setColumnCount(data->columnCount(i));
 106             rowCount = std::max(rowCount, data->columnCount(i));
 107         }
 108         setRowCount(rowCount);
 109     }
 110 }
 111
 112 void AnalysisDataAverageModule::frameStarted(const AnalysisDataFrameHeader& /*header*/) {}
 113
 114 void AnalysisDataAverageModule::pointsAdded(const AnalysisDataPointSetRef& points)
 115 {
 116     if (impl_->bDataSets_)
 117     {
 118         const int dataSet = points.dataSetIndex();
 119         for (int i = 0; i < points.columnCount(); ++i)
 120         {
 121             if (points.present(i))
 122             {
 123                 impl_->averagers_[0].addValue(dataSet, points.y(i));
 124             }
 125         }
 126     }
 127     else
 128     {
 129         impl_->averagers_[points.dataSetIndex()].addPoints(points);
 130     }
 131 }
 132
 133 void AnalysisDataAverageModule::frameFinished(const AnalysisDataFrameHeader& /*header*/) {}
 134
 135 void AnalysisDataAverageModule::dataFinished()
 136 {
 137     allocateValues();
 138     for (int i = 0; i < columnCount(); ++i)
 139     {
 140         impl_->averagers_[i].finish();
 141         int j = 0;
 142         for (; j < impl_->averagers_[i].columnCount(); ++j)
 143         {
 144             value(j, i).setValue(impl_->averagers_[i].average(j),
 145                                  std::sqrt(impl_->averagers_[i].variance(j)));
 146         }
 147         for (; j < rowCount(); ++j)
 148         {
 149             value(j, i).setValue(0.0, 0.0, false);
 150         }
 151     }
 152     valuesReady();
 153 }
 154
 155 real AnalysisDataAverageModule::average(int dataSet, int column) const
 156 {
 157     if (impl_->bDataSets_)
 158     {
 159         GMX_ASSERT(column == 0, "Column should be zero with setAverageDataSets(true)");
 160         std::swap(dataSet, column);
 161     }
 162     return value(column, dataSet).value();
 163 }
 164
 165 real AnalysisDataAverageModule::standardDeviation(int dataSet, int column) const
 166 {
 167     if (impl_->bDataSets_)
 168     {
 169         GMX_ASSERT(column == 0, "Column should be zero with setAverageDataSets(true)");
 170         std::swap(dataSet, column);
 171     }
 172     return value(column, dataSet).error();
 173 }
 174
 175 int AnalysisDataAverageModule::sampleCount(int dataSet, int column) const
 176 {
 177     if (impl_->bDataSets_)
 178     {
 179         GMX_ASSERT(column == 0, "Column should be zero with setAverageDataSets(true)");
 180         std::swap(dataSet, column);
 181     }
 182     return impl_->averagers_[dataSet].sampleCount(column);
 183 }
 184
 185
 186 /********************************************************************
 187  * AnalysisDataFrameAverageModule
 188  */
 189
 190 class AnalysisDataFrameAverageModule::Impl
 191 {
 192 public:
 193     //! Storage implementation object.
 194     AnalysisDataStorage storage_;
 195     //! Number of samples in a frame for each data set.
 196     std::vector<int> sampleCount_;
 197 };
 198
 199 AnalysisDataFrameAverageModule::AnalysisDataFrameAverageModule() : impl_(new Impl()) {}
 200
 201 AnalysisDataFrameAverageModule::~AnalysisDataFrameAverageModule() {}
 202
 203 int AnalysisDataFrameAverageModule::frameCount() const
 204 {
 205     return impl_->storage_.frameCount();
 206 }
 207
 208 int AnalysisDataFrameAverageModule::flags() const
 209 {
 210     return efAllowMultipoint | efAllowMulticolumn | efAllowMissing | efAllowMultipleDataSets;
 211 }
 212
 213 void AnalysisDataFrameAverageModule::dataStarted(AbstractAnalysisData* data)
 214 {
 215     setColumnCount(0, data->dataSetCount());
 216     impl_->sampleCount_.resize(data->dataSetCount());
 217     impl_->storage_.startDataStorage(this, &moduleManager());
 218 }
 219
 220 void AnalysisDataFrameAverageModule::frameStarted(const AnalysisDataFrameHeader& header)
 221 {
 222     AnalysisDataStorageFrame& frame = impl_->storage_.startFrame(header);
 223     for (int i = 0; i < columnCount(); ++i)
 224     {
 225         impl_->sampleCount_[i] = 0;
 226         frame.setValue(i, 0.0);
 227     }
 228 }
 229
 230 void AnalysisDataFrameAverageModule::pointsAdded(const AnalysisDataPointSetRef& points)
 231 {
 232     const int                 dataSet = points.dataSetIndex();
 233     AnalysisDataStorageFrame& frame   = impl_->storage_.currentFrame(points.frameIndex());
 234     for (int i = 0; i < points.columnCount(); ++i)
 235     {
 236         if (points.present(i))
 237         {
 238             // TODO: Consider using AnalysisDataFrameAverager
 239             const real y     = points.y(i);
 240             const real delta = y - frame.value(dataSet);
 241             impl_->sampleCount_[dataSet] += 1;
 242             frame.value(dataSet) += delta / impl_->sampleCount_[dataSet];
 243         }
 244     }
 245 }
 246
 247 void AnalysisDataFrameAverageModule::frameFinished(const AnalysisDataFrameHeader& header)
 248 {
 249     impl_->storage_.finishFrame(header.index());
 250 }
 251
 252 void AnalysisDataFrameAverageModule::dataFinished()
 253 {
 254     impl_->storage_.finishDataStorage();
 255 }
 256
 257 AnalysisDataFrameRef AnalysisDataFrameAverageModule::tryGetDataFrameInternal(int index) const
 258 {
 259     return impl_->storage_.tryGetDataFrame(index);
 260 }
 261
 262 bool AnalysisDataFrameAverageModule::requestStorageInternal(int nframes)
 263 {
 264     return impl_->storage_.requestStorage(nframes);
 265 }
 266
 267 } // namespace gmx