/*
+ * This file is part of the GROMACS molecular simulation package.
*
- * This source code is part of
+ * Copyright (c) 2010,2011,2012,2013,2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
*
- * G R O M A C S
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
*
- * GROningen MAchine for Chemical Simulations
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
*
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2009, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*
- * If you want to redistribute modifications, please consider that
- * scientific software is very special. Version control is crucial -
- * bugs must be traceable. We will be happy to consider code for
- * inclusion in the official distribution, but derived work must not
- * be called official GROMACS. Details are found in the README & COPYING
- * files - if they are missing, get the official version at www.gromacs.org.
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
*
* To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
- *
- * For more info, check our website at http://www.gromacs.org
+ * the research papers on the package. Check out http://www.gromacs.org.
*/
/*! \internal \file
* \brief
* Implements classes in histogram.h.
*
- * \author Teemu Murtola <teemu.murtola@cbr.su.se>
+ * \author Teemu Murtola <teemu.murtola@gmail.com>
* \ingroup module_analysisdata
*/
-#include "gromacs/analysisdata/modules/histogram.h"
+#include "gmxpre.h"
+
+#include "histogram.h"
#include <cmath>
#include <limits>
+#include <vector>
#include "gromacs/analysisdata/dataframe.h"
#include "gromacs/analysisdata/datastorage.h"
#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/gmxassert.h"
-#include "histogram-impl.h"
+#include "frameaverager.h"
-static const real UNDEFINED = std::numeric_limits<real>::max();
-static bool isDefined(real value)
+namespace
+{
+
+//! Value used to signify that a real-valued histogram setting is not set.
+const real UNDEFINED = std::numeric_limits<real>::max();
+//! Checks whether \p value is defined.
+bool isDefined(real value)
{
return value != UNDEFINED;
}
+} // namespace
+
namespace gmx
{
else
{
firstEdge_ = settings.min_;
- lastEdge_ = settings.max_;
+ lastEdge_ = settings.max_;
if (settings.binCount_ > 0)
{
binCount_ = settings.binCount_;
}
+/********************************************************************
+ * StaticAverageHistogram
+ */
+
+namespace
+{
+
+/*! \brief
+ * Represents copies of average histograms.
+ *
+ * Methods in AbstractAverageHistogram that return new histogram instances
+ * return objects of this class.
+ * Initialization of values is handled in those methods.
+ *
+ * \ingroup module_analysisdata
+ */
+class StaticAverageHistogram : public AbstractAverageHistogram
+{
+ public:
+ StaticAverageHistogram();
+ //! Creates an average histogram module with defined bin parameters.
+ explicit StaticAverageHistogram(const AnalysisHistogramSettings &settings);
+
+ // Copy and assign disallowed by base.
+};
+
+StaticAverageHistogram::StaticAverageHistogram()
+{
+}
+
+
+StaticAverageHistogram::StaticAverageHistogram(
+ const AnalysisHistogramSettings &settings)
+ : AbstractAverageHistogram(settings)
+{
+}
+
+} // namespace
+
+
/********************************************************************
* AbstractAverageHistogram
*/
}
AverageHistogramPointer dest(
- new internal::StaticAverageHistogram(
- histogramFromBins(xstart(), nbins, 2*xstep())
- .integerBins(bIntegerBins)));
+ new StaticAverageHistogram(
+ histogramFromBins(xstart(), nbins, 2*xstep())
+ .integerBins(bIntegerBins)));
dest->setColumnCount(columnCount());
dest->allocateValues();
for (int c = 0; c < columnCount(); ++c)
{
real v1, v2;
+ real e1, e2;
if (bFirstHalfBin)
{
- v1 = value(0, c);
+ v1 = value(0, c).value();
+ e1 = value(0, c).error();
v2 = 0;
+ e2 = 0;
}
else
{
- v1 = value(j, c);
- v2 = value(j + 1, c);
- }
- if (c == 1)
- {
- dest->setValue(i, c, sqrt(v1 * v1 + v2 * v2));
- }
- else
- {
- dest->setValue(i, c, v1 + v2);
+ v1 = value(j, c).value();
+ e1 = value(j, c).error();
+ v2 = value(j + 1, c).value();
+ e2 = value(j + 1, c).error();
}
+ dest->value(i, c).setValue(v1 + v2, std::sqrt(e1 * e1 + e2 * e2));
}
if (bFirstHalfBin)
{
AverageHistogramPointer
AbstractAverageHistogram::clone() const
{
- AverageHistogramPointer dest(new internal::StaticAverageHistogram());
+ AverageHistogramPointer dest(new StaticAverageHistogram());
copyContents(this, dest.get());
return dest;
}
void
AbstractAverageHistogram::normalizeProbability()
{
- real sum = 0;
- for (int i = 0; i < rowCount(); ++i)
+ for (int c = 0; c < columnCount(); ++c)
{
- sum += value(i, 0);
+ real sum = 0;
+ for (int i = 0; i < rowCount(); ++i)
+ {
+ sum += value(i, c).value();
+ }
+ if (sum > 0.0)
+ {
+ scaleSingle(c, 1.0 / (sum * xstep()));
+ }
}
- scale(1.0 / (sum * xstep()));
}
void
-AbstractAverageHistogram::scale(real norm)
+AbstractAverageHistogram::scaleSingle(int index, real factor)
{
for (int i = 0; i < rowCount(); ++i)
{
- value(i, 0) *= norm;
- value(i, 1) *= norm;
+ value(i, index).value() *= factor;
+ value(i, index).error() *= factor;
}
}
void
-AbstractAverageHistogram::scaleVector(real norm[])
+AbstractAverageHistogram::scaleAll(real factor)
{
- for (int i = 0; i < rowCount(); ++i)
+ for (int i = 0; i < columnCount(); ++i)
{
- value(i, 0) *= norm[i];
- value(i, 1) *= norm[i];
+ scaleSingle(i, factor);
+ }
+}
+
+
+void
+AbstractAverageHistogram::scaleAllByVector(real factor[])
+{
+ for (int c = 0; c < columnCount(); ++c)
+ {
+ for (int i = 0; i < rowCount(); ++i)
+ {
+ value(i, c).value() *= factor[i];
+ value(i, c).error() *= factor[i];
+ }
}
}
/********************************************************************
- * StaticAverageHistogram
+ * BasicAverageHistogramModule
*/
namespace internal
{
-StaticAverageHistogram::StaticAverageHistogram()
+/*! \internal
+ * \brief
+ * Implements average histogram module that averages per-frame histograms.
+ *
+ * This class is used for accumulating average histograms in per-frame
+ * histogram modules (those that use BasicHistogramImpl as their implementation
+ * class).
+ * There are two columns, first for the average and second for standard
+ * deviation.
+ *
+ * \ingroup module_analysisdata
+ */
+class BasicAverageHistogramModule : public AbstractAverageHistogram,
+ public AnalysisDataModuleSerial
{
-}
+ public:
+ BasicAverageHistogramModule();
+ //! Creates an average histogram module with defined bin parameters.
+ explicit BasicAverageHistogramModule(const AnalysisHistogramSettings &settings);
+ using AbstractAverageHistogram::init;
-StaticAverageHistogram::StaticAverageHistogram(
- const AnalysisHistogramSettings &settings)
- : AbstractAverageHistogram(settings)
-{
-}
+ virtual int flags() const;
+ virtual void dataStarted(AbstractAnalysisData *data);
+ virtual void frameStarted(const AnalysisDataFrameHeader &header);
+ virtual void pointsAdded(const AnalysisDataPointSetRef &points);
+ virtual void frameFinished(const AnalysisDataFrameHeader &header);
+ virtual void dataFinished();
-/********************************************************************
- * BasicAverageHistogramModule
- */
+ private:
+ //! Averaging helper objects for each input data set.
+ std::vector<AnalysisDataFrameAverager> averagers_;
+
+ // Copy and assign disallowed by base.
+};
BasicAverageHistogramModule::BasicAverageHistogramModule()
- : frameCount_(0)
{
- setColumnCount(2);
}
BasicAverageHistogramModule::BasicAverageHistogramModule(
const AnalysisHistogramSettings &settings)
- : AbstractAverageHistogram(settings), frameCount_(0)
+ : AbstractAverageHistogram(settings)
{
- setColumnCount(2);
}
int
BasicAverageHistogramModule::flags() const
{
- return efAllowMulticolumn;
+ return efAllowMulticolumn | efAllowMultipleDataSets;
}
void
BasicAverageHistogramModule::dataStarted(AbstractAnalysisData *data)
{
- GMX_RELEASE_ASSERT(rowCount() == data->columnCount(),
- "Inconsistent data sizes, something is wrong in the initialization");
- allocateValues();
+ setColumnCount(data->dataSetCount());
+ averagers_.resize(data->dataSetCount());
+ for (int i = 0; i < data->dataSetCount(); ++i)
+ {
+ GMX_RELEASE_ASSERT(rowCount() == data->columnCount(i),
+ "Inconsistent data sizes, something is wrong in the initialization");
+ averagers_[i].setColumnCount(data->columnCount(i));
+ }
}
void
BasicAverageHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
{
- int firstcol = points.firstColumn();
- for (int i = 0; i < points.columnCount(); ++i)
- {
- real y = points.y(i);
- value(firstcol + i, 0) += y;
- value(firstcol + i, 1) += y * y;
- }
+ averagers_[points.dataSetIndex()].addPoints(points);
}
void
BasicAverageHistogramModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
{
- ++frameCount_;
}
void
BasicAverageHistogramModule::dataFinished()
{
- for (int i = 0; i < rowCount(); ++i)
+ allocateValues();
+ for (int i = 0; i < columnCount(); ++i)
{
- real ave = value(i, 0) / frameCount_;
- real std = sqrt(value(i, 1) / frameCount_ - ave * ave);
- setValue(i, 0, ave);
- setValue(i, 1, std);
+ averagers_[i].finish();
+ for (int j = 0; j < rowCount(); ++j)
+ {
+ value(j, i).setValue(averagers_[i].average(j),
+ std::sqrt(averagers_[i].variance(j)));
+ }
}
}
* BasicHistogramImpl
*/
+/*! \internal \brief
+ * Private implementation class for AnalysisDataSimpleHistogramModule and
+ * AnalysisDataWeightedHistogramModule.
+ *
+ * \ingroup module_analysisdata
+ */
+class BasicHistogramImpl
+{
+ public:
+ //! Smart pointer to manage an BasicAverageHistogramModule object.
+ typedef boost::shared_ptr<BasicAverageHistogramModule>
+ BasicAverageHistogramModulePointer;
+
+ BasicHistogramImpl();
+ //! Creates an histogram impl with defined bin parameters.
+ explicit BasicHistogramImpl(const AnalysisHistogramSettings &settings);
+ ~BasicHistogramImpl();
+
+ /*! \brief
+ * (Re)initializes the histogram from settings.
+ */
+ void init(const AnalysisHistogramSettings &settings);
+ /*! \brief
+ * Initializes data storage frame when a new frame starts.
+ */
+ void initFrame(int dataSetCount, AnalysisDataStorageFrame *frame);
+
+ //! Storage implementation object.
+ AnalysisDataStorage storage_;
+ //! Settings for the histogram object.
+ AnalysisHistogramSettings settings_;
+ //! Averager module.
+ BasicAverageHistogramModulePointer averager_;
+};
+
BasicHistogramImpl::BasicHistogramImpl()
: averager_(new BasicAverageHistogramModule())
{
void
-BasicHistogramImpl::initFrame(AnalysisDataStorageFrame *frame)
+BasicHistogramImpl::initFrame(int dataSetCount, AnalysisDataStorageFrame *frame)
{
- for (int i = 0; i < frame->columnCount(); ++i)
+ for (int s = 0; s < dataSetCount; ++s)
{
- frame->setValue(i, 0.0);
+ frame->selectDataSet(s);
+ for (int i = 0; i < frame->columnCount(); ++i)
+ {
+ frame->setValue(i, 0.0);
+ }
}
+ frame->selectDataSet(0);
}
-} // namespace internal
+} // namespace internal
/********************************************************************
}
+int
+AnalysisDataSimpleHistogramModule::frameCount() const
+{
+ return impl_->storage_.frameCount();
+}
+
+
int
AnalysisDataSimpleHistogramModule::flags() const
{
- return efAllowMulticolumn | efAllowMultipoint;
+ return efAllowMulticolumn | efAllowMultipoint | efAllowMissing
+ | efAllowMultipleDataSets;
}
-void
-AnalysisDataSimpleHistogramModule::dataStarted(AbstractAnalysisData *data)
+bool
+AnalysisDataSimpleHistogramModule::parallelDataStarted(
+ AbstractAnalysisData *data,
+ const AnalysisDataParallelOptions &options)
{
addModule(impl_->averager_);
- setColumnCount(settings().binCount());
- notifyDataStart();
- impl_->storage_.startDataStorage(this);
+ setDataSetCount(data->dataSetCount());
+ for (int i = 0; i < data->dataSetCount(); ++i)
+ {
+ setColumnCount(i, settings().binCount());
+ }
+ impl_->storage_.startParallelDataStorage(this, &moduleManager(), options);
+ return true;
}
AnalysisDataSimpleHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
{
AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
- impl_->initFrame(&frame);
+ impl_->initFrame(dataSetCount(), &frame);
}
{
AnalysisDataStorageFrame &frame =
impl_->storage_.currentFrame(points.frameIndex());
+ frame.selectDataSet(points.dataSetIndex());
for (int i = 0; i < points.columnCount(); ++i)
{
- int bin = settings().findBin(points.y(i));
- if (bin != -1)
+ if (points.present(i))
{
- frame.value(bin) += 1;
+ const int bin = settings().findBin(points.y(i));
+ if (bin != -1)
+ {
+ frame.value(bin) += 1;
+ }
}
}
}
void
AnalysisDataSimpleHistogramModule::dataFinished()
{
- notifyDataFinish();
+ impl_->storage_.finishDataStorage();
}
}
+int
+AnalysisDataWeightedHistogramModule::frameCount() const
+{
+ return impl_->storage_.frameCount();
+}
+
+
int
AnalysisDataWeightedHistogramModule::flags() const
{
- return efAllowMulticolumn | efAllowMultipoint;
+ return efAllowMulticolumn | efAllowMultipoint | efAllowMultipleDataSets;
}
-void
-AnalysisDataWeightedHistogramModule::dataStarted(AbstractAnalysisData *data)
+bool
+AnalysisDataWeightedHistogramModule::parallelDataStarted(
+ AbstractAnalysisData *data,
+ const AnalysisDataParallelOptions &options)
{
addModule(impl_->averager_);
- setColumnCount(settings().binCount());
- notifyDataStart();
- impl_->storage_.startDataStorage(this);
+ setDataSetCount(data->dataSetCount());
+ for (int i = 0; i < data->dataSetCount(); ++i)
+ {
+ setColumnCount(i, settings().binCount());
+ }
+ impl_->storage_.startParallelDataStorage(this, &moduleManager(), options);
+ return true;
}
AnalysisDataWeightedHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
{
AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
- impl_->initFrame(&frame);
+ impl_->initFrame(dataSetCount(), &frame);
}
{
AnalysisDataStorageFrame &frame =
impl_->storage_.currentFrame(points.frameIndex());
+ frame.selectDataSet(points.dataSetIndex());
for (int i = 1; i < points.columnCount(); ++i)
{
frame.value(bin) += points.y(i);
void
AnalysisDataWeightedHistogramModule::dataFinished()
{
- notifyDataFinish();
+ impl_->storage_.finishDataStorage();
}
* AnalysisDataBinAverageModule
*/
+class AnalysisDataBinAverageModule::Impl
+{
+ public:
+ Impl() {}
+ explicit Impl(const AnalysisHistogramSettings &settings)
+ : settings_(settings)
+ {
+ }
+
+ //! Histogram settings.
+ AnalysisHistogramSettings settings_;
+ //! Averaging helper objects for each input data set.
+ std::vector<AnalysisDataFrameAverager> averagers_;
+};
+
AnalysisDataBinAverageModule::AnalysisDataBinAverageModule()
+ : impl_(new Impl())
{
setColumnCount(3);
}
AnalysisDataBinAverageModule::AnalysisDataBinAverageModule(
const AnalysisHistogramSettings &settings)
- : settings_(settings)
+ : impl_(new Impl(settings))
{
- setColumnCount(3);
setRowCount(settings.binCount());
setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
settings.binWidth());
void
AnalysisDataBinAverageModule::init(const AnalysisHistogramSettings &settings)
{
- settings_ = settings;
+ impl_->settings_ = settings;
setRowCount(settings.binCount());
setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
settings.binWidth());
}
+const AnalysisHistogramSettings &
+AnalysisDataBinAverageModule::settings() const
+{
+ return impl_->settings_;
+}
+
+
int
AnalysisDataBinAverageModule::flags() const
{
- return efAllowMulticolumn | efAllowMultipoint;
+ return efAllowMulticolumn | efAllowMultipoint | efAllowMultipleDataSets;
}
void
-AnalysisDataBinAverageModule::dataStarted(AbstractAnalysisData * /*data*/)
+AnalysisDataBinAverageModule::dataStarted(AbstractAnalysisData *data)
{
- allocateValues();
+ setColumnCount(data->dataSetCount());
+ impl_->averagers_.resize(data->dataSetCount());
+ for (int i = 0; i < data->dataSetCount(); ++i)
+ {
+ impl_->averagers_[i].setColumnCount(rowCount());
+ }
}
int bin = settings().findBin(points.y(0));
if (bin != -1)
{
+ AnalysisDataFrameAverager &averager = impl_->averagers_[points.dataSetIndex()];
for (int i = 1; i < points.columnCount(); ++i)
{
- real y = points.y(i);
- value(bin, 0) += y;
- value(bin, 1) += y * y;
+ averager.addValue(bin, points.y(i));
}
- value(bin, 2) += points.columnCount() - 1;
}
}
void
AnalysisDataBinAverageModule::dataFinished()
{
- for (int i = 0; i < settings().binCount(); ++i)
+ allocateValues();
+ for (int i = 0; i < columnCount(); ++i)
{
- real n = value(i, 2);
- if (n > 0)
+ AnalysisDataFrameAverager &averager = impl_->averagers_[i];
+ averager.finish();
+ for (int j = 0; j < rowCount(); ++j)
{
- real ave = value(i, 0) / n;
- real std = sqrt(value(i, 1) / n - ave * ave);
- setValue(i, 0, ave);
- setValue(i, 1, std);
+ value(j, i).setValue(averager.average(j),
+ std::sqrt(averager.variance(j)));
}
}
valuesReady();
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff