diff --git a/build_gcc_linux_install b/build_gcc_linux_install index 64925ef2e..f6b64353f 100755 --- a/build_gcc_linux_install +++ b/build_gcc_linux_install @@ -137,7 +137,7 @@ if [ "$res" = "y" ]; then fi #ENABLE_64_BIT_RNG_NUMBERS else - FLAG=$FLAG' '" -DENABLE_CXX11_RANDOM=false" + FLAG=$FLAG' '" -DENABLE_CXX11_RANDOM=false -DENABLE_64_BIT_RNG_NUMBERS=false" fi #ENABLE_CXX11_RANDOM color 5 "set $FLAG" diff --git a/src/eo/utils/eoRNG.h b/src/eo/utils/eoRNG.h index 02eda3efc..521633c4b 100755 --- a/src/eo/utils/eoRNG.h +++ b/src/eo/utils/eoRNG.h @@ -55,17 +55,15 @@ Old contact information: todos@geneura.ugr.es, http://geneura.ugr.es #include "../eoPersistent.h" #include "../eoObject.h" -#ifdef HAVE_RANDOM - #include // Mersenne Twister available since C++11 ! About random library : - // "this library allows to produce random numbers using combinations of generators and distributions". - // (see www.cplusplus.com/reference/random/) -#endif - #if (UINT_MAX == 0xFFFFFFFFU) // Compile time check (32-bit vs. 64-bit environment) #define ENV32BIT #endif #ifdef HAVE_RANDOM // If set to true (see CMake cache values) + #include // Mersenne Twister available since C++11 ! About random library : + // "this library allows to produce random numbers using combinations of generators and distributions". + // (see www.cplusplus.com/reference/random/) + #ifndef ENV32BIT // Only on 64-bit environment #ifdef WITH_64_BIT_RNG_NUMBERS // If set to true (see CMake cache values) typedef uint64_t uint_t; // The C++11 Mersenne Twister pseudo-random generator can generate 64-bits numbers ! @@ -150,7 +148,7 @@ public : */ eoRng(uint_t s) #ifdef HAVE_RANDOM - : seed(s), position(0) + : seed(s), position(0), discard (false) { generator.seed(s); // initialize the internal state value } @@ -181,9 +179,10 @@ public : void reseed(uint_t s) { #ifdef HAVE_RANDOM - seed = s; + seed = 2*s; position = 0; - generator.seed(s); // re-initialize the internal state value + discard = false; + generator.seed(seed); // re-initialize the internal state value #else initialize(2*s); #endif @@ -229,7 +228,7 @@ public : double uniform(double min, double max) { // random number between [min, max] #ifdef HAVE_RANDOM - increment_position(); + discard = true; std::uniform_real_distribution distribution(min, max); return distribution(generator); #else @@ -300,15 +299,9 @@ public : double normal(double mean, double stdev) { #ifdef HAVE_RANDOM - // Don't increment the position here ... + discard = true; std::normal_distribution distribution(mean, stdev); - double ret = distribution(generator); - // ... but reseed the generator and call discard() to go further into the state sequence. - generator.seed(seed); - generator.discard(position); // According to the C++11 random library documentation - // the function complexity will be linear in the - // number of equivalent advances - return ret; + return distribution(generator); #else return mean + normal(stdev); #endif @@ -322,7 +315,7 @@ public : double negexp(double mean) { #ifdef HAVE_RANDOM - increment_position(); + discard = true; std::exponential_distribution distribution(mean); return distribution(generator); #else @@ -404,10 +397,7 @@ public : #ifdef HAVE_RANDOM _is >> seed; _is >> position; - generator.seed(seed); - generator.discard(position); // According to the C++11 random library documentation - // the function complexity will be linear in the - // number of equivalent advances + discard = true; #else for (int i = 0; i < N; ++i) { @@ -448,9 +438,19 @@ public : void increment_position() { if (position >= generator.state_size) + { + generator.seed(seed); position = 1; - else - ++position; + discard = false; + return; + } + if (discard) + { + generator.seed(seed); + generator.discard(position); // linear complexity + discard = false; + } + ++position; } #else @@ -497,7 +497,9 @@ private: // (one of the two requested parameters for the serialization) unsigned position; // Position of the last generated number into the state sequence - // (one of the two requested parameters for the serialization) + // (one of the two requested parameters for the serialization) + + bool discard; // Call discard() after using a distribution #else diff --git a/test/eo/t-eoRNGcxx11.cpp b/test/eo/t-eoRNGcxx11.cpp index 91c7a16ec..07e0f5e12 100644 --- a/test/eo/t-eoRNGcxx11.cpp +++ b/test/eo/t-eoRNGcxx11.cpp @@ -41,12 +41,13 @@ void basic_tests() void test_function(const unsigned N, uint_t(eoRng::*ptr)()) { std::stringstream ss; - ss << rng; // Print eoRNG on stream - uint_t r1 = (rng.*ptr)(); + rng.rand(); for (unsigned i = 0; i < N; i++) (rng.*ptr)(); + ss << rng; // Print eoRNG on stream + uint_t r1 = rng.rand(); ss >> rng; // Read eoRNG from stream - uint_t r2 = (rng.*ptr)(); + uint_t r2 = rng.rand(); assert(r1 == r2); } @@ -54,12 +55,13 @@ void test_function(const unsigned N, uint_t(eoRng::*ptr)()) void test_function(const unsigned N, double(eoRng::*ptr)(double), const double m) { std::stringstream ss; - ss << rng; // Print eoRNG on stream - uint_t r1 = (rng.*ptr)(m); + rng.rand(); for (unsigned i = 0; i < N; i++) (rng.*ptr)(m); + ss << rng; // Print eoRNG on stream + uint_t r1 = rng.rand(); ss >> rng; // Read eoRNG from stream - uint_t r2 = (rng.*ptr)(m); + uint_t r2 = rng.rand(); assert(r1 == r2); } @@ -81,16 +83,16 @@ void serialization_tests(const unsigned N, const double d) uint_t s = static_cast(time(0)); rng.reseed(s); - test_function(N, ptr_rand); + test_function(N, ptr_rand); rng.reseed(s); - test_function(N, ptr_uniform, double(1.0)); + test_function(N, ptr_uniform, double(1.0)); rng.reseed(s); - test_function(N, ptr_normal, d); + test_function(N, ptr_normal, d); rng.reseed(s); - test_function(N, ptr_negexp, d); + test_function(N, ptr_negexp, d); std::cout << "ok" << std::endl; } @@ -145,7 +147,7 @@ void stat_tests(const unsigned N) { // Important : don't forget to reseed the generator s = static_cast(time(0) + i-1); // Chosen method to reseed generator (contestable): - // Add i-1 seconds to the current time based on the current system + // Add i-1 seconds to the current time based on the current system rng.reseed(s); x[i] = rng.normal(); @@ -176,7 +178,7 @@ void stat_tests(const unsigned N) double p_value = PyFloat_AsDouble(PyTuple_GetItem(pResult, 1)); std::cout << "Shapiro-Wilk test statistic:" << W << std::endl; std::cout << "p-value of the Shapiro-Wilk test (strong if < 0.05):" << p_value << std::endl; - //Py_DECREF(pResult); + Py_DECREF(pResult); } // Free allocated memory @@ -206,10 +208,10 @@ void stat_tests(const unsigned N) for (unsigned i = 0; i < N; i++) { - // Important : don't forget to reseed the generator - s = static_cast(time(0) + i-1); // Chosen method to reseed generator (contestable): - // Add i-1 seconds to the current time based on the current system - rng.reseed(s); + // // Important : don't forget to reseed the generator + // s = static_cast(time(0) + i-1); // Chosen method to reseed generator (contestable): + // // Add i-1 seconds to the current time based on the current system + // rng.reseed(s); // Observed frequencies obs_freq[i] = rng.uniform();