| /* boost random/discrete_distribution.hpp header file |
| * |
| * Copyright Steven Watanabe 2009-2011 |
| * 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) |
| * |
| * See http://www.boost.org for most recent version including documentation. |
| * |
| * $Id$ |
| */ |
| |
| #ifndef BOOST_RANDOM_DISCRETE_DISTRIBUTION_HPP_INCLUDED |
| #define BOOST_RANDOM_DISCRETE_DISTRIBUTION_HPP_INCLUDED |
| |
| #include <vector> |
| #include <limits> |
| #include <numeric> |
| #include <utility> |
| #include <iterator> |
| #include <boost/assert.hpp> |
| #include <boost/random/uniform_01.hpp> |
| #include <boost/random/uniform_int_distribution.hpp> |
| #include <boost/random/detail/config.hpp> |
| #include <boost/random/detail/operators.hpp> |
| #include <boost/random/detail/vector_io.hpp> |
| |
| #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST |
| #include <initializer_list> |
| #endif |
| |
| #include <boost/range/begin.hpp> |
| #include <boost/range/end.hpp> |
| |
| #include <boost/random/detail/disable_warnings.hpp> |
| |
| namespace boost { |
| namespace random { |
| namespace detail { |
| |
| template<class IntType, class WeightType> |
| struct integer_alias_table { |
| WeightType get_weight(IntType bin) const { |
| WeightType result = _average; |
| if(bin < _excess) ++result; |
| return result; |
| } |
| template<class Iter> |
| WeightType init_average(Iter begin, Iter end) { |
| WeightType weight_average = 0; |
| IntType excess = 0; |
| IntType n = 0; |
| // weight_average * n + excess == current partial sum |
| // This is a bit messy, but it's guaranteed not to overflow |
| for(Iter iter = begin; iter != end; ++iter) { |
| ++n; |
| if(*iter < weight_average) { |
| WeightType diff = weight_average - *iter; |
| weight_average -= diff / n; |
| if(diff % n > excess) { |
| --weight_average; |
| excess += n - diff % n; |
| } else { |
| excess -= diff % n; |
| } |
| } else { |
| WeightType diff = *iter - weight_average; |
| weight_average += diff / n; |
| if(diff % n < n - excess) { |
| excess += diff % n; |
| } else { |
| ++weight_average; |
| excess -= n - diff % n; |
| } |
| } |
| } |
| _alias_table.resize(static_cast<std::size_t>(n)); |
| _average = weight_average; |
| _excess = excess; |
| return weight_average; |
| } |
| void init_empty() |
| { |
| _alias_table.clear(); |
| _alias_table.push_back(std::make_pair(static_cast<WeightType>(1), |
| static_cast<IntType>(0))); |
| _average = static_cast<WeightType>(1); |
| _excess = static_cast<IntType>(0); |
| } |
| bool operator==(const integer_alias_table& other) const |
| { |
| return _alias_table == other._alias_table && |
| _average == other._average && _excess == other._excess; |
| } |
| static WeightType normalize(WeightType val, WeightType average) |
| { |
| return val; |
| } |
| static void normalize(std::vector<WeightType>&) {} |
| template<class URNG> |
| WeightType test(URNG &urng) const |
| { |
| return uniform_int_distribution<WeightType>(0, _average)(urng); |
| } |
| bool accept(IntType result, WeightType val) const |
| { |
| return result < _excess || val < _average; |
| } |
| static WeightType try_get_sum(const std::vector<WeightType>& weights) |
| { |
| WeightType result = static_cast<WeightType>(0); |
| for(typename std::vector<WeightType>::const_iterator |
| iter = weights.begin(), end = weights.end(); |
| iter != end; ++iter) |
| { |
| if((std::numeric_limits<WeightType>::max)() - result > *iter) { |
| return static_cast<WeightType>(0); |
| } |
| result += *iter; |
| } |
| return result; |
| } |
| template<class URNG> |
| static WeightType generate_in_range(URNG &urng, WeightType max) |
| { |
| return uniform_int_distribution<WeightType>( |
| static_cast<WeightType>(0), max-1)(urng); |
| } |
| typedef std::vector<std::pair<WeightType, IntType> > alias_table_t; |
| alias_table_t _alias_table; |
| WeightType _average; |
| IntType _excess; |
| }; |
| |
| template<class IntType, class WeightType> |
| struct real_alias_table { |
| WeightType get_weight(IntType) const |
| { |
| return WeightType(1.0); |
| } |
| template<class Iter> |
| WeightType init_average(Iter first, Iter last) |
| { |
| std::size_t size = std::distance(first, last); |
| WeightType weight_sum = |
| std::accumulate(first, last, static_cast<WeightType>(0)); |
| _alias_table.resize(size); |
| return weight_sum / size; |
| } |
| void init_empty() |
| { |
| _alias_table.clear(); |
| _alias_table.push_back(std::make_pair(static_cast<WeightType>(1), |
| static_cast<IntType>(0))); |
| } |
| bool operator==(const real_alias_table& other) const |
| { |
| return _alias_table == other._alias_table; |
| } |
| static WeightType normalize(WeightType val, WeightType average) |
| { |
| return val / average; |
| } |
| static void normalize(std::vector<WeightType>& weights) |
| { |
| WeightType sum = |
| std::accumulate(weights.begin(), weights.end(), |
| static_cast<WeightType>(0)); |
| for(typename std::vector<WeightType>::iterator |
| iter = weights.begin(), |
| end = weights.end(); |
| iter != end; ++iter) |
| { |
| *iter /= sum; |
| } |
| } |
| template<class URNG> |
| WeightType test(URNG &urng) const |
| { |
| return uniform_01<WeightType>()(urng); |
| } |
| bool accept(IntType, WeightType) const |
| { |
| return true; |
| } |
| static WeightType try_get_sum(const std::vector<WeightType>& weights) |
| { |
| return static_cast<WeightType>(1); |
| } |
| template<class URNG> |
| static WeightType generate_in_range(URNG &urng, WeightType) |
| { |
| return uniform_01<WeightType>()(urng); |
| } |
| typedef std::vector<std::pair<WeightType, IntType> > alias_table_t; |
| alias_table_t _alias_table; |
| }; |
| |
| template<bool IsIntegral> |
| struct select_alias_table; |
| |
| template<> |
| struct select_alias_table<true> { |
| template<class IntType, class WeightType> |
| struct apply { |
| typedef integer_alias_table<IntType, WeightType> type; |
| }; |
| }; |
| |
| template<> |
| struct select_alias_table<false> { |
| template<class IntType, class WeightType> |
| struct apply { |
| typedef real_alias_table<IntType, WeightType> type; |
| }; |
| }; |
| |
| } |
| |
| /** |
| * The class @c discrete_distribution models a \random_distribution. |
| * It produces integers in the range [0, n) with the probability |
| * of producing each value is specified by the parameters of the |
| * distribution. |
| */ |
| template<class IntType = int, class WeightType = double> |
| class discrete_distribution { |
| public: |
| typedef WeightType input_type; |
| typedef IntType result_type; |
| |
| class param_type { |
| public: |
| |
| typedef discrete_distribution distribution_type; |
| |
| /** |
| * Constructs a @c param_type object, representing a distribution |
| * with \f$p(0) = 1\f$ and \f$p(k|k>0) = 0\f$. |
| */ |
| param_type() : _probabilities(1, static_cast<WeightType>(1)) {} |
| /** |
| * If @c first == @c last, equivalent to the default constructor. |
| * Otherwise, the values of the range represent weights for the |
| * possible values of the distribution. |
| */ |
| template<class Iter> |
| param_type(Iter first, Iter last) : _probabilities(first, last) |
| { |
| normalize(); |
| } |
| #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST |
| /** |
| * If wl.size() == 0, equivalent to the default constructor. |
| * Otherwise, the values of the @c initializer_list represent |
| * weights for the possible values of the distribution. |
| */ |
| param_type(const std::initializer_list<WeightType>& wl) |
| : _probabilities(wl) |
| { |
| normalize(); |
| } |
| #endif |
| /** |
| * If the range is empty, equivalent to the default constructor. |
| * Otherwise, the elements of the range represent |
| * weights for the possible values of the distribution. |
| */ |
| template<class Range> |
| explicit param_type(const Range& range) |
| : _probabilities(boost::begin(range), boost::end(range)) |
| { |
| normalize(); |
| } |
| |
| /** |
| * If nw is zero, equivalent to the default constructor. |
| * Otherwise, the range of the distribution is [0, nw), |
| * and the weights are found by calling fw with values |
| * evenly distributed between \f$\mbox{xmin} + \delta/2\f$ and |
| * \f$\mbox{xmax} - \delta/2\f$, where |
| * \f$\delta = (\mbox{xmax} - \mbox{xmin})/\mbox{nw}\f$. |
| */ |
| template<class Func> |
| param_type(std::size_t nw, double xmin, double xmax, Func fw) |
| { |
| std::size_t n = (nw == 0) ? 1 : nw; |
| double delta = (xmax - xmin) / n; |
| BOOST_ASSERT(delta > 0); |
| for(std::size_t k = 0; k < n; ++k) { |
| _probabilities.push_back(fw(xmin + k*delta + delta/2)); |
| } |
| normalize(); |
| } |
| |
| /** |
| * Returns a vector containing the probabilities of each possible |
| * value of the distribution. |
| */ |
| std::vector<WeightType> probabilities() const |
| { |
| return _probabilities; |
| } |
| |
| /** Writes the parameters to a @c std::ostream. */ |
| BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm) |
| { |
| detail::print_vector(os, parm._probabilities); |
| return os; |
| } |
| |
| /** Reads the parameters from a @c std::istream. */ |
| BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm) |
| { |
| std::vector<WeightType> temp; |
| detail::read_vector(is, temp); |
| if(is) { |
| parm._probabilities.swap(temp); |
| } |
| return is; |
| } |
| |
| /** Returns true if the two sets of parameters are the same. */ |
| BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs) |
| { |
| return lhs._probabilities == rhs._probabilities; |
| } |
| /** Returns true if the two sets of parameters are different. */ |
| BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type) |
| private: |
| /// @cond show_private |
| friend class discrete_distribution; |
| explicit param_type(const discrete_distribution& dist) |
| : _probabilities(dist.probabilities()) |
| {} |
| void normalize() |
| { |
| impl_type::normalize(_probabilities); |
| } |
| std::vector<WeightType> _probabilities; |
| /// @endcond |
| }; |
| |
| /** |
| * Creates a new @c discrete_distribution object that has |
| * \f$p(0) = 1\f$ and \f$p(i|i>0) = 0\f$. |
| */ |
| discrete_distribution() |
| { |
| _impl.init_empty(); |
| } |
| /** |
| * Constructs a discrete_distribution from an iterator range. |
| * If @c first == @c last, equivalent to the default constructor. |
| * Otherwise, the values of the range represent weights for the |
| * possible values of the distribution. |
| */ |
| template<class Iter> |
| discrete_distribution(Iter first, Iter last) |
| { |
| init(first, last); |
| } |
| #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST |
| /** |
| * Constructs a @c discrete_distribution from a @c std::initializer_list. |
| * If the @c initializer_list is empty, equivalent to the default |
| * constructor. Otherwise, the values of the @c initializer_list |
| * represent weights for the possible values of the distribution. |
| * For example, given the distribution |
| * |
| * @code |
| * discrete_distribution<> dist{1, 4, 5}; |
| * @endcode |
| * |
| * The probability of a 0 is 1/10, the probability of a 1 is 2/5, |
| * the probability of a 2 is 1/2, and no other values are possible. |
| */ |
| discrete_distribution(std::initializer_list<WeightType> wl) |
| { |
| init(wl.begin(), wl.end()); |
| } |
| #endif |
| /** |
| * Constructs a discrete_distribution from a Boost.Range range. |
| * If the range is empty, equivalent to the default constructor. |
| * Otherwise, the values of the range represent weights for the |
| * possible values of the distribution. |
| */ |
| template<class Range> |
| explicit discrete_distribution(const Range& range) |
| { |
| init(boost::begin(range), boost::end(range)); |
| } |
| /** |
| * Constructs a discrete_distribution that approximates a function. |
| * If nw is zero, equivalent to the default constructor. |
| * Otherwise, the range of the distribution is [0, nw), |
| * and the weights are found by calling fw with values |
| * evenly distributed between \f$\mbox{xmin} + \delta/2\f$ and |
| * \f$\mbox{xmax} - \delta/2\f$, where |
| * \f$\delta = (\mbox{xmax} - \mbox{xmin})/\mbox{nw}\f$. |
| */ |
| template<class Func> |
| discrete_distribution(std::size_t nw, double xmin, double xmax, Func fw) |
| { |
| std::size_t n = (nw == 0) ? 1 : nw; |
| double delta = (xmax - xmin) / n; |
| BOOST_ASSERT(delta > 0); |
| std::vector<WeightType> weights; |
| for(std::size_t k = 0; k < n; ++k) { |
| weights.push_back(fw(xmin + k*delta + delta/2)); |
| } |
| init(weights.begin(), weights.end()); |
| } |
| /** |
| * Constructs a discrete_distribution from its parameters. |
| */ |
| explicit discrete_distribution(const param_type& parm) |
| { |
| param(parm); |
| } |
| |
| /** |
| * Returns a value distributed according to the parameters of the |
| * discrete_distribution. |
| */ |
| template<class URNG> |
| IntType operator()(URNG& urng) const |
| { |
| BOOST_ASSERT(!_impl._alias_table.empty()); |
| IntType result; |
| WeightType test; |
| do { |
| result = uniform_int_distribution<IntType>((min)(), (max)())(urng); |
| test = _impl.test(urng); |
| } while(!_impl.accept(result, test)); |
| if(test < _impl._alias_table[static_cast<std::size_t>(result)].first) { |
| return result; |
| } else { |
| return(_impl._alias_table[static_cast<std::size_t>(result)].second); |
| } |
| } |
| |
| /** |
| * Returns a value distributed according to the parameters |
| * specified by param. |
| */ |
| template<class URNG> |
| IntType operator()(URNG& urng, const param_type& parm) const |
| { |
| if(WeightType limit = impl_type::try_get_sum(parm._probabilities)) { |
| WeightType val = impl_type::generate_in_range(urng, limit); |
| WeightType sum = 0; |
| std::size_t result = 0; |
| for(typename std::vector<WeightType>::const_iterator |
| iter = parm._probabilities.begin(), |
| end = parm._probabilities.end(); |
| iter != end; ++iter, ++result) |
| { |
| sum += *iter; |
| if(sum > val) { |
| return result; |
| } |
| } |
| // This shouldn't be reachable, but round-off error |
| // can prevent any match from being found when val is |
| // very close to 1. |
| return static_cast<IntType>(parm._probabilities.size() - 1); |
| } else { |
| // WeightType is integral and sum(parm._probabilities) |
| // would overflow. Just use the easy solution. |
| return discrete_distribution(parm)(urng); |
| } |
| } |
| |
| /** Returns the smallest value that the distribution can produce. */ |
| result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; } |
| /** Returns the largest value that the distribution can produce. */ |
| result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const |
| { return static_cast<result_type>(_impl._alias_table.size() - 1); } |
| |
| /** |
| * Returns a vector containing the probabilities of each |
| * value of the distribution. For example, given |
| * |
| * @code |
| * discrete_distribution<> dist = { 1, 4, 5 }; |
| * std::vector<double> p = dist.param(); |
| * @endcode |
| * |
| * the vector, p will contain {0.1, 0.4, 0.5}. |
| * |
| * If @c WeightType is integral, then the weights |
| * will be returned unchanged. |
| */ |
| std::vector<WeightType> probabilities() const |
| { |
| std::vector<WeightType> result(_impl._alias_table.size(), static_cast<WeightType>(0)); |
| std::size_t i = 0; |
| for(typename impl_type::alias_table_t::const_iterator |
| iter = _impl._alias_table.begin(), |
| end = _impl._alias_table.end(); |
| iter != end; ++iter, ++i) |
| { |
| WeightType val = iter->first; |
| result[i] += val; |
| result[static_cast<std::size_t>(iter->second)] += _impl.get_weight(i) - val; |
| } |
| impl_type::normalize(result); |
| return(result); |
| } |
| |
| /** Returns the parameters of the distribution. */ |
| param_type param() const |
| { |
| return param_type(*this); |
| } |
| /** Sets the parameters of the distribution. */ |
| void param(const param_type& parm) |
| { |
| init(parm._probabilities.begin(), parm._probabilities.end()); |
| } |
| |
| /** |
| * Effects: Subsequent uses of the distribution do not depend |
| * on values produced by any engine prior to invoking reset. |
| */ |
| void reset() {} |
| |
| /** Writes a distribution to a @c std::ostream. */ |
| BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, discrete_distribution, dd) |
| { |
| os << dd.param(); |
| return os; |
| } |
| |
| /** Reads a distribution from a @c std::istream */ |
| BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, discrete_distribution, dd) |
| { |
| param_type parm; |
| if(is >> parm) { |
| dd.param(parm); |
| } |
| return is; |
| } |
| |
| /** |
| * Returns true if the two distributions will return the |
| * same sequence of values, when passed equal generators. |
| */ |
| BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(discrete_distribution, lhs, rhs) |
| { |
| return lhs._impl == rhs._impl; |
| } |
| /** |
| * Returns true if the two distributions may return different |
| * sequences of values, when passed equal generators. |
| */ |
| BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(discrete_distribution) |
| |
| private: |
| |
| /// @cond show_private |
| |
| template<class Iter> |
| void init(Iter first, Iter last, std::input_iterator_tag) |
| { |
| std::vector<WeightType> temp(first, last); |
| init(temp.begin(), temp.end()); |
| } |
| template<class Iter> |
| void init(Iter first, Iter last, std::forward_iterator_tag) |
| { |
| std::vector<std::pair<WeightType, IntType> > below_average; |
| std::vector<std::pair<WeightType, IntType> > above_average; |
| WeightType weight_average = _impl.init_average(first, last); |
| WeightType normalized_average = _impl.get_weight(0); |
| std::size_t i = 0; |
| for(; first != last; ++first, ++i) { |
| WeightType val = impl_type::normalize(*first, weight_average); |
| std::pair<WeightType, IntType> elem(val, static_cast<IntType>(i)); |
| if(val < normalized_average) { |
| below_average.push_back(elem); |
| } else { |
| above_average.push_back(elem); |
| } |
| } |
| |
| typename impl_type::alias_table_t::iterator |
| b_iter = below_average.begin(), |
| b_end = below_average.end(), |
| a_iter = above_average.begin(), |
| a_end = above_average.end() |
| ; |
| while(b_iter != b_end && a_iter != a_end) { |
| _impl._alias_table[static_cast<std::size_t>(b_iter->second)] = |
| std::make_pair(b_iter->first, a_iter->second); |
| a_iter->first -= (_impl.get_weight(b_iter->second) - b_iter->first); |
| if(a_iter->first < normalized_average) { |
| *b_iter = *a_iter++; |
| } else { |
| ++b_iter; |
| } |
| } |
| for(; b_iter != b_end; ++b_iter) { |
| _impl._alias_table[static_cast<std::size_t>(b_iter->second)].first = |
| _impl.get_weight(b_iter->second); |
| } |
| for(; a_iter != a_end; ++a_iter) { |
| _impl._alias_table[static_cast<std::size_t>(a_iter->second)].first = |
| _impl.get_weight(a_iter->second); |
| } |
| } |
| template<class Iter> |
| void init(Iter first, Iter last) |
| { |
| if(first == last) { |
| _impl.init_empty(); |
| } else { |
| typename std::iterator_traits<Iter>::iterator_category category; |
| init(first, last, category); |
| } |
| } |
| typedef typename detail::select_alias_table< |
| (::boost::is_integral<WeightType>::value) |
| >::template apply<IntType, WeightType>::type impl_type; |
| impl_type _impl; |
| /// @endcond |
| }; |
| |
| } |
| } |
| |
| #include <boost/random/detail/enable_warnings.hpp> |
| |
| #endif |