src/gromacs/random/tabulatednormaldistribution.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2015,2016, 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 Tabulated normal distribution
  38  *
  39  * A very fast normal distribution, but with limited resolution.
  40  *
  41  * \author Erik Lindahl <erik.lindahl@gmail.com>
  42  * \inpublicapi
  43  * \ingroup module_random
  44  */
  45
  46 #ifndef GMX_RANDOM_TABULATEDNORMALDISTRIBUTION_H
  47 #define GMX_RANDOM_TABULATEDNORMALDISTRIBUTION_H
  48
  49 #include <cmath>
  50
  51 #include <limits>
  52 #include <vector>
  53
  54 #include "gromacs/math/utilities.h"
  55 #include "gromacs/utility/basedefinitions.h"
  56 #include "gromacs/utility/classhelpers.h"
  57 #include "gromacs/utility/real.h"
  58
  59 namespace gmx
  60 {
  61
  62 /*! \brief Tabulated normal random distribution
  63  *
  64  *  Random distribution compatible with C++11 distributions - it can be
  65  *  used with any C++11 random engine.
  66  *
  67  *  \tparam RealType  Type of the return value. Float or double. Note that
  68  *                    GROMACS uses "real" type by default in contrast to the C++11
  69  *                    standard library, to avoid double/float conversions.
  70  *  \tparam tableBits Size of the table, specified in bits. The storage
  71  *                    space required is sizeof(RealType)*2^tableBits. To
  72  *                    keep things sane this is limited to 24 bits.
  73  *
  74  *  Some stochastic integrators depend on drawing a lot of normal
  75  *  distribution random numbers quickly, but in many cases the only
  76  *  important property is the distribution - given the noise in forces
  77  *  we do not need very high resolution.
  78  *  This distribution uses an internal table to return samples from a
  79  *  normal distribution with limited resolution. By default the table uses
  80  *  14 bits, but this is specified with a template parameter.
  81  *
  82  *  Since this distribution only uses tableBits bits per value generated,
  83  *  the values draw from the random engine are used for several results.
  84  *  To make sure you get a reproducible result when using counter-based
  85  *  random engines (such as ThreeFry2x64), remember to call the reset()
  86  *  method to cancel the internal memory of the distribution.
  87  *
  88  *  \note For modern NUMA systems, you likely want to use separate
  89  *        distributions for each thread, and make sure they are initialized
  90  *        on the CPU where they will run, so the table is placed in that
  91  *        NUMA memory pool.
  92  *  \note The finite table resolution means this distribution will NOT
  93  *        return arbitrarily small/large values, but with e.g. 14 bits
  94  *        the results are limited to roughly +/- 4 standard deviations.
  95  */
  96 template<class RealType = real, unsigned int tableBits = 14>
  97 class TabulatedNormalDistribution
  98 {
  99     static_assert(tableBits <= 24, "Normal distribution table is limited to 24bits (64MB in single precision)");
 100
 101     public:
 102         /*! \brief  Type of normal distribution results */
 103         typedef RealType result_type;
 104
 105         /*! \brief  Normal distribution parameter class (mean and stddev) */
 106         class param_type
 107         {
 108             public:
 109                 /*! \brief The type of distribution the parameters describe */
 110                 typedef TabulatedNormalDistribution distribution_type;
 111
 112                 /*! \brief Constructor. Default is classical distr. with mean 0, stddev 1.
 113                  *
 114                  * \param mean     Expectation value.
 115                  * \param stddev   Standard deviation.
 116                  *
 117                  */
 118                 explicit param_type(result_type mean = 0.0, result_type stddev = 1.0)
 119                     : mean_(mean), stddev_(stddev) {}
 120
 121                 /*! \brief Return mean parameter of normal distribution */
 122                 result_type
 123                 mean() const { return mean_; }
 124
 125                 /*! \brief Return standard deviation parameter of normal distribution */
 126                 result_type
 127                 stddev() const { return stddev_; }
 128
 129                 /*! \brief True if two sets of normal distributions parameters are identical
 130                  *
 131                  * \param x Instance to compare with.
 132                  */
 133                 bool
 134                 operator==(const param_type &x) const
 135                 {
 136                     return (mean_ == x.mean_ && stddev_ == x.stddev_);
 137                 }
 138
 139                 /*! \brief True if two sets of normal distributions parameters are different.
 140                  *
 141                  * \param x Instance to compare with.
 142                  */
 143                 bool
 144                 operator!=(const param_type &x) const { return !operator==(x); }
 145
 146             private:
 147                 /*! \brief Internal storage for mean of normal distribution */
 148                 result_type mean_;
 149                 /*! \brief Internal storage for standard deviation of normal distribution */
 150                 result_type stddev_;
 151         };
 152
 153     private:
 154
 155         /*! \brief Fill the table with values for the normal distribution
 156          *
 157          *  This routine returns a new reference to a std::vector allocated on
 158          *  the heap. It will only be called to generate the static
 159          *  constant table data at initialization time. While it is technically
 160          *  returning memory on the heap, this does not use/leak more memory
 161          *  than what we would use by having the same table as a static member
 162          *  of the class.
 163          */
 164         static
 165         std::vector<RealType> &
 166         // cppcheck-suppress unusedPrivateFunction
 167         makeTable()
 168         {
 169             /* Fill the table with the integral of a gaussian distribution:
 170              */
 171             std::size_t            tableSize        = 1 << tableBits;
 172             std::size_t            halfSize         = (1 << tableBits)/2;
 173             double                 invSize          = 1.0/tableSize;
 174             double                 factor           = std::sqrt(2.0*M_PI);
 175             double                 x                = 0.5*factor*invSize;
 176
 177             std::vector<RealType>* table  = new std::vector<RealType>(1ULL << tableBits);
 178
 179             for (std::size_t i = 0; i < halfSize; i++)
 180             {
 181                 if (i > 0)
 182                 {
 183                     double dx;
 184
 185                     if (i < halfSize-1)
 186                     {
 187                         double invNormal = factor*std::exp(0.5*x*x);
 188                         /* det is larger than 0 for all x, except the last */
 189                         double det = 1.0 - 2.0*invSize*x*invNormal;
 190                         dx = (1.0 - std::sqrt(det))/x;
 191                     }
 192                     else
 193                     {
 194                         dx = 1.0/x;
 195                     }
 196                     x = x + dx;
 197                 }
 198                 table->at(halfSize-1-i) = -x;
 199                 table->at(halfSize+i)   =  x;
 200             }
 201             return *table;
 202         }
 203
 204     public:
 205
 206         /*! \brief Construct new normal distribution with specified mean & stdddev.
 207          *
 208          *  \param mean    Mean value of tabulated normal distribution
 209          *  \param stddev  Standard deviation of tabulated normal distribution
 210          */
 211         explicit TabulatedNormalDistribution(result_type mean = 0.0, result_type stddev = 1.0 )
 212             : param_(param_type(mean, stddev)), savedRandomBits_(0), savedRandomBitsLeft_(0)
 213         {
 214         }
 215
 216         /*! \brief Construct new normal distribution from parameter type.
 217          *
 218          *  \param param Parameter class containing mean and standard deviation.
 219          */
 220         explicit TabulatedNormalDistribution(  const param_type &param )
 221             : param_(param), savedRandomBits_(0), savedRandomBitsLeft_(0)
 222         {
 223         }
 224
 225         /*! \brief Smallest value that can be generated in normal distrubiton.
 226          *
 227          * \note The smallest value is not -infinity with a table, but it
 228          *       depends on the table resolution. With 14 bits, this is roughly
 229          *       four standard deviations below the mean.
 230          */
 231         result_type
 232         min() const
 233         {
 234             return c_table_[0];
 235         }
 236
 237         /*! \brief Largest value that can be generated in normal distribution.
 238          *
 239          * \note The largest value is not infinity with a table, but it
 240          *       depends on the table resolution. With 14 bits, this is roughly
 241          *       four standard deviations above the mean.
 242          */
 243         result_type
 244         max() const
 245         {
 246             return c_table_[c_table_.size()-1];
 247         }
 248
 249         /*! \brief Mean of the present normal distribution */
 250         result_type
 251         mean() const
 252         {
 253             return param_.mean();
 254         }
 255
 256         /*! \brief Standard deviation of the present normal distribution */
 257
 258         result_type
 259         stddev() const
 260         {
 261             return param_.stddev();
 262         }
 263
 264         /*! \brief The parameter class (mean & stddev) of the normal distribution */
 265         param_type
 266         param() const
 267         {
 268             return param_;
 269         }
 270
 271         /*! \brief Clear all internal saved random bits from the random engine */
 272         void
 273         reset()
 274         {
 275             savedRandomBitsLeft_ = 0;
 276         }
 277
 278         /*! \brief Return normal distribution value specified by internal parameters.
 279          *
 280          * \tparam Rng   Random engine type used to provide uniform random bits.
 281          * \param  g     Random engine of class Rng. For normal GROMACS usage
 282          *               you likely want to use ThreeFry2x64.
 283          */
 284         template<class Rng>
 285         result_type
 286         operator()(Rng &g)
 287         {
 288             return (*this)(g, param_);
 289         }
 290
 291         /*! \brief Return normal distribution value specified by given parameters
 292          *
 293          * \tparam Rng   Random engine type used to provide uniform random bits.
 294          * \param  g     Random engine of class Rng. For normal GROMACS usage
 295          *               you likely want to use ThreeFry2x64.
 296          * \param  param Parameters used to specify normal distribution.
 297          */
 298         template<class Rng>
 299         result_type
 300         operator()(Rng &g, const param_type &param)
 301         {
 302             if (savedRandomBitsLeft_ < tableBits)
 303             {
 304                 // We do not know whether the generator g returns 64 or 32 bits,
 305                 // since g is not known when we construct this class.
 306                 // To keep things simple, we always draw one random number,
 307                 // store it in our 64-bit value, and set the number of active bits.
 308                 // For tableBits up to 16 this will be as efficient both with 32
 309                 // and 64 bit random engines when drawing multiple numbers
 310                 // (our default value is 14), It also avoids drawing multiple
 311                 // 32-bit random numbres even if we just call this routine for a
 312                 // single result.
 313                 savedRandomBits_     = static_cast<gmx_uint64_t>(g());
 314                 savedRandomBitsLeft_ = std::numeric_limits<typename Rng::result_type>::digits;
 315             }
 316             result_type value        = c_table_[savedRandomBits_ & ( (1ULL << tableBits) - 1 ) ];
 317             savedRandomBits_       >>= tableBits;
 318             savedRandomBitsLeft_    -= tableBits;
 319             return param.mean() + value * param.stddev();
 320         }
 321
 322         /*!\brief Check if two tabulated normal distributions have identical states.
 323          *
 324          * \param  x     Instance to compare with.
 325          */
 326         bool
 327         operator==(const TabulatedNormalDistribution<RealType, tableBits> &x) const
 328         {
 329             return (param_ == x.param_ &&
 330                     savedRandomBits_ == x.savedRandomBits_ &&
 331                     savedRandomBitsLeft_ == x.savedRandomBitsLeft_);
 332         }
 333
 334         /*!\brief Check if two tabulated normal distributions have different states.
 335          *
 336          * \param  x     Instance to compare with.
 337          */
 338         bool
 339         operator!=(const TabulatedNormalDistribution<RealType, tableBits> &x) const
 340         {
 341             return !operator==(x);
 342         }
 343
 344     private:
 345         /*! \brief Parameters of normal distribution (mean and stddev) */
 346         param_type                                                   param_;
 347         /*! \brief Array with tabluated values of normal distribution */
 348         static const std::vector<RealType>                           c_table_;
 349         /*! \brief Saved output from random engine, shifted tableBits right each time */
 350         gmx_uint64_t                                                 savedRandomBits_;
 351         /*! \brief Number of valid bits remaining i savedRandomBits_ */
 352         unsigned int                                                 savedRandomBitsLeft_;
 353
 354         GMX_DISALLOW_COPY_AND_ASSIGN(TabulatedNormalDistribution);
 355 };
 356
 357 // MSVC does not handle extern template class members correctly even in MSVC 2015,
 358 // so in that case we have to instantiate in every object using it. In addition,
 359 // doxygen is convinced this defines a function (which leads to crashes in our python
 360 // scripts), so to avoid confusion we hide it from doxygen too.
 361 #if (!defined(_MSC_VER) || defined(__INTEL_COMPILER)) && !defined(DOXYGEN)
 362 // Declaration of template specialization
 363 template<>
 364 const std::vector<real> TabulatedNormalDistribution<real, 14>::c_table_;
 365
 366 extern template
 367 const std::vector<real> TabulatedNormalDistribution<real, 14>::c_table_;
 368 #endif
 369
 370 // Instantiation for all tables without specialization
 371 template<class RealType, unsigned int tableBits>
 372 const std::vector<RealType> TabulatedNormalDistribution<RealType, tableBits>::c_table_ = TabulatedNormalDistribution<RealType, tableBits>::makeTable();
 373
 374 }      // namespace gmx
 375
 376 #endif // GMX_RANDOM_TABULATEDNORMALDISTRIBUTION_H