14#ifndef XTENSOR_RANDOM_HPP
15#define XTENSOR_RANDOM_HPP
23#include <xtl/xspan.hpp>
25#include "xbuilder.hpp"
26#include "xgenerator.hpp"
27#include "xindex_view.hpp"
30#include "xtensor_config.hpp"
41 using default_engine_type = std::mt19937;
42 using seed_type = default_engine_type::result_type;
44 default_engine_type& get_default_random_engine();
45 void seed(seed_type seed);
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());
50 template <
class T,
class S,
class E = random::default_engine_type>
54 T upper = (std::numeric_limits<T>::max)(),
55 E& engine = random::get_default_random_engine()
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());
61 template <
class T,
class S,
class D =
double,
class E = random::default_engine_type>
63 binomial(
const S& shape, T trials = 1, D prob = 0.5, E& engine = random::get_default_random_engine());
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());
68 template <
class T,
class S,
class D =
double,
class E = random::default_engine_type>
70 negative_binomial(
const S& shape, T k = 1, D prob = 0.5, E& engine = random::get_default_random_engine());
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());
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());
78 template <
class T,
class S,
class E = random::default_engine_type>
80 gamma(
const S& shape, T alpha = 1.0, T beta = 1.0, E& engine = random::get_default_random_engine());
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());
85 template <
class T,
class S,
class E = random::default_engine_type>
87 extreme_value(
const S& shape, T a = 0.0, T b = 1.0, E& engine = random::get_default_random_engine());
89 template <
class T,
class S,
class E = random::default_engine_type>
91 lognormal(
const S& shape, T
mean = 0, T std_dev = 1, E& engine = random::get_default_random_engine());
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());
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());
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());
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());
105 template <
class T,
class I, std::
size_t L,
class E = random::default_engine_type>
107 rand(
const I (&shape)[L], T lower = 0, T upper = 1, E& engine = random::get_default_random_engine());
109 template <
class T,
class I, std::
size_t L,
class E = random::default_engine_type>
113 T upper = (std::numeric_limits<T>::max)(),
114 E& engine = random::get_default_random_engine()
117 template <
class T,
class I, std::
size_t L,
class E = random::default_engine_type>
119 randn(
const I (&shape)[L], T
mean = 0, T std_dev = 1, E& engine = random::get_default_random_engine());
121 template <
class T,
class I, std::
size_t L,
class D =
double,
class E = random::default_engine_type>
123 binomial(
const I (&shape)[L], T trials = 1, D prob = 0.5, E& engine = random::get_default_random_engine());
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());
128 template <
class T,
class I, std::
size_t L,
class D =
double,
class E = random::default_engine_type>
129 auto negative_binomial(
133 E& engine = random::get_default_random_engine()
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());
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());
142 template <
class T,
class I, std::
size_t L,
class E = random::default_engine_type>
144 gamma(
const I (&shape)[L], T alpha = 1.0, T beta = 1.0, E& engine = random::get_default_random_engine());
146 template <
class T,
class I, std::
size_t L,
class E = random::default_engine_type>
148 weibull(
const I (&shape)[L], T a = 1.0, T b = 1.0, E& engine = random::get_default_random_engine());
150 template <
class T,
class I, std::
size_t L,
class E = random::default_engine_type>
152 extreme_value(
const I (&shape)[L], T a = 0.0, T b = 1.0, E& engine = random::get_default_random_engine());
154 template <
class T,
class I, std::
size_t L,
class E = random::default_engine_type>
159 E& engine = random::get_default_random_engine()
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());
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());
168 template <
class T,
class I, std::
size_t L,
class E = random::default_engine_type>
170 fisher_f(
const I (&shape)[L], T m = 1.0, T n = 1.0, E& engine = random::get_default_random_engine());
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());
175 template <
class T,
class E = random::default_engine_type>
176 void shuffle(xexpression<T>& e, E& engine = random::get_default_random_engine());
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());
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());
186 template <
class T,
class E = random::default_engine_type>
187 xtensor<typename T::value_type, 1> choice(
188 const xexpression<T>& e,
191 E& engine = random::get_default_random_engine()
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,
198 const xexpression<W>& weights,
200 E& engine = random::get_default_random_engine()
206 template <
class T,
class E,
class D>
209 using value_type = T;
211 random_impl(E& engine, D&& dist)
213 , m_dist(std::move(dist))
217 template <
class... Args>
218 inline value_type operator()(Args...)
const
220 return m_dist(m_engine);
224 inline value_type element(It, It)
const
226 return m_dist(m_engine);
230 inline void assign_to(xexpression<EX>& e)
const noexcept
233 auto& ed = e.derived_cast();
234 for (
auto&& el : ed.storage())
236 el = m_dist(m_engine);
252 inline default_engine_type& get_default_random_engine()
254 static default_engine_type mt;
262 inline void seed(seed_type seed)
264 get_default_random_engine().seed(seed);
279 template <
class T,
class S,
class E>
280 inline auto rand(
const S& shape, T lower, T upper, E& engine)
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)),
301 template <
class T,
class S,
class E>
302 inline auto randint(
const S& shape, T lower, T upper, E& engine)
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)),
324 template <
class T,
class S,
class E>
325 inline auto randn(
const S& shape, T
mean, T std_dev, E& engine)
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)),
347 template <
class T,
class S,
class D,
class E>
348 inline auto binomial(
const S& shape, T trials, D prob, E& engine)
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)),
369 template <
class T,
class S,
class D,
class E>
370 inline auto geometric(
const S& shape, D prob, E& engine)
372 std::geometric_distribution<T> dist(prob);
373 return detail::make_xgenerator(
374 detail::random_impl<T, E,
decltype(dist)>(engine, std::move(dist)),
393 template <
class T,
class S,
class D,
class E>
394 inline auto negative_binomial(
const S& shape, T k, D prob, E& engine)
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)),
414 template <
class T,
class S,
class D,
class E>
415 inline auto poisson(
const S& shape, D rate, E& engine)
417 std::poisson_distribution<T> dist(rate);
418 return detail::make_xgenerator(
419 detail::random_impl<T, E,
decltype(dist)>(engine, std::move(dist)),
435 template <
class T,
class S,
class E>
436 inline auto exponential(
const S& shape, T rate, E& engine)
438 std::exponential_distribution<T> dist(rate);
439 return detail::make_xgenerator(
440 detail::random_impl<T, E,
decltype(dist)>(engine, std::move(dist)),
457 template <
class T,
class S,
class E>
458 inline auto gamma(
const S& shape, T alpha, T beta, E& engine)
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)),
479 template <
class T,
class S,
class E>
480 inline auto weibull(
const S& shape, T a, T b, E& engine)
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)),
501 template <
class T,
class S,
class E>
502 inline auto extreme_value(
const S& shape, T a, T b, E& engine)
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)),
524 template <
class T,
class S,
class E>
525 inline auto lognormal(
const S& shape, T
mean, T std_dev, E& engine)
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)),
545 template <
class T,
class S,
class E>
546 inline auto chi_squared(
const S& shape, T deg, E& engine)
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)),
567 template <
class T,
class S,
class E>
568 inline auto cauchy(
const S& shape, T a, T b, E& engine)
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)),
590 template <
class T,
class S,
class E>
591 inline auto fisher_f(
const S& shape, T m, T n, E& engine)
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)),
612 template <
class T,
class S,
class E>
613 inline auto student_t(
const S& shape, T n, E& engine)
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)),
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)
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)),
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)
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)),
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)
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)),
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)
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)),
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)
665 std::geometric_distribution<T> dist(prob);
666 return detail::make_xgenerator(
667 detail::random_impl<T, E,
decltype(dist)>(engine, std::move(dist)),
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)
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)),
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)
685 std::poisson_distribution<T> dist(rate);
686 return detail::make_xgenerator(
687 detail::random_impl<T, E,
decltype(dist)>(engine, std::move(dist)),
692 template <
class T,
class I, std::
size_t L,
class E>
693 inline auto exponential(
const I (&shape)[L], T rate, E& engine)
695 std::exponential_distribution<T> dist(rate);
696 return detail::make_xgenerator(
697 detail::random_impl<T, E,
decltype(dist)>(engine, std::move(dist)),
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)
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)),
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)
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)),
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)
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)),
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)
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)),
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)
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)),
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)
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)),
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)
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)),
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)
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)),
789 template <
class T,
class E>
790 void shuffle(xexpression<T>& e, E& engine)
792 T& de = e.derived_cast();
794 if (de.dimension() == 1)
796 using size_type =
typename T::size_type;
797 auto first = de.begin();
798 auto last = de.end();
800 for (size_type i = std::size_t((last - first) - 1); i > 0; --i)
802 std::uniform_int_distribution<size_type> dist(0, i);
803 auto j = dist(engine);
805 swap(first[i], first[j]);
810 using size_type =
typename T::size_type;
813 for (size_type i = de.shape()[0] - 1; i > 0; --i)
815 std::uniform_int_distribution<size_type> dist(0, i);
816 size_type j = dist(engine);
838 template <
class T,
class E>
839 std::enable_if_t<xtl::is_integral<T>::value, xtensor<T, 1>> permutation(T e, E& engine)
842 shuffle(res, engine);
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)
850 using copy_type = std::decay_t<T>;
852 shuffle(res, engine);
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)
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;
883 auto dist = std::uniform_int_distribution<size_type>(0, de.size() - 1);
884 for (size_type i = 0; i < n; ++i)
886 result[i] = de.storage()[dist(engine)];
892 std::copy(de.storage().begin(), de.storage().begin() + n, result.begin());
894 for (
auto it = de.storage().begin() + n; it != de.storage().end(); ++it, ++i)
896 auto idx = std::uniform_int_distribution<size_type>(0, i)(engine);
899 result.storage()[idx] = *it;
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)
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));
945 std::is_floating_point<typename W::value_type>::value,
946 "Weight expression must be of floating point type"
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;
960 std::uniform_real_distribution<weight_type> weight_dist{0, wc[wc.size() - 1]};
961 for (
auto& x : result)
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()
973 xtensor<weight_type, 1> keys;
974 keys.resize({dweights.size()});
975 std::exponential_distribution<weight_type> randexp{weight_type(1)};
980 [&randexp, &engine](
auto w)
982 return w / randexp(engine);
987 xtensor<size_type, 1> indices = arange<size_type>(0, dweights.size());
992 [&keys](
auto i,
auto j)
994 return keys[i] > keys[j];
999 result =
index_view(de, xtl::span<size_type>{indices.data(), n});
auto cumsum(E &&e, std::ptrdiff_t axis)
Cumulative sum.
auto mean(E &&e, X &&axes, EVS es=EVS())
Mean of elements over given axes.
auto eval(T &&t) -> std::enable_if_t< detail::is_container< std::decay_t< T > >::value, T && >
Force evaluation of xexpression.
@ 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.
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,...
auto view(E &&e, S &&... slices)
Constructs and returns a view on the specified xexpression.