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