1 /* boost random/discrete_distribution.hpp header file
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)
8 * See http://www.boost.org for most recent version including documentation.
13 #ifndef BOOST_RANDOM_DISCRETE_DISTRIBUTION_HPP_INCLUDED
14 #define BOOST_RANDOM_DISCRETE_DISTRIBUTION_HPP_INCLUDED
21 #include <boost/assert.hpp>
22 #include <boost/random/uniform_01.hpp>
23 #include <boost/random/uniform_int_distribution.hpp>
24 #include <boost/random/detail/config.hpp>
25 #include <boost/random/detail/operators.hpp>
26 #include <boost/random/detail/vector_io.hpp>
28 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
29 #include <initializer_list>
32 #include <boost/range/begin.hpp>
33 #include <boost/range/end.hpp>
35 #include <boost/random/detail/disable_warnings.hpp>
41 template<class IntType, class WeightType>
42 struct integer_alias_table {
43 WeightType get_weight(IntType bin) const {
44 WeightType result = _average;
45 if(bin < _excess) ++result;
49 WeightType init_average(Iter begin, Iter end) {
50 WeightType weight_average = 0;
53 // weight_average * n + excess == current partial sum
54 // This is a bit messy, but it's guaranteed not to overflow
55 for(Iter iter = begin; iter != end; ++iter) {
57 if(*iter < weight_average) {
58 WeightType diff = weight_average - *iter;
59 weight_average -= diff / n;
60 if(diff % n > excess) {
62 excess += n - diff % n;
67 WeightType diff = *iter - weight_average;
68 weight_average += diff / n;
69 if(diff % n < n - excess) {
73 excess -= n - diff % n;
77 _alias_table.resize(static_cast<std::size_t>(n));
78 _average = weight_average;
80 return weight_average;
85 _alias_table.push_back(std::make_pair(static_cast<WeightType>(1),
86 static_cast<IntType>(0)));
87 _average = static_cast<WeightType>(1);
88 _excess = static_cast<IntType>(0);
90 bool operator==(const integer_alias_table& other) const
92 return _alias_table == other._alias_table &&
93 _average == other._average && _excess == other._excess;
95 static WeightType normalize(WeightType val, WeightType average)
99 static void normalize(std::vector<WeightType>&) {}
101 WeightType test(URNG &urng) const
103 return uniform_int_distribution<WeightType>(0, _average)(urng);
105 bool accept(IntType result, WeightType val) const
107 return result < _excess || val < _average;
109 static WeightType try_get_sum(const std::vector<WeightType>& weights)
111 WeightType result = static_cast<WeightType>(0);
112 for(typename std::vector<WeightType>::const_iterator
113 iter = weights.begin(), end = weights.end();
116 if((std::numeric_limits<WeightType>::max)() - result > *iter) {
117 return static_cast<WeightType>(0);
124 static WeightType generate_in_range(URNG &urng, WeightType max)
126 return uniform_int_distribution<WeightType>(
127 static_cast<WeightType>(0), max-1)(urng);
129 typedef std::vector<std::pair<WeightType, IntType> > alias_table_t;
130 alias_table_t _alias_table;
135 template<class IntType, class WeightType>
136 struct real_alias_table {
137 WeightType get_weight(IntType) const
139 return WeightType(1.0);
142 WeightType init_average(Iter first, Iter last)
144 std::size_t size = std::distance(first, last);
145 WeightType weight_sum =
146 std::accumulate(first, last, static_cast<WeightType>(0));
147 _alias_table.resize(size);
148 return weight_sum / size;
152 _alias_table.clear();
153 _alias_table.push_back(std::make_pair(static_cast<WeightType>(1),
154 static_cast<IntType>(0)));
156 bool operator==(const real_alias_table& other) const
158 return _alias_table == other._alias_table;
160 static WeightType normalize(WeightType val, WeightType average)
162 return val / average;
164 static void normalize(std::vector<WeightType>& weights)
167 std::accumulate(weights.begin(), weights.end(),
168 static_cast<WeightType>(0));
169 for(typename std::vector<WeightType>::iterator
170 iter = weights.begin(),
178 WeightType test(URNG &urng) const
180 return uniform_01<WeightType>()(urng);
182 bool accept(IntType, WeightType) const
186 static WeightType try_get_sum(const std::vector<WeightType>& weights)
188 return static_cast<WeightType>(1);
191 static WeightType generate_in_range(URNG &urng, WeightType)
193 return uniform_01<WeightType>()(urng);
195 typedef std::vector<std::pair<WeightType, IntType> > alias_table_t;
196 alias_table_t _alias_table;
199 template<bool IsIntegral>
200 struct select_alias_table;
203 struct select_alias_table<true> {
204 template<class IntType, class WeightType>
206 typedef integer_alias_table<IntType, WeightType> type;
211 struct select_alias_table<false> {
212 template<class IntType, class WeightType>
214 typedef real_alias_table<IntType, WeightType> type;
221 * The class @c discrete_distribution models a \random_distribution.
222 * It produces integers in the range [0, n) with the probability
223 * of producing each value is specified by the parameters of the
226 template<class IntType = int, class WeightType = double>
227 class discrete_distribution {
229 typedef WeightType input_type;
230 typedef IntType result_type;
235 typedef discrete_distribution distribution_type;
238 * Constructs a @c param_type object, representing a distribution
239 * with \f$p(0) = 1\f$ and \f$p(k|k>0) = 0\f$.
241 param_type() : _probabilities(1, static_cast<WeightType>(1)) {}
243 * If @c first == @c last, equivalent to the default constructor.
244 * Otherwise, the values of the range represent weights for the
245 * possible values of the distribution.
248 param_type(Iter first, Iter last) : _probabilities(first, last)
252 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
254 * If wl.size() == 0, equivalent to the default constructor.
255 * Otherwise, the values of the @c initializer_list represent
256 * weights for the possible values of the distribution.
258 param_type(const std::initializer_list<WeightType>& wl)
265 * If the range is empty, equivalent to the default constructor.
266 * Otherwise, the elements of the range represent
267 * weights for the possible values of the distribution.
269 template<class Range>
270 explicit param_type(const Range& range)
271 : _probabilities(boost::begin(range), boost::end(range))
277 * If nw is zero, equivalent to the default constructor.
278 * Otherwise, the range of the distribution is [0, nw),
279 * and the weights are found by calling fw with values
280 * evenly distributed between \f$\mbox{xmin} + \delta/2\f$ and
281 * \f$\mbox{xmax} - \delta/2\f$, where
282 * \f$\delta = (\mbox{xmax} - \mbox{xmin})/\mbox{nw}\f$.
285 param_type(std::size_t nw, double xmin, double xmax, Func fw)
287 std::size_t n = (nw == 0) ? 1 : nw;
288 double delta = (xmax - xmin) / n;
289 BOOST_ASSERT(delta > 0);
290 for(std::size_t k = 0; k < n; ++k) {
291 _probabilities.push_back(fw(xmin + k*delta + delta/2));
297 * Returns a vector containing the probabilities of each possible
298 * value of the distribution.
300 std::vector<WeightType> probabilities() const
302 return _probabilities;
305 /** Writes the parameters to a @c std::ostream. */
306 BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
308 detail::print_vector(os, parm._probabilities);
312 /** Reads the parameters from a @c std::istream. */
313 BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
315 std::vector<WeightType> temp;
316 detail::read_vector(is, temp);
318 parm._probabilities.swap(temp);
323 /** Returns true if the two sets of parameters are the same. */
324 BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
326 return lhs._probabilities == rhs._probabilities;
328 /** Returns true if the two sets of parameters are different. */
329 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
331 /// @cond show_private
332 friend class discrete_distribution;
333 explicit param_type(const discrete_distribution& dist)
334 : _probabilities(dist.probabilities())
338 impl_type::normalize(_probabilities);
340 std::vector<WeightType> _probabilities;
345 * Creates a new @c discrete_distribution object that has
346 * \f$p(0) = 1\f$ and \f$p(i|i>0) = 0\f$.
348 discrete_distribution()
353 * Constructs a discrete_distribution from an iterator range.
354 * If @c first == @c last, equivalent to the default constructor.
355 * Otherwise, the values of the range represent weights for the
356 * possible values of the distribution.
359 discrete_distribution(Iter first, Iter last)
363 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
365 * Constructs a @c discrete_distribution from a @c std::initializer_list.
366 * If the @c initializer_list is empty, equivalent to the default
367 * constructor. Otherwise, the values of the @c initializer_list
368 * represent weights for the possible values of the distribution.
369 * For example, given the distribution
372 * discrete_distribution<> dist{1, 4, 5};
375 * The probability of a 0 is 1/10, the probability of a 1 is 2/5,
376 * the probability of a 2 is 1/2, and no other values are possible.
378 discrete_distribution(std::initializer_list<WeightType> wl)
380 init(wl.begin(), wl.end());
384 * Constructs a discrete_distribution from a Boost.Range range.
385 * If the range is empty, equivalent to the default constructor.
386 * Otherwise, the values of the range represent weights for the
387 * possible values of the distribution.
389 template<class Range>
390 explicit discrete_distribution(const Range& range)
392 init(boost::begin(range), boost::end(range));
395 * Constructs a discrete_distribution that approximates a function.
396 * If nw is zero, equivalent to the default constructor.
397 * Otherwise, the range of the distribution is [0, nw),
398 * and the weights are found by calling fw with values
399 * evenly distributed between \f$\mbox{xmin} + \delta/2\f$ and
400 * \f$\mbox{xmax} - \delta/2\f$, where
401 * \f$\delta = (\mbox{xmax} - \mbox{xmin})/\mbox{nw}\f$.
404 discrete_distribution(std::size_t nw, double xmin, double xmax, Func fw)
406 std::size_t n = (nw == 0) ? 1 : nw;
407 double delta = (xmax - xmin) / n;
408 BOOST_ASSERT(delta > 0);
409 std::vector<WeightType> weights;
410 for(std::size_t k = 0; k < n; ++k) {
411 weights.push_back(fw(xmin + k*delta + delta/2));
413 init(weights.begin(), weights.end());
416 * Constructs a discrete_distribution from its parameters.
418 explicit discrete_distribution(const param_type& parm)
424 * Returns a value distributed according to the parameters of the
425 * discrete_distribution.
428 IntType operator()(URNG& urng) const
430 BOOST_ASSERT(!_impl._alias_table.empty());
434 result = uniform_int_distribution<IntType>((min)(), (max)())(urng);
435 test = _impl.test(urng);
436 } while(!_impl.accept(result, test));
437 if(test < _impl._alias_table[static_cast<std::size_t>(result)].first) {
440 return(_impl._alias_table[static_cast<std::size_t>(result)].second);
445 * Returns a value distributed according to the parameters
446 * specified by param.
449 IntType operator()(URNG& urng, const param_type& parm) const
451 if(WeightType limit = impl_type::try_get_sum(parm._probabilities)) {
452 WeightType val = impl_type::generate_in_range(urng, limit);
454 std::size_t result = 0;
455 for(typename std::vector<WeightType>::const_iterator
456 iter = parm._probabilities.begin(),
457 end = parm._probabilities.end();
458 iter != end; ++iter, ++result)
465 // This shouldn't be reachable, but round-off error
466 // can prevent any match from being found when val is
468 return static_cast<IntType>(parm._probabilities.size() - 1);
470 // WeightType is integral and sum(parm._probabilities)
471 // would overflow. Just use the easy solution.
472 return discrete_distribution(parm)(urng);
476 /** Returns the smallest value that the distribution can produce. */
477 result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; }
478 /** Returns the largest value that the distribution can produce. */
479 result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const
480 { return static_cast<result_type>(_impl._alias_table.size() - 1); }
483 * Returns a vector containing the probabilities of each
484 * value of the distribution. For example, given
487 * discrete_distribution<> dist = { 1, 4, 5 };
488 * std::vector<double> p = dist.param();
491 * the vector, p will contain {0.1, 0.4, 0.5}.
493 * If @c WeightType is integral, then the weights
494 * will be returned unchanged.
496 std::vector<WeightType> probabilities() const
498 std::vector<WeightType> result(_impl._alias_table.size(), static_cast<WeightType>(0));
500 for(typename impl_type::alias_table_t::const_iterator
501 iter = _impl._alias_table.begin(),
502 end = _impl._alias_table.end();
503 iter != end; ++iter, ++i)
505 WeightType val = iter->first;
507 result[static_cast<std::size_t>(iter->second)] += _impl.get_weight(i) - val;
509 impl_type::normalize(result);
513 /** Returns the parameters of the distribution. */
514 param_type param() const
516 return param_type(*this);
518 /** Sets the parameters of the distribution. */
519 void param(const param_type& parm)
521 init(parm._probabilities.begin(), parm._probabilities.end());
525 * Effects: Subsequent uses of the distribution do not depend
526 * on values produced by any engine prior to invoking reset.
530 /** Writes a distribution to a @c std::ostream. */
531 BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, discrete_distribution, dd)
537 /** Reads a distribution from a @c std::istream */
538 BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, discrete_distribution, dd)
548 * Returns true if the two distributions will return the
549 * same sequence of values, when passed equal generators.
551 BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(discrete_distribution, lhs, rhs)
553 return lhs._impl == rhs._impl;
556 * Returns true if the two distributions may return different
557 * sequences of values, when passed equal generators.
559 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(discrete_distribution)
563 /// @cond show_private
566 void init(Iter first, Iter last, std::input_iterator_tag)
568 std::vector<WeightType> temp(first, last);
569 init(temp.begin(), temp.end());
572 void init(Iter first, Iter last, std::forward_iterator_tag)
574 std::vector<std::pair<WeightType, IntType> > below_average;
575 std::vector<std::pair<WeightType, IntType> > above_average;
576 WeightType weight_average = _impl.init_average(first, last);
577 WeightType normalized_average = _impl.get_weight(0);
579 for(; first != last; ++first, ++i) {
580 WeightType val = impl_type::normalize(*first, weight_average);
581 std::pair<WeightType, IntType> elem(val, static_cast<IntType>(i));
582 if(val < normalized_average) {
583 below_average.push_back(elem);
585 above_average.push_back(elem);
589 typename impl_type::alias_table_t::iterator
590 b_iter = below_average.begin(),
591 b_end = below_average.end(),
592 a_iter = above_average.begin(),
593 a_end = above_average.end()
595 while(b_iter != b_end && a_iter != a_end) {
596 _impl._alias_table[static_cast<std::size_t>(b_iter->second)] =
597 std::make_pair(b_iter->first, a_iter->second);
598 a_iter->first -= (_impl.get_weight(b_iter->second) - b_iter->first);
599 if(a_iter->first < normalized_average) {
605 for(; b_iter != b_end; ++b_iter) {
606 _impl._alias_table[static_cast<std::size_t>(b_iter->second)].first =
607 _impl.get_weight(b_iter->second);
609 for(; a_iter != a_end; ++a_iter) {
610 _impl._alias_table[static_cast<std::size_t>(a_iter->second)].first =
611 _impl.get_weight(a_iter->second);
615 void init(Iter first, Iter last)
620 typename std::iterator_traits<Iter>::iterator_category category;
621 init(first, last, category);
624 typedef typename detail::select_alias_table<
625 (::boost::is_integral<WeightType>::value)
626 >::template apply<IntType, WeightType>::type impl_type;
634 #include <boost/random/detail/enable_warnings.hpp>