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