]>
Commit | Line | Data |
---|---|---|
1 | /* boost random/hyperexponential_distribution.hpp header file | |
2 | * | |
3 | * Copyright Marco Guazzone 2014 | |
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 | * Much of the code here taken by boost::math::hyperexponential_distribution. | |
11 | * To this end, we would like to thank Paul Bristow and John Maddock for their | |
12 | * valuable feedback. | |
13 | * | |
14 | * \author Marco Guazzone (marco.guazzone@gmail.com) | |
15 | */ | |
16 | ||
17 | #ifndef BOOST_RANDOM_HYPEREXPONENTIAL_DISTRIBUTION_HPP | |
18 | #define BOOST_RANDOM_HYPEREXPONENTIAL_DISTRIBUTION_HPP | |
19 | ||
20 | ||
21 | #include <boost/config.hpp> | |
22 | #include <boost/math/special_functions/fpclassify.hpp> | |
23 | #include <boost/random/detail/operators.hpp> | |
24 | #include <boost/random/detail/vector_io.hpp> | |
25 | #include <boost/random/discrete_distribution.hpp> | |
26 | #include <boost/random/exponential_distribution.hpp> | |
27 | #include <boost/range/begin.hpp> | |
28 | #include <boost/range/end.hpp> | |
29 | #include <boost/range/size.hpp> | |
30 | #include <boost/type_traits/has_pre_increment.hpp> | |
31 | #include <cassert> | |
32 | #include <cmath> | |
33 | #include <cstddef> | |
34 | #include <iterator> | |
35 | #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST | |
36 | # include <initializer_list> | |
37 | #endif // BOOST_NO_CXX11_HDR_INITIALIZER_LIST | |
38 | #include <iostream> | |
39 | #include <limits> | |
40 | #include <numeric> | |
41 | #include <vector> | |
42 | ||
43 | ||
44 | namespace boost { namespace random { | |
45 | ||
46 | namespace hyperexp_detail { | |
47 | ||
48 | template <typename T> | |
49 | std::vector<T>& normalize(std::vector<T>& v) | |
50 | { | |
51 | if (v.size() == 0) | |
52 | { | |
53 | return v; | |
54 | } | |
55 | ||
56 | const T sum = std::accumulate(v.begin(), v.end(), static_cast<T>(0)); | |
57 | T final_sum = 0; | |
58 | ||
59 | const typename std::vector<T>::iterator end = --v.end(); | |
60 | for (typename std::vector<T>::iterator it = v.begin(); | |
61 | it != end; | |
62 | ++it) | |
63 | { | |
64 | *it /= sum; | |
65 | final_sum += *it; | |
66 | } | |
67 | *end = 1-final_sum; // avoids round off errors thus ensuring the probabilities really sum to 1 | |
68 | ||
69 | return v; | |
70 | } | |
71 | ||
72 | template <typename RealT> | |
73 | bool check_probabilities(std::vector<RealT> const& probabilities) | |
74 | { | |
75 | const std::size_t n = probabilities.size(); | |
76 | RealT sum = 0; | |
77 | for (std::size_t i = 0; i < n; ++i) | |
78 | { | |
79 | if (probabilities[i] < 0 | |
80 | || probabilities[i] > 1 | |
81 | || !(boost::math::isfinite)(probabilities[i])) | |
82 | { | |
83 | return false; | |
84 | } | |
85 | sum += probabilities[i]; | |
86 | } | |
87 | ||
88 | //NOTE: the check below seems to fail on some architectures. | |
89 | // So we commented it. | |
90 | //// - We try to keep phase probabilities correctly normalized in the distribution constructors | |
91 | //// - However in practice we have to allow for a very slight divergence from a sum of exactly 1: | |
92 | ////if (std::abs(sum-1) > (std::numeric_limits<RealT>::epsilon()*2)) | |
93 | //// This is from Knuth "The Art of Computer Programming: Vol.2, 3rd Ed", and can be used to | |
94 | //// check is two numbers are approximately equal | |
95 | //const RealT one = 1; | |
96 | //const RealT tol = std::numeric_limits<RealT>::epsilon()*2.0; | |
97 | //if (std::abs(sum-one) > (std::max(std::abs(sum), std::abs(one))*tol)) | |
98 | //{ | |
99 | // return false; | |
100 | //} | |
101 | ||
102 | return true; | |
103 | } | |
104 | ||
105 | template <typename RealT> | |
106 | bool check_rates(std::vector<RealT> const& rates) | |
107 | { | |
108 | const std::size_t n = rates.size(); | |
109 | for (std::size_t i = 0; i < n; ++i) | |
110 | { | |
111 | if (rates[i] <= 0 | |
112 | || !(boost::math::isfinite)(rates[i])) | |
113 | { | |
114 | return false; | |
115 | } | |
116 | } | |
117 | return true; | |
118 | } | |
119 | ||
120 | template <typename RealT> | |
121 | bool check_params(std::vector<RealT> const& probabilities, std::vector<RealT> const& rates) | |
122 | { | |
123 | if (probabilities.size() != rates.size()) | |
124 | { | |
125 | return false; | |
126 | } | |
127 | ||
128 | return check_probabilities(probabilities) | |
129 | && check_rates(rates); | |
130 | } | |
131 | ||
132 | } // Namespace hyperexp_detail | |
133 | ||
134 | ||
135 | /** | |
136 | * The hyperexponential distribution is a real-valued continuous distribution | |
137 | * with two parameters, the <em>phase probability vector</em> \c probs and the | |
138 | * <em>rate vector</em> \c rates. | |
139 | * | |
140 | * A \f$k\f$-phase hyperexponential distribution is a mixture of \f$k\f$ | |
141 | * exponential distributions. | |
142 | * For this reason, it is also referred to as <em>mixed exponential | |
143 | * distribution</em> or <em>parallel \f$k\f$-phase exponential | |
144 | * distribution</em>. | |
145 | * | |
146 | * A \f$k\f$-phase hyperexponential distribution is characterized by two | |
147 | * parameters, namely a <em>phase probability vector</em> \f$\mathbf{\alpha}=(\alpha_1,\ldots,\alpha_k)\f$ and a <em>rate vector</em> \f$\mathbf{\lambda}=(\lambda_1,\ldots,\lambda_k)\f$. | |
148 | * | |
149 | * A \f$k\f$-phase hyperexponential distribution is frequently used in | |
150 | * <em>queueing theory</em> to model the distribution of the superposition of | |
151 | * \f$k\f$ independent events, like, for instance, the service time distribution | |
152 | * of a queueing station with \f$k\f$ servers in parallel where the \f$i\f$-th | |
153 | * server is chosen with probability \f$\alpha_i\f$ and its service time | |
154 | * distribution is an exponential distribution with rate \f$\lambda_i\f$ | |
155 | * (Allen,1990; Papadopolous et al.,1993; Trivedi,2002). | |
156 | * | |
157 | * For instance, CPUs service-time distribution in a computing system has often | |
158 | * been observed to possess such a distribution (Rosin,1965). | |
159 | * Also, the arrival of different types of customer to a single queueing station | |
160 | * is often modeled as a hyperexponential distribution (Papadopolous et al.,1993). | |
161 | * Similarly, if a product manufactured in several parallel assemply lines and | |
162 | * the outputs are merged, the failure density of the overall product is likely | |
163 | * to be hyperexponential (Trivedi,2002). | |
164 | * | |
165 | * Finally, since the hyperexponential distribution exhibits a high Coefficient | |
166 | * of Variation (CoV), that is a CoV > 1, it is especially suited to fit | |
167 | * empirical data with large CoV (Feitelson,2014; Wolski et al.,2013) and to | |
168 | * approximate <em>long-tail probability distributions</em> (Feldmann et al.,1998). | |
169 | * | |
170 | * See (Boost,2014) for more information and examples. | |
171 | * | |
172 | * A \f$k\f$-phase hyperexponential distribution has a probability density | |
173 | * function | |
174 | * \f[ | |
175 | * f(x) = \sum_{i=1}^k \alpha_i \lambda_i e^{-x\lambda_i} | |
176 | * \f] | |
177 | * where: | |
178 | * - \f$k\f$ is the <em>number of phases</em> and also the size of the input | |
179 | * vector parameters, | |
180 | * - \f$\mathbf{\alpha}=(\alpha_1,\ldots,\alpha_k)\f$ is the <em>phase probability | |
181 | * vector</em> parameter, and | |
182 | * - \f$\mathbf{\lambda}=(\lambda_1,\ldots,\lambda_k)\f$ is the <em>rate vector</em> | |
183 | * parameter. | |
184 | * . | |
185 | * | |
186 | * Given a \f$k\f$-phase hyperexponential distribution with phase probability | |
187 | * vector \f$\mathbf{\alpha}\f$ and rate vector \f$\mathbf{\lambda}\f$, the | |
188 | * random variate generation algorithm consists of the following steps (Tyszer,1999): | |
189 | * -# Generate a random variable \f$U\f$ uniformly distribution on the interval \f$(0,1)\f$. | |
190 | * -# Use \f$U\f$ to select the appropriate \f$\lambda_i\f$ (e.g., the | |
191 | * <em>alias method</em> can possibly be used for this step). | |
192 | * -# Generate an exponentially distributed random variable \f$X\f$ with rate parameter \f$\lambda_i\f$. | |
193 | * -# Return \f$X\f$. | |
194 | * . | |
195 | * | |
196 | * References: | |
197 | * -# A.O. Allen, <em>Probability, Statistics, and Queuing Theory with Computer Science Applications, Second Edition</em>, Academic Press, 1990. | |
198 | * -# Boost C++ Libraries, <em>Boost.Math / Statistical Distributions: Hyperexponential Distribution</em>, Online: http://www.boost.org/doc/libs/release/libs/math/doc/html/dist.html , 2014. | |
199 | * -# D.G. Feitelson, <em>Workload Modeling for Computer Systems Performance Evaluation</em>, Cambridge University Press, 2014 | |
200 | * -# A. Feldmann and W. Whitt, <em>Fitting mixtures of exponentials to long-tail distributions to analyze network performance models</em>, Performance Evaluation 31(3-4):245, doi:10.1016/S0166-5316(97)00003-5, 1998. | |
201 | * -# H.T. Papadopolous, C. Heavey and J. Browne, <em>Queueing Theory in Manufacturing Systems Analysis and Design</em>, Chapman & Hall/CRC, 1993, p. 35. | |
202 | * -# R.F. Rosin, <em>Determining a computing center environment</em>, Communications of the ACM 8(7):463-468, 1965. | |
203 | * -# K.S. Trivedi, <em>Probability and Statistics with Reliability, Queueing, and Computer Science Applications</em>, John Wiley & Sons, Inc., 2002. | |
204 | * -# J. Tyszer, <em>Object-Oriented Computer Simulation of Discrete-Event Systems</em>, Springer, 1999. | |
205 | * -# Wikipedia, <em>Hyperexponential Distribution</em>, Online: http://en.wikipedia.org/wiki/Hyperexponential_distribution , 2014. | |
206 | * -# Wolfram Mathematica, <em>Hyperexponential Distribution</em>, Online: http://reference.wolfram.com/language/ref/HyperexponentialDistribution.html , 2014. | |
207 | * . | |
208 | * | |
209 | * \author Marco Guazzone (marco.guazzone@gmail.com) | |
210 | */ | |
211 | template<class RealT = double> | |
212 | class hyperexponential_distribution | |
213 | { | |
214 | public: typedef RealT result_type; | |
215 | public: typedef RealT input_type; | |
216 | ||
217 | ||
218 | /** | |
219 | * The parameters of a hyperexponential distribution. | |
220 | * | |
221 | * Stores the <em>phase probability vector</em> and the <em>rate vector</em> | |
222 | * of the hyperexponential distribution. | |
223 | * | |
224 | * \author Marco Guazzone (marco.guazzone@gmail.com) | |
225 | */ | |
226 | public: class param_type | |
227 | { | |
228 | public: typedef hyperexponential_distribution distribution_type; | |
229 | ||
230 | /** | |
231 | * Constructs a \c param_type with the default parameters | |
232 | * of the distribution. | |
233 | */ | |
234 | public: param_type() | |
235 | : probs_(1, 1), | |
236 | rates_(1, 1) | |
237 | { | |
238 | } | |
239 | ||
240 | /** | |
241 | * Constructs a \c param_type from the <em>phase probability vector</em> | |
242 | * and <em>rate vector</em> parameters of the distribution. | |
243 | * | |
244 | * The <em>phase probability vector</em> parameter is given by the range | |
245 | * defined by [\a prob_first, \a prob_last) iterator pair, and the | |
246 | * <em>rate vector</em> parameter is given by the range defined by | |
247 | * [\a rate_first, \a rate_last) iterator pair. | |
248 | * | |
249 | * \tparam ProbIterT Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]). | |
250 | * \tparam RateIterT Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]). | |
251 | * | |
252 | * \param prob_first The iterator to the beginning of the range of non-negative real elements representing the phase probabilities; if elements don't sum to 1, they are normalized. | |
253 | * \param prob_last The iterator to the ending of the range of non-negative real elements representing the phase probabilities; if elements don't sum to 1, they are normalized. | |
254 | * \param rate_first The iterator to the beginning of the range of non-negative real elements representing the rates. | |
255 | * \param rate_last The iterator to the ending of the range of non-negative real elements representing the rates. | |
256 | * | |
257 | * References: | |
258 | * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014 | |
259 | * . | |
260 | */ | |
261 | public: template <typename ProbIterT, typename RateIterT> | |
262 | param_type(ProbIterT prob_first, ProbIterT prob_last, | |
263 | RateIterT rate_first, RateIterT rate_last) | |
264 | : probs_(prob_first, prob_last), | |
265 | rates_(rate_first, rate_last) | |
266 | { | |
267 | hyperexp_detail::normalize(probs_); | |
268 | ||
269 | assert( hyperexp_detail::check_params(probs_, rates_) ); | |
270 | } | |
271 | ||
272 | /** | |
273 | * Constructs a \c param_type from the <em>phase probability vector</em> | |
274 | * and <em>rate vector</em> parameters of the distribution. | |
275 | * | |
276 | * The <em>phase probability vector</em> parameter is given by the range | |
277 | * defined by \a prob_range, and the <em>rate vector</em> parameter is | |
278 | * given by the range defined by \a rate_range. | |
279 | * | |
280 | * \tparam ProbRangeT Must meet the requirements of <a href="boost:/libs/range/doc/html/range/concepts.html">Range</a> concept. | |
281 | * \tparam RateRangeT Must meet the requirements of <a href="boost:/libs/range/doc/html/range/concepts.html">Range</a> concept. | |
282 | * | |
283 | * \param prob_range The range of non-negative real elements representing the phase probabilities; if elements don't sum to 1, they are normalized. | |
284 | * \param rate_range The range of positive real elements representing the rates. | |
285 | * | |
286 | * \note | |
287 | * The final \c disable_if parameter is an implementation detail that | |
288 | * differentiates between this two argument constructor and the | |
289 | * iterator-based two argument constructor described below. | |
290 | */ | |
291 | // We SFINAE this out of existance if either argument type is | |
292 | // incrementable as in that case the type is probably an iterator: | |
293 | public: template <typename ProbRangeT, typename RateRangeT> | |
294 | param_type(ProbRangeT const& prob_range, | |
295 | RateRangeT const& rate_range, | |
296 | typename boost::disable_if_c<boost::has_pre_increment<ProbRangeT>::value || boost::has_pre_increment<RateRangeT>::value>::type* = 0) | |
297 | : probs_(boost::begin(prob_range), boost::end(prob_range)), | |
298 | rates_(boost::begin(rate_range), boost::end(rate_range)) | |
299 | { | |
300 | hyperexp_detail::normalize(probs_); | |
301 | ||
302 | assert( hyperexp_detail::check_params(probs_, rates_) ); | |
303 | } | |
304 | ||
305 | /** | |
306 | * Constructs a \c param_type from the <em>rate vector</em> parameter of | |
307 | * the distribution and with equal phase probabilities. | |
308 | * | |
309 | * The <em>rate vector</em> parameter is given by the range defined by | |
310 | * [\a rate_first, \a rate_last) iterator pair, and the <em>phase | |
311 | * probability vector</em> parameter is set to the equal phase | |
312 | * probabilities (i.e., to a vector of the same length \f$k\f$ of the | |
313 | * <em>rate vector</em> and with each element set to \f$1.0/k\f$). | |
314 | * | |
315 | * \tparam RateIterT Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]). | |
316 | * \tparam RateIterT2 Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]). | |
317 | * | |
318 | * \param rate_first The iterator to the beginning of the range of non-negative real elements representing the rates. | |
319 | * \param rate_last The iterator to the ending of the range of non-negative real elements representing the rates. | |
320 | * | |
321 | * \note | |
322 | * The final \c disable_if parameter is an implementation detail that | |
323 | * differentiates between this two argument constructor and the | |
324 | * range-based two argument constructor described above. | |
325 | * | |
326 | * References: | |
327 | * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014 | |
328 | * . | |
329 | */ | |
330 | // We SFINAE this out of existance if the argument type is | |
331 | // incrementable as in that case the type is probably an iterator. | |
332 | public: template <typename RateIterT> | |
333 | param_type(RateIterT rate_first, | |
334 | RateIterT rate_last, | |
335 | typename boost::enable_if_c<boost::has_pre_increment<RateIterT>::value>::type* = 0) | |
336 | : probs_(std::distance(rate_first, rate_last), 1), // will be normalized below | |
337 | rates_(rate_first, rate_last) | |
338 | { | |
339 | assert(probs_.size() == rates_.size()); | |
340 | } | |
341 | ||
342 | /** | |
343 | * Constructs a @c param_type from the "rates" parameters | |
344 | * of the distribution and with equal phase probabilities. | |
345 | * | |
346 | * The <em>rate vector</em> parameter is given by the range defined by | |
347 | * \a rate_range, and the <em>phase probability vector</em> parameter is | |
348 | * set to the equal phase probabilities (i.e., to a vector of the same | |
349 | * length \f$k\f$ of the <em>rate vector</em> and with each element set | |
350 | * to \f$1.0/k\f$). | |
351 | * | |
352 | * \tparam RateRangeT Must meet the requirements of <a href="boost:/libs/range/doc/html/range/concepts.html">Range</a> concept. | |
353 | * | |
354 | * \param rate_range The range of positive real elements representing the rates. | |
355 | */ | |
356 | public: template <typename RateRangeT> | |
357 | param_type(RateRangeT const& rate_range) | |
358 | : probs_(boost::size(rate_range), 1), // Will be normalized below | |
359 | rates_(boost::begin(rate_range), boost::end(rate_range)) | |
360 | { | |
361 | hyperexp_detail::normalize(probs_); | |
362 | ||
363 | assert( hyperexp_detail::check_params(probs_, rates_) ); | |
364 | } | |
365 | ||
366 | #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST | |
367 | /** | |
368 | * Constructs a \c param_type from the <em>phase probability vector</em> | |
369 | * and <em>rate vector</em> parameters of the distribution. | |
370 | * | |
371 | * The <em>phase probability vector</em> parameter is given by the | |
372 | * <em>brace-init-list</em> (ISO,2014,sec. 8.5.4 [dcl.init.list]) | |
373 | * defined by \a l1, and the <em>rate vector</em> parameter is given by the | |
374 | * <em>brace-init-list</em> (ISO,2014,sec. 8.5.4 [dcl.init.list]) | |
375 | * defined by \a l2. | |
376 | * | |
377 | * \param l1 The initializer list for inizializing the phase probability vector. | |
378 | * \param l2 The initializer list for inizializing the rate vector. | |
379 | * | |
380 | * References: | |
381 | * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014 | |
382 | * . | |
383 | */ | |
384 | public: param_type(std::initializer_list<RealT> l1, std::initializer_list<RealT> l2) | |
385 | : probs_(l1.begin(), l1.end()), | |
386 | rates_(l2.begin(), l2.end()) | |
387 | { | |
388 | hyperexp_detail::normalize(probs_); | |
389 | ||
390 | assert( hyperexp_detail::check_params(probs_, rates_) ); | |
391 | } | |
392 | ||
393 | /** | |
394 | * Constructs a \c param_type from the <em>rate vector</em> parameter | |
395 | * of the distribution and with equal phase probabilities. | |
396 | * | |
397 | * The <em>rate vector</em> parameter is given by the | |
398 | * <em>brace-init-list</em> (ISO,2014,sec. 8.5.4 [dcl.init.list]) | |
399 | * defined by \a l1, and the <em>phase probability vector</em> parameter is | |
400 | * set to the equal phase probabilities (i.e., to a vector of the same | |
401 | * length \f$k\f$ of the <em>rate vector</em> and with each element set | |
402 | * to \f$1.0/k\f$). | |
403 | * | |
404 | * \param l1 The initializer list for inizializing the rate vector. | |
405 | * | |
406 | * References: | |
407 | * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014 | |
408 | * . | |
409 | */ | |
410 | public: param_type(std::initializer_list<RealT> l1) | |
411 | : probs_(std::distance(l1.begin(), l1.end()), 1), // Will be normalized below | |
412 | rates_(l1.begin(), l1.end()) | |
413 | { | |
414 | hyperexp_detail::normalize(probs_); | |
415 | ||
416 | assert( hyperexp_detail::check_params(probs_, rates_) ); | |
417 | } | |
418 | #endif // BOOST_NO_CXX11_HDR_INITIALIZER_LIST | |
419 | ||
420 | /** | |
421 | * Gets the <em>phase probability vector</em> parameter of the distribtuion. | |
422 | * | |
423 | * \return The <em>phase probability vector</em> parameter of the distribution. | |
424 | * | |
425 | * \note | |
426 | * The returned probabilities are the normalized version of the ones | |
427 | * passed at construction time. | |
428 | */ | |
429 | public: std::vector<RealT> probabilities() const | |
430 | { | |
431 | return probs_; | |
432 | } | |
433 | ||
434 | /** | |
435 | * Gets the <em>rate vector</em> parameter of the distribtuion. | |
436 | * | |
437 | * \return The <em>rate vector</em> parameter of the distribution. | |
438 | */ | |
439 | public: std::vector<RealT> rates() const | |
440 | { | |
441 | return rates_; | |
442 | } | |
443 | ||
444 | /** Writes a \c param_type to a \c std::ostream. */ | |
445 | public: BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, param) | |
446 | { | |
447 | detail::print_vector(os, param.probs_); | |
448 | os << ' '; | |
449 | detail::print_vector(os, param.rates_); | |
450 | ||
451 | return os; | |
452 | } | |
453 | ||
454 | /** Reads a \c param_type from a \c std::istream. */ | |
455 | public: BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, param) | |
456 | { | |
457 | // NOTE: if \c std::ios_base::exceptions is set, the code below may | |
458 | // throw in case of a I/O failure. | |
459 | // To prevent leaving the state of \c param inconsistent: | |
460 | // - if an exception is thrown, the state of \c param is left | |
461 | // unchanged (i.e., is the same as the one at the beginning | |
462 | // of the function's execution), and | |
463 | // - the state of \c param only after reading the whole input. | |
464 | ||
465 | std::vector<RealT> probs; | |
466 | std::vector<RealT> rates; | |
467 | ||
468 | // Reads probability and rate vectors | |
469 | detail::read_vector(is, probs); | |
470 | if (!is) | |
471 | { | |
472 | return is; | |
473 | } | |
474 | is >> std::ws; | |
475 | detail::read_vector(is, rates); | |
476 | if (!is) | |
477 | { | |
478 | return is; | |
479 | } | |
480 | ||
481 | // Update the state of the param_type object | |
482 | if (probs.size() > 0) | |
483 | { | |
484 | param.probs_.swap(probs); | |
485 | probs.clear(); | |
486 | } | |
487 | if (rates.size() > 0) | |
488 | { | |
489 | param.rates_.swap(rates); | |
490 | rates.clear(); | |
491 | } | |
492 | ||
493 | bool fail = false; | |
494 | ||
495 | // Adjust vector sizes (if needed) | |
496 | if (param.probs_.size() != param.rates_.size() | |
497 | || param.probs_.size() == 0) | |
498 | { | |
499 | fail = true; | |
500 | ||
501 | const std::size_t np = param.probs_.size(); | |
502 | const std::size_t nr = param.rates_.size(); | |
503 | ||
504 | if (np > nr) | |
505 | { | |
506 | param.rates_.resize(np, 1); | |
507 | } | |
508 | else if (nr > np) | |
509 | { | |
510 | param.probs_.resize(nr, 1); | |
511 | } | |
512 | else | |
513 | { | |
514 | param.probs_.resize(1, 1); | |
515 | param.rates_.resize(1, 1); | |
516 | } | |
517 | } | |
518 | ||
519 | // Normalize probabilities | |
520 | // NOTE: this cannot be done earlier since the probability vector | |
521 | // can be changed due to size conformance | |
522 | hyperexp_detail::normalize(param.probs_); | |
523 | ||
524 | // Set the error state in the underlying stream in case of invalid input | |
525 | if (fail) | |
526 | { | |
527 | // This throws an exception if ios_base::exception(failbit) is enabled | |
528 | is.setstate(std::ios_base::failbit); | |
529 | } | |
530 | ||
531 | //post: vector size conformance | |
532 | assert(param.probs_.size() == param.rates_.size()); | |
533 | ||
534 | return is; | |
535 | } | |
536 | ||
537 | /** Returns true if the two sets of parameters are the same. */ | |
538 | public: BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs) | |
539 | { | |
540 | return lhs.probs_ == rhs.probs_ | |
541 | && lhs.rates_ == rhs.rates_; | |
542 | } | |
543 | ||
544 | /** Returns true if the two sets of parameters are the different. */ | |
545 | public: BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type) | |
546 | ||
547 | ||
548 | private: std::vector<RealT> probs_; ///< The <em>phase probability vector</em> parameter of the distribution | |
549 | private: std::vector<RealT> rates_; ///< The <em>rate vector</em> parameter of the distribution | |
550 | }; // param_type | |
551 | ||
552 | ||
553 | /** | |
554 | * Constructs a 1-phase \c hyperexponential_distribution (i.e., an | |
555 | * exponential distribution) with rate 1. | |
556 | */ | |
557 | public: hyperexponential_distribution() | |
558 | : dd_(std::vector<RealT>(1, 1)), | |
559 | rates_(1, 1) | |
560 | { | |
561 | // empty | |
562 | } | |
563 | ||
564 | /** | |
565 | * Constructs a \c hyperexponential_distribution from the <em>phase | |
566 | * probability vector</em> and <em>rate vector</em> parameters of the | |
567 | * distribution. | |
568 | * | |
569 | * The <em>phase probability vector</em> parameter is given by the range | |
570 | * defined by [\a prob_first, \a prob_last) iterator pair, and the | |
571 | * <em>rate vector</em> parameter is given by the range defined by | |
572 | * [\a rate_first, \a rate_last) iterator pair. | |
573 | * | |
574 | * \tparam ProbIterT Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]). | |
575 | * \tparam RateIterT Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]). | |
576 | * | |
577 | * \param prob_first The iterator to the beginning of the range of non-negative real elements representing the phase probabilities; if elements don't sum to 1, they are normalized. | |
578 | * \param prob_last The iterator to the ending of the range of non-negative real elements representing the phase probabilities; if elements don't sum to 1, they are normalized. | |
579 | * \param rate_first The iterator to the beginning of the range of non-negative real elements representing the rates. | |
580 | * \param rate_last The iterator to the ending of the range of non-negative real elements representing the rates. | |
581 | * | |
582 | * References: | |
583 | * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014 | |
584 | * . | |
585 | */ | |
586 | public: template <typename ProbIterT, typename RateIterT> | |
587 | hyperexponential_distribution(ProbIterT prob_first, ProbIterT prob_last, | |
588 | RateIterT rate_first, RateIterT rate_last) | |
589 | : dd_(prob_first, prob_last), | |
590 | rates_(rate_first, rate_last) | |
591 | { | |
592 | assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) ); | |
593 | } | |
594 | ||
595 | /** | |
596 | * Constructs a \c hyperexponential_distribution from the <em>phase | |
597 | * probability vector</em> and <em>rate vector</em> parameters of the | |
598 | * distribution. | |
599 | * | |
600 | * The <em>phase probability vector</em> parameter is given by the range | |
601 | * defined by \a prob_range, and the <em>rate vector</em> parameter is | |
602 | * given by the range defined by \a rate_range. | |
603 | * | |
604 | * \tparam ProbRangeT Must meet the requirements of <a href="boost:/libs/range/doc/html/range/concepts.html">Range</a> concept. | |
605 | * \tparam RateRangeT Must meet the requirements of <a href="boost:/libs/range/doc/html/range/concepts.html">Range</a> concept. | |
606 | * | |
607 | * \param prob_range The range of non-negative real elements representing the phase probabilities; if elements don't sum to 1, they are normalized. | |
608 | * \param rate_range The range of positive real elements representing the rates. | |
609 | * | |
610 | * \note | |
611 | * The final \c disable_if parameter is an implementation detail that | |
612 | * differentiates between this two argument constructor and the | |
613 | * iterator-based two argument constructor described below. | |
614 | */ | |
615 | // We SFINAE this out of existance if either argument type is | |
616 | // incrementable as in that case the type is probably an iterator: | |
617 | public: template <typename ProbRangeT, typename RateRangeT> | |
618 | hyperexponential_distribution(ProbRangeT const& prob_range, | |
619 | RateRangeT const& rate_range, | |
620 | typename boost::disable_if_c<boost::has_pre_increment<ProbRangeT>::value || boost::has_pre_increment<RateRangeT>::value>::type* = 0) | |
621 | : dd_(prob_range), | |
622 | rates_(boost::begin(rate_range), boost::end(rate_range)) | |
623 | { | |
624 | assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) ); | |
625 | } | |
626 | ||
627 | /** | |
628 | * Constructs a \c hyperexponential_distribution from the <em>rate | |
629 | * vector</em> parameter of the distribution and with equal phase | |
630 | * probabilities. | |
631 | * | |
632 | * The <em>rate vector</em> parameter is given by the range defined by | |
633 | * [\a rate_first, \a rate_last) iterator pair, and the <em>phase | |
634 | * probability vector</em> parameter is set to the equal phase | |
635 | * probabilities (i.e., to a vector of the same length \f$k\f$ of the | |
636 | * <em>rate vector</em> and with each element set to \f$1.0/k\f$). | |
637 | * | |
638 | * \tparam RateIterT Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]). | |
639 | * \tparam RateIterT2 Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]). | |
640 | * | |
641 | * \param rate_first The iterator to the beginning of the range of non-negative real elements representing the rates. | |
642 | * \param rate_last The iterator to the ending of the range of non-negative real elements representing the rates. | |
643 | * | |
644 | * \note | |
645 | * The final \c disable_if parameter is an implementation detail that | |
646 | * differentiates between this two argument constructor and the | |
647 | * range-based two argument constructor described above. | |
648 | * | |
649 | * References: | |
650 | * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014 | |
651 | * . | |
652 | */ | |
653 | // We SFINAE this out of existance if the argument type is | |
654 | // incrementable as in that case the type is probably an iterator. | |
655 | public: template <typename RateIterT> | |
656 | hyperexponential_distribution(RateIterT rate_first, | |
657 | RateIterT rate_last, | |
658 | typename boost::enable_if_c<boost::has_pre_increment<RateIterT>::value>::type* = 0) | |
659 | : dd_(std::vector<RealT>(std::distance(rate_first, rate_last), 1)), | |
660 | rates_(rate_first, rate_last) | |
661 | { | |
662 | assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) ); | |
663 | } | |
664 | ||
665 | /** | |
666 | * Constructs a @c param_type from the "rates" parameters | |
667 | * of the distribution and with equal phase probabilities. | |
668 | * | |
669 | * The <em>rate vector</em> parameter is given by the range defined by | |
670 | * \a rate_range, and the <em>phase probability vector</em> parameter is | |
671 | * set to the equal phase probabilities (i.e., to a vector of the same | |
672 | * length \f$k\f$ of the <em>rate vector</em> and with each element set | |
673 | * to \f$1.0/k\f$). | |
674 | * | |
675 | * \tparam RateRangeT Must meet the requirements of <a href="boost:/libs/range/doc/html/range/concepts.html">Range</a> concept. | |
676 | * | |
677 | * \param rate_range The range of positive real elements representing the rates. | |
678 | */ | |
679 | public: template <typename RateRangeT> | |
680 | hyperexponential_distribution(RateRangeT const& rate_range) | |
681 | : dd_(std::vector<RealT>(boost::size(rate_range), 1)), | |
682 | rates_(boost::begin(rate_range), boost::end(rate_range)) | |
683 | { | |
684 | assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) ); | |
685 | } | |
686 | ||
687 | /** | |
688 | * Constructs a \c hyperexponential_distribution from its parameters. | |
689 | * | |
690 | * \param param The parameters of the distribution. | |
691 | */ | |
692 | public: explicit hyperexponential_distribution(param_type const& param) | |
693 | : dd_(param.probabilities()), | |
694 | rates_(param.rates()) | |
695 | { | |
696 | assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) ); | |
697 | } | |
698 | ||
699 | #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST | |
700 | /** | |
701 | * Constructs a \c hyperexponential_distribution from the <em>phase | |
702 | * probability vector</em> and <em>rate vector</em> parameters of the | |
703 | * distribution. | |
704 | * | |
705 | * The <em>phase probability vector</em> parameter is given by the | |
706 | * <em>brace-init-list</em> (ISO,2014,sec. 8.5.4 [dcl.init.list]) | |
707 | * defined by \a l1, and the <em>rate vector</em> parameter is given by the | |
708 | * <em>brace-init-list</em> (ISO,2014,sec. 8.5.4 [dcl.init.list]) | |
709 | * defined by \a l2. | |
710 | * | |
711 | * \param l1 The initializer list for inizializing the phase probability vector. | |
712 | * \param l2 The initializer list for inizializing the rate vector. | |
713 | * | |
714 | * References: | |
715 | * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014 | |
716 | * . | |
717 | */ | |
718 | public: hyperexponential_distribution(std::initializer_list<RealT> const& l1, std::initializer_list<RealT> const& l2) | |
719 | : dd_(l1.begin(), l1.end()), | |
720 | rates_(l2.begin(), l2.end()) | |
721 | { | |
722 | assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) ); | |
723 | } | |
724 | ||
725 | /** | |
726 | * Constructs a \c hyperexponential_distribution from the <em>rate | |
727 | * vector</em> parameter of the distribution and with equal phase | |
728 | * probabilities. | |
729 | * | |
730 | * The <em>rate vector</em> parameter is given by the | |
731 | * <em>brace-init-list</em> (ISO,2014,sec. 8.5.4 [dcl.init.list]) | |
732 | * defined by \a l1, and the <em>phase probability vector</em> parameter is | |
733 | * set to the equal phase probabilities (i.e., to a vector of the same | |
734 | * length \f$k\f$ of the <em>rate vector</em> and with each element set | |
735 | * to \f$1.0/k\f$). | |
736 | * | |
737 | * \param l1 The initializer list for inizializing the rate vector. | |
738 | * | |
739 | * References: | |
740 | * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014 | |
741 | * . | |
742 | */ | |
743 | public: hyperexponential_distribution(std::initializer_list<RealT> const& l1) | |
744 | : dd_(std::vector<RealT>(std::distance(l1.begin(), l1.end()), 1)), | |
745 | rates_(l1.begin(), l1.end()) | |
746 | { | |
747 | assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) ); | |
748 | } | |
749 | #endif | |
750 | ||
751 | /** | |
752 | * Gets a random variate distributed according to the | |
753 | * hyperexponential distribution. | |
754 | * | |
755 | * \tparam URNG Must meet the requirements of \uniform_random_number_generator. | |
756 | * | |
757 | * \param urng A uniform random number generator object. | |
758 | * | |
759 | * \return A random variate distributed according to the hyperexponential distribution. | |
760 | */ | |
761 | public: template<class URNG>\ | |
762 | RealT operator()(URNG& urng) const | |
763 | { | |
764 | const int i = dd_(urng); | |
765 | ||
766 | return boost::random::exponential_distribution<RealT>(rates_[i])(urng); | |
767 | } | |
768 | ||
769 | /** | |
770 | * Gets a random variate distributed according to the hyperexponential | |
771 | * distribution with parameters specified by \c param. | |
772 | * | |
773 | * \tparam URNG Must meet the requirements of \uniform_random_number_generator. | |
774 | * | |
775 | * \param urng A uniform random number generator object. | |
776 | * \param param A distribution parameter object. | |
777 | * | |
778 | * \return A random variate distributed according to the hyperexponential distribution. | |
779 | * distribution with parameters specified by \c param. | |
780 | */ | |
781 | public: template<class URNG> | |
782 | RealT operator()(URNG& urng, const param_type& param) const | |
783 | { | |
784 | return hyperexponential_distribution(param)(urng); | |
785 | } | |
786 | ||
787 | /** Returns the number of phases of the distribution. */ | |
788 | public: std::size_t num_phases() const | |
789 | { | |
790 | return rates_.size(); | |
791 | } | |
792 | ||
793 | /** Returns the <em>phase probability vector</em> parameter of the distribution. */ | |
794 | public: std::vector<RealT> probabilities() const | |
795 | { | |
796 | return dd_.probabilities(); | |
797 | } | |
798 | ||
799 | /** Returns the <em>rate vector</em> parameter of the distribution. */ | |
800 | public: std::vector<RealT> rates() const | |
801 | { | |
802 | return rates_; | |
803 | } | |
804 | ||
805 | /** Returns the smallest value that the distribution can produce. */ | |
806 | public: RealT min BOOST_PREVENT_MACRO_SUBSTITUTION () const | |
807 | { | |
808 | return 0; | |
809 | } | |
810 | ||
811 | /** Returns the largest value that the distribution can produce. */ | |
812 | public: RealT max BOOST_PREVENT_MACRO_SUBSTITUTION () const | |
813 | { | |
814 | return std::numeric_limits<RealT>::infinity(); | |
815 | } | |
816 | ||
817 | /** Returns the parameters of the distribution. */ | |
818 | public: param_type param() const | |
819 | { | |
820 | std::vector<RealT> probs = dd_.probabilities(); | |
821 | ||
822 | return param_type(probs.begin(), probs.end(), rates_.begin(), rates_.end()); | |
823 | } | |
824 | ||
825 | /** Sets the parameters of the distribution. */ | |
826 | public: void param(param_type const& param) | |
827 | { | |
828 | dd_.param(typename boost::random::discrete_distribution<int,RealT>::param_type(param.probabilities())); | |
829 | rates_ = param.rates(); | |
830 | } | |
831 | ||
832 | /** | |
833 | * Effects: Subsequent uses of the distribution do not depend | |
834 | * on values produced by any engine prior to invoking reset. | |
835 | */ | |
836 | public: void reset() | |
837 | { | |
838 | // empty | |
839 | } | |
840 | ||
841 | /** Writes an @c hyperexponential_distribution to a @c std::ostream. */ | |
842 | public: BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, hyperexponential_distribution, hd) | |
843 | { | |
844 | os << hd.param(); | |
845 | return os; | |
846 | } | |
847 | ||
848 | /** Reads an @c hyperexponential_distribution from a @c std::istream. */ | |
849 | public: BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, hyperexponential_distribution, hd) | |
850 | { | |
851 | param_type param; | |
852 | if(is >> param) | |
853 | { | |
854 | hd.param(param); | |
855 | } | |
856 | return is; | |
857 | } | |
858 | ||
859 | /** | |
860 | * Returns true if the two instances of @c hyperexponential_distribution will | |
861 | * return identical sequences of values given equal generators. | |
862 | */ | |
863 | public: BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(hyperexponential_distribution, lhs, rhs) | |
864 | { | |
865 | return lhs.dd_ == rhs.dd_ | |
866 | && lhs.rates_ == rhs.rates_; | |
867 | } | |
868 | ||
869 | /** | |
870 | * Returns true if the two instances of @c hyperexponential_distribution will | |
871 | * return different sequences of values given equal generators. | |
872 | */ | |
873 | public: BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(hyperexponential_distribution) | |
874 | ||
875 | ||
876 | private: boost::random::discrete_distribution<int,RealT> dd_; ///< The \c discrete_distribution used to sample the phase probability and choose the rate | |
877 | private: std::vector<RealT> rates_; ///< The <em>rate vector</em> parameter of the distribution | |
878 | }; // hyperexponential_distribution | |
879 | ||
880 | }} // namespace boost::random | |
881 | ||
882 | ||
883 | #endif // BOOST_RANDOM_HYPEREXPONENTIAL_DISTRIBUTION_HPP |