boost/random/hyperexponential_distribution.hpp

0001 /* boost random/hyperexponential_distribution.hpp header file
0002  *
0003  * Copyright Marco Guazzone 2014
0004  * Distributed under the Boost Software License, Version 1.0. (See
0005  * accompanying file LICENSE_1_0.txt or copy at
0006  * http://www.boost.org/LICENSE_1_0.txt)
0007  *
0008  * See http://www.boost.org for most recent version including documentation.
0009  *
0010  * Much of the code here taken by boost::math::hyperexponential_distribution.
0011  * To this end, we would like to thank Paul Bristow and John Maddock for their
0012  * valuable feedback.
0013  *
0014  * \author Marco Guazzone (marco.guazzone@gmail.com)
0015  */
0016
0017 #ifndef BOOST_RANDOM_HYPEREXPONENTIAL_DISTRIBUTION_HPP
0018 #define BOOST_RANDOM_HYPEREXPONENTIAL_DISTRIBUTION_HPP
0019
0020
0021 #include <boost/config.hpp>
0022 #include <boost/core/cmath.hpp>
0023 #include <boost/random/detail/operators.hpp>
0024 #include <boost/random/detail/vector_io.hpp>
0025 #include <boost/random/discrete_distribution.hpp>
0026 #include <boost/random/exponential_distribution.hpp>
0027 #include <boost/range/begin.hpp>
0028 #include <boost/range/end.hpp>
0029 #include <boost/range/size.hpp>
0030 #include <boost/type_traits/has_pre_increment.hpp>
0031 #include <cassert>
0032 #include <cmath>
0033 #include <cstddef>
0034 #include <iterator>
0035 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
0036 # include <initializer_list>
0037 #endif // BOOST_NO_CXX11_HDR_INITIALIZER_LIST
0038 #include <iostream>
0039 #include <limits>
0040 #include <numeric>
0041 #include <vector>
0042
0043
0044 namespace boost { namespace random {
0045
0046 namespace hyperexp_detail {
0047
0048 template <typename T>
0049 std::vector<T>& normalize(std::vector<T>& v)
0050 {
0051     if (v.size() == 0)
0052     {
0053         return v;
0054     }
0055
0056     const T sum = std::accumulate(v.begin(), v.end(), static_cast<T>(0));
0057     T final_sum = 0;
0058
0059     const typename std::vector<T>::iterator end = --v.end();
0060     for (typename std::vector<T>::iterator it = v.begin();
0061          it != end;
0062          ++it)
0063     {
0064         *it /= sum;
0065         final_sum += *it;
0066     }
0067     *end = 1-final_sum; // avoids round off errors thus ensuring the probabilities really sum to 1
0068
0069     return v;
0070 }
0071
0072 template <typename RealT>
0073 bool check_probabilities(std::vector<RealT> const& probabilities)
0074 {
0075     const std::size_t n = probabilities.size();
0076     RealT sum = 0;
0077     for (std::size_t i = 0; i < n; ++i)
0078     {
0079         if (probabilities[i] < 0
0080             || probabilities[i] > 1
0081             || !(boost::core::isfinite)(probabilities[i]))
0082         {
0083             return false;
0084         }
0085         sum += probabilities[i];
0086     }
0087
0088     //NOTE: the check below seems to fail on some architectures.
0089     //      So we commented it.
0090     //// - We try to keep phase probabilities correctly normalized in the distribution constructors
0091     //// - However in practice we have to allow for a very slight divergence from a sum of exactly 1:
0092     ////if (std::abs(sum-1) > (std::numeric_limits<RealT>::epsilon()*2))
0093     //// This is from Knuth "The Art of Computer Programming: Vol.2, 3rd Ed", and can be used to
0094     //// check is two numbers are approximately equal
0095     //const RealT one = 1;
0096     //const RealT tol = std::numeric_limits<RealT>::epsilon()*2.0;
0097     //if (std::abs(sum-one) > (std::max(std::abs(sum), std::abs(one))*tol))
0098     //{
0099     //    return false;
0100     //}
0101
0102     return true;
0103 }
0104
0105 template <typename RealT>
0106 bool check_rates(std::vector<RealT> const& rates)
0107 {
0108     const std::size_t n = rates.size();
0109     for (std::size_t i = 0; i < n; ++i)
0110     {
0111         if (rates[i] <= 0
0112             || !(boost::core::isfinite)(rates[i]))
0113         {
0114             return false;
0115         }
0116     }
0117     return true;
0118 }
0119
0120 template <typename RealT>
0121 bool check_params(std::vector<RealT> const& probabilities, std::vector<RealT> const& rates)
0122 {
0123     if (probabilities.size() != rates.size())
0124     {
0125         return false;
0126     }
0127
0128     return check_probabilities(probabilities)
0129            && check_rates(rates);
0130 }
0131
0132 } // Namespace hyperexp_detail
0133
0134
0135 /**
0136  * The hyperexponential distribution is a real-valued continuous distribution
0137  * with two parameters, the <em>phase probability vector</em> \c probs and the
0138  * <em>rate vector</em> \c rates.
0139  *
0140  * A \f$k\f$-phase hyperexponential distribution is a mixture of \f$k\f$
0141  * exponential distributions.
0142  * For this reason, it is also referred to as <em>mixed exponential
0143  * distribution</em> or <em>parallel \f$k\f$-phase exponential
0144  * distribution</em>.
0145  *
0146  * A \f$k\f$-phase hyperexponential distribution is characterized by two
0147  * 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$.
0148  *
0149  * A \f$k\f$-phase hyperexponential distribution is frequently used in
0150  * <em>queueing theory</em> to model the distribution of the superposition of
0151  * \f$k\f$ independent events, like, for instance, the  service time distribution
0152  * of a queueing station with \f$k\f$ servers in parallel where the \f$i\f$-th
0153  * server is chosen with probability \f$\alpha_i\f$ and its service time
0154  * distribution is an exponential distribution with rate \f$\lambda_i\f$
0155  * (Allen,1990; Papadopolous et al.,1993; Trivedi,2002).
0156  *
0157  * For instance, CPUs service-time distribution in a computing system has often
0158  * been observed to possess such a distribution (Rosin,1965).
0159  * Also, the arrival of different types of customer to a single queueing station
0160  * is often modeled as a hyperexponential distribution (Papadopolous et al.,1993).
0161  * Similarly, if a product manufactured in several parallel assemply lines and
0162  * the outputs are merged, the failure density of the overall product is likely
0163  * to be hyperexponential (Trivedi,2002).
0164  *
0165  * Finally, since the hyperexponential distribution exhibits a high Coefficient
0166  * of Variation (CoV), that is a CoV > 1, it is especially suited to fit
0167  * empirical data with large CoV (Feitelson,2014; Wolski et al.,2013) and to
0168  * approximate <em>long-tail probability distributions</em> (Feldmann et al.,1998).
0169  *
0170  * See (Boost,2014) for more information and examples.
0171  *
0172  * A \f$k\f$-phase hyperexponential distribution has a probability density
0173  * function
0174  * \f[
0175  *  f(x) = \sum_{i=1}^k \alpha_i \lambda_i e^{-x\lambda_i}
0176  * \f]
0177  * where:
0178  * - \f$k\f$ is the <em>number of phases</em> and also the size of the input
0179  *   vector parameters,
0180  * - \f$\mathbf{\alpha}=(\alpha_1,\ldots,\alpha_k)\f$ is the <em>phase probability
0181  *   vector</em> parameter, and
0182  * - \f$\mathbf{\lambda}=(\lambda_1,\ldots,\lambda_k)\f$ is the <em>rate vector</em>
0183  *   parameter.
0184  * .
0185  *
0186  * Given a \f$k\f$-phase hyperexponential distribution with phase probability
0187  * vector \f$\mathbf{\alpha}\f$ and rate vector \f$\mathbf{\lambda}\f$, the
0188  * random variate generation algorithm consists of the following steps (Tyszer,1999):
0189  * -# Generate a random variable \f$U\f$ uniformly distribution on the interval \f$(0,1)\f$.
0190  * -# Use \f$U\f$ to select the appropriate \f$\lambda_i\f$ (e.g., the
0191  *  <em>alias method</em> can possibly be used for this step).
0192  * -# Generate an exponentially distributed random variable \f$X\f$ with rate parameter \f$\lambda_i\f$.
0193  * -# Return \f$X\f$.
0194  * .
0195  *
0196  * References:
0197  * -# A.O. Allen, <em>Probability, Statistics, and Queuing Theory with Computer Science Applications, Second Edition</em>, Academic Press, 1990.
0198  * -# 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.
0199  * -# D.G. Feitelson, <em>Workload Modeling for Computer Systems Performance Evaluation</em>, Cambridge University Press, 2014
0200  * -# 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.
0201  * -# H.T. Papadopolous, C. Heavey and J. Browne, <em>Queueing Theory in Manufacturing Systems Analysis and Design</em>, Chapman & Hall/CRC, 1993, p. 35.
0202  * -# R.F. Rosin, <em>Determining a computing center environment</em>, Communications of the ACM 8(7):463-468, 1965.
0203  * -# K.S. Trivedi, <em>Probability and Statistics with Reliability, Queueing, and Computer Science Applications</em>, John Wiley & Sons, Inc., 2002.
0204  * -# J. Tyszer, <em>Object-Oriented Computer Simulation of Discrete-Event Systems</em>, Springer, 1999.
0205  * -# Wikipedia, <em>Hyperexponential Distribution</em>, Online: http://en.wikipedia.org/wiki/Hyperexponential_distribution , 2014.
0206  * -# Wolfram Mathematica, <em>Hyperexponential Distribution</em>, Online: http://reference.wolfram.com/language/ref/HyperexponentialDistribution.html , 2014.
0207  * .
0208  *
0209  * \author Marco Guazzone (marco.guazzone@gmail.com)
0210  */
0211 template<class RealT = double>
0212 class hyperexponential_distribution
0213 {
0214     public: typedef RealT result_type;
0215     public: typedef RealT input_type;
0216
0217
0218     /**
0219      * The parameters of a hyperexponential distribution.
0220      *
0221      * Stores the <em>phase probability vector</em> and the <em>rate vector</em>
0222      * of the hyperexponential distribution.
0223      *
0224      * \author Marco Guazzone (marco.guazzone@gmail.com)
0225      */
0226     public: class param_type
0227     {
0228         public: typedef hyperexponential_distribution distribution_type;
0229
0230         /**
0231          * Constructs a \c param_type with the default parameters
0232          * of the distribution.
0233          */
0234         public: param_type()
0235         : probs_(1, 1),
0236           rates_(1, 1)
0237         {
0238         }
0239
0240         /**
0241          * Constructs a \c param_type from the <em>phase probability vector</em>
0242          * and <em>rate vector</em> parameters of the distribution.
0243          *
0244          * The <em>phase probability vector</em> parameter is given by the range
0245          * defined by [\a prob_first, \a prob_last) iterator pair, and the
0246          * <em>rate vector</em> parameter is given by the range defined by
0247          * [\a rate_first, \a rate_last) iterator pair.
0248          *
0249          * \tparam ProbIterT Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]).
0250          * \tparam RateIterT Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]).
0251          *
0252          * \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.
0253          * \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.
0254          * \param rate_first The iterator to the beginning of the range of non-negative real elements representing the rates.
0255          * \param rate_last The iterator to the ending of the range of non-negative real elements representing the rates.
0256          *
0257          * References:
0258          * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014
0259          * .
0260          */
0261         public: template <typename ProbIterT, typename RateIterT>
0262                 param_type(ProbIterT prob_first, ProbIterT prob_last,
0263                            RateIterT rate_first, RateIterT rate_last)
0264         : probs_(prob_first, prob_last),
0265           rates_(rate_first, rate_last)
0266         {
0267             hyperexp_detail::normalize(probs_);
0268
0269             assert( hyperexp_detail::check_params(probs_, rates_) );
0270         }
0271
0272         /**
0273          * Constructs a \c param_type from the <em>phase probability vector</em>
0274          * and <em>rate vector</em> parameters of the distribution.
0275          *
0276          * The <em>phase probability vector</em> parameter is given by the range
0277          * defined by \a prob_range, and the <em>rate vector</em> parameter is
0278          * given by the range defined by \a rate_range.
0279          *
0280          * \tparam ProbRangeT Must meet the requirements of <a href="boost:/libs/range/doc/html/range/concepts.html">Range</a> concept.
0281          * \tparam RateRangeT Must meet the requirements of <a href="boost:/libs/range/doc/html/range/concepts.html">Range</a> concept.
0282          *
0283          * \param prob_range The range of non-negative real elements representing the phase probabilities; if elements don't sum to 1, they are normalized.
0284          * \param rate_range The range of positive real elements representing the rates.
0285          *
0286          * \note
0287          *  The final \c disable_if parameter is an implementation detail that
0288          *  differentiates between this two argument constructor and the
0289          *  iterator-based two argument constructor described below.
0290          */
0291         //  We SFINAE this out of existance if either argument type is
0292         //  incrementable as in that case the type is probably an iterator:
0293         public: template <typename ProbRangeT, typename RateRangeT>
0294                 param_type(ProbRangeT const& prob_range,
0295                            RateRangeT const& rate_range,
0296                            typename boost::disable_if_c<boost::has_pre_increment<ProbRangeT>::value || boost::has_pre_increment<RateRangeT>::value>::type* = 0)
0297         : probs_(boost::begin(prob_range), boost::end(prob_range)),
0298           rates_(boost::begin(rate_range), boost::end(rate_range))
0299         {
0300             hyperexp_detail::normalize(probs_);
0301
0302             assert( hyperexp_detail::check_params(probs_, rates_) );
0303         }
0304
0305         /**
0306          * Constructs a \c param_type from the <em>rate vector</em> parameter of
0307          * the distribution and with equal phase probabilities.
0308          *
0309          * The <em>rate vector</em> parameter is given by the range defined by
0310          * [\a rate_first, \a rate_last) iterator pair, and the <em>phase
0311          * probability vector</em> parameter is set to the equal phase
0312          * probabilities (i.e., to a vector of the same length \f$k\f$ of the
0313          * <em>rate vector</em> and with each element set to \f$1.0/k\f$).
0314          *
0315          * \tparam RateIterT Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]).
0316          * \tparam RateIterT2 Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]).
0317          *
0318          * \param rate_first The iterator to the beginning of the range of non-negative real elements representing the rates.
0319          * \param rate_last The iterator to the ending of the range of non-negative real elements representing the rates.
0320          *
0321          * \note
0322          *  The final \c disable_if parameter is an implementation detail that
0323          *  differentiates between this two argument constructor and the
0324          *  range-based two argument constructor described above.
0325          *
0326          * References:
0327          * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014
0328          * .
0329          */
0330         //  We SFINAE this out of existance if the argument type is
0331         //  incrementable as in that case the type is probably an iterator.
0332         public: template <typename RateIterT>
0333                 param_type(RateIterT rate_first,
0334                            RateIterT rate_last,
0335                            typename boost::enable_if_c<boost::has_pre_increment<RateIterT>::value>::type* = 0)
0336         : probs_(std::distance(rate_first, rate_last), 1), // will be normalized below
0337           rates_(rate_first, rate_last)
0338         {
0339             assert(probs_.size() == rates_.size());
0340         }
0341
0342         /**
0343          * Constructs a @c param_type from the "rates" parameters
0344          * of the distribution and with equal phase probabilities.
0345          *
0346          * The <em>rate vector</em> parameter is given by the range defined by
0347          * \a rate_range, and the <em>phase probability vector</em> parameter is
0348          * set to the equal phase probabilities (i.e., to a vector of the same
0349          * length \f$k\f$ of the <em>rate vector</em> and with each element set
0350          * to \f$1.0/k\f$).
0351          *
0352          * \tparam RateRangeT Must meet the requirements of <a href="boost:/libs/range/doc/html/range/concepts.html">Range</a> concept.
0353          *
0354          * \param rate_range The range of positive real elements representing the rates.
0355          */
0356         public: template <typename RateRangeT>
0357                 param_type(RateRangeT const& rate_range)
0358         : probs_(boost::size(rate_range), 1), // Will be normalized below
0359           rates_(boost::begin(rate_range), boost::end(rate_range))
0360         {
0361             hyperexp_detail::normalize(probs_);
0362
0363             assert( hyperexp_detail::check_params(probs_, rates_) );
0364         }
0365
0366 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
0367         /**
0368          * Constructs a \c param_type from the <em>phase probability vector</em>
0369          * and <em>rate vector</em> parameters of the distribution.
0370          *
0371          * The <em>phase probability vector</em> parameter is given by the
0372          * <em>brace-init-list</em> (ISO,2014,sec. 8.5.4 [dcl.init.list])
0373          * defined by \a l1, and the <em>rate vector</em> parameter is given by the
0374          * <em>brace-init-list</em> (ISO,2014,sec. 8.5.4 [dcl.init.list])
0375          * defined by \a l2.
0376          *
0377          * \param l1 The initializer list for inizializing the phase probability vector.
0378          * \param l2 The initializer list for inizializing the rate vector.
0379          *
0380          * References:
0381          * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014
0382          * .
0383          */
0384         public: param_type(std::initializer_list<RealT> l1, std::initializer_list<RealT> l2)
0385         : probs_(l1.begin(), l1.end()),
0386           rates_(l2.begin(), l2.end())
0387         {
0388             hyperexp_detail::normalize(probs_);
0389
0390             assert( hyperexp_detail::check_params(probs_, rates_) );
0391         }
0392
0393         /**
0394          * Constructs a \c param_type from the <em>rate vector</em> parameter
0395          * of the distribution and with equal phase probabilities.
0396          *
0397          * The <em>rate vector</em> parameter is given by the
0398          * <em>brace-init-list</em> (ISO,2014,sec. 8.5.4 [dcl.init.list])
0399          * defined by \a l1, and the <em>phase probability vector</em> parameter is
0400          * set to the equal phase probabilities (i.e., to a vector of the same
0401          * length \f$k\f$ of the <em>rate vector</em> and with each element set
0402          * to \f$1.0/k\f$).
0403          *
0404          * \param l1 The initializer list for inizializing the rate vector.
0405          *
0406          * References:
0407          * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014
0408          * .
0409          */
0410         public: param_type(std::initializer_list<RealT> l1)
0411         : probs_(std::distance(l1.begin(), l1.end()), 1), // Will be normalized below
0412           rates_(l1.begin(), l1.end())
0413         {
0414             hyperexp_detail::normalize(probs_);
0415
0416             assert( hyperexp_detail::check_params(probs_, rates_) );
0417         }
0418 #endif // BOOST_NO_CXX11_HDR_INITIALIZER_LIST
0419
0420         /**
0421          * Gets the <em>phase probability vector</em> parameter of the distribtuion.
0422          *
0423          * \return The <em>phase probability vector</em> parameter of the distribution.
0424          *
0425          * \note
0426          *  The returned probabilities are the normalized version of the ones
0427          *  passed at construction time.
0428          */
0429         public: std::vector<RealT> probabilities() const
0430         {
0431             return probs_;
0432         }
0433
0434         /**
0435          * Gets the <em>rate vector</em> parameter of the distribtuion.
0436          *
0437          * \return The <em>rate vector</em> parameter of the distribution.
0438          */
0439         public: std::vector<RealT> rates() const
0440         {
0441             return rates_;
0442         }
0443
0444         /** Writes a \c param_type to a \c std::ostream. */
0445         public: BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, param)
0446         {
0447             detail::print_vector(os, param.probs_);
0448             os << ' ';
0449             detail::print_vector(os, param.rates_);
0450
0451             return os;
0452         }
0453
0454         /** Reads a \c param_type from a \c std::istream. */
0455         public: BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, param)
0456         {
0457             // NOTE: if \c std::ios_base::exceptions is set, the code below may
0458             //       throw in case of a I/O failure.
0459             //       To prevent leaving the state of \c param inconsistent:
0460             //       - if an exception is thrown, the state of \c param is left
0461             //         unchanged (i.e., is the same as the one at the beginning
0462             //         of the function's execution), and
0463             //       - the state of \c param only after reading the whole input.
0464
0465             std::vector<RealT> probs;
0466             std::vector<RealT> rates;
0467
0468             // Reads probability and rate vectors
0469             detail::read_vector(is, probs);
0470             if (!is)
0471             {
0472                 return is;
0473             }
0474             is >> std::ws;
0475             detail::read_vector(is, rates);
0476             if (!is)
0477             {
0478                 return is;
0479             }
0480
0481             // Update the state of the param_type object
0482             if (probs.size() > 0)
0483             {
0484                 param.probs_.swap(probs);
0485                 probs.clear();
0486             }
0487             if (rates.size() > 0)
0488             {
0489                 param.rates_.swap(rates);
0490                 rates.clear();
0491             }
0492
0493             // Adjust vector sizes (if needed)
0494             if (param.probs_.size() != param.rates_.size()
0495                 || param.probs_.size() == 0)
0496             {
0497                 const std::size_t np = param.probs_.size();
0498                 const std::size_t nr = param.rates_.size();
0499
0500                 if (np > nr)
0501                 {
0502                     param.rates_.resize(np, 1);
0503                 }
0504                 else if (nr > np)
0505                 {
0506                     param.probs_.resize(nr, 1);
0507                 }
0508                 else
0509                 {
0510                     param.probs_.resize(1, 1);
0511                     param.rates_.resize(1, 1);
0512                 }
0513             }
0514
0515             // Normalize probabilities
0516             // NOTE: this cannot be done earlier since the probability vector
0517             //       can be changed due to size conformance
0518             hyperexp_detail::normalize(param.probs_);
0519
0520             //post: vector size conformance
0521             assert(param.probs_.size() == param.rates_.size());
0522
0523             return is;
0524         }
0525
0526         /** Returns true if the two sets of parameters are the same. */
0527         public: BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
0528         {
0529             return lhs.probs_ == rhs.probs_
0530                    && lhs.rates_ == rhs.rates_;
0531         }
0532
0533         /** Returns true if the two sets of parameters are the different. */
0534         public: BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
0535
0536
0537         private: std::vector<RealT> probs_; ///< The <em>phase probability vector</em> parameter of the distribution
0538         private: std::vector<RealT> rates_; ///< The <em>rate vector</em> parameter of the distribution
0539     }; // param_type
0540
0541
0542     /**
0543      * Constructs a 1-phase \c hyperexponential_distribution (i.e., an
0544      * exponential distribution) with rate 1.
0545      */
0546     public: hyperexponential_distribution()
0547     : dd_(std::vector<RealT>(1, 1)),
0548       rates_(1, 1)
0549     {
0550         // empty
0551     }
0552
0553     /**
0554      * Constructs a \c hyperexponential_distribution from the <em>phase
0555      * probability vector</em> and <em>rate vector</em> parameters of the
0556      * distribution.
0557      *
0558      * The <em>phase probability vector</em> parameter is given by the range
0559      * defined by [\a prob_first, \a prob_last) iterator pair, and the
0560      * <em>rate vector</em> parameter is given by the range defined by
0561      * [\a rate_first, \a rate_last) iterator pair.
0562      *
0563      * \tparam ProbIterT Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]).
0564      * \tparam RateIterT Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]).
0565      *
0566      * \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.
0567      * \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.
0568      * \param rate_first The iterator to the beginning of the range of non-negative real elements representing the rates.
0569      * \param rate_last The iterator to the ending of the range of non-negative real elements representing the rates.
0570      *
0571      * References:
0572      * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014
0573      * .
0574      */
0575     public: template <typename ProbIterT, typename RateIterT>
0576             hyperexponential_distribution(ProbIterT prob_first, ProbIterT prob_last,
0577                                           RateIterT rate_first, RateIterT rate_last)
0578     : dd_(prob_first, prob_last),
0579       rates_(rate_first, rate_last)
0580     {
0581         assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) );
0582     }
0583
0584     /**
0585      * Constructs a \c hyperexponential_distribution from the <em>phase
0586      * probability vector</em> and <em>rate vector</em> parameters of the
0587      * distribution.
0588      *
0589      * The <em>phase probability vector</em> parameter is given by the range
0590      * defined by \a prob_range, and the <em>rate vector</em> parameter is
0591      * given by the range defined by \a rate_range.
0592      *
0593      * \tparam ProbRangeT Must meet the requirements of <a href="boost:/libs/range/doc/html/range/concepts.html">Range</a> concept.
0594      * \tparam RateRangeT Must meet the requirements of <a href="boost:/libs/range/doc/html/range/concepts.html">Range</a> concept.
0595      *
0596      * \param prob_range The range of non-negative real elements representing the phase probabilities; if elements don't sum to 1, they are normalized.
0597      * \param rate_range The range of positive real elements representing the rates.
0598      *
0599      * \note
0600      *  The final \c disable_if parameter is an implementation detail that
0601      *  differentiates between this two argument constructor and the
0602      *  iterator-based two argument constructor described below.
0603      */
0604     //  We SFINAE this out of existance if either argument type is
0605     //  incrementable as in that case the type is probably an iterator:
0606     public: template <typename ProbRangeT, typename RateRangeT>
0607             hyperexponential_distribution(ProbRangeT const& prob_range,
0608                                           RateRangeT const& rate_range,
0609                                           typename boost::disable_if_c<boost::has_pre_increment<ProbRangeT>::value || boost::has_pre_increment<RateRangeT>::value>::type* = 0)
0610     : dd_(prob_range),
0611       rates_(boost::begin(rate_range), boost::end(rate_range))
0612     {
0613         assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) );
0614     }
0615
0616     /**
0617      * Constructs a \c hyperexponential_distribution from the <em>rate
0618      * vector</em> parameter of the distribution and with equal phase
0619      * probabilities.
0620      *
0621      * The <em>rate vector</em> parameter is given by the range defined by
0622      * [\a rate_first, \a rate_last) iterator pair, and the <em>phase
0623      * probability vector</em> parameter is set to the equal phase
0624      * probabilities (i.e., to a vector of the same length \f$k\f$ of the
0625      * <em>rate vector</em> and with each element set to \f$1.0/k\f$).
0626      *
0627      * \tparam RateIterT Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]).
0628      * \tparam RateIterT2 Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]).
0629      *
0630      * \param rate_first The iterator to the beginning of the range of non-negative real elements representing the rates.
0631      * \param rate_last The iterator to the ending of the range of non-negative real elements representing the rates.
0632      *
0633      * \note
0634      *  The final \c disable_if parameter is an implementation detail that
0635      *  differentiates between this two argument constructor and the
0636      *  range-based two argument constructor described above.
0637      *
0638      * References:
0639      * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014
0640      * .
0641      */
0642     //  We SFINAE this out of existance if the argument type is
0643     //  incrementable as in that case the type is probably an iterator.
0644     public: template <typename RateIterT>
0645             hyperexponential_distribution(RateIterT rate_first,
0646                                           RateIterT rate_last,
0647                                           typename boost::enable_if_c<boost::has_pre_increment<RateIterT>::value>::type* = 0)
0648     : dd_(std::vector<RealT>(std::distance(rate_first, rate_last), 1)),
0649       rates_(rate_first, rate_last)
0650     {
0651         assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) );
0652     }
0653
0654     /**
0655      * Constructs a @c param_type from the "rates" parameters
0656      * of the distribution and with equal phase probabilities.
0657      *
0658      * The <em>rate vector</em> parameter is given by the range defined by
0659      * \a rate_range, and the <em>phase probability vector</em> parameter is
0660      * set to the equal phase probabilities (i.e., to a vector of the same
0661      * length \f$k\f$ of the <em>rate vector</em> and with each element set
0662      * to \f$1.0/k\f$).
0663      *
0664      * \tparam RateRangeT Must meet the requirements of <a href="boost:/libs/range/doc/html/range/concepts.html">Range</a> concept.
0665      *
0666      * \param rate_range The range of positive real elements representing the rates.
0667      */
0668     public: template <typename RateRangeT>
0669             hyperexponential_distribution(RateRangeT const& rate_range)
0670     : dd_(std::vector<RealT>(boost::size(rate_range), 1)),
0671       rates_(boost::begin(rate_range), boost::end(rate_range))
0672     {
0673         assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) );
0674     }
0675
0676     /**
0677      * Constructs a \c hyperexponential_distribution from its parameters.
0678      *
0679      * \param param The parameters of the distribution.
0680      */
0681     public: explicit hyperexponential_distribution(param_type const& param)
0682     : dd_(param.probabilities()),
0683       rates_(param.rates())
0684     {
0685         assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) );
0686     }
0687
0688 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
0689     /**
0690      * Constructs a \c hyperexponential_distribution from the <em>phase
0691      * probability vector</em> and <em>rate vector</em> parameters of the
0692      * distribution.
0693      *
0694      * The <em>phase probability vector</em> parameter is given by the
0695      * <em>brace-init-list</em> (ISO,2014,sec. 8.5.4 [dcl.init.list])
0696      * defined by \a l1, and the <em>rate vector</em> parameter is given by the
0697      * <em>brace-init-list</em> (ISO,2014,sec. 8.5.4 [dcl.init.list])
0698      * defined by \a l2.
0699      *
0700      * \param l1 The initializer list for inizializing the phase probability vector.
0701      * \param l2 The initializer list for inizializing the rate vector.
0702      *
0703      * References:
0704      * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014
0705      * .
0706      */
0707     public: hyperexponential_distribution(std::initializer_list<RealT> const& l1, std::initializer_list<RealT> const& l2)
0708     : dd_(l1.begin(), l1.end()),
0709       rates_(l2.begin(), l2.end())
0710     {
0711         assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) );
0712     }
0713
0714     /**
0715      * Constructs a \c hyperexponential_distribution from the <em>rate
0716      * vector</em> parameter of the distribution and with equal phase
0717      * probabilities.
0718      *
0719      * The <em>rate vector</em> parameter is given by the
0720      * <em>brace-init-list</em> (ISO,2014,sec. 8.5.4 [dcl.init.list])
0721      * defined by \a l1, and the <em>phase probability vector</em> parameter is
0722      * set to the equal phase probabilities (i.e., to a vector of the same
0723      * length \f$k\f$ of the <em>rate vector</em> and with each element set
0724      * to \f$1.0/k\f$).
0725      *
0726      * \param l1 The initializer list for inizializing the rate vector.
0727      *
0728      * References:
0729      * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014
0730      * .
0731      */
0732     public: hyperexponential_distribution(std::initializer_list<RealT> const& l1)
0733     : dd_(std::vector<RealT>(std::distance(l1.begin(), l1.end()), 1)),
0734       rates_(l1.begin(), l1.end())
0735     {
0736         assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) );
0737     }
0738 #endif
0739
0740     /**
0741      * Gets a random variate distributed according to the
0742      * hyperexponential distribution.
0743      *
0744      * \tparam URNG Must meet the requirements of \uniform_random_number_generator.
0745      *
0746      * \param urng A uniform random number generator object.
0747      *
0748      * \return A random variate distributed according to the hyperexponential distribution.
0749      */
0750     public: template<class URNG>\
0751             RealT operator()(URNG& urng) const
0752     {
0753         const int i = dd_(urng);
0754
0755         return boost::random::exponential_distribution<RealT>(rates_[i])(urng);
0756     }
0757
0758     /**
0759      * Gets a random variate distributed according to the hyperexponential
0760      * distribution with parameters specified by \c param.
0761      *
0762      * \tparam URNG Must meet the requirements of \uniform_random_number_generator.
0763      *
0764      * \param urng A uniform random number generator object.
0765      * \param param A distribution parameter object.
0766      *
0767      * \return A random variate distributed according to the hyperexponential distribution.
0768      *  distribution with parameters specified by \c param.
0769      */
0770     public: template<class URNG>
0771             RealT operator()(URNG& urng, const param_type& param) const
0772     {
0773         return hyperexponential_distribution(param)(urng);
0774     }
0775
0776     /** Returns the number of phases of the distribution. */
0777     public: std::size_t num_phases() const
0778     {
0779         return rates_.size();
0780     }
0781
0782     /** Returns the <em>phase probability vector</em> parameter of the distribution. */
0783     public: std::vector<RealT> probabilities() const
0784     {
0785         return dd_.probabilities();
0786     }
0787
0788     /** Returns the <em>rate vector</em> parameter of the distribution. */
0789     public: std::vector<RealT> rates() const
0790     {
0791         return rates_;
0792     }
0793
0794     /** Returns the smallest value that the distribution can produce. */
0795     public: RealT min BOOST_PREVENT_MACRO_SUBSTITUTION () const
0796     {
0797         return 0;
0798     }
0799
0800     /** Returns the largest value that the distribution can produce. */
0801     public: RealT max BOOST_PREVENT_MACRO_SUBSTITUTION () const
0802     {
0803         return std::numeric_limits<RealT>::infinity();
0804     }
0805
0806     /** Returns the parameters of the distribution. */
0807     public: param_type param() const
0808     {
0809         std::vector<RealT> probs = dd_.probabilities();
0810
0811         return param_type(probs.begin(), probs.end(), rates_.begin(), rates_.end());
0812     }
0813
0814     /** Sets the parameters of the distribution. */
0815     public: void param(param_type const& param)
0816     {
0817         dd_.param(typename boost::random::discrete_distribution<int,RealT>::param_type(param.probabilities()));
0818         rates_ = param.rates();
0819     }
0820
0821     /**
0822      * Effects: Subsequent uses of the distribution do not depend
0823      * on values produced by any engine prior to invoking reset.
0824      */
0825     public: void reset()
0826     {
0827         // empty
0828     }
0829
0830     /** Writes an @c hyperexponential_distribution to a @c std::ostream. */
0831     public: BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, hyperexponential_distribution, hd)
0832     {
0833         os << hd.param();
0834         return os;
0835     }
0836
0837     /** Reads an @c hyperexponential_distribution from a @c std::istream. */
0838     public: BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, hyperexponential_distribution, hd)
0839     {
0840         param_type param;
0841         if(is >> param)
0842         {
0843             hd.param(param);
0844         }
0845         return is;
0846     }
0847
0848     /**
0849      * Returns true if the two instances of @c hyperexponential_distribution will
0850      * return identical sequences of values given equal generators.
0851      */
0852     public: BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(hyperexponential_distribution, lhs, rhs)
0853     {
0854         return lhs.dd_ == rhs.dd_
0855                && lhs.rates_ == rhs.rates_;
0856     }
0857
0858     /**
0859      * Returns true if the two instances of @c hyperexponential_distribution will
0860      * return different sequences of values given equal generators.
0861      */
0862     public: BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(hyperexponential_distribution)
0863
0864
0865     private: boost::random::discrete_distribution<int,RealT> dd_; ///< The \c discrete_distribution used to sample the phase probability and choose the rate
0866     private: std::vector<RealT> rates_; ///< The <em>rate vector</em> parameter of the distribution
0867 }; // hyperexponential_distribution
0868
0869 }} // namespace boost::random
0870
0871
0872 #endif // BOOST_RANDOM_HYPEREXPONENTIAL_DISTRIBUTION_HPP