discrete_distribution.hpp 15 KB

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