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 by the GROMACS development team.
   5  * Copyright (c) 2019,2020, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 /*! \internal \file
  37  * \brief
  38  * Implements gmx::AnalysisDataLifetimeModule.
  39  *
  40  * \author Teemu Murtola <teemu.murtola@gmail.com>
  41  * \ingroup module_analysisdata
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "lifetime.h"
  46
  47 #include <cmath>
  48
  49 #include <algorithm>
  50 #include <deque>
  51 #include <vector>
  52
  53 #include "gromacs/analysisdata/dataframe.h"
  54 #include "gromacs/analysisdata/datastorage.h"
  55
  56 namespace gmx
  57 {
  58
  59 /********************************************************************
  60  * AnalysisDataLifetimeModule
  61  */
  62
  63 /*! \internal \brief
  64  * Private implementation class for AnalysisDataLifetimeModule.
  65  *
  66  * \ingroup module_analysisdata
  67  */
  68 class AnalysisDataLifetimeModule::Impl
  69 {
  70 public:
  71     //! Container type for storing a histogram during the calculation.
  72     typedef std::deque<int> LifetimeHistogram;
  73
  74     //! Initializes the implementation class with empty/default values.
  75     Impl() : firstx_(0.0), lastx_(0.0), frameCount_(0), bCumulative_(false) {}
  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() : impl_(new Impl()) {}
 119
 120 AnalysisDataLifetimeModule::~AnalysisDataLifetimeModule() {}
 121
 122 void AnalysisDataLifetimeModule::setCumulative(bool bCumulative)
 123 {
 124     impl_->bCumulative_ = bCumulative;
 125 }
 126
 127 int AnalysisDataLifetimeModule::flags() const
 128 {
 129     return efAllowMulticolumn | efAllowMissing | efAllowMultipleDataSets;
 130 }
 131
 132 void AnalysisDataLifetimeModule::dataStarted(AbstractAnalysisData* data)
 133 {
 134     impl_->currentLifetimes_.reserve(data->dataSetCount());
 135     impl_->lifetimeHistograms_.reserve(data->dataSetCount());
 136     for (int i = 0; i < data->dataSetCount(); ++i)
 137     {
 138         impl_->currentLifetimes_.emplace_back(data->columnCount(i), 0);
 139         impl_->lifetimeHistograms_.emplace_back();
 140     }
 141 }
 142
 143 void AnalysisDataLifetimeModule::frameStarted(const AnalysisDataFrameHeader& header)
 144 {
 145     if (header.index() == 0)
 146     {
 147         impl_->firstx_ = header.x();
 148     }
 149     impl_->lastx_ = header.x();
 150     ++impl_->frameCount_;
 151     // TODO: Check the input for even spacing.
 152 }
 153
 154 void AnalysisDataLifetimeModule::pointsAdded(const AnalysisDataPointSetRef& points)
 155 {
 156     const int dataSet = points.dataSetIndex();
 157     // This assumption is strictly not necessary, but this is how the
 158     // framework works currently, and makes the code below simpler.
 159     GMX_ASSERT(points.firstColumn() == 0
 160                        && points.lastColumn() == ssize(impl_->currentLifetimes_[dataSet]) - 1,
 161                "Point set should cover all columns");
 162     for (int i = 0; i < points.columnCount(); ++i)
 163     {
 164         // TODO: Perhaps add control over how this is determined?
 165         const bool bPresent = points.present(i) && points.y(i) > 0.0;
 166         if (bPresent)
 167         {
 168             ++impl_->currentLifetimes_[dataSet][i];
 169         }
 170         else if (impl_->currentLifetimes_[dataSet][i] > 0)
 171         {
 172             impl_->addLifetime(dataSet, impl_->currentLifetimes_[dataSet][i]);
 173             impl_->currentLifetimes_[dataSet][i] = 0;
 174         }
 175     }
 176 }
 177
 178 void AnalysisDataLifetimeModule::frameFinished(const AnalysisDataFrameHeader& /*header*/) {}
 179
 180 void AnalysisDataLifetimeModule::dataFinished()
 181 {
 182     // Need to process the elements present in the last frame explicitly.
 183     for (size_t i = 0; i < impl_->currentLifetimes_.size(); ++i)
 184     {
 185         for (size_t j = 0; j < impl_->currentLifetimes_[i].size(); ++j)
 186         {
 187             impl_->addLifetime(i, impl_->currentLifetimes_[i][j]);
 188         }
 189     }
 190     impl_->currentLifetimes_.clear();
 191
 192     if (impl_->bCumulative_)
 193     {
 194         // Sum up subintervals of longer intervals into the histograms
 195         // if explicitly requested.
 196         std::vector<Impl::LifetimeHistogram>::iterator histogram;
 197         for (histogram = impl_->lifetimeHistograms_.begin();
 198              histogram != impl_->lifetimeHistograms_.end(); ++histogram)
 199         {
 200             Impl::LifetimeHistogram::iterator shorter, longer;
 201             for (shorter = histogram->begin(); shorter != histogram->end(); ++shorter)
 202             {
 203                 int subIntervalCount = 2;
 204                 for (longer = shorter + 1; longer != histogram->end(); ++longer, ++subIntervalCount)
 205                 {
 206                     // Interval of length shorter contains (longer - shorter + 1)
 207                     // continuous intervals of length longer.
 208                     *shorter += subIntervalCount * (*longer);
 209                 }
 210             }
 211         }
 212     }
 213
 214     // X spacing is determined by averaging from the first and last frame
 215     // instead of first two frames to avoid rounding issues.
 216     const real spacing =
 217             (impl_->frameCount_ > 1) ? (impl_->lastx_ - impl_->firstx_) / (impl_->frameCount_ - 1) : 0.0;
 218     setXAxis(0.0, spacing);
 219
 220     // Determine output dimensionality to cover all the histograms.
 221     setColumnCount(impl_->lifetimeHistograms_.size());
 222     std::vector<Impl::LifetimeHistogram>::const_iterator histogram;
 223     size_t                                               maxLifetime = 1;
 224     for (histogram = impl_->lifetimeHistograms_.begin();
 225          histogram != impl_->lifetimeHistograms_.end(); ++histogram)
 226     {
 227         maxLifetime = std::max(maxLifetime, histogram->size());
 228     }
 229     setRowCount(maxLifetime);
 230
 231     // Fill up the output data from the histograms.
 232     allocateValues();
 233     int column = 0;
 234     for (histogram = impl_->lifetimeHistograms_.begin();
 235          histogram != impl_->lifetimeHistograms_.end(); ++histogram, ++column)
 236     {
 237         int                                     row = 0;
 238         Impl::LifetimeHistogram::const_iterator i;
 239         for (i = histogram->begin(); i != histogram->end(); ++i, ++row)
 240         {
 241             // Normalize by the number of frames, taking into account the
 242             // length of the interval (interval of length N cannot start in
 243             // N-1 last frames).  row is always smaller than frameCount_
 244             // because of the histograms have at most frameCount_ entries.
 245             const real normalized = *i / static_cast<real>(impl_->frameCount_ - row);
 246             value(row, column).setValue(normalized);
 247         }
 248         // Pad the rest of the histogram with zeros to match the longest
 249         // histogram.
 250         for (; row < rowCount(); ++row)
 251         {
 252             value(row, column).setValue(0.0);
 253         }
 254     }
 255     impl_->lifetimeHistograms_.clear();
 256     valuesReady();
 257 }
 258
 259 } // namespace gmx