piecewise_constant_distribution.hpp 17 KB

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