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