src/gromacs/analysisdata/modules/lifetime.cpp

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