2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2010,2011,2012,2013, by the GROMACS development team, led by
5 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
6 * others, as listed in the AUTHORS file in the top-level source
7 * 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 classes in histogram.h.
39 * \author Teemu Murtola <teemu.murtola@gmail.com>
40 * \ingroup module_analysisdata
42 #include "gromacs/analysisdata/modules/histogram.h"
48 #include "gromacs/analysisdata/dataframe.h"
49 #include "gromacs/analysisdata/datastorage.h"
50 #include "gromacs/utility/exceptions.h"
51 #include "gromacs/utility/gmxassert.h"
53 #include "frameaverager.h"
58 //! Value used to signify that a real-valued histogram setting is not set.
59 const real UNDEFINED = std::numeric_limits<real>::max();
60 //! Checks whether \p value is defined.
61 bool isDefined(real value)
63 return value != UNDEFINED;
71 /********************************************************************
72 * AnalysisHistogramSettingsInitializer
75 AnalysisHistogramSettingsInitializer::AnalysisHistogramSettingsInitializer()
76 : min_(UNDEFINED), max_(UNDEFINED), binWidth_(UNDEFINED),
77 binCount_(0), bIntegerBins_(false), bRoundRange_(false),
83 /********************************************************************
84 * AnalysisHistogramSettings
87 AnalysisHistogramSettings::AnalysisHistogramSettings()
88 : firstEdge_(0.0), lastEdge_(0.0), binWidth_(0.0), inverseBinWidth_(0.0),
89 binCount_(0), bAll_(false)
94 AnalysisHistogramSettings::AnalysisHistogramSettings(
95 const AnalysisHistogramSettingsInitializer &settings)
97 GMX_RELEASE_ASSERT(isDefined(settings.min_),
98 "Histogram start value must be defined");
99 GMX_RELEASE_ASSERT(!isDefined(settings.max_) || settings.max_ > settings.min_,
100 "Histogram end value must be larger than start value");
101 GMX_RELEASE_ASSERT(!isDefined(settings.binWidth_) || settings.binWidth_ > 0.0,
102 "Histogram bin width must be positive");
103 GMX_RELEASE_ASSERT(settings.binCount_ >= 0,
104 "Histogram bin count must be positive");
106 if (!isDefined(settings.max_))
108 GMX_RELEASE_ASSERT(isDefined(settings.binWidth_) && settings.binCount_ > 0,
109 "Not all required values provided");
110 GMX_RELEASE_ASSERT(!settings.bRoundRange_,
111 "Rounding only supported for min/max ranges");
113 firstEdge_ = settings.min_;
114 binCount_ = settings.binCount_;
115 binWidth_ = settings.binWidth_;
116 if (settings.bIntegerBins_)
118 firstEdge_ -= 0.5 * binWidth_;
120 lastEdge_ = firstEdge_ + binCount_ * binWidth_;
124 GMX_RELEASE_ASSERT(!(isDefined(settings.binWidth_) && settings.binCount_ > 0),
125 "Conflicting histogram bin specifications");
126 GMX_RELEASE_ASSERT(isDefined(settings.binWidth_) || settings.binCount_ > 0,
127 "Not all required values provided");
129 if (settings.bRoundRange_)
131 GMX_RELEASE_ASSERT(!settings.bIntegerBins_,
132 "Rounding and integer bins cannot be combined");
133 GMX_RELEASE_ASSERT(isDefined(settings.binWidth_),
134 "Rounding only makes sense with defined binwidth");
135 binWidth_ = settings.binWidth_;
136 firstEdge_ = binWidth_ * floor(settings.min_ / binWidth_);
137 lastEdge_ = binWidth_ * ceil(settings.max_ / binWidth_);
138 binCount_ = static_cast<int>((lastEdge_ - firstEdge_) / binWidth_ + 0.5);
142 firstEdge_ = settings.min_;
143 lastEdge_ = settings.max_;
144 if (settings.binCount_ > 0)
146 binCount_ = settings.binCount_;
147 if (settings.bIntegerBins_)
149 GMX_RELEASE_ASSERT(settings.binCount_ > 1,
150 "Bin count must be at least two with integer bins");
151 binWidth_ = (lastEdge_ - firstEdge_) / (binCount_ - 1);
152 firstEdge_ -= 0.5 * binWidth_;
153 lastEdge_ += 0.5 * binWidth_;
157 binWidth_ = (lastEdge_ - firstEdge_) / binCount_;
162 binWidth_ = settings.binWidth_;
163 binCount_ = static_cast<int>((lastEdge_ - firstEdge_) / binWidth_ + 0.5);
164 if (settings.bIntegerBins_)
166 firstEdge_ -= 0.5 * binWidth_;
169 lastEdge_ = firstEdge_ + binCount_ * binWidth_;
174 inverseBinWidth_ = 1.0 / binWidth_;
175 bAll_ = settings.bIncludeAll_;
180 AnalysisHistogramSettings::findBin(real y) const
184 return bAll_ ? 0 : -1;
186 int bin = static_cast<int>((y - firstEdge_) * inverseBinWidth_);
187 if (bin >= binCount_)
189 return bAll_ ? binCount_ - 1 : -1;
195 /********************************************************************
196 * StaticAverageHistogram
203 * Represents copies of average histograms.
205 * Methods in AbstractAverageHistogram that return new histogram instances
206 * return objects of this class.
207 * Initialization of values is handled in those methods.
209 * \ingroup module_analysisdata
211 class StaticAverageHistogram : public AbstractAverageHistogram
214 StaticAverageHistogram();
215 //! Creates an average histogram module with defined bin parameters.
216 explicit StaticAverageHistogram(const AnalysisHistogramSettings &settings);
218 // Copy and assign disallowed by base.
221 StaticAverageHistogram::StaticAverageHistogram()
226 StaticAverageHistogram::StaticAverageHistogram(
227 const AnalysisHistogramSettings &settings)
228 : AbstractAverageHistogram(settings)
235 /********************************************************************
236 * AbstractAverageHistogram
239 AbstractAverageHistogram::AbstractAverageHistogram()
244 AbstractAverageHistogram::AbstractAverageHistogram(
245 const AnalysisHistogramSettings &settings)
246 : settings_(settings)
248 setRowCount(settings.binCount());
249 setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
250 settings.binWidth());
254 AbstractAverageHistogram::~AbstractAverageHistogram()
260 AbstractAverageHistogram::init(const AnalysisHistogramSettings &settings)
262 settings_ = settings;
263 setRowCount(settings.binCount());
264 setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
265 settings.binWidth());
269 AverageHistogramPointer
270 AbstractAverageHistogram::resampleDoubleBinWidth(bool bIntegerBins) const
275 nbins = (rowCount() + 1) / 2;
279 nbins = rowCount() / 2;
282 AverageHistogramPointer dest(
283 new StaticAverageHistogram(
284 histogramFromBins(xstart(), nbins, 2*xstep())
285 .integerBins(bIntegerBins)));
286 dest->setColumnCount(columnCount());
287 dest->allocateValues();
290 for (i = j = 0; i < nbins; ++i)
292 const bool bFirstHalfBin = (bIntegerBins && i == 0);
293 for (int c = 0; c < columnCount(); ++c)
304 v2 = value(j + 1, c);
308 dest->setValue(i, c, sqrt(v1 * v1 + v2 * v2));
312 dest->setValue(i, c, v1 + v2);
328 AverageHistogramPointer
329 AbstractAverageHistogram::clone() const
331 AverageHistogramPointer dest(new StaticAverageHistogram());
332 copyContents(this, dest.get());
338 AbstractAverageHistogram::normalizeProbability()
341 for (int i = 0; i < rowCount(); ++i)
345 scale(1.0 / (sum * xstep()));
350 AbstractAverageHistogram::scale(real norm)
352 for (int i = 0; i < rowCount(); ++i)
361 AbstractAverageHistogram::scaleVector(real norm[])
363 for (int i = 0; i < rowCount(); ++i)
365 value(i, 0) *= norm[i];
366 value(i, 1) *= norm[i];
371 /********************************************************************
372 * BasicAverageHistogramModule
379 * Implements average histogram module that averages per-frame histograms.
381 * This class is used for accumulating average histograms in per-frame
382 * histogram modules (those that use BasicHistogramImpl as their implementation
384 * There are two columns, first for the average and second for standard
387 * \ingroup module_analysisdata
389 class BasicAverageHistogramModule : public AbstractAverageHistogram,
390 public AnalysisDataModuleInterface
393 BasicAverageHistogramModule();
394 //! Creates an average histogram module with defined bin parameters.
395 explicit BasicAverageHistogramModule(const AnalysisHistogramSettings &settings);
397 using AbstractAverageHistogram::init;
399 virtual int flags() const;
401 virtual void dataStarted(AbstractAnalysisData *data);
402 virtual void frameStarted(const AnalysisDataFrameHeader &header);
403 virtual void pointsAdded(const AnalysisDataPointSetRef &points);
404 virtual void frameFinished(const AnalysisDataFrameHeader &header);
405 virtual void dataFinished();
408 //! Averaging helper object.
409 AnalysisDataFrameAverager averager_;
411 // Copy and assign disallowed by base.
414 BasicAverageHistogramModule::BasicAverageHistogramModule()
420 BasicAverageHistogramModule::BasicAverageHistogramModule(
421 const AnalysisHistogramSettings &settings)
422 : AbstractAverageHistogram(settings)
429 BasicAverageHistogramModule::flags() const
431 return efAllowMulticolumn;
436 BasicAverageHistogramModule::dataStarted(AbstractAnalysisData *data)
438 GMX_RELEASE_ASSERT(rowCount() == data->columnCount(),
439 "Inconsistent data sizes, something is wrong in the initialization");
440 averager_.setColumnCount(rowCount());
445 BasicAverageHistogramModule::frameStarted(const AnalysisDataFrameHeader & /*header*/)
451 BasicAverageHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
453 averager_.addPoints(points);
458 BasicAverageHistogramModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
464 BasicAverageHistogramModule::dataFinished()
468 for (int i = 0; i < rowCount(); ++i)
470 setValue(i, 0, averager_.average(i));
471 setValue(i, 1, sqrt(averager_.variance(i)));
476 /********************************************************************
481 * Private implementation class for AnalysisDataSimpleHistogramModule and
482 * AnalysisDataWeightedHistogramModule.
484 * \ingroup module_analysisdata
486 class BasicHistogramImpl
489 //! Smart pointer to manage an BasicAverageHistogramModule object.
490 typedef boost::shared_ptr<BasicAverageHistogramModule>
491 BasicAverageHistogramModulePointer;
493 BasicHistogramImpl();
494 //! Creates an histogram impl with defined bin parameters.
495 explicit BasicHistogramImpl(const AnalysisHistogramSettings &settings);
496 ~BasicHistogramImpl();
499 * (Re)initializes the histogram from settings.
501 void init(const AnalysisHistogramSettings &settings);
503 * Initializes data storage frame when a new frame starts.
505 void initFrame(AnalysisDataStorageFrame *frame);
507 //! Storage implementation object.
508 AnalysisDataStorage storage_;
509 //! Settings for the histogram object.
510 AnalysisHistogramSettings settings_;
512 BasicAverageHistogramModulePointer averager_;
515 BasicHistogramImpl::BasicHistogramImpl()
516 : averager_(new BasicAverageHistogramModule())
521 BasicHistogramImpl::BasicHistogramImpl(const AnalysisHistogramSettings &settings)
522 : settings_(settings), averager_(new BasicAverageHistogramModule(settings))
527 BasicHistogramImpl::~BasicHistogramImpl()
532 void BasicHistogramImpl::init(const AnalysisHistogramSettings &settings)
534 settings_ = settings;
535 averager_->init(settings);
540 BasicHistogramImpl::initFrame(AnalysisDataStorageFrame *frame)
542 for (int i = 0; i < frame->columnCount(); ++i)
544 frame->setValue(i, 0.0);
548 } // namespace internal
551 /********************************************************************
552 * AnalysisDataSimpleHistogramModule
555 AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule()
556 : impl_(new internal::BasicHistogramImpl())
561 AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule(
562 const AnalysisHistogramSettings &settings)
563 : impl_(new internal::BasicHistogramImpl(settings))
568 AnalysisDataSimpleHistogramModule::~AnalysisDataSimpleHistogramModule()
573 void AnalysisDataSimpleHistogramModule::init(const AnalysisHistogramSettings &settings)
575 impl_->init(settings);
579 AbstractAverageHistogram &
580 AnalysisDataSimpleHistogramModule::averager()
582 return *impl_->averager_;
586 const AnalysisHistogramSettings &
587 AnalysisDataSimpleHistogramModule::settings() const
589 return impl_->settings_;
594 AnalysisDataSimpleHistogramModule::flags() const
596 return efAllowMulticolumn | efAllowMultipoint;
601 AnalysisDataSimpleHistogramModule::dataStarted(AbstractAnalysisData *data)
603 addModule(impl_->averager_);
604 setColumnCount(settings().binCount());
606 impl_->storage_.startDataStorage(this);
611 AnalysisDataSimpleHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
613 AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
614 impl_->initFrame(&frame);
619 AnalysisDataSimpleHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
621 AnalysisDataStorageFrame &frame =
622 impl_->storage_.currentFrame(points.frameIndex());
623 for (int i = 0; i < points.columnCount(); ++i)
625 int bin = settings().findBin(points.y(i));
628 frame.value(bin) += 1;
635 AnalysisDataSimpleHistogramModule::frameFinished(const AnalysisDataFrameHeader &header)
637 impl_->storage_.finishFrame(header.index());
642 AnalysisDataSimpleHistogramModule::dataFinished()
649 AnalysisDataSimpleHistogramModule::tryGetDataFrameInternal(int index) const
651 return impl_->storage_.tryGetDataFrame(index);
656 AnalysisDataSimpleHistogramModule::requestStorageInternal(int nframes)
658 return impl_->storage_.requestStorage(nframes);
662 /********************************************************************
663 * AnalysisDataWeightedHistogramModule
666 AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule()
667 : impl_(new internal::BasicHistogramImpl())
672 AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule(
673 const AnalysisHistogramSettings &settings)
674 : impl_(new internal::BasicHistogramImpl(settings))
679 AnalysisDataWeightedHistogramModule::~AnalysisDataWeightedHistogramModule()
684 void AnalysisDataWeightedHistogramModule::init(const AnalysisHistogramSettings &settings)
686 impl_->init(settings);
690 AbstractAverageHistogram &
691 AnalysisDataWeightedHistogramModule::averager()
693 return *impl_->averager_;
697 const AnalysisHistogramSettings &
698 AnalysisDataWeightedHistogramModule::settings() const
700 return impl_->settings_;
705 AnalysisDataWeightedHistogramModule::flags() const
707 return efAllowMulticolumn | efAllowMultipoint;
712 AnalysisDataWeightedHistogramModule::dataStarted(AbstractAnalysisData *data)
714 addModule(impl_->averager_);
715 setColumnCount(settings().binCount());
717 impl_->storage_.startDataStorage(this);
722 AnalysisDataWeightedHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
724 AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
725 impl_->initFrame(&frame);
730 AnalysisDataWeightedHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
732 if (points.firstColumn() != 0 || points.columnCount() < 2)
734 GMX_THROW(APIError("Invalid data layout"));
736 int bin = settings().findBin(points.y(0));
739 AnalysisDataStorageFrame &frame =
740 impl_->storage_.currentFrame(points.frameIndex());
741 for (int i = 1; i < points.columnCount(); ++i)
743 frame.value(bin) += points.y(i);
750 AnalysisDataWeightedHistogramModule::frameFinished(const AnalysisDataFrameHeader &header)
752 impl_->storage_.finishFrame(header.index());
757 AnalysisDataWeightedHistogramModule::dataFinished()
764 AnalysisDataWeightedHistogramModule::tryGetDataFrameInternal(int index) const
766 return impl_->storage_.tryGetDataFrame(index);
771 AnalysisDataWeightedHistogramModule::requestStorageInternal(int nframes)
773 return impl_->storage_.requestStorage(nframes);
777 /********************************************************************
778 * AnalysisDataBinAverageModule
781 class AnalysisDataBinAverageModule::Impl
785 explicit Impl(const AnalysisHistogramSettings &settings)
786 : settings_(settings)
790 //! Histogram settings.
791 AnalysisHistogramSettings settings_;
792 //! Averaging helper object.
793 AnalysisDataFrameAverager averager_;
796 AnalysisDataBinAverageModule::AnalysisDataBinAverageModule()
803 AnalysisDataBinAverageModule::AnalysisDataBinAverageModule(
804 const AnalysisHistogramSettings &settings)
805 : impl_(new Impl(settings))
808 setRowCount(settings.binCount());
809 setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
810 settings.binWidth());
814 AnalysisDataBinAverageModule::~AnalysisDataBinAverageModule()
820 AnalysisDataBinAverageModule::init(const AnalysisHistogramSettings &settings)
822 impl_->settings_ = settings;
823 setRowCount(settings.binCount());
824 setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
825 settings.binWidth());
829 const AnalysisHistogramSettings &
830 AnalysisDataBinAverageModule::settings() const
832 return impl_->settings_;
837 AnalysisDataBinAverageModule::flags() const
839 return efAllowMulticolumn | efAllowMultipoint;
844 AnalysisDataBinAverageModule::dataStarted(AbstractAnalysisData * /*data*/)
846 impl_->averager_.setColumnCount(rowCount());
851 AnalysisDataBinAverageModule::frameStarted(const AnalysisDataFrameHeader & /*header*/)
857 AnalysisDataBinAverageModule::pointsAdded(const AnalysisDataPointSetRef &points)
859 if (points.firstColumn() != 0 || points.columnCount() < 2)
861 GMX_THROW(APIError("Invalid data layout"));
863 int bin = settings().findBin(points.y(0));
866 for (int i = 1; i < points.columnCount(); ++i)
868 impl_->averager_.addValue(bin, points.y(i));
875 AnalysisDataBinAverageModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
881 AnalysisDataBinAverageModule::dataFinished()
883 impl_->averager_.finish();
885 for (int i = 0; i < settings().binCount(); ++i)
887 setValue(i, 0, impl_->averager_.average(i));
888 setValue(i, 1, std::sqrt(impl_->averager_.variance(i)));
889 setValue(i, 2, impl_->averager_.sampleCount(i));