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