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::flags() const
 621 {
 622     return efAllowMulticolumn | efAllowMultipoint | efAllowMultipleDataSets;
 623 }
 624
 625
 626 void
 627 AnalysisDataSimpleHistogramModule::dataStarted(AbstractAnalysisData *data)
 628 {
 629     addModule(impl_->averager_);
 630     setDataSetCount(data->dataSetCount());
 631     for (int i = 0; i < data->dataSetCount(); ++i)
 632     {
 633         setColumnCount(i, settings().binCount());
 634     }
 635     notifyDataStart();
 636     impl_->storage_.startDataStorage(this);
 637 }
 638
 639
 640 void
 641 AnalysisDataSimpleHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
 642 {
 643     AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
 644     impl_->initFrame(dataSetCount(), &frame);
 645 }
 646
 647
 648 void
 649 AnalysisDataSimpleHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
 650 {
 651     AnalysisDataStorageFrame &frame =
 652         impl_->storage_.currentFrame(points.frameIndex());
 653     frame.selectDataSet(points.dataSetIndex());
 654     for (int i = 0; i < points.columnCount(); ++i)
 655     {
 656         int bin = settings().findBin(points.y(i));
 657         if (bin != -1)
 658         {
 659             frame.value(bin) += 1;
 660         }
 661     }
 662 }
 663
 664
 665 void
 666 AnalysisDataSimpleHistogramModule::frameFinished(const AnalysisDataFrameHeader &header)
 667 {
 668     impl_->storage_.finishFrame(header.index());
 669 }
 670
 671
 672 void
 673 AnalysisDataSimpleHistogramModule::dataFinished()
 674 {
 675     notifyDataFinish();
 676 }
 677
 678
 679 AnalysisDataFrameRef
 680 AnalysisDataSimpleHistogramModule::tryGetDataFrameInternal(int index) const
 681 {
 682     return impl_->storage_.tryGetDataFrame(index);
 683 }
 684
 685
 686 bool
 687 AnalysisDataSimpleHistogramModule::requestStorageInternal(int nframes)
 688 {
 689     return impl_->storage_.requestStorage(nframes);
 690 }
 691
 692
 693 /********************************************************************
 694  * AnalysisDataWeightedHistogramModule
 695  */
 696
 697 AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule()
 698     : impl_(new internal::BasicHistogramImpl())
 699 {
 700 }
 701
 702
 703 AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule(
 704         const AnalysisHistogramSettings &settings)
 705     : impl_(new internal::BasicHistogramImpl(settings))
 706 {
 707 }
 708
 709
 710 AnalysisDataWeightedHistogramModule::~AnalysisDataWeightedHistogramModule()
 711 {
 712 }
 713
 714
 715 void AnalysisDataWeightedHistogramModule::init(const AnalysisHistogramSettings &settings)
 716 {
 717     impl_->init(settings);
 718 }
 719
 720
 721 AbstractAverageHistogram &
 722 AnalysisDataWeightedHistogramModule::averager()
 723 {
 724     return *impl_->averager_;
 725 }
 726
 727
 728 const AnalysisHistogramSettings &
 729 AnalysisDataWeightedHistogramModule::settings() const
 730 {
 731     return impl_->settings_;
 732 }
 733
 734
 735 int
 736 AnalysisDataWeightedHistogramModule::flags() const
 737 {
 738     return efAllowMulticolumn | efAllowMultipoint | efAllowMultipleDataSets;
 739 }
 740
 741
 742 void
 743 AnalysisDataWeightedHistogramModule::dataStarted(AbstractAnalysisData *data)
 744 {
 745     addModule(impl_->averager_);
 746     setDataSetCount(data->dataSetCount());
 747     for (int i = 0; i < data->dataSetCount(); ++i)
 748     {
 749         setColumnCount(i, settings().binCount());
 750     }
 751     notifyDataStart();
 752     impl_->storage_.startDataStorage(this);
 753 }
 754
 755
 756 void
 757 AnalysisDataWeightedHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
 758 {
 759     AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
 760     impl_->initFrame(dataSetCount(), &frame);
 761 }
 762
 763
 764 void
 765 AnalysisDataWeightedHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
 766 {
 767     if (points.firstColumn() != 0 || points.columnCount() < 2)
 768     {
 769         GMX_THROW(APIError("Invalid data layout"));
 770     }
 771     int bin = settings().findBin(points.y(0));
 772     if (bin != -1)
 773     {
 774         AnalysisDataStorageFrame &frame =
 775             impl_->storage_.currentFrame(points.frameIndex());
 776         frame.selectDataSet(points.dataSetIndex());
 777         for (int i = 1; i < points.columnCount(); ++i)
 778         {
 779             frame.value(bin) += points.y(i);
 780         }
 781     }
 782 }
 783
 784
 785 void
 786 AnalysisDataWeightedHistogramModule::frameFinished(const AnalysisDataFrameHeader &header)
 787 {
 788     impl_->storage_.finishFrame(header.index());
 789 }
 790
 791
 792 void
 793 AnalysisDataWeightedHistogramModule::dataFinished()
 794 {
 795     notifyDataFinish();
 796 }
 797
 798
 799 AnalysisDataFrameRef
 800 AnalysisDataWeightedHistogramModule::tryGetDataFrameInternal(int index) const
 801 {
 802     return impl_->storage_.tryGetDataFrame(index);
 803 }
 804
 805
 806 bool
 807 AnalysisDataWeightedHistogramModule::requestStorageInternal(int nframes)
 808 {
 809     return impl_->storage_.requestStorage(nframes);
 810 }
 811
 812
 813 /********************************************************************
 814  * AnalysisDataBinAverageModule
 815  */
 816
 817 class AnalysisDataBinAverageModule::Impl
 818 {
 819     public:
 820         Impl() {}
 821         explicit Impl(const AnalysisHistogramSettings &settings)
 822             : settings_(settings)
 823         {
 824         }
 825
 826         //! Histogram settings.
 827         AnalysisHistogramSettings               settings_;
 828         //! Averaging helper objects for each input data set.
 829         std::vector<AnalysisDataFrameAverager>  averagers_;
 830 };
 831
 832 AnalysisDataBinAverageModule::AnalysisDataBinAverageModule()
 833     : impl_(new Impl())
 834 {
 835     setColumnCount(3);
 836 }
 837
 838
 839 AnalysisDataBinAverageModule::AnalysisDataBinAverageModule(
 840         const AnalysisHistogramSettings &settings)
 841     : impl_(new Impl(settings))
 842 {
 843     setRowCount(settings.binCount());
 844     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
 845              settings.binWidth());
 846 }
 847
 848
 849 AnalysisDataBinAverageModule::~AnalysisDataBinAverageModule()
 850 {
 851 }
 852
 853
 854 void
 855 AnalysisDataBinAverageModule::init(const AnalysisHistogramSettings &settings)
 856 {
 857     impl_->settings_ = settings;
 858     setRowCount(settings.binCount());
 859     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
 860              settings.binWidth());
 861 }
 862
 863
 864 const AnalysisHistogramSettings &
 865 AnalysisDataBinAverageModule::settings() const
 866 {
 867     return impl_->settings_;
 868 }
 869
 870
 871 int
 872 AnalysisDataBinAverageModule::flags() const
 873 {
 874     return efAllowMulticolumn | efAllowMultipoint | efAllowMultipleDataSets;
 875 }
 876
 877
 878 void
 879 AnalysisDataBinAverageModule::dataStarted(AbstractAnalysisData *data)
 880 {
 881     setColumnCount(data->dataSetCount());
 882     impl_->averagers_.resize(data->dataSetCount());
 883     for (int i = 0; i < data->dataSetCount(); ++i)
 884     {
 885         impl_->averagers_[i].setColumnCount(rowCount());
 886     }
 887 }
 888
 889
 890 void
 891 AnalysisDataBinAverageModule::frameStarted(const AnalysisDataFrameHeader & /*header*/)
 892 {
 893 }
 894
 895
 896 void
 897 AnalysisDataBinAverageModule::pointsAdded(const AnalysisDataPointSetRef &points)
 898 {
 899     if (points.firstColumn() != 0 || points.columnCount() < 2)
 900     {
 901         GMX_THROW(APIError("Invalid data layout"));
 902     }
 903     int bin = settings().findBin(points.y(0));
 904     if (bin != -1)
 905     {
 906         AnalysisDataFrameAverager &averager = impl_->averagers_[points.dataSetIndex()];
 907         for (int i = 1; i < points.columnCount(); ++i)
 908         {
 909             averager.addValue(bin, points.y(i));
 910         }
 911     }
 912 }
 913
 914
 915 void
 916 AnalysisDataBinAverageModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
 917 {
 918 }
 919
 920
 921 void
 922 AnalysisDataBinAverageModule::dataFinished()
 923 {
 924     allocateValues();
 925     for (int i = 0; i < columnCount(); ++i)
 926     {
 927         AnalysisDataFrameAverager &averager = impl_->averagers_[i];
 928         averager.finish();
 929         for (int j = 0; j < rowCount(); ++j)
 930         {
 931             value(j, i).setValue(averager.average(j),
 932                                  std::sqrt(averager.variance(j)));
 933         }
 934     }
 935     valuesReady();
 936 }
 937
 938 } // namespace gmx