You can not select more than 25 topics Topics must start with a chinese character,a letter or number, can include dashes ('-') and can be up to 35 characters long.

discrete_distribution.h 9.6 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307
  1. // Copyright 2017 The Abseil Authors.
  2. //
  3. // Licensed under the Apache License, Version 2.0 (the "License");
  4. // you may not use this file except in compliance with the License.
  5. // You may obtain a copy of the License at
  6. //
  7. // https://www.apache.org/licenses/LICENSE-2.0
  8. //
  9. // Unless required by applicable law or agreed to in writing, software
  10. // distributed under the License is distributed on an "AS IS" BASIS,
  11. // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  12. // See the License for the specific language governing permissions and
  13. // limitations under the License.
  14. #ifndef ABSL_RANDOM_DISCRETE_DISTRIBUTION_H_
  15. #define ABSL_RANDOM_DISCRETE_DISTRIBUTION_H_
  16. #include <cassert>
  17. #include <cmath>
  18. #include <istream>
  19. #include <limits>
  20. #include <numeric>
  21. #include <type_traits>
  22. #include <utility>
  23. #include <vector>
  24. #include "absl/random/bernoulli_distribution.h"
  25. #include "absl/random/internal/iostream_state_saver.h"
  26. #include "absl/random/uniform_int_distribution.h"
  27. namespace absl
  28. {
  29. ABSL_NAMESPACE_BEGIN
  30. // absl::discrete_distribution
  31. //
  32. // A discrete distribution produces random integers i, where 0 <= i < n
  33. // distributed according to the discrete probability function:
  34. //
  35. // P(i|p0,...,pn−1)=pi
  36. //
  37. // This class is an implementation of discrete_distribution (see
  38. // [rand.dist.samp.discrete]).
  39. //
  40. // The algorithm used is Walker's Aliasing algorithm, described in Knuth, Vol 2.
  41. // absl::discrete_distribution takes O(N) time to precompute the probabilities
  42. // (where N is the number of possible outcomes in the distribution) at
  43. // construction, and then takes O(1) time for each variate generation. Many
  44. // other implementations also take O(N) time to construct an ordered sequence of
  45. // partial sums, plus O(log N) time per variate to binary search.
  46. //
  47. template<typename IntType = int>
  48. class discrete_distribution
  49. {
  50. public:
  51. using result_type = IntType;
  52. class param_type
  53. {
  54. public:
  55. using distribution_type = discrete_distribution;
  56. param_type()
  57. {
  58. init();
  59. }
  60. template<typename InputIterator>
  61. explicit param_type(InputIterator begin, InputIterator end) :
  62. p_(begin, end)
  63. {
  64. init();
  65. }
  66. explicit param_type(std::initializer_list<double> weights) :
  67. p_(weights)
  68. {
  69. init();
  70. }
  71. template<class UnaryOperation>
  72. explicit param_type(size_t nw, double xmin, double xmax, UnaryOperation fw)
  73. {
  74. if (nw > 0)
  75. {
  76. p_.reserve(nw);
  77. double delta = (xmax - xmin) / static_cast<double>(nw);
  78. assert(delta > 0);
  79. double t = delta * 0.5;
  80. for (size_t i = 0; i < nw; ++i)
  81. {
  82. p_.push_back(fw(xmin + i * delta + t));
  83. }
  84. }
  85. init();
  86. }
  87. const std::vector<double>& probabilities() const
  88. {
  89. return p_;
  90. }
  91. size_t n() const
  92. {
  93. return p_.size() - 1;
  94. }
  95. friend bool operator==(const param_type& a, const param_type& b)
  96. {
  97. return a.probabilities() == b.probabilities();
  98. }
  99. friend bool operator!=(const param_type& a, const param_type& b)
  100. {
  101. return !(a == b);
  102. }
  103. private:
  104. friend class discrete_distribution;
  105. void init();
  106. std::vector<double> p_; // normalized probabilities
  107. std::vector<std::pair<double, size_t>> q_; // (acceptance, alternate) pairs
  108. static_assert(std::is_integral<result_type>::value, "Class-template absl::discrete_distribution<> must be "
  109. "parameterized using an integral type.");
  110. };
  111. discrete_distribution() :
  112. param_()
  113. {
  114. }
  115. explicit discrete_distribution(const param_type& p) :
  116. param_(p)
  117. {
  118. }
  119. template<typename InputIterator>
  120. explicit discrete_distribution(InputIterator begin, InputIterator end) :
  121. param_(begin, end)
  122. {
  123. }
  124. explicit discrete_distribution(std::initializer_list<double> weights) :
  125. param_(weights)
  126. {
  127. }
  128. template<class UnaryOperation>
  129. explicit discrete_distribution(size_t nw, double xmin, double xmax, UnaryOperation fw) :
  130. param_(nw, xmin, xmax, std::move(fw))
  131. {
  132. }
  133. void reset()
  134. {
  135. }
  136. // generating functions
  137. template<typename URBG>
  138. result_type operator()(URBG& g)
  139. { // NOLINT(runtime/references)
  140. return (*this)(g, param_);
  141. }
  142. template<typename URBG>
  143. result_type operator()(URBG& g, // NOLINT(runtime/references)
  144. const param_type& p);
  145. const param_type& param() const
  146. {
  147. return param_;
  148. }
  149. void param(const param_type& p)
  150. {
  151. param_ = p;
  152. }
  153. result_type(min)() const
  154. {
  155. return 0;
  156. }
  157. result_type(max)() const
  158. {
  159. return static_cast<result_type>(param_.n());
  160. } // inclusive
  161. // NOTE [rand.dist.sample.discrete] returns a std::vector<double> not a
  162. // const std::vector<double>&.
  163. const std::vector<double>& probabilities() const
  164. {
  165. return param_.probabilities();
  166. }
  167. friend bool operator==(const discrete_distribution& a, const discrete_distribution& b)
  168. {
  169. return a.param_ == b.param_;
  170. }
  171. friend bool operator!=(const discrete_distribution& a, const discrete_distribution& b)
  172. {
  173. return a.param_ != b.param_;
  174. }
  175. private:
  176. param_type param_;
  177. };
  178. // --------------------------------------------------------------------------
  179. // Implementation details only below
  180. // --------------------------------------------------------------------------
  181. namespace random_internal
  182. {
  183. // Using the vector `*probabilities`, whose values are the weights or
  184. // probabilities of an element being selected, constructs the proportional
  185. // probabilities used by the discrete distribution. `*probabilities` will be
  186. // scaled, if necessary, so that its entries sum to a value sufficiently close
  187. // to 1.0.
  188. std::vector<std::pair<double, size_t>> InitDiscreteDistribution(
  189. std::vector<double>* probabilities
  190. );
  191. } // namespace random_internal
  192. template<typename IntType>
  193. void discrete_distribution<IntType>::param_type::init()
  194. {
  195. if (p_.empty())
  196. {
  197. p_.push_back(1.0);
  198. q_.emplace_back(1.0, 0);
  199. }
  200. else
  201. {
  202. assert(n() <= (std::numeric_limits<IntType>::max)());
  203. q_ = random_internal::InitDiscreteDistribution(&p_);
  204. }
  205. }
  206. template<typename IntType>
  207. template<typename URBG>
  208. typename discrete_distribution<IntType>::result_type
  209. discrete_distribution<IntType>::operator()(
  210. URBG& g, // NOLINT(runtime/references)
  211. const param_type& p
  212. )
  213. {
  214. const auto idx = absl::uniform_int_distribution<result_type>(0, p.n())(g);
  215. const auto& q = p.q_[idx];
  216. const bool selected = absl::bernoulli_distribution(q.first)(g);
  217. return selected ? idx : static_cast<result_type>(q.second);
  218. }
  219. template<typename CharT, typename Traits, typename IntType>
  220. std::basic_ostream<CharT, Traits>& operator<<(
  221. std::basic_ostream<CharT, Traits>& os, // NOLINT(runtime/references)
  222. const discrete_distribution<IntType>& x
  223. )
  224. {
  225. auto saver = random_internal::make_ostream_state_saver(os);
  226. const auto& probabilities = x.param().probabilities();
  227. os << probabilities.size();
  228. os.precision(random_internal::stream_precision_helper<double>::kPrecision);
  229. for (const auto& p : probabilities)
  230. {
  231. os << os.fill() << p;
  232. }
  233. return os;
  234. }
  235. template<typename CharT, typename Traits, typename IntType>
  236. std::basic_istream<CharT, Traits>& operator>>(
  237. std::basic_istream<CharT, Traits>& is, // NOLINT(runtime/references)
  238. discrete_distribution<IntType>& x
  239. )
  240. { // NOLINT(runtime/references)
  241. using param_type = typename discrete_distribution<IntType>::param_type;
  242. auto saver = random_internal::make_istream_state_saver(is);
  243. size_t n;
  244. std::vector<double> p;
  245. is >> n;
  246. if (is.fail())
  247. return is;
  248. if (n > 0)
  249. {
  250. p.reserve(n);
  251. for (IntType i = 0; i < n && !is.fail(); ++i)
  252. {
  253. auto tmp = random_internal::read_floating_point<double>(is);
  254. if (is.fail())
  255. return is;
  256. p.push_back(tmp);
  257. }
  258. }
  259. x.param(param_type(p.begin(), p.end()));
  260. return is;
  261. }
  262. ABSL_NAMESPACE_END
  263. } // namespace absl
  264. #endif // ABSL_RANDOM_DISCRETE_DISTRIBUTION_H_