/*
+ * 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
* \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"
+#include "frameaverager.h"
+
namespace
{
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)
+ {
+ 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()));
+ }
+ }
+}
+
+void
+AbstractAverageHistogram::makeCumulative()
+{
+ 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();
+ // Clear the error, as we don't cumulate that.
+ value(i, c).clear();
+ value(i, c).setValue(sum);
+ }
}
- scale(1.0 / (sum * xstep()));
+ 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)
{
- 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];
+ }
}
}
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:
- //! Number of frames accumulated so far.
- int frameCount_;
+ //! 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.
+/*! \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();
}
* 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();