src/gromacs/analysisdata/modules/histogram.cpp

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  10  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  11  * Copyright (c) 2001-2009, The GROMACS development team,
  12  * check out http://www.gromacs.org for more information.
  13
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * If you want to redistribute modifications, please consider that
  20  * scientific software is very special. Version control is crucial -
  21  * bugs must be traceable. We will be happy to consider code for
  22  * inclusion in the official distribution, but derived work must not
  23  * be called official GROMACS. Details are found in the README & COPYING
  24  * files - if they are missing, get the official version at www.gromacs.org.
  25  *
  26  * To help us fund GROMACS development, we humbly ask that you cite
  27  * the papers on the package - you can find them in the top README file.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  */
  31 /*! \internal \file
  32  * \brief
  33  * Implements classes in histogram.h.
  34  *
  35  * \author Teemu Murtola <teemu.murtola@cbr.su.se>
  36  * \ingroup module_analysisdata
  37  */
  38 #include "gromacs/analysisdata/modules/histogram.h"
  39
  40 #include <cmath>
  41
  42 #include <limits>
  43
  44 #include "gromacs/basicmath.h"
  45 #include "gromacs/analysisdata/dataframe.h"
  46 #include "gromacs/analysisdata/datastorage.h"
  47 #include "gromacs/fatalerror/exceptions.h"
  48 #include "gromacs/fatalerror/gmxassert.h"
  49
  50 #include "histogram-impl.h"
  51
  52 static const real UNDEFINED = std::numeric_limits<real>::max();
  53 static bool isDefined(real value)
  54 {
  55     return value != UNDEFINED;
  56 }
  57
  58 namespace gmx
  59 {
  60
  61 /********************************************************************
  62  * AnalysisHistogramSettingsInitializer
  63  */
  64
  65 AnalysisHistogramSettingsInitializer::AnalysisHistogramSettingsInitializer()
  66     : min_(UNDEFINED), max_(UNDEFINED), binWidth_(UNDEFINED),
  67       binCount_(0), bIntegerBins_(false), bRoundRange_(false),
  68       bIncludeAll_(false)
  69 {
  70 }
  71
  72
  73 /********************************************************************
  74  * AnalysisHistogramSettings
  75  */
  76
  77 AnalysisHistogramSettings::AnalysisHistogramSettings()
  78     : firstEdge_(0.0), lastEdge_(0.0), binWidth_(0.0), inverseBinWidth_(0.0),
  79       binCount_(0), bAll_(false)
  80 {
  81 }
  82
  83
  84 AnalysisHistogramSettings::AnalysisHistogramSettings(
  85         const AnalysisHistogramSettingsInitializer &settings)
  86 {
  87     GMX_RELEASE_ASSERT(isDefined(settings.min_),
  88                        "Histogram start value must be defined");
  89     GMX_RELEASE_ASSERT(!isDefined(settings.max_) || settings.max_ > settings.min_,
  90                        "Histogram end value must be larger than start value");
  91     GMX_RELEASE_ASSERT(!isDefined(settings.binWidth_) || settings.binWidth_ > 0.0,
  92                        "Histogram bin width must be positive");
  93     GMX_RELEASE_ASSERT(settings.binCount_ >= 0,
  94                        "Histogram bin count must be positive");
  95
  96     if (!isDefined(settings.max_))
  97     {
  98         GMX_RELEASE_ASSERT(isDefined(settings.binWidth_) && settings.binCount_ > 0,
  99                            "Not all required values provided");
 100         GMX_RELEASE_ASSERT(!settings.bRoundRange_,
 101                            "Rounding only supported for min/max ranges");
 102
 103         firstEdge_ = settings.min_;
 104         binCount_  = settings.binCount_;
 105         binWidth_  = settings.binWidth_;
 106         if (settings.bIntegerBins_)
 107         {
 108             firstEdge_ -= 0.5 * binWidth_;
 109         }
 110         lastEdge_ = firstEdge_ + binCount_ * binWidth_;
 111     }
 112     else
 113     {
 114         GMX_RELEASE_ASSERT(!(isDefined(settings.binWidth_) && settings.binCount_ > 0),
 115                            "Conflicting histogram bin specifications");
 116         GMX_RELEASE_ASSERT(isDefined(settings.binWidth_) || settings.binCount_ > 0,
 117                            "Not all required values provided");
 118
 119         if (settings.bRoundRange_)
 120         {
 121             GMX_RELEASE_ASSERT(!settings.bIntegerBins_,
 122                                "Rounding and integer bins cannot be combined");
 123             GMX_RELEASE_ASSERT(isDefined(settings.binWidth_),
 124                                "Rounding only makes sense with defined binwidth");
 125             binWidth_  = settings.binWidth_;
 126             firstEdge_ = binWidth_ * floor(settings.min_ / binWidth_);
 127             lastEdge_  = binWidth_ * ceil(settings.max_ / binWidth_);
 128             binCount_  = static_cast<int>((lastEdge_ - firstEdge_) / binWidth_ + 0.5);
 129         }
 130         else
 131         {
 132             firstEdge_     = settings.min_;
 133             lastEdge_     = settings.max_;
 134             if (settings.binCount_ > 0)
 135             {
 136                 binCount_ = settings.binCount_;
 137                 if (settings.bIntegerBins_)
 138                 {
 139                     GMX_RELEASE_ASSERT(settings.binCount_ > 1,
 140                                        "Bin count must be at least two with integer bins");
 141                     binWidth_   = (lastEdge_ - firstEdge_) / (binCount_ - 1);
 142                     firstEdge_ -= 0.5 * binWidth_;
 143                     lastEdge_  += 0.5 * binWidth_;
 144                 }
 145                 else
 146                 {
 147                     binWidth_ = (lastEdge_ - firstEdge_) / binCount_;
 148                 }
 149             }
 150             else
 151             {
 152                 binWidth_ = settings.binWidth_;
 153                 binCount_ = static_cast<int>((lastEdge_ - firstEdge_) / binWidth_ + 0.5);
 154                 if (settings.bIntegerBins_)
 155                 {
 156                     firstEdge_ -= 0.5 * binWidth_;
 157                     ++binCount_;
 158                 }
 159                 lastEdge_ = firstEdge_ + binCount_ * binWidth_;
 160             }
 161         }
 162     }
 163
 164     inverseBinWidth_ = 1.0 / binWidth_;
 165     bAll_            = settings.bIncludeAll_;
 166 }
 167
 168
 169 int
 170 AnalysisHistogramSettings::findBin(real y) const
 171 {
 172     if (y < firstEdge_)
 173     {
 174         return bAll_ ? 0 : -1;
 175     }
 176     int bin = static_cast<int>((y - firstEdge_) * inverseBinWidth_);
 177     if (bin >= binCount_)
 178     {
 179         return bAll_ ? binCount_ - 1 : -1;
 180     }
 181     return bin;
 182 }
 183
 184
 185 /********************************************************************
 186  * AbstractAverageHistogram
 187  */
 188
 189 AbstractAverageHistogram::AbstractAverageHistogram()
 190 {
 191 }
 192
 193
 194 AbstractAverageHistogram::AbstractAverageHistogram(
 195         const AnalysisHistogramSettings &settings)
 196     : settings_(settings)
 197 {
 198     setRowCount(settings.binCount());
 199     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
 200              settings.binWidth());
 201 }
 202
 203
 204 AbstractAverageHistogram::~AbstractAverageHistogram()
 205 {
 206 }
 207
 208
 209 void
 210 AbstractAverageHistogram::init(const AnalysisHistogramSettings &settings)
 211 {
 212     settings_ = settings;
 213     setRowCount(settings.binCount());
 214     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
 215              settings.binWidth());
 216 }
 217
 218
 219 AverageHistogramPointer
 220 AbstractAverageHistogram::resampleDoubleBinWidth(bool bIntegerBins) const
 221 {
 222     int nbins;
 223     if (bIntegerBins)
 224     {
 225         nbins = (rowCount() + 1) / 2;
 226     }
 227     else
 228     {
 229         nbins = rowCount() / 2;
 230     }
 231
 232     AverageHistogramPointer dest(
 233         new internal::StaticAverageHistogram(
 234             histogramFromBins(xstart(), nbins, 2*xstep())
 235                 .integerBins(bIntegerBins)));
 236     dest->setColumnCount(columnCount());
 237     dest->allocateValues();
 238
 239     int  i, j;
 240     for (i = j = 0; i < nbins; ++i)
 241     {
 242         const bool bFirstHalfBin = (bIntegerBins && i == 0);
 243         for (int c = 0; c < columnCount(); ++c)
 244         {
 245             real  v1, v2;
 246             if (bFirstHalfBin)
 247             {
 248                 v1 = value(0, c);
 249                 v2 = 0;
 250             }
 251             else
 252             {
 253                 v1 = value(j, c);
 254                 v2 = value(j + 1, c);
 255             }
 256             if (c == 1)
 257             {
 258                 dest->setValue(i, c, sqrt(v1 * v1 + v2 * v2));
 259             }
 260             else
 261             {
 262                 dest->setValue(i, c, v1 + v2);
 263             }
 264         }
 265         if (bFirstHalfBin)
 266         {
 267             ++j;
 268         }
 269         else
 270         {
 271             j += 2;
 272         }
 273     }
 274     return dest;
 275 }
 276
 277
 278 AverageHistogramPointer
 279 AbstractAverageHistogram::clone() const
 280 {
 281     AverageHistogramPointer dest(new internal::StaticAverageHistogram());
 282     copyContents(this, dest.get());
 283     return dest;
 284 }
 285
 286
 287 void
 288 AbstractAverageHistogram::normalizeProbability()
 289 {
 290     real sum = 0;
 291     for (int i = 0; i < rowCount(); ++i)
 292     {
 293         sum += value(i, 0);
 294     }
 295     scale(1.0 / (sum * xstep()));
 296 }
 297
 298
 299 void
 300 AbstractAverageHistogram::scale(real norm)
 301 {
 302     for (int i = 0; i < rowCount(); ++i)
 303     {
 304         value(i, 0) *= norm;
 305         value(i, 1) *= norm;
 306     }
 307 }
 308
 309
 310 void
 311 AbstractAverageHistogram::scaleVector(real norm[])
 312 {
 313     for (int i = 0; i < rowCount(); ++i)
 314     {
 315         value(i, 0) *= norm[i];
 316         value(i, 1) *= norm[i];
 317     }
 318 }
 319
 320
 321 /********************************************************************
 322  * StaticAverageHistogram
 323  */
 324
 325 namespace internal
 326 {
 327
 328 StaticAverageHistogram::StaticAverageHistogram()
 329 {
 330 }
 331
 332
 333 StaticAverageHistogram::StaticAverageHistogram(
 334         const AnalysisHistogramSettings &settings)
 335     : AbstractAverageHistogram(settings)
 336 {
 337 }
 338
 339
 340 /********************************************************************
 341  * BasicAverageHistogramModule
 342  */
 343
 344 BasicAverageHistogramModule::BasicAverageHistogramModule()
 345     : frameCount_(0)
 346 {
 347     setColumnCount(2);
 348 }
 349
 350
 351 BasicAverageHistogramModule::BasicAverageHistogramModule(
 352         const AnalysisHistogramSettings &settings)
 353     : AbstractAverageHistogram(settings), frameCount_(0)
 354 {
 355     setColumnCount(2);
 356 }
 357
 358
 359 int
 360 BasicAverageHistogramModule::flags() const
 361 {
 362     return efAllowMulticolumn;
 363 }
 364
 365
 366 void
 367 BasicAverageHistogramModule::dataStarted(AbstractAnalysisData *data)
 368 {
 369     GMX_RELEASE_ASSERT(rowCount() == data->columnCount(),
 370                        "Inconsistent data sizes, something is wrong in the initialization");
 371     allocateValues();
 372 }
 373
 374
 375 void
 376 BasicAverageHistogramModule::frameStarted(const AnalysisDataFrameHeader & /*header*/)
 377 {
 378 }
 379
 380
 381 void
 382 BasicAverageHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
 383 {
 384     int firstcol = points.firstColumn();
 385     for (int i = 0; i < points.columnCount(); ++i)
 386     {
 387         real y = points.y(i);
 388         value(firstcol + i, 0) += y;
 389         value(firstcol + i, 1) += y * y;
 390     }
 391 }
 392
 393
 394 void
 395 BasicAverageHistogramModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
 396 {
 397     ++frameCount_;
 398 }
 399
 400
 401 void
 402 BasicAverageHistogramModule::dataFinished()
 403 {
 404     for (int i = 0; i < rowCount(); ++i)
 405     {
 406         real ave = value(i, 0) / frameCount_;
 407         real std = sqrt(value(i, 1) / frameCount_ - ave * ave);
 408         setValue(i, 0, ave);
 409         setValue(i, 1, std);
 410     }
 411 }
 412
 413
 414 /********************************************************************
 415  * BasicHistogramImpl
 416  */
 417
 418 BasicHistogramImpl::BasicHistogramImpl()
 419     : averager_(new BasicAverageHistogramModule())
 420 {
 421 }
 422
 423
 424 BasicHistogramImpl::BasicHistogramImpl(const AnalysisHistogramSettings &settings)
 425     : settings_(settings), averager_(new BasicAverageHistogramModule(settings))
 426 {
 427 }
 428
 429
 430 BasicHistogramImpl::~BasicHistogramImpl()
 431 {
 432 }
 433
 434
 435 void BasicHistogramImpl::init(const AnalysisHistogramSettings &settings)
 436 {
 437     settings_ = settings;
 438     averager_->init(settings);
 439 }
 440
 441
 442 void
 443 BasicHistogramImpl::initFrame(AnalysisDataStorageFrame *frame)
 444 {
 445     for (int i = 0; i < frame->columnCount(); ++i)
 446     {
 447         frame->setValue(i, 0.0);
 448     }
 449 }
 450
 451 } // namespace internal
 452
 453
 454 /********************************************************************
 455  * AnalysisDataSimpleHistogramModule
 456  */
 457
 458 AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule()
 459     : impl_(new internal::BasicHistogramImpl())
 460 {
 461 }
 462
 463
 464 AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule(
 465         const AnalysisHistogramSettings &settings)
 466     : impl_(new internal::BasicHistogramImpl(settings))
 467 {
 468 }
 469
 470
 471 AnalysisDataSimpleHistogramModule::~AnalysisDataSimpleHistogramModule()
 472 {
 473 }
 474
 475
 476 void AnalysisDataSimpleHistogramModule::init(const AnalysisHistogramSettings &settings)
 477 {
 478     impl_->init(settings);
 479 }
 480
 481
 482 AbstractAverageHistogram &
 483 AnalysisDataSimpleHistogramModule::averager()
 484 {
 485     return *impl_->averager_;
 486 }
 487
 488
 489 const AnalysisHistogramSettings &
 490 AnalysisDataSimpleHistogramModule::settings() const
 491 {
 492     return impl_->settings_;
 493 }
 494
 495
 496 int
 497 AnalysisDataSimpleHistogramModule::flags() const
 498 {
 499     return efAllowMulticolumn | efAllowMultipoint;
 500 }
 501
 502
 503 void
 504 AnalysisDataSimpleHistogramModule::dataStarted(AbstractAnalysisData *data)
 505 {
 506     addModule(impl_->averager_);
 507     setColumnCount(settings().binCount());
 508     notifyDataStart();
 509     impl_->storage_.startDataStorage(this);
 510 }
 511
 512
 513 void
 514 AnalysisDataSimpleHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
 515 {
 516     AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
 517     impl_->initFrame(&frame);
 518 }
 519
 520
 521 void
 522 AnalysisDataSimpleHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
 523 {
 524     AnalysisDataStorageFrame &frame =
 525         impl_->storage_.currentFrame(points.frameIndex());
 526     for (int i = 0; i < points.columnCount(); ++i)
 527     {
 528         int bin = settings().findBin(points.y(i));
 529         if (bin != -1)
 530         {
 531             frame.value(bin) += 1;
 532         }
 533     }
 534 }
 535
 536
 537 void
 538 AnalysisDataSimpleHistogramModule::frameFinished(const AnalysisDataFrameHeader &header)
 539 {
 540     impl_->storage_.finishFrame(header.index());
 541 }
 542
 543
 544 void
 545 AnalysisDataSimpleHistogramModule::dataFinished()
 546 {
 547     notifyDataFinish();
 548 }
 549
 550
 551 AnalysisDataFrameRef
 552 AnalysisDataSimpleHistogramModule::tryGetDataFrameInternal(int index) const
 553 {
 554     return impl_->storage_.tryGetDataFrame(index);
 555 }
 556
 557
 558 bool
 559 AnalysisDataSimpleHistogramModule::requestStorageInternal(int nframes)
 560 {
 561     return impl_->storage_.requestStorage(nframes);
 562 }
 563
 564
 565 /********************************************************************
 566  * AnalysisDataWeightedHistogramModule
 567  */
 568
 569 AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule()
 570     : impl_(new internal::BasicHistogramImpl())
 571 {
 572 }
 573
 574
 575 AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule(
 576         const AnalysisHistogramSettings &settings)
 577     : impl_(new internal::BasicHistogramImpl(settings))
 578 {
 579 }
 580
 581
 582 AnalysisDataWeightedHistogramModule::~AnalysisDataWeightedHistogramModule()
 583 {
 584 }
 585
 586
 587 void AnalysisDataWeightedHistogramModule::init(const AnalysisHistogramSettings &settings)
 588 {
 589     impl_->init(settings);
 590 }
 591
 592
 593 AbstractAverageHistogram &
 594 AnalysisDataWeightedHistogramModule::averager()
 595 {
 596     return *impl_->averager_;
 597 }
 598
 599
 600 const AnalysisHistogramSettings &
 601 AnalysisDataWeightedHistogramModule::settings() const
 602 {
 603     return impl_->settings_;
 604 }
 605
 606
 607 int
 608 AnalysisDataWeightedHistogramModule::flags() const
 609 {
 610     return efAllowMulticolumn | efAllowMultipoint;
 611 }
 612
 613
 614 void
 615 AnalysisDataWeightedHistogramModule::dataStarted(AbstractAnalysisData *data)
 616 {
 617     addModule(impl_->averager_);
 618     setColumnCount(settings().binCount());
 619     notifyDataStart();
 620     impl_->storage_.startDataStorage(this);
 621 }
 622
 623
 624 void
 625 AnalysisDataWeightedHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
 626 {
 627     AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
 628     impl_->initFrame(&frame);
 629 }
 630
 631
 632 void
 633 AnalysisDataWeightedHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
 634 {
 635     if (points.firstColumn() != 0 || points.columnCount() < 2)
 636     {
 637         GMX_THROW(APIError("Invalid data layout"));
 638     }
 639     int bin = settings().findBin(points.y(0));
 640     if (bin != -1)
 641     {
 642         AnalysisDataStorageFrame &frame =
 643             impl_->storage_.currentFrame(points.frameIndex());
 644         for (int i = 1; i < points.columnCount(); ++i)
 645         {
 646             frame.value(bin) += points.y(i);
 647         }
 648     }
 649 }
 650
 651
 652 void
 653 AnalysisDataWeightedHistogramModule::frameFinished(const AnalysisDataFrameHeader &header)
 654 {
 655     impl_->storage_.finishFrame(header.index());
 656 }
 657
 658
 659 void
 660 AnalysisDataWeightedHistogramModule::dataFinished()
 661 {
 662     notifyDataFinish();
 663 }
 664
 665
 666 AnalysisDataFrameRef
 667 AnalysisDataWeightedHistogramModule::tryGetDataFrameInternal(int index) const
 668 {
 669     return impl_->storage_.tryGetDataFrame(index);
 670 }
 671
 672
 673 bool
 674 AnalysisDataWeightedHistogramModule::requestStorageInternal(int nframes)
 675 {
 676     return impl_->storage_.requestStorage(nframes);
 677 }
 678
 679
 680 /********************************************************************
 681  * AnalysisDataBinAverageModule
 682  */
 683
 684 AnalysisDataBinAverageModule::AnalysisDataBinAverageModule()
 685 {
 686     setColumnCount(3);
 687 }
 688
 689
 690 AnalysisDataBinAverageModule::AnalysisDataBinAverageModule(
 691         const AnalysisHistogramSettings &settings)
 692     : settings_(settings)
 693 {
 694     setColumnCount(3);
 695     setRowCount(settings.binCount());
 696     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
 697              settings.binWidth());
 698 }
 699
 700
 701 AnalysisDataBinAverageModule::~AnalysisDataBinAverageModule()
 702 {
 703 }
 704
 705
 706 void
 707 AnalysisDataBinAverageModule::init(const AnalysisHistogramSettings &settings)
 708 {
 709     settings_ = settings;
 710     setRowCount(settings.binCount());
 711     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
 712              settings.binWidth());
 713 }
 714
 715
 716 int
 717 AnalysisDataBinAverageModule::flags() const
 718 {
 719     return efAllowMulticolumn | efAllowMultipoint;
 720 }
 721
 722
 723 void
 724 AnalysisDataBinAverageModule::dataStarted(AbstractAnalysisData * /*data*/)
 725 {
 726     allocateValues();
 727 }
 728
 729
 730 void
 731 AnalysisDataBinAverageModule::frameStarted(const AnalysisDataFrameHeader & /*header*/)
 732 {
 733 }
 734
 735
 736 void
 737 AnalysisDataBinAverageModule::pointsAdded(const AnalysisDataPointSetRef &points)
 738 {
 739     if (points.firstColumn() != 0 || points.columnCount() < 2)
 740     {
 741         GMX_THROW(APIError("Invalid data layout"));
 742     }
 743     int bin = settings().findBin(points.y(0));
 744     if (bin != -1)
 745     {
 746         for (int i = 1; i < points.columnCount(); ++i)
 747         {
 748             real y = points.y(i);
 749             value(bin, 0) += y;
 750             value(bin, 1) += y * y;
 751         }
 752         value(bin, 2) += points.columnCount() - 1;
 753     }
 754 }
 755
 756
 757 void
 758 AnalysisDataBinAverageModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
 759 {
 760 }
 761
 762
 763 void
 764 AnalysisDataBinAverageModule::dataFinished()
 765 {
 766     for (int i = 0; i < settings().binCount(); ++i)
 767     {
 768         real n = value(i, 2);
 769         if (n > 0)
 770         {
 771             real ave = value(i, 0) / n;
 772             real std = sqrt(value(i, 1) / n - ave * ave);
 773             setValue(i, 0, ave);
 774             setValue(i, 1, std);
 775         }
 776     }
 777     valuesReady();
 778 }
 779
 780 } // namespace gmx