piecewise_constant_distribution.h 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302
  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. static_assert(__libcpp_random_is_valid_realtype<_RealType>::value,
  26. "RealType must be a supported floating-point type");
  27. public:
  28. // types
  29. typedef _RealType result_type;
  30. class _LIBCPP_TEMPLATE_VIS param_type {
  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, _InputIteratorW __f_w);
  39. #ifndef _LIBCPP_CXX03_LANG
  40. template <class _UnaryOperation>
  41. _LIBCPP_HIDE_FROM_ABI param_type(initializer_list<result_type> __bl, _UnaryOperation __fw);
  42. #endif // _LIBCPP_CXX03_LANG
  43. template <class _UnaryOperation>
  44. _LIBCPP_HIDE_FROM_ABI param_type(size_t __nw, result_type __xmin, result_type __xmax, _UnaryOperation __fw);
  45. _LIBCPP_HIDE_FROM_ABI param_type(param_type const&) = default;
  46. _LIBCPP_HIDE_FROM_ABI param_type& operator=(const param_type& __rhs);
  47. _LIBCPP_HIDE_FROM_ABI vector<result_type> intervals() const { return __b_; }
  48. _LIBCPP_HIDE_FROM_ABI vector<result_type> densities() const { return __densities_; }
  49. friend _LIBCPP_HIDE_FROM_ABI bool operator==(const param_type& __x, const param_type& __y) {
  50. return __x.__densities_ == __y.__densities_ && __x.__b_ == __y.__b_;
  51. }
  52. friend _LIBCPP_HIDE_FROM_ABI bool operator!=(const param_type& __x, const param_type& __y) { return !(__x == __y); }
  53. private:
  54. _LIBCPP_HIDE_FROM_ABI void __init();
  55. friend class piecewise_constant_distribution;
  56. template <class _CharT, class _Traits, class _RT>
  57. friend basic_ostream<_CharT, _Traits>&
  58. operator<<(basic_ostream<_CharT, _Traits>& __os, const piecewise_constant_distribution<_RT>& __x);
  59. template <class _CharT, class _Traits, class _RT>
  60. friend basic_istream<_CharT, _Traits>&
  61. operator>>(basic_istream<_CharT, _Traits>& __is, piecewise_constant_distribution<_RT>& __x);
  62. };
  63. private:
  64. param_type __p_;
  65. public:
  66. // constructor and reset functions
  67. _LIBCPP_HIDE_FROM_ABI piecewise_constant_distribution() {}
  68. template <class _InputIteratorB, class _InputIteratorW>
  69. _LIBCPP_HIDE_FROM_ABI
  70. piecewise_constant_distribution(_InputIteratorB __f_b, _InputIteratorB __l_b, _InputIteratorW __f_w)
  71. : __p_(__f_b, __l_b, __f_w) {}
  72. #ifndef _LIBCPP_CXX03_LANG
  73. template <class _UnaryOperation>
  74. _LIBCPP_HIDE_FROM_ABI piecewise_constant_distribution(initializer_list<result_type> __bl, _UnaryOperation __fw)
  75. : __p_(__bl, __fw) {}
  76. #endif // _LIBCPP_CXX03_LANG
  77. template <class _UnaryOperation>
  78. _LIBCPP_HIDE_FROM_ABI
  79. piecewise_constant_distribution(size_t __nw, result_type __xmin, result_type __xmax, _UnaryOperation __fw)
  80. : __p_(__nw, __xmin, __xmax, __fw) {}
  81. _LIBCPP_HIDE_FROM_ABI explicit piecewise_constant_distribution(const param_type& __p) : __p_(__p) {}
  82. _LIBCPP_HIDE_FROM_ABI void reset() {}
  83. // generating functions
  84. template <class _URNG>
  85. _LIBCPP_HIDE_FROM_ABI result_type operator()(_URNG& __g) {
  86. return (*this)(__g, __p_);
  87. }
  88. template <class _URNG>
  89. _LIBCPP_HIDE_FROM_ABI result_type operator()(_URNG& __g, const param_type& __p);
  90. // property functions
  91. _LIBCPP_HIDE_FROM_ABI vector<result_type> intervals() const { return __p_.intervals(); }
  92. _LIBCPP_HIDE_FROM_ABI vector<result_type> densities() const { return __p_.densities(); }
  93. _LIBCPP_HIDE_FROM_ABI param_type param() const { return __p_; }
  94. _LIBCPP_HIDE_FROM_ABI void param(const param_type& __p) { __p_ = __p; }
  95. _LIBCPP_HIDE_FROM_ABI result_type min() const { return __p_.__b_.front(); }
  96. _LIBCPP_HIDE_FROM_ABI result_type max() const { return __p_.__b_.back(); }
  97. friend _LIBCPP_HIDE_FROM_ABI bool
  98. operator==(const piecewise_constant_distribution& __x, const piecewise_constant_distribution& __y) {
  99. return __x.__p_ == __y.__p_;
  100. }
  101. friend _LIBCPP_HIDE_FROM_ABI bool
  102. operator!=(const piecewise_constant_distribution& __x, const piecewise_constant_distribution& __y) {
  103. return !(__x == __y);
  104. }
  105. template <class _CharT, class _Traits, class _RT>
  106. friend basic_ostream<_CharT, _Traits>&
  107. operator<<(basic_ostream<_CharT, _Traits>& __os, const piecewise_constant_distribution<_RT>& __x);
  108. template <class _CharT, class _Traits, class _RT>
  109. friend basic_istream<_CharT, _Traits>&
  110. operator>>(basic_istream<_CharT, _Traits>& __is, piecewise_constant_distribution<_RT>& __x);
  111. };
  112. template <class _RealType>
  113. typename piecewise_constant_distribution<_RealType>::param_type&
  114. piecewise_constant_distribution<_RealType>::param_type::operator=(const param_type& __rhs) {
  115. // These can throw
  116. __b_.reserve(__rhs.__b_.size());
  117. __densities_.reserve(__rhs.__densities_.size());
  118. __areas_.reserve(__rhs.__areas_.size());
  119. // These can not throw
  120. __b_ = __rhs.__b_;
  121. __densities_ = __rhs.__densities_;
  122. __areas_ = __rhs.__areas_;
  123. return *this;
  124. }
  125. template <class _RealType>
  126. void piecewise_constant_distribution<_RealType>::param_type::__init() {
  127. // __densities_ contains non-normalized areas
  128. result_type __total_area = std::accumulate(__densities_.begin(), __densities_.end(), result_type());
  129. for (size_t __i = 0; __i < __densities_.size(); ++__i)
  130. __densities_[__i] /= __total_area;
  131. // __densities_ contains normalized areas
  132. __areas_.assign(__densities_.size(), result_type());
  133. std::partial_sum(__densities_.begin(), __densities_.end() - 1, __areas_.begin() + 1);
  134. // __areas_ contains partial sums of normalized areas: [0, __densities_ - 1]
  135. __densities_.back() = 1 - __areas_.back(); // correct round off error
  136. for (size_t __i = 0; __i < __densities_.size(); ++__i)
  137. __densities_[__i] /= (__b_[__i + 1] - __b_[__i]);
  138. // __densities_ now contains __densities_
  139. }
  140. template <class _RealType>
  141. piecewise_constant_distribution<_RealType>::param_type::param_type() : __b_(2), __densities_(1, 1.0), __areas_(1, 0.0) {
  142. __b_[1] = 1;
  143. }
  144. template <class _RealType>
  145. template <class _InputIteratorB, class _InputIteratorW>
  146. piecewise_constant_distribution<_RealType>::param_type::param_type(
  147. _InputIteratorB __f_b, _InputIteratorB __l_b, _InputIteratorW __f_w)
  148. : __b_(__f_b, __l_b) {
  149. if (__b_.size() < 2) {
  150. __b_.resize(2);
  151. __b_[0] = 0;
  152. __b_[1] = 1;
  153. __densities_.assign(1, 1.0);
  154. __areas_.assign(1, 0.0);
  155. } else {
  156. __densities_.reserve(__b_.size() - 1);
  157. for (size_t __i = 0; __i < __b_.size() - 1; ++__i, ++__f_w)
  158. __densities_.push_back(*__f_w);
  159. __init();
  160. }
  161. }
  162. #ifndef _LIBCPP_CXX03_LANG
  163. template <class _RealType>
  164. template <class _UnaryOperation>
  165. piecewise_constant_distribution<_RealType>::param_type::param_type(
  166. initializer_list<result_type> __bl, _UnaryOperation __fw)
  167. : __b_(__bl.begin(), __bl.end()) {
  168. if (__b_.size() < 2) {
  169. __b_.resize(2);
  170. __b_[0] = 0;
  171. __b_[1] = 1;
  172. __densities_.assign(1, 1.0);
  173. __areas_.assign(1, 0.0);
  174. } else {
  175. __densities_.reserve(__b_.size() - 1);
  176. for (size_t __i = 0; __i < __b_.size() - 1; ++__i)
  177. __densities_.push_back(__fw((__b_[__i + 1] + __b_[__i]) * .5));
  178. __init();
  179. }
  180. }
  181. #endif // _LIBCPP_CXX03_LANG
  182. template <class _RealType>
  183. template <class _UnaryOperation>
  184. piecewise_constant_distribution<_RealType>::param_type::param_type(
  185. size_t __nw, result_type __xmin, result_type __xmax, _UnaryOperation __fw)
  186. : __b_(__nw == 0 ? 2 : __nw + 1) {
  187. size_t __n = __b_.size() - 1;
  188. result_type __d = (__xmax - __xmin) / __n;
  189. __densities_.reserve(__n);
  190. for (size_t __i = 0; __i < __n; ++__i) {
  191. __b_[__i] = __xmin + __i * __d;
  192. __densities_.push_back(__fw(__b_[__i] + __d * .5));
  193. }
  194. __b_[__n] = __xmax;
  195. __init();
  196. }
  197. template <class _RealType>
  198. template <class _URNG>
  199. _RealType piecewise_constant_distribution<_RealType>::operator()(_URNG& __g, const param_type& __p) {
  200. static_assert(__libcpp_random_is_valid_urng<_URNG>::value, "");
  201. typedef uniform_real_distribution<result_type> _Gen;
  202. result_type __u = _Gen()(__g);
  203. ptrdiff_t __k = std::upper_bound(__p.__areas_.begin(), __p.__areas_.end(), __u) - __p.__areas_.begin() - 1;
  204. return (__u - __p.__areas_[__k]) / __p.__densities_[__k] + __p.__b_[__k];
  205. }
  206. template <class _CharT, class _Traits, class _RT>
  207. _LIBCPP_HIDE_FROM_ABI basic_ostream<_CharT, _Traits>&
  208. operator<<(basic_ostream<_CharT, _Traits>& __os, const piecewise_constant_distribution<_RT>& __x) {
  209. __save_flags<_CharT, _Traits> __lx(__os);
  210. typedef basic_ostream<_CharT, _Traits> _OStream;
  211. __os.flags(_OStream::dec | _OStream::left | _OStream::fixed | _OStream::scientific);
  212. _CharT __sp = __os.widen(' ');
  213. __os.fill(__sp);
  214. size_t __n = __x.__p_.__b_.size();
  215. __os << __n;
  216. for (size_t __i = 0; __i < __n; ++__i)
  217. __os << __sp << __x.__p_.__b_[__i];
  218. __n = __x.__p_.__densities_.size();
  219. __os << __sp << __n;
  220. for (size_t __i = 0; __i < __n; ++__i)
  221. __os << __sp << __x.__p_.__densities_[__i];
  222. __n = __x.__p_.__areas_.size();
  223. __os << __sp << __n;
  224. for (size_t __i = 0; __i < __n; ++__i)
  225. __os << __sp << __x.__p_.__areas_[__i];
  226. return __os;
  227. }
  228. template <class _CharT, class _Traits, class _RT>
  229. _LIBCPP_HIDE_FROM_ABI basic_istream<_CharT, _Traits>&
  230. operator>>(basic_istream<_CharT, _Traits>& __is, piecewise_constant_distribution<_RT>& __x) {
  231. typedef piecewise_constant_distribution<_RT> _Eng;
  232. typedef typename _Eng::result_type result_type;
  233. __save_flags<_CharT, _Traits> __lx(__is);
  234. typedef basic_istream<_CharT, _Traits> _Istream;
  235. __is.flags(_Istream::dec | _Istream::skipws);
  236. size_t __n;
  237. __is >> __n;
  238. vector<result_type> __b(__n);
  239. for (size_t __i = 0; __i < __n; ++__i)
  240. __is >> __b[__i];
  241. __is >> __n;
  242. vector<result_type> __densities(__n);
  243. for (size_t __i = 0; __i < __n; ++__i)
  244. __is >> __densities[__i];
  245. __is >> __n;
  246. vector<result_type> __areas(__n);
  247. for (size_t __i = 0; __i < __n; ++__i)
  248. __is >> __areas[__i];
  249. if (!__is.fail()) {
  250. swap(__x.__p_.__b_, __b);
  251. swap(__x.__p_.__densities_, __densities);
  252. swap(__x.__p_.__areas_, __areas);
  253. }
  254. return __is;
  255. }
  256. _LIBCPP_END_NAMESPACE_STD
  257. _LIBCPP_POP_MACROS
  258. #endif // _LIBCPP___RANDOM_PIECEWISE_CONSTANT_DISTRIBUTION_H