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