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