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,
111 "Normal distribution table is limited to 24bits (64MB in single precision)");
114 /*! \brief Type of normal distribution results */
115 typedef RealType result_type;
117 /*! \brief Normal distribution parameter class (mean and stddev) */
121 /*! \brief The type of distribution the parameters describe */
122 typedef TabulatedNormalDistribution distribution_type;
124 /*! \brief Constructor. Default is classical distr. with mean 0, stddev 1.
126 * \param mean Expectation value.
127 * \param stddev Standard deviation.
130 explicit param_type(result_type mean = 0.0, result_type stddev = 1.0) :
136 /*! \brief Return mean parameter of normal distribution */
137 result_type mean() const { return mean_; }
139 /*! \brief Return standard deviation parameter of normal distribution */
140 result_type stddev() const { return stddev_; }
142 /*! \brief True if two sets of normal distributions parameters are identical
144 * \param x Instance to compare with.
146 bool operator==(const param_type& x) const
148 return (mean_ == x.mean_ && stddev_ == x.stddev_);
151 /*! \brief True if two sets of normal distributions parameters are different.
153 * \param x Instance to compare with.
155 bool 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.
172 static std::array<RealType, 1 << tableBits> makeTable()
174 /* Fill the table with the integral of a gaussian distribution, which
175 * corresponds to the inverse error function.
176 * We avoid integrating a gaussian numerically, since that leads to
177 * some loss-of-precision which also accumulates so it is worse for
178 * larger indices in the table. */
179 constexpr std::size_t tableSize = 1 << tableBits;
180 constexpr std::size_t halfSize = tableSize / 2;
181 constexpr double invHalfSize = 1.0 / halfSize;
183 std::array<RealType, tableSize> table;
185 // Fill in all but the extremal entries of the table
186 for (std::size_t i = 0; i < halfSize - 1; i++)
188 double r = (i + 0.5) * invHalfSize;
189 double x = std::sqrt(2.0) * erfinv(r);
191 table.at(halfSize - 1 - i) = -x;
192 table.at(halfSize + i) = x;
194 // We want to fill in the extremal table entries with
195 // values that make the total variance equal to 1, so
196 // measure the variance by summing the squares of the
197 // other values of the distribution, starting from the
199 double sumOfSquares = 0;
200 for (std::size_t i = 1; i < halfSize; i++)
202 double value = table.at(i);
203 sumOfSquares += value * value;
205 double missingVariance = 1.0 - 2.0 * sumOfSquares / tableSize;
206 GMX_RELEASE_ASSERT(missingVariance > 0,
207 "Incorrect computation of tabulated normal distribution");
208 double extremalValue = std::sqrt(0.5 * missingVariance * tableSize);
209 table.at(0) = -extremalValue;
210 table.back() = extremalValue;
216 /*! \brief Construct new normal distribution with specified mean & stdddev.
218 * \param mean Mean value of tabulated normal distribution
219 * \param stddev Standard deviation of tabulated normal distribution
221 explicit TabulatedNormalDistribution(result_type mean = 0.0, result_type stddev = 1.0) :
222 param_(param_type(mean, stddev)),
224 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& param) :
235 savedRandomBitsLeft_(0)
239 /*! \brief Smallest value that can be generated in normal distrubiton.
241 * \note The smallest value is not -infinity with a table, but it
242 * depends on the table resolution. With 14 bits, this is roughly
243 * four standard deviations below the mean.
245 result_type min() const { return c_table_[0]; }
247 /*! \brief Largest value that can be generated in normal distribution.
249 * \note The largest value is not infinity with a table, but it
250 * depends on the table resolution. With 14 bits, this is roughly
251 * four standard deviations above the mean.
253 result_type max() const { return c_table_[c_table_.size() - 1]; }
255 /*! \brief Mean of the present normal distribution */
256 result_type mean() const { return param_.mean(); }
258 /*! \brief Standard deviation of the present normal distribution */
260 result_type stddev() const { return param_.stddev(); }
262 /*! \brief The parameter class (mean & stddev) of the normal distribution */
263 param_type param() const { return param_; }
265 /*! \brief Clear all internal saved random bits from the random engine */
266 void reset() { savedRandomBitsLeft_ = 0; }
268 /*! \brief Return normal distribution value specified by internal parameters.
270 * \tparam Rng Random engine type used to provide uniform random bits.
271 * \param g Random engine of class Rng. For normal GROMACS usage
272 * you likely want to use ThreeFry2x64.
275 result_type operator()(Rng& g)
277 return (*this)(g, param_);
280 /*! \brief Return normal distribution value specified by given parameters
282 * \tparam Rng Random engine type used to provide uniform random bits.
283 * \param g Random engine of class Rng. For normal GROMACS usage
284 * you likely want to use ThreeFry2x64.
285 * \param param Parameters used to specify normal distribution.
288 result_type operator()(Rng& g, const param_type& param)
290 if (savedRandomBitsLeft_ < tableBits)
292 // We do not know whether the generator g returns 64 or 32 bits,
293 // since g is not known when we construct this class.
294 // To keep things simple, we always draw one random number,
295 // store it in our 64-bit value, and set the number of active bits.
296 // For tableBits up to 16 this will be as efficient both with 32
297 // and 64 bit random engines when drawing multiple numbers
298 // (our default value is
299 // c_TabulatedNormalDistributionDefaultBits == 14). It
300 // also avoids drawing multiple 32-bit random numbers
301 // even if we just call this routine for a single
303 savedRandomBits_ = static_cast<uint64_t>(g());
304 savedRandomBitsLeft_ = std::numeric_limits<typename Rng::result_type>::digits;
306 result_type value = c_table_[savedRandomBits_ & ((1ULL << tableBits) - 1)];
307 savedRandomBits_ >>= tableBits;
308 savedRandomBitsLeft_ -= tableBits;
309 return param.mean() + value * param.stddev();
312 /*!\brief Check if two tabulated normal distributions have identical states.
314 * \param x Instance to compare with.
316 bool operator==(const TabulatedNormalDistribution<RealType, tableBits>& x) const
318 return (param_ == x.param_ && savedRandomBits_ == x.savedRandomBits_
319 && savedRandomBitsLeft_ == x.savedRandomBitsLeft_);
322 /*!\brief Check if two tabulated normal distributions have different states.
324 * \param x Instance to compare with.
326 bool operator!=(const TabulatedNormalDistribution<RealType, tableBits>& x) const
328 return !operator==(x);
332 /*! \brief Parameters of normal distribution (mean and stddev) */
334 /*! \brief Array with tabluated values of normal distribution */
335 static const std::array<RealType, 1 << tableBits> c_table_;
336 /*! \brief Saved output from random engine, shifted tableBits right each time */
337 uint64_t savedRandomBits_;
338 /*! \brief Number of valid bits remaining i savedRandomBits_ */
339 unsigned int savedRandomBitsLeft_;
341 GMX_DISALLOW_COPY_AND_ASSIGN(TabulatedNormalDistribution);
344 // MSVC does not handle extern template class members correctly even in MSVC 2015,
345 // so in that case we have to instantiate in every object using it. In addition,
346 // doxygen is convinced this defines a function (which leads to crashes in our python
347 // scripts), so to avoid confusion we hide it from doxygen too.
348 #if !defined(_MSC_VER) && !defined(DOXYGEN)
349 // Declaration of template specialization
351 const std::array<real, 1 << detail::c_TabulatedNormalDistributionDefaultBits> TabulatedNormalDistribution<>::c_table_;
353 extern template const std::array<real, 1 << detail::c_TabulatedNormalDistributionDefaultBits>
354 TabulatedNormalDistribution<>::c_table_;
357 // Instantiation for all tables without specialization
358 template<class RealType, unsigned int tableBits>
359 const std::array<RealType, 1 << tableBits> TabulatedNormalDistribution<RealType, tableBits>::c_table_ =
360 TabulatedNormalDistribution<RealType, tableBits>::makeTable();
364 #endif // GMX_RANDOM_TABULATEDNORMALDISTRIBUTION_H