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),
 145             beta_(beta)
 146         {
 147         }
 148
 149         /*! \brief Return first parameter */
 150         result_type alpha() const { return alpha_; }
 151         /*! \brief Return second parameter */
 152         result_type beta() const { return beta_; }
 153
 154         /*! \brief True if two parameter sets will return the same distribution.
 155          *
 156          * \param x  Instance to compare with.
 157          */
 158         bool operator==(const param_type& x) const
 159         {
 160             return alpha_ == x.alpha_ && beta_ == x.beta_;
 161         }
 162
 163         /*! \brief True if two parameter sets will return different distributions
 164          *
 165          * \param x  Instance to compare with.
 166          */
 167         bool operator!=(const param_type& x) const { return !operator==(x); }
 168     };
 169
 170     /*! \brief Construct new distribution with given floating-point parameters.
 171      *
 172      * \param alpha  First parameter of gamma distribution
 173      * \param beta   Second parameter of gamma distribution
 174      */
 175     explicit GammaDistribution(result_type alpha = 1.0, result_type beta = 1.0) :
 176         param_(param_type(alpha, beta))
 177     {
 178     }
 179
 180     /*! \brief Construct new distribution from parameter class
 181      *
 182      * \param param  Parameter class as defined inside gmx::GammaDistribution.
 183      */
 184     explicit GammaDistribution(const param_type& param) : param_(param) {}
 185
 186     /*! \brief Flush all internal saved values  */
 187     void reset() {}
 188
 189     /*! \brief Return values from gamma distribution with internal parameters
 190      *
 191      *  \tparam Rng   Random engine class
 192      *
 193      *  \param  g     Random engine
 194      */
 195     template<class Rng>
 196     result_type operator()(Rng& g)
 197     {
 198         return (*this)(g, param_);
 199     }
 200
 201     /*! \brief Return value from gamma distribution with given parameters
 202      *
 203      *  \tparam Rng   Random engine class
 204      *
 205      *  \param  g     Random engine
 206      *  \param  param Parameters to use
 207      */
 208     template<class Rng>
 209     result_type operator()(Rng& g, const param_type& param)
 210     {
 211         result_type                          alpha = param.alpha();
 212         UniformRealDistribution<result_type> uniformDist(0, 1);
 213         ExponentialDistribution<result_type> expDist;
 214
 215         result_type x;
 216
 217         if (alpha == 1.0)
 218         {
 219             x = expDist(g);
 220         }
 221         else if (alpha > 1.0)
 222         {
 223             const result_type b = alpha - 1.0;
 224             const result_type c = 3.0 * alpha - result_type(0.75);
 225
 226             while (true)
 227             {
 228                 const result_type u = uniformDist(g);
 229                 const result_type v = uniformDist(g);
 230                 const result_type w = u * (1 - u);
 231
 232                 if (w != 0)
 233                 {
 234                     const result_type y = std::sqrt(c / w) * (u - result_type(0.5));
 235                     x                   = b + y;
 236
 237                     if (x >= 0)
 238                     {
 239                         const result_type z = 64 * w * w * w * v * v;
 240
 241                         if (z <= 1.0 - 2.0 * y * y / x)
 242                         {
 243                             break;
 244                         }
 245                         if (std::log(z) <= 2.0 * (b * std::log(x / b) - y))
 246                         {
 247                             break;
 248                         }
 249                     }
 250                 }
 251             }
 252         }
 253         else // __a < 1
 254         {
 255             while (true)
 256             {
 257                 const result_type u  = uniformDist(g);
 258                 const result_type es = expDist(g);
 259
 260                 if (u <= 1.0 - alpha)
 261                 {
 262                     x = std::pow(u, 1.0 / alpha);
 263
 264                     if (x <= es)
 265                     {
 266                         break;
 267                     }
 268                 }
 269                 else
 270                 {
 271                     const result_type e = -std::log((1.0 - u) / alpha);
 272                     x                   = std::pow(1.0 - alpha + alpha * e, 1.0 / alpha);
 273
 274                     if (x <= e + es)
 275                     {
 276                         break;
 277                     }
 278                 }
 279             }
 280         }
 281         return x * param.beta();
 282     }
 283
 284     /*! \brief Return the first parameter of gamma distribution */
 285     result_type alpha() const { return param_.alpha(); }
 286
 287     /*! \brief Return the second parameter of gamma distribution */
 288     result_type beta() const { return param_.beta(); }
 289
 290     /*! \brief Return the full parameter class of gamma distribution */
 291     param_type param() const { return param_; }
 292
 293     /*! \brief Smallest value that can be returned from gamma distribution */
 294     result_type min() const { return 0; }
 295
 296     /*! \brief Largest value that can be returned from gamma distribution */
 297     result_type max() const { return std::numeric_limits<result_type>::infinity(); }
 298
 299     /*! \brief True if two gamma distributions will produce the same values.
 300      *
 301      * \param  x     Instance to compare with.
 302      */
 303     bool operator==(const GammaDistribution& x) const { return param_ == x.param_; }
 304
 305     /*! \brief True if two gamma distributions will produce different values.
 306      *
 307      * \param  x     Instance to compare with.
 308      */
 309     bool operator!=(const GammaDistribution& x) const { return !operator==(x); }
 310
 311 private:
 312     /*! \brief Internal value for parameters, can be overridden at generation time. */
 313     param_type param_;
 314
 315     GMX_DISALLOW_COPY_AND_ASSIGN(GammaDistribution);
 316 };
 317
 318 } // namespace gmx
 319
 320 #endif // GMX_RANDOM_GAMMADISTRIBUTION_H