boost/random/piecewise_constant_distribution.hpp

   1 /* boost random/piecewise_constant_distribution.hpp header file
   2  *
   3  * Copyright Steven Watanabe 2011
   4  * Distributed under the Boost Software License, Version 1.0. (See
   5  * accompanying file LICENSE_1_0.txt or copy at
   6  * http://www.boost.org/LICENSE_1_0.txt)
   7  *
   8  * See http://www.boost.org for most recent version including documentation.
   9  *
  10  * $Id: piecewise_constant_distribution.hpp 85813 2013-09-21 20:17:00Z jewillco $
  11  */
  12
  13 #ifndef BOOST_RANDOM_PIECEWISE_CONSTANT_DISTRIBUTION_HPP_INCLUDED
  14 #define BOOST_RANDOM_PIECEWISE_CONSTANT_DISTRIBUTION_HPP_INCLUDED
  15
  16 #include <vector>
  17 #include <numeric>
  18 #include <boost/assert.hpp>
  19 #include <boost/random/uniform_real.hpp>
  20 #include <boost/random/discrete_distribution.hpp>
  21 #include <boost/random/detail/config.hpp>
  22 #include <boost/random/detail/operators.hpp>
  23 #include <boost/random/detail/vector_io.hpp>
  24
  25 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
  26 #include <initializer_list>
  27 #endif
  28
  29 #include <boost/range/begin.hpp>
  30 #include <boost/range/end.hpp>
  31
  32 namespace boost {
  33 namespace random {
  34
  35 /**
  36  * The class @c piecewise_constant_distribution models a \random_distribution.
  37  */
  38 template<class RealType = double, class WeightType = double>
  39 class piecewise_constant_distribution {
  40 public:
  41     typedef std::size_t input_type;
  42     typedef RealType result_type;
  43
  44     class param_type {
  45     public:
  46
  47         typedef piecewise_constant_distribution distribution_type;
  48
  49         /**
  50          * Constructs a @c param_type object, representing a distribution
  51          * that produces values uniformly distributed in the range [0, 1).
  52          */
  53         param_type()
  54         {
  55             _weights.push_back(WeightType(1));
  56             _intervals.push_back(RealType(0));
  57             _intervals.push_back(RealType(1));
  58         }
  59         /**
  60          * Constructs a @c param_type object from two iterator ranges
  61          * containing the interval boundaries and the interval weights.
  62          * If there are less than two boundaries, then this is equivalent to
  63          * the default constructor and creates a single interval, [0, 1).
  64          *
  65          * The values of the interval boundaries must be strictly
  66          * increasing, and the number of weights must be one less than
  67          * the number of interval boundaries.  If there are extra
  68          * weights, they are ignored.
  69          */
  70         template<class IntervalIter, class WeightIter>
  71         param_type(IntervalIter intervals_first, IntervalIter intervals_last,
  72                    WeightIter weight_first)
  73           : _intervals(intervals_first, intervals_last)
  74         {
  75             if(_intervals.size() < 2) {
  76                 _intervals.clear();
  77                 _intervals.push_back(RealType(0));
  78                 _intervals.push_back(RealType(1));
  79                 _weights.push_back(WeightType(1));
  80             } else {
  81                 _weights.reserve(_intervals.size() - 1);
  82                 for(std::size_t i = 0; i < _intervals.size() - 1; ++i) {
  83                     _weights.push_back(*weight_first++);
  84                 }
  85             }
  86         }
  87 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
  88         /**
  89          * Constructs a @c param_type object from an
  90          * initializer_list containing the interval boundaries
  91          * and a unary function specifying the weights.  Each
  92          * weight is determined by calling the function at the
  93          * midpoint of the corresponding interval.
  94          *
  95          * If the initializer_list contains less than two elements,
  96          * this is equivalent to the default constructor and the
  97          * distribution will produce values uniformly distributed
  98          * in the range [0, 1).
  99          */
 100         template<class T, class F>
 101         param_type(const std::initializer_list<T>& il, F f)
 102           : _intervals(il.begin(), il.end())
 103         {
 104             if(_intervals.size() < 2) {
 105                 _intervals.clear();
 106                 _intervals.push_back(RealType(0));
 107                 _intervals.push_back(RealType(1));
 108                 _weights.push_back(WeightType(1));
 109             } else {
 110                 _weights.reserve(_intervals.size() - 1);
 111                 for(std::size_t i = 0; i < _intervals.size() - 1; ++i) {
 112                     RealType midpoint = (_intervals[i] + _intervals[i + 1]) / 2;
 113                     _weights.push_back(f(midpoint));
 114                 }
 115             }
 116         }
 117 #endif
 118         /**
 119          * Constructs a @c param_type object from Boost.Range
 120          * ranges holding the interval boundaries and the weights.  If
 121          * there are less than two interval boundaries, this is equivalent
 122          * to the default constructor and the distribution will produce
 123          * values uniformly distributed in the range [0, 1).  The
 124          * number of weights must be one less than the number of
 125          * interval boundaries.
 126          */
 127         template<class IntervalRange, class WeightRange>
 128         param_type(const IntervalRange& intervals_arg,
 129                    const WeightRange& weights_arg)
 130           : _intervals(boost::begin(intervals_arg), boost::end(intervals_arg)),
 131             _weights(boost::begin(weights_arg), boost::end(weights_arg))
 132         {
 133             if(_intervals.size() < 2) {
 134                 _intervals.clear();
 135                 _intervals.push_back(RealType(0));
 136                 _intervals.push_back(RealType(1));
 137                 _weights.push_back(WeightType(1));
 138             }
 139         }
 140
 141         /**
 142          * Constructs the parameters for a distribution that approximates a
 143          * function.  The range of the distribution is [xmin, xmax).  This
 144          * range is divided into nw equally sized intervals and the weights
 145          * are found by calling the unary function f on the midpoints of the
 146          * intervals.
 147          */
 148         template<class F>
 149         param_type(std::size_t nw, RealType xmin, RealType xmax, F f)
 150         {
 151             std::size_t n = (nw == 0) ? 1 : nw;
 152             double delta = (xmax - xmin) / n;
 153             BOOST_ASSERT(delta > 0);
 154             for(std::size_t k = 0; k < n; ++k) {
 155                 _weights.push_back(f(xmin + k*delta + delta/2));
 156                 _intervals.push_back(xmin + k*delta);
 157             }
 158             _intervals.push_back(xmax);
 159         }
 160
 161         /**  Returns a vector containing the interval boundaries. */
 162         std::vector<RealType> intervals() const { return _intervals; }
 163
 164         /**
 165          * Returns a vector containing the probability densities
 166          * over all the intervals of the distribution.
 167          */
 168         std::vector<RealType> densities() const
 169         {
 170             RealType sum = std::accumulate(_weights.begin(), _weights.end(),
 171                                              static_cast<RealType>(0));
 172             std::vector<RealType> result;
 173             result.reserve(_weights.size());
 174             for(std::size_t i = 0; i < _weights.size(); ++i) {
 175                 RealType width = _intervals[i + 1] - _intervals[i];
 176                 result.push_back(_weights[i] / (sum * width));
 177             }
 178             return result;
 179         }
 180
 181         /** Writes the parameters to a @c std::ostream. */
 182         BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
 183         {
 184             detail::print_vector(os, parm._intervals);
 185             detail::print_vector(os, parm._weights);
 186             return os;
 187         }
 188
 189         /** Reads the parameters from a @c std::istream. */
 190         BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
 191         {
 192             std::vector<RealType> new_intervals;
 193             std::vector<WeightType> new_weights;
 194             detail::read_vector(is, new_intervals);
 195             detail::read_vector(is, new_weights);
 196             if(is) {
 197                 parm._intervals.swap(new_intervals);
 198                 parm._weights.swap(new_weights);
 199             }
 200             return is;
 201         }
 202
 203         /** Returns true if the two sets of parameters are the same. */
 204         BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
 205         {
 206             return lhs._intervals == rhs._intervals
 207                 && lhs._weights == rhs._weights;
 208         }
 209         /** Returns true if the two sets of parameters are different. */
 210         BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
 211
 212     private:
 213
 214         friend class piecewise_constant_distribution;
 215
 216         std::vector<RealType> _intervals;
 217         std::vector<WeightType> _weights;
 218     };
 219
 220     /**
 221      * Creates a new @c piecewise_constant_distribution with
 222      * a single interval, [0, 1).
 223      */
 224     piecewise_constant_distribution()
 225     {
 226         _intervals.push_back(RealType(0));
 227         _intervals.push_back(RealType(1));
 228     }
 229     /**
 230      * Constructs a piecewise_constant_distribution from two iterator ranges
 231      * containing the interval boundaries and the interval weights.
 232      * If there are less than two boundaries, then this is equivalent to
 233      * the default constructor and creates a single interval, [0, 1).
 234      *
 235      * The values of the interval boundaries must be strictly
 236      * increasing, and the number of weights must be one less than
 237      * the number of interval boundaries.  If there are extra
 238      * weights, they are ignored.
 239      *
 240      * For example,
 241      *
 242      * @code
 243      * double intervals[] = { 0.0, 1.0, 4.0 };
 244      * double weights[] = { 1.0, 1.0 };
 245      * piecewise_constant_distribution<> dist(
 246      *     &intervals[0], &intervals[0] + 3, &weights[0]);
 247      * @endcode
 248      *
 249      * The distribution has a 50% chance of producing a
 250      * value between 0 and 1 and a 50% chance of producing
 251      * a value between 1 and 4.
 252      */
 253     template<class IntervalIter, class WeightIter>
 254     piecewise_constant_distribution(IntervalIter first_interval,
 255                                     IntervalIter last_interval,
 256                                     WeightIter first_weight)
 257       : _intervals(first_interval, last_interval)
 258     {
 259         if(_intervals.size() < 2) {
 260             _intervals.clear();
 261             _intervals.push_back(RealType(0));
 262             _intervals.push_back(RealType(1));
 263         } else {
 264             std::vector<WeightType> actual_weights;
 265             actual_weights.reserve(_intervals.size() - 1);
 266             for(std::size_t i = 0; i < _intervals.size() - 1; ++i) {
 267                 actual_weights.push_back(*first_weight++);
 268             }
 269             typedef discrete_distribution<std::size_t, WeightType> bins_type;
 270             typename bins_type::param_type bins_param(actual_weights);
 271             _bins.param(bins_param);
 272         }
 273     }
 274 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
 275     /**
 276      * Constructs a piecewise_constant_distribution from an
 277      * initializer_list containing the interval boundaries
 278      * and a unary function specifying the weights.  Each
 279      * weight is determined by calling the function at the
 280      * midpoint of the corresponding interval.
 281      *
 282      * If the initializer_list contains less than two elements,
 283      * this is equivalent to the default constructor and the
 284      * distribution will produce values uniformly distributed
 285      * in the range [0, 1).
 286      */
 287     template<class T, class F>
 288     piecewise_constant_distribution(std::initializer_list<T> il, F f)
 289       : _intervals(il.begin(), il.end())
 290     {
 291         if(_intervals.size() < 2) {
 292             _intervals.clear();
 293             _intervals.push_back(RealType(0));
 294             _intervals.push_back(RealType(1));
 295         } else {
 296             std::vector<WeightType> actual_weights;
 297             actual_weights.reserve(_intervals.size() - 1);
 298             for(std::size_t i = 0; i < _intervals.size() - 1; ++i) {
 299                 RealType midpoint = (_intervals[i] + _intervals[i + 1]) / 2;
 300                 actual_weights.push_back(f(midpoint));
 301             }
 302             typedef discrete_distribution<std::size_t, WeightType> bins_type;
 303             typename bins_type::param_type bins_param(actual_weights);
 304             _bins.param(bins_param);
 305         }
 306     }
 307 #endif
 308     /**
 309      * Constructs a piecewise_constant_distribution from Boost.Range
 310      * ranges holding the interval boundaries and the weights.  If
 311      * there are less than two interval boundaries, this is equivalent
 312      * to the default constructor and the distribution will produce
 313      * values uniformly distributed in the range [0, 1).  The
 314      * number of weights must be one less than the number of
 315      * interval boundaries.
 316      */
 317     template<class IntervalsRange, class WeightsRange>
 318     piecewise_constant_distribution(const IntervalsRange& intervals_arg,
 319                                     const WeightsRange& weights_arg)
 320       : _bins(weights_arg),
 321         _intervals(boost::begin(intervals_arg), boost::end(intervals_arg))
 322     {
 323         if(_intervals.size() < 2) {
 324             _intervals.clear();
 325             _intervals.push_back(RealType(0));
 326             _intervals.push_back(RealType(1));
 327         }
 328     }
 329     /**
 330      * Constructs a piecewise_constant_distribution that approximates a
 331      * function.  The range of the distribution is [xmin, xmax).  This
 332      * range is divided into nw equally sized intervals and the weights
 333      * are found by calling the unary function f on the midpoints of the
 334      * intervals.
 335      */
 336     template<class F>
 337     piecewise_constant_distribution(std::size_t nw,
 338                                     RealType xmin,
 339                                     RealType xmax,
 340                                     F f)
 341       : _bins(nw, xmin, xmax, f)
 342     {
 343         if(nw == 0) { nw = 1; }
 344         RealType delta = (xmax - xmin) / nw;
 345         _intervals.reserve(nw + 1);
 346         for(std::size_t i = 0; i < nw; ++i) {
 347             _intervals.push_back(xmin + i * delta);
 348         }
 349         _intervals.push_back(xmax);
 350     }
 351     /**
 352      * Constructs a piecewise_constant_distribution from its parameters.
 353      */
 354     explicit piecewise_constant_distribution(const param_type& parm)
 355       : _bins(parm._weights),
 356         _intervals(parm._intervals)
 357     {
 358     }
 359
 360     /**
 361      * Returns a value distributed according to the parameters of the
 362      * piecewist_constant_distribution.
 363      */
 364     template<class URNG>
 365     RealType operator()(URNG& urng) const
 366     {
 367         std::size_t i = _bins(urng);
 368         return uniform_real<RealType>(_intervals[i], _intervals[i+1])(urng);
 369     }
 370
 371     /**
 372      * Returns a value distributed according to the parameters
 373      * specified by param.
 374      */
 375     template<class URNG>
 376     RealType operator()(URNG& urng, const param_type& parm) const
 377     {
 378         return piecewise_constant_distribution(parm)(urng);
 379     }
 380
 381     /** Returns the smallest value that the distribution can produce. */
 382     result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const
 383     { return _intervals.front(); }
 384     /** Returns the largest value that the distribution can produce. */
 385     result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const
 386     { return _intervals.back(); }
 387
 388     /**
 389      * Returns a vector containing the probability density
 390      * over each interval.
 391      */
 392     std::vector<RealType> densities() const
 393     {
 394         std::vector<RealType> result(_bins.probabilities());
 395         for(std::size_t i = 0; i < result.size(); ++i) {
 396             result[i] /= (_intervals[i+1] - _intervals[i]);
 397         }
 398         return(result);
 399     }
 400     /**  Returns a vector containing the interval boundaries. */
 401     std::vector<RealType> intervals() const { return _intervals; }
 402
 403     /** Returns the parameters of the distribution. */
 404     param_type param() const
 405     {
 406         return param_type(_intervals, _bins.probabilities());
 407     }
 408     /** Sets the parameters of the distribution. */
 409     void param(const param_type& parm)
 410     {
 411         std::vector<RealType> new_intervals(parm._intervals);
 412         typedef discrete_distribution<std::size_t, RealType> bins_type;
 413         typename bins_type::param_type bins_param(parm._weights);
 414         _bins.param(bins_param);
 415         _intervals.swap(new_intervals);
 416     }
 417
 418     /**
 419      * Effects: Subsequent uses of the distribution do not depend
 420      * on values produced by any engine prior to invoking reset.
 421      */
 422     void reset() { _bins.reset(); }
 423
 424     /** Writes a distribution to a @c std::ostream. */
 425     BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(
 426         os, piecewise_constant_distribution, pcd)
 427     {
 428         os << pcd.param();
 429         return os;
 430     }
 431
 432     /** Reads a distribution from a @c std::istream */
 433     BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(
 434         is, piecewise_constant_distribution, pcd)
 435     {
 436         param_type parm;
 437         if(is >> parm) {
 438             pcd.param(parm);
 439         }
 440         return is;
 441     }
 442
 443     /**
 444      * Returns true if the two distributions will return the
 445      * same sequence of values, when passed equal generators.
 446      */
 447     BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(
 448         piecewise_constant_distribution, lhs,  rhs)
 449     {
 450         return lhs._bins == rhs._bins && lhs._intervals == rhs._intervals;
 451     }
 452     /**
 453      * Returns true if the two distributions may return different
 454      * sequences of values, when passed equal generators.
 455      */
 456     BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(piecewise_constant_distribution)
 457
 458 private:
 459     discrete_distribution<std::size_t, WeightType> _bins;
 460     std::vector<RealType> _intervals;
 461 };
 462
 463 }
 464 }
 465
 466 #endif