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