piecewise_constant_distribution.hpp 17 KB

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