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          */
param_type()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>
param_type(IntervalIter intervals_first,IntervalIter intervals_last,WeightIter weight_first)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>
param_type(const std::initializer_list<T> & il,F 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>
param_type(const IntervalRange & intervals_arg,const WeightRange & weights_arg)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>
param_type(std::size_t nw,RealType xmin,RealType xmax,F 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. */
intervals() const172         std::vector<RealType> intervals() const { return _intervals; }
173 
174         /**
175          * Returns a vector containing the probability densities
176          * at all the interval boundaries.
177          */
densities() const178         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. */
BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os,param_type,parm)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. */
BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is,param_type,parm)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. */
BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type,lhs,rhs)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      */
piecewise_linear_distribution()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>
piecewise_linear_distribution(IntervalIter first_interval,IntervalIter last_interval,WeightIter first_weight)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>
piecewise_linear_distribution(std::initializer_list<T> il,F 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>
piecewise_linear_distribution(const IntervalsRange & intervals_arg,const WeightsRange & weights_arg)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>
piecewise_linear_distribution(std::size_t nw,RealType xmin,RealType xmax,F 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      */
piecewise_linear_distribution(const param_type & parm)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>
operator ()(URNG & urng) const372     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>
operator ()(URNG & urng,const param_type & parm) const392     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. */
BOOST_PREVENT_MACRO_SUBSTITUTION() const398     result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const
399     { return _intervals.front(); }
400     /** Returns the largest value that the distribution can produce. */
BOOST_PREVENT_MACRO_SUBSTITUTION() const401     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      */
densities() const408     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. */
intervals() const426     std::vector<RealType> intervals() const { return _intervals; }
427 
428     /** Returns the parameters of the distribution. */
param() const429     param_type param() const
430     {
431         return param_type(_intervals, _weights);
432     }
433     /** Sets the parameters of the distribution. */
param(const param_type & parm)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      */
reset()447     void reset() { _bins.reset(); }
448 
449     /** Writes a distribution to a @c std::ostream. */
BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os,piecewise_linear_distribution,pld)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 */
BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is,piecewise_linear_distribution,pld)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      */
BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(piecewise_linear_distribution,lhs,rhs)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 
init(const std::vector<RealType> & intervals_arg,const std::vector<RealType> & weights_arg)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 
init()505     void init()
506     {
507         init(_intervals, _weights);
508     }
509 
default_init()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