]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | /* boost random/piecewise_linear_distribution.hpp header file |
2 | * | |
3 | * Copyright Steven Watanabe 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$ | |
11 | */ | |
12 | ||
13 | #ifndef BOOST_RANDOM_PIECEWISE_LINEAR_DISTRIBUTION_HPP_INCLUDED | |
14 | #define BOOST_RANDOM_PIECEWISE_LINEAR_DISTRIBUTION_HPP_INCLUDED | |
15 | ||
16 | #include <vector> | |
17 | #include <algorithm> | |
18 | #include <cmath> | |
19 | #include <cstdlib> | |
20 | #include <boost/assert.hpp> | |
21 | #include <boost/random/uniform_real.hpp> | |
22 | #include <boost/random/discrete_distribution.hpp> | |
23 | #include <boost/random/detail/config.hpp> | |
24 | #include <boost/random/detail/operators.hpp> | |
25 | #include <boost/random/detail/vector_io.hpp> | |
26 | ||
27 | #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST | |
28 | #include <initializer_list> | |
29 | #endif | |
30 | ||
31 | #include <boost/range/begin.hpp> | |
32 | #include <boost/range/end.hpp> | |
33 | ||
34 | namespace boost { | |
35 | namespace random { | |
36 | ||
37 | /** | |
38 | * The class @c piecewise_linear_distribution models a \random_distribution. | |
39 | */ | |
40 | template<class RealType = double> | |
41 | class piecewise_linear_distribution { | |
42 | public: | |
43 | typedef std::size_t input_type; | |
44 | typedef RealType result_type; | |
45 | ||
46 | class param_type { | |
47 | public: | |
48 | ||
49 | typedef piecewise_linear_distribution distribution_type; | |
50 | ||
51 | /** | |
52 | * Constructs a @c param_type object, representing a distribution | |
53 | * that produces values uniformly distributed in the range [0, 1). | |
54 | */ | |
55 | param_type() | |
56 | { | |
57 | _weights.push_back(RealType(1)); | |
58 | _weights.push_back(RealType(1)); | |
59 | _intervals.push_back(RealType(0)); | |
60 | _intervals.push_back(RealType(1)); | |
61 | } | |
62 | /** | |
63 | * Constructs a @c param_type object from two iterator ranges | |
64 | * containing the interval boundaries and weights at the boundaries. | |
65 | * If there are fewer than two boundaries, then this is equivalent to | |
66 | * the default constructor and the distribution will produce values | |
67 | * uniformly distributed in the range [0, 1). | |
68 | * | |
69 | * The values of the interval boundaries must be strictly | |
70 | * increasing, and the number of weights must be the same as | |
71 | * the number of interval boundaries. If there are extra | |
72 | * weights, they are ignored. | |
73 | */ | |
74 | template<class IntervalIter, class WeightIter> | |
75 | param_type(IntervalIter intervals_first, IntervalIter intervals_last, | |
76 | WeightIter weight_first) | |
77 | : _intervals(intervals_first, intervals_last) | |
78 | { | |
79 | if(_intervals.size() < 2) { | |
80 | _intervals.clear(); | |
81 | _weights.push_back(RealType(1)); | |
82 | _weights.push_back(RealType(1)); | |
83 | _intervals.push_back(RealType(0)); | |
84 | _intervals.push_back(RealType(1)); | |
85 | } else { | |
86 | _weights.reserve(_intervals.size()); | |
87 | for(std::size_t i = 0; i < _intervals.size(); ++i) { | |
88 | _weights.push_back(*weight_first++); | |
89 | } | |
90 | } | |
91 | } | |
92 | #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST | |
93 | /** | |
94 | * Constructs a @c param_type object from an initializer_list | |
95 | * containing the interval boundaries and a unary function | |
96 | * specifying the weights at the boundaries. Each weight is | |
97 | * determined by calling the function at the corresponding point. | |
98 | * | |
99 | * If the initializer_list contains fewer than two elements, | |
100 | * this is equivalent to the default constructor and the | |
101 | * distribution will produce values uniformly distributed | |
102 | * in the range [0, 1). | |
103 | */ | |
104 | template<class T, class F> | |
105 | param_type(const std::initializer_list<T>& il, F f) | |
106 | : _intervals(il.begin(), il.end()) | |
107 | { | |
108 | if(_intervals.size() < 2) { | |
109 | _intervals.clear(); | |
110 | _weights.push_back(RealType(1)); | |
111 | _weights.push_back(RealType(1)); | |
112 | _intervals.push_back(RealType(0)); | |
113 | _intervals.push_back(RealType(1)); | |
114 | } else { | |
115 | _weights.reserve(_intervals.size()); | |
116 | for(typename std::vector<RealType>::const_iterator | |
117 | iter = _intervals.begin(), end = _intervals.end(); | |
118 | iter != end; ++iter) | |
119 | { | |
120 | _weights.push_back(f(*iter)); | |
121 | } | |
122 | } | |
123 | } | |
124 | #endif | |
125 | /** | |
126 | * Constructs a @c param_type object from Boost.Range ranges holding | |
127 | * the interval boundaries and the weights at the boundaries. If | |
128 | * there are fewer than two interval boundaries, this is equivalent | |
129 | * to the default constructor and the distribution will produce | |
130 | * values uniformly distributed in the range [0, 1). The | |
131 | * number of weights must be equal to the number of | |
132 | * interval boundaries. | |
133 | */ | |
134 | template<class IntervalRange, class WeightRange> | |
135 | param_type(const IntervalRange& intervals_arg, | |
136 | const WeightRange& weights_arg) | |
137 | : _intervals(boost::begin(intervals_arg), boost::end(intervals_arg)), | |
138 | _weights(boost::begin(weights_arg), boost::end(weights_arg)) | |
139 | { | |
140 | if(_intervals.size() < 2) { | |
141 | _weights.clear(); | |
142 | _weights.push_back(RealType(1)); | |
143 | _weights.push_back(RealType(1)); | |
144 | _intervals.clear(); | |
145 | _intervals.push_back(RealType(0)); | |
146 | _intervals.push_back(RealType(1)); | |
147 | } | |
148 | } | |
149 | ||
150 | /** | |
151 | * Constructs the parameters for a distribution that approximates a | |
152 | * function. The range of the distribution is [xmin, xmax). This | |
153 | * range is divided into nw equally sized intervals and the weights | |
154 | * are found by calling the unary function f on the boundaries of the | |
155 | * intervals. | |
156 | */ | |
157 | template<class F> | |
158 | param_type(std::size_t nw, RealType xmin, RealType xmax, F f) | |
159 | { | |
160 | std::size_t n = (nw == 0) ? 1 : nw; | |
161 | double delta = (xmax - xmin) / n; | |
162 | BOOST_ASSERT(delta > 0); | |
163 | for(std::size_t k = 0; k < n; ++k) { | |
164 | _weights.push_back(f(xmin + k*delta)); | |
165 | _intervals.push_back(xmin + k*delta); | |
166 | } | |
167 | _weights.push_back(f(xmax)); | |
168 | _intervals.push_back(xmax); | |
169 | } | |
170 | ||
171 | /** Returns a vector containing the interval boundaries. */ | |
172 | std::vector<RealType> intervals() const { return _intervals; } | |
173 | ||
174 | /** | |
175 | * Returns a vector containing the probability densities | |
176 | * at all the interval boundaries. | |
177 | */ | |
178 | std::vector<RealType> densities() const | |
179 | { | |
180 | RealType sum = static_cast<RealType>(0); | |
181 | for(std::size_t i = 0; i < _intervals.size() - 1; ++i) { | |
182 | RealType width = _intervals[i + 1] - _intervals[i]; | |
183 | sum += (_weights[i] + _weights[i + 1]) * width / 2; | |
184 | } | |
185 | std::vector<RealType> result; | |
186 | result.reserve(_weights.size()); | |
187 | for(typename std::vector<RealType>::const_iterator | |
188 | iter = _weights.begin(), end = _weights.end(); | |
189 | iter != end; ++iter) | |
190 | { | |
191 | result.push_back(*iter / sum); | |
192 | } | |
193 | return result; | |
194 | } | |
195 | ||
196 | /** Writes the parameters to a @c std::ostream. */ | |
197 | BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm) | |
198 | { | |
199 | detail::print_vector(os, parm._intervals); | |
200 | detail::print_vector(os, parm._weights); | |
201 | return os; | |
202 | } | |
203 | ||
204 | /** Reads the parameters from a @c std::istream. */ | |
205 | BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm) | |
206 | { | |
207 | std::vector<RealType> new_intervals; | |
208 | std::vector<RealType> new_weights; | |
209 | detail::read_vector(is, new_intervals); | |
210 | detail::read_vector(is, new_weights); | |
211 | if(is) { | |
212 | parm._intervals.swap(new_intervals); | |
213 | parm._weights.swap(new_weights); | |
214 | } | |
215 | return is; | |
216 | } | |
217 | ||
218 | /** Returns true if the two sets of parameters are the same. */ | |
219 | BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs) | |
220 | { | |
221 | return lhs._intervals == rhs._intervals | |
222 | && lhs._weights == rhs._weights; | |
223 | } | |
224 | /** Returns true if the two sets of parameters are different. */ | |
225 | BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type) | |
226 | ||
227 | private: | |
228 | friend class piecewise_linear_distribution; | |
229 | ||
230 | std::vector<RealType> _intervals; | |
231 | std::vector<RealType> _weights; | |
232 | }; | |
233 | ||
234 | /** | |
235 | * Creates a new @c piecewise_linear_distribution that | |
236 | * produces values uniformly distributed in the range [0, 1). | |
237 | */ | |
238 | piecewise_linear_distribution() | |
239 | { | |
240 | default_init(); | |
241 | } | |
242 | /** | |
243 | * Constructs a piecewise_linear_distribution from two iterator ranges | |
244 | * containing the interval boundaries and the weights at the boundaries. | |
245 | * If there are fewer than two boundaries, then this is equivalent to | |
246 | * the default constructor and creates a distribution that | |
247 | * produces values uniformly distributed in the range [0, 1). | |
248 | * | |
249 | * The values of the interval boundaries must be strictly | |
250 | * increasing, and the number of weights must be equal to | |
251 | * the number of interval boundaries. If there are extra | |
252 | * weights, they are ignored. | |
253 | * | |
254 | * For example, | |
255 | * | |
256 | * @code | |
257 | * double intervals[] = { 0.0, 1.0, 2.0 }; | |
258 | * double weights[] = { 0.0, 1.0, 0.0 }; | |
259 | * piecewise_constant_distribution<> dist( | |
260 | * &intervals[0], &intervals[0] + 3, &weights[0]); | |
261 | * @endcode | |
262 | * | |
263 | * produces a triangle distribution. | |
264 | */ | |
265 | template<class IntervalIter, class WeightIter> | |
266 | piecewise_linear_distribution(IntervalIter first_interval, | |
267 | IntervalIter last_interval, | |
268 | WeightIter first_weight) | |
269 | : _intervals(first_interval, last_interval) | |
270 | { | |
271 | if(_intervals.size() < 2) { | |
272 | default_init(); | |
273 | } else { | |
274 | _weights.reserve(_intervals.size()); | |
275 | for(std::size_t i = 0; i < _intervals.size(); ++i) { | |
276 | _weights.push_back(*first_weight++); | |
277 | } | |
278 | init(); | |
279 | } | |
280 | } | |
281 | #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST | |
282 | /** | |
283 | * Constructs a piecewise_linear_distribution from an | |
284 | * initializer_list containing the interval boundaries | |
285 | * and a unary function specifying the weights. Each | |
286 | * weight is determined by calling the function at the | |
287 | * corresponding interval boundary. | |
288 | * | |
289 | * If the initializer_list contains fewer than two elements, | |
290 | * this is equivalent to the default constructor and the | |
291 | * distribution will produce values uniformly distributed | |
292 | * in the range [0, 1). | |
293 | */ | |
294 | template<class T, class F> | |
295 | piecewise_linear_distribution(std::initializer_list<T> il, F f) | |
296 | : _intervals(il.begin(), il.end()) | |
297 | { | |
298 | if(_intervals.size() < 2) { | |
299 | default_init(); | |
300 | } else { | |
301 | _weights.reserve(_intervals.size()); | |
302 | for(typename std::vector<RealType>::const_iterator | |
303 | iter = _intervals.begin(), end = _intervals.end(); | |
304 | iter != end; ++iter) | |
305 | { | |
306 | _weights.push_back(f(*iter)); | |
307 | } | |
308 | init(); | |
309 | } | |
310 | } | |
311 | #endif | |
312 | /** | |
313 | * Constructs a piecewise_linear_distribution from Boost.Range | |
314 | * ranges holding the interval boundaries and the weights. If | |
315 | * there are fewer than two interval boundaries, this is equivalent | |
316 | * to the default constructor and the distribution will produce | |
317 | * values uniformly distributed in the range [0, 1). The | |
318 | * number of weights must be equal to the number of | |
319 | * interval boundaries. | |
320 | */ | |
321 | template<class IntervalsRange, class WeightsRange> | |
322 | piecewise_linear_distribution(const IntervalsRange& intervals_arg, | |
323 | const WeightsRange& weights_arg) | |
324 | : _intervals(boost::begin(intervals_arg), boost::end(intervals_arg)), | |
325 | _weights(boost::begin(weights_arg), boost::end(weights_arg)) | |
326 | { | |
327 | if(_intervals.size() < 2) { | |
328 | default_init(); | |
329 | } else { | |
330 | init(); | |
331 | } | |
332 | } | |
333 | /** | |
334 | * Constructs a piecewise_linear_distribution that approximates a | |
335 | * function. The range of the distribution is [xmin, xmax). This | |
336 | * range is divided into nw equally sized intervals and the weights | |
337 | * are found by calling the unary function f on the interval boundaries. | |
338 | */ | |
339 | template<class F> | |
340 | piecewise_linear_distribution(std::size_t nw, | |
341 | RealType xmin, | |
342 | RealType xmax, | |
343 | F f) | |
344 | { | |
345 | if(nw == 0) { nw = 1; } | |
346 | RealType delta = (xmax - xmin) / nw; | |
347 | _intervals.reserve(nw + 1); | |
348 | for(std::size_t i = 0; i < nw; ++i) { | |
349 | RealType x = xmin + i * delta; | |
350 | _intervals.push_back(x); | |
351 | _weights.push_back(f(x)); | |
352 | } | |
353 | _intervals.push_back(xmax); | |
354 | _weights.push_back(f(xmax)); | |
355 | init(); | |
356 | } | |
357 | /** | |
358 | * Constructs a piecewise_linear_distribution from its parameters. | |
359 | */ | |
360 | explicit piecewise_linear_distribution(const param_type& parm) | |
361 | : _intervals(parm._intervals), | |
362 | _weights(parm._weights) | |
363 | { | |
364 | init(); | |
365 | } | |
366 | ||
367 | /** | |
368 | * Returns a value distributed according to the parameters of the | |
369 | * piecewise_linear_distribution. | |
370 | */ | |
371 | template<class URNG> | |
372 | RealType operator()(URNG& urng) const | |
373 | { | |
374 | std::size_t i = _bins(urng); | |
375 | bool is_in_rectangle = (i % 2 == 0); | |
376 | i = i / 2; | |
377 | uniform_real<RealType> dist(_intervals[i], _intervals[i+1]); | |
378 | if(is_in_rectangle) { | |
379 | return dist(urng); | |
380 | } else if(_weights[i] < _weights[i+1]) { | |
381 | return (std::max)(dist(urng), dist(urng)); | |
382 | } else { | |
383 | return (std::min)(dist(urng), dist(urng)); | |
384 | } | |
385 | } | |
386 | ||
387 | /** | |
388 | * Returns a value distributed according to the parameters | |
389 | * specified by param. | |
390 | */ | |
391 | template<class URNG> | |
392 | RealType operator()(URNG& urng, const param_type& parm) const | |
393 | { | |
394 | return piecewise_linear_distribution(parm)(urng); | |
395 | } | |
396 | ||
397 | /** Returns the smallest value that the distribution can produce. */ | |
398 | result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const | |
399 | { return _intervals.front(); } | |
400 | /** Returns the largest value that the distribution can produce. */ | |
401 | result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const | |
402 | { return _intervals.back(); } | |
403 | ||
404 | /** | |
405 | * Returns a vector containing the probability densities | |
406 | * at the interval boundaries. | |
407 | */ | |
408 | std::vector<RealType> densities() const | |
409 | { | |
410 | RealType sum = static_cast<RealType>(0); | |
411 | for(std::size_t i = 0; i < _intervals.size() - 1; ++i) { | |
412 | RealType width = _intervals[i + 1] - _intervals[i]; | |
413 | sum += (_weights[i] + _weights[i + 1]) * width / 2; | |
414 | } | |
415 | std::vector<RealType> result; | |
416 | result.reserve(_weights.size()); | |
417 | for(typename std::vector<RealType>::const_iterator | |
418 | iter = _weights.begin(), end = _weights.end(); | |
419 | iter != end; ++iter) | |
420 | { | |
421 | result.push_back(*iter / sum); | |
422 | } | |
423 | return result; | |
424 | } | |
425 | /** Returns a vector containing the interval boundaries. */ | |
426 | std::vector<RealType> intervals() const { return _intervals; } | |
427 | ||
428 | /** Returns the parameters of the distribution. */ | |
429 | param_type param() const | |
430 | { | |
431 | return param_type(_intervals, _weights); | |
432 | } | |
433 | /** Sets the parameters of the distribution. */ | |
434 | void param(const param_type& parm) | |
435 | { | |
436 | std::vector<RealType> new_intervals(parm._intervals); | |
437 | std::vector<RealType> new_weights(parm._weights); | |
438 | init(new_intervals, new_weights); | |
439 | _intervals.swap(new_intervals); | |
440 | _weights.swap(new_weights); | |
441 | } | |
442 | ||
443 | /** | |
444 | * Effects: Subsequent uses of the distribution do not depend | |
445 | * on values produced by any engine prior to invoking reset. | |
446 | */ | |
447 | void reset() { _bins.reset(); } | |
448 | ||
449 | /** Writes a distribution to a @c std::ostream. */ | |
450 | BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR( | |
451 | os, piecewise_linear_distribution, pld) | |
452 | { | |
453 | os << pld.param(); | |
454 | return os; | |
455 | } | |
456 | ||
457 | /** Reads a distribution from a @c std::istream */ | |
458 | BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR( | |
459 | is, piecewise_linear_distribution, pld) | |
460 | { | |
461 | param_type parm; | |
462 | if(is >> parm) { | |
463 | pld.param(parm); | |
464 | } | |
465 | return is; | |
466 | } | |
467 | ||
468 | /** | |
469 | * Returns true if the two distributions will return the | |
470 | * same sequence of values, when passed equal generators. | |
471 | */ | |
472 | BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR( | |
473 | piecewise_linear_distribution, lhs, rhs) | |
474 | { | |
475 | return lhs._intervals == rhs._intervals && lhs._weights == rhs._weights; | |
476 | } | |
477 | /** | |
478 | * Returns true if the two distributions may return different | |
479 | * sequences of values, when passed equal generators. | |
480 | */ | |
481 | BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(piecewise_linear_distribution) | |
482 | ||
483 | private: | |
484 | ||
485 | /// @cond \show_private | |
486 | ||
487 | void init(const std::vector<RealType>& intervals_arg, | |
488 | const std::vector<RealType>& weights_arg) | |
489 | { | |
490 | using std::abs; | |
491 | std::vector<RealType> bin_weights; | |
492 | bin_weights.reserve((intervals_arg.size() - 1) * 2); | |
493 | for(std::size_t i = 0; i < intervals_arg.size() - 1; ++i) { | |
494 | RealType width = intervals_arg[i + 1] - intervals_arg[i]; | |
495 | RealType w1 = weights_arg[i]; | |
496 | RealType w2 = weights_arg[i + 1]; | |
497 | bin_weights.push_back((std::min)(w1, w2) * width); | |
498 | bin_weights.push_back(abs(w1 - w2) * width / 2); | |
499 | } | |
500 | typedef discrete_distribution<std::size_t, RealType> bins_type; | |
501 | typename bins_type::param_type bins_param(bin_weights); | |
502 | _bins.param(bins_param); | |
503 | } | |
504 | ||
505 | void init() | |
506 | { | |
507 | init(_intervals, _weights); | |
508 | } | |
509 | ||
510 | void default_init() | |
511 | { | |
512 | _intervals.clear(); | |
513 | _intervals.push_back(RealType(0)); | |
514 | _intervals.push_back(RealType(1)); | |
515 | _weights.clear(); | |
516 | _weights.push_back(RealType(1)); | |
517 | _weights.push_back(RealType(1)); | |
518 | init(); | |
519 | } | |
520 | ||
521 | discrete_distribution<std::size_t, RealType> _bins; | |
522 | std::vector<RealType> _intervals; | |
523 | std::vector<RealType> _weights; | |
524 | ||
525 | /// @endcond | |
526 | }; | |
527 | ||
528 | } | |
529 | } | |
530 | ||
531 | #endif |