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,2014,2015,2018,2019, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source 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 "gmxpre.h"
  43
  44 #include "histogram.h"
  45
  46 #include <cmath>
  47
  48 #include <limits>
  49 #include <vector>
  50
  51 #include "gromacs/analysisdata/dataframe.h"
  52 #include "gromacs/analysisdata/datastorage.h"
  53 #include "gromacs/analysisdata/framelocaldata.h"
  54 #include "gromacs/math/functions.h"
  55 #include "gromacs/utility/basedefinitions.h"
  56 #include "gromacs/utility/exceptions.h"
  57 #include "gromacs/utility/gmxassert.h"
  58
  59 #include "frameaverager.h"
  60
  61 namespace
  62 {
  63
  64 //! Value used to signify that a real-valued histogram setting is not set.
  65 const real UNDEFINED = std::numeric_limits<real>::max();
  66 //! Checks whether \p value is defined.
  67 bool isDefined(real value)
  68 {
  69     return value != UNDEFINED;
  70 }
  71
  72 } // namespace
  73
  74 namespace gmx
  75 {
  76
  77 /********************************************************************
  78  * AnalysisHistogramSettingsInitializer
  79  */
  80
  81 AnalysisHistogramSettingsInitializer::AnalysisHistogramSettingsInitializer() :
  82     min_(UNDEFINED),
  83     max_(UNDEFINED),
  84     binWidth_(UNDEFINED),
  85     binCount_(0),
  86     bIntegerBins_(false),
  87     bRoundRange_(false),
  88     bIncludeAll_(false)
  89 {
  90 }
  91
  92
  93 /********************************************************************
  94  * AnalysisHistogramSettings
  95  */
  96
  97 AnalysisHistogramSettings::AnalysisHistogramSettings() :
  98     firstEdge_(0.0),
  99     lastEdge_(0.0),
 100     binWidth_(0.0),
 101     inverseBinWidth_(0.0),
 102     binCount_(0),
 103     bAll_(false)
 104 {
 105 }
 106
 107
 108 AnalysisHistogramSettings::AnalysisHistogramSettings(const AnalysisHistogramSettingsInitializer& settings)
 109 {
 110     GMX_RELEASE_ASSERT(isDefined(settings.min_), "Histogram start value must be defined");
 111     GMX_RELEASE_ASSERT(!isDefined(settings.max_) || settings.max_ > settings.min_,
 112                        "Histogram end value must be larger than start value");
 113     GMX_RELEASE_ASSERT(!isDefined(settings.binWidth_) || settings.binWidth_ > 0.0,
 114                        "Histogram bin width must be positive");
 115     GMX_RELEASE_ASSERT(settings.binCount_ >= 0, "Histogram bin count must be positive");
 116
 117     if (!isDefined(settings.max_))
 118     {
 119         GMX_RELEASE_ASSERT(isDefined(settings.binWidth_) && settings.binCount_ > 0,
 120                            "Not all required values provided");
 121         GMX_RELEASE_ASSERT(!settings.bRoundRange_, "Rounding only supported for min/max ranges");
 122
 123         firstEdge_ = settings.min_;
 124         binCount_  = settings.binCount_;
 125         binWidth_  = settings.binWidth_;
 126         if (settings.bIntegerBins_)
 127         {
 128             firstEdge_ -= 0.5 * binWidth_;
 129         }
 130         lastEdge_ = firstEdge_ + binCount_ * binWidth_;
 131     }
 132     else
 133     {
 134         GMX_RELEASE_ASSERT(!(isDefined(settings.binWidth_) && settings.binCount_ > 0),
 135                            "Conflicting histogram bin specifications");
 136         GMX_RELEASE_ASSERT(isDefined(settings.binWidth_) || settings.binCount_ > 0,
 137                            "Not all required values provided");
 138
 139         if (settings.bRoundRange_)
 140         {
 141             GMX_RELEASE_ASSERT(!settings.bIntegerBins_,
 142                                "Rounding and integer bins cannot be combined");
 143             GMX_RELEASE_ASSERT(isDefined(settings.binWidth_),
 144                                "Rounding only makes sense with defined binwidth");
 145             binWidth_  = settings.binWidth_;
 146             firstEdge_ = binWidth_ * std::floor(settings.min_ / binWidth_);
 147             lastEdge_  = binWidth_ * std::ceil(settings.max_ / binWidth_);
 148             binCount_  = gmx::roundToInt((lastEdge_ - firstEdge_) / binWidth_);
 149         }
 150         else
 151         {
 152             firstEdge_ = settings.min_;
 153             lastEdge_  = settings.max_;
 154             if (settings.binCount_ > 0)
 155             {
 156                 binCount_ = settings.binCount_;
 157                 if (settings.bIntegerBins_)
 158                 {
 159                     GMX_RELEASE_ASSERT(settings.binCount_ > 1,
 160                                        "Bin count must be at least two with integer bins");
 161                     binWidth_ = (lastEdge_ - firstEdge_) / (binCount_ - 1);
 162                     firstEdge_ -= 0.5 * binWidth_;
 163                     lastEdge_ += 0.5 * binWidth_;
 164                 }
 165                 else
 166                 {
 167                     binWidth_ = (lastEdge_ - firstEdge_) / binCount_;
 168                 }
 169             }
 170             else
 171             {
 172                 binWidth_ = settings.binWidth_;
 173                 binCount_ = gmx::roundToInt((lastEdge_ - firstEdge_) / binWidth_);
 174                 if (settings.bIntegerBins_)
 175                 {
 176                     firstEdge_ -= 0.5 * binWidth_;
 177                     ++binCount_;
 178                 }
 179                 lastEdge_ = firstEdge_ + binCount_ * binWidth_;
 180             }
 181         }
 182     }
 183
 184     inverseBinWidth_ = 1.0 / binWidth_;
 185     bAll_            = settings.bIncludeAll_;
 186 }
 187
 188
 189 int AnalysisHistogramSettings::findBin(real y) const
 190 {
 191     if (y < firstEdge_)
 192     {
 193         return bAll_ ? 0 : -1;
 194     }
 195     int bin = static_cast<int>((y - firstEdge_) * inverseBinWidth_);
 196     if (bin >= binCount_)
 197     {
 198         return bAll_ ? binCount_ - 1 : -1;
 199     }
 200     return bin;
 201 }
 202
 203
 204 /********************************************************************
 205  * StaticAverageHistogram
 206  */
 207
 208 namespace
 209 {
 210
 211 /*! \brief
 212  * Represents copies of average histograms.
 213  *
 214  * Methods in AbstractAverageHistogram that return new histogram instances
 215  * return objects of this class.
 216  * Initialization of values is handled in those methods.
 217  *
 218  * \ingroup module_analysisdata
 219  */
 220 class StaticAverageHistogram : public AbstractAverageHistogram
 221 {
 222 public:
 223     StaticAverageHistogram();
 224     //! Creates an average histogram module with defined bin parameters.
 225     explicit StaticAverageHistogram(const AnalysisHistogramSettings& settings);
 226
 227     // Copy and assign disallowed by base.
 228 };
 229
 230 StaticAverageHistogram::StaticAverageHistogram() {}
 231
 232
 233 StaticAverageHistogram::StaticAverageHistogram(const AnalysisHistogramSettings& settings) :
 234     AbstractAverageHistogram(settings)
 235 {
 236 }
 237
 238 } // namespace
 239
 240
 241 /********************************************************************
 242  * AbstractAverageHistogram
 243  */
 244
 245 AbstractAverageHistogram::AbstractAverageHistogram() {}
 246
 247
 248 AbstractAverageHistogram::AbstractAverageHistogram(const AnalysisHistogramSettings& settings) :
 249     settings_(settings)
 250 {
 251     setRowCount(settings.binCount());
 252     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(), settings.binWidth());
 253 }
 254
 255
 256 AbstractAverageHistogram::~AbstractAverageHistogram() {}
 257
 258
 259 void AbstractAverageHistogram::init(const AnalysisHistogramSettings& settings)
 260 {
 261     settings_ = settings;
 262     setRowCount(settings.binCount());
 263     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(), settings.binWidth());
 264 }
 265
 266
 267 AverageHistogramPointer AbstractAverageHistogram::resampleDoubleBinWidth(bool bIntegerBins) const
 268 {
 269     int nbins;
 270     if (bIntegerBins)
 271     {
 272         nbins = (rowCount() + 1) / 2;
 273     }
 274     else
 275     {
 276         nbins = rowCount() / 2;
 277     }
 278
 279     AverageHistogramPointer dest(new StaticAverageHistogram(
 280             histogramFromBins(settings().firstEdge(), nbins, 2 * xstep()).integerBins(bIntegerBins)));
 281     dest->setColumnCount(columnCount());
 282     dest->allocateValues();
 283
 284     int i, j;
 285     for (i = j = 0; i < nbins; ++i)
 286     {
 287         const bool bFirstHalfBin = (bIntegerBins && i == 0);
 288         for (int c = 0; c < columnCount(); ++c)
 289         {
 290             real v1, v2;
 291             real e1, e2;
 292             if (bFirstHalfBin)
 293             {
 294                 v1 = value(0, c).value();
 295                 e1 = value(0, c).error();
 296                 v2 = 0;
 297                 e2 = 0;
 298             }
 299             else
 300             {
 301                 v1 = value(j, c).value();
 302                 e1 = value(j, c).error();
 303                 v2 = value(j + 1, c).value();
 304                 e2 = value(j + 1, c).error();
 305             }
 306             dest->value(i, c).setValue(v1 + v2, std::sqrt(e1 * e1 + e2 * e2));
 307         }
 308         if (bFirstHalfBin)
 309         {
 310             ++j;
 311         }
 312         else
 313         {
 314             j += 2;
 315         }
 316     }
 317     return dest;
 318 }
 319
 320
 321 AverageHistogramPointer AbstractAverageHistogram::clone() const
 322 {
 323     AverageHistogramPointer dest(new StaticAverageHistogram());
 324     copyContents(this, dest.get());
 325     dest->settings_ = settings_;
 326     return dest;
 327 }
 328
 329
 330 void AbstractAverageHistogram::normalizeProbability()
 331 {
 332     for (int c = 0; c < columnCount(); ++c)
 333     {
 334         double sum = 0;
 335         for (int i = 0; i < rowCount(); ++i)
 336         {
 337             sum += value(i, c).value();
 338         }
 339         if (sum > 0.0)
 340         {
 341             scaleSingle(c, 1.0 / (sum * xstep()));
 342         }
 343     }
 344 }
 345
 346 void AbstractAverageHistogram::makeCumulative()
 347 {
 348     for (int c = 0; c < columnCount(); ++c)
 349     {
 350         double sum = 0;
 351         for (int i = 0; i < rowCount(); ++i)
 352         {
 353             sum += value(i, c).value();
 354             // Clear the error, as we don't cumulate that.
 355             value(i, c).clear();
 356             value(i, c).setValue(sum);
 357         }
 358     }
 359     setXAxis(settings().firstEdge() + settings().binWidth(), settings().binWidth());
 360 }
 361
 362
 363 void AbstractAverageHistogram::scaleSingle(int index, real factor)
 364 {
 365     for (int i = 0; i < rowCount(); ++i)
 366     {
 367         value(i, index).value() *= factor;
 368         value(i, index).error() *= factor;
 369     }
 370 }
 371
 372
 373 void AbstractAverageHistogram::scaleAll(real factor)
 374 {
 375     for (int i = 0; i < columnCount(); ++i)
 376     {
 377         scaleSingle(i, factor);
 378     }
 379 }
 380
 381
 382 void AbstractAverageHistogram::scaleAllByVector(const real factor[])
 383 {
 384     for (int c = 0; c < columnCount(); ++c)
 385     {
 386         for (int i = 0; i < rowCount(); ++i)
 387         {
 388             value(i, c).value() *= factor[i];
 389             value(i, c).error() *= factor[i];
 390         }
 391     }
 392 }
 393
 394
 395 /********************************************************************
 396  * BasicAverageHistogramModule
 397  */
 398
 399 namespace internal
 400 {
 401
 402 /*! \internal
 403  * \brief
 404  * Implements average histogram module that averages per-frame histograms.
 405  *
 406  * This class is used for accumulating average histograms in per-frame
 407  * histogram modules (those that use BasicHistogramImpl as their implementation
 408  * class).
 409  * There are two columns, first for the average and second for standard
 410  * deviation.
 411  *
 412  * \ingroup module_analysisdata
 413  */
 414 class BasicAverageHistogramModule : public AbstractAverageHistogram, public AnalysisDataModuleSerial
 415 {
 416 public:
 417     BasicAverageHistogramModule();
 418     //! Creates an average histogram module with defined bin parameters.
 419     explicit BasicAverageHistogramModule(const AnalysisHistogramSettings& settings);
 420
 421     using AbstractAverageHistogram::init;
 422
 423     int flags() const override;
 424
 425     void dataStarted(AbstractAnalysisData* data) override;
 426     void frameStarted(const AnalysisDataFrameHeader& header) override;
 427     void pointsAdded(const AnalysisDataPointSetRef& points) override;
 428     void frameFinished(const AnalysisDataFrameHeader& header) override;
 429     void dataFinished() override;
 430
 431 private:
 432     //! Averaging helper objects for each input data set.
 433     std::vector<AnalysisDataFrameAverager> averagers_;
 434
 435     // Copy and assign disallowed by base.
 436 };
 437
 438 BasicAverageHistogramModule::BasicAverageHistogramModule() {}
 439
 440
 441 BasicAverageHistogramModule::BasicAverageHistogramModule(const AnalysisHistogramSettings& settings) :
 442     AbstractAverageHistogram(settings)
 443 {
 444 }
 445
 446
 447 int BasicAverageHistogramModule::flags() const
 448 {
 449     return efAllowMulticolumn | efAllowMultipleDataSets;
 450 }
 451
 452
 453 void 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 BasicAverageHistogramModule::frameStarted(const AnalysisDataFrameHeader& /*header*/) {}
 467
 468
 469 void BasicAverageHistogramModule::pointsAdded(const AnalysisDataPointSetRef& points)
 470 {
 471     averagers_[points.dataSetIndex()].addPoints(points);
 472 }
 473
 474
 475 void BasicAverageHistogramModule::frameFinished(const AnalysisDataFrameHeader& /*header*/) {}
 476
 477
 478 void BasicAverageHistogramModule::dataFinished()
 479 {
 480     allocateValues();
 481     for (int i = 0; i < columnCount(); ++i)
 482     {
 483         averagers_[i].finish();
 484         for (int j = 0; j < rowCount(); ++j)
 485         {
 486             value(j, i).setValue(averagers_[i].average(j), std::sqrt(averagers_[i].variance(j)));
 487         }
 488     }
 489 }
 490
 491
 492 /********************************************************************
 493  * BasicHistogramImpl
 494  */
 495
 496 /*! \internal
 497  * \brief
 498  * Base class for private implementation classes for histogram modules.
 499  *
 500  * Actual implementation classes are derived from this and add an accumulation
 501  * data member that is specific to the histogram type in question.
 502  * This is done like this to keep implementation details out of the header, and
 503  * to not unnecessarily duplicate code.
 504  *
 505  * \ingroup module_analysisdata
 506  */
 507 class BasicHistogramImpl
 508 {
 509 public:
 510     //! Smart pointer to manage an BasicAverageHistogramModule object.
 511     typedef std::shared_ptr<BasicAverageHistogramModule> BasicAverageHistogramModulePointer;
 512
 513     BasicHistogramImpl();
 514     //! Creates an histogram impl with defined bin parameters.
 515     explicit BasicHistogramImpl(const AnalysisHistogramSettings& settings);
 516     // Virtual only for simplicity.
 517     virtual ~BasicHistogramImpl();
 518
 519     /*! \brief
 520      * (Re)initializes the histogram from settings.
 521      */
 522     void init(const AnalysisHistogramSettings& settings);
 523
 524     //! Storage implementation object.
 525     AnalysisDataStorage storage_;
 526     //! Settings for the histogram object.
 527     AnalysisHistogramSettings settings_;
 528     //! Averager module.
 529     BasicAverageHistogramModulePointer averager_;
 530 };
 531
 532 BasicHistogramImpl::BasicHistogramImpl() : averager_(new BasicAverageHistogramModule()) {}
 533
 534
 535 BasicHistogramImpl::BasicHistogramImpl(const AnalysisHistogramSettings& settings) :
 536     settings_(settings),
 537     averager_(new BasicAverageHistogramModule(settings))
 538 {
 539 }
 540
 541
 542 BasicHistogramImpl::~BasicHistogramImpl() {}
 543
 544
 545 void BasicHistogramImpl::init(const AnalysisHistogramSettings& settings)
 546 {
 547     settings_ = settings;
 548     averager_->init(settings);
 549 }
 550
 551 } // namespace internal
 552
 553
 554 /********************************************************************
 555  * AnalysisDataSimpleHistogramModule
 556  */
 557
 558 /*! \internal \brief
 559  * Private implementation class for AnalysisDataSimpleHistogramModule.
 560  *
 561  * \ingroup module_analysisdata
 562  */
 563 class AnalysisDataSimpleHistogramModule::Impl : public internal::BasicHistogramImpl
 564 {
 565 public:
 566     //! Shorthand for the per-frame accumulation data structure type.
 567     typedef AnalysisDataFrameLocalData<int64_t> FrameLocalData;
 568
 569     Impl() {}
 570     //! Creates an histogram impl with defined bin parameters.
 571     explicit Impl(const AnalysisHistogramSettings& settings) : BasicHistogramImpl(settings) {}
 572
 573     //! Accumulates the histogram within a frame.
 574     FrameLocalData accumulator_;
 575 };
 576
 577 AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule() : impl_(new Impl()) {}
 578
 579
 580 AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule(const AnalysisHistogramSettings& settings) :
 581     impl_(new Impl(settings))
 582 {
 583 }
 584
 585
 586 AnalysisDataSimpleHistogramModule::~AnalysisDataSimpleHistogramModule() {}
 587
 588
 589 void AnalysisDataSimpleHistogramModule::init(const AnalysisHistogramSettings& settings)
 590 {
 591     impl_->init(settings);
 592 }
 593
 594
 595 AbstractAverageHistogram& AnalysisDataSimpleHistogramModule::averager()
 596 {
 597     return *impl_->averager_;
 598 }
 599
 600
 601 const AnalysisHistogramSettings& AnalysisDataSimpleHistogramModule::settings() const
 602 {
 603     return impl_->settings_;
 604 }
 605
 606
 607 int AnalysisDataSimpleHistogramModule::frameCount() const
 608 {
 609     return impl_->storage_.frameCount();
 610 }
 611
 612
 613 int AnalysisDataSimpleHistogramModule::flags() const
 614 {
 615     return efAllowMulticolumn | efAllowMultipoint | efAllowMissing | efAllowMultipleDataSets;
 616 }
 617
 618
 619 bool AnalysisDataSimpleHistogramModule::parallelDataStarted(AbstractAnalysisData*              data,
 620                                                             const AnalysisDataParallelOptions& options)
 621 {
 622     addModule(impl_->averager_);
 623     const int dataSetCount = data->dataSetCount();
 624     const int columnCount  = settings().binCount();
 625     setDataSetCount(dataSetCount);
 626     impl_->accumulator_.setDataSetCount(dataSetCount);
 627     for (int i = 0; i < dataSetCount; ++i)
 628     {
 629         setColumnCount(i, columnCount);
 630         impl_->accumulator_.setColumnCount(i, columnCount);
 631     }
 632     impl_->accumulator_.init(options);
 633     impl_->storage_.startParallelDataStorage(this, &moduleManager(), options);
 634     return true;
 635 }
 636
 637
 638 void AnalysisDataSimpleHistogramModule::frameStarted(const AnalysisDataFrameHeader& header)
 639 {
 640     impl_->accumulator_.frameData(header.index()).clear();
 641 }
 642
 643
 644 void AnalysisDataSimpleHistogramModule::pointsAdded(const AnalysisDataPointSetRef& points)
 645 {
 646     Impl::FrameLocalData::DataSetHandle handle =
 647             impl_->accumulator_.frameDataSet(points.frameIndex(), points.dataSetIndex());
 648     for (int i = 0; i < points.columnCount(); ++i)
 649     {
 650         if (points.present(i))
 651         {
 652             const int bin = settings().findBin(points.y(i));
 653             if (bin != -1)
 654             {
 655                 handle.value(bin) += 1;
 656             }
 657         }
 658     }
 659 }
 660
 661
 662 void AnalysisDataSimpleHistogramModule::frameFinished(const AnalysisDataFrameHeader& header)
 663 {
 664     Impl::FrameLocalData::FrameHandle handle      = impl_->accumulator_.frameData(header.index());
 665     AnalysisDataStorageFrame&         frame       = impl_->storage_.startFrame(header);
 666     const int                         columnCount = settings().binCount();
 667     for (int s = 0; s < dataSetCount(); ++s)
 668     {
 669         Impl::FrameLocalData::DataSetHandle dataSet = handle.dataSet(s);
 670         frame.selectDataSet(s);
 671         for (int i = 0; i < columnCount; ++i)
 672         {
 673             frame.setValue(i, dataSet.value(i));
 674         }
 675     }
 676     frame.finishFrame();
 677 }
 678
 679
 680 void AnalysisDataSimpleHistogramModule::frameFinishedSerial(int frameIndex)
 681 {
 682     impl_->storage_.finishFrameSerial(frameIndex);
 683 }
 684
 685
 686 void AnalysisDataSimpleHistogramModule::dataFinished()
 687 {
 688     impl_->storage_.finishDataStorage();
 689 }
 690
 691
 692 AnalysisDataFrameRef AnalysisDataSimpleHistogramModule::tryGetDataFrameInternal(int index) const
 693 {
 694     return impl_->storage_.tryGetDataFrame(index);
 695 }
 696
 697
 698 bool AnalysisDataSimpleHistogramModule::requestStorageInternal(int nframes)
 699 {
 700     return impl_->storage_.requestStorage(nframes);
 701 }
 702
 703
 704 /********************************************************************
 705  * AnalysisDataWeightedHistogramModule
 706  */
 707
 708 /*! \internal \brief
 709  * Private implementation class for AnalysisDataWeightedHistogramModule.
 710  *
 711  * \ingroup module_analysisdata
 712  */
 713 class AnalysisDataWeightedHistogramModule::Impl : public internal::BasicHistogramImpl
 714 {
 715 public:
 716     //! Shorthand for the per-frame accumulation data structure type.
 717     typedef AnalysisDataFrameLocalData<double> FrameLocalData;
 718
 719     Impl() {}
 720     //! Creates an histogram impl with defined bin parameters.
 721     explicit Impl(const AnalysisHistogramSettings& settings) : BasicHistogramImpl(settings) {}
 722
 723     //! Accumulates the histogram within a frame.
 724     FrameLocalData accumulator_;
 725 };
 726
 727 AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule() : impl_(new Impl()) {}
 728
 729
 730 AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule(const AnalysisHistogramSettings& settings) :
 731     impl_(new Impl(settings))
 732 {
 733 }
 734
 735
 736 AnalysisDataWeightedHistogramModule::~AnalysisDataWeightedHistogramModule() {}
 737
 738
 739 void AnalysisDataWeightedHistogramModule::init(const AnalysisHistogramSettings& settings)
 740 {
 741     impl_->init(settings);
 742 }
 743
 744
 745 AbstractAverageHistogram& AnalysisDataWeightedHistogramModule::averager()
 746 {
 747     return *impl_->averager_;
 748 }
 749
 750
 751 const AnalysisHistogramSettings& AnalysisDataWeightedHistogramModule::settings() const
 752 {
 753     return impl_->settings_;
 754 }
 755
 756
 757 int AnalysisDataWeightedHistogramModule::frameCount() const
 758 {
 759     return impl_->storage_.frameCount();
 760 }
 761
 762
 763 int AnalysisDataWeightedHistogramModule::flags() const
 764 {
 765     return efAllowMulticolumn | efAllowMultipoint | efAllowMultipleDataSets;
 766 }
 767
 768
 769 bool AnalysisDataWeightedHistogramModule::parallelDataStarted(AbstractAnalysisData* data,
 770                                                               const AnalysisDataParallelOptions& options)
 771 {
 772     addModule(impl_->averager_);
 773     const int dataSetCount = data->dataSetCount();
 774     const int columnCount  = settings().binCount();
 775     setDataSetCount(dataSetCount);
 776     impl_->accumulator_.setDataSetCount(dataSetCount);
 777     for (int i = 0; i < dataSetCount; ++i)
 778     {
 779         setColumnCount(i, columnCount);
 780         impl_->accumulator_.setColumnCount(i, columnCount);
 781     }
 782     impl_->accumulator_.init(options);
 783     impl_->storage_.startParallelDataStorage(this, &moduleManager(), options);
 784     return true;
 785 }
 786
 787
 788 void AnalysisDataWeightedHistogramModule::frameStarted(const AnalysisDataFrameHeader& header)
 789 {
 790     impl_->accumulator_.frameData(header.index()).clear();
 791 }
 792
 793
 794 void AnalysisDataWeightedHistogramModule::pointsAdded(const AnalysisDataPointSetRef& points)
 795 {
 796     if (points.firstColumn() != 0 || points.columnCount() < 2)
 797     {
 798         GMX_THROW(APIError("Invalid data layout"));
 799     }
 800     int bin = settings().findBin(points.y(0));
 801     if (bin != -1)
 802     {
 803         Impl::FrameLocalData::DataSetHandle handle =
 804                 impl_->accumulator_.frameDataSet(points.frameIndex(), points.dataSetIndex());
 805         for (int i = 1; i < points.columnCount(); ++i)
 806         {
 807             handle.value(bin) += points.y(i);
 808         }
 809     }
 810 }
 811
 812
 813 void AnalysisDataWeightedHistogramModule::frameFinished(const AnalysisDataFrameHeader& header)
 814 {
 815     Impl::FrameLocalData::FrameHandle handle      = impl_->accumulator_.frameData(header.index());
 816     AnalysisDataStorageFrame&         frame       = impl_->storage_.startFrame(header);
 817     const int                         columnCount = settings().binCount();
 818     for (int s = 0; s < dataSetCount(); ++s)
 819     {
 820         Impl::FrameLocalData::DataSetHandle dataSet = handle.dataSet(s);
 821         frame.selectDataSet(s);
 822         for (int i = 0; i < columnCount; ++i)
 823         {
 824             frame.setValue(i, dataSet.value(i));
 825         }
 826     }
 827     frame.finishFrame();
 828 }
 829
 830
 831 void AnalysisDataWeightedHistogramModule::frameFinishedSerial(int frameIndex)
 832 {
 833     impl_->storage_.finishFrameSerial(frameIndex);
 834 }
 835
 836
 837 void AnalysisDataWeightedHistogramModule::dataFinished()
 838 {
 839     impl_->storage_.finishDataStorage();
 840 }
 841
 842
 843 AnalysisDataFrameRef AnalysisDataWeightedHistogramModule::tryGetDataFrameInternal(int index) const
 844 {
 845     return impl_->storage_.tryGetDataFrame(index);
 846 }
 847
 848
 849 bool AnalysisDataWeightedHistogramModule::requestStorageInternal(int nframes)
 850 {
 851     return impl_->storage_.requestStorage(nframes);
 852 }
 853
 854
 855 /********************************************************************
 856  * AnalysisDataBinAverageModule
 857  */
 858
 859 class AnalysisDataBinAverageModule::Impl
 860 {
 861 public:
 862     Impl() {}
 863     explicit Impl(const AnalysisHistogramSettings& settings) : settings_(settings) {}
 864
 865     //! Histogram settings.
 866     AnalysisHistogramSettings settings_;
 867     //! Averaging helper objects for each input data set.
 868     std::vector<AnalysisDataFrameAverager> averagers_;
 869 };
 870
 871 AnalysisDataBinAverageModule::AnalysisDataBinAverageModule() : impl_(new Impl())
 872 {
 873     setColumnCount(3);
 874 }
 875
 876
 877 AnalysisDataBinAverageModule::AnalysisDataBinAverageModule(const AnalysisHistogramSettings& settings) :
 878     impl_(new Impl(settings))
 879 {
 880     setRowCount(settings.binCount());
 881     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(), settings.binWidth());
 882 }
 883
 884
 885 AnalysisDataBinAverageModule::~AnalysisDataBinAverageModule() {}
 886
 887
 888 void AnalysisDataBinAverageModule::init(const AnalysisHistogramSettings& settings)
 889 {
 890     impl_->settings_ = settings;
 891     setRowCount(settings.binCount());
 892     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(), settings.binWidth());
 893 }
 894
 895
 896 const AnalysisHistogramSettings& AnalysisDataBinAverageModule::settings() const
 897 {
 898     return impl_->settings_;
 899 }
 900
 901
 902 int AnalysisDataBinAverageModule::flags() const
 903 {
 904     return efAllowMulticolumn | efAllowMultipoint | efAllowMultipleDataSets;
 905 }
 906
 907
 908 void AnalysisDataBinAverageModule::dataStarted(AbstractAnalysisData* data)
 909 {
 910     setColumnCount(data->dataSetCount());
 911     impl_->averagers_.resize(data->dataSetCount());
 912     for (int i = 0; i < data->dataSetCount(); ++i)
 913     {
 914         impl_->averagers_[i].setColumnCount(rowCount());
 915     }
 916 }
 917
 918
 919 void AnalysisDataBinAverageModule::frameStarted(const AnalysisDataFrameHeader& /*header*/) {}
 920
 921
 922 void AnalysisDataBinAverageModule::pointsAdded(const AnalysisDataPointSetRef& points)
 923 {
 924     if (points.firstColumn() != 0 || points.columnCount() < 2)
 925     {
 926         GMX_THROW(APIError("Invalid data layout"));
 927     }
 928     int bin = settings().findBin(points.y(0));
 929     if (bin != -1)
 930     {
 931         AnalysisDataFrameAverager& averager = impl_->averagers_[points.dataSetIndex()];
 932         for (int i = 1; i < points.columnCount(); ++i)
 933         {
 934             averager.addValue(bin, points.y(i));
 935         }
 936     }
 937 }
 938
 939
 940 void AnalysisDataBinAverageModule::frameFinished(const AnalysisDataFrameHeader& /*header*/) {}
 941
 942
 943 void AnalysisDataBinAverageModule::dataFinished()
 944 {
 945     allocateValues();
 946     for (int i = 0; i < columnCount(); ++i)
 947     {
 948         AnalysisDataFrameAverager& averager = impl_->averagers_[i];
 949         averager.finish();
 950         for (int j = 0; j < rowCount(); ++j)
 951         {
 952             value(j, i).setValue(averager.average(j), std::sqrt(averager.variance(j)));
 953         }
 954     }
 955     valuesReady();
 956 }
 957
 958 } // namespace gmx