File indexing completed on 2025-01-30 09:59:13
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013 #ifndef BOOST_RANDOM_PIECEWISE_LINEAR_DISTRIBUTION_HPP_INCLUDED
0014 #define BOOST_RANDOM_PIECEWISE_LINEAR_DISTRIBUTION_HPP_INCLUDED
0015
0016 #include <vector>
0017 #include <algorithm>
0018 #include <cmath>
0019 #include <cstdlib>
0020 #include <boost/assert.hpp>
0021 #include <boost/random/uniform_real.hpp>
0022 #include <boost/random/discrete_distribution.hpp>
0023 #include <boost/random/detail/config.hpp>
0024 #include <boost/random/detail/operators.hpp>
0025 #include <boost/random/detail/vector_io.hpp>
0026
0027 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
0028 #include <initializer_list>
0029 #endif
0030
0031 #include <boost/range/begin.hpp>
0032 #include <boost/range/end.hpp>
0033
0034 namespace boost {
0035 namespace random {
0036
0037
0038
0039
0040 template<class RealType = double>
0041 class piecewise_linear_distribution {
0042 public:
0043 typedef std::size_t input_type;
0044 typedef RealType result_type;
0045
0046 class param_type {
0047 public:
0048
0049 typedef piecewise_linear_distribution distribution_type;
0050
0051
0052
0053
0054
0055 param_type()
0056 {
0057 _weights.push_back(RealType(1));
0058 _weights.push_back(RealType(1));
0059 _intervals.push_back(RealType(0));
0060 _intervals.push_back(RealType(1));
0061 }
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074 template<class IntervalIter, class WeightIter>
0075 param_type(IntervalIter intervals_first, IntervalIter intervals_last,
0076 WeightIter weight_first)
0077 : _intervals(intervals_first, intervals_last)
0078 {
0079 if(_intervals.size() < 2) {
0080 _intervals.clear();
0081 _weights.push_back(RealType(1));
0082 _weights.push_back(RealType(1));
0083 _intervals.push_back(RealType(0));
0084 _intervals.push_back(RealType(1));
0085 } else {
0086 _weights.reserve(_intervals.size());
0087 for(std::size_t i = 0; i < _intervals.size(); ++i) {
0088 _weights.push_back(*weight_first++);
0089 }
0090 }
0091 }
0092 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104 template<class T, class F>
0105 param_type(const std::initializer_list<T>& il, F f)
0106 : _intervals(il.begin(), il.end())
0107 {
0108 if(_intervals.size() < 2) {
0109 _intervals.clear();
0110 _weights.push_back(RealType(1));
0111 _weights.push_back(RealType(1));
0112 _intervals.push_back(RealType(0));
0113 _intervals.push_back(RealType(1));
0114 } else {
0115 _weights.reserve(_intervals.size());
0116 for(typename std::vector<RealType>::const_iterator
0117 iter = _intervals.begin(), end = _intervals.end();
0118 iter != end; ++iter)
0119 {
0120 _weights.push_back(f(*iter));
0121 }
0122 }
0123 }
0124 #endif
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134 template<class IntervalRange, class WeightRange>
0135 param_type(const IntervalRange& intervals_arg,
0136 const WeightRange& weights_arg)
0137 : _intervals(boost::begin(intervals_arg), boost::end(intervals_arg)),
0138 _weights(boost::begin(weights_arg), boost::end(weights_arg))
0139 {
0140 if(_intervals.size() < 2) {
0141 _weights.clear();
0142 _weights.push_back(RealType(1));
0143 _weights.push_back(RealType(1));
0144 _intervals.clear();
0145 _intervals.push_back(RealType(0));
0146 _intervals.push_back(RealType(1));
0147 }
0148 }
0149
0150
0151
0152
0153
0154
0155
0156
0157 template<class F>
0158 param_type(std::size_t nw, RealType xmin, RealType xmax, F f)
0159 {
0160 std::size_t n = (nw == 0) ? 1 : nw;
0161 double delta = (xmax - xmin) / n;
0162 BOOST_ASSERT(delta > 0);
0163 for(std::size_t k = 0; k < n; ++k) {
0164 _weights.push_back(f(xmin + k*delta));
0165 _intervals.push_back(xmin + k*delta);
0166 }
0167 _weights.push_back(f(xmax));
0168 _intervals.push_back(xmax);
0169 }
0170
0171
0172 std::vector<RealType> intervals() const { return _intervals; }
0173
0174
0175
0176
0177
0178 std::vector<RealType> densities() const
0179 {
0180 RealType sum = static_cast<RealType>(0);
0181 for(std::size_t i = 0; i < _intervals.size() - 1; ++i) {
0182 RealType width = _intervals[i + 1] - _intervals[i];
0183 sum += (_weights[i] + _weights[i + 1]) * width / 2;
0184 }
0185 std::vector<RealType> result;
0186 result.reserve(_weights.size());
0187 for(typename std::vector<RealType>::const_iterator
0188 iter = _weights.begin(), end = _weights.end();
0189 iter != end; ++iter)
0190 {
0191 result.push_back(*iter / sum);
0192 }
0193 return result;
0194 }
0195
0196
0197 BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
0198 {
0199 detail::print_vector(os, parm._intervals);
0200 detail::print_vector(os, parm._weights);
0201 return os;
0202 }
0203
0204
0205 BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
0206 {
0207 std::vector<RealType> new_intervals;
0208 std::vector<RealType> new_weights;
0209 detail::read_vector(is, new_intervals);
0210 detail::read_vector(is, new_weights);
0211 if(is) {
0212 parm._intervals.swap(new_intervals);
0213 parm._weights.swap(new_weights);
0214 }
0215 return is;
0216 }
0217
0218
0219 BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
0220 {
0221 return lhs._intervals == rhs._intervals
0222 && lhs._weights == rhs._weights;
0223 }
0224
0225 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
0226
0227 private:
0228 friend class piecewise_linear_distribution;
0229
0230 std::vector<RealType> _intervals;
0231 std::vector<RealType> _weights;
0232 };
0233
0234
0235
0236
0237
0238 piecewise_linear_distribution()
0239 {
0240 default_init();
0241 }
0242
0243
0244
0245
0246
0247
0248
0249
0250
0251
0252
0253
0254
0255
0256
0257
0258
0259
0260
0261
0262
0263
0264
0265 template<class IntervalIter, class WeightIter>
0266 piecewise_linear_distribution(IntervalIter first_interval,
0267 IntervalIter last_interval,
0268 WeightIter first_weight)
0269 : _intervals(first_interval, last_interval)
0270 {
0271 if(_intervals.size() < 2) {
0272 default_init();
0273 } else {
0274 _weights.reserve(_intervals.size());
0275 for(std::size_t i = 0; i < _intervals.size(); ++i) {
0276 _weights.push_back(*first_weight++);
0277 }
0278 init();
0279 }
0280 }
0281 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
0282
0283
0284
0285
0286
0287
0288
0289
0290
0291
0292
0293
0294 template<class T, class F>
0295 piecewise_linear_distribution(std::initializer_list<T> il, F f)
0296 : _intervals(il.begin(), il.end())
0297 {
0298 if(_intervals.size() < 2) {
0299 default_init();
0300 } else {
0301 _weights.reserve(_intervals.size());
0302 for(typename std::vector<RealType>::const_iterator
0303 iter = _intervals.begin(), end = _intervals.end();
0304 iter != end; ++iter)
0305 {
0306 _weights.push_back(f(*iter));
0307 }
0308 init();
0309 }
0310 }
0311 #endif
0312
0313
0314
0315
0316
0317
0318
0319
0320
0321 template<class IntervalsRange, class WeightsRange>
0322 piecewise_linear_distribution(const IntervalsRange& intervals_arg,
0323 const WeightsRange& weights_arg)
0324 : _intervals(boost::begin(intervals_arg), boost::end(intervals_arg)),
0325 _weights(boost::begin(weights_arg), boost::end(weights_arg))
0326 {
0327 if(_intervals.size() < 2) {
0328 default_init();
0329 } else {
0330 init();
0331 }
0332 }
0333
0334
0335
0336
0337
0338
0339 template<class F>
0340 piecewise_linear_distribution(std::size_t nw,
0341 RealType xmin,
0342 RealType xmax,
0343 F f)
0344 {
0345 if(nw == 0) { nw = 1; }
0346 RealType delta = (xmax - xmin) / nw;
0347 _intervals.reserve(nw + 1);
0348 for(std::size_t i = 0; i < nw; ++i) {
0349 RealType x = xmin + i * delta;
0350 _intervals.push_back(x);
0351 _weights.push_back(f(x));
0352 }
0353 _intervals.push_back(xmax);
0354 _weights.push_back(f(xmax));
0355 init();
0356 }
0357
0358
0359
0360 explicit piecewise_linear_distribution(const param_type& parm)
0361 : _intervals(parm._intervals),
0362 _weights(parm._weights)
0363 {
0364 init();
0365 }
0366
0367
0368
0369
0370
0371 template<class URNG>
0372 RealType operator()(URNG& urng) const
0373 {
0374 std::size_t i = _bins(urng);
0375 bool is_in_rectangle = (i % 2 == 0);
0376 i = i / 2;
0377 uniform_real<RealType> dist(_intervals[i], _intervals[i+1]);
0378 if(is_in_rectangle) {
0379 return dist(urng);
0380 } else if(_weights[i] < _weights[i+1]) {
0381 return (std::max)(dist(urng), dist(urng));
0382 } else {
0383 return (std::min)(dist(urng), dist(urng));
0384 }
0385 }
0386
0387
0388
0389
0390
0391 template<class URNG>
0392 RealType operator()(URNG& urng, const param_type& parm) const
0393 {
0394 return piecewise_linear_distribution(parm)(urng);
0395 }
0396
0397
0398 result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const
0399 { return _intervals.front(); }
0400
0401 result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const
0402 { return _intervals.back(); }
0403
0404
0405
0406
0407
0408 std::vector<RealType> densities() const
0409 {
0410 RealType sum = static_cast<RealType>(0);
0411 for(std::size_t i = 0; i < _intervals.size() - 1; ++i) {
0412 RealType width = _intervals[i + 1] - _intervals[i];
0413 sum += (_weights[i] + _weights[i + 1]) * width / 2;
0414 }
0415 std::vector<RealType> result;
0416 result.reserve(_weights.size());
0417 for(typename std::vector<RealType>::const_iterator
0418 iter = _weights.begin(), end = _weights.end();
0419 iter != end; ++iter)
0420 {
0421 result.push_back(*iter / sum);
0422 }
0423 return result;
0424 }
0425
0426 std::vector<RealType> intervals() const { return _intervals; }
0427
0428
0429 param_type param() const
0430 {
0431 return param_type(_intervals, _weights);
0432 }
0433
0434 void param(const param_type& parm)
0435 {
0436 std::vector<RealType> new_intervals(parm._intervals);
0437 std::vector<RealType> new_weights(parm._weights);
0438 init(new_intervals, new_weights);
0439 _intervals.swap(new_intervals);
0440 _weights.swap(new_weights);
0441 }
0442
0443
0444
0445
0446
0447 void reset() { _bins.reset(); }
0448
0449
0450 BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(
0451 os, piecewise_linear_distribution, pld)
0452 {
0453 os << pld.param();
0454 return os;
0455 }
0456
0457
0458 BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(
0459 is, piecewise_linear_distribution, pld)
0460 {
0461 param_type parm;
0462 if(is >> parm) {
0463 pld.param(parm);
0464 }
0465 return is;
0466 }
0467
0468
0469
0470
0471
0472 BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(
0473 piecewise_linear_distribution, lhs, rhs)
0474 {
0475 return lhs._intervals == rhs._intervals && lhs._weights == rhs._weights;
0476 }
0477
0478
0479
0480
0481 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(piecewise_linear_distribution)
0482
0483 private:
0484
0485
0486
0487 void init(const std::vector<RealType>& intervals_arg,
0488 const std::vector<RealType>& weights_arg)
0489 {
0490 using std::abs;
0491 std::vector<RealType> bin_weights;
0492 bin_weights.reserve((intervals_arg.size() - 1) * 2);
0493 for(std::size_t i = 0; i < intervals_arg.size() - 1; ++i) {
0494 RealType width = intervals_arg[i + 1] - intervals_arg[i];
0495 RealType w1 = weights_arg[i];
0496 RealType w2 = weights_arg[i + 1];
0497 bin_weights.push_back((std::min)(w1, w2) * width);
0498 bin_weights.push_back(abs(w1 - w2) * width / 2);
0499 }
0500 typedef discrete_distribution<std::size_t, RealType> bins_type;
0501 typename bins_type::param_type bins_param(bin_weights);
0502 _bins.param(bins_param);
0503 }
0504
0505 void init()
0506 {
0507 init(_intervals, _weights);
0508 }
0509
0510 void default_init()
0511 {
0512 _intervals.clear();
0513 _intervals.push_back(RealType(0));
0514 _intervals.push_back(RealType(1));
0515 _weights.clear();
0516 _weights.push_back(RealType(1));
0517 _weights.push_back(RealType(1));
0518 init();
0519 }
0520
0521 discrete_distribution<std::size_t, RealType> _bins;
0522 std::vector<RealType> _intervals;
0523 std::vector<RealType> _weights;
0524
0525
0526 };
0527
0528 }
0529 }
0530
0531 #endif