Mercurial > hg > vamp-build-and-test
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 |