piecewise_linear_distribution.h 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315
  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_LINEAR_DISTRIBUTION_H
  9. #define _LIBCPP___RANDOM_PIECEWISE_LINEAR_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 <cmath>
  15. #include <iosfwd>
  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_linear_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_linear_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_linear_distribution;
  56. template <class _CharT, class _Traits, class _RT>
  57. friend basic_ostream<_CharT, _Traits>&
  58. operator<<(basic_ostream<_CharT, _Traits>& __os, const piecewise_linear_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_linear_distribution<_RT>& __x);
  62. };
  63. private:
  64. param_type __p_;
  65. public:
  66. // constructor and reset functions
  67. _LIBCPP_HIDE_FROM_ABI piecewise_linear_distribution() {}
  68. template <class _InputIteratorB, class _InputIteratorW>
  69. _LIBCPP_HIDE_FROM_ABI
  70. piecewise_linear_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_linear_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_linear_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_linear_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_linear_distribution& __x, const piecewise_linear_distribution& __y) {
  99. return __x.__p_ == __y.__p_;
  100. }
  101. friend _LIBCPP_HIDE_FROM_ABI bool
  102. operator!=(const piecewise_linear_distribution& __x, const piecewise_linear_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_linear_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_linear_distribution<_RT>& __x);
  111. };
  112. template <class _RealType>
  113. typename piecewise_linear_distribution<_RealType>::param_type&
  114. piecewise_linear_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_linear_distribution<_RealType>::param_type::__init() {
  127. __areas_.assign(__densities_.size() - 1, result_type());
  128. result_type __sp = 0;
  129. for (size_t __i = 0; __i < __areas_.size(); ++__i) {
  130. __areas_[__i] = (__densities_[__i + 1] + __densities_[__i]) * (__b_[__i + 1] - __b_[__i]) * .5;
  131. __sp += __areas_[__i];
  132. }
  133. for (size_t __i = __areas_.size(); __i > 1;) {
  134. --__i;
  135. __areas_[__i] = __areas_[__i - 1] / __sp;
  136. }
  137. __areas_[0] = 0;
  138. for (size_t __i = 1; __i < __areas_.size(); ++__i)
  139. __areas_[__i] += __areas_[__i - 1];
  140. for (size_t __i = 0; __i < __densities_.size(); ++__i)
  141. __densities_[__i] /= __sp;
  142. }
  143. template <class _RealType>
  144. piecewise_linear_distribution<_RealType>::param_type::param_type() : __b_(2), __densities_(2, 1.0), __areas_(1, 0.0) {
  145. __b_[1] = 1;
  146. }
  147. template <class _RealType>
  148. template <class _InputIteratorB, class _InputIteratorW>
  149. piecewise_linear_distribution<_RealType>::param_type::param_type(
  150. _InputIteratorB __f_b, _InputIteratorB __l_b, _InputIteratorW __f_w)
  151. : __b_(__f_b, __l_b) {
  152. if (__b_.size() < 2) {
  153. __b_.resize(2);
  154. __b_[0] = 0;
  155. __b_[1] = 1;
  156. __densities_.assign(2, 1.0);
  157. __areas_.assign(1, 0.0);
  158. } else {
  159. __densities_.reserve(__b_.size());
  160. for (size_t __i = 0; __i < __b_.size(); ++__i, ++__f_w)
  161. __densities_.push_back(*__f_w);
  162. __init();
  163. }
  164. }
  165. #ifndef _LIBCPP_CXX03_LANG
  166. template <class _RealType>
  167. template <class _UnaryOperation>
  168. piecewise_linear_distribution<_RealType>::param_type::param_type(
  169. initializer_list<result_type> __bl, _UnaryOperation __fw)
  170. : __b_(__bl.begin(), __bl.end()) {
  171. if (__b_.size() < 2) {
  172. __b_.resize(2);
  173. __b_[0] = 0;
  174. __b_[1] = 1;
  175. __densities_.assign(2, 1.0);
  176. __areas_.assign(1, 0.0);
  177. } else {
  178. __densities_.reserve(__b_.size());
  179. for (size_t __i = 0; __i < __b_.size(); ++__i)
  180. __densities_.push_back(__fw(__b_[__i]));
  181. __init();
  182. }
  183. }
  184. #endif // _LIBCPP_CXX03_LANG
  185. template <class _RealType>
  186. template <class _UnaryOperation>
  187. piecewise_linear_distribution<_RealType>::param_type::param_type(
  188. size_t __nw, result_type __xmin, result_type __xmax, _UnaryOperation __fw)
  189. : __b_(__nw == 0 ? 2 : __nw + 1) {
  190. size_t __n = __b_.size() - 1;
  191. result_type __d = (__xmax - __xmin) / __n;
  192. __densities_.reserve(__b_.size());
  193. for (size_t __i = 0; __i < __n; ++__i) {
  194. __b_[__i] = __xmin + __i * __d;
  195. __densities_.push_back(__fw(__b_[__i]));
  196. }
  197. __b_[__n] = __xmax;
  198. __densities_.push_back(__fw(__b_[__n]));
  199. __init();
  200. }
  201. template <class _RealType>
  202. template <class _URNG>
  203. _RealType piecewise_linear_distribution<_RealType>::operator()(_URNG& __g, const param_type& __p) {
  204. static_assert(__libcpp_random_is_valid_urng<_URNG>::value, "");
  205. typedef uniform_real_distribution<result_type> _Gen;
  206. result_type __u = _Gen()(__g);
  207. ptrdiff_t __k = std::upper_bound(__p.__areas_.begin(), __p.__areas_.end(), __u) - __p.__areas_.begin() - 1;
  208. __u -= __p.__areas_[__k];
  209. const result_type __dk = __p.__densities_[__k];
  210. const result_type __dk1 = __p.__densities_[__k + 1];
  211. const result_type __deltad = __dk1 - __dk;
  212. const result_type __bk = __p.__b_[__k];
  213. if (__deltad == 0)
  214. return __u / __dk + __bk;
  215. const result_type __bk1 = __p.__b_[__k + 1];
  216. const result_type __deltab = __bk1 - __bk;
  217. return (__bk * __dk1 - __bk1 * __dk + std::sqrt(__deltab * (__deltab * __dk * __dk + 2 * __deltad * __u))) / __deltad;
  218. }
  219. template <class _CharT, class _Traits, class _RT>
  220. _LIBCPP_HIDE_FROM_ABI basic_ostream<_CharT, _Traits>&
  221. operator<<(basic_ostream<_CharT, _Traits>& __os, const piecewise_linear_distribution<_RT>& __x) {
  222. __save_flags<_CharT, _Traits> __lx(__os);
  223. typedef basic_ostream<_CharT, _Traits> _OStream;
  224. __os.flags(_OStream::dec | _OStream::left | _OStream::fixed | _OStream::scientific);
  225. _CharT __sp = __os.widen(' ');
  226. __os.fill(__sp);
  227. size_t __n = __x.__p_.__b_.size();
  228. __os << __n;
  229. for (size_t __i = 0; __i < __n; ++__i)
  230. __os << __sp << __x.__p_.__b_[__i];
  231. __n = __x.__p_.__densities_.size();
  232. __os << __sp << __n;
  233. for (size_t __i = 0; __i < __n; ++__i)
  234. __os << __sp << __x.__p_.__densities_[__i];
  235. __n = __x.__p_.__areas_.size();
  236. __os << __sp << __n;
  237. for (size_t __i = 0; __i < __n; ++__i)
  238. __os << __sp << __x.__p_.__areas_[__i];
  239. return __os;
  240. }
  241. template <class _CharT, class _Traits, class _RT>
  242. _LIBCPP_HIDE_FROM_ABI basic_istream<_CharT, _Traits>&
  243. operator>>(basic_istream<_CharT, _Traits>& __is, piecewise_linear_distribution<_RT>& __x) {
  244. typedef piecewise_linear_distribution<_RT> _Eng;
  245. typedef typename _Eng::result_type result_type;
  246. __save_flags<_CharT, _Traits> __lx(__is);
  247. typedef basic_istream<_CharT, _Traits> _Istream;
  248. __is.flags(_Istream::dec | _Istream::skipws);
  249. size_t __n;
  250. __is >> __n;
  251. vector<result_type> __b(__n);
  252. for (size_t __i = 0; __i < __n; ++__i)
  253. __is >> __b[__i];
  254. __is >> __n;
  255. vector<result_type> __densities(__n);
  256. for (size_t __i = 0; __i < __n; ++__i)
  257. __is >> __densities[__i];
  258. __is >> __n;
  259. vector<result_type> __areas(__n);
  260. for (size_t __i = 0; __i < __n; ++__i)
  261. __is >> __areas[__i];
  262. if (!__is.fail()) {
  263. swap(__x.__p_.__b_, __b);
  264. swap(__x.__p_.__densities_, __densities);
  265. swap(__x.__p_.__areas_, __areas);
  266. }
  267. return __is;
  268. }
  269. _LIBCPP_END_NAMESPACE_STD
  270. _LIBCPP_POP_MACROS
  271. #endif // _LIBCPP___RANDOM_PIECEWISE_LINEAR_DISTRIBUTION_H