boost/random/sobol.hpp

0001 /* boost random/sobol.hpp header file
0002  *
0003  * Copyright Justinas Vygintas Daugmaudis 2010-2018
0004  * Distributed under the Boost Software License, Version 1.0. (See
0005  * accompanying file LICENSE_1_0.txt or copy at
0006  * http://www.boost.org/LICENSE_1_0.txt)
0007  */
0008
0009 #ifndef BOOST_RANDOM_SOBOL_HPP
0010 #define BOOST_RANDOM_SOBOL_HPP
0011
0012 #include <boost/random/detail/sobol_table.hpp>
0013 #include <boost/random/detail/gray_coded_qrng.hpp>
0014 #include <boost/assert.hpp>
0015
0016 namespace boost {
0017 namespace random {
0018
0019 /** @cond */
0020 namespace qrng_detail {
0021
0022 // sobol_lattice sets up the random-number generator to produce a Sobol
0023 // sequence of at most max dims-dimensional quasi-random vectors.
0024 // Adapted from ACM TOMS algorithm 659, see
0025
0026 // http://doi.acm.org/10.1145/42288.214372
0027
0028 template<typename UIntType, unsigned w, typename SobolTables>
0029 struct sobol_lattice
0030 {
0031   typedef UIntType value_type;
0032
0033   BOOST_STATIC_ASSERT(w > 0u);
0034   BOOST_STATIC_CONSTANT(unsigned, bit_count = w);
0035
0036 private:
0037   typedef std::vector<value_type> container_type;
0038
0039 public:
0040   explicit sobol_lattice(std::size_t dimension)
0041   {
0042     resize(dimension);
0043   }
0044
0045   // default copy c-tor is fine
0046
0047   void resize(std::size_t dimension)
0048   {
0049     dimension_assert("Sobol", dimension, SobolTables::max_dimension);
0050
0051     // Initialize the bit array
0052     container_type cj(bit_count * dimension);
0053
0054     // Initialize direction table in dimension 0
0055     for (unsigned k = 0; k != bit_count; ++k)
0056       cj[dimension*k] = static_cast<value_type>(1);
0057
0058     // Initialize in remaining dimensions.
0059     for (std::size_t dim = 1; dim < dimension; ++dim)
0060     {
0061       const typename SobolTables::value_type poly = SobolTables::polynomial(dim-1);
0062       if (poly > (std::numeric_limits<value_type>::max)()) {
0063         boost::throw_exception( std::range_error("sobol: polynomial value outside the given value type range") );
0064       }
0065       const unsigned degree = qrng_detail::msb(poly); // integer log2(poly)
0066
0067       // set initial values of m from table
0068       for (unsigned k = 0; k != degree; ++k)
0069         cj[dimension*k + dim] = SobolTables::minit(dim-1, k);
0070
0071       // Calculate remaining elements for this dimension,
0072       // as explained in Bratley+Fox, section 2.
0073       for (unsigned j = degree; j < bit_count; ++j)
0074       {
0075         typename SobolTables::value_type p_i = poly;
0076         const std::size_t bit_offset = dimension*j + dim;
0077
0078         cj[bit_offset] = cj[dimension*(j-degree) + dim];
0079         for (unsigned k = 0; k != degree; ++k, p_i >>= 1)
0080         {
0081           int rem = degree - k;
0082           cj[bit_offset] ^= ((p_i & 1) * cj[dimension*(j-rem) + dim]) << rem;
0083         }
0084       }
0085     }
0086
0087     // Shift columns by appropriate power of 2.
0088     unsigned p = 1u;
0089     for (int j = bit_count-1-1; j >= 0; --j, ++p)
0090     {
0091       const std::size_t bit_offset = dimension * j;
0092       for (std::size_t dim = 0; dim != dimension; ++dim)
0093         cj[bit_offset + dim] <<= p;
0094     }
0095
0096     bits.swap(cj);
0097   }
0098
0099   typename container_type::const_iterator iter_at(std::size_t n) const
0100   {
0101     BOOST_ASSERT(!(n > bits.size()));
0102     return bits.begin() + n;
0103   }
0104
0105 private:
0106   container_type bits;
0107 };
0108
0109 } // namespace qrng_detail
0110
0111 typedef detail::qrng_tables::sobol default_sobol_table;
0112
0113 /** @endcond */
0114
0115 //!Instantiations of class template sobol_engine model a \quasi_random_number_generator.
0116 //!The sobol_engine uses the algorithm described in
0117 //! \blockquote
0118 //![Bratley+Fox, TOMS 14, 88 (1988)]
0119 //!and [Antonov+Saleev, USSR Comput. Maths. Math. Phys. 19, 252 (1980)]
0120 //! \endblockquote
0121 //!
0122 //!\attention sobol_engine skips trivial zeroes at the start of the sequence. For example, the beginning
0123 //!of the 2-dimensional Sobol sequence in @c uniform_01 distribution will look like this:
0124 //!\code{.cpp}
0125 //!0.5, 0.5,
0126 //!0.75, 0.25,
0127 //!0.25, 0.75,
0128 //!0.375, 0.375,
0129 //!0.875, 0.875,
0130 //!...
0131 //!\endcode
0132 //!
0133 //!In the following documentation @c X denotes the concrete class of the template
0134 //!sobol_engine returning objects of type @c UIntType, u and v are the values of @c X.
0135 //!
0136 //!Some member functions may throw exceptions of type @c std::range_error. This
0137 //!happens when the quasi-random domain is exhausted and the generator cannot produce
0138 //!any more values. The length of the low discrepancy sequence is given by \f$L=Dimension \times (2^{w} - 1)\f$.
0139 template<typename UIntType, unsigned w, typename SobolTables = default_sobol_table>
0140 class sobol_engine
0141   : public qrng_detail::gray_coded_qrng<
0142       qrng_detail::sobol_lattice<UIntType, w, SobolTables>
0143     >
0144 {
0145   typedef qrng_detail::sobol_lattice<UIntType, w, SobolTables> lattice_t;
0146   typedef qrng_detail::gray_coded_qrng<lattice_t> base_t;
0147
0148 public:
0149   //!Effects: Constructs the default `s`-dimensional Sobol quasi-random number generator.
0150   //!
0151   //!Throws: bad_alloc, invalid_argument, range_error.
0152   explicit sobol_engine(std::size_t s)
0153     : base_t(s)
0154   {}
0155
0156   // default copy c-tor is fine
0157
0158 #ifdef BOOST_RANDOM_DOXYGEN
0159   //=========================Doxygen needs this!==============================
0160   typedef UIntType result_type;
0161
0162   /** @copydoc boost::random::niederreiter_base2_engine::min() */
0163   static BOOST_CONSTEXPR result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
0164   { return (base_t::min)(); }
0165
0166   /** @copydoc boost::random::niederreiter_base2_engine::max() */
0167   static BOOST_CONSTEXPR result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
0168   { return (base_t::max)(); }
0169
0170   /** @copydoc boost::random::niederreiter_base2_engine::dimension() */
0171   std::size_t dimension() const { return base_t::dimension(); }
0172
0173   /** @copydoc boost::random::niederreiter_base2_engine::seed() */
0174   void seed()
0175   {
0176     base_t::seed();
0177   }
0178
0179   /** @copydoc boost::random::niederreiter_base2_engine::seed(UIntType) */
0180   void seed(UIntType init)
0181   {
0182     base_t::seed(init);
0183   }
0184
0185   /** @copydoc boost::random::niederreiter_base2_engine::operator()() */
0186   result_type operator()()
0187   {
0188     return base_t::operator()();
0189   }
0190
0191   /** @copydoc boost::random::niederreiter_base2_engine::discard(boost::uintmax_t) */
0192   void discard(boost::uintmax_t z)
0193   {
0194     base_t::discard(z);
0195   }
0196
0197   /** Returns true if the two generators will produce identical sequences of outputs. */
0198   BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(sobol_engine, x, y)
0199   { return static_cast<const base_t&>(x) == y; }
0200
0201   /** Returns true if the two generators will produce different sequences of outputs. */
0202   BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(sobol_engine)
0203
0204   /** Writes the textual representation of the generator to a @c std::ostream. */
0205   BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, sobol_engine, s)
0206   { return os << static_cast<const base_t&>(s); }
0207
0208   /** Reads the textual representation of the generator from a @c std::istream. */
0209   BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, sobol_engine, s)
0210   { return is >> static_cast<base_t&>(s); }
0211
0212 #endif // BOOST_RANDOM_DOXYGEN
0213 };
0214
0215 /**
0216  * @attention This specialization of \sobol_engine supports up to 3667 dimensions.
0217  *
0218  * Data on the primitive binary polynomials `a` and the corresponding starting values `m`
0219  * for Sobol sequences in up to 21201 dimensions was taken from
0220  *
0221  *  @blockquote
0222  *  S. Joe and F. Y. Kuo, Constructing Sobol sequences with better two-dimensional projections,
0223  *  SIAM J. Sci. Comput. 30, 2635-2654 (2008).
0224  *  @endblockquote
0225  *
0226  * 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
0227  *
0228  * For practical reasons the default table uses only the subset of binary polynomials `a` < 2<sup>16</sup>.
0229  *
0230  * However, it is possible to provide your own table to \sobol_engine should the default one be insufficient.
0231  */
0232 typedef sobol_engine<boost::uint_least64_t, 64u, default_sobol_table> sobol;
0233
0234 } // namespace random
0235
0236 } // namespace boost
0237
0238 #endif // BOOST_RANDOM_SOBOL_HPP