Fix histogram resampling output bin positions
[alexxy/gromacs.git] / src / gromacs / analysisdata / modules / histogram.h
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2010,2011,2012,2013,2014, by the GROMACS development team, led by
5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6  * and including many others, as listed in the AUTHORS file in the
7  * top-level source directory and at http://www.gromacs.org.
8  *
9  * GROMACS is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public License
11  * as published by the Free Software Foundation; either version 2.1
12  * of the License, or (at your option) any later version.
13  *
14  * GROMACS is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17  * Lesser General Public License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with GROMACS; if not, see
21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
23  *
24  * If you want to redistribute modifications to GROMACS, please
25  * consider that scientific software is very special. Version
26  * control is crucial - bugs must be traceable. We will be happy to
27  * consider code for inclusion in the official distribution, but
28  * derived work must not be called official GROMACS. Details are found
29  * in the README & COPYING files - if they are missing, get the
30  * official version at http://www.gromacs.org.
31  *
32  * To help us fund GROMACS development, we humbly ask that you cite
33  * the research papers on the package. Check out http://www.gromacs.org.
34  */
35 /*! \file
36  * \brief
37  * Declares analysis data modules for calculating histograms.
38  *
39  * \author Teemu Murtola <teemu.murtola@gmail.com>
40  * \inpublicapi
41  * \ingroup module_analysisdata
42  */
43 #ifndef GMX_ANALYSISDATA_MODULES_HISTOGRAM_H
44 #define GMX_ANALYSISDATA_MODULES_HISTOGRAM_H
45
46 #include <boost/shared_ptr.hpp>
47
48 #include "gromacs/analysisdata/abstractdata.h"
49 #include "gromacs/analysisdata/arraydata.h"
50 #include "gromacs/analysisdata/datamodule.h"
51
52 namespace gmx
53 {
54
55 class AnalysisHistogramSettings;
56
57 /*! \brief
58  * Provides "named parameter" idiom for constructing histograms.
59  *
60  * \see histogramFromBins()
61  * \see histogramFromRange()
62  *
63  * Methods in this class do not throw.
64  *
65  * \inpublicapi
66  * \ingroup module_analysisdata
67  */
68 class AnalysisHistogramSettingsInitializer
69 {
70     public:
71         /*! \brief
72          * Creates an empty initializer.
73          *
74          * Should not be called directly, but histogramFromRange() or
75          * histogramFromBins() should be used instead.
76          */
77         AnalysisHistogramSettingsInitializer();
78
79         /*! \brief
80          * Sets the first bin location.
81          *
82          * Typically should not be called directly, but through
83          * histogramFromBins().
84          */
85         AnalysisHistogramSettingsInitializer &start(real min)
86         { min_ = min; return *this; }
87         /*! \brief
88          * Sets the number of bins in the histogram.
89          *
90          * If only the first bin location is specified, this value is required
91          * (and automatically provided if histogramFromBins() is used).
92          * If both the first and last bins are specified, either this value or
93          * binWidth() is required.
94          */
95         AnalysisHistogramSettingsInitializer &binCount(int binCount)
96         { binCount_ = binCount; return *this; }
97         /*! \brief
98          * Sets the first and last bin locations.
99          *
100          * Typically should not be called directly, but through
101          * histogramFromRange().
102          */
103         AnalysisHistogramSettingsInitializer &range(real min, real max)
104         { min_ = min; max_ = max; return *this; }
105         /*! \brief
106          * Sets the bin width of the histogram.
107          *
108          * If only the first bin location is specified, this value is required
109          * (and automatically provided if histogramFromBins() is used).
110          * If both the first and last bins are specified, either this value or
111          * binCount() is required.
112          * If a bin width is provided with both first and last bin locations,
113          * and the given bin width does not divide the range exactly, the last
114          * bin location is adjusted to match.
115          */
116         AnalysisHistogramSettingsInitializer &binWidth(real binWidth)
117         { binWidth_ = binWidth; return *this; }
118         /*! \brief
119          * Indicate that first and last bin locations to specify bin centers.
120          *
121          * If set, the first and last bin locations are interpreted as bin
122          * centers.
123          * If not set (the default), the first and last bin locations are
124          * interpreted as the edges of the whole histogram.
125          *
126          * Cannot be specified together with roundRange().
127          */
128         AnalysisHistogramSettingsInitializer &integerBins(bool enabled = true)
129         { bIntegerBins_ = enabled; return *this; }
130         /*! \brief
131          * Round first and last bin locations.
132          *
133          * If set, the resulting histogram will cover the range specified, but
134          * the actual bin locations will be rounded such that the edges fall
135          * on multiples of the bin width.
136          * Only implemented when both first and last bin location and bin width
137          * are defined.
138          * Cannot be specified together with integerBins() or with binCount().
139          */
140         AnalysisHistogramSettingsInitializer &roundRange(bool enabled = true)
141         { bRoundRange_ = enabled; return *this; }
142         /*! \brief
143          * Sets the histogram to match all values.
144          *
145          * If set, the histogram behaves as if the bins at the ends extended to
146          * +-infinity.
147          */
148         AnalysisHistogramSettingsInitializer &includeAll(bool enabled = true)
149         { bIncludeAll_ = enabled; return *this; }
150
151     private:
152         real                    min_;
153         real                    max_;
154         real                    binWidth_;
155         int                     binCount_;
156         bool                    bIntegerBins_;
157         bool                    bRoundRange_;
158         bool                    bIncludeAll_;
159
160         friend class AnalysisHistogramSettings;
161 };
162
163 /*! \brief
164  * Initializes a histogram using a range and a bin width.
165  *
166  * Does not throw.
167  *
168  * \inpublicapi
169  */
170 inline AnalysisHistogramSettingsInitializer
171 histogramFromRange(real min, real max)
172 {
173     return AnalysisHistogramSettingsInitializer().range(min, max);
174 }
175
176 /*! \brief
177  * Initializes a histogram using bin width and the number of bins.
178  *
179  * Does not throw.
180  *
181  * \inpublicapi
182  */
183 inline AnalysisHistogramSettingsInitializer
184 histogramFromBins(real start, int nbins, real binwidth)
185 {
186     return AnalysisHistogramSettingsInitializer()
187                .start(start).binCount(nbins).binWidth(binwidth);
188 }
189
190
191 /*! \brief
192  * Contains parameters that specify histogram bin locations.
193  *
194  * Methods in this class do not throw.
195  *
196  * \inpublicapi
197  * \ingroup module_analysisdata
198  */
199 class AnalysisHistogramSettings
200 {
201     public:
202         //! Initializes undefined parameters.
203         AnalysisHistogramSettings();
204         /*! \brief
205          * Initializes parameters based on a named parameter object.
206          *
207          * This constructor is not explicit to allow initialization of
208          * histograms directly from AnalysisHistogramSettingsInitializer:
209          * \code
210            gmx::AnalysisDataSimpleHistogramModule *hist =
211                    new gmx::AnalysisDataSimpleHistogramModule(
212                            histogramFromRange(0.0, 5.0).binWidth(0.5));
213          * \endcode
214          */
215         AnalysisHistogramSettings(const AnalysisHistogramSettingsInitializer &settings);
216
217         //! Returns the left edge of the first bin.
218         real firstEdge() const { return firstEdge_; }
219         //! Returns the right edge of the first bin.
220         real lastEdge() const { return lastEdge_; }
221         //! Returns the number of bins in the histogram.
222         int binCount() const { return binCount_; }
223         //! Returns the width of a bin in the histogram.
224         real binWidth() const { return binWidth_; }
225         //! Whether values beyond the edges are mapped to the edge bins.
226         bool includeAll() const { return bAll_; }
227         //! Returns a zero-based bin index for a value, or -1 if not in range.
228         int findBin(real y) const;
229
230     private:
231         real                    firstEdge_;
232         real                    lastEdge_;
233         real                    binWidth_;
234         real                    inverseBinWidth_;
235         int                     binCount_;
236         bool                    bAll_;
237 };
238
239
240 class AbstractAverageHistogram;
241
242 //! Smart pointer to manage an AbstractAverageHistogram object.
243 typedef boost::shared_ptr<AbstractAverageHistogram>
244     AverageHistogramPointer;
245
246 /*! \brief
247  * Base class for representing histograms averaged over frames.
248  *
249  * The averaging module for a per-frame histogram is always created by the
250  * histogram module class (e.g., AnalysisDataSimpleHistogramModule), and can be
251  * accessed using, e.g., AnalysisDataSimpleHistogramModule::averager().
252  * The user can alter some properties of the average histogram directly, but
253  * the main use of the object is to postprocess the histogram once the
254  * calculation is finished.
255  *
256  * This class can represent multiple histograms in one object: each column in
257  * the data is an independent histogram.
258  * The X values correspond to center of the bins, except for a cumulative
259  * histogram made with makeCumulative().
260  *
261  * \inpublicapi
262  * \ingroup module_analysisdata
263  */
264 class AbstractAverageHistogram : public AbstractAnalysisArrayData
265 {
266     public:
267         virtual ~AbstractAverageHistogram();
268
269         //! Returns bin properties for the histogram.
270         const AnalysisHistogramSettings &settings() const { return settings_; }
271
272         /*! \brief
273          * Creates a copy of the histogram with double the bin width.
274          *
275          * \param[in] bIntegerBins If `true`, the first bin in the result will
276          *     cover the first bin from the source. Otherwise, the first bin
277          *     will cover first two bins from the source.
278          * \throws std::bad_alloc if out of memory.
279          *
280          * The caller is responsible of deleting the returned object.
281          */
282         AverageHistogramPointer resampleDoubleBinWidth(bool bIntegerBins) const;
283         /*! \brief
284          * Creates a deep copy of the histogram.
285          *
286          * \throws std::bad_alloc if out of memory.
287          *
288          * The returned histogram is not necessarily of the same dynamic type
289          * as the original object, but contains the same data from the point of
290          * view of the AbstractAverageHistogram interface.
291          *
292          * The caller is responsible of deleting the returned object.
293          */
294         AverageHistogramPointer clone() const;
295         //! Normalizes the histogram such that the integral over it is one.
296         void normalizeProbability();
297         /*! \brief
298          * Makes the histograms cumulative by summing up each bin to all bins
299          * after it.
300          *
301          * The X values in the data are adjusted such that they match the right
302          * edges of bins instead of bin centers.
303          */
304         void makeCumulative();
305         //! Scales a single histogram by a uniform scaling factor.
306         void scaleSingle(int index, real factor);
307         //! Scales all histograms by a uniform scaling factor.
308         void scaleAll(real factor);
309         //! Scales the value of each bin by a different scaling factor.
310         void scaleAllByVector(real factor[]);
311         /*! \brief
312          * Notifies attached modules of the histogram data.
313          *
314          * After this function has been called, it is no longer possible to
315          * alter the histogram.
316          */
317         void done() { AbstractAnalysisArrayData::valuesReady(); }
318
319     protected:
320         /*! \brief
321          * Creates a histogram module with undefined bins.
322          *
323          * Bin parameters must be defined with init() before data input is
324          * started.
325          */
326         AbstractAverageHistogram();
327         //! Creates a histogram module with defined bin parameters.
328         explicit AbstractAverageHistogram(const AnalysisHistogramSettings &settings);
329
330         /*! \brief
331          * (Re)initializes the histogram from settings.
332          */
333         void init(const AnalysisHistogramSettings &settings);
334
335     private:
336         AnalysisHistogramSettings  settings_;
337
338         // Copy and assign disallowed by base.
339 };
340
341
342 /*! \brief
343  * Data module for per-frame histograms.
344  *
345  * Output data contains the same number of frames and data sets as the input
346  * data.  Each frame contains the histogram(s) for the points in that frame.
347  * Each input data set is processed independently into the corresponding output
348  * data set.  Missing values are ignored.
349  * All input columns for a data set are averaged into the same histogram.
350  * The number of columns for all data sets equals the number of bins in the
351  * histogram.
352  *
353  * The histograms are accumulated as 64-bit integers within a frame and summed
354  * in double precision across frames, even if the output data is in single
355  * precision.
356  *
357  * \inpublicapi
358  * \ingroup module_analysisdata
359  */
360 class AnalysisDataSimpleHistogramModule : public AbstractAnalysisData,
361                                           public AnalysisDataModuleParallel
362 {
363     public:
364         /*! \brief
365          * Creates a histogram module with undefined bins.
366          *
367          * Bin parameters must be defined with init() before data input is
368          * started.
369          */
370         AnalysisDataSimpleHistogramModule();
371         //! Creates a histogram module with defined bin parameters.
372         explicit AnalysisDataSimpleHistogramModule(const AnalysisHistogramSettings &settings);
373         virtual ~AnalysisDataSimpleHistogramModule();
374
375         /*! \brief
376          * (Re)initializes the histogram from settings.
377          */
378         void init(const AnalysisHistogramSettings &settings);
379
380         /*! \brief
381          * Returns the average histogram over all frames.
382          *
383          * Can be called already before the histogram is calculated to
384          * customize the way the average histogram is calculated.
385          *
386          * \see AbstractAverageHistogram
387          */
388         AbstractAverageHistogram &averager();
389
390         //! Returns bin properties for the histogram.
391         const AnalysisHistogramSettings &settings() const;
392
393         virtual int frameCount() const;
394
395         virtual int flags() const;
396
397         virtual bool parallelDataStarted(
398             AbstractAnalysisData              *data,
399             const AnalysisDataParallelOptions &options);
400         virtual void frameStarted(const AnalysisDataFrameHeader &header);
401         virtual void pointsAdded(const AnalysisDataPointSetRef &points);
402         virtual void frameFinished(const AnalysisDataFrameHeader &header);
403         virtual void dataFinished();
404
405     private:
406         virtual AnalysisDataFrameRef tryGetDataFrameInternal(int index) const;
407         virtual bool requestStorageInternal(int nframes);
408
409         class Impl;
410
411         PrivateImplPointer<Impl> impl_;
412
413         // Copy and assign disallowed by base.
414 };
415
416
417 /*! \brief
418  * Data module for per-frame weighted histograms.
419  *
420  * Output data contains the same number of frames and data sets as the input
421  * data.  Each frame contains the histogram(s) for the points in that frame,
422  * interpreted such that the first column passed to pointsAdded() determines
423  * the bin and the rest give weights to be added to that bin (input data should
424  * have at least two colums, and at least two columns should be added at the
425  * same time).
426  * Each input data set is processed independently into the corresponding output
427  * data set.
428  * All input columns for a data set are averaged into the same histogram.
429  * The number of columns for all data sets equals the number of bins in the
430  * histogram.
431  *
432  * The histograms are accumulated in double precision, even if the output data
433  * is in single precision.
434  *
435  * \inpublicapi
436  * \ingroup module_analysisdata
437  */
438 class AnalysisDataWeightedHistogramModule : public AbstractAnalysisData,
439                                             public AnalysisDataModuleParallel
440 {
441     public:
442         //! \copydoc AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule()
443         AnalysisDataWeightedHistogramModule();
444         //! \copydoc AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule(const AnalysisHistogramSettings &)
445         explicit AnalysisDataWeightedHistogramModule(const AnalysisHistogramSettings &settings);
446         virtual ~AnalysisDataWeightedHistogramModule();
447
448         //! \copydoc AnalysisDataSimpleHistogramModule::init()
449         void init(const AnalysisHistogramSettings &settings);
450
451         //! \copydoc AnalysisDataSimpleHistogramModule::averager()
452         AbstractAverageHistogram &averager();
453
454         //! \copydoc AnalysisDataSimpleHistogramModule::settings()
455         const AnalysisHistogramSettings &settings() const;
456
457         virtual int frameCount() const;
458
459         virtual int flags() const;
460
461         virtual bool parallelDataStarted(
462             AbstractAnalysisData              *data,
463             const AnalysisDataParallelOptions &options);
464         virtual void frameStarted(const AnalysisDataFrameHeader &header);
465         virtual void pointsAdded(const AnalysisDataPointSetRef &points);
466         virtual void frameFinished(const AnalysisDataFrameHeader &header);
467         virtual void dataFinished();
468
469     private:
470         virtual AnalysisDataFrameRef tryGetDataFrameInternal(int index) const;
471         virtual bool requestStorageInternal(int nframes);
472
473         class Impl;
474
475         PrivateImplPointer<Impl> impl_;
476
477         // Copy and assign disallowed by base.
478 };
479
480
481 /*! \brief
482  * Data module for bin averages.
483  *
484  * Output data contains one row for each bin; see AbstractAverageHistogram.
485  * Output data contains one column for each input data set.
486  * The value in a column is the average over all frames of that data set for
487  * that bin.
488  * The input data is interpreted such that the first column passed to
489  * pointsAdded() determines the bin and the rest give values to be added to
490  * that bin (input data should have at least two colums, and at least two
491  * columns should be added at the same time).
492  * All input columns for a data set are averaged into the same histogram.
493  *
494  * \inpublicapi
495  * \ingroup module_analysisdata
496  */
497 class AnalysisDataBinAverageModule : public AbstractAnalysisArrayData,
498                                      public AnalysisDataModuleSerial
499 {
500     public:
501         //! \copydoc AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule()
502         AnalysisDataBinAverageModule();
503         //! \copydoc AnalysisDataSimpleHistogramModule::AnalysisDataSimpleHistogramModule(const AnalysisHistogramSettings &)
504         explicit AnalysisDataBinAverageModule(const AnalysisHistogramSettings &settings);
505         virtual ~AnalysisDataBinAverageModule();
506
507         //! \copydoc AnalysisDataSimpleHistogramModule::init()
508         void init(const AnalysisHistogramSettings &settings);
509
510         //! \copydoc AnalysisDataSimpleHistogramModule::settings()
511         const AnalysisHistogramSettings &settings() const;
512
513         virtual int flags() const;
514
515         virtual void dataStarted(AbstractAnalysisData *data);
516         virtual void frameStarted(const AnalysisDataFrameHeader &header);
517         virtual void pointsAdded(const AnalysisDataPointSetRef &points);
518         virtual void frameFinished(const AnalysisDataFrameHeader &header);
519         virtual void dataFinished();
520
521     private:
522         class Impl;
523
524         PrivateImplPointer<Impl>   impl_;
525
526         // Copy and assign disallowed by base.
527 };
528
529 //! Smart pointer to manage an AnalysisDataSimpleHistogramModule object.
530 typedef boost::shared_ptr<AnalysisDataSimpleHistogramModule>
531     AnalysisDataSimpleHistogramModulePointer;
532 //! Smart pointer to manage an AnalysisDataWeightedHistogramModule object.
533 typedef boost::shared_ptr<AnalysisDataWeightedHistogramModule>
534     AnalysisDataWeightedHistogramModulePointer;
535 //! Smart pointer to manage an AnalysisDataBinAverageModule object.
536 typedef boost::shared_ptr<AnalysisDataBinAverageModule>
537     AnalysisDataBinAverageModulePointer;
538
539 } // namespace gmx
540
541 #endif