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