/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2010,2011,2012,2013, by the GROMACS development team, led by
- * David van der Spoel, Berk Hess, Erik Lindahl, and including many
- * others, as listed in the AUTHORS file in the top-level source
- * directory and at http://www.gromacs.org.
+ * 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.
*
* GROMACS is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public License
* \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/analysisdata/framelocaldata.h"
+#include "gromacs/utility/basedefinitions.h"
#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/gmxassert.h"
namespace
{
-/*! \internal \brief
+/*! \brief
* Represents copies of average histograms.
*
* Methods in AbstractAverageHistogram that return new histogram instances
AverageHistogramPointer dest(
new StaticAverageHistogram(
- histogramFromBins(xstart(), nbins, 2*xstep())
+ histogramFromBins(settings().firstEdge(), 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 dest(new StaticAverageHistogram());
copyContents(this, dest.get());
+ dest->settings_ = settings_;
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);
+ double 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::makeCumulative()
+{
+ for (int c = 0; c < columnCount(); ++c)
+ {
+ double sum = 0;
+ for (int i = 0; i < rowCount(); ++i)
+ {
+ sum += value(i, c).value();
+ // Clear the error, as we don't cumulate that.
+ value(i, c).clear();
+ value(i, c).setValue(sum);
+ }
+ }
+ setXAxis(settings().firstEdge() + settings().binWidth(),
+ settings().binWidth());
}
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)
+ {
+ scaleSingle(i, factor);
+ }
+}
+
+
+void
+AbstractAverageHistogram::scaleAllByVector(real factor[])
+{
+ for (int c = 0; c < columnCount(); ++c)
{
- value(i, 0) *= norm[i];
- value(i, 1) *= norm[i];
+ for (int i = 0; i < rowCount(); ++i)
+ {
+ value(i, c).value() *= factor[i];
+ value(i, c).error() *= factor[i];
+ }
}
}
namespace internal
{
-/*! \internal \brief
+/*! \internal
+ * \brief
* Implements average histogram module that averages per-frame histograms.
*
* This class is used for accumulating average histograms in per-frame
* \ingroup module_analysisdata
*/
class BasicAverageHistogramModule : public AbstractAverageHistogram,
- public AnalysisDataModuleInterface
+ public AnalysisDataModuleSerial
{
public:
BasicAverageHistogramModule();
virtual void dataFinished();
private:
- //! Averaging helper object.
- AnalysisDataFrameAverager averager_;
+ //! Averaging helper objects for each input data set.
+ std::vector<AnalysisDataFrameAverager> averagers_;
// Copy and assign disallowed by base.
};
BasicAverageHistogramModule::BasicAverageHistogramModule()
{
- setColumnCount(2);
}
const AnalysisHistogramSettings &settings)
: 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");
- averager_.setColumnCount(rowCount());
+ 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)
{
- averager_.addPoints(points);
+ averagers_[points.dataSetIndex()].addPoints(points);
}
void
BasicAverageHistogramModule::dataFinished()
{
- averager_.finish();
allocateValues();
- for (int i = 0; i < rowCount(); ++i)
+ for (int i = 0; i < columnCount(); ++i)
{
- setValue(i, 0, averager_.average(i));
- setValue(i, 1, sqrt(averager_.variance(i)));
+ 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.
+/*! \internal
+ * \brief
+ * Base class for private implementation classes for histogram modules.
+ *
+ * Actual implementation classes are derived from this and add an accumulation
+ * data member that is specific to the histogram type in question.
+ * This is done like this to keep implementation details out of the header, and
+ * to not unnecessarily duplicate code.
*
* \ingroup module_analysisdata
*/
BasicHistogramImpl();
//! Creates an histogram impl with defined bin parameters.
explicit BasicHistogramImpl(const AnalysisHistogramSettings &settings);
- ~BasicHistogramImpl();
+ // Virtual only for simplicity.
+ virtual ~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(AnalysisDataStorageFrame *frame);
//! Storage implementation object.
AnalysisDataStorage storage_;
averager_->init(settings);
}
-
-void
-BasicHistogramImpl::initFrame(AnalysisDataStorageFrame *frame)
-{
- for (int i = 0; i < frame->columnCount(); ++i)
- {
- frame->setValue(i, 0.0);
- }
-}
-
} // namespace internal
* AnalysisDataSimpleHistogramModule
*/
+/*! \internal \brief
+ * Private implementation class for AnalysisDataSimpleHistogramModule.
+ *
+ * \ingroup module_analysisdata
+ */
+class AnalysisDataSimpleHistogramModule::Impl : public internal::BasicHistogramImpl
+{
+ public:
+ //! Shorthand for the per-frame accumulation data structure type.
+ typedef AnalysisDataFrameLocalData<gmx_int64_t> FrameLocalData;
+
+ Impl() {}
+ //! Creates an histogram impl with defined bin parameters.
+ explicit Impl(const AnalysisHistogramSettings &settings)
+ : BasicHistogramImpl(settings)
+ {
+ }
+
+ //! Accumulates the histogram within a frame.
+ FrameLocalData accumulator_;
+};
+
AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule()
- : impl_(new internal::BasicHistogramImpl())
+ : impl_(new Impl())
{
}
AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule(
const AnalysisHistogramSettings &settings)
- : impl_(new internal::BasicHistogramImpl(settings))
+ : impl_(new Impl(settings))
{
}
}
+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);
+ const int dataSetCount = data->dataSetCount();
+ const int columnCount = settings().binCount();
+ setDataSetCount(dataSetCount);
+ impl_->accumulator_.setDataSetCount(dataSetCount);
+ for (int i = 0; i < dataSetCount; ++i)
+ {
+ setColumnCount(i, columnCount);
+ impl_->accumulator_.setColumnCount(i, columnCount);
+ }
+ impl_->accumulator_.init(options);
+ impl_->storage_.startParallelDataStorage(this, &moduleManager(), options);
+ return true;
}
void
AnalysisDataSimpleHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
{
- AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
- impl_->initFrame(&frame);
+ impl_->accumulator_.frameData(header.index()).clear();
}
void
AnalysisDataSimpleHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
{
- AnalysisDataStorageFrame &frame =
- impl_->storage_.currentFrame(points.frameIndex());
+ Impl::FrameLocalData::DataSetHandle handle
+ = impl_->accumulator_.frameDataSet(points.frameIndex(), 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)
+ {
+ handle.value(bin) += 1;
+ }
}
}
}
void
AnalysisDataSimpleHistogramModule::frameFinished(const AnalysisDataFrameHeader &header)
{
- impl_->storage_.finishFrame(header.index());
+ Impl::FrameLocalData::FrameHandle handle
+ = impl_->accumulator_.frameData(header.index());
+ AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
+ const int columnCount = settings().binCount();
+ for (int s = 0; s < dataSetCount(); ++s)
+ {
+ Impl::FrameLocalData::DataSetHandle dataSet = handle.dataSet(s);
+ frame.selectDataSet(s);
+ for (int i = 0; i < columnCount; ++i)
+ {
+ frame.setValue(i, dataSet.value(i));
+ }
+ }
+ frame.finishFrame();
}
void
AnalysisDataSimpleHistogramModule::dataFinished()
{
- notifyDataFinish();
+ impl_->storage_.finishDataStorage();
}
* AnalysisDataWeightedHistogramModule
*/
+/*! \internal \brief
+ * Private implementation class for AnalysisDataWeightedHistogramModule.
+ *
+ * \ingroup module_analysisdata
+ */
+class AnalysisDataWeightedHistogramModule::Impl : public internal::BasicHistogramImpl
+{
+ public:
+ //! Shorthand for the per-frame accumulation data structure type.
+ typedef AnalysisDataFrameLocalData<double> FrameLocalData;
+
+ Impl() {}
+ //! Creates an histogram impl with defined bin parameters.
+ explicit Impl(const AnalysisHistogramSettings &settings)
+ : BasicHistogramImpl(settings)
+ {
+ }
+
+ //! Accumulates the histogram within a frame.
+ FrameLocalData accumulator_;
+};
+
AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule()
- : impl_(new internal::BasicHistogramImpl())
+ : impl_(new Impl())
{
}
AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule(
const AnalysisHistogramSettings &settings)
- : impl_(new internal::BasicHistogramImpl(settings))
+ : impl_(new Impl(settings))
{
}
}
+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);
+ const int dataSetCount = data->dataSetCount();
+ const int columnCount = settings().binCount();
+ setDataSetCount(dataSetCount);
+ impl_->accumulator_.setDataSetCount(dataSetCount);
+ for (int i = 0; i < dataSetCount; ++i)
+ {
+ setColumnCount(i, columnCount);
+ impl_->accumulator_.setColumnCount(i, columnCount);
+ }
+ impl_->accumulator_.init(options);
+ impl_->storage_.startParallelDataStorage(this, &moduleManager(), options);
+ return true;
}
void
AnalysisDataWeightedHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
{
- AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
- impl_->initFrame(&frame);
+ impl_->accumulator_.frameData(header.index()).clear();
}
int bin = settings().findBin(points.y(0));
if (bin != -1)
{
- AnalysisDataStorageFrame &frame =
- impl_->storage_.currentFrame(points.frameIndex());
+ Impl::FrameLocalData::DataSetHandle handle
+ = impl_->accumulator_.frameDataSet(points.frameIndex(), points.dataSetIndex());
for (int i = 1; i < points.columnCount(); ++i)
{
- frame.value(bin) += points.y(i);
+ handle.value(bin) += points.y(i);
}
}
}
void
AnalysisDataWeightedHistogramModule::frameFinished(const AnalysisDataFrameHeader &header)
{
- impl_->storage_.finishFrame(header.index());
+ Impl::FrameLocalData::FrameHandle handle
+ = impl_->accumulator_.frameData(header.index());
+ AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
+ const int columnCount = settings().binCount();
+ for (int s = 0; s < dataSetCount(); ++s)
+ {
+ Impl::FrameLocalData::DataSetHandle dataSet = handle.dataSet(s);
+ frame.selectDataSet(s);
+ for (int i = 0; i < columnCount; ++i)
+ {
+ frame.setValue(i, dataSet.value(i));
+ }
+ }
+ frame.finishFrame();
}
void
AnalysisDataWeightedHistogramModule::dataFinished()
{
- notifyDataFinish();
+ impl_->storage_.finishDataStorage();
}
}
//! Histogram settings.
- AnalysisHistogramSettings settings_;
- //! Averaging helper object.
- AnalysisDataFrameAverager averager_;
+ AnalysisHistogramSettings settings_;
+ //! Averaging helper objects for each input data set.
+ std::vector<AnalysisDataFrameAverager> averagers_;
};
AnalysisDataBinAverageModule::AnalysisDataBinAverageModule()
const AnalysisHistogramSettings &settings)
: impl_(new Impl(settings))
{
- setColumnCount(3);
setRowCount(settings.binCount());
setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
settings.binWidth());
int
AnalysisDataBinAverageModule::flags() const
{
- return efAllowMulticolumn | efAllowMultipoint;
+ return efAllowMulticolumn | efAllowMultipoint | efAllowMultipleDataSets;
}
void
-AnalysisDataBinAverageModule::dataStarted(AbstractAnalysisData * /*data*/)
+AnalysisDataBinAverageModule::dataStarted(AbstractAnalysisData *data)
{
- impl_->averager_.setColumnCount(rowCount());
+ 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)
{
- impl_->averager_.addValue(bin, points.y(i));
+ averager.addValue(bin, points.y(i));
}
}
}
void
AnalysisDataBinAverageModule::dataFinished()
{
- impl_->averager_.finish();
allocateValues();
- for (int i = 0; i < settings().binCount(); ++i)
+ for (int i = 0; i < columnCount(); ++i)
{
- setValue(i, 0, impl_->averager_.average(i));
- setValue(i, 1, std::sqrt(impl_->averager_.variance(i)));
- setValue(i, 2, impl_->averager_.sampleCount(i));
+ AnalysisDataFrameAverager &averager = impl_->averagers_[i];
+ averager.finish();
+ for (int j = 0; j < rowCount(); ++j)
+ {
+ value(j, i).setValue(averager.average(j),
+ std::sqrt(averager.variance(j)));
+ }
}
valuesReady();
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff