src/gromacs/random/gammadistribution.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2015,2016,2019,2021, 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
  36 /*! \file
  37  * \brief The gamma distribution
  38  *
  39  * Portable version of the gamma distribution that generates the same sequence
  40  * on all platforms.
  41  *
  42  * \note The gamma distribution is broken in some standard library headers
  43  * (including those shipped with gcc-4.9), and it is not guaranteed to
  44  * generate the same result on stdlibc++ and libc++. Use this one instead so
  45  * our unit tests produce the same values on all platforms.
  46  *
  47  * \author Erik Lindahl <erik.lindahl@gmail.com>
  48  * \inpublicapi
  49  * \ingroup module_random
  50  */
  51
  52 #ifndef GMX_RANDOM_GAMMADISTRIBUTION_H
  53 #define GMX_RANDOM_GAMMADISTRIBUTION_H
  54
  55 #include <cmath>
  56
  57 #include <limits>
  58 #include <memory>
  59
  60 #include "gromacs/random/exponentialdistribution.h"
  61 #include "gromacs/random/uniformrealdistribution.h"
  62 #include "gromacs/utility/classhelpers.h"
  63
  64 /*
  65  * The workaround implementation for the broken std::gamma_distribution in the
  66  * gcc-4.6 headers has been modified from the LLVM libcxx headers, distributed
  67  * under the MIT license:
  68  *
  69  * Copyright (c) The LLVM compiler infrastructure
  70  *
  71  * Permission is hereby granted, free of charge, to any person obtaining a copy
  72  * of this software and associated documentation files (the "Software"), to deal
  73  * in the Software without restriction, including without limitation the rights
  74  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  75  * copies of the Software, and to permit persons to whom the Software is
  76  * furnished to do so, subject to the following conditions:
  77  *
  78  * The above copyright notice and this permission notice shall be included in
  79  * all copies or substantial portions of the Software.
  80  *
  81  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  82  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  83  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  84  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  85  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  86  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  87  * THE SOFTWARE.
  88  */
  89
  90 namespace gmx
  91 {
  92
  93 /*! \brief Gamma distribution
  94  *
  95  *  The C++ standard library does provide a gamma distribution, but when
  96  *  using libstdc++-4.4.7 with at least gcc-4.6 the headers
  97  *  produce errors. Even for newer compilers, libstdc++ and libc++ appear to
  98  *  use different algorithms to generate it, which means their values differ
  99  *  in contrast to the uniform and normal distributions where they are
 100  *  identical. To avoid both compiler bugs and make it easier to use
 101  *  GROMACS unit tests that depend on random numbers, we have our
 102  *  own implementation.
 103  *
 104  *  Be warned that the gamma distribution works like the standard
 105  *  normal distribution and keeps drawing values from the random engine
 106  *  in a loop, so you want to make sure you use a random stream with a
 107  *  very large margin to make sure you do not run out of random numbers
 108  *  in an unlucky case (which will lead to an exception with the GROMACS
 109  *  default random engine).
 110  *
 111  *  The gamma distribution is defined as
 112  *
 113  * \f[
 114  *     p(x|\alpha,\beta) = \frac{1}{\Gamma(\alpha)\beta^{alpha}} x^{\alpha - 1}
 115  * e^{-\frac{x}{\beta}}, x\geq 0 \f]
 116  *
 117  * \tparam RealType Floating-point type, real by default in GROMACS.
 118  */
 119 template<class RealType = real>
 120 class GammaDistribution
 121 {
 122 public:
 123     /*! \brief Type of values returned */
 124     typedef RealType result_type;
 125
 126     /*! \brief Gamma distribution parameters */
 127     class param_type
 128     {
 129         /*! \brief First parameter of gamma distribution */
 130         result_type alpha_;
 131         /*! \brief Second parameter of gamma distribution */
 132         result_type beta_;
 133
 134     public:
 135         /*! \brief Reference back to the distribution class */
 136         typedef GammaDistribution distribution_type;
 137
 138         /*! \brief Construct parameter block
 139          *
 140          * \param alpha  First parameter of gamma distribution
 141          * \param beta   Second parameter of gamma distribution
 142          */
 143         explicit param_type(result_type alpha = 1.0, result_type beta = 1.0) :
 144             alpha_(alpha), beta_(beta)
 145         {
 146         }
 147
 148         /*! \brief Return first parameter */
 149         result_type alpha() const { return alpha_; }
 150         /*! \brief Return second parameter */
 151         result_type beta() const { return beta_; }
 152
 153         /*! \brief True if two parameter sets will return the same distribution.
 154          *
 155          * \param x  Instance to compare with.
 156          */
 157         bool operator==(const param_type& x) const
 158         {
 159             return alpha_ == x.alpha_ && beta_ == x.beta_;
 160         }
 161
 162         /*! \brief True if two parameter sets will return different distributions
 163          *
 164          * \param x  Instance to compare with.
 165          */
 166         bool operator!=(const param_type& x) const { return !operator==(x); }
 167     };
 168
 169     /*! \brief Construct new distribution with given floating-point parameters.
 170      *
 171      * \param alpha  First parameter of gamma distribution
 172      * \param beta   Second parameter of gamma distribution
 173      */
 174     explicit GammaDistribution(result_type alpha = 1.0, result_type beta = 1.0) :
 175         param_(param_type(alpha, beta))
 176     {
 177     }
 178
 179     /*! \brief Construct new distribution from parameter class
 180      *
 181      * \param param  Parameter class as defined inside gmx::GammaDistribution.
 182      */
 183     explicit GammaDistribution(const param_type& param) : param_(param) {}
 184
 185     /*! \brief Flush all internal saved values  */
 186     void reset() {}
 187
 188     /*! \brief Return values from gamma distribution with internal parameters
 189      *
 190      *  \tparam Rng   Random engine class
 191      *
 192      *  \param  g     Random engine
 193      */
 194     template<class Rng>
 195     result_type operator()(Rng& g)
 196     {
 197         return (*this)(g, param_);
 198     }
 199
 200     /*! \brief Return value from gamma distribution with given parameters
 201      *
 202      *  \tparam Rng   Random engine class
 203      *
 204      *  \param  g     Random engine
 205      *  \param  param Parameters to use
 206      */
 207     template<class Rng>
 208     result_type operator()(Rng& g, const param_type& param)
 209     {
 210         result_type                          alpha = param.alpha();
 211         UniformRealDistribution<result_type> uniformDist(0, 1);
 212         ExponentialDistribution<result_type> expDist;
 213
 214         result_type x;
 215
 216         if (alpha == 1.0)
 217         {
 218             x = expDist(g);
 219         }
 220         else if (alpha > 1.0)
 221         {
 222             const result_type b = alpha - 1.0;
 223             const result_type c = 3.0 * alpha - result_type(0.75);
 224
 225             while (true)
 226             {
 227                 const result_type u = uniformDist(g);
 228                 const result_type v = uniformDist(g);
 229                 const result_type w = u * (1 - u);
 230
 231                 if (w != 0)
 232                 {
 233                     const result_type y = std::sqrt(c / w) * (u - result_type(0.5));
 234                     x                   = b + y;
 235
 236                     if (x >= 0)
 237                     {
 238                         const result_type z = 64 * w * w * w * v * v;
 239
 240                         if (z <= 1.0 - 2.0 * y * y / x)
 241                         {
 242                             break;
 243                         }
 244                         if (std::log(z) <= 2.0 * (b * std::log(x / b) - y))
 245                         {
 246                             break;
 247                         }
 248                     }
 249                 }
 250             }
 251         }
 252         else // __a < 1
 253         {
 254             while (true)
 255             {
 256                 const result_type u  = uniformDist(g);
 257                 const result_type es = expDist(g);
 258
 259                 if (u <= 1.0 - alpha)
 260                 {
 261                     x = std::pow(u, 1.0 / alpha);
 262
 263                     if (x <= es)
 264                     {
 265                         break;
 266                     }
 267                 }
 268                 else
 269                 {
 270                     const result_type e = -std::log((1.0 - u) / alpha);
 271                     x                   = std::pow(1.0 - alpha + alpha * e, 1.0 / alpha);
 272
 273                     if (x <= e + es)
 274                     {
 275                         break;
 276                     }
 277                 }
 278             }
 279         }
 280         return x * param.beta();
 281     }
 282
 283     /*! \brief Return the first parameter of gamma distribution */
 284     result_type alpha() const { return param_.alpha(); }
 285
 286     /*! \brief Return the second parameter of gamma distribution */
 287     result_type beta() const { return param_.beta(); }
 288
 289     /*! \brief Return the full parameter class of gamma distribution */
 290     param_type param() const { return param_; }
 291
 292     /*! \brief Smallest value that can be returned from gamma distribution */
 293     result_type min() const { return 0; }
 294
 295     /*! \brief Largest value that can be returned from gamma distribution */
 296     result_type max() const { return std::numeric_limits<result_type>::infinity(); }
 297
 298     /*! \brief True if two gamma distributions will produce the same values.
 299      *
 300      * \param  x     Instance to compare with.
 301      */
 302     bool operator==(const GammaDistribution& x) const { return param_ == x.param_; }
 303
 304     /*! \brief True if two gamma distributions will produce different values.
 305      *
 306      * \param  x     Instance to compare with.
 307      */
 308     bool operator!=(const GammaDistribution& x) const { return !operator==(x); }
 309
 310 private:
 311     /*! \brief Internal value for parameters, can be overridden at generation time. */
 312     param_type param_;
 313
 314     GMX_DISALLOW_COPY_AND_ASSIGN(GammaDistribution);
 315 };
 316
 317 } // namespace gmx
 318
 319 #endif // GMX_RANDOM_GAMMADISTRIBUTION_H