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