Back to home page

EIC code displayed by LXR

 
 

    


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

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