More precise analysisdata histogram accumulation
[alexxy/gromacs.git] / src / gromacs / analysisdata / modules / histogram.cpp
index aafd293be5814d08dee47498828c2a76efce49ad..24f67024a74217c8391166d56b187ce72c51080b 100644 (file)
@@ -1,59 +1,75 @@
 /*
+ * 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/fatalerror/exceptions.h"
-#include "gromacs/fatalerror/gmxassert.h"
+#include "gromacs/analysisdata/framelocaldata.h"
+#include "gromacs/utility/basedefinitions.h"
+#include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/gmxassert.h"
+
+#include "frameaverager.h"
 
-#include "histogram-impl.h"
+namespace
+{
 
-static const real UNDEFINED = std::numeric_limits<real>::max();
-static bool isDefined(real value)
+//! 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
 {
 
@@ -129,7 +145,7 @@ AnalysisHistogramSettings::AnalysisHistogramSettings(
         else
         {
             firstEdge_     = settings.min_;
-            lastEdge_     = settings.max_;
+            lastEdge_      = settings.max_;
             if (settings.binCount_ > 0)
             {
                 binCount_ = settings.binCount_;
@@ -181,6 +197,46 @@ AnalysisHistogramSettings::findBin(real y) const
 }
 
 
+/********************************************************************
+ * 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
  */
@@ -229,9 +285,9 @@ AbstractAverageHistogram::resampleDoubleBinWidth(bool bIntegerBins) const
     }
 
     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();
 
@@ -242,24 +298,22 @@ AbstractAverageHistogram::resampleDoubleBinWidth(bool bIntegerBins) const
         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)
         {
@@ -277,7 +331,7 @@ AbstractAverageHistogram::resampleDoubleBinWidth(bool bIntegerBins) const
 AverageHistogramPointer
 AbstractAverageHistogram::clone() const
 {
-    AverageHistogramPointer dest(new internal::StaticAverageHistogram());
+    AverageHistogramPointer dest(new StaticAverageHistogram());
     copyContents(this, dest.get());
     return dest;
 }
@@ -286,88 +340,130 @@ AbstractAverageHistogram::clone() const
 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::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];
+        }
     }
 }
 
 
 /********************************************************************
- * 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));
+    }
 }
 
 
@@ -380,32 +476,28 @@ BasicAverageHistogramModule::frameStarted(const AnalysisDataFrameHeader & /*head
 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)));
+        }
     }
 }
 
@@ -414,6 +506,43 @@ BasicAverageHistogramModule::dataFinished()
  * BasicHistogramImpl
  */
 
+/*! \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
+ */
+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);
+        // Virtual only for simplicity.
+        virtual ~BasicHistogramImpl();
+
+        /*! \brief
+         * (Re)initializes the histogram from settings.
+         */
+        void init(const AnalysisHistogramSettings &settings);
+
+        //! Storage implementation object.
+        AnalysisDataStorage                  storage_;
+        //! Settings for the histogram object.
+        AnalysisHistogramSettings            settings_;
+        //! Averager module.
+        BasicAverageHistogramModulePointer   averager_;
+};
+
 BasicHistogramImpl::BasicHistogramImpl()
     : averager_(new BasicAverageHistogramModule())
 {
@@ -437,32 +566,44 @@ void BasicHistogramImpl::init(const AnalysisHistogramSettings &settings)
     averager_->init(settings);
 }
 
-
-void
-BasicHistogramImpl::initFrame(AnalysisDataStorageFrame *frame)
-{
-    for (int i = 0; i < frame->columnCount(); ++i)
-    {
-        frame->setValue(i, 0.0);
-    }
-}
-
-} // namespace internal
+}   // 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))
 {
 }
 
@@ -492,42 +633,63 @@ AnalysisDataSimpleHistogramModule::settings() const
 }
 
 
+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;
+            }
         }
     }
 }
@@ -536,14 +698,27 @@ AnalysisDataSimpleHistogramModule::pointsAdded(const AnalysisDataPointSetRef &po
 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();
 }
 
 
@@ -565,15 +740,37 @@ AnalysisDataSimpleHistogramModule::requestStorageInternal(int nframes)
  * 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))
 {
 }
 
@@ -603,28 +800,45 @@ AnalysisDataWeightedHistogramModule::settings() const
 }
 
 
+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();
 }
 
 
@@ -638,11 +852,11 @@ AnalysisDataWeightedHistogramModule::pointsAdded(const AnalysisDataPointSetRef &
     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);
         }
     }
 }
@@ -651,14 +865,27 @@ AnalysisDataWeightedHistogramModule::pointsAdded(const AnalysisDataPointSetRef &
 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();
 }
 
 
@@ -680,7 +907,23 @@ AnalysisDataWeightedHistogramModule::requestStorageInternal(int nframes)
  * 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);
 }
@@ -688,9 +931,8 @@ AnalysisDataBinAverageModule::AnalysisDataBinAverageModule()
 
 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());
@@ -705,24 +947,36 @@ AnalysisDataBinAverageModule::~AnalysisDataBinAverageModule()
 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());
+    }
 }
 
 
@@ -742,13 +996,11 @@ AnalysisDataBinAverageModule::pointsAdded(const AnalysisDataPointSetRef &points)
     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;
     }
 }
 
@@ -762,15 +1014,15 @@ AnalysisDataBinAverageModule::frameFinished(const AnalysisDataFrameHeader & /*he
 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();