discrete_distribution.hpp 5.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160
  1. //---------------------------------------------------------------------------//
  2. // Copyright (c) 2014 Roshan <thisisroshansmail@gmail.com>
  3. //
  4. // Distributed under the Boost Software License, Version 1.0
  5. // See accompanying file LICENSE_1_0.txt or copy at
  6. // http://www.boost.org/LICENSE_1_0.txt
  7. //
  8. // See http://boostorg.github.com/compute for more information.
  9. //---------------------------------------------------------------------------//
  10. #ifndef BOOST_COMPUTE_RANDOM_DISCRETE_DISTRIBUTION_HPP
  11. #define BOOST_COMPUTE_RANDOM_DISCRETE_DISTRIBUTION_HPP
  12. #include <numeric>
  13. #include <boost/config.hpp>
  14. #include <boost/type_traits.hpp>
  15. #include <boost/static_assert.hpp>
  16. #include <boost/compute/command_queue.hpp>
  17. #include <boost/compute/function.hpp>
  18. #include <boost/compute/algorithm/accumulate.hpp>
  19. #include <boost/compute/algorithm/copy.hpp>
  20. #include <boost/compute/algorithm/transform.hpp>
  21. #include <boost/compute/detail/literal.hpp>
  22. #include <boost/compute/types/fundamental.hpp>
  23. namespace boost {
  24. namespace compute {
  25. /// \class discrete_distribution
  26. /// \brief Produces random integers on the interval [0, n), where
  27. /// probability of each integer is given by the weight of the ith
  28. /// integer divided by the sum of all weights.
  29. ///
  30. /// The following example shows how to setup a discrete distribution to
  31. /// produce 0 and 1 with equal probability
  32. ///
  33. /// \snippet test/test_discrete_distribution.cpp generate
  34. ///
  35. template<class IntType = uint_>
  36. class discrete_distribution
  37. {
  38. public:
  39. typedef IntType result_type;
  40. /// Creates a new discrete distribution with a single weight p = { 1 }.
  41. /// This distribution produces only zeroes.
  42. discrete_distribution()
  43. : m_probabilities(1, double(1)),
  44. m_scanned_probabilities(1, double(1))
  45. {
  46. }
  47. /// Creates a new discrete distribution with weights given by
  48. /// the range [\p first, \p last).
  49. template<class InputIterator>
  50. discrete_distribution(InputIterator first, InputIterator last)
  51. : m_probabilities(first, last),
  52. m_scanned_probabilities(std::distance(first, last))
  53. {
  54. if(first != last) {
  55. // after this m_scanned_probabilities.back() is a sum of all
  56. // weights from the range [first, last)
  57. std::partial_sum(first, last, m_scanned_probabilities.begin());
  58. std::vector<double>::iterator i = m_probabilities.begin();
  59. std::vector<double>::iterator j = m_scanned_probabilities.begin();
  60. for(; i != m_probabilities.end(); ++i, ++j)
  61. {
  62. // dividing each weight by sum of all weights to
  63. // get probabilities
  64. *i = *i / m_scanned_probabilities.back();
  65. // dividing each partial sum of weights by sum of
  66. // all weights to get partial sums of probabilities
  67. *j = *j / m_scanned_probabilities.back();
  68. }
  69. }
  70. else {
  71. m_probabilities.push_back(double(1));
  72. m_scanned_probabilities.push_back(double(1));
  73. }
  74. }
  75. /// Destroys the discrete_distribution object.
  76. ~discrete_distribution()
  77. {
  78. }
  79. /// Returns the probabilities
  80. ::std::vector<double> probabilities() const
  81. {
  82. return m_probabilities;
  83. }
  84. /// Returns the minimum potentially generated value.
  85. result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const
  86. {
  87. return result_type(0);
  88. }
  89. /// Returns the maximum potentially generated value.
  90. result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const
  91. {
  92. size_t type_max = static_cast<size_t>(
  93. (std::numeric_limits<result_type>::max)()
  94. );
  95. if(m_probabilities.size() - 1 > type_max) {
  96. return (std::numeric_limits<result_type>::max)();
  97. }
  98. return static_cast<result_type>(m_probabilities.size() - 1);
  99. }
  100. /// Generates uniformly distributed integers and stores
  101. /// them to the range [\p first, \p last).
  102. template<class OutputIterator, class Generator>
  103. void generate(OutputIterator first,
  104. OutputIterator last,
  105. Generator &generator,
  106. command_queue &queue)
  107. {
  108. std::string source = "inline IntType scale_random(uint x)\n";
  109. source = source +
  110. "{\n" +
  111. "float rno = convert_float(x) / UINT_MAX;\n";
  112. for(size_t i = 0; i < m_scanned_probabilities.size() - 1; i++)
  113. {
  114. source = source +
  115. "if(rno <= " + detail::make_literal<float>(m_scanned_probabilities[i]) + ")\n" +
  116. " return " + detail::make_literal(i) + ";\n";
  117. }
  118. source = source +
  119. "return " + detail::make_literal(m_scanned_probabilities.size() - 1) + ";\n" +
  120. "}\n";
  121. BOOST_COMPUTE_FUNCTION(IntType, scale_random, (const uint_ x), {});
  122. scale_random.set_source(source);
  123. scale_random.define("IntType", type_name<IntType>());
  124. generator.generate(first, last, scale_random, queue);
  125. }
  126. private:
  127. ::std::vector<double> m_probabilities;
  128. ::std::vector<double> m_scanned_probabilities;
  129. BOOST_STATIC_ASSERT_MSG(
  130. boost::is_integral<IntType>::value,
  131. "Template argument must be integral"
  132. );
  133. };
  134. } // end compute namespace
  135. } // end boost namespace
  136. #endif // BOOST_COMPUTE_RANDOM_UNIFORM_INT_DISTRIBUTION_HPP