3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
9 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
10 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
11 * Copyright (c) 2001-2009, The GROMACS development team,
12 * check out http://www.gromacs.org for more information.
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
19 * If you want to redistribute modifications, please consider that
20 * scientific software is very special. Version control is crucial -
21 * bugs must be traceable. We will be happy to consider code for
22 * inclusion in the official distribution, but derived work must not
23 * be called official GROMACS. Details are found in the README & COPYING
24 * files - if they are missing, get the official version at www.gromacs.org.
26 * To help us fund GROMACS development, we humbly ask that you cite
27 * the papers on the package - you can find them in the top README file.
29 * For more info, check our website at http://www.gromacs.org
33 * Implements classes in histogram.h.
35 * \author Teemu Murtola <teemu.murtola@cbr.su.se>
36 * \ingroup module_analysisdata
38 #include "gromacs/analysisdata/modules/histogram.h"
44 #include "gromacs/analysisdata/dataframe.h"
45 #include "gromacs/analysisdata/datastorage.h"
46 #include "gromacs/utility/exceptions.h"
47 #include "gromacs/utility/gmxassert.h"
49 static const real UNDEFINED = std::numeric_limits<real>::max();
50 static bool isDefined(real value)
52 return value != UNDEFINED;
58 /********************************************************************
59 * AnalysisHistogramSettingsInitializer
62 AnalysisHistogramSettingsInitializer::AnalysisHistogramSettingsInitializer()
63 : min_(UNDEFINED), max_(UNDEFINED), binWidth_(UNDEFINED),
64 binCount_(0), bIntegerBins_(false), bRoundRange_(false),
70 /********************************************************************
71 * AnalysisHistogramSettings
74 AnalysisHistogramSettings::AnalysisHistogramSettings()
75 : firstEdge_(0.0), lastEdge_(0.0), binWidth_(0.0), inverseBinWidth_(0.0),
76 binCount_(0), bAll_(false)
81 AnalysisHistogramSettings::AnalysisHistogramSettings(
82 const AnalysisHistogramSettingsInitializer &settings)
84 GMX_RELEASE_ASSERT(isDefined(settings.min_),
85 "Histogram start value must be defined");
86 GMX_RELEASE_ASSERT(!isDefined(settings.max_) || settings.max_ > settings.min_,
87 "Histogram end value must be larger than start value");
88 GMX_RELEASE_ASSERT(!isDefined(settings.binWidth_) || settings.binWidth_ > 0.0,
89 "Histogram bin width must be positive");
90 GMX_RELEASE_ASSERT(settings.binCount_ >= 0,
91 "Histogram bin count must be positive");
93 if (!isDefined(settings.max_))
95 GMX_RELEASE_ASSERT(isDefined(settings.binWidth_) && settings.binCount_ > 0,
96 "Not all required values provided");
97 GMX_RELEASE_ASSERT(!settings.bRoundRange_,
98 "Rounding only supported for min/max ranges");
100 firstEdge_ = settings.min_;
101 binCount_ = settings.binCount_;
102 binWidth_ = settings.binWidth_;
103 if (settings.bIntegerBins_)
105 firstEdge_ -= 0.5 * binWidth_;
107 lastEdge_ = firstEdge_ + binCount_ * binWidth_;
111 GMX_RELEASE_ASSERT(!(isDefined(settings.binWidth_) && settings.binCount_ > 0),
112 "Conflicting histogram bin specifications");
113 GMX_RELEASE_ASSERT(isDefined(settings.binWidth_) || settings.binCount_ > 0,
114 "Not all required values provided");
116 if (settings.bRoundRange_)
118 GMX_RELEASE_ASSERT(!settings.bIntegerBins_,
119 "Rounding and integer bins cannot be combined");
120 GMX_RELEASE_ASSERT(isDefined(settings.binWidth_),
121 "Rounding only makes sense with defined binwidth");
122 binWidth_ = settings.binWidth_;
123 firstEdge_ = binWidth_ * floor(settings.min_ / binWidth_);
124 lastEdge_ = binWidth_ * ceil(settings.max_ / binWidth_);
125 binCount_ = static_cast<int>((lastEdge_ - firstEdge_) / binWidth_ + 0.5);
129 firstEdge_ = settings.min_;
130 lastEdge_ = settings.max_;
131 if (settings.binCount_ > 0)
133 binCount_ = settings.binCount_;
134 if (settings.bIntegerBins_)
136 GMX_RELEASE_ASSERT(settings.binCount_ > 1,
137 "Bin count must be at least two with integer bins");
138 binWidth_ = (lastEdge_ - firstEdge_) / (binCount_ - 1);
139 firstEdge_ -= 0.5 * binWidth_;
140 lastEdge_ += 0.5 * binWidth_;
144 binWidth_ = (lastEdge_ - firstEdge_) / binCount_;
149 binWidth_ = settings.binWidth_;
150 binCount_ = static_cast<int>((lastEdge_ - firstEdge_) / binWidth_ + 0.5);
151 if (settings.bIntegerBins_)
153 firstEdge_ -= 0.5 * binWidth_;
156 lastEdge_ = firstEdge_ + binCount_ * binWidth_;
161 inverseBinWidth_ = 1.0 / binWidth_;
162 bAll_ = settings.bIncludeAll_;
167 AnalysisHistogramSettings::findBin(real y) const
171 return bAll_ ? 0 : -1;
173 int bin = static_cast<int>((y - firstEdge_) * inverseBinWidth_);
174 if (bin >= binCount_)
176 return bAll_ ? binCount_ - 1 : -1;
182 /********************************************************************
183 * StaticAverageHistogram
190 * Represents copies of average histograms.
192 * Methods in AbstractAverageHistogram that return new histogram instances
193 * return objects of this class.
194 * Initialization of values is handled in those methods.
196 * \ingroup module_analysisdata
198 class StaticAverageHistogram : public AbstractAverageHistogram
201 StaticAverageHistogram();
202 //! Creates an average histogram module with defined bin parameters.
203 explicit StaticAverageHistogram(const AnalysisHistogramSettings &settings);
205 // Copy and assign disallowed by base.
208 StaticAverageHistogram::StaticAverageHistogram()
213 StaticAverageHistogram::StaticAverageHistogram(
214 const AnalysisHistogramSettings &settings)
215 : AbstractAverageHistogram(settings)
222 /********************************************************************
223 * AbstractAverageHistogram
226 AbstractAverageHistogram::AbstractAverageHistogram()
231 AbstractAverageHistogram::AbstractAverageHistogram(
232 const AnalysisHistogramSettings &settings)
233 : settings_(settings)
235 setRowCount(settings.binCount());
236 setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
237 settings.binWidth());
241 AbstractAverageHistogram::~AbstractAverageHistogram()
247 AbstractAverageHistogram::init(const AnalysisHistogramSettings &settings)
249 settings_ = settings;
250 setRowCount(settings.binCount());
251 setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
252 settings.binWidth());
256 AverageHistogramPointer
257 AbstractAverageHistogram::resampleDoubleBinWidth(bool bIntegerBins) const
262 nbins = (rowCount() + 1) / 2;
266 nbins = rowCount() / 2;
269 AverageHistogramPointer dest(
270 new StaticAverageHistogram(
271 histogramFromBins(xstart(), nbins, 2*xstep())
272 .integerBins(bIntegerBins)));
273 dest->setColumnCount(columnCount());
274 dest->allocateValues();
277 for (i = j = 0; i < nbins; ++i)
279 const bool bFirstHalfBin = (bIntegerBins && i == 0);
280 for (int c = 0; c < columnCount(); ++c)
291 v2 = value(j + 1, c);
295 dest->setValue(i, c, sqrt(v1 * v1 + v2 * v2));
299 dest->setValue(i, c, v1 + v2);
315 AverageHistogramPointer
316 AbstractAverageHistogram::clone() const
318 AverageHistogramPointer dest(new StaticAverageHistogram());
319 copyContents(this, dest.get());
325 AbstractAverageHistogram::normalizeProbability()
328 for (int i = 0; i < rowCount(); ++i)
332 scale(1.0 / (sum * xstep()));
337 AbstractAverageHistogram::scale(real norm)
339 for (int i = 0; i < rowCount(); ++i)
348 AbstractAverageHistogram::scaleVector(real norm[])
350 for (int i = 0; i < rowCount(); ++i)
352 value(i, 0) *= norm[i];
353 value(i, 1) *= norm[i];
358 /********************************************************************
359 * BasicAverageHistogramModule
366 * Implements average histogram module that averages per-frame histograms.
368 * This class is used for accumulating average histograms in per-frame
369 * histogram modules (those that use BasicHistogramImpl as their implementation
371 * There are two columns, first for the average and second for standard
374 * \ingroup module_analysisdata
376 class BasicAverageHistogramModule : public AbstractAverageHistogram,
377 public AnalysisDataModuleInterface
380 BasicAverageHistogramModule();
381 //! Creates an average histogram module with defined bin parameters.
382 explicit BasicAverageHistogramModule(const AnalysisHistogramSettings &settings);
384 using AbstractAverageHistogram::init;
386 virtual int flags() const;
388 virtual void dataStarted(AbstractAnalysisData *data);
389 virtual void frameStarted(const AnalysisDataFrameHeader &header);
390 virtual void pointsAdded(const AnalysisDataPointSetRef &points);
391 virtual void frameFinished(const AnalysisDataFrameHeader &header);
392 virtual void dataFinished();
395 //! Number of frames accumulated so far.
398 // Copy and assign disallowed by base.
401 BasicAverageHistogramModule::BasicAverageHistogramModule()
408 BasicAverageHistogramModule::BasicAverageHistogramModule(
409 const AnalysisHistogramSettings &settings)
410 : AbstractAverageHistogram(settings), frameCount_(0)
417 BasicAverageHistogramModule::flags() const
419 return efAllowMulticolumn;
424 BasicAverageHistogramModule::dataStarted(AbstractAnalysisData *data)
426 GMX_RELEASE_ASSERT(rowCount() == data->columnCount(),
427 "Inconsistent data sizes, something is wrong in the initialization");
433 BasicAverageHistogramModule::frameStarted(const AnalysisDataFrameHeader & /*header*/)
439 BasicAverageHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
441 int firstcol = points.firstColumn();
442 for (int i = 0; i < points.columnCount(); ++i)
444 real y = points.y(i);
445 value(firstcol + i, 0) += y;
446 value(firstcol + i, 1) += y * y;
452 BasicAverageHistogramModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
459 BasicAverageHistogramModule::dataFinished()
461 for (int i = 0; i < rowCount(); ++i)
463 real ave = value(i, 0) / frameCount_;
464 real std = sqrt(value(i, 1) / frameCount_ - ave * ave);
471 /********************************************************************
476 * Private implementation class for AnalysisDataSimpleHistogramModule and
477 * AnalysisDataWeightedHistogramModule.
479 * \ingroup module_analysisdata
481 class BasicHistogramImpl
484 //! Smart pointer to manage an BasicAverageHistogramModule object.
485 typedef boost::shared_ptr<BasicAverageHistogramModule>
486 BasicAverageHistogramModulePointer;
488 BasicHistogramImpl();
489 //! Creates an histogram impl with defined bin parameters.
490 explicit BasicHistogramImpl(const AnalysisHistogramSettings &settings);
491 ~BasicHistogramImpl();
494 * (Re)initializes the histogram from settings.
496 void init(const AnalysisHistogramSettings &settings);
498 * Initializes data storage frame when a new frame starts.
500 void initFrame(AnalysisDataStorageFrame *frame);
502 //! Storage implementation object.
503 AnalysisDataStorage storage_;
504 //! Settings for the histogram object.
505 AnalysisHistogramSettings settings_;
507 BasicAverageHistogramModulePointer averager_;
510 BasicHistogramImpl::BasicHistogramImpl()
511 : averager_(new BasicAverageHistogramModule())
516 BasicHistogramImpl::BasicHistogramImpl(const AnalysisHistogramSettings &settings)
517 : settings_(settings), averager_(new BasicAverageHistogramModule(settings))
522 BasicHistogramImpl::~BasicHistogramImpl()
527 void BasicHistogramImpl::init(const AnalysisHistogramSettings &settings)
529 settings_ = settings;
530 averager_->init(settings);
535 BasicHistogramImpl::initFrame(AnalysisDataStorageFrame *frame)
537 for (int i = 0; i < frame->columnCount(); ++i)
539 frame->setValue(i, 0.0);
543 } // namespace internal
546 /********************************************************************
547 * AnalysisDataSimpleHistogramModule
550 AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule()
551 : impl_(new internal::BasicHistogramImpl())
556 AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule(
557 const AnalysisHistogramSettings &settings)
558 : impl_(new internal::BasicHistogramImpl(settings))
563 AnalysisDataSimpleHistogramModule::~AnalysisDataSimpleHistogramModule()
568 void AnalysisDataSimpleHistogramModule::init(const AnalysisHistogramSettings &settings)
570 impl_->init(settings);
574 AbstractAverageHistogram &
575 AnalysisDataSimpleHistogramModule::averager()
577 return *impl_->averager_;
581 const AnalysisHistogramSettings &
582 AnalysisDataSimpleHistogramModule::settings() const
584 return impl_->settings_;
589 AnalysisDataSimpleHistogramModule::flags() const
591 return efAllowMulticolumn | efAllowMultipoint;
596 AnalysisDataSimpleHistogramModule::dataStarted(AbstractAnalysisData *data)
598 addModule(impl_->averager_);
599 setColumnCount(settings().binCount());
601 impl_->storage_.startDataStorage(this);
606 AnalysisDataSimpleHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
608 AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
609 impl_->initFrame(&frame);
614 AnalysisDataSimpleHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
616 AnalysisDataStorageFrame &frame =
617 impl_->storage_.currentFrame(points.frameIndex());
618 for (int i = 0; i < points.columnCount(); ++i)
620 int bin = settings().findBin(points.y(i));
623 frame.value(bin) += 1;
630 AnalysisDataSimpleHistogramModule::frameFinished(const AnalysisDataFrameHeader &header)
632 impl_->storage_.finishFrame(header.index());
637 AnalysisDataSimpleHistogramModule::dataFinished()
644 AnalysisDataSimpleHistogramModule::tryGetDataFrameInternal(int index) const
646 return impl_->storage_.tryGetDataFrame(index);
651 AnalysisDataSimpleHistogramModule::requestStorageInternal(int nframes)
653 return impl_->storage_.requestStorage(nframes);
657 /********************************************************************
658 * AnalysisDataWeightedHistogramModule
661 AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule()
662 : impl_(new internal::BasicHistogramImpl())
667 AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule(
668 const AnalysisHistogramSettings &settings)
669 : impl_(new internal::BasicHistogramImpl(settings))
674 AnalysisDataWeightedHistogramModule::~AnalysisDataWeightedHistogramModule()
679 void AnalysisDataWeightedHistogramModule::init(const AnalysisHistogramSettings &settings)
681 impl_->init(settings);
685 AbstractAverageHistogram &
686 AnalysisDataWeightedHistogramModule::averager()
688 return *impl_->averager_;
692 const AnalysisHistogramSettings &
693 AnalysisDataWeightedHistogramModule::settings() const
695 return impl_->settings_;
700 AnalysisDataWeightedHistogramModule::flags() const
702 return efAllowMulticolumn | efAllowMultipoint;
707 AnalysisDataWeightedHistogramModule::dataStarted(AbstractAnalysisData *data)
709 addModule(impl_->averager_);
710 setColumnCount(settings().binCount());
712 impl_->storage_.startDataStorage(this);
717 AnalysisDataWeightedHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
719 AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
720 impl_->initFrame(&frame);
725 AnalysisDataWeightedHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
727 if (points.firstColumn() != 0 || points.columnCount() < 2)
729 GMX_THROW(APIError("Invalid data layout"));
731 int bin = settings().findBin(points.y(0));
734 AnalysisDataStorageFrame &frame =
735 impl_->storage_.currentFrame(points.frameIndex());
736 for (int i = 1; i < points.columnCount(); ++i)
738 frame.value(bin) += points.y(i);
745 AnalysisDataWeightedHistogramModule::frameFinished(const AnalysisDataFrameHeader &header)
747 impl_->storage_.finishFrame(header.index());
752 AnalysisDataWeightedHistogramModule::dataFinished()
759 AnalysisDataWeightedHistogramModule::tryGetDataFrameInternal(int index) const
761 return impl_->storage_.tryGetDataFrame(index);
766 AnalysisDataWeightedHistogramModule::requestStorageInternal(int nframes)
768 return impl_->storage_.requestStorage(nframes);
772 /********************************************************************
773 * AnalysisDataBinAverageModule
776 AnalysisDataBinAverageModule::AnalysisDataBinAverageModule()
782 AnalysisDataBinAverageModule::AnalysisDataBinAverageModule(
783 const AnalysisHistogramSettings &settings)
784 : settings_(settings)
787 setRowCount(settings.binCount());
788 setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
789 settings.binWidth());
793 AnalysisDataBinAverageModule::~AnalysisDataBinAverageModule()
799 AnalysisDataBinAverageModule::init(const AnalysisHistogramSettings &settings)
801 settings_ = settings;
802 setRowCount(settings.binCount());
803 setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
804 settings.binWidth());
809 AnalysisDataBinAverageModule::flags() const
811 return efAllowMulticolumn | efAllowMultipoint;
816 AnalysisDataBinAverageModule::dataStarted(AbstractAnalysisData * /*data*/)
823 AnalysisDataBinAverageModule::frameStarted(const AnalysisDataFrameHeader & /*header*/)
829 AnalysisDataBinAverageModule::pointsAdded(const AnalysisDataPointSetRef &points)
831 if (points.firstColumn() != 0 || points.columnCount() < 2)
833 GMX_THROW(APIError("Invalid data layout"));
835 int bin = settings().findBin(points.y(0));
838 for (int i = 1; i < points.columnCount(); ++i)
840 real y = points.y(i);
842 value(bin, 1) += y * y;
844 value(bin, 2) += points.columnCount() - 1;
850 AnalysisDataBinAverageModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
856 AnalysisDataBinAverageModule::dataFinished()
858 for (int i = 0; i < settings().binCount(); ++i)
860 real n = value(i, 2);
863 real ave = value(i, 0) / n;
864 real std = sqrt(value(i, 1) / n - ave * ave);