]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | ||
32 | namespace boost { | |
33 | namespace random { | |
34 | ||
35 | /** | |
36 | * The class @c piecewise_constant_distribution models a \random_distribution. | |
37 | */ | |
38 | template<class RealType = double, class WeightType = double> | |
39 | class piecewise_constant_distribution { | |
40 | public: | |
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 | ||
458 | private: | |
459 | discrete_distribution<std::size_t, WeightType> _bins; | |
460 | std::vector<RealType> _intervals; | |
461 | }; | |
462 | ||
463 | } | |
464 | } | |
465 | ||
466 | #endif |