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