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