comparison DEPENDENCIES/generic/include/boost/random/discrete_distribution.hpp @ 16:2665513ce2d3

Add boost headers
author Chris Cannam
date Tue, 05 Aug 2014 11:11:38 +0100
parents
children c530137014c0
comparison
equal deleted inserted replaced
15:663ca0da4350 16:2665513ce2d3
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