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 |
532 | |