src/gromacs/analysisdata/modules/lifetime.cpp

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