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$ 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 */ param_type()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> param_type(IntervalIter intervals_first,IntervalIter intervals_last,WeightIter weight_first)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> param_type(const std::initializer_list<T> & il,F 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> param_type(const IntervalRange & intervals_arg,const WeightRange & weights_arg)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> param_type(std::size_t nw,RealType xmin,RealType xmax,F 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. */ intervals() const162 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 */ densities() const168 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. */ BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os,param_type,parm)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. */ BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is,param_type,parm)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. */ BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type,lhs,rhs)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 */ piecewise_constant_distribution()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> piecewise_constant_distribution(IntervalIter first_interval,IntervalIter last_interval,WeightIter first_weight)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> piecewise_constant_distribution(std::initializer_list<T> il,F 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> piecewise_constant_distribution(const IntervalsRange & intervals_arg,const WeightsRange & weights_arg)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> piecewise_constant_distribution(std::size_t nw,RealType xmin,RealType xmax,F 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 */ piecewise_constant_distribution(const param_type & parm)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> operator ()(URNG & urng) const365 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> operator ()(URNG & urng,const param_type & parm) const376 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. */ BOOST_PREVENT_MACRO_SUBSTITUTION() const382 result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const 383 { return _intervals.front(); } 384 /** Returns the largest value that the distribution can produce. */ BOOST_PREVENT_MACRO_SUBSTITUTION() const385 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 */ densities() const392 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. */ intervals() const401 std::vector<RealType> intervals() const { return _intervals; } 402 403 /** Returns the parameters of the distribution. */ param() const404 param_type param() const 405 { 406 return param_type(_intervals, _bins.probabilities()); 407 } 408 /** Sets the parameters of the distribution. */ param(const param_type & parm)409 void param(const param_type& parm) 410 { 411 std::vector<RealType> new_intervals(parm._intervals); 412 typedef discrete_distribution<std::size_t, WeightType> 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 */ reset()422 void reset() { _bins.reset(); } 423 424 /** Writes a distribution to a @c std::ostream. */ BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os,piecewise_constant_distribution,pcd)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 */ BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is,piecewise_constant_distribution,pcd)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 */ BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(piecewise_constant_distribution,lhs,rhs)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 467