src/gromacs/analysisdata/modules/histogram.cpp

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