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