xtensor
Loading...
Searching...
No Matches
xrandom.hpp
1/***************************************************************************
2 * Copyright (c) Johan Mabille, Sylvain Corlay and Wolf Vollprecht *
3 * Copyright (c) QuantStack *
4 * *
5 * Distributed under the terms of the BSD 3-Clause License. *
6 * *
7 * The full license is in the file LICENSE, distributed with this software. *
8 ****************************************************************************/
9
14#ifndef XTENSOR_RANDOM_HPP
15#define XTENSOR_RANDOM_HPP
16
17#include <algorithm>
18#include <functional>
19#include <random>
20#include <type_traits>
21#include <utility>
22
23#include <xtl/xspan.hpp>
24
25#include "xbuilder.hpp"
26#include "xgenerator.hpp"
27#include "xindex_view.hpp"
28#include "xmath.hpp"
29#include "xtensor.hpp"
30#include "xtensor_config.hpp"
31#include "xview.hpp"
32
33namespace xt
34{
35 /*********************
36 * Random generators *
37 *********************/
38
39 namespace random
40 {
41 using default_engine_type = std::mt19937;
42 using seed_type = default_engine_type::result_type;
43
44 default_engine_type& get_default_random_engine();
45 void seed(seed_type seed);
46
47 template <class T, class S, class E = random::default_engine_type>
48 auto rand(const S& shape, T lower = 0, T upper = 1, E& engine = random::get_default_random_engine());
49
50 template <class T, class S, class E = random::default_engine_type>
51 auto randint(
52 const S& shape,
53 T lower = 0,
54 T upper = (std::numeric_limits<T>::max)(),
55 E& engine = random::get_default_random_engine()
56 );
57
58 template <class T, class S, class E = random::default_engine_type>
59 auto randn(const S& shape, T mean = 0, T std_dev = 1, E& engine = random::get_default_random_engine());
60
61 template <class T, class S, class D = double, class E = random::default_engine_type>
62 auto
63 binomial(const S& shape, T trials = 1, D prob = 0.5, E& engine = random::get_default_random_engine());
64
65 template <class T, class S, class D = double, class E = random::default_engine_type>
66 auto geometric(const S& shape, D prob = 0.5, E& engine = random::get_default_random_engine());
67
68 template <class T, class S, class D = double, class E = random::default_engine_type>
69 auto
70 negative_binomial(const S& shape, T k = 1, D prob = 0.5, E& engine = random::get_default_random_engine());
71
72 template <class T, class S, class D = double, class E = random::default_engine_type>
73 auto poisson(const S& shape, D rate = 1.0, E& engine = random::get_default_random_engine());
74
75 template <class T, class S, class E = random::default_engine_type>
76 auto exponential(const S& shape, T rate = 1.0, E& engine = random::get_default_random_engine());
77
78 template <class T, class S, class E = random::default_engine_type>
79 auto
80 gamma(const S& shape, T alpha = 1.0, T beta = 1.0, E& engine = random::get_default_random_engine());
81
82 template <class T, class S, class E = random::default_engine_type>
83 auto weibull(const S& shape, T a = 1.0, T b = 1.0, E& engine = random::get_default_random_engine());
84
85 template <class T, class S, class E = random::default_engine_type>
86 auto
87 extreme_value(const S& shape, T a = 0.0, T b = 1.0, E& engine = random::get_default_random_engine());
88
89 template <class T, class S, class E = random::default_engine_type>
90 auto
91 lognormal(const S& shape, T mean = 0, T std_dev = 1, E& engine = random::get_default_random_engine());
92
93 template <class T, class S, class E = random::default_engine_type>
94 auto chi_squared(const S& shape, T deg = 1.0, E& engine = random::get_default_random_engine());
95
96 template <class T, class S, class E = random::default_engine_type>
97 auto cauchy(const S& shape, T a = 0.0, T b = 1.0, E& engine = random::get_default_random_engine());
98
99 template <class T, class S, class E = random::default_engine_type>
100 auto fisher_f(const S& shape, T m = 1.0, T n = 1.0, E& engine = random::get_default_random_engine());
101
102 template <class T, class S, class E = random::default_engine_type>
103 auto student_t(const S& shape, T n = 1.0, E& engine = random::get_default_random_engine());
104
105 template <class T, class I, std::size_t L, class E = random::default_engine_type>
106 auto
107 rand(const I (&shape)[L], T lower = 0, T upper = 1, E& engine = random::get_default_random_engine());
108
109 template <class T, class I, std::size_t L, class E = random::default_engine_type>
110 auto randint(
111 const I (&shape)[L],
112 T lower = 0,
113 T upper = (std::numeric_limits<T>::max)(),
114 E& engine = random::get_default_random_engine()
115 );
116
117 template <class T, class I, std::size_t L, class E = random::default_engine_type>
118 auto
119 randn(const I (&shape)[L], T mean = 0, T std_dev = 1, E& engine = random::get_default_random_engine());
120
121 template <class T, class I, std::size_t L, class D = double, class E = random::default_engine_type>
122 auto
123 binomial(const I (&shape)[L], T trials = 1, D prob = 0.5, E& engine = random::get_default_random_engine());
124
125 template <class T, class I, std::size_t L, class D = double, class E = random::default_engine_type>
126 auto geometric(const I (&shape)[L], D prob = 0.5, E& engine = random::get_default_random_engine());
127
128 template <class T, class I, std::size_t L, class D = double, class E = random::default_engine_type>
129 auto negative_binomial(
130 const I (&shape)[L],
131 T k = 1,
132 D prob = 0.5,
133 E& engine = random::get_default_random_engine()
134 );
135
136 template <class T, class I, std::size_t L, class D = double, class E = random::default_engine_type>
137 auto poisson(const I (&shape)[L], D rate = 1.0, E& engine = random::get_default_random_engine());
138
139 template <class T, class I, std::size_t L, class E = random::default_engine_type>
140 auto exponential(const I (&shape)[L], T rate = 1.0, E& engine = random::get_default_random_engine());
141
142 template <class T, class I, std::size_t L, class E = random::default_engine_type>
143 auto
144 gamma(const I (&shape)[L], T alpha = 1.0, T beta = 1.0, E& engine = random::get_default_random_engine());
145
146 template <class T, class I, std::size_t L, class E = random::default_engine_type>
147 auto
148 weibull(const I (&shape)[L], T a = 1.0, T b = 1.0, E& engine = random::get_default_random_engine());
149
150 template <class T, class I, std::size_t L, class E = random::default_engine_type>
151 auto
152 extreme_value(const I (&shape)[L], T a = 0.0, T b = 1.0, E& engine = random::get_default_random_engine());
153
154 template <class T, class I, std::size_t L, class E = random::default_engine_type>
155 auto lognormal(
156 const I (&shape)[L],
157 T mean = 0.0,
158 T std_dev = 1.0,
159 E& engine = random::get_default_random_engine()
160 );
161
162 template <class T, class I, std::size_t L, class E = random::default_engine_type>
163 auto chi_squared(const I (&shape)[L], T deg = 1.0, E& engine = random::get_default_random_engine());
164
165 template <class T, class I, std::size_t L, class E = random::default_engine_type>
166 auto cauchy(const I (&shape)[L], T a = 0.0, T b = 1.0, E& engine = random::get_default_random_engine());
167
168 template <class T, class I, std::size_t L, class E = random::default_engine_type>
169 auto
170 fisher_f(const I (&shape)[L], T m = 1.0, T n = 1.0, E& engine = random::get_default_random_engine());
171
172 template <class T, class I, std::size_t L, class E = random::default_engine_type>
173 auto student_t(const I (&shape)[L], T n = 1.0, E& engine = random::get_default_random_engine());
174
175 template <class T, class E = random::default_engine_type>
176 void shuffle(xexpression<T>& e, E& engine = random::get_default_random_engine());
177
178 template <class T, class E = random::default_engine_type>
179 std::enable_if_t<xtl::is_integral<T>::value, xtensor<T, 1>>
180 permutation(T e, E& engine = random::get_default_random_engine());
181
182 template <class T, class E = random::default_engine_type>
183 std::enable_if_t<is_xexpression<std::decay_t<T>>::value, std::decay_t<T>>
184 permutation(T&& e, E& engine = random::get_default_random_engine());
185
186 template <class T, class E = random::default_engine_type>
187 xtensor<typename T::value_type, 1> choice(
188 const xexpression<T>& e,
189 std::size_t n,
190 bool replace = true,
191 E& engine = random::get_default_random_engine()
192 );
193
194 template <class T, class W, class E = random::default_engine_type>
195 xtensor<typename T::value_type, 1> choice(
196 const xexpression<T>& e,
197 std::size_t n,
198 const xexpression<W>& weights,
199 bool replace = true,
200 E& engine = random::get_default_random_engine()
201 );
202 }
203
204 namespace detail
205 {
206 template <class T, class E, class D>
207 struct random_impl
208 {
209 using value_type = T;
210
211 random_impl(E& engine, D&& dist)
212 : m_engine(engine)
213 , m_dist(std::move(dist))
214 {
215 }
216
217 template <class... Args>
218 inline value_type operator()(Args...) const
219 {
220 return m_dist(m_engine);
221 }
222
223 template <class It>
224 inline value_type element(It, It) const
225 {
226 return m_dist(m_engine);
227 }
228
229 template <class EX>
230 inline void assign_to(xexpression<EX>& e) const noexcept
231 {
232 // Note: we're not going row/col major here
233 auto& ed = e.derived_cast();
234 for (auto&& el : ed.storage())
235 {
236 el = m_dist(m_engine);
237 }
238 }
239
240 private:
241
242 E& m_engine;
243 mutable D m_dist;
244 };
245 }
246
247 namespace random
248 {
252 inline default_engine_type& get_default_random_engine()
253 {
254 static default_engine_type mt;
255 return mt;
256 }
257
262 inline void seed(seed_type seed)
263 {
264 get_default_random_engine().seed(seed);
265 }
266
279 template <class T, class S, class E>
280 inline auto rand(const S& shape, T lower, T upper, E& engine)
281 {
282 std::uniform_real_distribution<T> dist(lower, upper);
283 return detail::make_xgenerator(
284 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
285 shape
286 );
287 }
288
301 template <class T, class S, class E>
302 inline auto randint(const S& shape, T lower, T upper, E& engine)
303 {
304 std::uniform_int_distribution<T> dist(lower, T(upper - 1));
305 return detail::make_xgenerator(
306 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
307 shape
308 );
309 }
310
324 template <class T, class S, class E>
325 inline auto randn(const S& shape, T mean, T std_dev, E& engine)
326 {
327 std::normal_distribution<T> dist(mean, std_dev);
328 return detail::make_xgenerator(
329 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
330 shape
331 );
332 }
333
347 template <class T, class S, class D, class E>
348 inline auto binomial(const S& shape, T trials, D prob, E& engine)
349 {
350 std::binomial_distribution<T> dist(trials, prob);
351 return detail::make_xgenerator(
352 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
353 shape
354 );
355 }
356
369 template <class T, class S, class D, class E>
370 inline auto geometric(const S& shape, D prob, E& engine)
371 {
372 std::geometric_distribution<T> dist(prob);
373 return detail::make_xgenerator(
374 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
375 shape
376 );
377 }
378
393 template <class T, class S, class D, class E>
394 inline auto negative_binomial(const S& shape, T k, D prob, E& engine)
395 {
396 std::negative_binomial_distribution<T> dist(k, prob);
397 return detail::make_xgenerator(
398 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
399 shape
400 );
401 }
402
414 template <class T, class S, class D, class E>
415 inline auto poisson(const S& shape, D rate, E& engine)
416 {
417 std::poisson_distribution<T> dist(rate);
418 return detail::make_xgenerator(
419 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
420 shape
421 );
422 }
423
435 template <class T, class S, class E>
436 inline auto exponential(const S& shape, T rate, E& engine)
437 {
438 std::exponential_distribution<T> dist(rate);
439 return detail::make_xgenerator(
440 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
441 shape
442 );
443 }
444
457 template <class T, class S, class E>
458 inline auto gamma(const S& shape, T alpha, T beta, E& engine)
459 {
460 std::gamma_distribution<T> dist(alpha, beta);
461 return detail::make_xgenerator(
462 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
463 shape
464 );
465 }
466
479 template <class T, class S, class E>
480 inline auto weibull(const S& shape, T a, T b, E& engine)
481 {
482 std::weibull_distribution<T> dist(a, b);
483 return detail::make_xgenerator(
484 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
485 shape
486 );
487 }
488
501 template <class T, class S, class E>
502 inline auto extreme_value(const S& shape, T a, T b, E& engine)
503 {
504 std::extreme_value_distribution<T> dist(a, b);
505 return detail::make_xgenerator(
506 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
507 shape
508 );
509 }
510
524 template <class T, class S, class E>
525 inline auto lognormal(const S& shape, T mean, T std_dev, E& engine)
526 {
527 std::lognormal_distribution<T> dist(mean, std_dev);
528 return detail::make_xgenerator(
529 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
530 shape
531 );
532 }
533
545 template <class T, class S, class E>
546 inline auto chi_squared(const S& shape, T deg, E& engine)
547 {
548 std::chi_squared_distribution<T> dist(deg);
549 return detail::make_xgenerator(
550 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
551 shape
552 );
553 }
554
567 template <class T, class S, class E>
568 inline auto cauchy(const S& shape, T a, T b, E& engine)
569 {
570 std::cauchy_distribution<T> dist(a, b);
571 return detail::make_xgenerator(
572 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
573 shape
574 );
575 }
576
590 template <class T, class S, class E>
591 inline auto fisher_f(const S& shape, T m, T n, E& engine)
592 {
593 std::fisher_f_distribution<T> dist(m, n);
594 return detail::make_xgenerator(
595 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
596 shape
597 );
598 }
599
612 template <class T, class S, class E>
613 inline auto student_t(const S& shape, T n, E& engine)
614 {
615 std::student_t_distribution<T> dist(n);
616 return detail::make_xgenerator(
617 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
618 shape
619 );
620 }
621
622 template <class T, class I, std::size_t L, class E>
623 inline auto rand(const I (&shape)[L], T lower, T upper, E& engine)
624 {
625 std::uniform_real_distribution<T> dist(lower, upper);
626 return detail::make_xgenerator(
627 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
628 shape
629 );
630 }
631
632 template <class T, class I, std::size_t L, class E>
633 inline auto randint(const I (&shape)[L], T lower, T upper, E& engine)
634 {
635 std::uniform_int_distribution<T> dist(lower, T(upper - 1));
636 return detail::make_xgenerator(
637 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
638 shape
639 );
640 }
641
642 template <class T, class I, std::size_t L, class E>
643 inline auto randn(const I (&shape)[L], T mean, T std_dev, E& engine)
644 {
645 std::normal_distribution<T> dist(mean, std_dev);
646 return detail::make_xgenerator(
647 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
648 shape
649 );
650 }
651
652 template <class T, class I, std::size_t L, class D, class E>
653 inline auto binomial(const I (&shape)[L], T trials, D prob, E& engine)
654 {
655 std::binomial_distribution<T> dist(trials, prob);
656 return detail::make_xgenerator(
657 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
658 shape
659 );
660 }
661
662 template <class T, class I, std::size_t L, class D, class E>
663 inline auto geometric(const I (&shape)[L], D prob, E& engine)
664 {
665 std::geometric_distribution<T> dist(prob);
666 return detail::make_xgenerator(
667 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
668 shape
669 );
670 }
671
672 template <class T, class I, std::size_t L, class D, class E>
673 inline auto negative_binomial(const I (&shape)[L], T k, D prob, E& engine)
674 {
675 std::negative_binomial_distribution<T> dist(k, prob);
676 return detail::make_xgenerator(
677 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
678 shape
679 );
680 }
681
682 template <class T, class I, std::size_t L, class D, class E>
683 inline auto poisson(const I (&shape)[L], D rate, E& engine)
684 {
685 std::poisson_distribution<T> dist(rate);
686 return detail::make_xgenerator(
687 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
688 shape
689 );
690 }
691
692 template <class T, class I, std::size_t L, class E>
693 inline auto exponential(const I (&shape)[L], T rate, E& engine)
694 {
695 std::exponential_distribution<T> dist(rate);
696 return detail::make_xgenerator(
697 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
698 shape
699 );
700 }
701
702 template <class T, class I, std::size_t L, class E>
703 inline auto gamma(const I (&shape)[L], T alpha, T beta, E& engine)
704 {
705 std::gamma_distribution<T> dist(alpha, beta);
706 return detail::make_xgenerator(
707 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
708 shape
709 );
710 }
711
712 template <class T, class I, std::size_t L, class E>
713 inline auto weibull(const I (&shape)[L], T a, T b, E& engine)
714 {
715 std::weibull_distribution<T> dist(a, b);
716 return detail::make_xgenerator(
717 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
718 shape
719 );
720 }
721
722 template <class T, class I, std::size_t L, class E>
723 inline auto extreme_value(const I (&shape)[L], T a, T b, E& engine)
724 {
725 std::extreme_value_distribution<T> dist(a, b);
726 return detail::make_xgenerator(
727 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
728 shape
729 );
730 }
731
732 template <class T, class I, std::size_t L, class E>
733 inline auto lognormal(const I (&shape)[L], T mean, T std_dev, E& engine)
734 {
735 std::lognormal_distribution<T> dist(mean, std_dev);
736 return detail::make_xgenerator(
737 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
738 shape
739 );
740 }
741
742 template <class T, class I, std::size_t L, class E>
743 inline auto chi_squared(const I (&shape)[L], T deg, E& engine)
744 {
745 std::chi_squared_distribution<T> dist(deg);
746 return detail::make_xgenerator(
747 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
748 shape
749 );
750 }
751
752 template <class T, class I, std::size_t L, class E>
753 inline auto cauchy(const I (&shape)[L], T a, T b, E& engine)
754 {
755 std::cauchy_distribution<T> dist(a, b);
756 return detail::make_xgenerator(
757 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
758 shape
759 );
760 }
761
762 template <class T, class I, std::size_t L, class E>
763 inline auto fisher_f(const I (&shape)[L], T m, T n, E& engine)
764 {
765 std::fisher_f_distribution<T> dist(m, n);
766 return detail::make_xgenerator(
767 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
768 shape
769 );
770 }
771
772 template <class T, class I, std::size_t L, class E>
773 inline auto student_t(const I (&shape)[L], T n, E& engine)
774 {
775 std::student_t_distribution<T> dist(n);
776 return detail::make_xgenerator(
777 detail::random_impl<T, E, decltype(dist)>(engine, std::move(dist)),
778 shape
779 );
780 }
781
789 template <class T, class E>
790 void shuffle(xexpression<T>& e, E& engine)
791 {
792 T& de = e.derived_cast();
793
794 if (de.dimension() == 1)
795 {
796 using size_type = typename T::size_type;
797 auto first = de.begin();
798 auto last = de.end();
799
800 for (size_type i = std::size_t((last - first) - 1); i > 0; --i)
801 {
802 std::uniform_int_distribution<size_type> dist(0, i);
803 auto j = dist(engine);
804 using std::swap;
805 swap(first[i], first[j]);
806 }
807 }
808 else
809 {
810 using size_type = typename T::size_type;
811 decltype(auto) buf = empty_like(view(de, 0));
812
813 for (size_type i = de.shape()[0] - 1; i > 0; --i)
814 {
815 std::uniform_int_distribution<size_type> dist(0, i);
816 size_type j = dist(engine);
817
818 buf = view(de, j);
819 view(de, j) = view(de, i);
820 view(de, i) = buf;
821 }
822 }
823 }
824
838 template <class T, class E>
839 std::enable_if_t<xtl::is_integral<T>::value, xtensor<T, 1>> permutation(T e, E& engine)
840 {
842 shuffle(res, engine);
843 return res;
844 }
845
847 template <class T, class E>
848 std::enable_if_t<is_xexpression<std::decay_t<T>>::value, std::decay_t<T>> permutation(T&& e, E& engine)
849 {
850 using copy_type = std::decay_t<T>;
851 copy_type res = e;
852 shuffle(res, engine);
853 return res;
854 }
855
857
869 template <class T, class E>
870 xtensor<typename T::value_type, 1>
871 choice(const xexpression<T>& e, std::size_t n, bool replace, E& engine)
872 {
873 const auto& de = e.derived_cast();
874 XTENSOR_ASSERT((de.dimension() == 1));
875 XTENSOR_ASSERT((replace || n <= de.size()));
876 using result_type = xtensor<typename T::value_type, 1>;
877 using size_type = typename result_type::size_type;
878 result_type result;
879 result.resize({n});
880
881 if (replace)
882 {
883 auto dist = std::uniform_int_distribution<size_type>(0, de.size() - 1);
884 for (size_type i = 0; i < n; ++i)
885 {
886 result[i] = de.storage()[dist(engine)];
887 }
888 }
889 else
890 {
891 // Naive resevoir sampling without weighting:
892 std::copy(de.storage().begin(), de.storage().begin() + n, result.begin());
893 size_type i = n;
894 for (auto it = de.storage().begin() + n; it != de.storage().end(); ++it, ++i)
895 {
896 auto idx = std::uniform_int_distribution<size_type>(0, i)(engine);
897 if (idx < n)
898 {
899 result.storage()[idx] = *it;
900 }
901 }
902 }
903 return result;
904 }
905
933 template <class T, class W, class E>
934 xtensor<typename T::value_type, 1>
935 choice(const xexpression<T>& e, std::size_t n, const xexpression<W>& weights, bool replace, E& engine)
936 {
937 const auto& de = e.derived_cast();
938 const auto& dweights = weights.derived_cast();
939 XTENSOR_ASSERT((de.dimension() == 1));
940 XTENSOR_ASSERT((replace || n <= de.size()));
941 XTENSOR_ASSERT((de.size() == dweights.size()));
942 XTENSOR_ASSERT((de.dimension() == dweights.dimension()));
943 XTENSOR_ASSERT(xt::all(dweights >= 0));
944 static_assert(
945 std::is_floating_point<typename W::value_type>::value,
946 "Weight expression must be of floating point type"
947 );
948 using result_type = xtensor<typename T::value_type, 1>;
949 using size_type = typename result_type::size_type;
950 using weight_type = typename W::value_type;
951 result_type result;
952 result.resize({n});
953
954 if (replace)
955 {
956 // Sample u uniformly in the range [0, sum(weights)[
957 // The index idx of the sampled element is such that weight_cumul[idx - 1] <= u <
958 // weight_cumul[idx]. Where weight_cumul[-1] is implicitly 0, as the empty sum.
959 const auto wc = eval(cumsum(dweights));
960 std::uniform_real_distribution<weight_type> weight_dist{0, wc[wc.size() - 1]};
961 for (auto& x : result)
962 {
963 const auto u = weight_dist(engine);
964 const auto idx = static_cast<size_type>(
965 std::upper_bound(wc.cbegin(), wc.cend(), u) - wc.cbegin()
966 );
967 x = de[idx];
968 }
969 }
970 else
971 {
972 // Compute (modified) keys as weight/randexp(1).
973 xtensor<weight_type, 1> keys;
974 keys.resize({dweights.size()});
975 std::exponential_distribution<weight_type> randexp{weight_type(1)};
976 std::transform(
977 dweights.cbegin(),
978 dweights.cend(),
979 keys.begin(),
980 [&randexp, &engine](auto w)
981 {
982 return w / randexp(engine);
983 }
984 );
985
986 // Find indexes for the n biggest key
987 xtensor<size_type, 1> indices = arange<size_type>(0, dweights.size());
988 std::partial_sort(
989 indices.begin(),
990 indices.begin() + n,
991 indices.end(),
992 [&keys](auto i, auto j)
993 {
994 return keys[i] > keys[j];
995 }
996 );
997
998 // Return samples with the n biggest keys
999 result = index_view(de, xtl::span<size_type>{indices.data(), n});
1000 }
1001 return result;
1002 }
1003
1004 }
1005}
1006
1007#endif
auto cumsum(E &&e, std::ptrdiff_t axis)
Cumulative sum.
Definition xmath.hpp:2286
auto mean(E &&e, X &&axes, EVS es=EVS())
Mean of elements over given axes.
Definition xmath.hpp:1935
auto eval(T &&t) -> std::enable_if_t< detail::is_container< std::decay_t< T > >::value, T && >
Force evaluation of xexpression.
Definition xeval.hpp:46
@ weibull
Method 6 of (Hyndman and Fan, 1996) with alpha=0 and beta=0.
standard mathematical functions for xexpressions
auto all() noexcept
Returns a slice representing a full dimension, to be used as an argument of view function.
Definition xslice.hpp:234
auto index_view(E &&e, I &&indices) noexcept
creates an indexview from a container of indices.
auto empty_like(const xexpression< E > &e)
Create a xcontainer (xarray, xtensor or xtensor_fixed) with uninitialized values of the same shape,...
Definition xbuilder.hpp:121
auto view(E &&e, S &&... slices)
Constructs and returns a view on the specified xexpression.
Definition xview.hpp:1834