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