Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-30 09:59:13

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