9
3

sobol.hpp 7.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237
  1. /* boost random/sobol.hpp header file
  2. *
  3. * Copyright Justinas Vygintas Daugmaudis 2010-2018
  4. * Distributed under the Boost Software License, Version 1.0. (See
  5. * accompanying file LICENSE_1_0.txt or copy at
  6. * http://www.boost.org/LICENSE_1_0.txt)
  7. */
  8. #ifndef BOOST_RANDOM_SOBOL_HPP
  9. #define BOOST_RANDOM_SOBOL_HPP
  10. #include <boost/random/detail/sobol_table.hpp>
  11. #include <boost/random/detail/gray_coded_qrng.hpp>
  12. namespace boost {
  13. namespace random {
  14. /** @cond */
  15. namespace qrng_detail {
  16. // sobol_lattice sets up the random-number generator to produce a Sobol
  17. // sequence of at most max dims-dimensional quasi-random vectors.
  18. // Adapted from ACM TOMS algorithm 659, see
  19. // http://doi.acm.org/10.1145/42288.214372
  20. template<typename UIntType, unsigned w, typename SobolTables>
  21. struct sobol_lattice
  22. {
  23. typedef UIntType value_type;
  24. BOOST_STATIC_ASSERT(w > 0u);
  25. BOOST_STATIC_CONSTANT(unsigned, bit_count = w);
  26. private:
  27. typedef std::vector<value_type> container_type;
  28. public:
  29. explicit sobol_lattice(std::size_t dimension)
  30. {
  31. resize(dimension);
  32. }
  33. // default copy c-tor is fine
  34. void resize(std::size_t dimension)
  35. {
  36. dimension_assert("Sobol", dimension, SobolTables::max_dimension);
  37. // Initialize the bit array
  38. container_type cj(bit_count * dimension);
  39. // Initialize direction table in dimension 0
  40. for (unsigned k = 0; k != bit_count; ++k)
  41. cj[dimension*k] = static_cast<value_type>(1);
  42. // Initialize in remaining dimensions.
  43. for (std::size_t dim = 1; dim < dimension; ++dim)
  44. {
  45. const typename SobolTables::value_type poly = SobolTables::polynomial(dim-1);
  46. if (poly > std::numeric_limits<value_type>::max()) {
  47. boost::throw_exception( std::range_error("sobol: polynomial value outside the given value type range") );
  48. }
  49. const unsigned degree = multiprecision::msb(poly); // integer log2(poly)
  50. // set initial values of m from table
  51. for (unsigned k = 0; k != degree; ++k)
  52. cj[dimension*k + dim] = SobolTables::minit(dim-1, k);
  53. // Calculate remaining elements for this dimension,
  54. // as explained in Bratley+Fox, section 2.
  55. for (unsigned j = degree; j < bit_count; ++j)
  56. {
  57. typename SobolTables::value_type p_i = poly;
  58. const std::size_t bit_offset = dimension*j + dim;
  59. cj[bit_offset] = cj[dimension*(j-degree) + dim];
  60. for (unsigned k = 0; k != degree; ++k, p_i >>= 1)
  61. {
  62. int rem = degree - k;
  63. cj[bit_offset] ^= ((p_i & 1) * cj[dimension*(j-rem) + dim]) << rem;
  64. }
  65. }
  66. }
  67. // Shift columns by appropriate power of 2.
  68. unsigned p = 1u;
  69. for (int j = bit_count-1-1; j >= 0; --j, ++p)
  70. {
  71. const std::size_t bit_offset = dimension * j;
  72. for (std::size_t dim = 0; dim != dimension; ++dim)
  73. cj[bit_offset + dim] <<= p;
  74. }
  75. bits.swap(cj);
  76. }
  77. typename container_type::const_iterator iter_at(std::size_t n) const
  78. {
  79. BOOST_ASSERT(!(n > bits.size()));
  80. return bits.begin() + n;
  81. }
  82. private:
  83. container_type bits;
  84. };
  85. } // namespace qrng_detail
  86. typedef detail::qrng_tables::sobol default_sobol_table;
  87. /** @endcond */
  88. //!Instantiations of class template sobol_engine model a \quasi_random_number_generator.
  89. //!The sobol_engine uses the algorithm described in
  90. //! \blockquote
  91. //![Bratley+Fox, TOMS 14, 88 (1988)]
  92. //!and [Antonov+Saleev, USSR Comput. Maths. Math. Phys. 19, 252 (1980)]
  93. //! \endblockquote
  94. //!
  95. //!\attention sobol_engine skips trivial zeroes at the start of the sequence. For example, the beginning
  96. //!of the 2-dimensional Sobol sequence in @c uniform_01 distribution will look like this:
  97. //!\code{.cpp}
  98. //!0.5, 0.5,
  99. //!0.75, 0.25,
  100. //!0.25, 0.75,
  101. //!0.375, 0.375,
  102. //!0.875, 0.875,
  103. //!...
  104. //!\endcode
  105. //!
  106. //!In the following documentation @c X denotes the concrete class of the template
  107. //!sobol_engine returning objects of type @c UIntType, u and v are the values of @c X.
  108. //!
  109. //!Some member functions may throw exceptions of type @c std::range_error. This
  110. //!happens when the quasi-random domain is exhausted and the generator cannot produce
  111. //!any more values. The length of the low discrepancy sequence is given by \f$L=Dimension \times (2^{w} - 1)\f$.
  112. template<typename UIntType, unsigned w, typename SobolTables = default_sobol_table>
  113. class sobol_engine
  114. : public qrng_detail::gray_coded_qrng<
  115. qrng_detail::sobol_lattice<UIntType, w, SobolTables>
  116. >
  117. {
  118. typedef qrng_detail::sobol_lattice<UIntType, w, SobolTables> lattice_t;
  119. typedef qrng_detail::gray_coded_qrng<lattice_t> base_t;
  120. public:
  121. //!Effects: Constructs the default `s`-dimensional Sobol quasi-random number generator.
  122. //!
  123. //!Throws: bad_alloc, invalid_argument, range_error.
  124. explicit sobol_engine(std::size_t s)
  125. : base_t(s)
  126. {}
  127. // default copy c-tor is fine
  128. #ifdef BOOST_RANDOM_DOXYGEN
  129. //=========================Doxygen needs this!==============================
  130. typedef UIntType result_type;
  131. /** @copydoc boost::random::niederreiter_base2_engine::min() */
  132. static BOOST_CONSTEXPR result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
  133. { return (base_t::min)(); }
  134. /** @copydoc boost::random::niederreiter_base2_engine::max() */
  135. static BOOST_CONSTEXPR result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
  136. { return (base_t::max)(); }
  137. /** @copydoc boost::random::niederreiter_base2_engine::dimension() */
  138. std::size_t dimension() const { return base_t::dimension(); }
  139. /** @copydoc boost::random::niederreiter_base2_engine::seed() */
  140. void seed()
  141. {
  142. base_t::seed();
  143. }
  144. /** @copydoc boost::random::niederreiter_base2_engine::seed(UIntType) */
  145. void seed(UIntType init)
  146. {
  147. base_t::seed(init);
  148. }
  149. /** @copydoc boost::random::niederreiter_base2_engine::operator()() */
  150. result_type operator()()
  151. {
  152. return base_t::operator()();
  153. }
  154. /** @copydoc boost::random::niederreiter_base2_engine::discard(boost::uintmax_t) */
  155. void discard(boost::uintmax_t z)
  156. {
  157. base_t::discard(z);
  158. }
  159. /** Returns true if the two generators will produce identical sequences of outputs. */
  160. BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(sobol_engine, x, y)
  161. { return static_cast<const base_t&>(x) == y; }
  162. /** Returns true if the two generators will produce different sequences of outputs. */
  163. BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(sobol_engine)
  164. /** Writes the textual representation of the generator to a @c std::ostream. */
  165. BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, sobol_engine, s)
  166. { return os << static_cast<const base_t&>(s); }
  167. /** Reads the textual representation of the generator from a @c std::istream. */
  168. BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, sobol_engine, s)
  169. { return is >> static_cast<base_t&>(s); }
  170. #endif // BOOST_RANDOM_DOXYGEN
  171. };
  172. /**
  173. * @attention This specialization of \sobol_engine supports up to 3667 dimensions.
  174. *
  175. * Data on the primitive binary polynomials `a` and the corresponding starting values `m`
  176. * for Sobol sequences in up to 21201 dimensions was taken from
  177. *
  178. * @blockquote
  179. * S. Joe and F. Y. Kuo, Constructing Sobol sequences with better two-dimensional projections,
  180. * SIAM J. Sci. Comput. 30, 2635-2654 (2008).
  181. * @endblockquote
  182. *
  183. * See the original tables up to dimension 21201: https://web.archive.org/web/20170802022909/http://web.maths.unsw.edu.au/~fkuo/sobol/new-joe-kuo-6.21201
  184. *
  185. * For practical reasons the default table uses only the subset of binary polynomials `a` < 2<sup>16</sup>.
  186. *
  187. * However, it is possible to provide your own table to \sobol_engine should the default one be insufficient.
  188. */
  189. typedef sobol_engine<boost::uint_least64_t, 64u, default_sobol_table> sobol;
  190. } // namespace random
  191. } // namespace boost
  192. #endif // BOOST_RANDOM_SOBOL_HPP