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