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(settings().firstEdge(), 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     dest->settings_ = settings_;
 337     return dest;
 338 }
 339
 340
 341 void
 342 AbstractAverageHistogram::normalizeProbability()
 343 {
 344     for (int c = 0; c < columnCount(); ++c)
 345     {
 346         double sum = 0;
 347         for (int i = 0; i < rowCount(); ++i)
 348         {
 349             sum += value(i, c).value();
 350         }
 351         if (sum > 0.0)
 352         {
 353             scaleSingle(c, 1.0 / (sum * xstep()));
 354         }
 355     }
 356 }
 357
 358 void
 359 AbstractAverageHistogram::makeCumulative()
 360 {
 361     for (int c = 0; c < columnCount(); ++c)
 362     {
 363         double sum = 0;
 364         for (int i = 0; i < rowCount(); ++i)
 365         {
 366             sum += value(i, c).value();
 367             // Clear the error, as we don't cumulate that.
 368             value(i, c).clear();
 369             value(i, c).setValue(sum);
 370         }
 371     }
 372     setXAxis(settings().firstEdge() + settings().binWidth(),
 373              settings().binWidth());
 374 }
 375
 376
 377 void
 378 AbstractAverageHistogram::scaleSingle(int index, real factor)
 379 {
 380     for (int i = 0; i < rowCount(); ++i)
 381     {
 382         value(i, index).value() *= factor;
 383         value(i, index).error() *= factor;
 384     }
 385 }
 386
 387
 388 void
 389 AbstractAverageHistogram::scaleAll(real factor)
 390 {
 391     for (int i = 0; i < columnCount(); ++i)
 392     {
 393         scaleSingle(i, factor);
 394     }
 395 }
 396
 397
 398 void
 399 AbstractAverageHistogram::scaleAllByVector(real factor[])
 400 {
 401     for (int c = 0; c < columnCount(); ++c)
 402     {
 403         for (int i = 0; i < rowCount(); ++i)
 404         {
 405             value(i, c).value() *= factor[i];
 406             value(i, c).error() *= factor[i];
 407         }
 408     }
 409 }
 410
 411
 412 /********************************************************************
 413  * BasicAverageHistogramModule
 414  */
 415
 416 namespace internal
 417 {
 418
 419 /*! \internal
 420  * \brief
 421  * Implements average histogram module that averages per-frame histograms.
 422  *
 423  * This class is used for accumulating average histograms in per-frame
 424  * histogram modules (those that use BasicHistogramImpl as their implementation
 425  * class).
 426  * There are two columns, first for the average and second for standard
 427  * deviation.
 428  *
 429  * \ingroup module_analysisdata
 430  */
 431 class BasicAverageHistogramModule : public AbstractAverageHistogram,
 432                                     public AnalysisDataModuleSerial
 433 {
 434     public:
 435         BasicAverageHistogramModule();
 436         //! Creates an average histogram module with defined bin parameters.
 437         explicit BasicAverageHistogramModule(const AnalysisHistogramSettings &settings);
 438
 439         using AbstractAverageHistogram::init;
 440
 441         virtual int flags() const;
 442
 443         virtual void dataStarted(AbstractAnalysisData *data);
 444         virtual void frameStarted(const AnalysisDataFrameHeader &header);
 445         virtual void pointsAdded(const AnalysisDataPointSetRef &points);
 446         virtual void frameFinished(const AnalysisDataFrameHeader &header);
 447         virtual void dataFinished();
 448
 449     private:
 450         //! Averaging helper objects for each input data set.
 451         std::vector<AnalysisDataFrameAverager> averagers_;
 452
 453         // Copy and assign disallowed by base.
 454 };
 455
 456 BasicAverageHistogramModule::BasicAverageHistogramModule()
 457 {
 458 }
 459
 460
 461 BasicAverageHistogramModule::BasicAverageHistogramModule(
 462         const AnalysisHistogramSettings &settings)
 463     : AbstractAverageHistogram(settings)
 464 {
 465 }
 466
 467
 468 int
 469 BasicAverageHistogramModule::flags() const
 470 {
 471     return efAllowMulticolumn | efAllowMultipleDataSets;
 472 }
 473
 474
 475 void
 476 BasicAverageHistogramModule::dataStarted(AbstractAnalysisData *data)
 477 {
 478     setColumnCount(data->dataSetCount());
 479     averagers_.resize(data->dataSetCount());
 480     for (int i = 0; i < data->dataSetCount(); ++i)
 481     {
 482         GMX_RELEASE_ASSERT(rowCount() == data->columnCount(i),
 483                            "Inconsistent data sizes, something is wrong in the initialization");
 484         averagers_[i].setColumnCount(data->columnCount(i));
 485     }
 486 }
 487
 488
 489 void
 490 BasicAverageHistogramModule::frameStarted(const AnalysisDataFrameHeader & /*header*/)
 491 {
 492 }
 493
 494
 495 void
 496 BasicAverageHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
 497 {
 498     averagers_[points.dataSetIndex()].addPoints(points);
 499 }
 500
 501
 502 void
 503 BasicAverageHistogramModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
 504 {
 505 }
 506
 507
 508 void
 509 BasicAverageHistogramModule::dataFinished()
 510 {
 511     allocateValues();
 512     for (int i = 0; i < columnCount(); ++i)
 513     {
 514         averagers_[i].finish();
 515         for (int j = 0; j < rowCount(); ++j)
 516         {
 517             value(j, i).setValue(averagers_[i].average(j),
 518                                  std::sqrt(averagers_[i].variance(j)));
 519         }
 520     }
 521 }
 522
 523
 524 /********************************************************************
 525  * BasicHistogramImpl
 526  */
 527
 528 /*! \internal
 529  * \brief
 530  * Base class for private implementation classes for histogram modules.
 531  *
 532  * Actual implementation classes are derived from this and add an accumulation
 533  * data member that is specific to the histogram type in question.
 534  * This is done like this to keep implementation details out of the header, and
 535  * to not unnecessarily duplicate code.
 536  *
 537  * \ingroup module_analysisdata
 538  */
 539 class BasicHistogramImpl
 540 {
 541     public:
 542         //! Smart pointer to manage an BasicAverageHistogramModule object.
 543         typedef boost::shared_ptr<BasicAverageHistogramModule>
 544             BasicAverageHistogramModulePointer;
 545
 546         BasicHistogramImpl();
 547         //! Creates an histogram impl with defined bin parameters.
 548         explicit BasicHistogramImpl(const AnalysisHistogramSettings &settings);
 549         // Virtual only for simplicity.
 550         virtual ~BasicHistogramImpl();
 551
 552         /*! \brief
 553          * (Re)initializes the histogram from settings.
 554          */
 555         void init(const AnalysisHistogramSettings &settings);
 556
 557         //! Storage implementation object.
 558         AnalysisDataStorage                  storage_;
 559         //! Settings for the histogram object.
 560         AnalysisHistogramSettings            settings_;
 561         //! Averager module.
 562         BasicAverageHistogramModulePointer   averager_;
 563 };
 564
 565 BasicHistogramImpl::BasicHistogramImpl()
 566     : averager_(new BasicAverageHistogramModule())
 567 {
 568 }
 569
 570
 571 BasicHistogramImpl::BasicHistogramImpl(const AnalysisHistogramSettings &settings)
 572     : settings_(settings), averager_(new BasicAverageHistogramModule(settings))
 573 {
 574 }
 575
 576
 577 BasicHistogramImpl::~BasicHistogramImpl()
 578 {
 579 }
 580
 581
 582 void BasicHistogramImpl::init(const AnalysisHistogramSettings &settings)
 583 {
 584     settings_ = settings;
 585     averager_->init(settings);
 586 }
 587
 588 }   // namespace internal
 589
 590
 591 /********************************************************************
 592  * AnalysisDataSimpleHistogramModule
 593  */
 594
 595 /*! \internal \brief
 596  * Private implementation class for AnalysisDataSimpleHistogramModule.
 597  *
 598  * \ingroup module_analysisdata
 599  */
 600 class AnalysisDataSimpleHistogramModule::Impl : public internal::BasicHistogramImpl
 601 {
 602     public:
 603         //! Shorthand for the per-frame accumulation data structure type.
 604         typedef AnalysisDataFrameLocalData<gmx_int64_t> FrameLocalData;
 605
 606         Impl() {}
 607         //! Creates an histogram impl with defined bin parameters.
 608         explicit Impl(const AnalysisHistogramSettings &settings)
 609             : BasicHistogramImpl(settings)
 610         {
 611         }
 612
 613         //! Accumulates the histogram within a frame.
 614         FrameLocalData  accumulator_;
 615 };
 616
 617 AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule()
 618     : impl_(new Impl())
 619 {
 620 }
 621
 622
 623 AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule(
 624         const AnalysisHistogramSettings &settings)
 625     : impl_(new Impl(settings))
 626 {
 627 }
 628
 629
 630 AnalysisDataSimpleHistogramModule::~AnalysisDataSimpleHistogramModule()
 631 {
 632 }
 633
 634
 635 void AnalysisDataSimpleHistogramModule::init(const AnalysisHistogramSettings &settings)
 636 {
 637     impl_->init(settings);
 638 }
 639
 640
 641 AbstractAverageHistogram &
 642 AnalysisDataSimpleHistogramModule::averager()
 643 {
 644     return *impl_->averager_;
 645 }
 646
 647
 648 const AnalysisHistogramSettings &
 649 AnalysisDataSimpleHistogramModule::settings() const
 650 {
 651     return impl_->settings_;
 652 }
 653
 654
 655 int
 656 AnalysisDataSimpleHistogramModule::frameCount() const
 657 {
 658     return impl_->storage_.frameCount();
 659 }
 660
 661
 662 int
 663 AnalysisDataSimpleHistogramModule::flags() const
 664 {
 665     return efAllowMulticolumn | efAllowMultipoint | efAllowMissing
 666            | efAllowMultipleDataSets;
 667 }
 668
 669
 670 bool
 671 AnalysisDataSimpleHistogramModule::parallelDataStarted(
 672         AbstractAnalysisData              *data,
 673         const AnalysisDataParallelOptions &options)
 674 {
 675     addModule(impl_->averager_);
 676     const int dataSetCount = data->dataSetCount();
 677     const int columnCount  = settings().binCount();
 678     setDataSetCount(dataSetCount);
 679     impl_->accumulator_.setDataSetCount(dataSetCount);
 680     for (int i = 0; i < dataSetCount; ++i)
 681     {
 682         setColumnCount(i, columnCount);
 683         impl_->accumulator_.setColumnCount(i, columnCount);
 684     }
 685     impl_->accumulator_.init(options);
 686     impl_->storage_.startParallelDataStorage(this, &moduleManager(), options);
 687     return true;
 688 }
 689
 690
 691 void
 692 AnalysisDataSimpleHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
 693 {
 694     impl_->accumulator_.frameData(header.index()).clear();
 695 }
 696
 697
 698 void
 699 AnalysisDataSimpleHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
 700 {
 701     Impl::FrameLocalData::DataSetHandle handle
 702         = impl_->accumulator_.frameDataSet(points.frameIndex(), points.dataSetIndex());
 703     for (int i = 0; i < points.columnCount(); ++i)
 704     {
 705         if (points.present(i))
 706         {
 707             const int bin = settings().findBin(points.y(i));
 708             if (bin != -1)
 709             {
 710                 handle.value(bin) += 1;
 711             }
 712         }
 713     }
 714 }
 715
 716
 717 void
 718 AnalysisDataSimpleHistogramModule::frameFinished(const AnalysisDataFrameHeader &header)
 719 {
 720     Impl::FrameLocalData::FrameHandle  handle
 721         = impl_->accumulator_.frameData(header.index());
 722     AnalysisDataStorageFrame          &frame = impl_->storage_.startFrame(header);
 723     const int columnCount                    = settings().binCount();
 724     for (int s = 0; s < dataSetCount(); ++s)
 725     {
 726         Impl::FrameLocalData::DataSetHandle dataSet = handle.dataSet(s);
 727         frame.selectDataSet(s);
 728         for (int i = 0; i < columnCount; ++i)
 729         {
 730             frame.setValue(i, dataSet.value(i));
 731         }
 732     }
 733     frame.finishFrame();
 734 }
 735
 736
 737 void
 738 AnalysisDataSimpleHistogramModule::dataFinished()
 739 {
 740     impl_->storage_.finishDataStorage();
 741 }
 742
 743
 744 AnalysisDataFrameRef
 745 AnalysisDataSimpleHistogramModule::tryGetDataFrameInternal(int index) const
 746 {
 747     return impl_->storage_.tryGetDataFrame(index);
 748 }
 749
 750
 751 bool
 752 AnalysisDataSimpleHistogramModule::requestStorageInternal(int nframes)
 753 {
 754     return impl_->storage_.requestStorage(nframes);
 755 }
 756
 757
 758 /********************************************************************
 759  * AnalysisDataWeightedHistogramModule
 760  */
 761
 762 /*! \internal \brief
 763  * Private implementation class for AnalysisDataWeightedHistogramModule.
 764  *
 765  * \ingroup module_analysisdata
 766  */
 767 class AnalysisDataWeightedHistogramModule::Impl : public internal::BasicHistogramImpl
 768 {
 769     public:
 770         //! Shorthand for the per-frame accumulation data structure type.
 771         typedef AnalysisDataFrameLocalData<double> FrameLocalData;
 772
 773         Impl() {}
 774         //! Creates an histogram impl with defined bin parameters.
 775         explicit Impl(const AnalysisHistogramSettings &settings)
 776             : BasicHistogramImpl(settings)
 777         {
 778         }
 779
 780         //! Accumulates the histogram within a frame.
 781         FrameLocalData  accumulator_;
 782 };
 783
 784 AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule()
 785     : impl_(new Impl())
 786 {
 787 }
 788
 789
 790 AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule(
 791         const AnalysisHistogramSettings &settings)
 792     : impl_(new Impl(settings))
 793 {
 794 }
 795
 796
 797 AnalysisDataWeightedHistogramModule::~AnalysisDataWeightedHistogramModule()
 798 {
 799 }
 800
 801
 802 void AnalysisDataWeightedHistogramModule::init(const AnalysisHistogramSettings &settings)
 803 {
 804     impl_->init(settings);
 805 }
 806
 807
 808 AbstractAverageHistogram &
 809 AnalysisDataWeightedHistogramModule::averager()
 810 {
 811     return *impl_->averager_;
 812 }
 813
 814
 815 const AnalysisHistogramSettings &
 816 AnalysisDataWeightedHistogramModule::settings() const
 817 {
 818     return impl_->settings_;
 819 }
 820
 821
 822 int
 823 AnalysisDataWeightedHistogramModule::frameCount() const
 824 {
 825     return impl_->storage_.frameCount();
 826 }
 827
 828
 829 int
 830 AnalysisDataWeightedHistogramModule::flags() const
 831 {
 832     return efAllowMulticolumn | efAllowMultipoint | efAllowMultipleDataSets;
 833 }
 834
 835
 836 bool
 837 AnalysisDataWeightedHistogramModule::parallelDataStarted(
 838         AbstractAnalysisData              *data,
 839         const AnalysisDataParallelOptions &options)
 840 {
 841     addModule(impl_->averager_);
 842     const int dataSetCount = data->dataSetCount();
 843     const int columnCount  = settings().binCount();
 844     setDataSetCount(dataSetCount);
 845     impl_->accumulator_.setDataSetCount(dataSetCount);
 846     for (int i = 0; i < dataSetCount; ++i)
 847     {
 848         setColumnCount(i, columnCount);
 849         impl_->accumulator_.setColumnCount(i, columnCount);
 850     }
 851     impl_->accumulator_.init(options);
 852     impl_->storage_.startParallelDataStorage(this, &moduleManager(), options);
 853     return true;
 854 }
 855
 856
 857 void
 858 AnalysisDataWeightedHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
 859 {
 860     impl_->accumulator_.frameData(header.index()).clear();
 861 }
 862
 863
 864 void
 865 AnalysisDataWeightedHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
 866 {
 867     if (points.firstColumn() != 0 || points.columnCount() < 2)
 868     {
 869         GMX_THROW(APIError("Invalid data layout"));
 870     }
 871     int bin = settings().findBin(points.y(0));
 872     if (bin != -1)
 873     {
 874         Impl::FrameLocalData::DataSetHandle  handle
 875             = impl_->accumulator_.frameDataSet(points.frameIndex(), points.dataSetIndex());
 876         for (int i = 1; i < points.columnCount(); ++i)
 877         {
 878             handle.value(bin) += points.y(i);
 879         }
 880     }
 881 }
 882
 883
 884 void
 885 AnalysisDataWeightedHistogramModule::frameFinished(const AnalysisDataFrameHeader &header)
 886 {
 887     Impl::FrameLocalData::FrameHandle  handle
 888         = impl_->accumulator_.frameData(header.index());
 889     AnalysisDataStorageFrame          &frame = impl_->storage_.startFrame(header);
 890     const int columnCount                    = settings().binCount();
 891     for (int s = 0; s < dataSetCount(); ++s)
 892     {
 893         Impl::FrameLocalData::DataSetHandle dataSet = handle.dataSet(s);
 894         frame.selectDataSet(s);
 895         for (int i = 0; i < columnCount; ++i)
 896         {
 897             frame.setValue(i, dataSet.value(i));
 898         }
 899     }
 900     frame.finishFrame();
 901 }
 902
 903
 904 void
 905 AnalysisDataWeightedHistogramModule::dataFinished()
 906 {
 907     impl_->storage_.finishDataStorage();
 908 }
 909
 910
 911 AnalysisDataFrameRef
 912 AnalysisDataWeightedHistogramModule::tryGetDataFrameInternal(int index) const
 913 {
 914     return impl_->storage_.tryGetDataFrame(index);
 915 }
 916
 917
 918 bool
 919 AnalysisDataWeightedHistogramModule::requestStorageInternal(int nframes)
 920 {
 921     return impl_->storage_.requestStorage(nframes);
 922 }
 923
 924
 925 /********************************************************************
 926  * AnalysisDataBinAverageModule
 927  */
 928
 929 class AnalysisDataBinAverageModule::Impl
 930 {
 931     public:
 932         Impl() {}
 933         explicit Impl(const AnalysisHistogramSettings &settings)
 934             : settings_(settings)
 935         {
 936         }
 937
 938         //! Histogram settings.
 939         AnalysisHistogramSettings               settings_;
 940         //! Averaging helper objects for each input data set.
 941         std::vector<AnalysisDataFrameAverager>  averagers_;
 942 };
 943
 944 AnalysisDataBinAverageModule::AnalysisDataBinAverageModule()
 945     : impl_(new Impl())
 946 {
 947     setColumnCount(3);
 948 }
 949
 950
 951 AnalysisDataBinAverageModule::AnalysisDataBinAverageModule(
 952         const AnalysisHistogramSettings &settings)
 953     : impl_(new Impl(settings))
 954 {
 955     setRowCount(settings.binCount());
 956     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
 957              settings.binWidth());
 958 }
 959
 960
 961 AnalysisDataBinAverageModule::~AnalysisDataBinAverageModule()
 962 {
 963 }
 964
 965
 966 void
 967 AnalysisDataBinAverageModule::init(const AnalysisHistogramSettings &settings)
 968 {
 969     impl_->settings_ = settings;
 970     setRowCount(settings.binCount());
 971     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
 972              settings.binWidth());
 973 }
 974
 975
 976 const AnalysisHistogramSettings &
 977 AnalysisDataBinAverageModule::settings() const
 978 {
 979     return impl_->settings_;
 980 }
 981
 982
 983 int
 984 AnalysisDataBinAverageModule::flags() const
 985 {
 986     return efAllowMulticolumn | efAllowMultipoint | efAllowMultipleDataSets;
 987 }
 988
 989
 990 void
 991 AnalysisDataBinAverageModule::dataStarted(AbstractAnalysisData *data)
 992 {
 993     setColumnCount(data->dataSetCount());
 994     impl_->averagers_.resize(data->dataSetCount());
 995     for (int i = 0; i < data->dataSetCount(); ++i)
 996     {
 997         impl_->averagers_[i].setColumnCount(rowCount());
 998     }
 999 }
1000
1001
1002 void
1003 AnalysisDataBinAverageModule::frameStarted(const AnalysisDataFrameHeader & /*header*/)
1004 {
1005 }
1006
1007
1008 void
1009 AnalysisDataBinAverageModule::pointsAdded(const AnalysisDataPointSetRef &points)
1010 {
1011     if (points.firstColumn() != 0 || points.columnCount() < 2)
1012     {
1013         GMX_THROW(APIError("Invalid data layout"));
1014     }
1015     int bin = settings().findBin(points.y(0));
1016     if (bin != -1)
1017     {
1018         AnalysisDataFrameAverager &averager = impl_->averagers_[points.dataSetIndex()];
1019         for (int i = 1; i < points.columnCount(); ++i)
1020         {
1021             averager.addValue(bin, points.y(i));
1022         }
1023     }
1024 }
1025
1026
1027 void
1028 AnalysisDataBinAverageModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
1029 {
1030 }
1031
1032
1033 void
1034 AnalysisDataBinAverageModule::dataFinished()
1035 {
1036     allocateValues();
1037     for (int i = 0; i < columnCount(); ++i)
1038     {
1039         AnalysisDataFrameAverager &averager = impl_->averagers_[i];
1040         averager.finish();
1041         for (int j = 0; j < rowCount(); ++j)
1042         {
1043             value(j, i).setValue(averager.average(j),
1044                                  std::sqrt(averager.variance(j)));
1045         }
1046     }
1047     valuesReady();
1048 }
1049
1050 } // namespace gmx