1 /* boost random/discrete_distribution.hpp header file 2 * 3 * Copyright Steven Watanabe 2009-2011 4 * Distributed under the Boost Software License, Version 1.0. (See 5 * accompanying file LICENSE_1_0.txt or copy at 6 * http://www.boost.org/LICENSE_1_0.txt) 7 * 8 * See http://www.boost.org for most recent version including documentation. 9 * 10 * $Id: discrete_distribution.hpp 71018 2011-04-05 21:27:52Z steven_watanabe $ 11 */ 12 13 #ifndef BOOST_RANDOM_DISCRETE_DISTRIBUTION_HPP_INCLUDED 14 #define BOOST_RANDOM_DISCRETE_DISTRIBUTION_HPP_INCLUDED 15 16 #include <vector> 17 #include <limits> 18 #include <numeric> 19 #include <utility> 20 #include <iterator> 21 #include <boost/assert.hpp> 22 #include <boost/random/uniform_01.hpp> 23 #include <boost/random/uniform_int.hpp> 24 #include <boost/random/detail/config.hpp> 25 #include <boost/random/detail/operators.hpp> 26 #include <boost/random/detail/vector_io.hpp> 27 28 #ifndef BOOST_NO_INITIALIZER_LISTS 29 #include <initializer_list> 30 #endif 31 32 #include <boost/range/begin.hpp> 33 #include <boost/range/end.hpp> 34 35 #include <boost/random/detail/disable_warnings.hpp> 36 37 namespace boost { 38 namespace random { 39 40 /** 41 * The class @c discrete_distribution models a \random_distribution. 42 * It produces integers in the range [0, n) with the probability 43 * of producing each value is specified by the parameters of the 44 * distribution. 45 */ 46 template<class IntType = int, class WeightType = double> 47 class discrete_distribution { 48 public: 49 typedef WeightType input_type; 50 typedef IntType result_type; 51 52 class param_type { 53 public: 54 55 typedef discrete_distribution distribution_type; 56 57 /** 58 * Constructs a @c param_type object, representing a distribution 59 * with \f$p(0) = 1\f$ and \f$p(k|k>0) = 0\f$. 60 */ param_type()61 param_type() : _probabilities(1, static_cast<WeightType>(1)) {} 62 /** 63 * If @c first == @c last, equivalent to the default constructor. 64 * Otherwise, the values of the range represent weights for the 65 * possible values of the distribution. 66 */ 67 template<class Iter> param_type(Iter first,Iter last)68 param_type(Iter first, Iter last) : _probabilities(first, last) 69 { 70 normalize(); 71 } 72 #ifndef BOOST_NO_INITIALIZER_LISTS 73 /** 74 * If wl.size() == 0, equivalent to the default constructor. 75 * Otherwise, the values of the @c initializer_list represent 76 * weights for the possible values of the distribution. 77 */ param_type(const std::initializer_list<WeightType> & wl)78 param_type(const std::initializer_list<WeightType>& wl) 79 : _probabilities(wl) 80 { 81 normalize(); 82 } 83 #endif 84 /** 85 * If the range is empty, equivalent to the default constructor. 86 * Otherwise, the elements of the range represent 87 * weights for the possible values of the distribution. 88 */ 89 template<class Range> param_type(const Range & range)90 explicit param_type(const Range& range) 91 : _probabilities(boost::begin(range), boost::end(range)) 92 { 93 normalize(); 94 } 95 96 /** 97 * If nw is zero, equivalent to the default constructor. 98 * Otherwise, the range of the distribution is [0, nw), 99 * and the weights are found by calling fw with values 100 * evenly distributed between \f$\mbox{xmin} + \delta/2\f$ and 101 * \f$\mbox{xmax} - \delta/2\f$, where 102 * \f$\delta = (\mbox{xmax} - \mbox{xmin})/\mbox{nw}\f$. 103 */ 104 template<class Func> param_type(std::size_t nw,double xmin,double xmax,Func fw)105 param_type(std::size_t nw, double xmin, double xmax, Func fw) 106 { 107 std::size_t n = (nw == 0) ? 1 : nw; 108 double delta = (xmax - xmin) / n; 109 BOOST_ASSERT(delta > 0); 110 for(std::size_t k = 0; k < n; ++k) { 111 _probabilities.push_back(fw(xmin + k*delta + delta/2)); 112 } 113 normalize(); 114 } 115 116 /** 117 * Returns a vector containing the probabilities of each possible 118 * value of the distribution. 119 */ probabilities() const120 std::vector<WeightType> probabilities() const 121 { 122 return _probabilities; 123 } 124 125 /** Writes the parameters to a @c std::ostream. */ BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os,param_type,parm)126 BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm) 127 { 128 detail::print_vector(os, parm._probabilities); 129 return os; 130 } 131 132 /** Reads the parameters from a @c std::istream. */ BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is,param_type,parm)133 BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm) 134 { 135 std::vector<WeightType> temp; 136 detail::read_vector(is, temp); 137 if(is) { 138 parm._probabilities.swap(temp); 139 } 140 return is; 141 } 142 143 /** Returns true if the two sets of parameters are the same. */ BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type,lhs,rhs)144 BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs) 145 { 146 return lhs._probabilities == rhs._probabilities; 147 } 148 /** Returns true if the two sets of parameters are different. */ 149 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type) 150 private: 151 /// @cond show_private 152 friend class discrete_distribution; param_type(const discrete_distribution & dist)153 explicit param_type(const discrete_distribution& dist) 154 : _probabilities(dist.probabilities()) 155 {} normalize()156 void normalize() 157 { 158 WeightType sum = 159 std::accumulate(_probabilities.begin(), _probabilities.end(), 160 static_cast<WeightType>(0)); 161 for(typename std::vector<WeightType>::iterator 162 iter = _probabilities.begin(), 163 end = _probabilities.end(); 164 iter != end; ++iter) 165 { 166 *iter /= sum; 167 } 168 } 169 std::vector<WeightType> _probabilities; 170 /// @endcond 171 }; 172 173 /** 174 * Creates a new @c discrete_distribution object that has 175 * \f$p(0) = 1\f$ and \f$p(i|i>0) = 0\f$. 176 */ discrete_distribution()177 discrete_distribution() 178 { 179 _alias_table.push_back(std::make_pair(static_cast<WeightType>(1), 180 static_cast<IntType>(0))); 181 } 182 /** 183 * Constructs a discrete_distribution from an iterator range. 184 * If @c first == @c last, equivalent to the default constructor. 185 * Otherwise, the values of the range represent weights for the 186 * possible values of the distribution. 187 */ 188 template<class Iter> discrete_distribution(Iter first,Iter last)189 discrete_distribution(Iter first, Iter last) 190 { 191 init(first, last); 192 } 193 #ifndef BOOST_NO_INITIALIZER_LISTS 194 /** 195 * Constructs a @c discrete_distribution from a @c std::initializer_list. 196 * If the @c initializer_list is empty, equivalent to the default 197 * constructor. Otherwise, the values of the @c initializer_list 198 * represent weights for the possible values of the distribution. 199 * For example, given the distribution 200 * 201 * @code 202 * discrete_distribution<> dist{1, 4, 5}; 203 * @endcode 204 * 205 * The probability of a 0 is 1/10, the probability of a 1 is 2/5, 206 * the probability of a 2 is 1/2, and no other values are possible. 207 */ discrete_distribution(std::initializer_list<WeightType> wl)208 discrete_distribution(std::initializer_list<WeightType> wl) 209 { 210 init(wl.begin(), wl.end()); 211 } 212 #endif 213 /** 214 * Constructs a discrete_distribution from a Boost.Range range. 215 * If the range is empty, equivalent to the default constructor. 216 * Otherwise, the values of the range represent weights for the 217 * possible values of the distribution. 218 */ 219 template<class Range> discrete_distribution(const Range & range)220 explicit discrete_distribution(const Range& range) 221 { 222 init(boost::begin(range), boost::end(range)); 223 } 224 /** 225 * Constructs a discrete_distribution that approximates a function. 226 * If nw is zero, equivalent to the default constructor. 227 * Otherwise, the range of the distribution is [0, nw), 228 * and the weights are found by calling fw with values 229 * evenly distributed between \f$\mbox{xmin} + \delta/2\f$ and 230 * \f$\mbox{xmax} - \delta/2\f$, where 231 * \f$\delta = (\mbox{xmax} - \mbox{xmin})/\mbox{nw}\f$. 232 */ 233 template<class Func> discrete_distribution(std::size_t nw,double xmin,double xmax,Func fw)234 discrete_distribution(std::size_t nw, double xmin, double xmax, Func fw) 235 { 236 std::size_t n = (nw == 0) ? 1 : nw; 237 double delta = (xmax - xmin) / n; 238 BOOST_ASSERT(delta > 0); 239 std::vector<WeightType> weights; 240 for(std::size_t k = 0; k < n; ++k) { 241 weights.push_back(fw(xmin + k*delta + delta/2)); 242 } 243 init(weights.begin(), weights.end()); 244 } 245 /** 246 * Constructs a discrete_distribution from its parameters. 247 */ discrete_distribution(const param_type & parm)248 explicit discrete_distribution(const param_type& parm) 249 { 250 param(parm); 251 } 252 253 /** 254 * Returns a value distributed according to the parameters of the 255 * discrete_distribution. 256 */ 257 template<class URNG> operator ()(URNG & urng) const258 IntType operator()(URNG& urng) const 259 { 260 BOOST_ASSERT(!_alias_table.empty()); 261 WeightType test = uniform_01<WeightType>()(urng); 262 IntType result = uniform_int<IntType>((min)(), (max)())(urng); 263 if(test < _alias_table[result].first) { 264 return result; 265 } else { 266 return(_alias_table[result].second); 267 } 268 } 269 270 /** 271 * Returns a value distributed according to the parameters 272 * specified by param. 273 */ 274 template<class URNG> operator ()(URNG & urng,const param_type & parm) const275 IntType operator()(URNG& urng, const param_type& parm) const 276 { 277 while(true) { 278 WeightType val = uniform_01<WeightType>()(urng); 279 WeightType sum = 0; 280 std::size_t result = 0; 281 for(typename std::vector<WeightType>::const_iterator 282 iter = parm._probabilities.begin(), 283 end = parm._probabilities.end(); 284 iter != end; ++iter, ++result) 285 { 286 sum += *iter; 287 if(sum > val) { 288 return result; 289 } 290 } 291 } 292 } 293 294 /** Returns the smallest value that the distribution can produce. */ BOOST_PREVENT_MACRO_SUBSTITUTION() const295 result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; } 296 /** Returns the largest value that the distribution can produce. */ BOOST_PREVENT_MACRO_SUBSTITUTION() const297 result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const 298 { return static_cast<result_type>(_alias_table.size() - 1); } 299 300 /** 301 * Returns a vector containing the probabilities of each 302 * value of the distribution. For example, given 303 * 304 * @code 305 * discrete_distribution<> dist = { 1, 4, 5 }; 306 * std::vector<double> p = dist.param(); 307 * @endcode 308 * 309 * the vector, p will contain {0.1, 0.4, 0.5}. 310 */ probabilities() const311 std::vector<WeightType> probabilities() const 312 { 313 std::vector<WeightType> result(_alias_table.size()); 314 const WeightType mean = 315 static_cast<WeightType>(1) / _alias_table.size(); 316 std::size_t i = 0; 317 for(typename alias_table_t::const_iterator 318 iter = _alias_table.begin(), 319 end = _alias_table.end(); 320 iter != end; ++iter, ++i) 321 { 322 WeightType val = iter->first * mean; 323 result[i] += val; 324 result[iter->second] += mean - val; 325 } 326 return(result); 327 } 328 329 /** Returns the parameters of the distribution. */ param() const330 param_type param() const 331 { 332 return param_type(*this); 333 } 334 /** Sets the parameters of the distribution. */ param(const param_type & parm)335 void param(const param_type& parm) 336 { 337 init(parm._probabilities.begin(), parm._probabilities.end()); 338 } 339 340 /** 341 * Effects: Subsequent uses of the distribution do not depend 342 * on values produced by any engine prior to invoking reset. 343 */ reset()344 void reset() {} 345 346 /** Writes a distribution to a @c std::ostream. */ BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os,discrete_distribution,dd)347 BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, discrete_distribution, dd) 348 { 349 os << dd.param(); 350 return os; 351 } 352 353 /** Reads a distribution from a @c std::istream */ BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is,discrete_distribution,dd)354 BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, discrete_distribution, dd) 355 { 356 param_type parm; 357 if(is >> parm) { 358 dd.param(parm); 359 } 360 return is; 361 } 362 363 /** 364 * Returns true if the two distributions will return the 365 * same sequence of values, when passed equal generators. 366 */ BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(discrete_distribution,lhs,rhs)367 BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(discrete_distribution, lhs, rhs) 368 { 369 return lhs._alias_table == rhs._alias_table; 370 } 371 /** 372 * Returns true if the two distributions may return different 373 * sequences of values, when passed equal generators. 374 */ 375 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(discrete_distribution) 376 377 private: 378 379 /// @cond show_private 380 381 template<class Iter> init(Iter first,Iter last,std::input_iterator_tag)382 void init(Iter first, Iter last, std::input_iterator_tag) 383 { 384 std::vector<WeightType> temp(first, last); 385 init(temp.begin(), temp.end()); 386 } 387 template<class Iter> init(Iter first,Iter last,std::forward_iterator_tag)388 void init(Iter first, Iter last, std::forward_iterator_tag) 389 { 390 std::vector<std::pair<WeightType, IntType> > below_average; 391 std::vector<std::pair<WeightType, IntType> > above_average; 392 std::size_t size = std::distance(first, last); 393 WeightType weight_sum = 394 std::accumulate(first, last, static_cast<WeightType>(0)); 395 WeightType weight_average = weight_sum / size; 396 std::size_t i = 0; 397 for(; first != last; ++first, ++i) { 398 WeightType val = *first / weight_average; 399 std::pair<WeightType, IntType> elem(val, static_cast<IntType>(i)); 400 if(val < static_cast<WeightType>(1)) { 401 below_average.push_back(elem); 402 } else { 403 above_average.push_back(elem); 404 } 405 } 406 407 _alias_table.resize(size); 408 typename alias_table_t::iterator 409 b_iter = below_average.begin(), 410 b_end = below_average.end(), 411 a_iter = above_average.begin(), 412 a_end = above_average.end() 413 ; 414 while(b_iter != b_end && a_iter != a_end) { 415 _alias_table[b_iter->second] = 416 std::make_pair(b_iter->first, a_iter->second); 417 a_iter->first -= (static_cast<WeightType>(1) - b_iter->first); 418 if(a_iter->first < static_cast<WeightType>(1)) { 419 *b_iter = *a_iter++; 420 } else { 421 ++b_iter; 422 } 423 } 424 for(; b_iter != b_end; ++b_iter) { 425 _alias_table[b_iter->second].first = static_cast<WeightType>(1); 426 } 427 for(; a_iter != a_end; ++a_iter) { 428 _alias_table[a_iter->second].first = static_cast<WeightType>(1); 429 } 430 } 431 template<class Iter> init(Iter first,Iter last)432 void init(Iter first, Iter last) 433 { 434 if(first == last) { 435 _alias_table.clear(); 436 _alias_table.push_back(std::make_pair(static_cast<WeightType>(1), 437 static_cast<IntType>(0))); 438 } else { 439 typename std::iterator_traits<Iter>::iterator_category category; 440 init(first, last, category); 441 } 442 } 443 typedef std::vector<std::pair<WeightType, IntType> > alias_table_t; 444 alias_table_t _alias_table; 445 /// @endcond 446 }; 447 448 } 449 } 450 451 #include <boost/random/detail/enable_warnings.hpp> 452 453 #endif 454