Квази-Генераторы Случайных Чисел


Я хотел бы передать кусок кода, который я надеюсь скоро быть в целом полезным для всех, кто программирует на C++: набор квази-случайных чисел генераторы , предложенные для того , чтобы увеличить.случайные. Чтобы сохранить объем несколько ограничен, я только представить последовательность Соболь за комментарий (хотя есть два других последовательностей в ПР), и, надеюсь, советы будут достаточно общими, чтобы устранить любые другие проблемы.

Во-первых, использование:

#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());
}

Далее базовый класс:

/* 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

Другой базовый класс:

/* 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

И 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_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


267
5
задан 16 февраля 2018 в 05:02 Источник Поделиться
Комментарии