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
34namespace boost {
35namespace random {
36
37/**
38 * The class @c piecewise_linear_distribution models a \random_distribution.
39 */
40template<class RealType = double>
41class piecewise_linear_distribution {
42public:
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
483private:
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