5
\$\begingroup\$

I would like to communicate a piece of code which I hope will soon be broadly useful for everyone who programs in C++: A set of quasi-random number generators proposed for addition to boost.random. To keep the scope somewhat bounded, I will only submit the Sobol sequence for review (though there are two other sequences in the PR), and hopefully the advice will be generic enough to fix any other problems encountered.

First, usage:

#include <boost/random/sobol.hpp>
#include <boost/random/uniform_01.hpp>
#include <boost/random/variate_generator.hpp>

int main()
{
    static const std::size_t dimension = 4;

    // Create a generator
    typedef boost::variate_generator<boost::random::sobol64&, boost::uniform_01<double> > quasi_random_gen_t;

    // Initialize the engine to draw randomness out of thin air
    boost::random::sobol64 engine(dimension);

    // Glue the engine and the distribution together
    quasi_random_gen_t gen(engine, boost::uniform_01<double>());

    std::vector<double> sample(dimension);

    // At this point you can use std::generate, generate member f-n, etc.
    std::generate(sample.begin(), sample.end(), gen);
    engine.generate(sample.begin(), sample.end());
}

Next, base class:

/* boost random/detail/gray_coded_qrng_base.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_DETAIL_GRAY_CODED_QRNG_BASE_HPP
#define BOOST_RANDOM_DETAIL_GRAY_CODED_QRNG_BASE_HPP

#include <boost/random/detail/qrng_base.hpp>

#include <boost/multiprecision/integer.hpp> // lsb

//!\file
//!Describes the gray-coded quasi-random number generator base class template.

namespace boost {
namespace random {

namespace detail {

template<typename DerivedT, typename LatticeT, typename SizeT>
class gray_coded_qrng_base : public qrng_base<DerivedT, LatticeT, SizeT>
{
private:
  typedef gray_coded_qrng_base<DerivedT, LatticeT, SizeT> self_t;
  typedef qrng_base<DerivedT, LatticeT, SizeT> base_t;

  // The base needs to access modifying member f-ns, and we
  // don't want these functions to be available for the public use
  friend class qrng_base<DerivedT, LatticeT, SizeT>;

public:
  typedef typename base_t::size_type size_type;
  typedef typename LatticeT::value_type result_type;

  explicit gray_coded_qrng_base(std::size_t dimension)
    : base_t(dimension)
  {}

  // default copy c-tor is fine

  // default assignment operator is fine

  void seed()
  {
    base_t::set_zero();
    update_quasi(0);
  }

  void seed(size_type init)
  {
    this->curr_elem = 0;
    if (init != this->seq_count)
    {
      base_t::set_zero();

      this->seq_count = init++;
      init ^= (init >> 1);
      for (unsigned r = 0; init != 0; ++r, init >>= 1)
      {
        if (init & 1)
          update_quasi(r);
      }
    }
  }

private:
  void compute_seq(size_type cnt)
  {
    // Find the position of the least-significant zero in sequence count.
    // This is the bit that changes in the Gray-code representation as
    // the count is advanced.
    update_quasi(multiprecision::lsb(~cnt));
  }

  void update_quasi(unsigned r)
  {
    // Calculate the next state.
    for (std::size_t i = 0; i != this->dimension(); ++i)
      this->quasi_state[i] ^= this->lattice(r, i);
  }
};

}} // namespace detail::random

} // namespace boost

#endif // BOOST_RANDOM_DETAIL_GRAY_CODED_QRNG_BASE_HPP

Another base class:

/* boost random/detail/quasi_random_number_generator_base.hpp header file
 *
 * Copyright Justinas Vygintas Daugmaudis 2010-2017
 * 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_DETAIL_QRNG_BASE_HPP
#define BOOST_RANDOM_DETAIL_QRNG_BASE_HPP

#include <istream>
#include <ostream>

#include <stdexcept>
#include <vector>
#include <sstream>

#include <boost/random/detail/operators.hpp>

#include <boost/throw_exception.hpp>

//!\file
//!Describes the quasi-random number generator base class template.

namespace boost {
namespace random {

namespace detail {

template<typename DerivedT, typename LatticeT, typename SizeT>
class qrng_base
{
public:
  typedef SizeT size_type;
  typedef typename LatticeT::value_type result_type;

  explicit qrng_base(std::size_t dimension)
    // Guard against invalid dimensions before creating the lattice
    : lattice(prevent_zero_dimension(dimension))
    , quasi_state(dimension)
  {
    derived().seed();
  }

  // default copy c-tor is fine

  // default assignment operator is fine

  //!Returns: The dimension of of the quasi-random domain.
  //!
  //!Throws: nothing.
  std::size_t dimension() const { return quasi_state.size(); }

  //!Requirements: *this is mutable.
  //!
  //!Returns: Returns a successive element of an s-dimensional
  //!(s = X::dimension()) vector at each invocation. When all elements are
  //!exhausted, X::operator() begins anew with the starting element of a
  //!subsequent s-dimensional vector.
  //!
  //!Throws: overflow_error.
  result_type operator()()
  {
    return curr_elem != dimension() ? load_cached(): next_state();
  }

  //!Fills a range with quasi-random values.
  template<typename Iter> void generate(Iter first, Iter last)
  {
    for (; first != last; ++first)
      *first = this->operator()();
  }

  //!Requirements: *this is mutable.
  //!
  //!Effects: Advances *this state as if z consecutive
  //!X::operator() invocations were executed.
  //!
  //!Throws: overflow_error.
  void discard(size_type z)
  {
    const std::size_t dimension_value = dimension();

    std::size_t vec_n  = z / dimension_value;
    std::size_t elem_n = z - vec_n * dimension_value; // z % Dimension
    std::size_t vec_offset = vec_n + (curr_elem + elem_n) / dimension_value;
    // Discards vec_offset consecutive s-dimensional vectors
    discard_vector(vec_offset);
    // Sets up the proper position of the element-to-read
    curr_elem += (z - dimension_value * vec_offset);
  }

  //!Writes a @c DerivedT to a @c std::ostream.
  BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, DerivedT, s)
  {
    os << s.dimension() << " " << s.seq_count << " " << s.curr_elem;
    return os;
  }

  //!Reads a @c DerivedT from a @c std::istream.
  BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, DerivedT, s)
  {
    std::size_t dim, seed, z;
    if (is >> dim >> std::ws >> seed >> std::ws >> z) // initialize iff success!
    {
      if (s.dimension() != dim)
      {
        prevent_zero_dimension(dim);
        s.lattice.resize(dim);
        s.quasi_state.resize(dim);
      }
      // Fast-forward to the correct state
      s.seed(seed);
      s.discard(z);
    }
    return is;
  }

  //!Returns true if the two generators will produce identical sequences.
  BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(DerivedT, x, y)
  {
    const std::size_t dimension_value = x.dimension();

    // Note that two generators with different seq_counts and curr_elems can
    // produce the same sequence because the generator triple
    // (D, S, D) is equivalent to (D, S + 1, 0), where D is dimension, S -- seq_count,
    // and the last one is curr_elem.

    return (dimension_value == y.dimension()) &&
      (x.seq_count + (x.curr_elem / dimension_value) == y.seq_count + (y.curr_elem / dimension_value)) &&
      (x.curr_elem % dimension_value == y.curr_elem % dimension_value);
  }

  //!Returns true if the two generators will produce different sequences,
  BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(DerivedT)

protected:
  DerivedT& derived() throw()
  {
    return *static_cast<DerivedT * const>(this);
  }

  void set_zero()
  {
    curr_elem = 0;
    seq_count = 0;
    std::fill(quasi_state.begin(), quasi_state.end(), result_type /*zero*/());
  }

private:
  inline static std::size_t prevent_zero_dimension(std::size_t dimension)
  {
    if (dimension == 0)
      boost::throw_exception( std::invalid_argument("qrng_base: zero dimension") );
    return dimension;
  }

  // Load the result from the saved state.
  result_type load_cached()
  {
    return quasi_state[curr_elem++];
  }

  result_type next_state()
  {
    size_type new_seq = seq_count;
    if (++new_seq > seq_count)
    {
      derived().compute_seq(new_seq);
      seq_count = new_seq;
      curr_elem = 0;
      return load_cached();
    }
    boost::throw_exception( std::overflow_error("qrng_base: next_state") );
  }

  // Discards z consecutive s-dimensional vectors,
  // and preserves the position of the element-to-read
  void discard_vector(size_type z)
  {
    size_type inc_seq_count = seq_count + z;
    // Here we check that no overflow occurs before we
    // begin seeding the new value
    if (inc_seq_count > seq_count)
    {
      std::size_t tmp = curr_elem;

      derived().seed(inc_seq_count);

      curr_elem = tmp;
    }
    else if (inc_seq_count < seq_count) // Increment overflowed?
    {
      boost::throw_exception( std::overflow_error("discard_vector") );
    }
  }

protected:
  LatticeT lattice;
  std::size_t curr_elem;
  size_type seq_count;
  std::vector<result_type> quasi_state;
};

inline void dimension_assert(const char* generator, std::size_t dim, std::size_t maxdim)
{
  if (dim > maxdim)
  {
    std::ostringstream os;
    os << "The " << generator << " quasi-random number generator only supports up to "
      << maxdim << " dimensions.";
    boost::throw_exception( std::invalid_argument(os.str()) );
  }
}

}} // namespace detail::random

} // namespace boost

#endif // BOOST_RANDOM_DETAIL_QRNG_BASE_HPP

And the sobol.hpp file:

/* 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_data.hpp>
#include <boost/random/detail/gray_coded_qrng_base.hpp>

#include <limits>

#include <boost/cstdint.hpp>

#include <boost/multi_array.hpp>

//!\file
//!Describes the quasi-random number generator class template sobol_engine.
//!
//!\b Note: it is especially useful in conjunction with class template uniform_real.

namespace boost {
namespace random {

/** @cond */
namespace 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 IntType, typename SobolTables>
struct sobol_lattice
{
  typedef IntType value_type;

  BOOST_STATIC_CONSTANT(unsigned, bit_count = std::numeric_limits<IntType>::digits);

  // default copy c-tor is fine

  explicit sobol_lattice(std::size_t dimension)
  {
    resize(dimension);
  }

  void resize(std::size_t dimension)
  {
    detail::dimension_assert("Sobol", dimension, SobolTables::max_dimension);

    // Initialize the bit array
    bits.resize(boost::extents[bit_count][dimension]);

    // Initialize direction table in dimension 0
    for (unsigned k = 0; k != bit_count; ++k)
      bits[k][0] = static_cast<value_type>(1);

    // Initialize in remaining dimensions.
    for (std::size_t dim = 1; dim < dimension; ++dim)
    {
      const unsigned int 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::detail::find_msb(poly); // integer log2(poly)

      // set initial values of m from table
      for (unsigned k = 0; k != degree; ++k)
        bits[k][dim] = SobolTables::minit(k, dim-1);

      // Calculate remaining elements for this dimension,
      // as explained in Bratley+Fox, section 2.
      for (unsigned j = degree; j < bit_count; ++j)
      {
        unsigned int p_i = poly;
        bits[j][dim] = bits[j - degree][dim];
        for (unsigned k = 0; k != degree; ++k)
        {
          int rem = degree - k;
          bits[j][dim] ^= ((p_i & 1) * bits[j-rem][dim]) << rem;
          p_i >>= 1;
        }
      }
    }

    // Shift columns by appropriate power of 2.
    value_type p = static_cast<value_type>(1);
    for (int j = bit_count-1-1; j >= 0; --j, ++p)
      for (std::size_t dim = 0; dim != dimension; ++dim)
        bits[j][dim] <<= p;
  }

  value_type operator()(int i, int j) const
  {
    return bits[i][j];
  }

private:
  boost::multi_array<value_type, 2> bits;
};

} // namespace detail
/** @endcond */

//!class template sobol_engine implements a quasi-random number generator as described in
//! \blockquote
//![Bratley+Fox, TOMS 14, 88 (1988)]
//!and [Antonov+Saleev, USSR Comput. Maths. Math. Phys. 19, 252 (1980)]
//! \endblockquote
//!
//!In the following documentation @c X denotes the concrete class of the template
//!sobol_engine returning objects of type @c IntType, u and v are the values of @c X.
//!
//!Some member functions may throw exceptions of type @c std::overflow_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}\f$,
//!where `w`= std::numeric_limits<IntType>::digits.
template<typename IntType, typename SobolTables>
class sobol_engine : public detail::gray_coded_qrng_base<
                                sobol_engine<IntType, SobolTables>
                              , detail::sobol_lattice<IntType, SobolTables>
                              , IntType
                              >
{
  typedef sobol_engine<IntType, SobolTables> self_t;
  typedef detail::sobol_lattice<IntType, SobolTables> lattice_t;
  typedef detail::gray_coded_qrng_base<self_t, lattice_t, IntType> base_t;

public:
  typedef typename base_t::result_type result_type;

  /** @copydoc boost::random::niederreiter_base2_engine::min() */
  static BOOST_CONSTEXPR result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
  { return 0; }

  /** @copydoc boost::random::niederreiter_base2_engine::max() */
  static BOOST_CONSTEXPR result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
  { return (std::numeric_limits<result_type>::max)(); }

  //!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!==============================

  /** @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(IntType) */
  void seed(IntType 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(IntType) */
  void discard(IntType z)
  {
    base_t::discard(z);
  }
#endif // BOOST_RANDOM_DOXYGEN
};

/**
 * @attention This specialization of \sobol_engine supports up to 3667 dimensions.
 *
 * Data on the primitive binary polynomials \f$a\f$ and the corresponding starting values \f$m\f$,
 * for Sobol sequences in up to 21201 dimensions, 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
 *
 * For practical reasons the default table uses only the subset of binary polynomials \f$a\f$ that
 * satisfy \f$\forall a < 2^{16}\f$.
 *
 * However, it is possible to provide your own table should the default one be insufficient.
 */
typedef sobol_engine<boost::uint_least64_t,
  #ifdef BOOST_RANDOM_DOXYGEN
    sobol-tables
  #else
    detail::qrng_tables::sobol
  #endif
> sobol;

} // namespace random

} // namespace boost

#endif // BOOST_RANDOM_SOBOL_HPP
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.