boost/random/sobol.hpp
/* boost random/sobol.hpp header file * * Copyright Justinas Vygintas Daugmaudis 2010-2018 * Distributed under the Boost Software License, Version 1.0. (See * accompanying file LICENSE_1_0.txt or copy at * http://www.boost.org/LICENSE_1_0.txt) */ #ifndef BOOST_RANDOM_SOBOL_HPP #define BOOST_RANDOM_SOBOL_HPP #include <boost/random/detail/sobol_table.hpp> #include <boost/random/detail/gray_coded_qrng.hpp> namespace boost { namespace random { /** @cond */ namespace qrng_detail { // sobol_lattice sets up the random-number generator to produce a Sobol // sequence of at most max dims-dimensional quasi-random vectors. // Adapted from ACM TOMS algorithm 659, see // http://doi.acm.org/10.1145/42288.214372 template<typename UIntType, unsigned w, typename SobolTables> struct sobol_lattice { typedef UIntType value_type; BOOST_STATIC_ASSERT(w > 0u); BOOST_STATIC_CONSTANT(unsigned, bit_count = w); private: typedef std::vector<value_type> container_type; public: explicit sobol_lattice(std::size_t dimension) { resize(dimension); } // default copy c-tor is fine void resize(std::size_t dimension) { dimension_assert("Sobol", dimension, SobolTables::max_dimension); // Initialize the bit array container_type cj(bit_count * dimension); // Initialize direction table in dimension 0 for (unsigned k = 0; k != bit_count; ++k) cj[dimension*k] = static_cast<value_type>(1); // Initialize in remaining dimensions. for (std::size_t dim = 1; dim < dimension; ++dim) { const typename SobolTables::value_type poly = SobolTables::polynomial(dim-1); if (poly > std::numeric_limits<value_type>::max()) { boost::throw_exception( std::range_error("sobol: polynomial value outside the given value type range") ); } const unsigned degree = multiprecision::msb(poly); // integer log2(poly) // set initial values of m from table for (unsigned k = 0; k != degree; ++k) cj[dimension*k + dim] = SobolTables::minit(dim-1, k); // Calculate remaining elements for this dimension, // as explained in Bratley+Fox, section 2. for (unsigned j = degree; j < bit_count; ++j) { typename SobolTables::value_type p_i = poly; const std::size_t bit_offset = dimension*j + dim; cj[bit_offset] = cj[dimension*(j-degree) + dim]; for (unsigned k = 0; k != degree; ++k, p_i >>= 1) { int rem = degree - k; cj[bit_offset] ^= ((p_i & 1) * cj[dimension*(j-rem) + dim]) << rem; } } } // Shift columns by appropriate power of 2. unsigned p = 1u; for (int j = bit_count-1-1; j >= 0; --j, ++p) { const std::size_t bit_offset = dimension * j; for (std::size_t dim = 0; dim != dimension; ++dim) cj[bit_offset + dim] <<= p; } bits.swap(cj); } typename container_type::const_iterator iter_at(std::size_t n) const { BOOST_ASSERT(!(n > bits.size())); return bits.begin() + n; } private: container_type bits; }; } // namespace qrng_detail typedef detail::qrng_tables::sobol default_sobol_table; /** @endcond */ //!Instantiations of class template sobol_engine model a \quasi_random_number_generator. //!The sobol_engine uses the algorithm described in //! \blockquote //![Bratley+Fox, TOMS 14, 88 (1988)] //!and [Antonov+Saleev, USSR Comput. Maths. Math. Phys. 19, 252 (1980)] //! \endblockquote //! //!\attention sobol_engine skips trivial zeroes at the start of the sequence. For example, the beginning //!of the 2-dimensional Sobol sequence in @c uniform_01 distribution will look like this: //!\code{.cpp} //!0.5, 0.5, //!0.75, 0.25, //!0.25, 0.75, //!0.375, 0.375, //!0.875, 0.875, //!... //!\endcode //! //!In the following documentation @c X denotes the concrete class of the template //!sobol_engine returning objects of type @c UIntType, u and v are the values of @c X. //! //!Some member functions may throw exceptions of type @c std::range_error. This //!happens when the quasi-random domain is exhausted and the generator cannot produce //!any more values. The length of the low discrepancy sequence is given by \f$L=Dimension \times (2^{w} - 1)\f$. template<typename UIntType, unsigned w, typename SobolTables = default_sobol_table> class sobol_engine : public qrng_detail::gray_coded_qrng< qrng_detail::sobol_lattice<UIntType, w, SobolTables> > { typedef qrng_detail::sobol_lattice<UIntType, w, SobolTables> lattice_t; typedef qrng_detail::gray_coded_qrng<lattice_t> base_t; public: //!Effects: Constructs the default `s`-dimensional Sobol quasi-random number generator. //! //!Throws: bad_alloc, invalid_argument, range_error. explicit sobol_engine(std::size_t s) : base_t(s) {} // default copy c-tor is fine #ifdef BOOST_RANDOM_DOXYGEN //=========================Doxygen needs this!============================== typedef UIntType result_type; /** @copydoc boost::random::niederreiter_base2_engine::min() */ static BOOST_CONSTEXPR result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () { return (base_t::min)(); } /** @copydoc boost::random::niederreiter_base2_engine::max() */ static BOOST_CONSTEXPR result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () { return (base_t::max)(); } /** @copydoc boost::random::niederreiter_base2_engine::dimension() */ std::size_t dimension() const { return base_t::dimension(); } /** @copydoc boost::random::niederreiter_base2_engine::seed() */ void seed() { base_t::seed(); } /** @copydoc boost::random::niederreiter_base2_engine::seed(UIntType) */ void seed(UIntType init) { base_t::seed(init); } /** @copydoc boost::random::niederreiter_base2_engine::operator()() */ result_type operator()() { return base_t::operator()(); } /** @copydoc boost::random::niederreiter_base2_engine::discard(boost::uintmax_t) */ void discard(boost::uintmax_t z) { base_t::discard(z); } /** Returns true if the two generators will produce identical sequences of outputs. */ BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(sobol_engine, x, y) { return static_cast<const base_t&>(x) == y; } /** Returns true if the two generators will produce different sequences of outputs. */ BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(sobol_engine) /** Writes the textual representation of the generator to a @c std::ostream. */ BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, sobol_engine, s) { return os << static_cast<const base_t&>(s); } /** Reads the textual representation of the generator from a @c std::istream. */ BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, sobol_engine, s) { return is >> static_cast<base_t&>(s); } #endif // BOOST_RANDOM_DOXYGEN }; /** * @attention This specialization of \sobol_engine supports up to 3667 dimensions. * * Data on the primitive binary polynomials `a` and the corresponding starting values `m` * for Sobol sequences in up to 21201 dimensions was taken from * * @blockquote * S. Joe and F. Y. Kuo, Constructing Sobol sequences with better two-dimensional projections, * SIAM J. Sci. Comput. 30, 2635-2654 (2008). * @endblockquote * * 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 * * For practical reasons the default table uses only the subset of binary polynomials `a` < 2<sup>16</sup>. * * However, it is possible to provide your own table to \sobol_engine should the default one be insufficient. */ typedef sobol_engine<boost::uint_least64_t, 64u, default_sobol_table> sobol; } // namespace random } // namespace boost #endif // BOOST_RANDOM_SOBOL_HPP