src/gromacs/analysisdata/modules/lifetime.cpp

   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  * 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::AnalysisDataLifetimeModule.
  38  *
  39  * \author Teemu Murtola <teemu.murtola@gmail.com>
  40  * \ingroup module_analysisdata
  41  */
  42 #include "lifetime.h"
  43
  44 #include <cmath>
  45
  46 #include <deque>
  47 #include <vector>
  48
  49 #include "gromacs/analysisdata/dataframe.h"
  50 #include "gromacs/analysisdata/datastorage.h"
  51
  52 namespace gmx
  53 {
  54
  55 /********************************************************************
  56  * AnalysisDataLifetimeModule
  57  */
  58
  59 /*! \internal \brief
  60  * Private implementation class for AnalysisDataLifetimeModule.
  61  *
  62  * \ingroup module_analysisdata
  63  */
  64 class AnalysisDataLifetimeModule::Impl
  65 {
  66     public:
  67         //! Container type for storing a histogram during the calculation.
  68         typedef std::deque<int> LifetimeHistogram;
  69
  70         //! Initializes the implementation class with empty/default values.
  71         Impl() : firstx_(0.0), lastx_(0.0), frameCount_(0), bCumulative_(false)
  72         {
  73         }
  74
  75         /*! \brief
  76          * Increments a lifetime histogram with a single lifetime.
  77          *
  78          * \param[in] dataSet   Index of the histogram to increment.
  79          * \param[in] lifetime  Lifetime to add to the histogram.
  80          */
  81         void addLifetime(int dataSet, int lifetime)
  82         {
  83             if (lifetime > 0)
  84             {
  85                 LifetimeHistogram &histogram = lifetimeHistograms_[dataSet];
  86                 if (histogram.size() < static_cast<unsigned>(lifetime))
  87                 {
  88                     histogram.resize(lifetime, 0);
  89                 }
  90                 ++histogram[lifetime - 1];
  91             }
  92         }
  93
  94         //! X value of the first frame (used for determining output spacing).
  95         real                            firstx_;
  96         //! X value of the last frame (used for determining output spacing).
  97         real                            lastx_;
  98         //! Total number of frames (used for normalization and output spacing).
  99         int                             frameCount_;
 100         //! Whether to add subintervals of longer intervals explicitly.
 101         bool                            bCumulative_;
 102         /*! \brief
 103          * Length of current continuously present interval for each data column.
 104          *
 105          * While frame N has been processed, stores the length of an interval
 106          * for each data column where that column has been continuously present
 107          * up to and including frame N.
 108          */
 109         std::vector<std::vector<int> >  currentLifetimes_;
 110         /*! \brief
 111          * Accumulated lifetime histograms for each data set.
 112          */
 113         std::vector<LifetimeHistogram>  lifetimeHistograms_;
 114 };
 115
 116 AnalysisDataLifetimeModule::AnalysisDataLifetimeModule()
 117     : impl_(new Impl())
 118 {
 119 }
 120
 121 AnalysisDataLifetimeModule::~AnalysisDataLifetimeModule()
 122 {
 123 }
 124
 125 void AnalysisDataLifetimeModule::setCumulative(bool bCumulative)
 126 {
 127     impl_->bCumulative_ = bCumulative;
 128 }
 129
 130 int AnalysisDataLifetimeModule::flags() const
 131 {
 132     return efAllowMulticolumn | efAllowMissing | efAllowMultipleDataSets;
 133 }
 134
 135 void
 136 AnalysisDataLifetimeModule::dataStarted(AbstractAnalysisData *data)
 137 {
 138     impl_->currentLifetimes_.reserve(data->dataSetCount());
 139     impl_->lifetimeHistograms_.reserve(data->dataSetCount());
 140     for (int i = 0; i < data->dataSetCount(); ++i)
 141     {
 142         impl_->currentLifetimes_.push_back(std::vector<int>(data->columnCount(i), 0));
 143         impl_->lifetimeHistograms_.push_back(std::deque<int>());
 144     }
 145 }
 146
 147 void
 148 AnalysisDataLifetimeModule::frameStarted(const AnalysisDataFrameHeader &header)
 149 {
 150     if (header.index() == 0)
 151     {
 152         impl_->firstx_ = header.x();
 153     }
 154     impl_->lastx_ = header.x();
 155     ++impl_->frameCount_;
 156     // TODO: Check the input for even spacing.
 157 }
 158
 159 void
 160 AnalysisDataLifetimeModule::pointsAdded(const AnalysisDataPointSetRef &points)
 161 {
 162     const int dataSet = points.dataSetIndex();
 163     // This assumption is strictly not necessary, but this is how the
 164     // framework works currently, and makes the code below simpler.
 165     GMX_ASSERT(points.firstColumn() == 0
 166                && points.lastColumn() == static_cast<int>(impl_->currentLifetimes_[dataSet].size()) - 1,
 167                "Point set should cover all columns");
 168     for (int i = 0; i < points.columnCount(); ++i)
 169     {
 170         // TODO: Perhaps add control over how this is determined?
 171         const bool bPresent = points.present(i) && points.y(i) > 0.0;
 172         if (bPresent)
 173         {
 174             ++impl_->currentLifetimes_[dataSet][i];
 175         }
 176         else if (impl_->currentLifetimes_[dataSet][i] > 0)
 177         {
 178             impl_->addLifetime(dataSet, impl_->currentLifetimes_[dataSet][i]);
 179             impl_->currentLifetimes_[dataSet][i] = 0;
 180         }
 181     }
 182 }
 183
 184 void
 185 AnalysisDataLifetimeModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
 186 {
 187 }
 188
 189 void
 190 AnalysisDataLifetimeModule::dataFinished()
 191 {
 192     // Need to process the elements present in the last frame explicitly.
 193     for (size_t i = 0; i < impl_->currentLifetimes_.size(); ++i)
 194     {
 195         for (size_t j = 0; j < impl_->currentLifetimes_[i].size(); ++j)
 196         {
 197             impl_->addLifetime(i, impl_->currentLifetimes_[i][j]);
 198         }
 199     }
 200     impl_->currentLifetimes_.clear();
 201
 202     if (impl_->bCumulative_)
 203     {
 204         // Sum up subintervals of longer intervals into the histograms
 205         // if explicitly requested.
 206         std::vector<Impl::LifetimeHistogram>::iterator histogram;
 207         for (histogram  = impl_->lifetimeHistograms_.begin();
 208              histogram != impl_->lifetimeHistograms_.end();
 209              ++histogram)
 210         {
 211             Impl::LifetimeHistogram::iterator shorter, longer;
 212             for (shorter = histogram->begin(); shorter != histogram->end(); ++shorter)
 213             {
 214                 int subIntervalCount = 2;
 215                 for (longer = shorter + 1; longer != histogram->end();
 216                      ++longer, ++subIntervalCount)
 217                 {
 218                     // Interval of length shorter contains (longer - shorter + 1)
 219                     // continuous intervals of length longer.
 220                     *shorter += subIntervalCount * (*longer);
 221                 }
 222             }
 223         }
 224     }
 225
 226     // X spacing is determined by averaging from the first and last frame
 227     // instead of first two frames to avoid rounding issues.
 228     const real spacing =
 229         (impl_->frameCount_ > 1)
 230         ? (impl_->lastx_ - impl_->firstx_) / (impl_->frameCount_ - 1)
 231         : 0.0;
 232     setXAxis(0.0, spacing);
 233
 234     // Determine output dimensionality to cover all the histograms.
 235     setColumnCount(impl_->lifetimeHistograms_.size());
 236     std::vector<Impl::LifetimeHistogram>::const_iterator histogram;
 237     size_t maxLifetime = 1;
 238     for (histogram  = impl_->lifetimeHistograms_.begin();
 239          histogram != impl_->lifetimeHistograms_.end();
 240          ++histogram)
 241     {
 242         maxLifetime = std::max(maxLifetime, histogram->size());
 243     }
 244     setRowCount(maxLifetime);
 245
 246     // Fill up the output data from the histograms.
 247     allocateValues();
 248     int column = 0;
 249     for (histogram  = impl_->lifetimeHistograms_.begin();
 250          histogram != impl_->lifetimeHistograms_.end();
 251          ++histogram, ++column)
 252     {
 253         int row = 0;
 254         Impl::LifetimeHistogram::const_iterator i;
 255         for (i = histogram->begin(); i != histogram->end(); ++i, ++row)
 256         {
 257             // Normalize by the number of frames, taking into account the
 258             // length of the interval (interval of length N cannot start in
 259             // N-1 last frames).  row is always smaller than frameCount_
 260             // because of the histograms have at most frameCount_ entries.
 261             const real normalized = *i / static_cast<real>(impl_->frameCount_ - row);
 262             value(row, column).setValue(normalized);
 263         }
 264         // Pad the rest of the histogram with zeros to match the longest
 265         // histogram.
 266         for (; row < rowCount(); ++row)
 267         {
 268             value(row, column).setValue(0.0);
 269         }
 270     }
 271     impl_->lifetimeHistograms_.clear();
 272     valuesReady();
 273 }
 274
 275 } // namespace gmx