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