piecewise_constant_distribution.h 12 KB

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