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