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_CONSTANT_DISTRIBUTION_H
10 #define _LIBCPP___RANDOM_PIECEWISE_CONSTANT_DISTRIBUTION_H
11 
12 #include <__algorithm/upper_bound.h>
13 #include <__config>
14 #include <__random/uniform_real_distribution.h>
15 #include <iosfwd>
16 #include <numeric>
17 #include <vector>
18 
19 #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
20 #pragma GCC system_header
21 #endif
22 
23 _LIBCPP_PUSH_MACROS
24 #include <__undef_macros>
25 
26 _LIBCPP_BEGIN_NAMESPACE_STD
27 
28 template<class _RealType = double>
29 class _LIBCPP_TEMPLATE_VIS piecewise_constant_distribution
30 {
31 public:
32     // types
33     typedef _RealType result_type;
34 
35     class _LIBCPP_TEMPLATE_VIS param_type
36     {
37         vector<result_type> __b_;
38         vector<result_type> __densities_;
39         vector<result_type> __areas_;
40     public:
41         typedef piecewise_constant_distribution distribution_type;
42 
43         param_type();
44         template<class _InputIteratorB, class _InputIteratorW>
45             param_type(_InputIteratorB __fB, _InputIteratorB __lB,
46                        _InputIteratorW __fW);
47 #ifndef _LIBCPP_CXX03_LANG
48         template<class _UnaryOperation>
49             param_type(initializer_list<result_type> __bl, _UnaryOperation __fw);
50 #endif // _LIBCPP_CXX03_LANG
51         template<class _UnaryOperation>
52             param_type(size_t __nw, result_type __xmin, result_type __xmax,
53                        _UnaryOperation __fw);
54         param_type(param_type const&) = default;
55         param_type & operator=(const param_type& __rhs);
56 
57         _LIBCPP_INLINE_VISIBILITY
58         vector<result_type> intervals() const {return __b_;}
59         _LIBCPP_INLINE_VISIBILITY
60         vector<result_type> densities() const {return __densities_;}
61 
62         friend _LIBCPP_INLINE_VISIBILITY
63             bool operator==(const param_type& __x, const param_type& __y)
64             {return __x.__densities_ == __y.__densities_ && __x.__b_ == __y.__b_;}
65         friend _LIBCPP_INLINE_VISIBILITY
66             bool operator!=(const param_type& __x, const param_type& __y)
67             {return !(__x == __y);}
68 
69     private:
70         void __init();
71 
72         friend class piecewise_constant_distribution;
73 
74         template <class _CharT, class _Traits, class _RT>
75         friend
76         basic_ostream<_CharT, _Traits>&
77         operator<<(basic_ostream<_CharT, _Traits>& __os,
78                    const piecewise_constant_distribution<_RT>& __x);
79 
80         template <class _CharT, class _Traits, class _RT>
81         friend
82         basic_istream<_CharT, _Traits>&
83         operator>>(basic_istream<_CharT, _Traits>& __is,
84                    piecewise_constant_distribution<_RT>& __x);
85     };
86 
87 private:
88     param_type __p_;
89 
90 public:
91     // constructor and reset functions
92     _LIBCPP_INLINE_VISIBILITY
93     piecewise_constant_distribution() {}
94     template<class _InputIteratorB, class _InputIteratorW>
95         _LIBCPP_INLINE_VISIBILITY
96         piecewise_constant_distribution(_InputIteratorB __fB,
97                                         _InputIteratorB __lB,
98                                         _InputIteratorW __fW)
99         : __p_(__fB, __lB, __fW) {}
100 
101 #ifndef _LIBCPP_CXX03_LANG
102     template<class _UnaryOperation>
103         _LIBCPP_INLINE_VISIBILITY
104         piecewise_constant_distribution(initializer_list<result_type> __bl,
105                                         _UnaryOperation __fw)
106         : __p_(__bl, __fw) {}
107 #endif // _LIBCPP_CXX03_LANG
108 
109     template<class _UnaryOperation>
110         _LIBCPP_INLINE_VISIBILITY
111         piecewise_constant_distribution(size_t __nw, result_type __xmin,
112                                         result_type __xmax, _UnaryOperation __fw)
113         : __p_(__nw, __xmin, __xmax, __fw) {}
114 
115     _LIBCPP_INLINE_VISIBILITY
116     explicit piecewise_constant_distribution(const param_type& __p)
117         : __p_(__p) {}
118 
119     _LIBCPP_INLINE_VISIBILITY
120     void reset() {}
121 
122     // generating functions
123     template<class _URNG>
124         _LIBCPP_INLINE_VISIBILITY
125         result_type operator()(_URNG& __g)
126         {return (*this)(__g, __p_);}
127     template<class _URNG> result_type operator()(_URNG& __g, const param_type& __p);
128 
129     // property functions
130     _LIBCPP_INLINE_VISIBILITY
131     vector<result_type> intervals() const {return __p_.intervals();}
132     _LIBCPP_INLINE_VISIBILITY
133     vector<result_type> densities() const {return __p_.densities();}
134 
135     _LIBCPP_INLINE_VISIBILITY
136     param_type param() const {return __p_;}
137     _LIBCPP_INLINE_VISIBILITY
138     void param(const param_type& __p) {__p_ = __p;}
139 
140     _LIBCPP_INLINE_VISIBILITY
141     result_type min() const {return __p_.__b_.front();}
142     _LIBCPP_INLINE_VISIBILITY
143     result_type max() const {return __p_.__b_.back();}
144 
145     friend _LIBCPP_INLINE_VISIBILITY
146         bool operator==(const piecewise_constant_distribution& __x,
147                         const piecewise_constant_distribution& __y)
148         {return __x.__p_ == __y.__p_;}
149     friend _LIBCPP_INLINE_VISIBILITY
150         bool operator!=(const piecewise_constant_distribution& __x,
151                            const piecewise_constant_distribution& __y)
152         {return !(__x == __y);}
153 
154     template <class _CharT, class _Traits, class _RT>
155     friend
156     basic_ostream<_CharT, _Traits>&
157     operator<<(basic_ostream<_CharT, _Traits>& __os,
158                const piecewise_constant_distribution<_RT>& __x);
159 
160     template <class _CharT, class _Traits, class _RT>
161     friend
162     basic_istream<_CharT, _Traits>&
163     operator>>(basic_istream<_CharT, _Traits>& __is,
164                piecewise_constant_distribution<_RT>& __x);
165 };
166 
167 template<class _RealType>
168 typename piecewise_constant_distribution<_RealType>::param_type &
169 piecewise_constant_distribution<_RealType>::param_type::operator=
170                                                        (const param_type& __rhs)
171 {
172 //  These can throw
173     __b_.reserve        (__rhs.__b_.size ());
174     __densities_.reserve(__rhs.__densities_.size());
175     __areas_.reserve    (__rhs.__areas_.size());
176 
177 //  These can not throw
178     __b_         = __rhs.__b_;
179     __densities_ = __rhs.__densities_;
180     __areas_     =  __rhs.__areas_;
181     return *this;
182 }
183 
184 template<class _RealType>
185 void
186 piecewise_constant_distribution<_RealType>::param_type::__init()
187 {
188     // __densities_ contains non-normalized areas
189     result_type __total_area = _VSTD::accumulate(__densities_.begin(),
190                                                 __densities_.end(),
191                                                 result_type());
192     for (size_t __i = 0; __i < __densities_.size(); ++__i)
193         __densities_[__i] /= __total_area;
194     // __densities_ contains normalized areas
195     __areas_.assign(__densities_.size(), result_type());
196     _VSTD::partial_sum(__densities_.begin(), __densities_.end() - 1,
197                                                           __areas_.begin() + 1);
198     // __areas_ contains partial sums of normalized areas: [0, __densities_ - 1]
199     __densities_.back() = 1 - __areas_.back();  // correct round off error
200     for (size_t __i = 0; __i < __densities_.size(); ++__i)
201         __densities_[__i] /= (__b_[__i+1] - __b_[__i]);
202     // __densities_ now contains __densities_
203 }
204 
205 template<class _RealType>
206 piecewise_constant_distribution<_RealType>::param_type::param_type()
207     : __b_(2),
208       __densities_(1, 1.0),
209       __areas_(1, 0.0)
210 {
211     __b_[1] = 1;
212 }
213 
214 template<class _RealType>
215 template<class _InputIteratorB, class _InputIteratorW>
216 piecewise_constant_distribution<_RealType>::param_type::param_type(
217         _InputIteratorB __fB, _InputIteratorB __lB, _InputIteratorW __fW)
218     : __b_(__fB, __lB)
219 {
220     if (__b_.size() < 2)
221     {
222         __b_.resize(2);
223         __b_[0] = 0;
224         __b_[1] = 1;
225         __densities_.assign(1, 1.0);
226         __areas_.assign(1, 0.0);
227     }
228     else
229     {
230         __densities_.reserve(__b_.size() - 1);
231         for (size_t __i = 0; __i < __b_.size() - 1; ++__i, ++__fW)
232             __densities_.push_back(*__fW);
233         __init();
234     }
235 }
236 
237 #ifndef _LIBCPP_CXX03_LANG
238 
239 template<class _RealType>
240 template<class _UnaryOperation>
241 piecewise_constant_distribution<_RealType>::param_type::param_type(
242         initializer_list<result_type> __bl, _UnaryOperation __fw)
243     : __b_(__bl.begin(), __bl.end())
244 {
245     if (__b_.size() < 2)
246     {
247         __b_.resize(2);
248         __b_[0] = 0;
249         __b_[1] = 1;
250         __densities_.assign(1, 1.0);
251         __areas_.assign(1, 0.0);
252     }
253     else
254     {
255         __densities_.reserve(__b_.size() - 1);
256         for (size_t __i = 0; __i < __b_.size() - 1; ++__i)
257             __densities_.push_back(__fw((__b_[__i+1] + __b_[__i])*.5));
258         __init();
259     }
260 }
261 
262 #endif // _LIBCPP_CXX03_LANG
263 
264 template<class _RealType>
265 template<class _UnaryOperation>
266 piecewise_constant_distribution<_RealType>::param_type::param_type(
267         size_t __nw, result_type __xmin, result_type __xmax, _UnaryOperation __fw)
268     : __b_(__nw == 0 ? 2 : __nw + 1)
269 {
270     size_t __n = __b_.size() - 1;
271     result_type __d = (__xmax - __xmin) / __n;
272     __densities_.reserve(__n);
273     for (size_t __i = 0; __i < __n; ++__i)
274     {
275         __b_[__i] = __xmin + __i * __d;
276         __densities_.push_back(__fw(__b_[__i] + __d*.5));
277     }
278     __b_[__n] = __xmax;
279     __init();
280 }
281 
282 template<class _RealType>
283 template<class _URNG>
284 _RealType
285 piecewise_constant_distribution<_RealType>::operator()(_URNG& __g, const param_type& __p)
286 {
287     typedef uniform_real_distribution<result_type> _Gen;
288     result_type __u = _Gen()(__g);
289     ptrdiff_t __k = _VSTD::upper_bound(__p.__areas_.begin(), __p.__areas_.end(),
290                                       __u) - __p.__areas_.begin() - 1;
291     return (__u - __p.__areas_[__k]) / __p.__densities_[__k] + __p.__b_[__k];
292 }
293 
294 template <class _CharT, class _Traits, class _RT>
295 basic_ostream<_CharT, _Traits>&
296 operator<<(basic_ostream<_CharT, _Traits>& __os,
297            const piecewise_constant_distribution<_RT>& __x)
298 {
299     __save_flags<_CharT, _Traits> __lx(__os);
300     typedef basic_ostream<_CharT, _Traits> _OStream;
301     __os.flags(_OStream::dec | _OStream::left | _OStream::fixed |
302                _OStream::scientific);
303     _CharT __sp = __os.widen(' ');
304     __os.fill(__sp);
305     size_t __n = __x.__p_.__b_.size();
306     __os << __n;
307     for (size_t __i = 0; __i < __n; ++__i)
308         __os << __sp << __x.__p_.__b_[__i];
309     __n = __x.__p_.__densities_.size();
310     __os << __sp << __n;
311     for (size_t __i = 0; __i < __n; ++__i)
312         __os << __sp << __x.__p_.__densities_[__i];
313     __n = __x.__p_.__areas_.size();
314     __os << __sp << __n;
315     for (size_t __i = 0; __i < __n; ++__i)
316         __os << __sp << __x.__p_.__areas_[__i];
317     return __os;
318 }
319 
320 template <class _CharT, class _Traits, class _RT>
321 basic_istream<_CharT, _Traits>&
322 operator>>(basic_istream<_CharT, _Traits>& __is,
323            piecewise_constant_distribution<_RT>& __x)
324 {
325     typedef piecewise_constant_distribution<_RT> _Eng;
326     typedef typename _Eng::result_type result_type;
327     __save_flags<_CharT, _Traits> __lx(__is);
328     typedef basic_istream<_CharT, _Traits> _Istream;
329     __is.flags(_Istream::dec | _Istream::skipws);
330     size_t __n;
331     __is >> __n;
332     vector<result_type> __b(__n);
333     for (size_t __i = 0; __i < __n; ++__i)
334         __is >> __b[__i];
335     __is >> __n;
336     vector<result_type> __densities(__n);
337     for (size_t __i = 0; __i < __n; ++__i)
338         __is >> __densities[__i];
339     __is >> __n;
340     vector<result_type> __areas(__n);
341     for (size_t __i = 0; __i < __n; ++__i)
342         __is >> __areas[__i];
343     if (!__is.fail())
344     {
345         swap(__x.__p_.__b_, __b);
346         swap(__x.__p_.__densities_, __densities);
347         swap(__x.__p_.__areas_, __areas);
348     }
349     return __is;
350 }
351 
352 _LIBCPP_END_NAMESPACE_STD
353 
354 _LIBCPP_POP_MACROS
355 
356 #endif // _LIBCPP___RANDOM_PIECEWISE_CONSTANT_DISTRIBUTION_H
357