boost/random/piecewise_linear_distribution.hpp

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