piecewise_linear_distribution.h 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375
  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 <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_linear_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_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,
  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_linear_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_linear_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_linear_distribution<_RT>& __x);
  72. };
  73. private:
  74. param_type __p_;
  75. public:
  76. // constructor and reset functions
  77. _LIBCPP_INLINE_VISIBILITY
  78. piecewise_linear_distribution() {}
  79. template<class _InputIteratorB, class _InputIteratorW>
  80. _LIBCPP_INLINE_VISIBILITY
  81. piecewise_linear_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_linear_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_linear_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_linear_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_linear_distribution& __x,
  124. const piecewise_linear_distribution& __y)
  125. {return __x.__p_ == __y.__p_;}
  126. friend _LIBCPP_INLINE_VISIBILITY
  127. bool operator!=(const piecewise_linear_distribution& __x,
  128. const piecewise_linear_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_linear_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_linear_distribution<_RT>& __x);
  140. };
  141. template<class _RealType>
  142. typename piecewise_linear_distribution<_RealType>::param_type &
  143. piecewise_linear_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_linear_distribution<_RealType>::param_type::__init()
  159. {
  160. __areas_.assign(__densities_.size() - 1, result_type());
  161. result_type __sp = 0;
  162. for (size_t __i = 0; __i < __areas_.size(); ++__i)
  163. {
  164. __areas_[__i] = (__densities_[__i+1] + __densities_[__i]) *
  165. (__b_[__i+1] - __b_[__i]) * .5;
  166. __sp += __areas_[__i];
  167. }
  168. for (size_t __i = __areas_.size(); __i > 1;)
  169. {
  170. --__i;
  171. __areas_[__i] = __areas_[__i-1] / __sp;
  172. }
  173. __areas_[0] = 0;
  174. for (size_t __i = 1; __i < __areas_.size(); ++__i)
  175. __areas_[__i] += __areas_[__i-1];
  176. for (size_t __i = 0; __i < __densities_.size(); ++__i)
  177. __densities_[__i] /= __sp;
  178. }
  179. template<class _RealType>
  180. piecewise_linear_distribution<_RealType>::param_type::param_type()
  181. : __b_(2),
  182. __densities_(2, 1.0),
  183. __areas_(1, 0.0)
  184. {
  185. __b_[1] = 1;
  186. }
  187. template<class _RealType>
  188. template<class _InputIteratorB, class _InputIteratorW>
  189. piecewise_linear_distribution<_RealType>::param_type::param_type(
  190. _InputIteratorB __f_b, _InputIteratorB __l_b, _InputIteratorW __f_w)
  191. : __b_(__f_b, __l_b)
  192. {
  193. if (__b_.size() < 2)
  194. {
  195. __b_.resize(2);
  196. __b_[0] = 0;
  197. __b_[1] = 1;
  198. __densities_.assign(2, 1.0);
  199. __areas_.assign(1, 0.0);
  200. }
  201. else
  202. {
  203. __densities_.reserve(__b_.size());
  204. for (size_t __i = 0; __i < __b_.size(); ++__i, ++__f_w)
  205. __densities_.push_back(*__f_w);
  206. __init();
  207. }
  208. }
  209. #ifndef _LIBCPP_CXX03_LANG
  210. template<class _RealType>
  211. template<class _UnaryOperation>
  212. piecewise_linear_distribution<_RealType>::param_type::param_type(
  213. initializer_list<result_type> __bl, _UnaryOperation __fw)
  214. : __b_(__bl.begin(), __bl.end())
  215. {
  216. if (__b_.size() < 2)
  217. {
  218. __b_.resize(2);
  219. __b_[0] = 0;
  220. __b_[1] = 1;
  221. __densities_.assign(2, 1.0);
  222. __areas_.assign(1, 0.0);
  223. }
  224. else
  225. {
  226. __densities_.reserve(__b_.size());
  227. for (size_t __i = 0; __i < __b_.size(); ++__i)
  228. __densities_.push_back(__fw(__b_[__i]));
  229. __init();
  230. }
  231. }
  232. #endif // _LIBCPP_CXX03_LANG
  233. template<class _RealType>
  234. template<class _UnaryOperation>
  235. piecewise_linear_distribution<_RealType>::param_type::param_type(
  236. size_t __nw, result_type __xmin, result_type __xmax, _UnaryOperation __fw)
  237. : __b_(__nw == 0 ? 2 : __nw + 1)
  238. {
  239. size_t __n = __b_.size() - 1;
  240. result_type __d = (__xmax - __xmin) / __n;
  241. __densities_.reserve(__b_.size());
  242. for (size_t __i = 0; __i < __n; ++__i)
  243. {
  244. __b_[__i] = __xmin + __i * __d;
  245. __densities_.push_back(__fw(__b_[__i]));
  246. }
  247. __b_[__n] = __xmax;
  248. __densities_.push_back(__fw(__b_[__n]));
  249. __init();
  250. }
  251. template<class _RealType>
  252. template<class _URNG>
  253. _RealType
  254. piecewise_linear_distribution<_RealType>::operator()(_URNG& __g, const param_type& __p)
  255. {
  256. static_assert(__libcpp_random_is_valid_urng<_URNG>::value, "");
  257. typedef uniform_real_distribution<result_type> _Gen;
  258. result_type __u = _Gen()(__g);
  259. ptrdiff_t __k = _VSTD::upper_bound(__p.__areas_.begin(), __p.__areas_.end(),
  260. __u) - __p.__areas_.begin() - 1;
  261. __u -= __p.__areas_[__k];
  262. const result_type __dk = __p.__densities_[__k];
  263. const result_type __dk1 = __p.__densities_[__k+1];
  264. const result_type __deltad = __dk1 - __dk;
  265. const result_type __bk = __p.__b_[__k];
  266. if (__deltad == 0)
  267. return __u / __dk + __bk;
  268. const result_type __bk1 = __p.__b_[__k+1];
  269. const result_type __deltab = __bk1 - __bk;
  270. return (__bk * __dk1 - __bk1 * __dk +
  271. _VSTD::sqrt(__deltab * (__deltab * __dk * __dk + 2 * __deltad * __u))) /
  272. __deltad;
  273. }
  274. template <class _CharT, class _Traits, class _RT>
  275. _LIBCPP_HIDE_FROM_ABI basic_ostream<_CharT, _Traits>&
  276. operator<<(basic_ostream<_CharT, _Traits>& __os,
  277. const piecewise_linear_distribution<_RT>& __x)
  278. {
  279. __save_flags<_CharT, _Traits> __lx(__os);
  280. typedef basic_ostream<_CharT, _Traits> _OStream;
  281. __os.flags(_OStream::dec | _OStream::left | _OStream::fixed |
  282. _OStream::scientific);
  283. _CharT __sp = __os.widen(' ');
  284. __os.fill(__sp);
  285. size_t __n = __x.__p_.__b_.size();
  286. __os << __n;
  287. for (size_t __i = 0; __i < __n; ++__i)
  288. __os << __sp << __x.__p_.__b_[__i];
  289. __n = __x.__p_.__densities_.size();
  290. __os << __sp << __n;
  291. for (size_t __i = 0; __i < __n; ++__i)
  292. __os << __sp << __x.__p_.__densities_[__i];
  293. __n = __x.__p_.__areas_.size();
  294. __os << __sp << __n;
  295. for (size_t __i = 0; __i < __n; ++__i)
  296. __os << __sp << __x.__p_.__areas_[__i];
  297. return __os;
  298. }
  299. template <class _CharT, class _Traits, class _RT>
  300. _LIBCPP_HIDE_FROM_ABI basic_istream<_CharT, _Traits>&
  301. operator>>(basic_istream<_CharT, _Traits>& __is,
  302. piecewise_linear_distribution<_RT>& __x)
  303. {
  304. typedef piecewise_linear_distribution<_RT> _Eng;
  305. typedef typename _Eng::result_type result_type;
  306. __save_flags<_CharT, _Traits> __lx(__is);
  307. typedef basic_istream<_CharT, _Traits> _Istream;
  308. __is.flags(_Istream::dec | _Istream::skipws);
  309. size_t __n;
  310. __is >> __n;
  311. vector<result_type> __b(__n);
  312. for (size_t __i = 0; __i < __n; ++__i)
  313. __is >> __b[__i];
  314. __is >> __n;
  315. vector<result_type> __densities(__n);
  316. for (size_t __i = 0; __i < __n; ++__i)
  317. __is >> __densities[__i];
  318. __is >> __n;
  319. vector<result_type> __areas(__n);
  320. for (size_t __i = 0; __i < __n; ++__i)
  321. __is >> __areas[__i];
  322. if (!__is.fail())
  323. {
  324. swap(__x.__p_.__b_, __b);
  325. swap(__x.__p_.__densities_, __densities);
  326. swap(__x.__p_.__areas_, __areas);
  327. }
  328. return __is;
  329. }
  330. _LIBCPP_END_NAMESPACE_STD
  331. _LIBCPP_POP_MACROS
  332. #endif // _LIBCPP___RANDOM_PIECEWISE_LINEAR_DISTRIBUTION_H