src/gmxlib/statistics/histogram.c

   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 /*! \page histograms Histogram calculation
  32  *
  33  * Functions to calculate histograms are defined in histogram.h.
  34  * The principle is simple: histogram calculation is set up with
  35  * gmx_histogram_create() and gmx_histogram_set_*(), the histogram is sampled
  36  * using the provided functions, and gmx_histogram_finish() is called to
  37  * finish the sampling.
  38  * Various post-processing functions can then be used to normalize the
  39  * histogram in various ways and to write it into a file.
  40  * gmx_histogram_get_*() can be used to access data from the histogram.
  41  * gmx_histogram_free() can be used to free the memory allocated for a
  42  * histogram.
  43  *
  44  *
  45  * \section histogram_sampling Initialization and sampling
  46  *
  47  * A histogram calculation is initialized by calling gmx_histogram_create()
  48  * or gmx_histogram_create_range().
  49  * The first function takes the type of the histogram (see \ref e_histogram_t)
  50  * and the number of bins, and allocates the necessary memory.
  51  * The bin locations can be initialized with gmx_histogram_set_integerbins()
  52  * and one of gmx_histogram_set_binwidth() and gmx_histogram_set_range().
  53  * Only uniformly spaced bins are currently supported.
  54  * The second initialization function takes the beginning and end of the range
  55  * and the bin width.
  56  * The treatment of values that fall outside the histogram range can be
  57  * set with gmx_histogram_set_all().
  58  * gmx_histogram_set_blocksize() can be used to specify a block size for error
  59  * estimates. If not called, no error estimates are calculated.
  60  * Otherwise, a histogram is calculated for each block of given number of
  61  * frames, and the error estimated as the error of the mean of these block
  62  * histograms.
  63  * If the number of frames is not divisible by the block size,
  64  * the last block is discarded.
  65  * gmx_histogram_set_block_output() can be used to write the individual
  66  * block histograms to a file (it is currently not very flexible).
  67  *
  68  * After initialization, three sets of sampling functions are provided
  69  * to sample the differen types of histograms:
  70  *  - Use gmx_histogram_increment() or gmx_histogram_increment_bin() to
  71  *    calculate simple histograms (\ref HIST_SIMPLE histograms).
  72  *  - Use gmx_histogram_add() or gmx_histogram_add_to_bin() to calculate
  73  *    histograms where different points can have different weight
  74  *    (\ref HIST_WEIGHT histograms).
  75  *  - Use gmx_histogram_add_item() or gmx_histogram_add_item_to_bin() to
  76  *    calculate averages within each bin
  77  *    (\ref HIST_BINAVER histograms).
  78  *
  79  * The functions aboge that take bin indices directly do no range checking,
  80  * i.e., one must ensure that the bin number is between 0 and the number of
  81  * bins in the histogram. In contrast, the functions that use real values to
  82  * locate the correct bin use range checking, and silently discard values
  83  * outside the histogram (unless gmx_histogram_set_all() has been called).
  84  * gmx_histogram_find_bin() can be used to find the bin number corresponding
  85  * to a value, adhering to the gmx_histogram_set_all() setting.
  86  *
  87  * After a frame is sampled, gmx_histogram_finish_frame() needs to be
  88  * called.
  89  * After all data is sampled, gmx_histogram_finish() should be called.
  90  * After these calls, the histogram will be normalized such that the
  91  * sum over all bins equals the average number of samples per frame, and
  92  * the error field contains the error of the mean at each bin.
  93  *
  94  *
  95  * \section histogram_processing Post-processing
  96  *
  97  * The following functions can be used to process the histograms:
  98  *  - gmx_histogram_clone() makes a deep copy.
  99  *  - gmx_histogram_resample_dblbw() also makes a deep copy,
 100  *    but uses a binwidth that is double of the input one
 101  *    (see below for usage).
 102  *  - gmx_histogram_normalize_prob() normalizes a histogram such that its
 103  *    integral becomes one.
 104  *  - gmx_histogram_scale() and gmx_histogram_scale_vec() can be used
 105  *    to scale a histogram uniformly or non-uniformly, respectively.
 106  *
 107  * There are also functions to write histograms to a file:
 108  *  - gmx_histogram_write() writes a single histogram, optionally with
 109  *    errors.
 110  *  - gmx_histogram_write_array() writes multiple histograms with identical
 111  *    bins. Values and/or errors can be written.
 112  *  - gmx_histogram_write_cum_array() writes cumulative histograms for
 113  *    multiple histograms with identical bins.
 114  *
 115  * gmx_histogram_resample_dblbw() is useful if a histogram and its
 116  * cumulative sum are required at same points.
 117  * One can then sample the histogram with half the desired binwidth, and
 118  * then use gmx_histogram_resample_dblbw() to construct two different
 119  * versions of the histogram (by changing \p bIntegerBins), one suitable for
 120  * writing out the histogram and the other for writing the cumulative
 121  * distribution.
 122  *
 123  * Access to the histogram data is also provided through the following
 124  * functions:
 125  *  - gmx_histogram_get_nbins()
 126  *  - gmx_histogram_get_binwidth()
 127  *  - gmx_histogram_get_value()
 128  *  - gmx_histogram_get_bin_value()
 129  *  - gmx_histogram_get_values()
 130  *  - gmx_histogram_get_errors()
 131  *
 132  * gmx_histogram_free() can be used to free memory associated with a
 133  * histogram when it is no longer needed.
 134  */
 135 /*! \internal \file
 136  * \brief Implementation of functions in histogram.h.
 137  */
 138 #ifdef HAVE_CONFIG_H
 139 #include <config.h>
 140 #endif
 141
 142 #include <math.h>
 143 #include <string.h>
 144
 145 #include <smalloc.h>
 146 #include <vec.h>
 147
 148 #include <histogram.h>
 149
 150 /*! \internal \brief
 151  * Stores data for a histogram.
 152  */
 153 struct gmx_histogram_t
 154 {
 155     /** The left edge of the first bin. */
 156     real                 start;
 157     /** The right edge of the last bin. */
 158     real                 end;
 159     /** Bin width. */
 160     real                 binwidth;
 161     /** Number of bins. */
 162     int                  nbins;
 163     /** The final histogram. */
 164     double              *hist;
 165     /** Standard deviation of each bin (calculated from block histograms). */
 166     double              *histerr;
 167
 168     /** Histogram type. */
 169     e_histogram_t        type;
 170     /** Histogram flags; */
 171     int                  flags;
 172     /** Block size for averaging. */
 173     int                  bsize;
 174     /** Output file for block histograms (can be NULL). */
 175     FILE                *blockfp;
 176
 177     /** Inverse bin width. */
 178     real                 invbw;
 179     /** Current histogram value. */
 180     double              *chist;
 181     /** Current histogram count. */
 182     int                 *cn;
 183     /** Number of frames read for the current block so far. */
 184     int                  nframes;
 185     /** Number of blocks averaged so far. */
 186     int                  nblocks;
 187 };
 188
 189 /*!
 190  * \param[out] hp      Histogram to create.
 191  * \param[in]  type    Histogram type to create;
 192  * \param[in]  nbins   Number of bins.
 193  * \returns    0 on success, a non-zero error code on error.
 194  *
 195  * Initialized the histogram structure \p h with the provided values,
 196  * allocating the necessary memory.
 197  */
 198 int
 199 gmx_histogram_create(gmx_histogram_t **hp, e_histogram_t type, int nbins)
 200 {
 201     gmx_histogram_t *h;
 202
 203     if (nbins <= 0)
 204     {
 205         *hp = NULL;
 206         gmx_incons("number of histogram bins <= 0");
 207         return EINVAL;
 208     }
 209
 210     snew(h, 1);
 211     h->start    = 0;
 212     h->end      = 0;
 213     h->binwidth = 0;
 214     h->nbins    = nbins;
 215     h->hist     = NULL;
 216     h->histerr  = NULL;
 217
 218     h->type     = type;
 219     h->flags    = 0;
 220
 221     h->bsize    = 0;
 222     h->blockfp  = NULL;
 223     h->invbw    = 0;
 224     h->chist    = NULL;
 225     h->cn       = NULL;
 226     h->nframes  = 0;
 227     h->nblocks  = 0;
 228
 229     snew(h->hist,    nbins + 1);
 230     snew(h->histerr, nbins + 1);
 231     snew(h->cn,      nbins + 1);
 232     if (type != HIST_SIMPLE)
 233     {
 234         snew(h->chist, nbins + 1);
 235     }
 236     gmx_histogram_clear(h);
 237
 238     *hp = h;
 239     return 0;
 240 }
 241
 242 /*!
 243  * \param[out] hp      Histogram to create.
 244  * \param[in]  type    Histogram type to create;
 245  * \param[in]  start   Left edge/center of the first bin
 246  *   (depending on \p bIntegerBins).
 247  * \param[in]  end     Right edge/center of the last bin
 248  *   (depending on \p bIntegerBins).
 249  * \param[in]  binw    Bin width.
 250  * \param[in]  bIntegerBins If true, histogram bins are centered at
 251  *   \c start+n*binwidth, otherwise the centers are at
 252  *   \c start+(n+0.5)*binwidth.
 253  * \returns    0 on success, a non-zero error code on error.
 254  *
 255  * Initialized the histogram structure \p h with the provided values,
 256  * allocating the necessary memory.
 257  */
 258 int
 259 gmx_histogram_create_range(gmx_histogram_t **hp, e_histogram_t type,
 260                            real start, real end, real binw, gmx_bool bIntegerBins)
 261 {
 262     gmx_histogram_t *h;
 263     int              nbins;
 264     int              rc;
 265
 266     *hp = NULL;
 267     if (start >= end || binw <= 0)
 268     {
 269         gmx_incons("histogram left edge larger than right edge or bin width <= 0");
 270         return EINVAL;
 271     }
 272
 273     /* Adjust the end points and calculate the number of bins */
 274     if (bIntegerBins)
 275     {
 276         nbins = ceil((end - start) / binw) + 1;
 277         end   = start + (nbins - 1) * binw;
 278     }
 279     else
 280     {
 281         start = binw * floor(start / binw);
 282         end   = binw * ceil(end / binw);
 283         if (start != 0)
 284         {
 285             start -= binw;
 286         }
 287         end += binw;
 288         nbins = (int)((end - start) / binw + 0.5);
 289     }
 290     /* Create the histogram */
 291     rc = gmx_histogram_create(&h, type, nbins);
 292     if (rc != 0)
 293     {
 294         return rc;
 295     }
 296     /* Set it up */
 297     gmx_histogram_set_integerbins(h, bIntegerBins);
 298     gmx_histogram_set_binwidth(h, start, binw);
 299
 300     *hp = h;
 301     return 0;
 302 }
 303
 304 /*!
 305  * \param[in,out] h       Histogram to clear.
 306  *
 307  * Histograms are automatically cleared when initialized; you only need to
 308  * call this function if you want to reuse a histogram structure that has
 309  * already been used for sampling.
 310  */
 311 void
 312 gmx_histogram_clear(gmx_histogram_t *h)
 313 {
 314     int              i;
 315
 316     if (h->nbins <= 0)
 317     {
 318         return;
 319     }
 320     for (i = 0; i <= h->nbins; ++i)
 321     {
 322         h->hist[i]    = 0;
 323         h->histerr[i] = 0;
 324         if (h->chist)
 325         {
 326             h->chist[i] = 0;
 327         }
 328         h->cn[i]      = 0;
 329     }
 330     h->nframes      = 0;
 331     h->nblocks      = 0;
 332 }
 333
 334 /*!
 335  * \param[in] h  Histogram to free.
 336  *
 337  * The pointer \p h is invalid after the call.
 338  */
 339 void
 340 gmx_histogram_free(gmx_histogram_t *h)
 341 {
 342     sfree(h->chist);
 343     sfree(h->cn);
 344     sfree(h->hist);
 345     sfree(h->histerr);
 346     sfree(h);
 347 }
 348
 349 /*!
 350  * \param[in,out] h      Histogram data structure.
 351  * \param[in]     start  Left edge/center of the first bin
 352  *   (depending on gmx_histogram_set_integerbins()).
 353  * \param[in]     binw   Bin width.
 354  * \returns       0 on success, a non-zero error code on error.
 355  */
 356 int
 357 gmx_histogram_set_binwidth(gmx_histogram_t *h, real start, real binw)
 358 {
 359     if (binw <= 0)
 360     {
 361         gmx_incons("histogram binwidth <= 0");
 362         return EINVAL;
 363     }
 364     if (h->flags & HIST_INTEGERBINS)
 365     {
 366         start   -= 0.5*binw;
 367     }
 368     h->start     = start;
 369     h->binwidth  = binw;
 370     h->end       = start + h->nbins * binw;
 371     h->invbw     = 1.0/binw;
 372     h->flags    |= HIST_INITBW;
 373     return 0;
 374 }
 375
 376 /*!
 377  * \param[in,out] h      Histogram data structure.
 378  * \param[in]     start  Left edge/center of the first bin
 379  *   (depending on gmx_histogram_set_integerbins()).
 380  * \param[in]     end    Right edge/center of the last bin
 381  *   (depending on gmx_histogram_set_integerbins()).
 382  * \returns       0 on success, a non-zero error code on error.
 383  */
 384 int
 385 gmx_histogram_set_range(gmx_histogram_t *h, real start, real end)
 386 {
 387     if (start >= end)
 388     {
 389         gmx_incons("histogram left edge larger than right edge");
 390         return EINVAL;
 391     }
 392     h->start         = start;
 393     h->end           = end;
 394     if (h->flags & HIST_INTEGERBINS)
 395     {
 396         h->binwidth  = (end - start) / (h->nbins - 1);
 397         start       -= 0.5*h->binwidth;
 398         end         += 0.5*h->binwidth;
 399     }
 400     else
 401     {
 402         h->binwidth  = (end - start) / h->nbins;
 403     }
 404     h->invbw         = 1.0/h->binwidth;
 405     h->flags        &= ~HIST_INITBW;
 406     return 0;
 407 }
 408
 409 /*!
 410  * \param[in,out] h      Histogram data structure.
 411  * \param[in]     bIntegerBins If true, histogram bins are centered at
 412  *   \c start+n*binwidth, otherwise the centers are at
 413  *   \c start+(n+0.5)*binwidth.
 414  */
 415 void
 416 gmx_histogram_set_integerbins(gmx_histogram_t *h, gmx_bool bIntegerBins)
 417 {
 418     /* Adjust the ranges if they have been initialized */
 419     if (h->start < h->end)
 420     {
 421         if (h->flags & HIST_INTEGERBINS)
 422         {
 423             h->start   += 0.5*h->binwidth;
 424             if (h->flags & HIST_INITBW)
 425             {
 426                 h->end += 0.5*h->binwidth;
 427             }
 428             else
 429             {
 430                 h->end -= 0.5*h->binwidth;
 431             }
 432         }
 433         if (bIntegerBins)
 434         {
 435             h->start   -= 0.5*h->binwidth;
 436             if (h->flags & HIST_INITBW)
 437             {
 438                 h->end -= 0.5*h->binwidth;
 439             }
 440             else
 441             {
 442                 h->end += 0.5*h->binwidth;
 443             }
 444         }
 445     }
 446     if (bIntegerBins)
 447     {
 448         h->flags |= HIST_INTEGERBINS;
 449     }
 450     else
 451     {
 452         h->flags &= ~HIST_INTEGERBINS;
 453     }
 454 }
 455
 456 /*!
 457  * \param[in,out] h     Histogram data structure.
 458  * \param[in]     bAll  If true, values outside the histogram edges are added
 459  *   to the bins at the edges.
 460  *
 461  * \p bAll can be used to avoid rounding errors in cases where the histogram
 462  * spans the full range of possible values. If set, the values that are at
 463  * the exact maximum are still correctly included.
 464  */
 465 void
 466 gmx_histogram_set_all(gmx_histogram_t *h, gmx_bool bAll)
 467 {
 468     if (bAll)
 469     {
 470         h->flags |= HIST_ALL;
 471     }
 472     else
 473     {
 474         h->flags &= ~HIST_ALL;
 475     }
 476 }
 477
 478 /*!
 479  * \param[in,out] h      Histogram data structure.
 480  * \param[in]     bsize  Block size for error estimates.
 481  * \returns       0 on success, a non-zero error code on error.
 482  *
 483  * If \p bsize is zero, no block averaging or error estimates are done.
 484  * This is also the case if this function is not called.
 485  */
 486 int
 487 gmx_histogram_set_blocksize(gmx_histogram_t *h, int bsize)
 488 {
 489     if (bsize < 0)
 490     {
 491         gmx_incons("histogram block size < 0");
 492         return EINVAL;
 493     }
 494     h->bsize = bsize;
 495     return 0;
 496 }
 497
 498 /*!
 499  * \param[in,out] h  Histogram data structure.
 500  * \param[in]     fp File for block output.
 501  * \returns       0 on success, a non-zero error code on error.
 502  *
 503  * Sets a file into which each block histogram (the histogram after each
 504  * \c gmx_histogram_t::bsize frames) is written.
 505  * All histograms are written to the same file, separated by empty lines.
 506  * If this function is not called, the block histograms are only used for
 507  * error estimation, not written to a file.
 508  */
 509 int
 510 gmx_histogram_set_block_output(gmx_histogram_t *h, FILE *fp)
 511 {
 512     if (h->bsize <= 0)
 513     {
 514         gmx_incons("histogram block size not set but output initialized");
 515         return EINVAL;
 516     }
 517     h->blockfp = fp;
 518     return 0;
 519 }
 520
 521 /*!
 522  * \param[in] h     Histogram data.
 523  * \param[in] pos   Position.
 524  * \returns   Bin index in \p h corresponding to \p pos, or -1 if there is no
 525  *    such bin.
 526  *
 527  * If gmx_histogram_set_all() has been called with TRUE, values outside the
 528  * histogram are mapped to the bins at the edges.
 529  */
 530 int
 531 gmx_histogram_find_bin(gmx_histogram_t *h, real pos)
 532 {
 533     if (pos < h->start)
 534     {
 535         if (h->flags & HIST_ALL)
 536         {
 537             return 0;
 538         }
 539         else
 540         {
 541             return -1;
 542         }
 543     }
 544     else if (pos >= h->end)
 545     {
 546         if (h->flags & HIST_ALL)
 547         {
 548             return h->nbins - 1;
 549         }
 550         else
 551         {
 552             return -1;
 553         }
 554     }
 555     /* Calculate the bin index if we are inside the box */
 556     return (int)((pos - h->start)*h->invbw);
 557 }
 558
 559 /*!
 560  * \param[in] h     Histogram data.
 561  * \returns   Number of bins in \p h.
 562  */
 563 int
 564 gmx_histogram_get_nbins(gmx_histogram_t *h)
 565 {
 566     return h->nbins;
 567 }
 568
 569 /*!
 570  * \param[in] h     Histogram data.
 571  * \returns   Bin width of \p h, 0 if no binwidth has been set.
 572  */
 573 real
 574 gmx_histogram_get_binwidth(gmx_histogram_t *h)
 575 {
 576     return h->binwidth;
 577 }
 578
 579 /*!
 580  * \param[in]  h     Histogram data.
 581  * \param[in]  pos   Position.
 582  * \param[out] value Pointer to receive the value (can be NULL).
 583  * \param[out] err   Pointer to receive the value (can be NULL).
 584  *
 585  * If \p pos is outside the range of the histogram, zeros are returned.
 586  */
 587 void
 588 gmx_histogram_get_value(gmx_histogram_t *h, real pos, double *value, double *err)
 589 {
 590     int     bin;
 591     double  v, e;
 592
 593     if (pos < h->start || pos > h->end)
 594     {
 595         v = e = 0;
 596     }
 597     else
 598     {
 599         bin = gmx_histogram_find_bin(h, pos);
 600         if (bin < 0)
 601         {
 602             v = e = 0;
 603         }
 604         else
 605         {
 606             v = h->hist[bin];
 607             e = h->histerr[bin];
 608         }
 609     }
 610     if (value)
 611     {
 612         *value = v;
 613     }
 614     if (err)
 615     {
 616         *err = e;
 617     }
 618 }
 619
 620 /*!
 621  * \param[in]  h     Histogram data.
 622  * \param[in]  bin   Bin number.
 623  * \param[out] value Pointer to receive the value (can be NULL).
 624  * \param[out] err   Pointer to receive the value (can be NULL).
 625  *
 626  * If \p bin is outside the valid range, zeros are returned.
 627  */
 628 void
 629 gmx_histogram_get_bin_value(gmx_histogram_t *h, int bin, double *value, double *err)
 630 {
 631     double  v, e;
 632
 633     if (bin < 0 || bin >= h->nbins)
 634     {
 635         v = e = 0;
 636     }
 637     else
 638     {
 639         v = h->hist[bin];
 640         e = h->histerr[bin];
 641     }
 642     if (value)
 643     {
 644         *value = v;
 645     }
 646     if (err)
 647     {
 648         *err = e;
 649     }
 650 }
 651
 652 /*!
 653  * \param[in]  h     Histogram data.
 654  * \returns    Pointer to an array of values for \p h.
 655  *
 656  * The returned array has one element for each bin of \p h.
 657  * The returned pointer should not be freed.
 658  */
 659 double *
 660 gmx_histogram_get_values(gmx_histogram_t *h)
 661 {
 662     return h->hist;
 663 }
 664
 665 /*!
 666  * \param[in]  h     Histogram data.
 667  * \returns    Pointer to an array of errors for \p h.
 668  *
 669  * The returned array has one element for each bin of \p h.
 670  * The returned pointer should not be freed.
 671  */
 672 double *
 673 gmx_histogram_get_errors(gmx_histogram_t *h)
 674 {
 675     return h->histerr;
 676 }
 677
 678 /*!
 679  * \param[in,out] h     Histogram data.
 680  * \param[in]     pos   Position.
 681  */
 682 void
 683 gmx_histogram_increment(gmx_histogram_t *h, real pos)
 684 {
 685     int bin = gmx_histogram_find_bin(h, pos);
 686     if (bin < 0)
 687     {
 688         return;
 689     }
 690     h->cn[bin]++;
 691 }
 692
 693 /*!
 694  * \param[in,out] h     Histogram data.
 695  * \param[in]     pos   Position.
 696  * \param[in]     value Value to add.
 697  */
 698 void
 699 gmx_histogram_add(gmx_histogram_t *h, real pos, double value)
 700 {
 701     int bin = gmx_histogram_find_bin(h, pos);
 702     if (bin < 0)
 703     {
 704         return;
 705     }
 706     h->chist[bin] += value;
 707 }
 708
 709 /*!
 710  * \param[in,out] h     Histogram data.
 711  * \param[in]     pos   Position.
 712  * \param[in]     value Value to add.
 713  */
 714 void
 715 gmx_histogram_add_item(gmx_histogram_t *h, real pos, double value)
 716 {
 717     int bin = gmx_histogram_find_bin(h, pos);
 718     if (bin < 0)
 719     {
 720         return;
 721     }
 722     h->chist[bin] += value;
 723     h->cn[bin]++;
 724 }
 725
 726 /*!
 727  * \param[in,out] h     Histogram data.
 728  * \param[in]     bin   Bin number.
 729  *
 730  * No checks for out-of-bound errors are done: \p bin should be >=0 and
 731  * < \p h->nbins.
 732  */
 733 void
 734 gmx_histogram_increment_bin(gmx_histogram_t *h, int bin)
 735 {
 736     h->cn[bin]++;
 737 }
 738
 739 /*!
 740  * \param[in,out] h     Histogram data.
 741  * \param[in]     bin   Bin number.
 742  * \param[in]     value Value to add.
 743  *
 744  * No checks for out-of-bound errors are done: \p bin should be >=0 and
 745  * < \p h->nbins.
 746  */
 747 void
 748 gmx_histogram_add_to_bin(gmx_histogram_t *h, int bin, double value)
 749 {
 750     h->chist[bin] += value;
 751 }
 752
 753 /*!
 754  * \param[in,out] h     Histogram data.
 755  * \param[in]     bin   Bin number.
 756  * \param[in]     value Value to add.
 757  *
 758  * No checks for out-of-bound errors are done: \p bin should be >=0 and
 759  * < \p h->nbins.
 760  */
 761 void
 762 gmx_histogram_add_item_to_bin(gmx_histogram_t *h, int bin, double value)
 763 {
 764     h->chist[bin] += value;
 765     h->cn[bin]++;
 766 }
 767
 768 /*! \brief
 769  * Processes a sampled block histogram.
 770  *
 771  * \param[in,out] h    Histogram data structure.
 772  */
 773 static void
 774 finish_histogram_block(gmx_histogram_t *h)
 775 {
 776     int   i;
 777     real  v;
 778
 779     if (h->nframes == 0)
 780     {
 781         return;
 782     }
 783
 784     if (h->flags & HIST_ALL)
 785     {
 786         if (h->chist)
 787         {
 788             h->chist[h->nbins-1] += h->chist[h->nbins];
 789         }
 790         h->cn[h->nbins-1] += h->cn[h->nbins];
 791     }
 792
 793     for (i = 0; i <= h->nbins; ++i)
 794     {
 795         if (h->chist)
 796         {
 797             v = h->chist[i] / (h->cn[i] > 0 ? h->cn[i] : h->nframes);
 798         }
 799         else
 800         {
 801             v = ((real)h->cn[i]) / h->nframes;
 802         }
 803         if (h->blockfp)
 804         {
 805             fprintf(h->blockfp, "%10g %10g\n", h->start + (i+0.5)*h->binwidth, v);
 806         }
 807         h->hist[i]    += v;
 808         h->histerr[i] += sqr(v);
 809         if (h->chist)
 810         {
 811             h->chist[i] = 0;
 812         }
 813         h->cn[i] = 0;
 814     }
 815     if (h->blockfp)
 816     {
 817         fprintf(h->blockfp, "\n");
 818     }
 819     h->nblocks++;
 820     h->nframes = 0;
 821 }
 822
 823 /*!
 824  * \param[in,out] h    Histogram data structure.
 825  *
 826  * Should be called after each frame of data.
 827  *
 828  * \see gmx_histogram_finish()
 829  */
 830 void
 831 gmx_histogram_finish_frame(gmx_histogram_t *h)
 832 {
 833     h->nframes++;
 834     if (h->nframes == h->bsize)
 835     {
 836         finish_histogram_block(h);
 837     }
 838 }
 839
 840 /*!
 841  * \param[in,out] h    Histogram data structure.
 842  *
 843  * Normalizes the histogram.
 844  * Should be called after all frames have been sampled and before any
 845  * post-processing of the histogram.
 846  * If block size has not been set with gmx_histogram_set_blocksize(),
 847  * this function does no normalization, but it should still be called,
 848  * otherwise the \c gmx_histogram_t::hist array will contain only zeros.
 849  *
 850  * \see gmx_histogram_finish_frame()
 851  */
 852 void
 853 gmx_histogram_finish(gmx_histogram_t *h)
 854 {
 855     int  i;
 856
 857     if (h->nframes > 0 || h->bsize == 0)
 858     {
 859         if (h->nframes < h->bsize)
 860         {
 861             fprintf(stderr, "Last block smaller (%d frames) than the target size (%d frames) skipped \n",
 862                     h->nframes, h->bsize);
 863         }
 864         else
 865         {
 866             finish_histogram_block(h);
 867         }
 868     }
 869     if (h->nblocks == 0)
 870     {
 871         return;
 872     }
 873
 874     for (i = 0; i <= h->nbins; ++i)
 875     {
 876         h->hist[i]    /= h->nblocks;
 877         h->histerr[i] /= h->nblocks;
 878         h->histerr[i]  = sqrt((h->histerr[i] - sqr(h->hist[i])) / h->nblocks);
 879     }
 880 }
 881
 882 /*!
 883  * \param[out] destp        Destination histogram.
 884  * \param[in]  src          Source histogram.
 885  * \param[in]  bIntegerBins Control bin center position.
 886  *
 887  * Uses \p src to create a new histogram that has double the binwidth.
 888  * If \p bIntegerBins is TRUE, the first new bin will be centered at
 889  * \c src->start, otherwise the left edge will be at \c src->start.
 890  *
 891  * This function is mostly useful if one needs to sample both the
 892  * histogram and the cumulative histogram at same points.
 893  * To achieve this, first sample a histogram with half the desired
 894  * binwidth, and then use this function to create two different verions
 895  * of it with different values of \p bIntegerBins.
 896  * gmx_histogram_write_array() and gmx_histogram_write_cum_array() can
 897  * now used to write out the values correctly at identical values.
 898  *
 899  * Should be called only after gmx_histogram_finish() has been called.
 900  * \p src should have been sampled without gmx_histogram_set_integerbins().
 901  */
 902 void
 903 gmx_histogram_resample_dblbw(gmx_histogram_t **destp, gmx_histogram_t *src,
 904                              gmx_bool bIntegerBins)
 905 {
 906     gmx_histogram_t *dest;
 907     int              i, j;
 908     real             v, ve;
 909
 910     gmx_histogram_create(destp, src->type, src->nbins / 2);
 911     dest = *destp;
 912     sfree(dest->chist); dest->chist = NULL;
 913     sfree(dest->cn);    dest->cn    = NULL;
 914     gmx_histogram_set_binwidth(dest, src->start, src->binwidth * 2);
 915     gmx_histogram_set_integerbins(dest, bIntegerBins);
 916
 917     for (i = j = 0; i < dest->nbins; ++i)
 918     {
 919         if (bIntegerBins && i == 0)
 920         {
 921             v  = src->hist[0];
 922             ve = sqr(src->histerr[0]);
 923             ++j;
 924         }
 925         else
 926         {
 927             v  = src->hist[j]         + src->hist[j+1];
 928             ve = sqr(src->histerr[j]) + sqr(src->histerr[j+1]);
 929             j += 2;
 930         }
 931         ve = sqrt(ve);
 932         dest->hist[i]    = v;
 933         dest->histerr[i] = ve;
 934     }
 935     dest->hist[dest->nbins]    = 0;
 936     dest->histerr[dest->nbins] = 0;
 937 }
 938
 939 /*!
 940  * \param[out] destp Destination histogram.
 941  * \param[in]  src   Source histogram.
 942  *
 943  * Makes a clone of \p src for post-processing.
 944  */
 945 void
 946 gmx_histogram_clone(gmx_histogram_t **destp, gmx_histogram_t *src)
 947 {
 948     gmx_histogram_t *dest;
 949
 950     snew(dest, 1);
 951     memcpy(dest, src, sizeof(*dest));
 952
 953     /* These are not needed in post-processing */
 954     dest->blockfp = NULL;
 955     dest->chist   = NULL;
 956     dest->cn      = NULL;
 957
 958     /* Make a deep copy of the actual histograms */
 959     snew(dest->hist,    src->nbins+1);
 960     snew(dest->histerr, src->nbins+1);
 961     memcpy(dest->hist,    src->hist,    (src->nbins+1)*sizeof(real));
 962     memcpy(dest->histerr, src->histerr, (src->nbins+1)*sizeof(real));
 963
 964     *destp = dest;
 965 }
 966
 967 /*!
 968  * \param[in,out] h  Histogram to normalize.
 969  *
 970  * Normalizes the histogram such that its integral equals one.
 971  */
 972 void
 973 gmx_histogram_normalize_prob(gmx_histogram_t *h)
 974 {
 975     int  i;
 976     real sum;
 977     real normfac;
 978
 979     sum = 0;
 980     for (i = 0; i <= h->nbins; ++i)
 981     {
 982         sum += h->hist[i];
 983     }
 984
 985     normfac = h->invbw / sum;
 986     gmx_histogram_scale(h, normfac);
 987 }
 988
 989 /*!
 990  * \param[in,out] h    Histogram to normalize.
 991  * \param[in]     norm Scaling factor.
 992  *
 993  * All bin values are multiplied by \p norm.
 994  */
 995 void
 996 gmx_histogram_scale(gmx_histogram_t *h, real norm)
 997 {
 998     int  i;
 999
1000     for (i = 0; i <= h->nbins; ++i)
1001     {
1002         h->hist[i]    *= norm;
1003         h->histerr[i] *= norm;
1004     }
1005 }
1006
1007 /*!
1008  * \param[in,out] h    Histogram to normalize.
1009  * \param[in]     norm Scaling vector.
1010  *
1011  * The i'th bin is multiplied by \p norm[i].
1012  */
1013 void
1014 gmx_histogram_scale_vec(gmx_histogram_t *h, real norm[])
1015 {
1016     int  i;
1017
1018     for (i = 0; i < h->nbins; ++i)
1019     {
1020         h->hist[i]    *= norm[i];
1021         h->histerr[i] *= norm[i];
1022     }
1023     h->hist[h->nbins]    *= norm[h->nbins-1];
1024     h->histerr[h->nbins] *= norm[h->nbins-1];
1025 }
1026
1027 /*! \brief
1028  * Makes some checks on output histograms and finds the maximum number of bins.
1029  *
1030  * \param[in]  n     Number of histograms in \p h.
1031  * \param[in]  h     Array of histograms.
1032  * \param[out] nbins Pointer to a value that will receive the maximum number
1033  *   of bins in \p h.
1034  */
1035 static void
1036 prepare_output(int n, gmx_histogram_t *h[], int *nbins)
1037 {
1038     int  j;
1039
1040     *nbins = 0;
1041     for (j = 0; j < n; ++j)
1042     {
1043         if (!gmx_within_tol(h[j]->start, h[0]->start, GMX_REAL_EPS))
1044         {
1045             fprintf(stderr, "gmx_ana_histogram_write: histogram start values not identical\n");
1046         }
1047         if (!gmx_within_tol(h[j]->binwidth, h[0]->binwidth, GMX_REAL_EPS))
1048         {
1049             fprintf(stderr, "gmx_ana_histogram_write: bin widths not identical\n");
1050         }
1051         if (*nbins < h[j]->nbins)
1052         {
1053             *nbins = h[j]->nbins;
1054         }
1055     }
1056 }
1057
1058 /*!
1059  * \param[in] fp           Output file.
1060  * \param[in] h            Array of histograms to write.
1061  * \param[in] bErrors      If TRUE, histogram errors are written.
1062  *
1063  * Convenience wrapper for gmx_histogram_write_array() for writing a single
1064  * histogram.
1065  *
1066  * \see gmx_histogram_write_array()
1067  */
1068 void
1069 gmx_histogram_write(FILE *fp, gmx_histogram_t *h, gmx_bool bErrors)
1070 {
1071     gmx_histogram_write_array(fp, 1, &h, TRUE, bErrors);
1072 }
1073
1074 /*!
1075  * \param[in] fp           Output file.
1076  * \param[in] n            Number of histograms in \p h.
1077  * \param[in] h            Array of histograms to write.
1078  * \param[in] bValue       If TRUE, histogram values are written.
1079  * \param[in] bErrors      If TRUE, histogram errors are written.
1080  *
1081  * All the histograms in the array \p h should have the same bin widths and
1082  * left edges, otherwise the behavior is undefined.
1083  * The output format is one bin per line. The first number on the line is the
1084  * center of the bin, followed by one or two numbers for each histogram
1085  * (depending on \p bValue and \p bErrors).
1086  * If both \p bValue and \p bErrors are both TRUE, the values are written
1087  * before the errors.
1088  */
1089 void
1090 gmx_histogram_write_array(FILE *fp, int n, gmx_histogram_t *h[],
1091                           gmx_bool bValue, gmx_bool bErrors)
1092 {
1093     int           i, j, nbins;
1094
1095     prepare_output(n, h, &nbins);
1096     for (i = 0; i < nbins; ++i)
1097     {
1098         fprintf(fp, "%10g", h[0]->start + (i+0.5)*h[0]->binwidth);
1099         for (j = 0; j < n; ++j)
1100         {
1101             if (bValue)
1102             {
1103                 fprintf(fp, " %10g", i < h[j]->nbins ? h[j]->hist[i] : 0.0);
1104             }
1105             if (bErrors)
1106             {
1107                 fprintf(fp, " %10g", i < h[j]->nbins ? h[j]->histerr[i] : 0.0);
1108             }
1109         }
1110         fprintf(fp, "\n");
1111     }
1112     fprintf(fp, "\n");
1113 }
1114
1115 /*!
1116  * \param[in] fp           Output file.
1117  * \param[in] n            Number of histograms in \p h.
1118  * \param[in] h            Array of histograms to write.
1119  * \param[in] bValue       If TRUE, histogram values are written.
1120  * \param[in] bErrors      If TRUE, histogram errors are written.
1121  *
1122  * Works as gmx_histogram_write_array(), but writes the cumulative
1123  * histograms.
1124  * The first column in output will be the right edges of the bins.
1125  *
1126  * \note
1127  * Error output is not currently supported (zeros will be written if asked).
1128  *
1129  * \see gmx_histogram_write_array()
1130  */
1131 void
1132 gmx_histogram_write_cum_array(FILE *fp, int n, gmx_histogram_t *h[],
1133                                   gmx_bool bValue, gmx_bool bErrors)
1134 {
1135     int           i, j, nbins;
1136     double       *sum;
1137
1138     prepare_output(n, h, &nbins);
1139     snew(sum, n);
1140
1141     fprintf(fp, "%10g", h[0]->start);
1142     for (j = 0; j < n; ++j)
1143     {
1144         if (bValue)
1145         {
1146             fprintf(fp, " %10g", 0.0);
1147         }
1148         if (bErrors)
1149         {
1150             fprintf(fp, " %10g", 0.0);
1151         }
1152     }
1153     fprintf(fp, "\n");
1154     for (i = 0; i < nbins; ++i)
1155     {
1156         fprintf(fp, "%10g", h[0]->start + (i+1)*h[0]->binwidth);
1157         for (j = 0; j < n; ++j)
1158         {
1159             sum[j] += i < h[j]->nbins ? h[j]->hist[i] : 0.0;
1160             if (bValue)
1161             {
1162                 fprintf(fp, " %10g", sum[j]);
1163             }
1164             /* TODO: Error output not implemented */
1165             if (bErrors)
1166             {
1167                 fprintf(fp, " %10g", 0.0);
1168             }
1169         }
1170         fprintf(fp, "\n");
1171     }
1172     fprintf(fp, "\n");
1173
1174     sfree(sum);
1175 }