diff --git a/include/armadillo_bits/arma_rng.hpp b/include/armadillo_bits/arma_rng.hpp index a569172e..76946a05 100644 --- a/include/armadillo_bits/arma_rng.hpp +++ b/include/armadillo_bits/arma_rng.hpp @@ -23,21 +23,20 @@ #endif -// workaround for issue on macOS 11 and/or AppleClang 12.0 -// see https://gitlab.com/conradsnicta/armadillo-code/-/issues/173 -// the workaround is here instead of CMakeLists.txt -// to ensure that the armadillo runtime library has arma_rng_cxx11_instance -// for already compiled programs running on earlier versions of macOS -#if defined(__APPLE__) || defined(__apple_build_version__) - #if !defined(ARMA_DONT_DISABLE_EXTERN_RNG) - #undef ARMA_USE_EXTERN_RNG - #endif -#endif +// // workaround for issue on macOS 11 and/or AppleClang 12.0 +// // see https://gitlab.com/conradsnicta/armadillo-code/-/issues/173 +// // the workaround is here instead of CMakeLists.txt +// // to ensure that the armadillo runtime library has arma_rng_cxx11_instance +// // for already compiled programs running on earlier versions of macOS +// #if defined(__APPLE__) || defined(__apple_build_version__) +// #if !defined(ARMA_DONT_DISABLE_EXTERN_RNG) +// #undef ARMA_USE_EXTERN_RNG +// #endif +// #endif #if defined(ARMA_USE_EXTERN_RNG) - extern thread_local arma_rng_cxx11 arma_rng_cxx11_instance; - // namespace { thread_local arma_rng_cxx11 arma_rng_cxx11_instance; } + extern thread_local std::mt19937_64 mt19937_64_instance; #endif @@ -46,11 +45,11 @@ class arma_rng public: #if defined(ARMA_RNG_ALT) - typedef arma_rng_alt::seed_type seed_type; + typedef arma_rng_alt::seed_type seed_type; #elif defined(ARMA_USE_EXTERN_RNG) - typedef arma_rng_cxx11::seed_type seed_type; + typedef std::mt19937_64::result_type seed_type; #else - typedef arma_rng_cxx98::seed_type seed_type; + typedef arma_rng_cxx98::seed_type seed_type; #endif #if defined(ARMA_RNG_ALT) @@ -67,6 +66,7 @@ class arma_rng template struct randi; template struct randu; template struct randn; + template struct randg; }; @@ -81,7 +81,7 @@ arma_rng::set_seed(const arma_rng::seed_type val) } #elif defined(ARMA_USE_EXTERN_RNG) { - arma_rng_cxx11_instance.set_seed(val); + mt19937_64_instance.seed(val); } #else { @@ -173,10 +173,14 @@ arma_rng::set_seed_random() +// + + + template struct arma_rng::randi { - arma_inline + inline operator eT () { #if defined(ARMA_RNG_ALT) @@ -185,7 +189,9 @@ struct arma_rng::randi } #elif defined(ARMA_USE_EXTERN_RNG) { - return eT( arma_rng_cxx11_instance.randi_val() ); + constexpr double scale = double(std::numeric_limits::max()) / double(std::mt19937_64::max()); + + return eT( double(mt19937_64_instance()) * scale ); } #else { @@ -206,7 +212,7 @@ struct arma_rng::randi } #elif defined(ARMA_USE_EXTERN_RNG) { - return arma_rng_cxx11::randi_max_val(); + return std::numeric_limits::max(); } #else { @@ -227,7 +233,23 @@ struct arma_rng::randi } #elif defined(ARMA_USE_EXTERN_RNG) { - arma_rng_cxx11_instance.randi_fill(mem, N, a, b); + if(N == uword(1)) + { + std::uniform_int_distribution local_i_distr(a, b); + + mem[0] = eT(local_i_distr(mt19937_64_instance)); + + return; + } + + std::mt19937_64 local_engine; + std::uniform_int_distribution local_i_distr(a, b); + + typedef typename std::mt19937_64::result_type seed_type; + + local_engine.seed( seed_type(mt19937_64_instance()) ); + + for(uword i=0; i struct arma_rng::randu { - arma_inline + inline operator eT () { #if defined(ARMA_RNG_ALT) @@ -260,7 +286,9 @@ struct arma_rng::randu } #elif defined(ARMA_USE_EXTERN_RNG) { - return eT( arma_rng_cxx11_instance.randu_val() ); + constexpr double scale = double(1.0) / double(std::mt19937_64::max()); + + return eT( double(mt19937_64_instance()) * scale ); } #else { @@ -275,9 +303,22 @@ struct arma_rng::randu void fill(eT* mem, const uword N) { - #if defined(ARMA_RNG_ALT) || defined(ARMA_USE_EXTERN_RNG) + #if defined(ARMA_RNG_ALT) { - for(uword i=0; i < N; ++i) { mem[i] = eT( arma_rng::randu() ); } + for(uword i=0; i < N; ++i) { mem[i] = eT( arma_rng_alt::randu_val() ); } + } + #elif defined(ARMA_USE_EXTERN_RNG) + { + if(N == uword(1)) { mem[0] = eT( arma_rng::randu() ); return; } + + typedef typename std::mt19937_64::result_type seed_type; + + std::mt19937_64 local_engine; + std::uniform_real_distribution local_u_distr; + + local_engine.seed( seed_type(mt19937_64_instance()) ); + + for(uword i=0; i < N; ++i) { mem[i] = eT( local_u_distr(local_engine) ); } } #else { @@ -316,13 +357,40 @@ struct arma_rng::randu< std::complex > void fill(std::complex* mem, const uword N) { - #if defined(ARMA_RNG_ALT) || defined(ARMA_USE_EXTERN_RNG) + #if defined(ARMA_RNG_ALT) { for(uword i=0; i < N; ++i) + { + const T a = T( arma_rng_alt::randu_val() ); + const T b = T( arma_rng_alt::randu_val() ); + + mem[i] = std::complex(a, b); + } + } + #elif defined(ARMA_USE_EXTERN_RNG) + { + if(N == uword(1)) { const T a = T( arma_rng::randu() ); const T b = T( arma_rng::randu() ); + mem[0] = std::complex(a, b); + + return; + } + + typedef typename std::mt19937_64::result_type seed_type; + + std::mt19937_64 local_engine; + std::uniform_real_distribution local_u_distr; + + local_engine.seed( seed_type(mt19937_64_instance()) ); + + for(uword i=0; i < N; ++i) + { + const T a = T( local_u_distr(local_engine) ); + const T b = T( local_u_distr(local_engine) ); + mem[i] = std::complex(a, b); } } @@ -359,6 +427,10 @@ struct arma_rng::randu< std::complex > +// + + + template struct arma_rng::randn { @@ -371,7 +443,9 @@ struct arma_rng::randn } #elif defined(ARMA_USE_EXTERN_RNG) { - return eT( arma_rng_cxx11_instance.randn_val() ); + std::normal_distribution local_n_distr; + + return eT( local_n_distr(mt19937_64_instance) ); } #else { @@ -392,7 +466,10 @@ struct arma_rng::randn } #elif defined(ARMA_USE_EXTERN_RNG) { - arma_rng_cxx11_instance.randn_dual_val(out1, out2); + std::normal_distribution local_n_distr; + + out1 = eT( local_n_distr(mt19937_64_instance) ); + out2 = eT( local_n_distr(mt19937_64_instance) ); } #else { @@ -407,19 +484,20 @@ struct arma_rng::randn void fill_simple(eT* mem, const uword N) { - #if defined(ARMA_RNG_ALT) || defined(ARMA_USE_EXTERN_RNG) + #if defined(ARMA_RNG_ALT) { - uword i, j; + for(uword i=0; i < N; ++i) { mem[i] = eT( arma_rng_alt::randn_val() ); } + } + #elif defined(ARMA_USE_EXTERN_RNG) + { + typedef typename std::mt19937_64::result_type seed_type; - for(i=0, j=1; j < N; i+=2, j+=2) - { - arma_rng::randn::dual_val( mem[i], mem[j] ); - } + std::mt19937_64 local_engine; + std::normal_distribution local_n_distr; - if(i < N) - { - mem[i] = eT( arma_rng::randn() ); - } + local_engine.seed( seed_type(mt19937_64_instance()) ); + + for(uword i=0; i < N; ++i) { mem[i] = eT( local_n_distr(local_engine) ); } } #else { @@ -539,10 +617,27 @@ struct arma_rng::randn< std::complex > void fill_simple(std::complex* mem, const uword N) { - #if defined(ARMA_RNG_ALT) || defined(ARMA_USE_EXTERN_RNG) + #if defined(ARMA_RNG_ALT) { for(uword i=0; i < N; ++i) { mem[i] = std::complex( arma_rng::randn< std::complex >() ); } } + #elif defined(ARMA_USE_EXTERN_RNG) + { + typedef typename std::mt19937_64::result_type seed_type; + + std::mt19937_64 local_engine; + std::normal_distribution local_n_distr; + + local_engine.seed( seed_type(mt19937_64_instance()) ); + + for(uword i=0; i < N; ++i) + { + const T a = T( local_n_distr(local_engine) ); + const T b = T( local_n_distr(local_engine) ); + + mem[i] = std::complex(a,b); + } + } #else { if(N == uword(1)) @@ -640,4 +735,111 @@ struct arma_rng::randn< std::complex > +// + + + +template +struct arma_rng::randg + { + inline + static + void + fill_simple(eT* mem, const uword N, const double a, const double b) + { + #if defined(ARMA_USE_EXTERN_RNG) + { + if(N == uword(1)) + { + std::gamma_distribution local_g_distr(a,b); + + mem[0] = eT(local_g_distr(mt19937_64_instance)); + + return; + } + + typedef typename std::mt19937_64::result_type seed_type; + + std::mt19937_64 local_engine; + std::gamma_distribution local_g_distr(a,b); + + local_engine.seed( seed_type(mt19937_64_instance()) ); + + for(uword i=0; i local_g_distr(a,b); + + local_engine.seed( seed_type(std::rand()) ); + + for(uword i=0; i::fill_simple(mem, N, a, b); return; } + + typedef std::mt19937_64 motor_type; + typedef std::mt19937_64::result_type ovum_type; + typedef std::gamma_distribution distr_type; + + const uword n_threads = uword( mp_thread_limit::get() ); + + std::vector g_motor(n_threads); + std::vector g_distr(n_threads); + + const distr_type g_distr_base(a,b); + + for(uword t=0; t < n_threads; ++t) + { + motor_type& g_motor_t = g_motor[t]; + distr_type& g_distr_t = g_distr[t]; + + g_motor_t.seed( ovum_type(t) + ovum_type((*this).randi_val()) ); + + g_distr_t.param( g_distr_base.param() ); + } + + const uword chunk_size = N / n_threads; + + #pragma omp parallel for schedule(static) num_threads(int(n_threads)) + for(uword t=0; t < n_threads; ++t) + { + const uword start = (t+0) * chunk_size; + const uword endp1 = (t+1) * chunk_size; + + motor_type& g_motor_t = g_motor[t]; + distr_type& g_distr_t = g_distr[t]; + + for(uword i=start; i < endp1; ++i) { mem[i] = eT( g_distr_t(g_motor_t)); } + } + + motor_type& g_motor_0 = g_motor[0]; + distr_type& g_distr_0 = g_distr[0]; + + for(uword i=(n_threads*chunk_size); i < N; ++i) { mem[i] = eT( g_distr_0(g_motor_0)); } + } + #else + { + arma_rng::randg::fill_simple(mem, N, a, b); + } + #endif + } + + }; + + + //! @} diff --git a/include/armadillo_bits/fn_randg.hpp b/include/armadillo_bits/fn_randg.hpp index bcb4f6fa..ad344044 100644 --- a/include/armadillo_bits/fn_randg.hpp +++ b/include/armadillo_bits/fn_randg.hpp @@ -28,6 +28,8 @@ randg(const uword n_rows, const uword n_cols, const distr_param& param = distr_p arma_extra_debug_sigprint(); arma_ignore(junk); + typedef typename obj_type::elem_type eT; + if(is_Col::value) { arma_debug_check( (n_cols != 1), "randg(): incompatible size" ); @@ -62,21 +64,7 @@ randg(const uword n_rows, const uword n_cols, const distr_param& param = distr_p arma_debug_check( ((a <= double(0)) || (b <= double(0))), "randg(): a and b must be greater than zero" ); - #if defined(ARMA_USE_EXTERN_RNG) - { - arma_rng_cxx11_instance.randg_fill(out.memptr(), out.n_elem, a, b); - } - #else - { - arma_rng_cxx11 local_arma_rng_cxx11_instance; - - typedef typename arma_rng_cxx11::seed_type seed_type; - - local_arma_rng_cxx11_instance.set_seed( seed_type(arma_rng::randi()) ); - - local_arma_rng_cxx11_instance.randg_fill(out.memptr(), out.n_elem, a, b); - } - #endif + arma_rng::randg::fill(out.memptr(), out.n_elem, a, b); return out; } @@ -107,14 +95,10 @@ randg(const uword n_elem, const distr_param& param = distr_param(), const arma_e arma_ignore(junk1); arma_ignore(junk2); - if(is_Row::value) - { - return randg(1, n_elem, param); - } - else - { - return randg(n_elem, 1, param); - } + const uword n_rows = (is_Row::value) ? uword(1) : n_elem; + const uword n_cols = (is_Row::value) ? n_elem : uword(1); + + return randg(n_rows, n_cols, param); } @@ -187,6 +171,8 @@ randg(const uword n_rows, const uword n_cols, const uword n_slices, const distr_ arma_extra_debug_sigprint(); arma_ignore(junk); + typedef typename cube_type::elem_type eT; + cube_type out(n_rows, n_cols, n_slices); double a; @@ -211,21 +197,7 @@ randg(const uword n_rows, const uword n_cols, const uword n_slices, const distr_ arma_debug_check( ((a <= double(0)) || (b <= double(0))), "randg(): a and b must be greater than zero" ); - #if defined(ARMA_USE_EXTERN_RNG) - { - arma_rng_cxx11_instance.randg_fill(out.memptr(), out.n_elem, a, b); - } - #else - { - arma_rng_cxx11 local_arma_rng_cxx11_instance; - - typedef typename arma_rng_cxx11::seed_type seed_type; - - local_arma_rng_cxx11_instance.set_seed( seed_type(arma_rng::randi()) ); - - local_arma_rng_cxx11_instance.randg_fill(out.memptr(), out.n_elem, a, b); - } - #endif + arma_rng::randg::fill(out.memptr(), out.n_elem, a, b); return out; } diff --git a/include/armadillo_bits/fn_randn.hpp b/include/armadillo_bits/fn_randn.hpp index cb9dab6b..369133a0 100644 --- a/include/armadillo_bits/fn_randn.hpp +++ b/include/armadillo_bits/fn_randn.hpp @@ -63,14 +63,10 @@ randn(const uword n_elem, const arma_empty_class junk1 = arma_empty_class(), con arma_ignore(junk1); arma_ignore(junk2); - if(is_Row::value) - { - return Gen(1, n_elem); - } - else - { - return Gen(n_elem, 1); - } + const uword n_rows = (is_Row::value) ? uword(1) : n_elem; + const uword n_cols = (is_Row::value) ? n_elem : uword(1); + + return Gen(n_rows, n_cols); } diff --git a/include/armadillo_bits/fn_randu.hpp b/include/armadillo_bits/fn_randu.hpp index ca4084d6..93c31451 100644 --- a/include/armadillo_bits/fn_randu.hpp +++ b/include/armadillo_bits/fn_randu.hpp @@ -63,14 +63,10 @@ randu(const uword n_elem, const arma_empty_class junk1 = arma_empty_class(), con arma_ignore(junk1); arma_ignore(junk2); - if(is_Row::value) - { - return Gen(1, n_elem); - } - else - { - return Gen(n_elem, 1); - } + const uword n_rows = (is_Row::value) ? uword(1) : n_elem; + const uword n_cols = (is_Row::value) ? n_elem : uword(1); + + return Gen(n_rows, n_cols); } diff --git a/src/wrapper1.cpp b/src/wrapper1.cpp index c3cb806b..8a3fffe9 100644 --- a/src/wrapper1.cpp +++ b/src/wrapper1.cpp @@ -34,8 +34,13 @@ namespace arma { + // NOTE: arma_rng_cxx11_instance is kept only for compatibility with earlier versions of armadillo + // TODO: remove arma_rng_cxx11_instance when the major version is bumped + #include "armadillo_bits/arma_rng_cxx11.hpp" - thread_local arma_rng_cxx11 arma_rng_cxx11_instance; + thread_local arma_rng_cxx11 arma_rng_cxx11_instance; + + thread_local std::mt19937_64 mt19937_64_instance; } #endif