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