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