Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-09-17 08:49:57

0001 /* boost random/piecewise_linear_distribution.hpp header file
0002  *
0003  * Copyright Steven Watanabe 2011
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  * $Id$
0011  */
0012 
0013 #ifndef BOOST_RANDOM_PIECEWISE_LINEAR_DISTRIBUTION_HPP_INCLUDED
0014 #define BOOST_RANDOM_PIECEWISE_LINEAR_DISTRIBUTION_HPP_INCLUDED
0015 
0016 #include <vector>
0017 #include <algorithm>
0018 #include <cmath>
0019 #include <cstdlib>
0020 #include <boost/assert.hpp>
0021 #include <boost/random/uniform_real.hpp>
0022 #include <boost/random/discrete_distribution.hpp>
0023 #include <boost/random/detail/config.hpp>
0024 #include <boost/random/detail/operators.hpp>
0025 #include <boost/random/detail/vector_io.hpp>
0026 
0027 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
0028 #include <initializer_list>
0029 #endif
0030 
0031 namespace boost {
0032 namespace random {
0033 
0034 /**
0035  * The class @c piecewise_linear_distribution models a \random_distribution.
0036  */
0037 template<class RealType = double>
0038 class piecewise_linear_distribution {
0039 public:
0040     typedef std::size_t input_type;
0041     typedef RealType result_type;
0042 
0043     class param_type {
0044     public:
0045 
0046         typedef piecewise_linear_distribution distribution_type;
0047 
0048         /**
0049          * Constructs a @c param_type object, representing a distribution
0050          * that produces values uniformly distributed in the range [0, 1).
0051          */
0052         param_type()
0053         {
0054             _weights.push_back(RealType(1));
0055             _weights.push_back(RealType(1));
0056             _intervals.push_back(RealType(0));
0057             _intervals.push_back(RealType(1));
0058         }
0059         /**
0060          * Constructs a @c param_type object from two iterator ranges
0061          * containing the interval boundaries and weights at the boundaries.
0062          * If there are fewer than two boundaries, then this is equivalent to
0063          * the default constructor and the distribution will produce values
0064          * uniformly distributed in the range [0, 1).
0065          *
0066          * The values of the interval boundaries must be strictly
0067          * increasing, and the number of weights must be the same as
0068          * the number of interval boundaries.  If there are extra
0069          * weights, they are ignored.
0070          */
0071         template<class IntervalIter, class WeightIter>
0072         param_type(IntervalIter intervals_first, IntervalIter intervals_last,
0073                    WeightIter weight_first)
0074           : _intervals(intervals_first, intervals_last)
0075         {
0076             if(_intervals.size() < 2) {
0077                 _intervals.clear();
0078                 _weights.push_back(RealType(1));
0079                 _weights.push_back(RealType(1));
0080                 _intervals.push_back(RealType(0));
0081                 _intervals.push_back(RealType(1));
0082             } else {
0083                 _weights.reserve(_intervals.size());
0084                 for(std::size_t i = 0; i < _intervals.size(); ++i) {
0085                     _weights.push_back(*weight_first++);
0086                 }
0087             }
0088         }
0089 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
0090         /**
0091          * Constructs a @c param_type object from an initializer_list
0092          * containing the interval boundaries and a unary function
0093          * specifying the weights at the boundaries.  Each weight is
0094          * determined by calling the function at the corresponding point.
0095          *
0096          * If the initializer_list contains fewer than two elements,
0097          * this is equivalent to the default constructor and the
0098          * distribution will produce values uniformly distributed
0099          * in the range [0, 1).
0100          */
0101         template<class T, class F>
0102         param_type(const std::initializer_list<T>& il, F f)
0103           : _intervals(il.begin(), il.end())
0104         {
0105             if(_intervals.size() < 2) {
0106                 _intervals.clear();
0107                 _weights.push_back(RealType(1));
0108                 _weights.push_back(RealType(1));
0109                 _intervals.push_back(RealType(0));
0110                 _intervals.push_back(RealType(1));
0111             } else {
0112                 _weights.reserve(_intervals.size());
0113                 for(typename std::vector<RealType>::const_iterator
0114                     iter = _intervals.begin(), end = _intervals.end();
0115                     iter != end; ++iter)
0116                 {
0117                     _weights.push_back(f(*iter));
0118                 }
0119             }
0120         }
0121 #endif
0122         /**
0123          * Constructs a @c param_type object from Boost.Range ranges holding
0124          * the interval boundaries and the weights at the boundaries.  If
0125          * there are fewer than two interval boundaries, this is equivalent
0126          * to the default constructor and the distribution will produce
0127          * values uniformly distributed in the range [0, 1).  The
0128          * number of weights must be equal to the number of
0129          * interval boundaries.
0130          */
0131         template<class IntervalRange, class WeightRange>
0132         param_type(const IntervalRange& intervals_arg,
0133                    const WeightRange& weights_arg)
0134           : _intervals(std::begin(intervals_arg), std::end(intervals_arg)),
0135             _weights(std::begin(weights_arg), std::end(weights_arg))
0136         {
0137             if(_intervals.size() < 2) {
0138                 _weights.clear();
0139                 _weights.push_back(RealType(1));
0140                 _weights.push_back(RealType(1));
0141                 _intervals.clear();
0142                 _intervals.push_back(RealType(0));
0143                 _intervals.push_back(RealType(1));
0144             }
0145         }
0146 
0147         /**
0148          * Constructs the parameters for a distribution that approximates a
0149          * function.  The range of the distribution is [xmin, xmax).  This
0150          * range is divided into nw equally sized intervals and the weights
0151          * are found by calling the unary function f on the boundaries of the
0152          * intervals.
0153          */
0154         template<class F>
0155         param_type(std::size_t nw, RealType xmin, RealType xmax, F f)
0156         {
0157             std::size_t n = (nw == 0) ? 1 : nw;
0158             double delta = (xmax - xmin) / n;
0159             BOOST_ASSERT(delta > 0);
0160             for(std::size_t k = 0; k < n; ++k) {
0161                 _weights.push_back(f(xmin + k*delta));
0162                 _intervals.push_back(xmin + k*delta);
0163             }
0164             _weights.push_back(f(xmax));
0165             _intervals.push_back(xmax);
0166         }
0167 
0168         /**  Returns a vector containing the interval boundaries. */
0169         std::vector<RealType> intervals() const { return _intervals; }
0170 
0171         /**
0172          * Returns a vector containing the probability densities
0173          * at all the interval boundaries.
0174          */
0175         std::vector<RealType> densities() const
0176         {
0177             RealType sum = static_cast<RealType>(0);
0178             for(std::size_t i = 0; i < _intervals.size() - 1; ++i) {
0179                 RealType width = _intervals[i + 1] - _intervals[i];
0180                 sum += (_weights[i] + _weights[i + 1]) * width / 2;
0181             }
0182             std::vector<RealType> result;
0183             result.reserve(_weights.size());
0184             for(typename std::vector<RealType>::const_iterator
0185                 iter = _weights.begin(), end = _weights.end();
0186                 iter != end; ++iter)
0187             {
0188                 result.push_back(*iter / sum);
0189             }
0190             return result;
0191         }
0192 
0193         /** Writes the parameters to a @c std::ostream. */
0194         BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
0195         {
0196             detail::print_vector(os, parm._intervals);
0197             detail::print_vector(os, parm._weights);
0198             return os;
0199         }
0200 
0201         /** Reads the parameters from a @c std::istream. */
0202         BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
0203         {
0204             std::vector<RealType> new_intervals;
0205             std::vector<RealType> new_weights;
0206             detail::read_vector(is, new_intervals);
0207             detail::read_vector(is, new_weights);
0208             if(is) {
0209                 parm._intervals.swap(new_intervals);
0210                 parm._weights.swap(new_weights);
0211             }
0212             return is;
0213         }
0214 
0215         /** Returns true if the two sets of parameters are the same. */
0216         BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
0217         {
0218             return lhs._intervals == rhs._intervals
0219                 && lhs._weights == rhs._weights;
0220         }
0221         /** Returns true if the two sets of parameters are different. */
0222         BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
0223 
0224     private:
0225         friend class piecewise_linear_distribution;
0226 
0227         std::vector<RealType> _intervals;
0228         std::vector<RealType> _weights;
0229     };
0230 
0231     /**
0232      * Creates a new @c piecewise_linear_distribution that
0233      * produces values uniformly distributed in the range [0, 1).
0234      */
0235     piecewise_linear_distribution()
0236     {
0237         default_init();
0238     }
0239     /**
0240      * Constructs a piecewise_linear_distribution from two iterator ranges
0241      * containing the interval boundaries and the weights at the boundaries.
0242      * If there are fewer than two boundaries, then this is equivalent to
0243      * the default constructor and creates a distribution that
0244      * produces values uniformly distributed in the range [0, 1).
0245      *
0246      * The values of the interval boundaries must be strictly
0247      * increasing, and the number of weights must be equal to
0248      * the number of interval boundaries.  If there are extra
0249      * weights, they are ignored.
0250      *
0251      * For example,
0252      *
0253      * @code
0254      * double intervals[] = { 0.0, 1.0, 2.0 };
0255      * double weights[] = { 0.0, 1.0, 0.0 };
0256      * piecewise_constant_distribution<> dist(
0257      *     &intervals[0], &intervals[0] + 3, &weights[0]);
0258      * @endcode
0259      *
0260      * produces a triangle distribution.
0261      */
0262     template<class IntervalIter, class WeightIter>
0263     piecewise_linear_distribution(IntervalIter first_interval,
0264                                   IntervalIter last_interval,
0265                                   WeightIter first_weight)
0266       : _intervals(first_interval, last_interval)
0267     {
0268         if(_intervals.size() < 2) {
0269             default_init();
0270         } else {
0271             _weights.reserve(_intervals.size());
0272             for(std::size_t i = 0; i < _intervals.size(); ++i) {
0273                 _weights.push_back(*first_weight++);
0274             }
0275             init();
0276         }
0277     }
0278 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
0279     /**
0280      * Constructs a piecewise_linear_distribution from an
0281      * initializer_list containing the interval boundaries
0282      * and a unary function specifying the weights.  Each
0283      * weight is determined by calling the function at the
0284      * corresponding interval boundary.
0285      *
0286      * If the initializer_list contains fewer than two elements,
0287      * this is equivalent to the default constructor and the
0288      * distribution will produce values uniformly distributed
0289      * in the range [0, 1).
0290      */
0291     template<class T, class F>
0292     piecewise_linear_distribution(std::initializer_list<T> il, F f)
0293       : _intervals(il.begin(), il.end())
0294     {
0295         if(_intervals.size() < 2) {
0296             default_init();
0297         } else {
0298             _weights.reserve(_intervals.size());
0299             for(typename std::vector<RealType>::const_iterator
0300                 iter = _intervals.begin(), end = _intervals.end();
0301                 iter != end; ++iter)
0302             {
0303                 _weights.push_back(f(*iter));
0304             }
0305             init();
0306         }
0307     }
0308 #endif
0309     /**
0310      * Constructs a piecewise_linear_distribution from Boost.Range
0311      * ranges holding the interval boundaries and the weights.  If
0312      * there are fewer than two interval boundaries, this is equivalent
0313      * to the default constructor and the distribution will produce
0314      * values uniformly distributed in the range [0, 1).  The
0315      * number of weights must be equal to the number of
0316      * interval boundaries.
0317      */
0318     template<class IntervalsRange, class WeightsRange>
0319     piecewise_linear_distribution(const IntervalsRange& intervals_arg,
0320                                   const WeightsRange& weights_arg)
0321       : _intervals(std::begin(intervals_arg), std::end(intervals_arg)),
0322         _weights(std::begin(weights_arg), std::end(weights_arg))
0323     {
0324         if(_intervals.size() < 2) {
0325             default_init();
0326         } else {
0327             init();
0328         }
0329     }
0330     /**
0331      * Constructs a piecewise_linear_distribution that approximates a
0332      * function.  The range of the distribution is [xmin, xmax).  This
0333      * range is divided into nw equally sized intervals and the weights
0334      * are found by calling the unary function f on the interval boundaries.
0335      */
0336     template<class F>
0337     piecewise_linear_distribution(std::size_t nw,
0338                                   RealType xmin,
0339                                   RealType xmax,
0340                                   F f)
0341     {
0342         if(nw == 0) { nw = 1; }
0343         RealType delta = (xmax - xmin) / nw;
0344         _intervals.reserve(nw + 1);
0345         for(std::size_t i = 0; i < nw; ++i) {
0346             RealType x = xmin + i * delta;
0347             _intervals.push_back(x);
0348             _weights.push_back(f(x));
0349         }
0350         _intervals.push_back(xmax);
0351         _weights.push_back(f(xmax));
0352         init();
0353     }
0354     /**
0355      * Constructs a piecewise_linear_distribution from its parameters.
0356      */
0357     explicit piecewise_linear_distribution(const param_type& parm)
0358       : _intervals(parm._intervals),
0359         _weights(parm._weights)
0360     {
0361         init();
0362     }
0363 
0364     /**
0365      * Returns a value distributed according to the parameters of the
0366      * piecewise_linear_distribution.
0367      */
0368     template<class URNG>
0369     RealType operator()(URNG& urng) const
0370     {
0371         std::size_t i = _bins(urng);
0372         bool is_in_rectangle = (i % 2 == 0);
0373         i = i / 2;
0374         uniform_real<RealType> dist(_intervals[i], _intervals[i+1]);
0375         if(is_in_rectangle) {
0376             return dist(urng);
0377         } else if(_weights[i] < _weights[i+1]) {
0378             return (std::max)(dist(urng), dist(urng));
0379         } else {
0380             return (std::min)(dist(urng), dist(urng));
0381         }
0382     }
0383 
0384     /**
0385      * Returns a value distributed according to the parameters
0386      * specified by param.
0387      */
0388     template<class URNG>
0389     RealType operator()(URNG& urng, const param_type& parm) const
0390     {
0391         return piecewise_linear_distribution(parm)(urng);
0392     }
0393 
0394     /** Returns the smallest value that the distribution can produce. */
0395     result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const
0396     { return _intervals.front(); }
0397     /** Returns the largest value that the distribution can produce. */
0398     result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const
0399     { return _intervals.back(); }
0400 
0401     /**
0402      * Returns a vector containing the probability densities
0403      * at the interval boundaries.
0404      */
0405     std::vector<RealType> densities() const
0406     {
0407         RealType sum = static_cast<RealType>(0);
0408         for(std::size_t i = 0; i < _intervals.size() - 1; ++i) {
0409             RealType width = _intervals[i + 1] - _intervals[i];
0410             sum += (_weights[i] + _weights[i + 1]) * width / 2;
0411         }
0412         std::vector<RealType> result;
0413         result.reserve(_weights.size());
0414         for(typename std::vector<RealType>::const_iterator
0415             iter = _weights.begin(), end = _weights.end();
0416             iter != end; ++iter)
0417         {
0418             result.push_back(*iter / sum);
0419         }
0420         return result;
0421     }
0422     /**  Returns a vector containing the interval boundaries. */
0423     std::vector<RealType> intervals() const { return _intervals; }
0424 
0425     /** Returns the parameters of the distribution. */
0426     param_type param() const
0427     {
0428         return param_type(_intervals, _weights);
0429     }
0430     /** Sets the parameters of the distribution. */
0431     void param(const param_type& parm)
0432     {
0433         std::vector<RealType> new_intervals(parm._intervals);
0434         std::vector<RealType> new_weights(parm._weights);
0435         init(new_intervals, new_weights);
0436         _intervals.swap(new_intervals);
0437         _weights.swap(new_weights);
0438     }
0439 
0440     /**
0441      * Effects: Subsequent uses of the distribution do not depend
0442      * on values produced by any engine prior to invoking reset.
0443      */
0444     void reset() { _bins.reset(); }
0445 
0446     /** Writes a distribution to a @c std::ostream. */
0447     BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(
0448         os, piecewise_linear_distribution, pld)
0449     {
0450         os << pld.param();
0451         return os;
0452     }
0453 
0454     /** Reads a distribution from a @c std::istream */
0455     BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(
0456         is, piecewise_linear_distribution, pld)
0457     {
0458         param_type parm;
0459         if(is >> parm) {
0460             pld.param(parm);
0461         }
0462         return is;
0463     }
0464 
0465     /**
0466      * Returns true if the two distributions will return the
0467      * same sequence of values, when passed equal generators.
0468      */
0469     BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(
0470         piecewise_linear_distribution, lhs,  rhs)
0471     {
0472         return lhs._intervals == rhs._intervals && lhs._weights == rhs._weights;
0473     }
0474     /**
0475      * Returns true if the two distributions may return different
0476      * sequences of values, when passed equal generators.
0477      */
0478     BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(piecewise_linear_distribution)
0479 
0480 private:
0481 
0482     /// @cond \show_private
0483 
0484     void init(const std::vector<RealType>& intervals_arg,
0485               const std::vector<RealType>& weights_arg)
0486     {
0487         using std::abs;
0488         std::vector<RealType> bin_weights;
0489         bin_weights.reserve((intervals_arg.size() - 1) * 2);
0490         for(std::size_t i = 0; i < intervals_arg.size() - 1; ++i) {
0491             RealType width = intervals_arg[i + 1] - intervals_arg[i];
0492             RealType w1 = weights_arg[i];
0493             RealType w2 = weights_arg[i + 1];
0494             bin_weights.push_back((std::min)(w1, w2) * width);
0495             bin_weights.push_back(abs(w1 - w2) * width / 2);
0496         }
0497         typedef discrete_distribution<std::size_t, RealType> bins_type;
0498         typename bins_type::param_type bins_param(bin_weights);
0499         _bins.param(bins_param);
0500     }
0501 
0502     void init()
0503     {
0504         init(_intervals, _weights);
0505     }
0506 
0507     void default_init()
0508     {
0509         _intervals.clear();
0510         _intervals.push_back(RealType(0));
0511         _intervals.push_back(RealType(1));
0512         _weights.clear();
0513         _weights.push_back(RealType(1));
0514         _weights.push_back(RealType(1));
0515         init();
0516     }
0517 
0518     discrete_distribution<std::size_t, RealType> _bins;
0519     std::vector<RealType> _intervals;
0520     std::vector<RealType> _weights;
0521 
0522     /// @endcond
0523 };
0524 
0525 }
0526 }
0527 
0528 #endif