2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
38 * Implements gmx::AnalysisDataLifetimeModule.
40 * \author Teemu Murtola <teemu.murtola@gmail.com>
41 * \ingroup module_analysisdata
53 #include "gromacs/analysisdata/dataframe.h"
54 #include "gromacs/analysisdata/datastorage.h"
59 /********************************************************************
60 * AnalysisDataLifetimeModule
64 * Private implementation class for AnalysisDataLifetimeModule.
66 * \ingroup module_analysisdata
68 class AnalysisDataLifetimeModule::Impl
71 //! Container type for storing a histogram during the calculation.
72 typedef std::deque
<int> LifetimeHistogram
;
74 //! Initializes the implementation class with empty/default values.
75 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() : impl_(new Impl()) {}
120 AnalysisDataLifetimeModule::~AnalysisDataLifetimeModule() {}
122 void AnalysisDataLifetimeModule::setCumulative(bool bCumulative
)
124 impl_
->bCumulative_
= bCumulative
;
127 int AnalysisDataLifetimeModule::flags() const
129 return efAllowMulticolumn
| efAllowMissing
| efAllowMultipleDataSets
;
132 void AnalysisDataLifetimeModule::dataStarted(AbstractAnalysisData
* data
)
134 impl_
->currentLifetimes_
.reserve(data
->dataSetCount());
135 impl_
->lifetimeHistograms_
.reserve(data
->dataSetCount());
136 for (int i
= 0; i
< data
->dataSetCount(); ++i
)
138 impl_
->currentLifetimes_
.emplace_back(data
->columnCount(i
), 0);
139 impl_
->lifetimeHistograms_
.emplace_back();
143 void AnalysisDataLifetimeModule::frameStarted(const AnalysisDataFrameHeader
& header
)
145 if (header
.index() == 0)
147 impl_
->firstx_
= header
.x();
149 impl_
->lastx_
= header
.x();
150 ++impl_
->frameCount_
;
151 // TODO: Check the input for even spacing.
154 void AnalysisDataLifetimeModule::pointsAdded(const AnalysisDataPointSetRef
& points
)
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
)
164 // TODO: Perhaps add control over how this is determined?
165 const bool bPresent
= points
.present(i
) && points
.y(i
) > 0.0;
168 ++impl_
->currentLifetimes_
[dataSet
][i
];
170 else if (impl_
->currentLifetimes_
[dataSet
][i
] > 0)
172 impl_
->addLifetime(dataSet
, impl_
->currentLifetimes_
[dataSet
][i
]);
173 impl_
->currentLifetimes_
[dataSet
][i
] = 0;
178 void AnalysisDataLifetimeModule::frameFinished(const AnalysisDataFrameHeader
& /*header*/) {}
180 void AnalysisDataLifetimeModule::dataFinished()
182 // Need to process the elements present in the last frame explicitly.
183 for (size_t i
= 0; i
< impl_
->currentLifetimes_
.size(); ++i
)
185 for (size_t j
= 0; j
< impl_
->currentLifetimes_
[i
].size(); ++j
)
187 impl_
->addLifetime(i
, impl_
->currentLifetimes_
[i
][j
]);
190 impl_
->currentLifetimes_
.clear();
192 if (impl_
->bCumulative_
)
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
)
200 Impl::LifetimeHistogram::iterator shorter
, longer
;
201 for (shorter
= histogram
->begin(); shorter
!= histogram
->end(); ++shorter
)
203 int subIntervalCount
= 2;
204 for (longer
= shorter
+ 1; longer
!= histogram
->end(); ++longer
, ++subIntervalCount
)
206 // Interval of length shorter contains (longer - shorter + 1)
207 // continuous intervals of length longer.
208 *shorter
+= subIntervalCount
* (*longer
);
214 // X spacing is determined by averaging from the first and last frame
215 // instead of first two frames to avoid rounding issues.
217 (impl_
->frameCount_
> 1) ? (impl_
->lastx_
- impl_
->firstx_
) / (impl_
->frameCount_
- 1) : 0.0;
218 setXAxis(0.0, spacing
);
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
)
227 maxLifetime
= std::max(maxLifetime
, histogram
->size());
229 setRowCount(maxLifetime
);
231 // Fill up the output data from the histograms.
234 for (histogram
= impl_
->lifetimeHistograms_
.begin();
235 histogram
!= impl_
->lifetimeHistograms_
.end(); ++histogram
, ++column
)
238 Impl::LifetimeHistogram::const_iterator i
;
239 for (i
= histogram
->begin(); i
!= histogram
->end(); ++i
, ++row
)
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
);
248 // Pad the rest of the histogram with zeros to match the longest
250 for (; row
< rowCount(); ++row
)
252 value(row
, column
).setValue(0.0);
255 impl_
->lifetimeHistograms_
.clear();