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