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     real minx = xstart();
 233     AverageHistogramPointer dest(
 234         new internal::StaticAverageHistogram(
 235             histogramFromBins(xstart(), nbins, 2*xstep())
 236                 .integerBins(bIntegerBins)));
 237     dest->setColumnCount(columnCount());
 238     dest->allocateValues();
 239
 240     int  i, j;
 241     for (i = j = 0; i < nbins; ++i)
 242     {
 243         const bool bFirstHalfBin = (bIntegerBins && i == 0);
 244         for (int c = 0; c < columnCount(); ++c)
 245         {
 246             real  v1, v2;
 247             if (bFirstHalfBin)
 248             {
 249                 v1 = value(0, c);
 250                 v2 = 0;
 251             }
 252             else
 253             {
 254                 v1 = value(j, c);
 255                 v2 = value(j + 1, c);
 256             }
 257             if (c == 1)
 258             {
 259                 dest->setValue(i, c, sqrt(v1 * v1 + v2 * v2));
 260             }
 261             else
 262             {
 263                 dest->setValue(i, c, v1 + v2);
 264             }
 265         }
 266         if (bFirstHalfBin)
 267         {
 268             ++j;
 269         }
 270         else
 271         {
 272             j += 2;
 273         }
 274     }
 275     return dest;
 276 }
 277
 278
 279 AverageHistogramPointer
 280 AbstractAverageHistogram::clone() const
 281 {
 282     AverageHistogramPointer dest(new internal::StaticAverageHistogram());
 283     copyContents(this, dest.get());
 284     return dest;
 285 }
 286
 287
 288 void
 289 AbstractAverageHistogram::normalizeProbability()
 290 {
 291     real sum = 0;
 292     for (int i = 0; i < rowCount(); ++i)
 293     {
 294         sum += value(i, 0);
 295     }
 296     scale(1.0 / (sum * xstep()));
 297 }
 298
 299
 300 void
 301 AbstractAverageHistogram::scale(real norm)
 302 {
 303     for (int i = 0; i < rowCount(); ++i)
 304     {
 305         value(i, 0) *= norm;
 306         value(i, 1) *= norm;
 307     }
 308 }
 309
 310
 311 void
 312 AbstractAverageHistogram::scaleVector(real norm[])
 313 {
 314     for (int i = 0; i < rowCount(); ++i)
 315     {
 316         value(i, 0) *= norm[i];
 317         value(i, 1) *= norm[i];
 318     }
 319 }
 320
 321
 322 /********************************************************************
 323  * StaticAverageHistogram
 324  */
 325
 326 namespace internal
 327 {
 328
 329 StaticAverageHistogram::StaticAverageHistogram()
 330 {
 331 }
 332
 333
 334 StaticAverageHistogram::StaticAverageHistogram(
 335         const AnalysisHistogramSettings &settings)
 336     : AbstractAverageHistogram(settings)
 337 {
 338 }
 339
 340
 341 /********************************************************************
 342  * BasicAverageHistogramModule
 343  */
 344
 345 BasicAverageHistogramModule::BasicAverageHistogramModule()
 346     : frameCount_(0)
 347 {
 348     setColumnCount(2);
 349 }
 350
 351
 352 BasicAverageHistogramModule::BasicAverageHistogramModule(
 353         const AnalysisHistogramSettings &settings)
 354     : AbstractAverageHistogram(settings), frameCount_(0)
 355 {
 356     setColumnCount(2);
 357 }
 358
 359
 360 int
 361 BasicAverageHistogramModule::flags() const
 362 {
 363     return efAllowMulticolumn;
 364 }
 365
 366
 367 void
 368 BasicAverageHistogramModule::dataStarted(AbstractAnalysisData *data)
 369 {
 370     GMX_RELEASE_ASSERT(rowCount() == data->columnCount(),
 371                        "Inconsistent data sizes, something is wrong in the initialization");
 372     allocateValues();
 373 }
 374
 375
 376 void
 377 BasicAverageHistogramModule::frameStarted(const AnalysisDataFrameHeader & /*header*/)
 378 {
 379 }
 380
 381
 382 void
 383 BasicAverageHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
 384 {
 385     int firstcol = points.firstColumn();
 386     for (int i = 0; i < points.columnCount(); ++i)
 387     {
 388         real y = points.y(i);
 389         value(firstcol + i, 0) += y;
 390         value(firstcol + i, 1) += y * y;
 391     }
 392 }
 393
 394
 395 void
 396 BasicAverageHistogramModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
 397 {
 398     ++frameCount_;
 399 }
 400
 401
 402 void
 403 BasicAverageHistogramModule::dataFinished()
 404 {
 405     for (int i = 0; i < rowCount(); ++i)
 406     {
 407         real ave = value(i, 0) / frameCount_;
 408         real std = sqrt(value(i, 1) / frameCount_ - ave * ave);
 409         setValue(i, 0, ave);
 410         setValue(i, 1, std);
 411     }
 412 }
 413
 414
 415 /********************************************************************
 416  * BasicHistogramImpl
 417  */
 418
 419 BasicHistogramImpl::BasicHistogramImpl()
 420     : averager_(new BasicAverageHistogramModule())
 421 {
 422 }
 423
 424
 425 BasicHistogramImpl::BasicHistogramImpl(const AnalysisHistogramSettings &settings)
 426     : settings_(settings), averager_(new BasicAverageHistogramModule(settings))
 427 {
 428 }
 429
 430
 431 BasicHistogramImpl::~BasicHistogramImpl()
 432 {
 433 }
 434
 435
 436 void BasicHistogramImpl::init(const AnalysisHistogramSettings &settings)
 437 {
 438     settings_ = settings;
 439     averager_->init(settings);
 440 }
 441
 442
 443 void
 444 BasicHistogramImpl::initFrame(AnalysisDataStorageFrame *frame)
 445 {
 446     for (int i = 0; i < frame->columnCount(); ++i)
 447     {
 448         frame->setValue(i, 0.0);
 449     }
 450 }
 451
 452 } // namespace internal
 453
 454
 455 /********************************************************************
 456  * AnalysisDataSimpleHistogramModule
 457  */
 458
 459 AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule()
 460     : impl_(new internal::BasicHistogramImpl())
 461 {
 462 }
 463
 464
 465 AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule(
 466         const AnalysisHistogramSettings &settings)
 467     : impl_(new internal::BasicHistogramImpl(settings))
 468 {
 469 }
 470
 471
 472 AnalysisDataSimpleHistogramModule::~AnalysisDataSimpleHistogramModule()
 473 {
 474 }
 475
 476
 477 void AnalysisDataSimpleHistogramModule::init(const AnalysisHistogramSettings &settings)
 478 {
 479     impl_->init(settings);
 480 }
 481
 482
 483 AbstractAverageHistogram &
 484 AnalysisDataSimpleHistogramModule::averager()
 485 {
 486     return *impl_->averager_;
 487 }
 488
 489
 490 const AnalysisHistogramSettings &
 491 AnalysisDataSimpleHistogramModule::settings() const
 492 {
 493     return impl_->settings_;
 494 }
 495
 496
 497 int
 498 AnalysisDataSimpleHistogramModule::flags() const
 499 {
 500     return efAllowMulticolumn | efAllowMultipoint;
 501 }
 502
 503
 504 void
 505 AnalysisDataSimpleHistogramModule::dataStarted(AbstractAnalysisData *data)
 506 {
 507     addModule(impl_->averager_);
 508     setColumnCount(settings().binCount());
 509     notifyDataStart();
 510     impl_->storage_.startDataStorage(this);
 511 }
 512
 513
 514 void
 515 AnalysisDataSimpleHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
 516 {
 517     AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
 518     impl_->initFrame(&frame);
 519 }
 520
 521
 522 void
 523 AnalysisDataSimpleHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
 524 {
 525     AnalysisDataStorageFrame &frame =
 526         impl_->storage_.currentFrame(points.frameIndex());
 527     for (int i = 0; i < points.columnCount(); ++i)
 528     {
 529         int bin = settings().findBin(points.y(i));
 530         if (bin != -1)
 531         {
 532             frame.value(bin) += 1;
 533         }
 534     }
 535 }
 536
 537
 538 void
 539 AnalysisDataSimpleHistogramModule::frameFinished(const AnalysisDataFrameHeader &header)
 540 {
 541     impl_->storage_.finishFrame(header.index());
 542 }
 543
 544
 545 void
 546 AnalysisDataSimpleHistogramModule::dataFinished()
 547 {
 548     notifyDataFinish();
 549 }
 550
 551
 552 AnalysisDataFrameRef
 553 AnalysisDataSimpleHistogramModule::tryGetDataFrameInternal(int index) const
 554 {
 555     return impl_->storage_.tryGetDataFrame(index);
 556 }
 557
 558
 559 bool
 560 AnalysisDataSimpleHistogramModule::requestStorageInternal(int nframes)
 561 {
 562     return impl_->storage_.requestStorage(nframes);
 563 }
 564
 565
 566 /********************************************************************
 567  * AnalysisDataWeightedHistogramModule
 568  */
 569
 570 AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule()
 571     : impl_(new internal::BasicHistogramImpl())
 572 {
 573 }
 574
 575
 576 AnalysisDataWeightedHistogramModule::AnalysisDataWeightedHistogramModule(
 577         const AnalysisHistogramSettings &settings)
 578     : impl_(new internal::BasicHistogramImpl(settings))
 579 {
 580 }
 581
 582
 583 AnalysisDataWeightedHistogramModule::~AnalysisDataWeightedHistogramModule()
 584 {
 585 }
 586
 587
 588 void AnalysisDataWeightedHistogramModule::init(const AnalysisHistogramSettings &settings)
 589 {
 590     impl_->init(settings);
 591 }
 592
 593
 594 AbstractAverageHistogram &
 595 AnalysisDataWeightedHistogramModule::averager()
 596 {
 597     return *impl_->averager_;
 598 }
 599
 600
 601 const AnalysisHistogramSettings &
 602 AnalysisDataWeightedHistogramModule::settings() const
 603 {
 604     return impl_->settings_;
 605 }
 606
 607
 608 int
 609 AnalysisDataWeightedHistogramModule::flags() const
 610 {
 611     return efAllowMulticolumn | efAllowMultipoint;
 612 }
 613
 614
 615 void
 616 AnalysisDataWeightedHistogramModule::dataStarted(AbstractAnalysisData *data)
 617 {
 618     addModule(impl_->averager_);
 619     setColumnCount(settings().binCount());
 620     notifyDataStart();
 621     impl_->storage_.startDataStorage(this);
 622 }
 623
 624
 625 void
 626 AnalysisDataWeightedHistogramModule::frameStarted(const AnalysisDataFrameHeader &header)
 627 {
 628     AnalysisDataStorageFrame &frame = impl_->storage_.startFrame(header);
 629     impl_->initFrame(&frame);
 630 }
 631
 632
 633 void
 634 AnalysisDataWeightedHistogramModule::pointsAdded(const AnalysisDataPointSetRef &points)
 635 {
 636     if (points.firstColumn() != 0 || points.columnCount() < 2)
 637     {
 638         GMX_THROW(APIError("Invalid data layout"));
 639     }
 640     int bin = settings().findBin(points.y(0));
 641     if (bin != -1)
 642     {
 643         AnalysisDataStorageFrame &frame =
 644             impl_->storage_.currentFrame(points.frameIndex());
 645         for (int i = 1; i < points.columnCount(); ++i)
 646         {
 647             frame.value(bin) += points.y(i);
 648         }
 649     }
 650 }
 651
 652
 653 void
 654 AnalysisDataWeightedHistogramModule::frameFinished(const AnalysisDataFrameHeader &header)
 655 {
 656     impl_->storage_.finishFrame(header.index());
 657 }
 658
 659
 660 void
 661 AnalysisDataWeightedHistogramModule::dataFinished()
 662 {
 663     notifyDataFinish();
 664 }
 665
 666
 667 AnalysisDataFrameRef
 668 AnalysisDataWeightedHistogramModule::tryGetDataFrameInternal(int index) const
 669 {
 670     return impl_->storage_.tryGetDataFrame(index);
 671 }
 672
 673
 674 bool
 675 AnalysisDataWeightedHistogramModule::requestStorageInternal(int nframes)
 676 {
 677     return impl_->storage_.requestStorage(nframes);
 678 }
 679
 680
 681 /********************************************************************
 682  * AnalysisDataBinAverageModule
 683  */
 684
 685 AnalysisDataBinAverageModule::AnalysisDataBinAverageModule()
 686 {
 687     setColumnCount(3);
 688 }
 689
 690
 691 AnalysisDataBinAverageModule::AnalysisDataBinAverageModule(
 692         const AnalysisHistogramSettings &settings)
 693     : settings_(settings)
 694 {
 695     setColumnCount(3);
 696     setRowCount(settings.binCount());
 697     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
 698              settings.binWidth());
 699 }
 700
 701
 702 AnalysisDataBinAverageModule::~AnalysisDataBinAverageModule()
 703 {
 704 }
 705
 706
 707 void
 708 AnalysisDataBinAverageModule::init(const AnalysisHistogramSettings &settings)
 709 {
 710     settings_ = settings;
 711     setRowCount(settings.binCount());
 712     setXAxis(settings.firstEdge() + 0.5 * settings.binWidth(),
 713              settings.binWidth());
 714 }
 715
 716
 717 int
 718 AnalysisDataBinAverageModule::flags() const
 719 {
 720     return efAllowMulticolumn | efAllowMultipoint;
 721 }
 722
 723
 724 void
 725 AnalysisDataBinAverageModule::dataStarted(AbstractAnalysisData * /*data*/)
 726 {
 727     allocateValues();
 728 }
 729
 730
 731 void
 732 AnalysisDataBinAverageModule::frameStarted(const AnalysisDataFrameHeader & /*header*/)
 733 {
 734 }
 735
 736
 737 void
 738 AnalysisDataBinAverageModule::pointsAdded(const AnalysisDataPointSetRef &points)
 739 {
 740     if (points.firstColumn() != 0 || points.columnCount() < 2)
 741     {
 742         GMX_THROW(APIError("Invalid data layout"));
 743     }
 744     int bin = settings().findBin(points.y(0));
 745     if (bin != -1)
 746     {
 747         for (int i = 1; i < points.columnCount(); ++i)
 748         {
 749             real y = points.y(i);
 750             value(bin, 0) += y;
 751             value(bin, 1) += y * y;
 752         }
 753         value(bin, 2) += points.columnCount() - 1;
 754     }
 755 }
 756
 757
 758 void
 759 AnalysisDataBinAverageModule::frameFinished(const AnalysisDataFrameHeader & /*header*/)
 760 {
 761 }
 762
 763
 764 void
 765 AnalysisDataBinAverageModule::dataFinished()
 766 {
 767     for (int i = 0; i < settings().binCount(); ++i)
 768     {
 769         real n = value(i, 2);
 770         if (n > 0)
 771         {
 772             real ave = value(i, 0) / n;
 773             real std = sqrt(value(i, 1) / n - ave * ave);
 774             setValue(i, 0, ave);
 775             setValue(i, 1, std);
 776         }
 777     }
 778     valuesReady();
 779 }
 780
 781 } // namespace gmx