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