Chris@16: /* boost random/piecewise_linear_distribution.hpp header file
Chris@16:  *
Chris@16:  * Copyright Steven Watanabe 2011
Chris@16:  * Distributed under the Boost Software License, Version 1.0. (See
Chris@16:  * accompanying file LICENSE_1_0.txt or copy at
Chris@16:  * http://www.boost.org/LICENSE_1_0.txt)
Chris@16:  *
Chris@16:  * See http://www.boost.org for most recent version including documentation.
Chris@16:  *
Chris@101:  * $Id$
Chris@16:  */
Chris@16: 
Chris@16: #ifndef BOOST_RANDOM_PIECEWISE_LINEAR_DISTRIBUTION_HPP_INCLUDED
Chris@16: #define BOOST_RANDOM_PIECEWISE_LINEAR_DISTRIBUTION_HPP_INCLUDED
Chris@16: 
Chris@16: #include <vector>
Chris@16: #include <algorithm>
Chris@16: #include <cmath>
Chris@16: #include <cstdlib>
Chris@16: #include <boost/assert.hpp>
Chris@16: #include <boost/random/uniform_real.hpp>
Chris@16: #include <boost/random/discrete_distribution.hpp>
Chris@16: #include <boost/random/detail/config.hpp>
Chris@16: #include <boost/random/detail/operators.hpp>
Chris@16: #include <boost/random/detail/vector_io.hpp>
Chris@16: 
Chris@16: #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
Chris@16: #include <initializer_list>
Chris@16: #endif
Chris@16: 
Chris@16: #include <boost/range/begin.hpp>
Chris@16: #include <boost/range/end.hpp>
Chris@16: 
Chris@16: namespace boost {
Chris@16: namespace random {
Chris@16: 
Chris@16: /**
Chris@16:  * The class @c piecewise_linear_distribution models a \random_distribution.
Chris@16:  */
Chris@16: template<class RealType = double>
Chris@16: class piecewise_linear_distribution {
Chris@16: public:
Chris@16:     typedef std::size_t input_type;
Chris@16:     typedef RealType result_type;
Chris@16: 
Chris@16:     class param_type {
Chris@16:     public:
Chris@16: 
Chris@16:         typedef piecewise_linear_distribution distribution_type;
Chris@16: 
Chris@16:         /**
Chris@16:          * Constructs a @c param_type object, representing a distribution
Chris@16:          * that produces values uniformly distributed in the range [0, 1).
Chris@16:          */
Chris@16:         param_type()
Chris@16:         {
Chris@16:             _weights.push_back(RealType(1));
Chris@16:             _weights.push_back(RealType(1));
Chris@16:             _intervals.push_back(RealType(0));
Chris@16:             _intervals.push_back(RealType(1));
Chris@16:         }
Chris@16:         /**
Chris@16:          * Constructs a @c param_type object from two iterator ranges
Chris@16:          * containing the interval boundaries and weights at the boundaries.
Chris@16:          * If there are fewer than two boundaries, then this is equivalent to
Chris@16:          * the default constructor and the distribution will produce values
Chris@16:          * uniformly distributed in the range [0, 1).
Chris@16:          *
Chris@16:          * The values of the interval boundaries must be strictly
Chris@16:          * increasing, and the number of weights must be the same as
Chris@16:          * the number of interval boundaries.  If there are extra
Chris@16:          * weights, they are ignored.
Chris@16:          */
Chris@16:         template<class IntervalIter, class WeightIter>
Chris@16:         param_type(IntervalIter intervals_first, IntervalIter intervals_last,
Chris@16:                    WeightIter weight_first)
Chris@16:           : _intervals(intervals_first, intervals_last)
Chris@16:         {
Chris@16:             if(_intervals.size() < 2) {
Chris@16:                 _intervals.clear();
Chris@16:                 _weights.push_back(RealType(1));
Chris@16:                 _weights.push_back(RealType(1));
Chris@16:                 _intervals.push_back(RealType(0));
Chris@16:                 _intervals.push_back(RealType(1));
Chris@16:             } else {
Chris@16:                 _weights.reserve(_intervals.size());
Chris@16:                 for(std::size_t i = 0; i < _intervals.size(); ++i) {
Chris@16:                     _weights.push_back(*weight_first++);
Chris@16:                 }
Chris@16:             }
Chris@16:         }
Chris@16: #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
Chris@16:         /**
Chris@16:          * Constructs a @c param_type object from an initializer_list
Chris@16:          * containing the interval boundaries and a unary function
Chris@16:          * specifying the weights at the boundaries.  Each weight is
Chris@16:          * determined by calling the function at the corresponding point.
Chris@16:          *
Chris@16:          * If the initializer_list contains fewer than two elements,
Chris@16:          * this is equivalent to the default constructor and the
Chris@16:          * distribution will produce values uniformly distributed
Chris@16:          * in the range [0, 1).
Chris@16:          */
Chris@16:         template<class T, class F>
Chris@16:         param_type(const std::initializer_list<T>& il, F f)
Chris@16:           : _intervals(il.begin(), il.end())
Chris@16:         {
Chris@16:             if(_intervals.size() < 2) {
Chris@16:                 _intervals.clear();
Chris@16:                 _weights.push_back(RealType(1));
Chris@16:                 _weights.push_back(RealType(1));
Chris@16:                 _intervals.push_back(RealType(0));
Chris@16:                 _intervals.push_back(RealType(1));
Chris@16:             } else {
Chris@16:                 _weights.reserve(_intervals.size());
Chris@16:                 for(typename std::vector<RealType>::const_iterator
Chris@16:                     iter = _intervals.begin(), end = _intervals.end();
Chris@16:                     iter != end; ++iter)
Chris@16:                 {
Chris@16:                     _weights.push_back(f(*iter));
Chris@16:                 }
Chris@16:             }
Chris@16:         }
Chris@16: #endif
Chris@16:         /**
Chris@16:          * Constructs a @c param_type object from Boost.Range ranges holding
Chris@16:          * the interval boundaries and the weights at the boundaries.  If
Chris@16:          * there are fewer than two interval boundaries, this is equivalent
Chris@16:          * to the default constructor and the distribution will produce
Chris@16:          * values uniformly distributed in the range [0, 1).  The
Chris@16:          * number of weights must be equal to the number of
Chris@16:          * interval boundaries.
Chris@16:          */
Chris@16:         template<class IntervalRange, class WeightRange>
Chris@16:         param_type(const IntervalRange& intervals_arg,
Chris@16:                    const WeightRange& weights_arg)
Chris@16:           : _intervals(boost::begin(intervals_arg), boost::end(intervals_arg)),
Chris@16:             _weights(boost::begin(weights_arg), boost::end(weights_arg))
Chris@16:         {
Chris@16:             if(_intervals.size() < 2) {
Chris@16:                 _weights.clear();
Chris@16:                 _weights.push_back(RealType(1));
Chris@16:                 _weights.push_back(RealType(1));
Chris@16:                 _intervals.clear();
Chris@16:                 _intervals.push_back(RealType(0));
Chris@16:                 _intervals.push_back(RealType(1));
Chris@16:             }
Chris@16:         }
Chris@16: 
Chris@16:         /**
Chris@16:          * Constructs the parameters for a distribution that approximates a
Chris@16:          * function.  The range of the distribution is [xmin, xmax).  This
Chris@16:          * range is divided into nw equally sized intervals and the weights
Chris@16:          * are found by calling the unary function f on the boundaries of the
Chris@16:          * intervals.
Chris@16:          */
Chris@16:         template<class F>
Chris@16:         param_type(std::size_t nw, RealType xmin, RealType xmax, F f)
Chris@16:         {
Chris@16:             std::size_t n = (nw == 0) ? 1 : nw;
Chris@16:             double delta = (xmax - xmin) / n;
Chris@16:             BOOST_ASSERT(delta > 0);
Chris@16:             for(std::size_t k = 0; k < n; ++k) {
Chris@16:                 _weights.push_back(f(xmin + k*delta));
Chris@16:                 _intervals.push_back(xmin + k*delta);
Chris@16:             }
Chris@16:             _weights.push_back(f(xmax));
Chris@16:             _intervals.push_back(xmax);
Chris@16:         }
Chris@16: 
Chris@16:         /**  Returns a vector containing the interval boundaries. */
Chris@16:         std::vector<RealType> intervals() const { return _intervals; }
Chris@16: 
Chris@16:         /**
Chris@16:          * Returns a vector containing the probability densities
Chris@16:          * at all the interval boundaries.
Chris@16:          */
Chris@16:         std::vector<RealType> densities() const
Chris@16:         {
Chris@16:             RealType sum = static_cast<RealType>(0);
Chris@16:             for(std::size_t i = 0; i < _intervals.size() - 1; ++i) {
Chris@16:                 RealType width = _intervals[i + 1] - _intervals[i];
Chris@16:                 sum += (_weights[i] + _weights[i + 1]) * width / 2;
Chris@16:             }
Chris@16:             std::vector<RealType> result;
Chris@16:             result.reserve(_weights.size());
Chris@16:             for(typename std::vector<RealType>::const_iterator
Chris@16:                 iter = _weights.begin(), end = _weights.end();
Chris@16:                 iter != end; ++iter)
Chris@16:             {
Chris@16:                 result.push_back(*iter / sum);
Chris@16:             }
Chris@16:             return result;
Chris@16:         }
Chris@16: 
Chris@16:         /** Writes the parameters to a @c std::ostream. */
Chris@16:         BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
Chris@16:         {
Chris@16:             detail::print_vector(os, parm._intervals);
Chris@16:             detail::print_vector(os, parm._weights);
Chris@16:             return os;
Chris@16:         }
Chris@16:         
Chris@16:         /** Reads the parameters from a @c std::istream. */
Chris@16:         BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
Chris@16:         {
Chris@16:             std::vector<RealType> new_intervals;
Chris@16:             std::vector<RealType> new_weights;
Chris@16:             detail::read_vector(is, new_intervals);
Chris@16:             detail::read_vector(is, new_weights);
Chris@16:             if(is) {
Chris@16:                 parm._intervals.swap(new_intervals);
Chris@16:                 parm._weights.swap(new_weights);
Chris@16:             }
Chris@16:             return is;
Chris@16:         }
Chris@16: 
Chris@16:         /** Returns true if the two sets of parameters are the same. */
Chris@16:         BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
Chris@16:         {
Chris@16:             return lhs._intervals == rhs._intervals
Chris@16:                 && lhs._weights == rhs._weights;
Chris@16:         }
Chris@16:         /** Returns true if the two sets of parameters are different. */
Chris@16:         BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
Chris@16: 
Chris@16:     private:
Chris@16:         friend class piecewise_linear_distribution;
Chris@16: 
Chris@16:         std::vector<RealType> _intervals;
Chris@16:         std::vector<RealType> _weights;
Chris@16:     };
Chris@16: 
Chris@16:     /**
Chris@16:      * Creates a new @c piecewise_linear_distribution that
Chris@16:      * produces values uniformly distributed in the range [0, 1).
Chris@16:      */
Chris@16:     piecewise_linear_distribution()
Chris@16:     {
Chris@16:         default_init();
Chris@16:     }
Chris@16:     /**
Chris@16:      * Constructs a piecewise_linear_distribution from two iterator ranges
Chris@16:      * containing the interval boundaries and the weights at the boundaries.
Chris@16:      * If there are fewer than two boundaries, then this is equivalent to
Chris@16:      * the default constructor and creates a distribution that
Chris@16:      * produces values uniformly distributed in the range [0, 1).
Chris@16:      *
Chris@16:      * The values of the interval boundaries must be strictly
Chris@16:      * increasing, and the number of weights must be equal to
Chris@16:      * the number of interval boundaries.  If there are extra
Chris@16:      * weights, they are ignored.
Chris@16:      *
Chris@16:      * For example,
Chris@16:      *
Chris@16:      * @code
Chris@16:      * double intervals[] = { 0.0, 1.0, 2.0 };
Chris@16:      * double weights[] = { 0.0, 1.0, 0.0 };
Chris@16:      * piecewise_constant_distribution<> dist(
Chris@16:      *     &intervals[0], &intervals[0] + 3, &weights[0]);
Chris@16:      * @endcode
Chris@16:      *
Chris@16:      * produces a triangle distribution.
Chris@16:      */
Chris@16:     template<class IntervalIter, class WeightIter>
Chris@16:     piecewise_linear_distribution(IntervalIter first_interval,
Chris@16:                                   IntervalIter last_interval,
Chris@16:                                   WeightIter first_weight)
Chris@16:       : _intervals(first_interval, last_interval)
Chris@16:     {
Chris@16:         if(_intervals.size() < 2) {
Chris@16:             default_init();
Chris@16:         } else {
Chris@16:             _weights.reserve(_intervals.size());
Chris@16:             for(std::size_t i = 0; i < _intervals.size(); ++i) {
Chris@16:                 _weights.push_back(*first_weight++);
Chris@16:             }
Chris@16:             init();
Chris@16:         }
Chris@16:     }
Chris@16: #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
Chris@16:     /**
Chris@16:      * Constructs a piecewise_linear_distribution from an
Chris@16:      * initializer_list containing the interval boundaries
Chris@16:      * and a unary function specifying the weights.  Each
Chris@16:      * weight is determined by calling the function at the
Chris@16:      * corresponding interval boundary.
Chris@16:      *
Chris@16:      * If the initializer_list contains fewer than two elements,
Chris@16:      * this is equivalent to the default constructor and the
Chris@16:      * distribution will produce values uniformly distributed
Chris@16:      * in the range [0, 1).
Chris@16:      */
Chris@16:     template<class T, class F>
Chris@16:     piecewise_linear_distribution(std::initializer_list<T> il, F f)
Chris@16:       : _intervals(il.begin(), il.end())
Chris@16:     {
Chris@16:         if(_intervals.size() < 2) {
Chris@16:             default_init();
Chris@16:         } else {
Chris@16:             _weights.reserve(_intervals.size());
Chris@16:             for(typename std::vector<RealType>::const_iterator
Chris@16:                 iter = _intervals.begin(), end = _intervals.end();
Chris@16:                 iter != end; ++iter)
Chris@16:             {
Chris@16:                 _weights.push_back(f(*iter));
Chris@16:             }
Chris@16:             init();
Chris@16:         }
Chris@16:     }
Chris@16: #endif
Chris@16:     /**
Chris@16:      * Constructs a piecewise_linear_distribution from Boost.Range
Chris@16:      * ranges holding the interval boundaries and the weights.  If
Chris@16:      * there are fewer than two interval boundaries, this is equivalent
Chris@16:      * to the default constructor and the distribution will produce
Chris@16:      * values uniformly distributed in the range [0, 1).  The
Chris@16:      * number of weights must be equal to the number of
Chris@16:      * interval boundaries.
Chris@16:      */
Chris@16:     template<class IntervalsRange, class WeightsRange>
Chris@16:     piecewise_linear_distribution(const IntervalsRange& intervals_arg,
Chris@16:                                   const WeightsRange& weights_arg)
Chris@16:       : _intervals(boost::begin(intervals_arg), boost::end(intervals_arg)),
Chris@16:         _weights(boost::begin(weights_arg), boost::end(weights_arg))
Chris@16:     {
Chris@16:         if(_intervals.size() < 2) {
Chris@16:             default_init();
Chris@16:         } else {
Chris@16:             init();
Chris@16:         }
Chris@16:     }
Chris@16:     /**
Chris@16:      * Constructs a piecewise_linear_distribution that approximates a
Chris@16:      * function.  The range of the distribution is [xmin, xmax).  This
Chris@16:      * range is divided into nw equally sized intervals and the weights
Chris@16:      * are found by calling the unary function f on the interval boundaries.
Chris@16:      */
Chris@16:     template<class F>
Chris@16:     piecewise_linear_distribution(std::size_t nw,
Chris@16:                                   RealType xmin,
Chris@16:                                   RealType xmax,
Chris@16:                                   F f)
Chris@16:     {
Chris@16:         if(nw == 0) { nw = 1; }
Chris@16:         RealType delta = (xmax - xmin) / nw;
Chris@16:         _intervals.reserve(nw + 1);
Chris@16:         for(std::size_t i = 0; i < nw; ++i) {
Chris@16:             RealType x = xmin + i * delta;
Chris@16:             _intervals.push_back(x);
Chris@16:             _weights.push_back(f(x));
Chris@16:         }
Chris@16:         _intervals.push_back(xmax);
Chris@16:         _weights.push_back(f(xmax));
Chris@16:         init();
Chris@16:     }
Chris@16:     /**
Chris@16:      * Constructs a piecewise_linear_distribution from its parameters.
Chris@16:      */
Chris@16:     explicit piecewise_linear_distribution(const param_type& parm)
Chris@16:       : _intervals(parm._intervals),
Chris@16:         _weights(parm._weights)
Chris@16:     {
Chris@16:         init();
Chris@16:     }
Chris@16: 
Chris@16:     /**
Chris@16:      * Returns a value distributed according to the parameters of the
Chris@16:      * piecewise_linear_distribution.
Chris@16:      */
Chris@16:     template<class URNG>
Chris@16:     RealType operator()(URNG& urng) const
Chris@16:     {
Chris@16:         std::size_t i = _bins(urng);
Chris@16:         bool is_in_rectangle = (i % 2 == 0);
Chris@16:         i = i / 2;
Chris@16:         uniform_real<RealType> dist(_intervals[i], _intervals[i+1]);
Chris@16:         if(is_in_rectangle) {
Chris@16:             return dist(urng);
Chris@16:         } else if(_weights[i] < _weights[i+1]) {
Chris@16:             return (std::max)(dist(urng), dist(urng));
Chris@16:         } else {
Chris@16:             return (std::min)(dist(urng), dist(urng));
Chris@16:         }
Chris@16:     }
Chris@16:     
Chris@16:     /**
Chris@16:      * Returns a value distributed according to the parameters
Chris@16:      * specified by param.
Chris@16:      */
Chris@16:     template<class URNG>
Chris@16:     RealType operator()(URNG& urng, const param_type& parm) const
Chris@16:     {
Chris@16:         return piecewise_linear_distribution(parm)(urng);
Chris@16:     }
Chris@16:     
Chris@16:     /** Returns the smallest value that the distribution can produce. */
Chris@16:     result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const
Chris@16:     { return _intervals.front(); }
Chris@16:     /** Returns the largest value that the distribution can produce. */
Chris@16:     result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const
Chris@16:     { return _intervals.back(); }
Chris@16: 
Chris@16:     /**
Chris@16:      * Returns a vector containing the probability densities
Chris@16:      * at the interval boundaries.
Chris@16:      */
Chris@16:     std::vector<RealType> densities() const
Chris@16:     {
Chris@16:         RealType sum = static_cast<RealType>(0);
Chris@16:         for(std::size_t i = 0; i < _intervals.size() - 1; ++i) {
Chris@16:             RealType width = _intervals[i + 1] - _intervals[i];
Chris@16:             sum += (_weights[i] + _weights[i + 1]) * width / 2;
Chris@16:         }
Chris@16:         std::vector<RealType> result;
Chris@16:         result.reserve(_weights.size());
Chris@16:         for(typename std::vector<RealType>::const_iterator
Chris@16:             iter = _weights.begin(), end = _weights.end();
Chris@16:             iter != end; ++iter)
Chris@16:         {
Chris@16:             result.push_back(*iter / sum);
Chris@16:         }
Chris@16:         return result;
Chris@16:     }
Chris@16:     /**  Returns a vector containing the interval boundaries. */
Chris@16:     std::vector<RealType> intervals() const { return _intervals; }
Chris@16: 
Chris@16:     /** Returns the parameters of the distribution. */
Chris@16:     param_type param() const
Chris@16:     {
Chris@16:         return param_type(_intervals, _weights);
Chris@16:     }
Chris@16:     /** Sets the parameters of the distribution. */
Chris@16:     void param(const param_type& parm)
Chris@16:     {
Chris@16:         std::vector<RealType> new_intervals(parm._intervals);
Chris@16:         std::vector<RealType> new_weights(parm._weights);
Chris@16:         init(new_intervals, new_weights);
Chris@16:         _intervals.swap(new_intervals);
Chris@16:         _weights.swap(new_weights);
Chris@16:     }
Chris@16:     
Chris@16:     /**
Chris@16:      * Effects: Subsequent uses of the distribution do not depend
Chris@16:      * on values produced by any engine prior to invoking reset.
Chris@16:      */
Chris@16:     void reset() { _bins.reset(); }
Chris@16: 
Chris@16:     /** Writes a distribution to a @c std::ostream. */
Chris@16:     BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(
Chris@16:         os, piecewise_linear_distribution, pld)
Chris@16:     {
Chris@16:         os << pld.param();
Chris@16:         return os;
Chris@16:     }
Chris@16: 
Chris@16:     /** Reads a distribution from a @c std::istream */
Chris@16:     BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(
Chris@16:         is, piecewise_linear_distribution, pld)
Chris@16:     {
Chris@16:         param_type parm;
Chris@16:         if(is >> parm) {
Chris@16:             pld.param(parm);
Chris@16:         }
Chris@16:         return is;
Chris@16:     }
Chris@16: 
Chris@16:     /**
Chris@16:      * Returns true if the two distributions will return the
Chris@16:      * same sequence of values, when passed equal generators.
Chris@16:      */
Chris@16:     BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(
Chris@16:         piecewise_linear_distribution, lhs,  rhs)
Chris@16:     {
Chris@16:         return lhs._intervals == rhs._intervals && lhs._weights == rhs._weights;
Chris@16:     }
Chris@16:     /**
Chris@16:      * Returns true if the two distributions may return different
Chris@16:      * sequences of values, when passed equal generators.
Chris@16:      */
Chris@16:     BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(piecewise_linear_distribution)
Chris@16: 
Chris@16: private:
Chris@16: 
Chris@16:     /// @cond \show_private
Chris@16: 
Chris@16:     void init(const std::vector<RealType>& intervals_arg,
Chris@16:               const std::vector<RealType>& weights_arg)
Chris@16:     {
Chris@16:         std::vector<RealType> bin_weights;
Chris@16:         bin_weights.reserve((intervals_arg.size() - 1) * 2);
Chris@16:         for(std::size_t i = 0; i < intervals_arg.size() - 1; ++i) {
Chris@16:             RealType width = intervals_arg[i + 1] - intervals_arg[i];
Chris@16:             RealType w1 = weights_arg[i];
Chris@16:             RealType w2 = weights_arg[i + 1];
Chris@16:             bin_weights.push_back((std::min)(w1, w2) * width);
Chris@16:             bin_weights.push_back(std::abs(w1 - w2) * width / 2);
Chris@16:         }
Chris@16:         typedef discrete_distribution<std::size_t, RealType> bins_type;
Chris@16:         typename bins_type::param_type bins_param(bin_weights);
Chris@16:         _bins.param(bins_param);
Chris@16:     }
Chris@16: 
Chris@16:     void init()
Chris@16:     {
Chris@16:         init(_intervals, _weights);
Chris@16:     }
Chris@16: 
Chris@16:     void default_init()
Chris@16:     {
Chris@16:         _intervals.clear();
Chris@16:         _intervals.push_back(RealType(0));
Chris@16:         _intervals.push_back(RealType(1));
Chris@16:         _weights.clear();
Chris@16:         _weights.push_back(RealType(1));
Chris@16:         _weights.push_back(RealType(1));
Chris@16:         init();
Chris@16:     }
Chris@16: 
Chris@16:     discrete_distribution<std::size_t, RealType> _bins;
Chris@16:     std::vector<RealType> _intervals;
Chris@16:     std::vector<RealType> _weights;
Chris@16: 
Chris@16:     /// @endcond
Chris@16: };
Chris@16: 
Chris@16: }
Chris@16: }
Chris@16: 
Chris@16: #endif