src/gromacs/random/tabulatednormaldistribution.h

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