2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2015,2016,2018,2019, 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.
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.
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.
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.
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.
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.
37 * \brief Tabulated normal distribution
39 * A very fast normal distribution, but with limited resolution.
41 * \author Erik Lindahl <erik.lindahl@gmail.com>
43 * \ingroup module_random
46 #ifndef GMX_RANDOM_TABULATEDNORMALDISTRIBUTION_H
47 #define GMX_RANDOM_TABULATEDNORMALDISTRIBUTION_H
54 #include "gromacs/math/functions.h"
55 #include "gromacs/math/utilities.h"
56 #include "gromacs/utility/basedefinitions.h"
57 #include "gromacs/utility/classhelpers.h"
58 #include "gromacs/utility/gmxassert.h"
59 #include "gromacs/utility/real.h"
67 //! Number of bits that determines the resolution of the lookup table for the normal distribution.
68 constexpr int c_TabulatedNormalDistributionDefaultBits = 14;
72 /*! \brief Tabulated normal random distribution
74 * Random distribution compatible with C++11 distributions - it can be
75 * used with any C++11 random engine.
77 * \tparam RealType Type of the return value. Float or double. Note that
78 * GROMACS uses "real" type by default in contrast to the C++11
79 * standard library, to avoid double/float conversions.
80 * \tparam tableBits Size of the table, specified in bits. The storage
81 * space required is sizeof(RealType)*2^tableBits. To
82 * keep things sane this is limited to 24 bits.
84 * Some stochastic integrators depend on drawing a lot of normal
85 * distribution random numbers quickly, but in many cases the only
86 * important property is the distribution - given the noise in forces
87 * we do not need very high resolution.
88 * This distribution uses an internal table to return samples from a
89 * normal distribution with limited resolution. By default the table
90 * uses c_TabulatedNormalDistributionDefaultBits bits, but this is
91 * specified with a template parameter.
93 * Since this distribution only uses tableBits bits per value generated,
94 * the values draw from the random engine are used for several results.
95 * To make sure you get a reproducible result when using counter-based
96 * random engines (such as ThreeFry2x64), remember to call the reset()
97 * method to cancel the internal memory of the distribution.
99 * \note For modern NUMA systems, you likely want to use separate
100 * distributions for each thread, and make sure they are initialized
101 * on the CPU where they will run, so the table is placed in that
103 * \note The finite table resolution means this distribution will NOT
104 * return arbitrarily small/large values, but with e.g. 14 bits
105 * the results are limited to roughly +/- 4 standard deviations.
107 template<class RealType = real, unsigned int tableBits = detail::c_TabulatedNormalDistributionDefaultBits>
108 class TabulatedNormalDistribution
110 static_assert(tableBits <= 24, "Normal distribution table is limited to 24bits (64MB in single precision)");
113 /*! \brief Type of normal distribution results */
114 typedef RealType result_type;
116 /*! \brief Normal distribution parameter class (mean and stddev) */
120 /*! \brief The type of distribution the parameters describe */
121 typedef TabulatedNormalDistribution distribution_type;
123 /*! \brief Constructor. Default is classical distr. with mean 0, stddev 1.
125 * \param mean Expectation value.
126 * \param stddev Standard deviation.
129 explicit param_type(result_type mean = 0.0, result_type stddev = 1.0)
130 : mean_(mean), stddev_(stddev) {}
132 /*! \brief Return mean parameter of normal distribution */
134 mean() const { return mean_; }
136 /*! \brief Return standard deviation parameter of normal distribution */
138 stddev() const { return stddev_; }
140 /*! \brief True if two sets of normal distributions parameters are identical
142 * \param x Instance to compare with.
145 operator==(const param_type &x) const
147 return (mean_ == x.mean_ && stddev_ == x.stddev_);
150 /*! \brief True if two sets of normal distributions parameters are different.
152 * \param x Instance to compare with.
155 operator!=(const param_type &x) const { return !operator==(x); }
158 /*! \brief Internal storage for mean of normal distribution */
160 /*! \brief Internal storage for standard deviation of normal distribution */
164 /*! \brief Fill the table with values for the normal distribution
166 * This routine returns a new a std::array with the table data.
168 * This routine is used to help construct objects of this class,
169 * and is exposed only to permit testing. Normal code should not
170 * need to call this function.
173 std::array<RealType, 1<<tableBits>
176 /* Fill the table with the integral of a gaussian distribution, which
177 * corresponds to the inverse error function.
178 * We avoid integrating a gaussian numerically, since that leads to
179 * some loss-of-precision which also accumulates so it is worse for
180 * larger indices in the table. */
181 constexpr std::size_t tableSize = 1 << tableBits;
182 constexpr std::size_t halfSize = tableSize/2;
183 constexpr double invHalfSize = 1.0/halfSize;
185 std::array<RealType, tableSize> table;
187 // Fill in all but the extremal entries of the table
188 for (std::size_t i = 0; i < halfSize-1; i++)
190 double r = (i + 0.5) * invHalfSize;
191 double x = std::sqrt(2.0) * erfinv(r);
193 table.at(halfSize-1-i) = -x;
194 table.at(halfSize+i) = x;
196 // We want to fill in the extremal table entries with
197 // values that make the total variance equal to 1, so
198 // measure the variance by summing the squares of the
199 // other values of the distribution, starting from the
201 double sumOfSquares = 0;
202 for (std::size_t i = 1; i < halfSize; i++)
204 double value = table.at(i);
205 sumOfSquares += value * value;
207 double missingVariance = 1.0 - 2.0*sumOfSquares/tableSize;
208 GMX_RELEASE_ASSERT(missingVariance > 0, "Incorrect computation of tabulated normal distribution");
209 double extremalValue = std::sqrt(0.5*missingVariance*tableSize);
210 table.at(0) = -extremalValue;
211 table.back() = extremalValue;
218 /*! \brief Construct new normal distribution with specified mean & stdddev.
220 * \param mean Mean value of tabulated normal distribution
221 * \param stddev Standard deviation of tabulated normal distribution
223 explicit TabulatedNormalDistribution(result_type mean = 0.0, result_type stddev = 1.0 )
224 : param_(param_type(mean, stddev)), savedRandomBits_(0), savedRandomBitsLeft_(0)
228 /*! \brief Construct new normal distribution from parameter type.
230 * \param param Parameter class containing mean and standard deviation.
232 explicit TabulatedNormalDistribution( const param_type ¶m )
233 : param_(param), savedRandomBits_(0), savedRandomBitsLeft_(0)
237 /*! \brief Smallest value that can be generated in normal distrubiton.
239 * \note The smallest 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 below the mean.
249 /*! \brief Largest value that can be generated in normal distribution.
251 * \note The largest value is not infinity with a table, but it
252 * depends on the table resolution. With 14 bits, this is roughly
253 * four standard deviations above the mean.
258 return c_table_[c_table_.size()-1];
261 /*! \brief Mean of the present normal distribution */
265 return param_.mean();
268 /*! \brief Standard deviation of the present normal distribution */
273 return param_.stddev();
276 /*! \brief The parameter class (mean & stddev) of the normal distribution */
283 /*! \brief Clear all internal saved random bits from the random engine */
287 savedRandomBitsLeft_ = 0;
290 /*! \brief Return normal distribution value specified by internal parameters.
292 * \tparam Rng Random engine type used to provide uniform random bits.
293 * \param g Random engine of class Rng. For normal GROMACS usage
294 * you likely want to use ThreeFry2x64.
300 return (*this)(g, param_);
303 /*! \brief Return normal distribution value specified by given parameters
305 * \tparam Rng Random engine type used to provide uniform random bits.
306 * \param g Random engine of class Rng. For normal GROMACS usage
307 * you likely want to use ThreeFry2x64.
308 * \param param Parameters used to specify normal distribution.
312 operator()(Rng &g, const param_type ¶m)
314 if (savedRandomBitsLeft_ < tableBits)
316 // We do not know whether the generator g returns 64 or 32 bits,
317 // since g is not known when we construct this class.
318 // To keep things simple, we always draw one random number,
319 // store it in our 64-bit value, and set the number of active bits.
320 // For tableBits up to 16 this will be as efficient both with 32
321 // and 64 bit random engines when drawing multiple numbers
322 // (our default value is
323 // c_TabulatedNormalDistributionDefaultBits == 14). It
324 // also avoids drawing multiple 32-bit random numbers
325 // even if we just call this routine for a single
327 savedRandomBits_ = static_cast<uint64_t>(g());
328 savedRandomBitsLeft_ = std::numeric_limits<typename Rng::result_type>::digits;
330 result_type value = c_table_[savedRandomBits_ & ( (1ULL << tableBits) - 1 ) ];
331 savedRandomBits_ >>= tableBits;
332 savedRandomBitsLeft_ -= tableBits;
333 return param.mean() + value * param.stddev();
336 /*!\brief Check if two tabulated normal distributions have identical states.
338 * \param x Instance to compare with.
341 operator==(const TabulatedNormalDistribution<RealType, tableBits> &x) const
343 return (param_ == x.param_ &&
344 savedRandomBits_ == x.savedRandomBits_ &&
345 savedRandomBitsLeft_ == x.savedRandomBitsLeft_);
348 /*!\brief Check if two tabulated normal distributions have different states.
350 * \param x Instance to compare with.
353 operator!=(const TabulatedNormalDistribution<RealType, tableBits> &x) const
355 return !operator==(x);
359 /*! \brief Parameters of normal distribution (mean and stddev) */
361 /*! \brief Array with tabluated values of normal distribution */
362 static const std::array<RealType, 1 << tableBits> c_table_;
363 /*! \brief Saved output from random engine, shifted tableBits right each time */
364 uint64_t savedRandomBits_;
365 /*! \brief Number of valid bits remaining i savedRandomBits_ */
366 unsigned int savedRandomBitsLeft_;
368 GMX_DISALLOW_COPY_AND_ASSIGN(TabulatedNormalDistribution);
371 // MSVC does not handle extern template class members correctly even in MSVC 2015,
372 // so in that case we have to instantiate in every object using it. In addition,
373 // doxygen is convinced this defines a function (which leads to crashes in our python
374 // scripts), so to avoid confusion we hide it from doxygen too.
375 #if !defined(_MSC_VER) && !defined(DOXYGEN)
376 // Declaration of template specialization
378 const std::array<real, 1<<detail::c_TabulatedNormalDistributionDefaultBits> TabulatedNormalDistribution<>::c_table_;
381 const std::array<real, 1<<detail::c_TabulatedNormalDistributionDefaultBits> TabulatedNormalDistribution<>::c_table_;
384 // Instantiation for all tables without specialization
385 template<class RealType, unsigned int tableBits>
386 const std::array<RealType, 1<<tableBits> TabulatedNormalDistribution<RealType, tableBits>::c_table_ = TabulatedNormalDistribution<RealType, tableBits>::makeTable();
390 #endif // GMX_RANDOM_TABULATEDNORMALDISTRIBUTION_H