src/gromacs/analysisdata/modules/histogram.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2010,2011,2012,2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Implements classes in histogram.h.
  38  *
  39  * \author Teemu Murtola <teemu.murtola@gmail.com>
  40  * \ingroup module_analysisdata
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "histogram.h"
  45
  46 #include <cmath>
  47
  48 #include <limits>
  49 #include <vector>
  50
  51 #include "gromacs/analysisdata/dataframe.h"
  52 #include "gromacs/analysisdata/datastorage.h"
  53 #include "gromacs/utility/exceptions.h"
  54 #include "gromacs/utility/gmxassert.h"
  55
  56 #include "frameaverager.h"
  57
  58 namespace
  59 {
  60
  61 //! Value used to signify that a real-valued histogram setting is not set.
  62 const real UNDEFINED = std::numeric_limits<real>::max();
  63 //! Checks whether \p value is defined.
  64 bool isDefined(real value)
  65 {
  66     return value != UNDEFINED;
  67 }
  68
  69 } // namespace
  70
  71 namespace gmx
  72 {
  73
  74 /********************************************************************
  75  * AnalysisHistogramSettingsInitializer
  76  */
  77
  78 AnalysisHistogramSettingsInitializer::AnalysisHistogramSettingsInitializer()
  79     : min_(UNDEFINED), max_(UNDEFINED), binWidth_(UNDEFINED),
  80       binCount_(0), bIntegerBins_(false), bRoundRange_(false),
  81       bIncludeAll_(false)
  82 {
  83 }
  84
  85
  86 /********************************************************************
  87  * AnalysisHistogramSettings
  88  */
  89
  90 AnalysisHistogramSettings::AnalysisHistogramSettings()
  91     : firstEdge_(0.0), lastEdge_(0.0), binWidth_(0.0), inverseBinWidth_(0.0),
  92       binCount_(0), bAll_(false)
  93 {
  94 }
  95
  96
  97 AnalysisHistogramSettings::AnalysisHistogramSettings(
  98         const AnalysisHistogramSettingsInitializer &settings)
  99 {
 100     GMX_RELEASE_ASSERT(isDefined(settings.min_),
 101                        "Histogram start value must be defined");
 102     GMX_RELEASE_ASSERT(!isDefined(settings.max_) || settings.max_ > settings.min_,
 103                        "Histogram end value must be larger than start value");
 104     GMX_RELEASE_ASSERT(!isDefined(settings.binWidth_) || settings.binWidth_ > 0.0,
 105                        "Histogram bin width must be positive");
 106     GMX_RELEASE_ASSERT(settings.binCount_ >= 0,
 107                        "Histogram bin count must be positive");
 108
 109     if (!isDefined(settings.max_))
 110     {
 111         GMX_RELEASE_ASSERT(isDefined(settings.binWidth_) && settings.binCount_ > 0,
 112                            "Not all required values provided");
 113         GMX_RELEASE_ASSERT(!settings.bRoundRange_,
 114                            "Rounding only supported for min/max ranges");
 115
 116         firstEdge_ = settings.min_;
 117         binCount_  = settings.binCount_;
 118         binWidth_  = settings.binWidth_;
 119         if (settings.bIntegerBins_)
 120         {
 121             firstEdge_ -= 0.5 * binWidth_;
 122         }
 123         lastEdge_ = firstEdge_ + binCount_ * binWidth_;
 124     }
 125     else
 126     {
 127         GMX_RELEASE_ASSERT(!(isDefined(settings.binWidth_) && settings.binCount_ > 0),
 128                            "Conflicting histogram bin specifications");
 129         GMX_RELEASE_ASSERT(isDefined(settings.binWidth_) || settings.binCount_ > 0,
 130                            "Not all required values provided");
 131
 132         if (settings.bRoundRange_)
 133         {
 134             GMX_RELEASE_ASSERT(!settings.bIntegerBins_,
 135                                "Rounding and integer bins cannot be combined");
 136             GMX_RELEASE_ASSERT(isDefined(settings.binWidth_),
 137                                "Rounding only makes sense with defined binwidth");
 138             binWidth_  = settings.binWidth_;
 139             firstEdge_ = binWidth_ * floor(settings.min_ / binWidth_);
 140             lastEdge_  = binWidth_ * ceil(settings.max_ / binWidth_);
 141             binCount_  = static_cast<int>((lastEdge_ - firstEdge_) / binWidth_ + 0.5);
 142         }
 143         else
 144         {
 145             firstEdge_     = settings.min_;
 146             lastEdge_      = settings.max_;
 147             if (settings.binCount_ > 0)
 148             {
 149                 binCount_ = settings.binCount_;
 150                 if (settings.bIntegerBins_)
 151                 {
 152                     GMX_RELEASE_ASSERT(settings.binCount_ > 1,
 153                                        "Bin count must be at least two with integer bins");
 154                     binWidth_   = (lastEdge_ - firstEdge_) / (binCount_ - 1);
 155                     firstEdge_ -= 0.5 * binWidth_;
 156                     lastEdge_  += 0.5 * binWidth_;
 157                 }
 158                 else
 159                 {
 160                     binWidth_ = (lastEdge_ - firstEdge_) / binCount_;
 161                 }
 162             }
 163             else
 164             {
 165                 binWidth_ = settings.binWidth_;
 166                 binCount_ = static_cast<int>((lastEdge_ - firstEdge_) / binWidth_ + 0.5);
 167                 if (settings.bIntegerBins_)
 168                 {
 169                     firstEdge_ -= 0.5 * binWidth_;
 170                     ++binCount_;
 171                 }
 172                 lastEdge_ = firstEdge_ + binCount_ * binWidth_;
 173             }
 174         }
 175     }
 176
 177     inverseBinWidth_ = 1.0 / binWidth_;
 178     bAll_            = settings.bIncludeAll_;
 179 }
 180
 181
 182 int
 183 AnalysisHistogramSettings::findBin(real y) const
 184 {
 185     if (y < firstEdge_)
 186     {
 187         return bAll_ ? 0 : -1;
 188     }
 189     int bin = static_cast<int>((y - firstEdge_) * inverseBinWidth_);
 190     if (bin >= binCount_)
 191     {
 192         return bAll_ ? binCount_ - 1 : -1;
 193     }
 194     return bin;
 195 }
 196
 197
 198 /********************************************************************
 199  * StaticAverageHistogram
 200  */
 201
 202 namespace
 203 {
 204
 205 /*! \brief
 206  * Represents copies of average histograms.
 207  *
 208  * Methods in AbstractAverageHistogram that return new histogram instances
 209  * return objects of this class.
 210  * Initialization of values is handled in those methods.
 211  *
 212  * \ingroup module_analysisdata
 213  */
 214 class StaticAverageHistogram : public AbstractAverageHistogram
 215 {
 216     public:
 217         StaticAverageHistogram();
 218         //! Creates an average histogram module with defined bin parameters.
 219         explicit StaticAverageHistogram(const AnalysisHistogramSettings &settings);
 220
 221         // Copy and assign disallowed by base.
 222 };
 223
 224 StaticAverageHistogram::StaticAverageHistogram()
 225 {
 226 }
 227
 228
 229 StaticAverageHistogram::StaticAverageHistogram(
 230         const AnalysisHistogramSettings &settings)
 231     : AbstractAverageHistogram(settings)
 232 {
 233 }
 234
 235 }   // namespace
 236
 237
 238 /********************************************************************
 239  * AbstractAverageHistogram
 240  */
 241
 242 AbstractAverageHistogram::AbstractAverageHistogram()
 243 {
 244 }
 245
 246
 247 AbstractAverageHistogram::AbstractAverageHistogram(
 248         const AnalysisHistogramSettings &settings)
 249     : settings_(settings)
 250 {
 251     setRowCount(settings.binCount());
 252     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
 253              settings.binWidth());
 254 }
 255
 256
 257 AbstractAverageHistogram::~AbstractAverageHistogram()
 258 {
 259 }
 260
 261
 262 void
 263 AbstractAverageHistogram::init(const AnalysisHistogramSettings &settings)
 264 {
 265     settings_ = settings;
 266     setRowCount(settings.binCount());
 267     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
 268              settings.binWidth());
 269 }
 270
 271
 272 AverageHistogramPointer
 273 AbstractAverageHistogram::resampleDoubleBinWidth(bool bIntegerBins) const
 274 {
 275     int nbins;
 276     if (bIntegerBins)
 277     {
 278         nbins = (rowCount() + 1) / 2;
 279     }
 280     else
 281     {
 282         nbins = rowCount() / 2;
 283     }
 284
 285     AverageHistogramPointer dest(
 286             new StaticAverageHistogram(
 287                     histogramFromBins(xstart(), nbins, 2*xstep())
 288                         .integerBins(bIntegerBins)));
 289     dest->setColumnCount(columnCount());
 290     dest->allocateValues();
 291
 292     int  i, j;
 293     for (i = j = 0; i < nbins; ++i)
 294     {
 295         const bool bFirstHalfBin = (bIntegerBins && i == 0);
 296         for (int c = 0; c < columnCount(); ++c)
 297         {
 298             real  v1, v2;
 299             real  e1, e2;
 300             if (bFirstHalfBin)
 301             {
 302                 v1 = value(0, c).value();
 303                 e1 = value(0, c).error();
 304                 v2 = 0;
 305                 e2 = 0;
 306             }
 307             else
 308             {
 309                 v1 = value(j, c).value();
 310                 e1 = value(j, c).error();
 311                 v2 = value(j + 1, c).value();
 312                 e2 = value(j + 1, c).error();
 313             }
 314             dest->value(i, c).setValue(v1 + v2, std::sqrt(e1 * e1 + e2 * e2));
 315         }
 316         if (bFirstHalfBin)
 317         {
 318             ++j;
 319         }
 320         else
 321         {
 322             j += 2;
 323         }
 324     }
 325     return dest;
 326 }
 327
 328
 329 AverageHistogramPointer
 330 AbstractAverageHistogram::clone() const
 331 {
 332     AverageHistogramPointer dest(new StaticAverageHistogram());
 333     copyContents(this, dest.get());
 334     return dest;
 335 }
 336
 337
 338 void
 339 AbstractAverageHistogram::normalizeProbability()
 340 {
 341     for (int c = 0; c < columnCount(); ++c)
 342     {
 343         real sum = 0;
 344         for (int i = 0; i < rowCount(); ++i)
 345         {
 346             sum += value(i, c).value();
 347         }
 348         if (sum > 0.0)
 349         {
 350             scaleSingle(c, 1.0 / (sum * xstep()));
 351         }
 352     }
 353 }
 354
 355
 356 void
 357 AbstractAverageHistogram::scaleSingle(int index, real factor)
 358 {
 359     for (int i = 0; i < rowCount(); ++i)
 360     {
 361         value(i, index).value() *= factor;
 362         value(i, index).error() *= factor;
 363     }
 364 }
 365
 366
 367 void
 368 AbstractAverageHistogram::scaleAll(real factor)
 369 {
 370     for (int i = 0; i < columnCount(); ++i)
 371     {
 372         scaleSingle(i, factor);
 373     }
 374 }
 375
 376
 377 void
 378 AbstractAverageHistogram::scaleAllByVector(real factor[])
 379 {
 380     for (int c = 0; c < columnCount(); ++c)
 381     {
 382         for (int i = 0; i < rowCount(); ++i)
 383         {
 384             value(i, c).value() *= factor[i];
 385             value(i, c).error() *= factor[i];
 386         }
 387     }
 388 }
 389
 390
 391 /********************************************************************
 392  * BasicAverageHistogramModule
 393  */
 394
 395 namespace internal
 396 {
 397
 398 /*! \internal
 399  * \brief
 400  * Implements average histogram module that averages per-frame histograms.
 401  *
 402  * This class is used for accumulating average histograms in per-frame
 403  * histogram modules (those that use BasicHistogramImpl as their implementation
 404  * class).
 405  * There are two columns, first for the average and second for standard
 406  * deviation.
 407  *
 408  * \ingroup module_analysisdata
 409  */
 410 class BasicAverageHistogramModule : public AbstractAverageHistogram,
 411                                     public AnalysisDataModuleSerial
 412 {
 413     public:
 414         BasicAverageHistogramModule();
 415         //! Creates an average histogram module with defined bin parameters.
 416         explicit BasicAverageHistogramModule(const AnalysisHistogramSettings &settings);
 417
 418         using AbstractAverageHistogram::init;
 419
 420         virtual int flags() const;
 421
 422         virtual void dataStarted(AbstractAnalysisData *data);
 423         virtual void frameStarted(const AnalysisDataFrameHeader &header);
 424         virtual void pointsAdded(const AnalysisDataPointSetRef &points);
 425         virtual void frameFinished(const AnalysisDataFrameHeader &header);
 426         virtual void dataFinished();
 427
 428     private:
 429         //! Averaging helper objects for each input data set.
 430         std::vector<AnalysisDataFrameAverager> averagers_;
 431
 432         // Copy and assign disallowed by base.
 433 };
 434
 435 BasicAverageHistogramModule::BasicAverageHistogramModule()
 436 {
 437 }
 438
 439
 440 BasicAverageHistogramModule::BasicAverageHistogramModule(
 441         const AnalysisHistogramSettings &settings)
 442     : AbstractAverageHistogram(settings)
 443 {
 444 }
 445
 446
 447 int
 448 BasicAverageHistogramModule::flags() const
 449 {
 450     return efAllowMulticolumn | efAllowMultipleDataSets;
 451 }
 452
 453
 454 void
 455 BasicAverageHistogramModule::dataStarted(AbstractAnalysisData *data)
 456 {
 457     setColumnCount(data->dataSetCount());
 458     averagers_.resize(data->dataSetCount());
 459     for (int i = 0; i < data->dataSetCount(); ++i)
 460     {
 461         GMX_RELEASE_ASSERT(rowCount() == data->columnCount(i),
 462                            "Inconsistent data sizes, something is wrong in the initialization");
 463         averagers_[i].setColumnCount(data->columnCount(i));
 464     }
 465 }
 466
 467
 468 void
 469 BasicAverageHistogramModule::frameStarted(const AnalysisDataFrameHeader & /*header*/)
 470 {
 471 }
 472
 473
 474 void
 475 BasicAverageHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
 476 {
 477     averagers_[points.dataSetIndex()].addPoints(points);
 478 }
 479
 480
 481 void
 482 BasicAverageHistogramModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
 483 {
 484 }
 485
 486
 487 void
 488 BasicAverageHistogramModule::dataFinished()
 489 {
 490     allocateValues();
 491     for (int i = 0; i < columnCount(); ++i)
 492     {
 493         averagers_[i].finish();
 494         for (int j = 0; j < rowCount(); ++j)
 495         {
 496             value(j, i).setValue(averagers_[i].average(j),
 497                                  std::sqrt(averagers_[i].variance(j)));
 498         }
 499     }
 500 }
 501
 502
 503 /********************************************************************
 504  * BasicHistogramImpl
 505  */
 506
 507 /*! \internal \brief
 508  * Private implementation class for AnalysisDataSimpleHistogramModule and
 509  * AnalysisDataWeightedHistogramModule.
 510  *
 511  * \ingroup module_analysisdata
 512  */
 513 class BasicHistogramImpl
 514 {
 515     public:
 516         //! Smart pointer to manage an BasicAverageHistogramModule object.
 517         typedef boost::shared_ptr<BasicAverageHistogramModule>
 518             BasicAverageHistogramModulePointer;
 519
 520         BasicHistogramImpl();
 521         //! Creates an histogram impl with defined bin parameters.
 522         explicit BasicHistogramImpl(const AnalysisHistogramSettings &settings);
 523         ~BasicHistogramImpl();
 524
 525         /*! \brief
 526          * (Re)initializes the histogram from settings.
 527          */
 528         void init(const AnalysisHistogramSettings &settings);
 529         /*! \brief
 530          * Initializes data storage frame when a new frame starts.
 531          */
 532         void initFrame(int dataSetCount, AnalysisDataStorageFrame *frame);
 533
 534         //! Storage implementation object.
 535         AnalysisDataStorage                  storage_;
 536         //! Settings for the histogram object.
 537         AnalysisHistogramSettings            settings_;
 538         //! Averager module.
 539         BasicAverageHistogramModulePointer   averager_;
 540 };
 541
 542 BasicHistogramImpl::BasicHistogramImpl()
 543     : averager_(new BasicAverageHistogramModule())
 544 {
 545 }
 546
 547
 548 BasicHistogramImpl::BasicHistogramImpl(const AnalysisHistogramSettings &settings)
 549     : settings_(settings), averager_(new BasicAverageHistogramModule(settings))
 550 {
 551 }
 552
 553
 554 BasicHistogramImpl::~BasicHistogramImpl()
 555 {
 556 }
 557
 558
 559 void BasicHistogramImpl::init(const AnalysisHistogramSettings &settings)
 560 {
 561     settings_ = settings;
 562     averager_->init(settings);
 563 }
 564
 565
 566 void
 567 BasicHistogramImpl::initFrame(int dataSetCount, AnalysisDataStorageFrame *frame)
 568 {
 569     for (int s = 0; s < dataSetCount; ++s)
 570     {
 571         frame->selectDataSet(s);
 572         for (int i = 0; i < frame->columnCount(); ++i)
 573         {
 574             frame->setValue(i, 0.0);
 575         }
 576     }
 577     frame->selectDataSet(0);
 578 }
 579
 580 }   // namespace internal
 581
 582
 583 /********************************************************************
 584  * AnalysisDataSimpleHistogramModule
 585  */
 586
 587 AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule()
 588     : impl_(new internal::BasicHistogramImpl())
 589 {
 590 }
 591
 592
 593 AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule(
 594         const AnalysisHistogramSettings &settings)
 595     : impl_(new internal::BasicHistogramImpl(settings))
 596 {
 597 }
 598
 599
 600 AnalysisDataSimpleHistogramModule::~AnalysisDataSimpleHistogramModule()
 601 {
 602 }
 603
 604
 605 void AnalysisDataSimpleHistogramModule::init(const AnalysisHistogramSettings &settings)
 606 {
 607     impl_->init(settings);
 608 }
 609
 610
 611 AbstractAverageHistogram &
 612 AnalysisDataSimpleHistogramModule::averager()
 613 {
 614     return *impl_->averager_;
 615 }
 616
 617
 618 const AnalysisHistogramSettings &
 619 AnalysisDataSimpleHistogramModule::settings() const
 620 {
 621     return impl_->settings_;
 622 }
 623
 624
 625 int
 626 AnalysisDataSimpleHistogramModule::frameCount() const
 627 {
 628     return impl_->storage_.frameCount();
 629 }
 630
 631
 632 int
 633 AnalysisDataSimpleHistogramModule::flags() const
 634 {
 635     return efAllowMulticolumn | efAllowMultipoint | efAllowMissing
 636            | efAllowMultipleDataSets;
 637 }
 638
 639
 640 bool
 641 AnalysisDataSimpleHistogramModule::parallelDataStarted(
 642         AbstractAnalysisData              *data,
 643         const AnalysisDataParallelOptions &options)
 644 {
 645     addModule(impl_->averager_);
 646     setDataSetCount(data->dataSetCount());
 647     for (int i = 0; i < data->dataSetCount(); ++i)
 648     {
 649         setColumnCount(i, settings().binCount());
 650     }
 651     impl_->storage_.startParallelDataStorage(this, &moduleManager(), options);
 652     return true;
 653 }
 654
 655
 656 void
 657 AnalysisDataSimpleHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
 658 {
 659     AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
 660     impl_->initFrame(dataSetCount(), &frame);
 661 }
 662
 663
 664 void
 665 AnalysisDataSimpleHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
 666 {
 667     AnalysisDataStorageFrame &frame =
 668         impl_->storage_.currentFrame(points.frameIndex());
 669     frame.selectDataSet(points.dataSetIndex());
 670     for (int i = 0; i < points.columnCount(); ++i)
 671     {
 672         if (points.present(i))
 673         {
 674             const int bin = settings().findBin(points.y(i));
 675             if (bin != -1)
 676             {
 677                 frame.value(bin) += 1;
 678             }
 679         }
 680     }
 681 }
 682
 683
 684 void
 685 AnalysisDataSimpleHistogramModule::frameFinished(const AnalysisDataFrameHeader &header)
 686 {
 687     impl_->storage_.finishFrame(header.index());
 688 }
 689
 690
 691 void
 692 AnalysisDataSimpleHistogramModule::dataFinished()
 693 {
 694     impl_->storage_.finishDataStorage();
 695 }
 696
 697
 698 AnalysisDataFrameRef
 699 AnalysisDataSimpleHistogramModule::tryGetDataFrameInternal(int index) const
 700 {
 701     return impl_->storage_.tryGetDataFrame(index);
 702 }
 703
 704
 705 bool
 706 AnalysisDataSimpleHistogramModule::requestStorageInternal(int nframes)
 707 {
 708     return impl_->storage_.requestStorage(nframes);
 709 }
 710
 711
 712 /********************************************************************
 713  * AnalysisDataWeightedHistogramModule
 714  */
 715
 716 AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule()
 717     : impl_(new internal::BasicHistogramImpl())
 718 {
 719 }
 720
 721
 722 AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule(
 723         const AnalysisHistogramSettings &settings)
 724     : impl_(new internal::BasicHistogramImpl(settings))
 725 {
 726 }
 727
 728
 729 AnalysisDataWeightedHistogramModule::~AnalysisDataWeightedHistogramModule()
 730 {
 731 }
 732
 733
 734 void AnalysisDataWeightedHistogramModule::init(const AnalysisHistogramSettings &settings)
 735 {
 736     impl_->init(settings);
 737 }
 738
 739
 740 AbstractAverageHistogram &
 741 AnalysisDataWeightedHistogramModule::averager()
 742 {
 743     return *impl_->averager_;
 744 }
 745
 746
 747 const AnalysisHistogramSettings &
 748 AnalysisDataWeightedHistogramModule::settings() const
 749 {
 750     return impl_->settings_;
 751 }
 752
 753
 754 int
 755 AnalysisDataWeightedHistogramModule::frameCount() const
 756 {
 757     return impl_->storage_.frameCount();
 758 }
 759
 760
 761 int
 762 AnalysisDataWeightedHistogramModule::flags() const
 763 {
 764     return efAllowMulticolumn | efAllowMultipoint | efAllowMultipleDataSets;
 765 }
 766
 767
 768 bool
 769 AnalysisDataWeightedHistogramModule::parallelDataStarted(
 770         AbstractAnalysisData              *data,
 771         const AnalysisDataParallelOptions &options)
 772 {
 773     addModule(impl_->averager_);
 774     setDataSetCount(data->dataSetCount());
 775     for (int i = 0; i < data->dataSetCount(); ++i)
 776     {
 777         setColumnCount(i, settings().binCount());
 778     }
 779     impl_->storage_.startParallelDataStorage(this, &moduleManager(), options);
 780     return true;
 781 }
 782
 783
 784 void
 785 AnalysisDataWeightedHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
 786 {
 787     AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
 788     impl_->initFrame(dataSetCount(), &frame);
 789 }
 790
 791
 792 void
 793 AnalysisDataWeightedHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
 794 {
 795     if (points.firstColumn() != 0 || points.columnCount() < 2)
 796     {
 797         GMX_THROW(APIError("Invalid data layout"));
 798     }
 799     int bin = settings().findBin(points.y(0));
 800     if (bin != -1)
 801     {
 802         AnalysisDataStorageFrame &frame =
 803             impl_->storage_.currentFrame(points.frameIndex());
 804         frame.selectDataSet(points.dataSetIndex());
 805         for (int i = 1; i < points.columnCount(); ++i)
 806         {
 807             frame.value(bin) += points.y(i);
 808         }
 809     }
 810 }
 811
 812
 813 void
 814 AnalysisDataWeightedHistogramModule::frameFinished(const AnalysisDataFrameHeader &header)
 815 {
 816     impl_->storage_.finishFrame(header.index());
 817 }
 818
 819
 820 void
 821 AnalysisDataWeightedHistogramModule::dataFinished()
 822 {
 823     impl_->storage_.finishDataStorage();
 824 }
 825
 826
 827 AnalysisDataFrameRef
 828 AnalysisDataWeightedHistogramModule::tryGetDataFrameInternal(int index) const
 829 {
 830     return impl_->storage_.tryGetDataFrame(index);
 831 }
 832
 833
 834 bool
 835 AnalysisDataWeightedHistogramModule::requestStorageInternal(int nframes)
 836 {
 837     return impl_->storage_.requestStorage(nframes);
 838 }
 839
 840
 841 /********************************************************************
 842  * AnalysisDataBinAverageModule
 843  */
 844
 845 class AnalysisDataBinAverageModule::Impl
 846 {
 847     public:
 848         Impl() {}
 849         explicit Impl(const AnalysisHistogramSettings &settings)
 850             : settings_(settings)
 851         {
 852         }
 853
 854         //! Histogram settings.
 855         AnalysisHistogramSettings               settings_;
 856         //! Averaging helper objects for each input data set.
 857         std::vector<AnalysisDataFrameAverager>  averagers_;
 858 };
 859
 860 AnalysisDataBinAverageModule::AnalysisDataBinAverageModule()
 861     : impl_(new Impl())
 862 {
 863     setColumnCount(3);
 864 }
 865
 866
 867 AnalysisDataBinAverageModule::AnalysisDataBinAverageModule(
 868         const AnalysisHistogramSettings &settings)
 869     : impl_(new Impl(settings))
 870 {
 871     setRowCount(settings.binCount());
 872     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
 873              settings.binWidth());
 874 }
 875
 876
 877 AnalysisDataBinAverageModule::~AnalysisDataBinAverageModule()
 878 {
 879 }
 880
 881
 882 void
 883 AnalysisDataBinAverageModule::init(const AnalysisHistogramSettings &settings)
 884 {
 885     impl_->settings_ = settings;
 886     setRowCount(settings.binCount());
 887     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
 888              settings.binWidth());
 889 }
 890
 891
 892 const AnalysisHistogramSettings &
 893 AnalysisDataBinAverageModule::settings() const
 894 {
 895     return impl_->settings_;
 896 }
 897
 898
 899 int
 900 AnalysisDataBinAverageModule::flags() const
 901 {
 902     return efAllowMulticolumn | efAllowMultipoint | efAllowMultipleDataSets;
 903 }
 904
 905
 906 void
 907 AnalysisDataBinAverageModule::dataStarted(AbstractAnalysisData *data)
 908 {
 909     setColumnCount(data->dataSetCount());
 910     impl_->averagers_.resize(data->dataSetCount());
 911     for (int i = 0; i < data->dataSetCount(); ++i)
 912     {
 913         impl_->averagers_[i].setColumnCount(rowCount());
 914     }
 915 }
 916
 917
 918 void
 919 AnalysisDataBinAverageModule::frameStarted(const AnalysisDataFrameHeader & /*header*/)
 920 {
 921 }
 922
 923
 924 void
 925 AnalysisDataBinAverageModule::pointsAdded(const AnalysisDataPointSetRef &points)
 926 {
 927     if (points.firstColumn() != 0 || points.columnCount() < 2)
 928     {
 929         GMX_THROW(APIError("Invalid data layout"));
 930     }
 931     int bin = settings().findBin(points.y(0));
 932     if (bin != -1)
 933     {
 934         AnalysisDataFrameAverager &averager = impl_->averagers_[points.dataSetIndex()];
 935         for (int i = 1; i < points.columnCount(); ++i)
 936         {
 937             averager.addValue(bin, points.y(i));
 938         }
 939     }
 940 }
 941
 942
 943 void
 944 AnalysisDataBinAverageModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
 945 {
 946 }
 947
 948
 949 void
 950 AnalysisDataBinAverageModule::dataFinished()
 951 {
 952     allocateValues();
 953     for (int i = 0; i < columnCount(); ++i)
 954     {
 955         AnalysisDataFrameAverager &averager = impl_->averagers_[i];
 956         averager.finish();
 957         for (int j = 0; j < rowCount(); ++j)
 958         {
 959             value(j, i).setValue(averager.average(j),
 960                                  std::sqrt(averager.variance(j)));
 961         }
 962     }
 963     valuesReady();
 964 }
 965
 966 } // namespace gmx