piecewise_constant_distribution.h 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356
  1. //===----------------------------------------------------------------------===//
  2. //
  3. // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
  4. // See https://llvm.org/LICENSE.txt for license information.
  5. // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
  6. //
  7. //===----------------------------------------------------------------------===//
  8. #ifndef _LIBCPP___RANDOM_PIECEWISE_CONSTANT_DISTRIBUTION_H
  9. #define _LIBCPP___RANDOM_PIECEWISE_CONSTANT_DISTRIBUTION_H
  10. #include <__algorithm/upper_bound.h>
  11. #include <__config>
  12. #include <__random/uniform_real_distribution.h>
  13. #include <iosfwd>
  14. #include <numeric>
  15. #include <vector>
  16. #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
  17. # pragma GCC system_header
  18. #endif
  19. _LIBCPP_PUSH_MACROS
  20. #include <__undef_macros>
  21. _LIBCPP_BEGIN_NAMESPACE_STD
  22. template<class _RealType = double>
  23. class _LIBCPP_TEMPLATE_VIS piecewise_constant_distribution
  24. {
  25. public:
  26. // types
  27. typedef _RealType result_type;
  28. class _LIBCPP_TEMPLATE_VIS param_type
  29. {
  30. vector<result_type> __b_;
  31. vector<result_type> __densities_;
  32. vector<result_type> __areas_;
  33. public:
  34. typedef piecewise_constant_distribution distribution_type;
  35. param_type();
  36. template<class _InputIteratorB, class _InputIteratorW>
  37. param_type(_InputIteratorB __fB, _InputIteratorB __lB,
  38. _InputIteratorW __fW);
  39. #ifndef _LIBCPP_CXX03_LANG
  40. template<class _UnaryOperation>
  41. param_type(initializer_list<result_type> __bl, _UnaryOperation __fw);
  42. #endif // _LIBCPP_CXX03_LANG
  43. template<class _UnaryOperation>
  44. param_type(size_t __nw, result_type __xmin, result_type __xmax,
  45. _UnaryOperation __fw);
  46. param_type(param_type const&) = default;
  47. param_type & operator=(const param_type& __rhs);
  48. _LIBCPP_INLINE_VISIBILITY
  49. vector<result_type> intervals() const {return __b_;}
  50. _LIBCPP_INLINE_VISIBILITY
  51. vector<result_type> densities() const {return __densities_;}
  52. friend _LIBCPP_INLINE_VISIBILITY
  53. bool operator==(const param_type& __x, const param_type& __y)
  54. {return __x.__densities_ == __y.__densities_ && __x.__b_ == __y.__b_;}
  55. friend _LIBCPP_INLINE_VISIBILITY
  56. bool operator!=(const param_type& __x, const param_type& __y)
  57. {return !(__x == __y);}
  58. private:
  59. void __init();
  60. friend class piecewise_constant_distribution;
  61. template <class _CharT, class _Traits, class _RT>
  62. friend
  63. basic_ostream<_CharT, _Traits>&
  64. operator<<(basic_ostream<_CharT, _Traits>& __os,
  65. const piecewise_constant_distribution<_RT>& __x);
  66. template <class _CharT, class _Traits, class _RT>
  67. friend
  68. basic_istream<_CharT, _Traits>&
  69. operator>>(basic_istream<_CharT, _Traits>& __is,
  70. piecewise_constant_distribution<_RT>& __x);
  71. };
  72. private:
  73. param_type __p_;
  74. public:
  75. // constructor and reset functions
  76. _LIBCPP_INLINE_VISIBILITY
  77. piecewise_constant_distribution() {}
  78. template<class _InputIteratorB, class _InputIteratorW>
  79. _LIBCPP_INLINE_VISIBILITY
  80. piecewise_constant_distribution(_InputIteratorB __fB,
  81. _InputIteratorB __lB,
  82. _InputIteratorW __fW)
  83. : __p_(__fB, __lB, __fW) {}
  84. #ifndef _LIBCPP_CXX03_LANG
  85. template<class _UnaryOperation>
  86. _LIBCPP_INLINE_VISIBILITY
  87. piecewise_constant_distribution(initializer_list<result_type> __bl,
  88. _UnaryOperation __fw)
  89. : __p_(__bl, __fw) {}
  90. #endif // _LIBCPP_CXX03_LANG
  91. template<class _UnaryOperation>
  92. _LIBCPP_INLINE_VISIBILITY
  93. piecewise_constant_distribution(size_t __nw, result_type __xmin,
  94. result_type __xmax, _UnaryOperation __fw)
  95. : __p_(__nw, __xmin, __xmax, __fw) {}
  96. _LIBCPP_INLINE_VISIBILITY
  97. explicit piecewise_constant_distribution(const param_type& __p)
  98. : __p_(__p) {}
  99. _LIBCPP_INLINE_VISIBILITY
  100. void reset() {}
  101. // generating functions
  102. template<class _URNG>
  103. _LIBCPP_INLINE_VISIBILITY
  104. result_type operator()(_URNG& __g)
  105. {return (*this)(__g, __p_);}
  106. template<class _URNG> result_type operator()(_URNG& __g, const param_type& __p);
  107. // property functions
  108. _LIBCPP_INLINE_VISIBILITY
  109. vector<result_type> intervals() const {return __p_.intervals();}
  110. _LIBCPP_INLINE_VISIBILITY
  111. vector<result_type> densities() const {return __p_.densities();}
  112. _LIBCPP_INLINE_VISIBILITY
  113. param_type param() const {return __p_;}
  114. _LIBCPP_INLINE_VISIBILITY
  115. void param(const param_type& __p) {__p_ = __p;}
  116. _LIBCPP_INLINE_VISIBILITY
  117. result_type min() const {return __p_.__b_.front();}
  118. _LIBCPP_INLINE_VISIBILITY
  119. result_type max() const {return __p_.__b_.back();}
  120. friend _LIBCPP_INLINE_VISIBILITY
  121. bool operator==(const piecewise_constant_distribution& __x,
  122. const piecewise_constant_distribution& __y)
  123. {return __x.__p_ == __y.__p_;}
  124. friend _LIBCPP_INLINE_VISIBILITY
  125. bool operator!=(const piecewise_constant_distribution& __x,
  126. const piecewise_constant_distribution& __y)
  127. {return !(__x == __y);}
  128. template <class _CharT, class _Traits, class _RT>
  129. friend
  130. basic_ostream<_CharT, _Traits>&
  131. operator<<(basic_ostream<_CharT, _Traits>& __os,
  132. const piecewise_constant_distribution<_RT>& __x);
  133. template <class _CharT, class _Traits, class _RT>
  134. friend
  135. basic_istream<_CharT, _Traits>&
  136. operator>>(basic_istream<_CharT, _Traits>& __is,
  137. piecewise_constant_distribution<_RT>& __x);
  138. };
  139. template<class _RealType>
  140. typename piecewise_constant_distribution<_RealType>::param_type &
  141. piecewise_constant_distribution<_RealType>::param_type::operator=
  142. (const param_type& __rhs)
  143. {
  144. // These can throw
  145. __b_.reserve (__rhs.__b_.size ());
  146. __densities_.reserve(__rhs.__densities_.size());
  147. __areas_.reserve (__rhs.__areas_.size());
  148. // These can not throw
  149. __b_ = __rhs.__b_;
  150. __densities_ = __rhs.__densities_;
  151. __areas_ = __rhs.__areas_;
  152. return *this;
  153. }
  154. template<class _RealType>
  155. void
  156. piecewise_constant_distribution<_RealType>::param_type::__init()
  157. {
  158. // __densities_ contains non-normalized areas
  159. result_type __total_area = _VSTD::accumulate(__densities_.begin(),
  160. __densities_.end(),
  161. result_type());
  162. for (size_t __i = 0; __i < __densities_.size(); ++__i)
  163. __densities_[__i] /= __total_area;
  164. // __densities_ contains normalized areas
  165. __areas_.assign(__densities_.size(), result_type());
  166. _VSTD::partial_sum(__densities_.begin(), __densities_.end() - 1,
  167. __areas_.begin() + 1);
  168. // __areas_ contains partial sums of normalized areas: [0, __densities_ - 1]
  169. __densities_.back() = 1 - __areas_.back(); // correct round off error
  170. for (size_t __i = 0; __i < __densities_.size(); ++__i)
  171. __densities_[__i] /= (__b_[__i+1] - __b_[__i]);
  172. // __densities_ now contains __densities_
  173. }
  174. template<class _RealType>
  175. piecewise_constant_distribution<_RealType>::param_type::param_type()
  176. : __b_(2),
  177. __densities_(1, 1.0),
  178. __areas_(1, 0.0)
  179. {
  180. __b_[1] = 1;
  181. }
  182. template<class _RealType>
  183. template<class _InputIteratorB, class _InputIteratorW>
  184. piecewise_constant_distribution<_RealType>::param_type::param_type(
  185. _InputIteratorB __fB, _InputIteratorB __lB, _InputIteratorW __fW)
  186. : __b_(__fB, __lB)
  187. {
  188. if (__b_.size() < 2)
  189. {
  190. __b_.resize(2);
  191. __b_[0] = 0;
  192. __b_[1] = 1;
  193. __densities_.assign(1, 1.0);
  194. __areas_.assign(1, 0.0);
  195. }
  196. else
  197. {
  198. __densities_.reserve(__b_.size() - 1);
  199. for (size_t __i = 0; __i < __b_.size() - 1; ++__i, ++__fW)
  200. __densities_.push_back(*__fW);
  201. __init();
  202. }
  203. }
  204. #ifndef _LIBCPP_CXX03_LANG
  205. template<class _RealType>
  206. template<class _UnaryOperation>
  207. piecewise_constant_distribution<_RealType>::param_type::param_type(
  208. initializer_list<result_type> __bl, _UnaryOperation __fw)
  209. : __b_(__bl.begin(), __bl.end())
  210. {
  211. if (__b_.size() < 2)
  212. {
  213. __b_.resize(2);
  214. __b_[0] = 0;
  215. __b_[1] = 1;
  216. __densities_.assign(1, 1.0);
  217. __areas_.assign(1, 0.0);
  218. }
  219. else
  220. {
  221. __densities_.reserve(__b_.size() - 1);
  222. for (size_t __i = 0; __i < __b_.size() - 1; ++__i)
  223. __densities_.push_back(__fw((__b_[__i+1] + __b_[__i])*.5));
  224. __init();
  225. }
  226. }
  227. #endif // _LIBCPP_CXX03_LANG
  228. template<class _RealType>
  229. template<class _UnaryOperation>
  230. piecewise_constant_distribution<_RealType>::param_type::param_type(
  231. size_t __nw, result_type __xmin, result_type __xmax, _UnaryOperation __fw)
  232. : __b_(__nw == 0 ? 2 : __nw + 1)
  233. {
  234. size_t __n = __b_.size() - 1;
  235. result_type __d = (__xmax - __xmin) / __n;
  236. __densities_.reserve(__n);
  237. for (size_t __i = 0; __i < __n; ++__i)
  238. {
  239. __b_[__i] = __xmin + __i * __d;
  240. __densities_.push_back(__fw(__b_[__i] + __d*.5));
  241. }
  242. __b_[__n] = __xmax;
  243. __init();
  244. }
  245. template<class _RealType>
  246. template<class _URNG>
  247. _RealType
  248. piecewise_constant_distribution<_RealType>::operator()(_URNG& __g, const param_type& __p)
  249. {
  250. typedef uniform_real_distribution<result_type> _Gen;
  251. result_type __u = _Gen()(__g);
  252. ptrdiff_t __k = _VSTD::upper_bound(__p.__areas_.begin(), __p.__areas_.end(),
  253. __u) - __p.__areas_.begin() - 1;
  254. return (__u - __p.__areas_[__k]) / __p.__densities_[__k] + __p.__b_[__k];
  255. }
  256. template <class _CharT, class _Traits, class _RT>
  257. basic_ostream<_CharT, _Traits>&
  258. operator<<(basic_ostream<_CharT, _Traits>& __os,
  259. const piecewise_constant_distribution<_RT>& __x)
  260. {
  261. __save_flags<_CharT, _Traits> __lx(__os);
  262. typedef basic_ostream<_CharT, _Traits> _OStream;
  263. __os.flags(_OStream::dec | _OStream::left | _OStream::fixed |
  264. _OStream::scientific);
  265. _CharT __sp = __os.widen(' ');
  266. __os.fill(__sp);
  267. size_t __n = __x.__p_.__b_.size();
  268. __os << __n;
  269. for (size_t __i = 0; __i < __n; ++__i)
  270. __os << __sp << __x.__p_.__b_[__i];
  271. __n = __x.__p_.__densities_.size();
  272. __os << __sp << __n;
  273. for (size_t __i = 0; __i < __n; ++__i)
  274. __os << __sp << __x.__p_.__densities_[__i];
  275. __n = __x.__p_.__areas_.size();
  276. __os << __sp << __n;
  277. for (size_t __i = 0; __i < __n; ++__i)
  278. __os << __sp << __x.__p_.__areas_[__i];
  279. return __os;
  280. }
  281. template <class _CharT, class _Traits, class _RT>
  282. basic_istream<_CharT, _Traits>&
  283. operator>>(basic_istream<_CharT, _Traits>& __is,
  284. piecewise_constant_distribution<_RT>& __x)
  285. {
  286. typedef piecewise_constant_distribution<_RT> _Eng;
  287. typedef typename _Eng::result_type result_type;
  288. __save_flags<_CharT, _Traits> __lx(__is);
  289. typedef basic_istream<_CharT, _Traits> _Istream;
  290. __is.flags(_Istream::dec | _Istream::skipws);
  291. size_t __n;
  292. __is >> __n;
  293. vector<result_type> __b(__n);
  294. for (size_t __i = 0; __i < __n; ++__i)
  295. __is >> __b[__i];
  296. __is >> __n;
  297. vector<result_type> __densities(__n);
  298. for (size_t __i = 0; __i < __n; ++__i)
  299. __is >> __densities[__i];
  300. __is >> __n;
  301. vector<result_type> __areas(__n);
  302. for (size_t __i = 0; __i < __n; ++__i)
  303. __is >> __areas[__i];
  304. if (!__is.fail())
  305. {
  306. swap(__x.__p_.__b_, __b);
  307. swap(__x.__p_.__densities_, __densities);
  308. swap(__x.__p_.__areas_, __areas);
  309. }
  310. return __is;
  311. }
  312. _LIBCPP_END_NAMESPACE_STD
  313. _LIBCPP_POP_MACROS
  314. #endif // _LIBCPP___RANDOM_PIECEWISE_CONSTANT_DISTRIBUTION_H