2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2009,2010,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.
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.
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.
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.
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.
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.
35 /*! \page histograms Histogram calculation
37 * Functions to calculate histograms are defined in histogram.h.
38 * The principle is simple: histogram calculation is set up with
39 * gmx_histogram_create() and gmx_histogram_set_*(), the histogram is sampled
40 * using the provided functions, and gmx_histogram_finish() is called to
41 * finish the sampling.
42 * Various post-processing functions can then be used to normalize the
43 * histogram in various ways and to write it into a file.
44 * gmx_histogram_get_*() can be used to access data from the histogram.
45 * gmx_histogram_free() can be used to free the memory allocated for a
49 * \section histogram_sampling Initialization and sampling
51 * A histogram calculation is initialized by calling gmx_histogram_create()
52 * or gmx_histogram_create_range().
53 * The first function takes the type of the histogram (see \ref e_histogram_t)
54 * and the number of bins, and allocates the necessary memory.
55 * The bin locations can be initialized with gmx_histogram_set_integerbins()
56 * and one of gmx_histogram_set_binwidth() and gmx_histogram_set_range().
57 * Only uniformly spaced bins are currently supported.
58 * The second initialization function takes the beginning and end of the range
60 * The treatment of values that fall outside the histogram range can be
61 * set with gmx_histogram_set_all().
62 * gmx_histogram_set_blocksize() can be used to specify a block size for error
63 * estimates. If not called, no error estimates are calculated.
64 * Otherwise, a histogram is calculated for each block of given number of
65 * frames, and the error estimated as the error of the mean of these block
67 * If the number of frames is not divisible by the block size,
68 * the last block is discarded.
69 * gmx_histogram_set_block_output() can be used to write the individual
70 * block histograms to a file (it is currently not very flexible).
72 * After initialization, three sets of sampling functions are provided
73 * to sample the differen types of histograms:
74 * - Use gmx_histogram_increment() or gmx_histogram_increment_bin() to
75 * calculate simple histograms (\ref HIST_SIMPLE histograms).
76 * - Use gmx_histogram_add() or gmx_histogram_add_to_bin() to calculate
77 * histograms where different points can have different weight
78 * (\ref HIST_WEIGHT histograms).
79 * - Use gmx_histogram_add_item() or gmx_histogram_add_item_to_bin() to
80 * calculate averages within each bin
81 * (\ref HIST_BINAVER histograms).
83 * The functions aboge that take bin indices directly do no range checking,
84 * i.e., one must ensure that the bin number is between 0 and the number of
85 * bins in the histogram. In contrast, the functions that use real values to
86 * locate the correct bin use range checking, and silently discard values
87 * outside the histogram (unless gmx_histogram_set_all() has been called).
88 * gmx_histogram_find_bin() can be used to find the bin number corresponding
89 * to a value, adhering to the gmx_histogram_set_all() setting.
91 * After a frame is sampled, gmx_histogram_finish_frame() needs to be
93 * After all data is sampled, gmx_histogram_finish() should be called.
94 * After these calls, the histogram will be normalized such that the
95 * sum over all bins equals the average number of samples per frame, and
96 * the error field contains the error of the mean at each bin.
99 * \section histogram_processing Post-processing
101 * The following functions can be used to process the histograms:
102 * - gmx_histogram_clone() makes a deep copy.
103 * - gmx_histogram_resample_dblbw() also makes a deep copy,
104 * but uses a binwidth that is double of the input one
105 * (see below for usage).
106 * - gmx_histogram_normalize_prob() normalizes a histogram such that its
107 * integral becomes one.
108 * - gmx_histogram_scale() and gmx_histogram_scale_vec() can be used
109 * to scale a histogram uniformly or non-uniformly, respectively.
111 * There are also functions to write histograms to a file:
112 * - gmx_histogram_write() writes a single histogram, optionally with
114 * - gmx_histogram_write_array() writes multiple histograms with identical
115 * bins. Values and/or errors can be written.
116 * - gmx_histogram_write_cum_array() writes cumulative histograms for
117 * multiple histograms with identical bins.
119 * gmx_histogram_resample_dblbw() is useful if a histogram and its
120 * cumulative sum are required at same points.
121 * One can then sample the histogram with half the desired binwidth, and
122 * then use gmx_histogram_resample_dblbw() to construct two different
123 * versions of the histogram (by changing \p bIntegerBins), one suitable for
124 * writing out the histogram and the other for writing the cumulative
127 * Access to the histogram data is also provided through the following
129 * - gmx_histogram_get_nbins()
130 * - gmx_histogram_get_binwidth()
131 * - gmx_histogram_get_value()
132 * - gmx_histogram_get_bin_value()
133 * - gmx_histogram_get_values()
134 * - gmx_histogram_get_errors()
136 * gmx_histogram_free() can be used to free memory associated with a
137 * histogram when it is no longer needed.
140 * \brief Implementation of functions in histogram.h.
152 #include "histogram.h"
155 * Stores data for a histogram.
157 struct gmx_histogram_t
159 /** The left edge of the first bin. */
161 /** The right edge of the last bin. */
165 /** Number of bins. */
167 /** The final histogram. */
169 /** Standard deviation of each bin (calculated from block histograms). */
172 /** Histogram type. */
174 /** Histogram flags; */
176 /** Block size for averaging. */
178 /** Output file for block histograms (can be NULL). */
181 /** Inverse bin width. */
183 /** Current histogram value. */
185 /** Current histogram count. */
187 /** Number of frames read for the current block so far. */
189 /** Number of blocks averaged so far. */
194 * \param[out] hp Histogram to create.
195 * \param[in] type Histogram type to create;
196 * \param[in] nbins Number of bins.
197 * \returns 0 on success, a non-zero error code on error.
199 * Initialized the histogram structure \p h with the provided values,
200 * allocating the necessary memory.
203 gmx_histogram_create(gmx_histogram_t **hp, e_histogram_t type, int nbins)
210 gmx_incons("number of histogram bins <= 0");
233 snew(h->hist, nbins + 1);
234 snew(h->histerr, nbins + 1);
235 snew(h->cn, nbins + 1);
236 if (type != HIST_SIMPLE)
238 snew(h->chist, nbins + 1);
240 gmx_histogram_clear(h);
247 * \param[out] hp Histogram to create.
248 * \param[in] type Histogram type to create;
249 * \param[in] start Left edge/center of the first bin
250 * (depending on \p bIntegerBins).
251 * \param[in] end Right edge/center of the last bin
252 * (depending on \p bIntegerBins).
253 * \param[in] binw Bin width.
254 * \param[in] bIntegerBins If true, histogram bins are centered at
255 * \c start+n*binwidth, otherwise the centers are at
256 * \c start+(n+0.5)*binwidth.
257 * \returns 0 on success, a non-zero error code on error.
259 * Initialized the histogram structure \p h with the provided values,
260 * allocating the necessary memory.
263 gmx_histogram_create_range(gmx_histogram_t **hp, e_histogram_t type,
264 real start, real end, real binw, gmx_bool bIntegerBins)
271 if (start >= end || binw <= 0)
273 gmx_incons("histogram left edge larger than right edge or bin width <= 0");
277 /* Adjust the end points and calculate the number of bins */
280 nbins = ceil((end - start) / binw) + 1;
281 end = start + (nbins - 1) * binw;
285 start = binw * floor(start / binw);
286 end = binw * ceil(end / binw);
292 nbins = (int)((end - start) / binw + 0.5);
294 /* Create the histogram */
295 rc = gmx_histogram_create(&h, type, nbins);
301 gmx_histogram_set_integerbins(h, bIntegerBins);
302 gmx_histogram_set_binwidth(h, start, binw);
309 * \param[in,out] h Histogram to clear.
311 * Histograms are automatically cleared when initialized; you only need to
312 * call this function if you want to reuse a histogram structure that has
313 * already been used for sampling.
316 gmx_histogram_clear(gmx_histogram_t *h)
324 for (i = 0; i <= h->nbins; ++i)
339 * \param[in] h Histogram to free.
341 * The pointer \p h is invalid after the call.
344 gmx_histogram_free(gmx_histogram_t *h)
354 * \param[in,out] h Histogram data structure.
355 * \param[in] start Left edge/center of the first bin
356 * (depending on gmx_histogram_set_integerbins()).
357 * \param[in] binw Bin width.
358 * \returns 0 on success, a non-zero error code on error.
361 gmx_histogram_set_binwidth(gmx_histogram_t *h, real start, real binw)
365 gmx_incons("histogram binwidth <= 0");
368 if (h->flags & HIST_INTEGERBINS)
374 h->end = start + h->nbins * binw;
376 h->flags |= HIST_INITBW;
381 * \param[in,out] h Histogram data structure.
382 * \param[in] start Left edge/center of the first bin
383 * (depending on gmx_histogram_set_integerbins()).
384 * \param[in] end Right edge/center of the last bin
385 * (depending on gmx_histogram_set_integerbins()).
386 * \returns 0 on success, a non-zero error code on error.
389 gmx_histogram_set_range(gmx_histogram_t *h, real start, real end)
393 gmx_incons("histogram left edge larger than right edge");
398 if (h->flags & HIST_INTEGERBINS)
400 h->binwidth = (end - start) / (h->nbins - 1);
401 start -= 0.5*h->binwidth;
402 end += 0.5*h->binwidth;
406 h->binwidth = (end - start) / h->nbins;
408 h->invbw = 1.0/h->binwidth;
409 h->flags &= ~HIST_INITBW;
414 * \param[in,out] h Histogram data structure.
415 * \param[in] bIntegerBins If true, histogram bins are centered at
416 * \c start+n*binwidth, otherwise the centers are at
417 * \c start+(n+0.5)*binwidth.
420 gmx_histogram_set_integerbins(gmx_histogram_t *h, gmx_bool bIntegerBins)
422 /* Adjust the ranges if they have been initialized */
423 if (h->start < h->end)
425 if (h->flags & HIST_INTEGERBINS)
427 h->start += 0.5*h->binwidth;
428 if (h->flags & HIST_INITBW)
430 h->end += 0.5*h->binwidth;
434 h->end -= 0.5*h->binwidth;
439 h->start -= 0.5*h->binwidth;
440 if (h->flags & HIST_INITBW)
442 h->end -= 0.5*h->binwidth;
446 h->end += 0.5*h->binwidth;
452 h->flags |= HIST_INTEGERBINS;
456 h->flags &= ~HIST_INTEGERBINS;
461 * \param[in,out] h Histogram data structure.
462 * \param[in] bAll If true, values outside the histogram edges are added
463 * to the bins at the edges.
465 * \p bAll can be used to avoid rounding errors in cases where the histogram
466 * spans the full range of possible values. If set, the values that are at
467 * the exact maximum are still correctly included.
470 gmx_histogram_set_all(gmx_histogram_t *h, gmx_bool bAll)
474 h->flags |= HIST_ALL;
478 h->flags &= ~HIST_ALL;
483 * \param[in,out] h Histogram data structure.
484 * \param[in] bsize Block size for error estimates.
485 * \returns 0 on success, a non-zero error code on error.
487 * If \p bsize is zero, no block averaging or error estimates are done.
488 * This is also the case if this function is not called.
491 gmx_histogram_set_blocksize(gmx_histogram_t *h, int bsize)
495 gmx_incons("histogram block size < 0");
503 * \param[in,out] h Histogram data structure.
504 * \param[in] fp File for block output.
505 * \returns 0 on success, a non-zero error code on error.
507 * Sets a file into which each block histogram (the histogram after each
508 * \c gmx_histogram_t::bsize frames) is written.
509 * All histograms are written to the same file, separated by empty lines.
510 * If this function is not called, the block histograms are only used for
511 * error estimation, not written to a file.
514 gmx_histogram_set_block_output(gmx_histogram_t *h, FILE *fp)
518 gmx_incons("histogram block size not set but output initialized");
526 * \param[in] h Histogram data.
527 * \param[in] pos Position.
528 * \returns Bin index in \p h corresponding to \p pos, or -1 if there is no
531 * If gmx_histogram_set_all() has been called with TRUE, values outside the
532 * histogram are mapped to the bins at the edges.
535 gmx_histogram_find_bin(gmx_histogram_t *h, real pos)
539 if (h->flags & HIST_ALL)
548 else if (pos >= h->end)
550 if (h->flags & HIST_ALL)
559 /* Calculate the bin index if we are inside the box */
560 return (int)((pos - h->start)*h->invbw);
564 * \param[in] h Histogram data.
565 * \returns Number of bins in \p h.
568 gmx_histogram_get_nbins(gmx_histogram_t *h)
574 * \param[in] h Histogram data.
575 * \returns Bin width of \p h, 0 if no binwidth has been set.
578 gmx_histogram_get_binwidth(gmx_histogram_t *h)
584 * \param[in] h Histogram data.
585 * \param[in] pos Position.
586 * \param[out] value Pointer to receive the value (can be NULL).
587 * \param[out] err Pointer to receive the value (can be NULL).
589 * If \p pos is outside the range of the histogram, zeros are returned.
592 gmx_histogram_get_value(gmx_histogram_t *h, real pos, double *value, double *err)
597 if (pos < h->start || pos > h->end)
603 bin = gmx_histogram_find_bin(h, pos);
625 * \param[in] h Histogram data.
626 * \param[in] bin Bin number.
627 * \param[out] value Pointer to receive the value (can be NULL).
628 * \param[out] err Pointer to receive the value (can be NULL).
630 * If \p bin is outside the valid range, zeros are returned.
633 gmx_histogram_get_bin_value(gmx_histogram_t *h, int bin, double *value, double *err)
637 if (bin < 0 || bin >= h->nbins)
657 * \param[in] h Histogram data.
658 * \returns Pointer to an array of values for \p h.
660 * The returned array has one element for each bin of \p h.
661 * The returned pointer should not be freed.
664 gmx_histogram_get_values(gmx_histogram_t *h)
670 * \param[in] h Histogram data.
671 * \returns Pointer to an array of errors for \p h.
673 * The returned array has one element for each bin of \p h.
674 * The returned pointer should not be freed.
677 gmx_histogram_get_errors(gmx_histogram_t *h)
683 * \param[in,out] h Histogram data.
684 * \param[in] pos Position.
687 gmx_histogram_increment(gmx_histogram_t *h, real pos)
689 int bin = gmx_histogram_find_bin(h, pos);
698 * \param[in,out] h Histogram data.
699 * \param[in] pos Position.
700 * \param[in] value Value to add.
703 gmx_histogram_add(gmx_histogram_t *h, real pos, double value)
705 int bin = gmx_histogram_find_bin(h, pos);
710 h->chist[bin] += value;
714 * \param[in,out] h Histogram data.
715 * \param[in] pos Position.
716 * \param[in] value Value to add.
719 gmx_histogram_add_item(gmx_histogram_t *h, real pos, double value)
721 int bin = gmx_histogram_find_bin(h, pos);
726 h->chist[bin] += value;
731 * \param[in,out] h Histogram data.
732 * \param[in] bin Bin number.
734 * No checks for out-of-bound errors are done: \p bin should be >=0 and
738 gmx_histogram_increment_bin(gmx_histogram_t *h, int bin)
744 * \param[in,out] h Histogram data.
745 * \param[in] bin Bin number.
746 * \param[in] value Value to add.
748 * No checks for out-of-bound errors are done: \p bin should be >=0 and
752 gmx_histogram_add_to_bin(gmx_histogram_t *h, int bin, double value)
754 h->chist[bin] += value;
758 * \param[in,out] h Histogram data.
759 * \param[in] bin Bin number.
760 * \param[in] value Value to add.
762 * No checks for out-of-bound errors are done: \p bin should be >=0 and
766 gmx_histogram_add_item_to_bin(gmx_histogram_t *h, int bin, double value)
768 h->chist[bin] += value;
773 * Processes a sampled block histogram.
775 * \param[in,out] h Histogram data structure.
778 finish_histogram_block(gmx_histogram_t *h)
788 if (h->flags & HIST_ALL)
792 h->chist[h->nbins-1] += h->chist[h->nbins];
794 h->cn[h->nbins-1] += h->cn[h->nbins];
797 for (i = 0; i <= h->nbins; ++i)
801 v = h->chist[i] / (h->cn[i] > 0 ? h->cn[i] : h->nframes);
805 v = ((real)h->cn[i]) / h->nframes;
809 fprintf(h->blockfp, "%10g %10g\n", h->start + (i+0.5)*h->binwidth, v);
812 h->histerr[i] += sqr(v);
821 fprintf(h->blockfp, "\n");
828 * \param[in,out] h Histogram data structure.
830 * Should be called after each frame of data.
832 * \see gmx_histogram_finish()
835 gmx_histogram_finish_frame(gmx_histogram_t *h)
838 if (h->nframes == h->bsize)
840 finish_histogram_block(h);
845 * \param[in,out] h Histogram data structure.
847 * Normalizes the histogram.
848 * Should be called after all frames have been sampled and before any
849 * post-processing of the histogram.
850 * If block size has not been set with gmx_histogram_set_blocksize(),
851 * this function does no normalization, but it should still be called,
852 * otherwise the \c gmx_histogram_t::hist array will contain only zeros.
854 * \see gmx_histogram_finish_frame()
857 gmx_histogram_finish(gmx_histogram_t *h)
861 if (h->nframes > 0 || h->bsize == 0)
863 if (h->nframes < h->bsize)
865 fprintf(stderr, "Last block smaller (%d frames) than the target size (%d frames) skipped \n",
866 h->nframes, h->bsize);
870 finish_histogram_block(h);
878 for (i = 0; i <= h->nbins; ++i)
880 h->hist[i] /= h->nblocks;
881 h->histerr[i] /= h->nblocks;
882 h->histerr[i] = sqrt((h->histerr[i] - sqr(h->hist[i])) / h->nblocks);
887 * \param[out] destp Destination histogram.
888 * \param[in] src Source histogram.
889 * \param[in] bIntegerBins Control bin center position.
891 * Uses \p src to create a new histogram that has double the binwidth.
892 * If \p bIntegerBins is TRUE, the first new bin will be centered at
893 * \c src->start, otherwise the left edge will be at \c src->start.
895 * This function is mostly useful if one needs to sample both the
896 * histogram and the cumulative histogram at same points.
897 * To achieve this, first sample a histogram with half the desired
898 * binwidth, and then use this function to create two different verions
899 * of it with different values of \p bIntegerBins.
900 * gmx_histogram_write_array() and gmx_histogram_write_cum_array() can
901 * now used to write out the values correctly at identical values.
903 * Should be called only after gmx_histogram_finish() has been called.
904 * \p src should have been sampled without gmx_histogram_set_integerbins().
907 gmx_histogram_resample_dblbw(gmx_histogram_t **destp, gmx_histogram_t *src,
908 gmx_bool bIntegerBins)
910 gmx_histogram_t *dest;
914 gmx_histogram_create(destp, src->type, src->nbins / 2);
916 sfree(dest->chist); dest->chist = NULL;
917 sfree(dest->cn); dest->cn = NULL;
918 gmx_histogram_set_binwidth(dest, src->start, src->binwidth * 2);
919 gmx_histogram_set_integerbins(dest, bIntegerBins);
921 for (i = j = 0; i < dest->nbins; ++i)
923 if (bIntegerBins && i == 0)
926 ve = sqr(src->histerr[0]);
931 v = src->hist[j] + src->hist[j+1];
932 ve = sqr(src->histerr[j]) + sqr(src->histerr[j+1]);
937 dest->histerr[i] = ve;
939 dest->hist[dest->nbins] = 0;
940 dest->histerr[dest->nbins] = 0;
944 * \param[out] destp Destination histogram.
945 * \param[in] src Source histogram.
947 * Makes a clone of \p src for post-processing.
950 gmx_histogram_clone(gmx_histogram_t **destp, gmx_histogram_t *src)
952 gmx_histogram_t *dest;
955 memcpy(dest, src, sizeof(*dest));
957 /* These are not needed in post-processing */
958 dest->blockfp = NULL;
962 /* Make a deep copy of the actual histograms */
963 snew(dest->hist, src->nbins+1);
964 snew(dest->histerr, src->nbins+1);
965 memcpy(dest->hist, src->hist, (src->nbins+1)*sizeof(real));
966 memcpy(dest->histerr, src->histerr, (src->nbins+1)*sizeof(real));
972 * \param[in,out] h Histogram to normalize.
974 * Normalizes the histogram such that its integral equals one.
977 gmx_histogram_normalize_prob(gmx_histogram_t *h)
984 for (i = 0; i <= h->nbins; ++i)
989 normfac = h->invbw / sum;
990 gmx_histogram_scale(h, normfac);
994 * \param[in,out] h Histogram to normalize.
995 * \param[in] norm Scaling factor.
997 * All bin values are multiplied by \p norm.
1000 gmx_histogram_scale(gmx_histogram_t *h, real norm)
1004 for (i = 0; i <= h->nbins; ++i)
1007 h->histerr[i] *= norm;
1012 * \param[in,out] h Histogram to normalize.
1013 * \param[in] norm Scaling vector.
1015 * The i'th bin is multiplied by \p norm[i].
1018 gmx_histogram_scale_vec(gmx_histogram_t *h, real norm[])
1022 for (i = 0; i < h->nbins; ++i)
1024 h->hist[i] *= norm[i];
1025 h->histerr[i] *= norm[i];
1027 h->hist[h->nbins] *= norm[h->nbins-1];
1028 h->histerr[h->nbins] *= norm[h->nbins-1];
1032 * Makes some checks on output histograms and finds the maximum number of bins.
1034 * \param[in] n Number of histograms in \p h.
1035 * \param[in] h Array of histograms.
1036 * \param[out] nbins Pointer to a value that will receive the maximum number
1040 prepare_output(int n, gmx_histogram_t *h[], int *nbins)
1045 for (j = 0; j < n; ++j)
1047 if (!gmx_within_tol(h[j]->start, h[0]->start, GMX_REAL_EPS))
1049 fprintf(stderr, "gmx_ana_histogram_write: histogram start values not identical\n");
1051 if (!gmx_within_tol(h[j]->binwidth, h[0]->binwidth, GMX_REAL_EPS))
1053 fprintf(stderr, "gmx_ana_histogram_write: bin widths not identical\n");
1055 if (*nbins < h[j]->nbins)
1057 *nbins = h[j]->nbins;
1063 * \param[in] fp Output file.
1064 * \param[in] h Array of histograms to write.
1065 * \param[in] bErrors If TRUE, histogram errors are written.
1067 * Convenience wrapper for gmx_histogram_write_array() for writing a single
1070 * \see gmx_histogram_write_array()
1073 gmx_histogram_write(FILE *fp, gmx_histogram_t *h, gmx_bool bErrors)
1075 gmx_histogram_write_array(fp, 1, &h, TRUE, bErrors);
1079 * \param[in] fp Output file.
1080 * \param[in] n Number of histograms in \p h.
1081 * \param[in] h Array of histograms to write.
1082 * \param[in] bValue If TRUE, histogram values are written.
1083 * \param[in] bErrors If TRUE, histogram errors are written.
1085 * All the histograms in the array \p h should have the same bin widths and
1086 * left edges, otherwise the behavior is undefined.
1087 * The output format is one bin per line. The first number on the line is the
1088 * center of the bin, followed by one or two numbers for each histogram
1089 * (depending on \p bValue and \p bErrors).
1090 * If both \p bValue and \p bErrors are both TRUE, the values are written
1091 * before the errors.
1094 gmx_histogram_write_array(FILE *fp, int n, gmx_histogram_t *h[],
1095 gmx_bool bValue, gmx_bool bErrors)
1099 prepare_output(n, h, &nbins);
1100 for (i = 0; i < nbins; ++i)
1102 fprintf(fp, "%10g", h[0]->start + (i+0.5)*h[0]->binwidth);
1103 for (j = 0; j < n; ++j)
1107 fprintf(fp, " %10g", i < h[j]->nbins ? h[j]->hist[i] : 0.0);
1111 fprintf(fp, " %10g", i < h[j]->nbins ? h[j]->histerr[i] : 0.0);
1120 * \param[in] fp Output file.
1121 * \param[in] n Number of histograms in \p h.
1122 * \param[in] h Array of histograms to write.
1123 * \param[in] bValue If TRUE, histogram values are written.
1124 * \param[in] bErrors If TRUE, histogram errors are written.
1126 * Works as gmx_histogram_write_array(), but writes the cumulative
1128 * The first column in output will be the right edges of the bins.
1131 * Error output is not currently supported (zeros will be written if asked).
1133 * \see gmx_histogram_write_array()
1136 gmx_histogram_write_cum_array(FILE *fp, int n, gmx_histogram_t *h[],
1137 gmx_bool bValue, gmx_bool bErrors)
1142 prepare_output(n, h, &nbins);
1145 fprintf(fp, "%10g", h[0]->start);
1146 for (j = 0; j < n; ++j)
1150 fprintf(fp, " %10g", 0.0);
1154 fprintf(fp, " %10g", 0.0);
1158 for (i = 0; i < nbins; ++i)
1160 fprintf(fp, "%10g", h[0]->start + (i+1)*h[0]->binwidth);
1161 for (j = 0; j < n; ++j)
1163 sum[j] += i < h[j]->nbins ? h[j]->hist[i] : 0.0;
1166 fprintf(fp, " %10g", sum[j]);
1168 /* TODO: Error output not implemented */
1171 fprintf(fp, " %10g", 0.0);