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 29 template<class _RealType = double> 30 class _LIBCPP_TEMPLATE_VIS piecewise_constant_distribution 31 { 32 public: 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 88 private: 89 param_type __p_; 90 91 public: 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 168 template<class _RealType> 169 typename piecewise_constant_distribution<_RealType>::param_type & 170 piecewise_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 185 template<class _RealType> 186 void 187 piecewise_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 206 template<class _RealType> 207 piecewise_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 215 template<class _RealType> 216 template<class _InputIteratorB, class _InputIteratorW> 217 piecewise_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 240 template<class _RealType> 241 template<class _UnaryOperation> 242 piecewise_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 265 template<class _RealType> 266 template<class _UnaryOperation> 267 piecewise_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 283 template<class _RealType> 284 template<class _URNG> 285 _RealType 286 piecewise_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 296 template <class _CharT, class _Traits, class _RT> 297 basic_ostream<_CharT, _Traits>& 298 operator<<(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 322 template <class _CharT, class _Traits, class _RT> 323 basic_istream<_CharT, _Traits>& 324 operator>>(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