src/gromacs/random/normaldistribution.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 normal distribution
  38  *
  39  * Portable version of the normal distribution that generates the same sequence
  40  * on all platforms. Since stdlibc++ and libc++ provide different sequences
  41  * we prefer this one so unit tests produce the same values on all platforms.
  42  *
  43  * \author Erik Lindahl <erik.lindahl@gmail.com>
  44  * \inpublicapi
  45  * \ingroup module_random
  46  */
  47
  48 #ifndef GMX_RANDOM_NORMALDISTRIBUTION_H
  49 #define GMX_RANDOM_NORMALDISTRIBUTION_H
  50
  51 #include <cmath>
  52
  53 #include <limits>
  54 #include <memory>
  55
  56 #include "gromacs/random/uniformrealdistribution.h"
  57 #include "gromacs/utility/classhelpers.h"
  58
  59 /*
  60  * The portable version of the normal distribution (to make sure we get the same
  61  * values on all platforms) has been modified from the LLVM libcxx headers,
  62  * distributed under the MIT license:
  63  *
  64  * Copyright (c) The LLVM compiler infrastructure
  65  *
  66  * Permission is hereby granted, free of charge, to any person obtaining a copy
  67  * of this software and associated documentation files (the "Software"), to deal
  68  * in the Software without restriction, including without limitation the rights
  69  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  70  * copies of the Software, and to permit persons to whom the Software is
  71  * furnished to do so, subject to the following conditions:
  72  *
  73  * The above copyright notice and this permission notice shall be included in
  74  * all copies or substantial portions of the Software.
  75  *
  76  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  77  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  78  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  79  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  80  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  81  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  82  * THE SOFTWARE.
  83  */
  84
  85 namespace gmx
  86 {
  87
  88 /*! \brief Normal distribution
  89  *
  90  *  The C++ standard library does provide a normal distribution, but even
  91  *  though they all sample from the normal distribution different standard
  92  *  library implementations appear to return different sequences of numbers
  93  *  for the same random number generator. To make it easier to use GROMACS
  94  *  unit tests that depend on random numbers we have our own implementation.
  95  *
  96  *  Be warned that the normal distribution draws values from the random engine
  97  *  in a loop, so you want to make sure you use a random stream with a
  98  *  very large margin to make sure you do not run out of random numbers
  99  *  in an unlucky case (which will lead to an exception with the GROMACS
 100  *  default random engine).
 101  *
 102  * \tparam RealType Floating-point type, real by default in GROMACS.
 103  */
 104 template<class RealType = real>
 105 class NormalDistribution
 106 {
 107 public:
 108     /*! \brief Type of values returned */
 109     typedef RealType result_type;
 110
 111     /*! \brief Normal distribution parameters */
 112     class param_type
 113     {
 114         /*! \brief Mean of normal distribution */
 115         result_type mean_;
 116         /*! \brief Standard deviation of distribution */
 117         result_type stddev_;
 118
 119     public:
 120         /*! \brief Reference back to the distribution class */
 121         typedef NormalDistribution distribution_type;
 122
 123         /*! \brief Construct parameter block
 124          *
 125          * \param mean     Mean of normal distribution
 126          * \param stddev   Standard deviation of normal distribution
 127          */
 128         explicit param_type(result_type mean = 0.0, result_type stddev = 1.0) :
 129             mean_(mean),
 130             stddev_(stddev)
 131         {
 132         }
 133
 134         /*! \brief Return first parameter */
 135         result_type mean() const { return mean_; }
 136         /*! \brief Return second parameter */
 137         result_type stddev() const { return stddev_; }
 138
 139         /*! \brief True if two parameter sets will return the same normal distribution.
 140          *
 141          * \param x  Instance to compare with.
 142          */
 143         bool operator==(const param_type& x) const
 144         {
 145             return mean_ == x.mean_ && stddev_ == x.stddev_;
 146         }
 147
 148         /*! \brief True if two parameter sets will return different normal distributions
 149          *
 150          * \param x  Instance to compare with.
 151          */
 152         bool operator!=(const param_type& x) const { return !operator==(x); }
 153     };
 154
 155 public:
 156     /*! \brief Construct new distribution with given floating-point parameters.
 157      *
 158      * \param mean     Mean of normal distribution
 159      * \param stddev   Standard deviation of normal distribution
 160      */
 161     explicit NormalDistribution(result_type mean = 0.0, result_type stddev = 1.0) :
 162         param_(param_type(mean, stddev)),
 163         hot_(false),
 164         saved_(0)
 165     {
 166     }
 167
 168     /*! \brief Construct new distribution from parameter class
 169      *
 170      * \param param  Parameter class as defined inside gmx::NormalDistribution.
 171      */
 172     explicit NormalDistribution(const param_type& param) : param_(param), hot_(false), saved_(0) {}
 173
 174     /*! \brief Flush all internal saved values  */
 175     void reset() { hot_ = false; }
 176
 177     /*! \brief Return values from normal distribution with internal parameters
 178      *
 179      *  \tparam Rng   Random engine class
 180      *
 181      *  \param  g     Random engine
 182      */
 183     template<class Rng>
 184     result_type operator()(Rng& g)
 185     {
 186         return (*this)(g, param_);
 187     }
 188
 189     /*! \brief Return value from normal distribution with given parameters
 190      *
 191      *  \tparam Rng   Random engine class
 192      *
 193      *  \param  g     Random engine
 194      *  \param  param Parameters to use
 195      */
 196     template<class Rng>
 197     result_type operator()(Rng& g, const param_type& param)
 198     {
 199         result_type result;
 200
 201         if (hot_)
 202         {
 203             hot_   = false;
 204             result = saved_;
 205         }
 206         else
 207         {
 208             UniformRealDistribution<result_type> uniformDist(-1.0, 1.0);
 209             result_type                          u;
 210             result_type                          v;
 211             result_type                          s;
 212
 213             do
 214             {
 215                 u = uniformDist(g);
 216                 v = uniformDist(g);
 217                 s = u * u + v * v;
 218             } while (s > 1.0 || s == 0.0);
 219
 220             s      = std::sqrt(-2.0 * std::log(s) / s);
 221             saved_ = v * s;
 222             hot_   = true;
 223             result = u * s;
 224         }
 225         return result * param.stddev() + param.mean();
 226     }
 227
 228     /*! \brief Return the mean of the normal distribution */
 229     result_type mean() const { return param_.mean(); }
 230
 231     /*! \brief Return the standard deviation of the normal distribution */
 232     result_type stddev() const { return param_.stddev(); }
 233
 234     /*! \brief Return the full parameter class of the normal distribution */
 235     param_type param() const { return param_; }
 236
 237     /*! \brief Smallest value that can be returned from normal distribution */
 238     result_type min() const { return -std::numeric_limits<result_type>::infinity(); }
 239
 240     /*! \brief Largest value that can be returned from normal distribution */
 241     result_type max() const { return std::numeric_limits<result_type>::infinity(); }
 242
 243     /*! \brief True if two normal distributions will produce the same values.
 244      *
 245      * \param  x     Instance to compare with.
 246      */
 247     bool operator==(const NormalDistribution& x) const
 248     {
 249         /* Equal if: Params are identical, and saved-state is identical,
 250          * and if we have something saved, it must be identical.
 251          */
 252         return param_ == x.param_ && hot_ == x.hot_ && (!hot_ || saved_ == x.saved_);
 253     }
 254
 255     /*! \brief True if two normal distributions will produce different values.
 256      *
 257      * \param  x     Instance to compare with.
 258      */
 259     bool operator!=(const NormalDistribution& x) const { return !operator==(x); }
 260
 261 private:
 262     /*! \brief Internal value for parameters, can be overridden at generation time. */
 263     param_type param_;
 264     /*! \brief True if there is a saved result to return */
 265     bool hot_;
 266     /*! \brief The saved result to return - only valid if hot_ is true */
 267     result_type saved_;
 268
 269     GMX_DISALLOW_COPY_AND_ASSIGN(NormalDistribution);
 270 };
 271
 272
 273 } // namespace gmx
 274
 275 #endif // GMX_RANDOM_NORMALDISTRIBUTION_H