boost/random/discrete_distribution.hpp

   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 85813 2013-09-21 20:17:00Z jewillco $
  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_CXX11_HDR_INITIALIZER_LIST
  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          */
  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>
  68         param_type(Iter first, Iter last) : _probabilities(first, last)
  69         {
  70             normalize();
  71         }
  72 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
  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          */
  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>
  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>
 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          */
 120         std::vector<WeightType> probabilities() const
 121         {
 122             return _probabilities;
 123         }
 124
 125         /** Writes the parameters to a @c std::ostream. */
 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. */
 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. */
 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;
 153         explicit param_type(const discrete_distribution& dist)
 154           : _probabilities(dist.probabilities())
 155         {}
 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      */
 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>
 189     discrete_distribution(Iter first, Iter last)
 190     {
 191         init(first, last);
 192     }
 193 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
 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      */
 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>
 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>
 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      */
 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>
 258     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>
 275     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. */
 295     result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; }
 296     /** Returns the largest value that the distribution can produce. */
 297     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      */
 311     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. */
 330     param_type param() const
 331     {
 332         return param_type(*this);
 333     }
 334     /** Sets the parameters of the distribution. */
 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      */
 344     void reset() {}
 345
 346     /** Writes a distribution to a @c std::ostream. */
 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 */
 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      */
 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>
 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>
 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>
 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