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