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