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
44namespace boost { namespace random {
45
46namespace hyperexp_detail {
47
48template <typename T>
49std::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
72template <typename RealT>
73bool 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
105template <typename RealT>
106bool 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
120template <typename RealT>
121bool 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 */
211template<class RealT = double>
212class 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
884