piecewise_linear_distribution.h 12 KB

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