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