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 /*! \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         if (sum > 0.0)
 347         {
 348             scaleSingle(c, 1.0 / (sum * xstep()));
 349         }
 350     }
 351 }
 352
 353
 354 void
 355 AbstractAverageHistogram::scaleSingle(int index, real factor)
 356 {
 357     for (int i = 0; i < rowCount(); ++i)
 358     {
 359         value(i, index).value() *= factor;
 360         value(i, index).error() *= factor;
 361     }
 362 }
 363
 364
 365 void
 366 AbstractAverageHistogram::scaleAll(real factor)
 367 {
 368     for (int i = 0; i < columnCount(); ++i)
 369     {
 370         scaleSingle(i, factor);
 371     }
 372 }
 373
 374
 375 void
 376 AbstractAverageHistogram::scaleAllByVector(real factor[])
 377 {
 378     for (int c = 0; c < columnCount(); ++c)
 379     {
 380         for (int i = 0; i < rowCount(); ++i)
 381         {
 382             value(i, c).value() *= factor[i];
 383             value(i, c).error() *= factor[i];
 384         }
 385     }
 386 }
 387
 388
 389 /********************************************************************
 390  * BasicAverageHistogramModule
 391  */
 392
 393 namespace internal
 394 {
 395
 396 /*! \internal
 397  * \brief
 398  * Implements average histogram module that averages per-frame histograms.
 399  *
 400  * This class is used for accumulating average histograms in per-frame
 401  * histogram modules (those that use BasicHistogramImpl as their implementation
 402  * class).
 403  * There are two columns, first for the average and second for standard
 404  * deviation.
 405  *
 406  * \ingroup module_analysisdata
 407  */
 408 class BasicAverageHistogramModule : public AbstractAverageHistogram,
 409                                     public AnalysisDataModuleSerial
 410 {
 411     public:
 412         BasicAverageHistogramModule();
 413         //! Creates an average histogram module with defined bin parameters.
 414         explicit BasicAverageHistogramModule(const AnalysisHistogramSettings &settings);
 415
 416         using AbstractAverageHistogram::init;
 417
 418         virtual int flags() const;
 419
 420         virtual void dataStarted(AbstractAnalysisData *data);
 421         virtual void frameStarted(const AnalysisDataFrameHeader &header);
 422         virtual void pointsAdded(const AnalysisDataPointSetRef &points);
 423         virtual void frameFinished(const AnalysisDataFrameHeader &header);
 424         virtual void dataFinished();
 425
 426     private:
 427         //! Averaging helper objects for each input data set.
 428         std::vector<AnalysisDataFrameAverager> averagers_;
 429
 430         // Copy and assign disallowed by base.
 431 };
 432
 433 BasicAverageHistogramModule::BasicAverageHistogramModule()
 434 {
 435 }
 436
 437
 438 BasicAverageHistogramModule::BasicAverageHistogramModule(
 439         const AnalysisHistogramSettings &settings)
 440     : AbstractAverageHistogram(settings)
 441 {
 442 }
 443
 444
 445 int
 446 BasicAverageHistogramModule::flags() const
 447 {
 448     return efAllowMulticolumn | efAllowMultipleDataSets;
 449 }
 450
 451
 452 void
 453 BasicAverageHistogramModule::dataStarted(AbstractAnalysisData *data)
 454 {
 455     setColumnCount(data->dataSetCount());
 456     averagers_.resize(data->dataSetCount());
 457     for (int i = 0; i < data->dataSetCount(); ++i)
 458     {
 459         GMX_RELEASE_ASSERT(rowCount() == data->columnCount(i),
 460                            "Inconsistent data sizes, something is wrong in the initialization");
 461         averagers_[i].setColumnCount(data->columnCount(i));
 462     }
 463 }
 464
 465
 466 void
 467 BasicAverageHistogramModule::frameStarted(const AnalysisDataFrameHeader & /*header*/)
 468 {
 469 }
 470
 471
 472 void
 473 BasicAverageHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
 474 {
 475     averagers_[points.dataSetIndex()].addPoints(points);
 476 }
 477
 478
 479 void
 480 BasicAverageHistogramModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
 481 {
 482 }
 483
 484
 485 void
 486 BasicAverageHistogramModule::dataFinished()
 487 {
 488     allocateValues();
 489     for (int i = 0; i < columnCount(); ++i)
 490     {
 491         averagers_[i].finish();
 492         for (int j = 0; j < rowCount(); ++j)
 493         {
 494             value(j, i).setValue(averagers_[i].average(j),
 495                                  std::sqrt(averagers_[i].variance(j)));
 496         }
 497     }
 498 }
 499
 500
 501 /********************************************************************
 502  * BasicHistogramImpl
 503  */
 504
 505 /*! \internal \brief
 506  * Private implementation class for AnalysisDataSimpleHistogramModule and
 507  * AnalysisDataWeightedHistogramModule.
 508  *
 509  * \ingroup module_analysisdata
 510  */
 511 class BasicHistogramImpl
 512 {
 513     public:
 514         //! Smart pointer to manage an BasicAverageHistogramModule object.
 515         typedef boost::shared_ptr<BasicAverageHistogramModule>
 516             BasicAverageHistogramModulePointer;
 517
 518         BasicHistogramImpl();
 519         //! Creates an histogram impl with defined bin parameters.
 520         explicit BasicHistogramImpl(const AnalysisHistogramSettings &settings);
 521         ~BasicHistogramImpl();
 522
 523         /*! \brief
 524          * (Re)initializes the histogram from settings.
 525          */
 526         void init(const AnalysisHistogramSettings &settings);
 527         /*! \brief
 528          * Initializes data storage frame when a new frame starts.
 529          */
 530         void initFrame(int dataSetCount, AnalysisDataStorageFrame *frame);
 531
 532         //! Storage implementation object.
 533         AnalysisDataStorage                  storage_;
 534         //! Settings for the histogram object.
 535         AnalysisHistogramSettings            settings_;
 536         //! Averager module.
 537         BasicAverageHistogramModulePointer   averager_;
 538 };
 539
 540 BasicHistogramImpl::BasicHistogramImpl()
 541     : averager_(new BasicAverageHistogramModule())
 542 {
 543 }
 544
 545
 546 BasicHistogramImpl::BasicHistogramImpl(const AnalysisHistogramSettings &settings)
 547     : settings_(settings), averager_(new BasicAverageHistogramModule(settings))
 548 {
 549 }
 550
 551
 552 BasicHistogramImpl::~BasicHistogramImpl()
 553 {
 554 }
 555
 556
 557 void BasicHistogramImpl::init(const AnalysisHistogramSettings &settings)
 558 {
 559     settings_ = settings;
 560     averager_->init(settings);
 561 }
 562
 563
 564 void
 565 BasicHistogramImpl::initFrame(int dataSetCount, AnalysisDataStorageFrame *frame)
 566 {
 567     for (int s = 0; s < dataSetCount; ++s)
 568     {
 569         frame->selectDataSet(s);
 570         for (int i = 0; i < frame->columnCount(); ++i)
 571         {
 572             frame->setValue(i, 0.0);
 573         }
 574     }
 575     frame->selectDataSet(0);
 576 }
 577
 578 }   // namespace internal
 579
 580
 581 /********************************************************************
 582  * AnalysisDataSimpleHistogramModule
 583  */
 584
 585 AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule()
 586     : impl_(new internal::BasicHistogramImpl())
 587 {
 588 }
 589
 590
 591 AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule(
 592         const AnalysisHistogramSettings &settings)
 593     : impl_(new internal::BasicHistogramImpl(settings))
 594 {
 595 }
 596
 597
 598 AnalysisDataSimpleHistogramModule::~AnalysisDataSimpleHistogramModule()
 599 {
 600 }
 601
 602
 603 void AnalysisDataSimpleHistogramModule::init(const AnalysisHistogramSettings &settings)
 604 {
 605     impl_->init(settings);
 606 }
 607
 608
 609 AbstractAverageHistogram &
 610 AnalysisDataSimpleHistogramModule::averager()
 611 {
 612     return *impl_->averager_;
 613 }
 614
 615
 616 const AnalysisHistogramSettings &
 617 AnalysisDataSimpleHistogramModule::settings() const
 618 {
 619     return impl_->settings_;
 620 }
 621
 622
 623 int
 624 AnalysisDataSimpleHistogramModule::frameCount() const
 625 {
 626     return impl_->storage_.frameCount();
 627 }
 628
 629
 630 int
 631 AnalysisDataSimpleHistogramModule::flags() const
 632 {
 633     return efAllowMulticolumn | efAllowMultipoint | efAllowMissing
 634            | efAllowMultipleDataSets;
 635 }
 636
 637
 638 bool
 639 AnalysisDataSimpleHistogramModule::parallelDataStarted(
 640         AbstractAnalysisData              *data,
 641         const AnalysisDataParallelOptions &options)
 642 {
 643     addModule(impl_->averager_);
 644     setDataSetCount(data->dataSetCount());
 645     for (int i = 0; i < data->dataSetCount(); ++i)
 646     {
 647         setColumnCount(i, settings().binCount());
 648     }
 649     impl_->storage_.startParallelDataStorage(this, &moduleManager(), options);
 650     return true;
 651 }
 652
 653
 654 void
 655 AnalysisDataSimpleHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
 656 {
 657     AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
 658     impl_->initFrame(dataSetCount(), &frame);
 659 }
 660
 661
 662 void
 663 AnalysisDataSimpleHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
 664 {
 665     AnalysisDataStorageFrame &frame =
 666         impl_->storage_.currentFrame(points.frameIndex());
 667     frame.selectDataSet(points.dataSetIndex());
 668     for (int i = 0; i < points.columnCount(); ++i)
 669     {
 670         if (points.present(i))
 671         {
 672             const int bin = settings().findBin(points.y(i));
 673             if (bin != -1)
 674             {
 675                 frame.value(bin) += 1;
 676             }
 677         }
 678     }
 679 }
 680
 681
 682 void
 683 AnalysisDataSimpleHistogramModule::frameFinished(const AnalysisDataFrameHeader &header)
 684 {
 685     impl_->storage_.finishFrame(header.index());
 686 }
 687
 688
 689 void
 690 AnalysisDataSimpleHistogramModule::dataFinished()
 691 {
 692     impl_->storage_.finishDataStorage();
 693 }
 694
 695
 696 AnalysisDataFrameRef
 697 AnalysisDataSimpleHistogramModule::tryGetDataFrameInternal(int index) const
 698 {
 699     return impl_->storage_.tryGetDataFrame(index);
 700 }
 701
 702
 703 bool
 704 AnalysisDataSimpleHistogramModule::requestStorageInternal(int nframes)
 705 {
 706     return impl_->storage_.requestStorage(nframes);
 707 }
 708
 709
 710 /********************************************************************
 711  * AnalysisDataWeightedHistogramModule
 712  */
 713
 714 AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule()
 715     : impl_(new internal::BasicHistogramImpl())
 716 {
 717 }
 718
 719
 720 AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule(
 721         const AnalysisHistogramSettings &settings)
 722     : impl_(new internal::BasicHistogramImpl(settings))
 723 {
 724 }
 725
 726
 727 AnalysisDataWeightedHistogramModule::~AnalysisDataWeightedHistogramModule()
 728 {
 729 }
 730
 731
 732 void AnalysisDataWeightedHistogramModule::init(const AnalysisHistogramSettings &settings)
 733 {
 734     impl_->init(settings);
 735 }
 736
 737
 738 AbstractAverageHistogram &
 739 AnalysisDataWeightedHistogramModule::averager()
 740 {
 741     return *impl_->averager_;
 742 }
 743
 744
 745 const AnalysisHistogramSettings &
 746 AnalysisDataWeightedHistogramModule::settings() const
 747 {
 748     return impl_->settings_;
 749 }
 750
 751
 752 int
 753 AnalysisDataWeightedHistogramModule::frameCount() const
 754 {
 755     return impl_->storage_.frameCount();
 756 }
 757
 758
 759 int
 760 AnalysisDataWeightedHistogramModule::flags() const
 761 {
 762     return efAllowMulticolumn | efAllowMultipoint | efAllowMultipleDataSets;
 763 }
 764
 765
 766 bool
 767 AnalysisDataWeightedHistogramModule::parallelDataStarted(
 768         AbstractAnalysisData              *data,
 769         const AnalysisDataParallelOptions &options)
 770 {
 771     addModule(impl_->averager_);
 772     setDataSetCount(data->dataSetCount());
 773     for (int i = 0; i < data->dataSetCount(); ++i)
 774     {
 775         setColumnCount(i, settings().binCount());
 776     }
 777     impl_->storage_.startParallelDataStorage(this, &moduleManager(), options);
 778     return true;
 779 }
 780
 781
 782 void
 783 AnalysisDataWeightedHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
 784 {
 785     AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
 786     impl_->initFrame(dataSetCount(), &frame);
 787 }
 788
 789
 790 void
 791 AnalysisDataWeightedHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
 792 {
 793     if (points.firstColumn() != 0 || points.columnCount() < 2)
 794     {
 795         GMX_THROW(APIError("Invalid data layout"));
 796     }
 797     int bin = settings().findBin(points.y(0));
 798     if (bin != -1)
 799     {
 800         AnalysisDataStorageFrame &frame =
 801             impl_->storage_.currentFrame(points.frameIndex());
 802         frame.selectDataSet(points.dataSetIndex());
 803         for (int i = 1; i < points.columnCount(); ++i)
 804         {
 805             frame.value(bin) += points.y(i);
 806         }
 807     }
 808 }
 809
 810
 811 void
 812 AnalysisDataWeightedHistogramModule::frameFinished(const AnalysisDataFrameHeader &header)
 813 {
 814     impl_->storage_.finishFrame(header.index());
 815 }
 816
 817
 818 void
 819 AnalysisDataWeightedHistogramModule::dataFinished()
 820 {
 821     impl_->storage_.finishDataStorage();
 822 }
 823
 824
 825 AnalysisDataFrameRef
 826 AnalysisDataWeightedHistogramModule::tryGetDataFrameInternal(int index) const
 827 {
 828     return impl_->storage_.tryGetDataFrame(index);
 829 }
 830
 831
 832 bool
 833 AnalysisDataWeightedHistogramModule::requestStorageInternal(int nframes)
 834 {
 835     return impl_->storage_.requestStorage(nframes);
 836 }
 837
 838
 839 /********************************************************************
 840  * AnalysisDataBinAverageModule
 841  */
 842
 843 class AnalysisDataBinAverageModule::Impl
 844 {
 845     public:
 846         Impl() {}
 847         explicit Impl(const AnalysisHistogramSettings &settings)
 848             : settings_(settings)
 849         {
 850         }
 851
 852         //! Histogram settings.
 853         AnalysisHistogramSettings               settings_;
 854         //! Averaging helper objects for each input data set.
 855         std::vector<AnalysisDataFrameAverager>  averagers_;
 856 };
 857
 858 AnalysisDataBinAverageModule::AnalysisDataBinAverageModule()
 859     : impl_(new Impl())
 860 {
 861     setColumnCount(3);
 862 }
 863
 864
 865 AnalysisDataBinAverageModule::AnalysisDataBinAverageModule(
 866         const AnalysisHistogramSettings &settings)
 867     : impl_(new Impl(settings))
 868 {
 869     setRowCount(settings.binCount());
 870     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
 871              settings.binWidth());
 872 }
 873
 874
 875 AnalysisDataBinAverageModule::~AnalysisDataBinAverageModule()
 876 {
 877 }
 878
 879
 880 void
 881 AnalysisDataBinAverageModule::init(const AnalysisHistogramSettings &settings)
 882 {
 883     impl_->settings_ = settings;
 884     setRowCount(settings.binCount());
 885     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
 886              settings.binWidth());
 887 }
 888
 889
 890 const AnalysisHistogramSettings &
 891 AnalysisDataBinAverageModule::settings() const
 892 {
 893     return impl_->settings_;
 894 }
 895
 896
 897 int
 898 AnalysisDataBinAverageModule::flags() const
 899 {
 900     return efAllowMulticolumn | efAllowMultipoint | efAllowMultipleDataSets;
 901 }
 902
 903
 904 void
 905 AnalysisDataBinAverageModule::dataStarted(AbstractAnalysisData *data)
 906 {
 907     setColumnCount(data->dataSetCount());
 908     impl_->averagers_.resize(data->dataSetCount());
 909     for (int i = 0; i < data->dataSetCount(); ++i)
 910     {
 911         impl_->averagers_[i].setColumnCount(rowCount());
 912     }
 913 }
 914
 915
 916 void
 917 AnalysisDataBinAverageModule::frameStarted(const AnalysisDataFrameHeader & /*header*/)
 918 {
 919 }
 920
 921
 922 void
 923 AnalysisDataBinAverageModule::pointsAdded(const AnalysisDataPointSetRef &points)
 924 {
 925     if (points.firstColumn() != 0 || points.columnCount() < 2)
 926     {
 927         GMX_THROW(APIError("Invalid data layout"));
 928     }
 929     int bin = settings().findBin(points.y(0));
 930     if (bin != -1)
 931     {
 932         AnalysisDataFrameAverager &averager = impl_->averagers_[points.dataSetIndex()];
 933         for (int i = 1; i < points.columnCount(); ++i)
 934         {
 935             averager.addValue(bin, points.y(i));
 936         }
 937     }
 938 }
 939
 940
 941 void
 942 AnalysisDataBinAverageModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
 943 {
 944 }
 945
 946
 947 void
 948 AnalysisDataBinAverageModule::dataFinished()
 949 {
 950     allocateValues();
 951     for (int i = 0; i < columnCount(); ++i)
 952     {
 953         AnalysisDataFrameAverager &averager = impl_->averagers_[i];
 954         averager.finish();
 955         for (int j = 0; j < rowCount(); ++j)
 956         {
 957             value(j, i).setValue(averager.average(j),
 958                                  std::sqrt(averager.variance(j)));
 959         }
 960     }
 961     valuesReady();
 962 }
 963
 964 } // namespace gmx