2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
37 * Implements gmx::AnalysisDataLifetimeModule.
39 * \author Teemu Murtola <teemu.murtola@gmail.com>
40 * \ingroup module_analysisdata
51 #include "gromacs/analysisdata/dataframe.h"
52 #include "gromacs/analysisdata/datastorage.h"
57 /********************************************************************
58 * AnalysisDataLifetimeModule
62 * Private implementation class for AnalysisDataLifetimeModule.
64 * \ingroup module_analysisdata
66 class AnalysisDataLifetimeModule::Impl
69 //! Container type for storing a histogram during the calculation.
70 typedef std::deque<int> LifetimeHistogram;
72 //! Initializes the implementation class with empty/default values.
73 Impl() : firstx_(0.0), lastx_(0.0), frameCount_(0), bCumulative_(false)
78 * Increments a lifetime histogram with a single lifetime.
80 * \param[in] dataSet Index of the histogram to increment.
81 * \param[in] lifetime Lifetime to add to the histogram.
83 void addLifetime(int dataSet, int lifetime)
87 LifetimeHistogram &histogram = lifetimeHistograms_[dataSet];
88 if (histogram.size() < static_cast<unsigned>(lifetime))
90 histogram.resize(lifetime, 0);
92 ++histogram[lifetime - 1];
96 //! X value of the first frame (used for determining output spacing).
98 //! X value of the last frame (used for determining output spacing).
100 //! Total number of frames (used for normalization and output spacing).
102 //! Whether to add subintervals of longer intervals explicitly.
105 * Length of current continuously present interval for each data column.
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.
111 std::vector<std::vector<int> > currentLifetimes_;
113 * Accumulated lifetime histograms for each data set.
115 std::vector<LifetimeHistogram> lifetimeHistograms_;
118 AnalysisDataLifetimeModule::AnalysisDataLifetimeModule()
123 AnalysisDataLifetimeModule::~AnalysisDataLifetimeModule()
127 void AnalysisDataLifetimeModule::setCumulative(bool bCumulative)
129 impl_->bCumulative_ = bCumulative;
132 int AnalysisDataLifetimeModule::flags() const
134 return efAllowMulticolumn | efAllowMissing | efAllowMultipleDataSets;
138 AnalysisDataLifetimeModule::dataStarted(AbstractAnalysisData *data)
140 impl_->currentLifetimes_.reserve(data->dataSetCount());
141 impl_->lifetimeHistograms_.reserve(data->dataSetCount());
142 for (int i = 0; i < data->dataSetCount(); ++i)
144 impl_->currentLifetimes_.push_back(std::vector<int>(data->columnCount(i), 0));
145 impl_->lifetimeHistograms_.push_back(std::deque<int>());
150 AnalysisDataLifetimeModule::frameStarted(const AnalysisDataFrameHeader &header)
152 if (header.index() == 0)
154 impl_->firstx_ = header.x();
156 impl_->lastx_ = header.x();
157 ++impl_->frameCount_;
158 // TODO: Check the input for even spacing.
162 AnalysisDataLifetimeModule::pointsAdded(const AnalysisDataPointSetRef &points)
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)
172 // TODO: Perhaps add control over how this is determined?
173 const bool bPresent = points.present(i) && points.y(i) > 0.0;
176 ++impl_->currentLifetimes_[dataSet][i];
178 else if (impl_->currentLifetimes_[dataSet][i] > 0)
180 impl_->addLifetime(dataSet, impl_->currentLifetimes_[dataSet][i]);
181 impl_->currentLifetimes_[dataSet][i] = 0;
187 AnalysisDataLifetimeModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
192 AnalysisDataLifetimeModule::dataFinished()
194 // Need to process the elements present in the last frame explicitly.
195 for (size_t i = 0; i < impl_->currentLifetimes_.size(); ++i)
197 for (size_t j = 0; j < impl_->currentLifetimes_[i].size(); ++j)
199 impl_->addLifetime(i, impl_->currentLifetimes_[i][j]);
202 impl_->currentLifetimes_.clear();
204 if (impl_->bCumulative_)
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();
213 Impl::LifetimeHistogram::iterator shorter, longer;
214 for (shorter = histogram->begin(); shorter != histogram->end(); ++shorter)
216 int subIntervalCount = 2;
217 for (longer = shorter + 1; longer != histogram->end();
218 ++longer, ++subIntervalCount)
220 // Interval of length shorter contains (longer - shorter + 1)
221 // continuous intervals of length longer.
222 *shorter += subIntervalCount * (*longer);
228 // X spacing is determined by averaging from the first and last frame
229 // instead of first two frames to avoid rounding issues.
231 (impl_->frameCount_ > 1)
232 ? (impl_->lastx_ - impl_->firstx_) / (impl_->frameCount_ - 1)
234 setXAxis(0.0, spacing);
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();
244 maxLifetime = std::max(maxLifetime, histogram->size());
246 setRowCount(maxLifetime);
248 // Fill up the output data from the histograms.
251 for (histogram = impl_->lifetimeHistograms_.begin();
252 histogram != impl_->lifetimeHistograms_.end();
253 ++histogram, ++column)
256 Impl::LifetimeHistogram::const_iterator i;
257 for (i = histogram->begin(); i != histogram->end(); ++i, ++row)
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);
266 // Pad the rest of the histogram with zeros to match the longest
268 for (; row < rowCount(); ++row)
270 value(row, column).setValue(0.0);
273 impl_->lifetimeHistograms_.clear();