From 1353157cb6733ac43d972bcb1bbaa2aa7080f320 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ad=C3=A8le=20Harrissart?= Date: Sun, 28 Sep 2014 23:08:53 +0200 Subject: [PATCH] Using library in order to re-implement the methods of the class eoRng (issue #24). Add two new CMake Cache Values: ENABLE_CXX11_RANDOM & ENABLE_64_BIT_RNG_NUMBERS. Add a test file; Python chi-square & shapiro-wilk tests: done. --- CMakeLists.txt | 3 +- build_gcc_linux_install | 40 ++- cmake/Config.cmake | 29 +- src/eo/config.h.cmake | 6 + src/eo/pyeo/random_numbers.cpp | 10 +- src/eo/utils/eoRNG.cpp | 10 +- src/eo/utils/eoRNG.h | 600 ++++++++++++++++++++------------- test/eo/CMakeLists.txt | 8 +- test/eo/t-eoRNGcxx11.cpp | 292 ++++++++++++++++ 9 files changed, 744 insertions(+), 254 deletions(-) create mode 100644 test/eo/t-eoRNGcxx11.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index a92068c4f..a031a5230 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -109,11 +109,10 @@ if(ENABLE_CMAKE_TESTING) add_subdirectory(test) endif(ENABLE_CMAKE_TESTING) - ###################################################################################### ### 6) Packaging : only in release ! ###################################################################################### if("${CMAKE_BUILD_TYPE}" STREQUAL "Release") include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/Package.cmake) -endif() +endif() \ No newline at end of file diff --git a/build_gcc_linux_install b/build_gcc_linux_install index 8c8d308b5..64925ef2e 100755 --- a/build_gcc_linux_install +++ b/build_gcc_linux_install @@ -14,7 +14,7 @@ color() { mkdir -pv build cd build -color 5 '\nDefine CMAKE cache values.\n' +color 5 'Define CMAKE cache values.' ################################## # choose modules for compilation # @@ -121,16 +121,38 @@ else FLAG=$FLAG' '" -DENABLE_OPENMP=false" fi #ENABLE_OPENMP -color 5 "\nset $FLAG\n" +color 1 'Do you want to enable C++11 random numbers (no by default)?(y/n)' +read res -cmake $FLAG .. +if [ "$res" = "y" ]; then + FLAG=$FLAG' '" -DENABLE_CXX11_RANDOM=true" + + color 1 'Do you want to enable 64-bit random numbers (no by default)?(y/n)' + read res + + if [ "$res" = "y" ]; then + FLAG=$FLAG' '" -DENABLE_64_BIT_RNG_NUMBERS=true" + else + FLAG=$FLAG' '" -DENABLE_64_BIT_RNG_NUMBERS=false" + fi #ENABLE_64_BIT_RNG_NUMBERS + +else + FLAG=$FLAG' '" -DENABLE_CXX11_RANDOM=false" +fi #ENABLE_CXX11_RANDOM + +color 5 "set $FLAG" + +FLAG2="-DENABLE_CMAKE_TESTING=false" +FLAG2=$FLAG2' '" -DENABLE_CMAKE_EXAMPLE=false" + +cmake $FLAG $FLAG2 .. ############### # compile doc # # default: no # ############### -color 1 '\nDo you want to compile the documentation (be patient it can be a bit long)?(y/n)' +color 1 'Do you want to compile the documentation (be patient it can be a bit long)?(y/n)' read res if [ "$res" = "y" ]; then @@ -143,14 +165,14 @@ fi # default path: /usr/local/include/paradiseo # ############################################## -color 1 '\nDo you want to install the sources on your laptop?(y/n) +color 1 'Do you want to install the sources on your laptop?(y/n) If "no" you will not be able to compile tests or tutorials which require the project library installation on your laptop.' read res if [ "$res" = "y" ]; then - color 5 '\nOk. For your information we let you the possibility + color 5 'Ok. For your information we let you the possibility to change the install path on your laptop but also to personalize the project name in order to avoid some name conflicts with headers. @@ -176,7 +198,7 @@ Go see the script line 166-170 to see what it looks like.' # default: no # ################# - color 1 '\nDo you want to compile tests?(y/n)' + color 1 'Do you want to compile tests?(y/n)' read test_compile if [ "$test_compile" = "y" ]; then @@ -207,7 +229,7 @@ Go see the script line 166-170 to see what it looks like.' fi - color 5 "\nset $FLAG\n" + color 5 "set $FLAG" cmake $FLAG .. if [ "$test_compile" = "y" ] || [ "$tutorial_compile" = "y" ]; then @@ -248,3 +270,5 @@ else make fi + +cd .. diff --git a/cmake/Config.cmake b/cmake/Config.cmake index a1a709c6e..278298caf 100755 --- a/cmake/Config.cmake +++ b/cmake/Config.cmake @@ -54,7 +54,6 @@ elseif(INSTALL_TYPE STREQUAL min OR NOT DEFINED INSTALL_TYPE) set(ENABLE_CMAKE_TESTING "false" CACHE BOOL "ParadisEO tests") endif() - ###################################################################################### ### 2) Define profiling flags ###################################################################################### @@ -75,11 +74,37 @@ if(ENABLE_CMAKE_TESTING) endif(ENABLE_CMAKE_TESTING) ###################################################################################### -### 5) Build examples ? +### 4) Build examples ? ###################################################################################### set(ENABLE_CMAKE_EXAMPLE "true" CACHE BOOL "ParadisEO examples") +###################################################################################### +### 5) Random numbers +###################################################################################### + +set(ENABLE_CXX11_RANDOM "false" CACHE BOOL "For C++11 random numbers") +set(ENABLE_64_BIT_RNG_NUMBERS "false" CACHE BOOL "For 64-bit random numbers") + +# For C++11 random numbers +if(ENABLE_CXX11_RANDOM) + + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DHAVE_RANDOM") + + if (UNIX OR CMAKE_COMPILER_IS_GNUCXX) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=gnu++11") + endif (UNIX OR CMAKE_COMPILER_IS_GNUCXX) + if (APPLE) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -mmacosx-version-min=10.7") + endif (APPLE) + + # For 64-bit random numbers + if(ENABLE_64_BIT_RNG_NUMBERS) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DWITH_64_BIT_RNG_NUMBERS") + endif(ENABLE_64_BIT_RNG_NUMBERS) + +endif(ENABLE_CXX11_RANDOM) + ###################################################################################### ### 6) Determine prefix for installation ###################################################################################### diff --git a/src/eo/config.h.cmake b/src/eo/config.h.cmake index 24568be10..d2055bcc0 100755 --- a/src/eo/config.h.cmake +++ b/src/eo/config.h.cmake @@ -6,6 +6,12 @@ /* gnuplot graphical display */ #cmakedefine HAVE_GNUPLOT +/* Enable C++11 random library */ +#cmakedefine HAVE_RANDOM + +/* Enable 64-bit random numbers using C++11 random library */ +#cmakedefine WITH_64_BIT_RNG_NUMBERS + /* Define to 1 if you have the header file. */ #cmakedefine HAVE_INTTYPES_H diff --git a/src/eo/pyeo/random_numbers.cpp b/src/eo/pyeo/random_numbers.cpp index 01bb883b4..996d9f9c6 100755 --- a/src/eo/pyeo/random_numbers.cpp +++ b/src/eo/pyeo/random_numbers.cpp @@ -28,10 +28,12 @@ using namespace boost::python; using namespace boost::python; -// initialize static constants -const uint32_t eoRng::K(0x9908B0DFU); -const int eoRng::M(397); -const int eoRng::N(624); +#ifndef HAVE_RANDOM + // initialize static constants + const uint32_t eoRng::K(0x9908B0DFU); + const int eoRng::M(397); + const int eoRng::N(624); +#endif HAVE_RANDOM namespace eo { diff --git a/src/eo/utils/eoRNG.cpp b/src/eo/utils/eoRNG.cpp index b538499f1..2a9fd5f4c 100755 --- a/src/eo/utils/eoRNG.cpp +++ b/src/eo/utils/eoRNG.cpp @@ -6,10 +6,12 @@ #include #include "eoRNG.h" -// initialize static constants -const uint32_t eoRng::K(0x9908B0DFU); -const int eoRng::M(397); -const int eoRng::N(624); +#ifndef HAVE_RANDOM + // initialize static constants + const uint32_t eoRng::K(0x9908B0DFU); + const int eoRng::M(397); + const int eoRng::N(624); +#endif namespace eo { diff --git a/src/eo/utils/eoRNG.h b/src/eo/utils/eoRNG.h index 6e9af74f3..02eda3efc 100755 --- a/src/eo/utils/eoRNG.h +++ b/src/eo/utils/eoRNG.h @@ -27,9 +27,9 @@ Old contact information: todos@geneura.ugr.es, http://geneura.ugr.es /** @addtogroup Random * @{ - * */ + * */ -# if (defined _MSC_VER) +#if (defined _MSC_VER) /** uint32_t is an unsigned integer type capable of holding 32 bits. In the applicatione here exactly 32 are used. @@ -37,22 +37,48 @@ Old contact information: todos@geneura.ugr.es, http://geneura.ugr.es optimization levels so feel free to try your options and see what's best for you. */ -typedef unsigned long uint32_t; + typedef unsigned long uint32_t; #else -#if (! defined __sun) -// The C99-standard defines uint32_t to be declared in stdint.h, but some -// systems don't have that and implement it in inttypes.h. -#include -#else -#include -#endif + #if (! defined __sun) + // The C99-standard defines uint32_t to be declared in stdint.h, but some + // systems don't have that and implement it in inttypes.h. + #include + #else + #include + #endif #endif #include #include +#include // std::accumulate (see roulette_wheel()) + #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) + #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 ! + #else + typedef uint32_t uint_t; + #endif + #else + typedef uint32_t uint_t; + #endif +#else + typedef uint32_t uint_t; +#endif + /** Random Number Generator @@ -122,18 +148,25 @@ public : @see reseed for details on usage of the seeding value. */ - eoRng(uint32_t s) - : state(0), next(0), left(-1), cached(false) - { - state = new uint32_t[N+1]; - initialize(2*s); - } + eoRng(uint_t s) + #ifdef HAVE_RANDOM + : seed(s), position(0) + { + generator.seed(s); // initialize the internal state value + } + #else + : state(0), next(0), left(-1), cached(false) + { + state = new uint_t[N+1]; + initialize(2*s); + } - /** Destructor */ - ~eoRng() - { - delete [] state; - } + /** Destructor */ + ~eoRng() + { + delete [] state; + } + #endif /** Re-initializes the Random Number Generator. @@ -143,11 +176,17 @@ public : Manually divide the seed by 2 if you want to re-run old runs - @version MS. 5 Oct. 2001 + @version MS. 5 Oct. 2001eoRNG */ - void reseed(uint32_t s) + void reseed(uint_t s) { - initialize(2*s); + #ifdef HAVE_RANDOM + seed = s; + position = 0; + generator.seed(s); // re-initialize the internal state value + #else + initialize(2*s); + #endif } /* FIXME remove in next release @@ -174,7 +213,11 @@ public : */ double uniform(double m = 1.0) { // random number between [0, m] - return m * double(rand()) / double(1.0 + rand_max()); + #ifdef HAVE_RANDOM + return uniform(0.0, m); + #else + return m * double(rand()) / double(1.0 + rand_max()); + #endif } /** Random number from unifom distribution @@ -185,7 +228,13 @@ public : */ double uniform(double min, double max) { // random number between [min, max] - return min + uniform(max - min); + #ifdef HAVE_RANDOM + increment_position(); + std::uniform_real_distribution distribution(min, max); + return distribution(generator); + #else + return min + uniform(max - min); + #endif } /** Random integer number from unifom distribution @@ -193,14 +242,14 @@ public : @param m Define interval for random number to [0, m) @return random integer in the range [0, m) */ - uint32_t random(uint32_t m) + uint_t random(uint_t m) { // C++ Standard (4.9 Floatingintegral conversions [conv.fpint]) // defines floating point to integer conversion as truncation // ("rounding towards zero"): "An rvalue of a floating point type // can be converted to an rvalue of an integer type. The conversion // truncates; that is, the fractional part is discarded" - return uint32_t(uniform() * double(m)); + return uint_t(uniform() * double(m)); } /** Biased coin toss @@ -232,7 +281,13 @@ public : @return Random Gaussian deviate */ double normal(double stdev) - { return stdev * normal(); } + { + #ifdef HAVE_RANDOM + return normal(0.0, stdev); + #else + return stdev * normal(); + #endif + } /** Gaussian deviate @@ -243,19 +298,21 @@ public : @return Random Gaussian deviate */ double normal(double mean, double stdev) - { return mean + normal(stdev); } - - /** - * @brief Forgets the last cached value of normal(), so as to be able to perform some repeatable calls to normal(). - * - * As normal() stores a cached value for performance purposes, sequences of pseudo random numbers can't be repeated - * when reseeding, since the cached value can be yield before a number is generated. To avoid that, this method - * allows one to clean the cache and force to regenerate a new pseudo random number. - */ - void clearCache() - { - cached = false; - } + { + #ifdef HAVE_RANDOM + // Don't increment the position here ... + 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; + #else + return mean + normal(stdev); + #endif + } /** Random numbers using a negative exponential distribution @@ -264,19 +321,15 @@ public : */ double negexp(double mean) { - return -mean*log(double(rand()) / rand_max()); + #ifdef HAVE_RANDOM + increment_position(); + std::exponential_distribution distribution(mean); + return distribution(generator); + #else + return -mean*log(double(rand()) / rand_max()); + #endif } - /** - rand() returns a random number in the range [0, rand_max) - */ - uint32_t rand(); - - /** - rand_max() the maximum returned by rand() - */ - uint32_t rand_max() const { return uint32_t(0xffffffff); } - /** Roulette wheel selection roulette_wheel(vec, total = 0) does a roulette wheel selection @@ -287,10 +340,8 @@ public : int roulette_wheel(const std::vector& vec, TYPE total = 0) { if (total == 0) - { // count - for (unsigned i = 0; i < vec.size(); ++i) - total += vec[i]; - } + total = std::accumulate(vec.begin(), vec.end(), total); // count + double fortune = uniform() * total; int i = 0; while (fortune >= 0) @@ -300,7 +351,6 @@ public : return --i; } - /** Randomly select element from vector. @return Uniformly chosen element from the vector. @@ -330,121 +380,98 @@ public : */ void printOn(std::ostream& _os) const { - for (int i = 0; i < N; ++i) - { - _os << state[i] << ' '; - } - _os << int(next - state) << ' '; - _os << left << ' ' << cached << ' ' << cacheValue; + #ifdef HAVE_RANDOM + // print the two requested parameters on the stream + _os << seed << ' '; + _os << position << ' '; + #else + for (int i = 0; i < N; ++i) + { + _os << state[i] << ' '; + } + _os << int(next - state) << ' '; + _os << left << ' ' << cached << ' ' << cacheValue; + #endif } - /** @brief Read RNG - - @param _is Stream to read RNG from + /** @brief Read RNG. + If the C++11 random library is used the complexity + will be "linear in the number of equivalent advances". + @param _is Stream to read RNG from */ - void readFrom(std::istream& _is) + void readFrom(std::istream& _is) { - for (int i = 0; i < N; ++i) - { - _is >> state[i]; - } + #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 + #else + for (int i = 0; i < N; ++i) + { + _is >> state[i]; + } + int n; + _is >> n; + next = state + n; - int n; - _is >> n; - next = state + n; - - _is >> left; - _is >> cached; - _is >> cacheValue; + _is >> left; + _is >> cached; + _is >> cacheValue; + #endif } + /** + rand() returns a random number in the range [0, rand_max) + */ + uint_t rand(); + + /** + rand_max() the maximum returned by rand() + */ + uint_t rand_max() const + { + #ifdef HAVE_RANDOM + return generator.max(); + #else + return uint_t(0xffffffff); + #endif + } + + #ifdef HAVE_RANDOM + + /** + * @brief Increment the position of the last generated number into the state sequence + */ + void increment_position() + { + if (position >= generator.state_size) + position = 1; + else + ++position; + } + + #else + /** + * @brief Forgets the last cached value of normal(), so as to be able to perform some repeatable calls to normal(). + * + * As normal() stores a cached value for performance purposes, sequences of pseudo random numbers can't be repeated + * when reseeding, since the cached value can be yield before a number is generated. To avoid that, this method + * allows one to clean the cache and force to regenerate a new pseudo random number. + */ + void clearCache() + { + cached = false; + } + + #endif + std::string className() const { return "Mersenne-Twister"; } private: - uint32_t restart(); - - /* @brief Initialize state - - We initialize state[0..(N-1)] via the generator - - x_new = (69069 * x_old) mod 2^32 - - from Line 15 of Table 1, p. 106, Sec. 3.3.4 of Knuth's _The Art of Computer - Programming_, Volume 2, 3rd ed. - - Notes (SJC): I do not know what the initial state requirements of the - Mersenne Twister are, but it seems this seeding generator could be better. - It achieves the maximum period for its modulus (2^30) iff x_initial is odd - (p. 20-21, Sec. 3.2.1.2, Knuth); if x_initial can be even, you have - sequences like 0, 0, 0, ...; 2^31, 2^31, 2^31, ...; 2^30, 2^30, 2^30, ...; - 2^29, 2^29 + 2^31, 2^29, 2^29 + 2^31, ..., etc. so I force seed to be odd - below. - - Even if x_initial is odd, if x_initial is 1 mod 4 then - - the lowest bit of x is always 1, - the next-to-lowest bit of x is always 0, - the 2nd-from-lowest bit of x alternates ... 0 1 0 1 0 1 0 1 ... , - the 3rd-from-lowest bit of x 4-cycles ... 0 1 1 0 0 1 1 0 ... , - the 4th-from-lowest bit of x has the 8-cycle ... 0 0 0 1 1 1 1 0 ... , - ... - - and if x_initial is 3 mod 4 then - - the lowest bit of x is always 1, - the next-to-lowest bit of x is always 1, - the 2nd-from-lowest bit of x alternates ... 0 1 0 1 0 1 0 1 ... , - the 3rd-from-lowest bit of x 4-cycles ... 0 0 1 1 0 0 1 1 ... , - the 4th-from-lowest bit of x has the 8-cycle ... 0 0 1 1 1 1 0 0 ... , - ... - - The generator's potency (min. s>=0 with (69069-1)^s = 0 mod 2^32) is 16, - which seems to be alright by p. 25, Sec. 3.2.1.3 of Knuth. It also does well - in the dimension 2..5 spectral tests, but it could be better in dimension 6 - (Line 15, Table 1, p. 106, Sec. 3.3.4, Knuth). - - Note that the random number user does not see the values generated here - directly since restart() will always munge them first, so maybe none of all - of this matters. In fact, the seed values made here could even be - extra-special desirable if the Mersenne Twister theory says so-- that's why - the only change I made is to restrict to odd seeds. - */ - void initialize(uint32_t seed); - - /** @brief Array for the state */ - uint32_t *state; - - /** Pointer to next available random number */ - uint32_t *next; - - /** Number of random numbers currently left */ - int left; - - /** @brief Is there a valid cached value for the normal distribution? */ - bool cached; - - /** @brief Cached value for normal distribution? */ - double cacheValue; - - /** @brief Size of the state-vector */ - static const int N; - - /** Internal constant */ - static const int M; - - /** @brief Magic constant */ - static const uint32_t K; - - - /** @brief Copy constructor - - Private copy ctor and assignment operator to make sure that nobody - accidentally copies the random number generator. If you want similar RNG's, - make two RNG's and initialize them with the same seed. - - As it cannot be called, we do not provide an implementation. - */ eoRng(const eoRng&); /** @brief Assignment operator @@ -452,12 +479,116 @@ private: @see Copy constructor eoRng(const eoRng&). */ eoRng& operator=(const eoRng&); + + #ifdef HAVE_RANDOM // If set to true (see CMake cache values) + #ifndef ENV32BIT // Only on 64-bit environment + #ifdef WITH_64_BIT_RNG_NUMBERS // If set to true (see CMake cache values) + std::mt19937_64 generator; // Mersenne Twister pseudo-random generator of 64-bits numbers with a state of 19937 bits + // 312 elements in the state sequence + #else + std::mt19937 generator; + #endif + #else + std::mt19937 generator; // Mersenne Twister pseudo-random generator of 32-bits numbers with a state of 19937 bits + // 624 elements in the state sequence + #endif + + uint_t seed; // Last value which was used to initialize the internal state value of the generator + // (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) + + #else + + uint_t restart(); + + /* @brief Initialize state + + We initialize state[0..(N-1)] via the generator + + x_new = (69069 * x_old) mod 2^32 + + from Line 15 of Table 1, p. 106, Sec. 3.3.4 of Knuth's _The Art of Computer + Programming_, Volume 2, 3rd ed. + + Notes (SJC): I do not know what the initial state requirements of the + Mersenne Twister are, but it seems this seeding generator could be better. + It achieves the maximum period for its modulus (2^30) iff x_initial is odd + (p. 20-21, Sec. 3.2.1.2, Knuth); if x_initial can be even, you have + sequences like 0, 0, 0, ...; 2^31, 2^31, 2^31, ...; 2^30, 2^30, 2^30, ...; + 2^29, 2^29 + 2^31, 2^29, 2^29 + 2^31, ..., etc. so I force seed to be odd + below. + + Even if x_initial is odd, if x_initial is 1 mod 4 then + + the lowest bit of x is always 1, + the next-to-lowest bit of x is always 0, + the 2nd-from-lowest bit of x alternates ... 0 1 0 1 0 1 0 1 ... , + the 3rd-from-lowest bit of x 4-cycles ... 0 1 1 0 0 1 1 0 ... , + the 4th-from-lowest bit of x has the 8-cycle ... 0 0 0 1 1 1 1 0 ... , + ... + + and if x_initial is 3 mod 4 then + + the lowest bit of x is always 1, + the next-to-lowest bit of x is always 1, + the 2nd-from-lowest bit of x alternates ... 0 1 0 1 0 1 0 1 ... , + the 3rd-from-lowest bit of x 4-cycles ... 0 0 1 1 0 0 1 1 ... , + the 4th-from-lowest bit of x has the 8-cycle ... 0 0 1 1 1 1 0 0 ... , + ... + + The generator's potency (min. s>=0 with (69069-1)^s = 0 mod 2^32) is 16, + which seems to be alright by p. 25, Sec. 3.2.1.3 of Knuth. It also does well + in the dimension 2..5 spectral tests, but it could be better in dimension 6 + (Line 15, Table 1, p. 106, Sec. 3.3.4, Knuth). + + Note that the random number user does not see the values generated here + directly since restart() will always munge them first, so maybe none of all + of this matters. In fact, the seed values made here could even be + extra-special desirable if the Mersenne Twister theory says so-- that's why + the only change I made is to restrict to odd seeds. + */ + void initialize(uint_t seed); + + /** @brief Array for the state */ + uint_t *state; + + /** Pointer to next available random number */ + uint_t *next; + + /** Number of random numbers currently left */ + int left; + + /** @brief Is there a valid cached value for the normal distribution? */ + bool cached; + + /** @brief Cached value for normal distribution? */ + double cacheValue; + + /** @brief Size of the state-vector */ + static const int N; + + /** Internal constant */ + static const int M; + + /** @brief Magic constant */ + static const uint_t K; + + + /** @brief Copy constructor + + Private copy ctor and assignment operator to make sure that nobody + accidentally copies the random number generator. If you want similar RNG's, + make two RNG's and initialize them with the same seed. + + As it cannot be called, we do not provide an implementation. + */ + #endif }; /** @example t-eoRNG.cpp */ - - namespace eo { /** The one and only global eoRng object */ @@ -468,83 +599,87 @@ using eo::rng; /** @} */ +#ifndef HAVE_RANDOM + + // Implementation of some eoRng members.... Don't mind the mess, it does work. + + #define hiBit(u) ((u) & 0x80000000U) // mask all but highest bit of u + #define loBit(u) ((u) & 0x00000001U) // mask all but lowest bit of u + #define loBits(u) ((u) & 0x7FFFFFFFU) // mask the highest bit of u + #define mixBits(u, v) (hiBit(u)|loBits(v)) // move hi bit of u to hi bit of v -// Implementation of some eoRng members.... Don't mind the mess, it does work. + inline void eoRng::initialize(uint_t seed) + { + left = -1; -#define hiBit(u) ((u) & 0x80000000U) // mask all but highest bit of u -#define loBit(u) ((u) & 0x00000001U) // mask all but lowest bit of u -#define loBits(u) ((u) & 0x7FFFFFFFU) // mask the highest bit of u -#define mixBits(u, v) (hiBit(u)|loBits(v)) // move hi bit of u to hi bit of v + register uint_t x = (seed | 1U) & 0xFFFFFFFFU, *s = state; + register int j; + + for(left=0, *s++=x, j=N; --j; + *s++ = (x*=69069U) & 0xFFFFFFFFU) ; + } -inline void eoRng::initialize(uint32_t seed) + inline uint_t eoRng::restart() + { + register uint_t *p0=state, *p2=state+2, *pM=state+M, s0, s1; + register int j; + + left=N-1, next=state+1; + + for(s0=state[0], s1=state[1], j=N-M+1; --j; s0=s1, s1=*p2++) + *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U); + + for(pM=state, j=M; --j; s0=s1, s1=*p2++) + *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U); + + s1=state[0], *p0 = *pM ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U); + s1 ^= (s1 >> 11); + s1 ^= (s1 << 7) & 0x9D2C5680U; + s1 ^= (s1 << 15) & 0xEFC60000U; + return(s1 ^ (s1 >> 18)); + } + +#endif + +inline uint_t eoRng::rand() { - left = -1; - - register uint32_t x = (seed | 1U) & 0xFFFFFFFFU, *s = state; - register int j; - - for(left=0, *s++=x, j=N; --j; - *s++ = (x*=69069U) & 0xFFFFFFFFU) ; + #ifdef HAVE_RANDOM + increment_position(); + return generator(); + #else + if(--left < 0) + return(restart()); + uint_t y = *next++; + y ^= (y >> 11); + y ^= (y << 7) & 0x9D2C5680U; + y ^= (y << 15) & 0xEFC60000U; + return(y ^ (y >> 18)); + #endif } - - -inline uint32_t eoRng::restart() -{ - register uint32_t *p0=state, *p2=state+2, *pM=state+M, s0, s1; - register int j; - - left=N-1, next=state+1; - - for(s0=state[0], s1=state[1], j=N-M+1; --j; s0=s1, s1=*p2++) - *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U); - - for(pM=state, j=M; --j; s0=s1, s1=*p2++) - *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U); - - s1=state[0], *p0 = *pM ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U); - s1 ^= (s1 >> 11); - s1 ^= (s1 << 7) & 0x9D2C5680U; - s1 ^= (s1 << 15) & 0xEFC60000U; - return(s1 ^ (s1 >> 18)); -} - - - -inline uint32_t eoRng::rand() -{ - if(--left < 0) - return(restart()); - uint32_t y = *next++; - y ^= (y >> 11); - y ^= (y << 7) & 0x9D2C5680U; - y ^= (y << 15) & 0xEFC60000U; - return(y ^ (y >> 18)); -} - - - inline double eoRng::normal() -{ - if (cached) { - cached = false; - return cacheValue; + { + #ifdef HAVE_RANDOM + return eoRng::normal(0.0, 1.0); + #else + if (cached) { + cached = false; + return cacheValue; + } + double rSquare, var1, var2; + do { + var1 = 2.0 * uniform() - 1.0; + var2 = 2.0 * uniform() - 1.0; + rSquare = var1 * var1 + var2 * var2; + } while (rSquare >= 1.0 || rSquare == 0.0); + double factor = sqrt(-2.0 * log(rSquare) / rSquare); + cacheValue = var1 * factor; + cached = true; + return (var2 * factor); + #endif } - double rSquare, var1, var2; - do { - var1 = 2.0 * uniform() - 1.0; - var2 = 2.0 * uniform() - 1.0; - rSquare = var1 * var1 + var2 * var2; - } while (rSquare >= 1.0 || rSquare == 0.0); - double factor = sqrt(-2.0 * log(rSquare) / rSquare); - cacheValue = var1 * factor; - cached = true; - return (var2 * factor); -} - - namespace eo { @@ -599,7 +734,6 @@ namespace eo inline double normal() { return rng.normal(); } } - #endif @@ -609,4 +743,4 @@ namespace eo // c-file-offsets: ((c . 0)) // c-file-style: "Stroustrup" // fill-column: 80 -// End: +// End: \ No newline at end of file diff --git a/test/eo/CMakeLists.txt b/test/eo/CMakeLists.txt index 0dff4c9ce..f4b03c406 100755 --- a/test/eo/CMakeLists.txt +++ b/test/eo/CMakeLists.txt @@ -71,12 +71,18 @@ set (TEST_LIST t-eoParser ) +# For C++11 random numbers +if(ENABLE_CXX11_RANDOM) + set (TEST_LIST ${TEST_LIST} t-eoRNGcxx11) # tests the eoRng class update (thanks to the new C++11 random library) + + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -I /usr/include/python2.7 -lpython2.7") # Install developpement package for Python and add Python library link + +endif(ENABLE_CXX11_RANDOM) foreach (test ${TEST_LIST}) set ("T_${test}_SOURCES" "${test}.cpp") endforeach (test) - if(ENABLE_MINIMAL_CMAKE_TESTING) set (MIN_TEST_LIST t-eoEasyPSO) diff --git a/test/eo/t-eoRNGcxx11.cpp b/test/eo/t-eoRNGcxx11.cpp new file mode 100644 index 000000000..91c7a16ec --- /dev/null +++ b/test/eo/t-eoRNGcxx11.cpp @@ -0,0 +1,292 @@ +#include +#include + +#include // std::stringstream +#include // assertions are used by serialization tests +#include +#include +#include // in order to reseed generator +#include // statistical tests using Python +#include // statistical tests using Python +#include + +void basic_tests() +{ + std::cout << "basic_tests" << std::endl; + + std::cout << "rand():" << rng.rand() << std::endl; + std::cout << "rand_max():" << rng.rand_max() << std::endl; + + uint_t s = static_cast(time(0)); + rng.reseed(s); + + std::cout << "uniform():" << rng.uniform() << std::endl; + std::cout << "uniform(2.4):" << rng.uniform(2.4) << std::endl; + std::cout << "uniform(0.0, 20.0):" << rng.uniform(0.0, 20.0) << std::endl; + + std::cout << "flip:" << rng.flip() << std::endl; + + std::cout << "normal():" << rng.normal() << std::endl; + std::cout << "normal(6.7):" << rng.normal(6.7) << std::endl; + std::cout << "normal(1000, 0.1):" << rng.normal(1000, 0.1) << std::endl; + + std::cout << "negexp:" << rng.negexp(4.7) << std::endl; + + int tab[] = {3, 46, 6, 89, 50, 78}; + std::vector v (tab, tab + sizeof(tab) / sizeof(int)); + std::cout << "roulette_wheel(v):" << rng.roulette_wheel(v) << std::endl; + std::cout << "choice(v):" << rng.choice(v) << std::endl; +} + +void test_function(const unsigned N, uint_t(eoRng::*ptr)()) +{ + std::stringstream ss; + ss << rng; // Print eoRNG on stream + uint_t r1 = (rng.*ptr)(); + for (unsigned i = 0; i < N; i++) + (rng.*ptr)(); + ss >> rng; // Read eoRNG from stream + uint_t r2 = (rng.*ptr)(); + + assert(r1 == r2); +} + +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); + for (unsigned i = 0; i < N; i++) + (rng.*ptr)(m); + ss >> rng; // Read eoRNG from stream + uint_t r2 = (rng.*ptr)(m); + + assert(r1 == r2); +} + +/** Test the serialization + * The serialization allows the user to stop the pseudo-random generator + * and save its state in a stream in order to continue the generation + * from this exact state later. + */ +void serialization_tests(const unsigned N, const double d) +{ + std::cout << "serialization_test with N=" << N << " and d=" << d << ":" << std::endl; + + uint_t(eoRng::*ptr_rand)() = &eoRng::rand; + double(eoRng::*ptr_uniform)(double) = &eoRng::uniform; + double(eoRng::*ptr_normal)(double) = &eoRng::normal; + double(eoRng::*ptr_negexp)(double) = &eoRng::negexp; + + uint_t s = static_cast(time(0)); + + rng.reseed(s); + test_function(N, ptr_rand); + + rng.reseed(s); + test_function(N, ptr_uniform, double(1.0)); + + rng.reseed(s); + test_function(N, ptr_normal, d); + + rng.reseed(s); + test_function(N, ptr_negexp, d); + + std::cout << "ok" << std::endl; +} + +/** Statistical tests using Python + * + * scipy.stats.shapiro: + * "Perform the Shapiro-Wilk test for normality". + * "The Shapiro-Wilk test tests the null hypothesis that the data was drawn from a normal distribution." + * See documentation for more details : http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.shapiro.html + * This test is based on the 'normal()' eoRng-implemented method. + * + * scipy.stats.chisquare: + * "Calculates a one-way chi square test." + * "The chi square test tests the null hypothesis that the categorical data has the given frequencies." + * See documentation for more details : http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chisquare.html + * This test is based on the 'uniform()' eoRng-implemented method. + */ +void stat_tests(const unsigned N) +{ + // Initialize the Python interpreter + Py_Initialize(); + // Create Python objects which will be used later + PyObject *pName, *pModule, *pDict, *pFunc, *pArgs_shapi, *pArgs_csquare, *f_obs, *f_expect, *pX, *pResult = NULL; + // Convert the module name into a Python string + pName = PyString_FromString("scipy.stats"); + // Import the Python module + pModule = PyImport_Import(pName); + // Import Numpty too + pName = PyString_FromString("numpy"); + PyImport_Import(pName); + // Create a dictionary for the contents of the module + pDict = PyModule_GetDict(pModule); + // In order to reseed generator + uint_t s; + // For Python N-arrays constructions + import_array(); + npy_intp dims[1] = {N}; + + std::cout << "Python Shapiro-Wilk test based on normal() distribution " << N << " values:" << std::endl; + + // Memory allocation requested (if N too big) + double* x = NULL; + x = (double*) malloc(N*sizeof(double)); + if (x == NULL) + { + std::cout << "Memory allocation failed" << std::endl; + exit(-1); + } + + 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); + + x[i] = rng.normal(); + } + + // Create a Python tuple to hold the method arguments + pArgs_shapi = PyTuple_New(1); + + pX = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, x); + PyTuple_SetItem(pArgs_shapi, 0, pX); + + // Get the method from the dictionary + pFunc = PyDict_GetItemString(pDict, "shapiro"); // scipy.stats.shapiro(x, a=None, reta=False) + + // Call the Python method + pResult = PyObject_CallObject(pFunc, pArgs_shapi); + + // Print a message if the method call failed + if (pResult == NULL) + { + std::cout << "Python Shapiro-Wilk test method call has failed. See shapirowilk_test()." << std::endl; + PyErr_Print(); + exit(-1); + } + else + { + double W = PyFloat_AsDouble(PyTuple_GetItem(pResult, 0)); + 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); + } + + // Free allocated memory + free(x); + + std::cout << "Python chi square test based on uniform() distribution with " << N << " values:" << std::endl; + + uint_t max = rng.rand_max(); + double v_exp = 1.0/max; + + // Memory allocation requested (if N too big) + double* expect_freq = NULL; + double* obs_freq = NULL; + + expect_freq = (double*) malloc(N*sizeof(double)); + if (expect_freq == NULL) + { + std::cout << "First memory allocation failed" << std::endl; + exit(-1); + } + obs_freq = (double*) malloc(N*sizeof(double)); + if (obs_freq == NULL) + { + std::cout << "Second memory allocation failed" << std::endl; + exit(-1); + } + + 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); + + // Observed frequencies + obs_freq[i] = rng.uniform(); + + // Expected frequencies + expect_freq[i] = v_exp; + } + + f_obs = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, obs_freq); + f_expect = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, expect_freq); + + // Create Python tuple to hold the method arguments + pArgs_csquare = PyTuple_New(2); + + PyTuple_SetItem(pArgs_csquare, 0, f_obs); + PyTuple_SetItem(pArgs_csquare, 1, f_expect); + + // Get the method from the dictionary + pFunc = PyDict_GetItemString(pDict, "chisquare"); // scipy.stats.chisquare(f_obs, f_exp=None, ddof=0, axis=0) + + // Call the Python method + pResult = PyObject_CallObject(pFunc, pArgs_csquare); + + // Print a message if the method call failed + if (pResult == NULL) // it should be two float numbers (chisq := the chi square test statistic, p := the p-value of the test) + { + std::cout << "Python chi square test method call has failed. See chisquare_test()." << std::endl; + PyErr_Print(); + } + else + { + double chisq = PyFloat_AsDouble(PyTuple_GetItem(pResult, 0)); + double p_value = PyFloat_AsDouble(PyTuple_GetItem(pResult, 1)); + std::cout << "chi square test statistic:" << chisq << std::endl; + std::cout << "p-value of the chi square test (strong if < 0.05):" << p_value << std::endl; + Py_DECREF(pResult); + } + + // Free allocated memory + free(expect_freq); + free(obs_freq); + + // Clean up + // Desallocation Python function call + // Python objects are : PyObject *pName, *pModule, *pDict, *pFunc, *pArgs_shapi, *pArgs_csquare, *f_obs, *f_expect, *pX, *pResult; + Py_DECREF(pName); + Py_DECREF(pModule); + //Py_DECREF(pDict); // Causes a 'Fatal Python error: deletion of interned string failed. Aborted' + Py_DECREF(pFunc); + Py_DECREF(pArgs_shapi); // Causes a 'Segmentation fault' error. + Py_DECREF(pArgs_csquare); + //Py_DECREF(pX); // Causes a 'Segmentation fault' error. + //Py_DECREF(f_obs); // Causes a 'Segmentation fault' error. + //Py_DECREF(f_expect); // Causes a 'Segmentation fault' error. + + // Destroy the Python interpreter + Py_Finalize(); +} + +int main() +{ + basic_tests(); + + rng.reseed(static_cast(time(0))); + constexpr double min = std::numeric_limits::min(); + constexpr double max = std::numeric_limits::max(); + double d = rng.uniform(min, max); + serialization_tests(635, d); // Reminder: + // Mersenne Twister pseudo-random generator of 32-bits numbers with a state of 19937 bits + // 624 elements in the state sequence + // Mersenne Twister pseudo-random generator of 64-bits numbers with a state of 19937 bits + // 312 elements in the state sequence + + stat_tests(2500); // Warning: limit for N (Shapiro-Wilk test) + // if N > 5000: warnings.warn("p-value may not be accurate for N > 5000.") + // See https://github.com/scipy/scipy/blob/master/scipy/stats/morestats.py for more details. + // Personal warning with API C-Python : causes valgrind errors + + return 0; +} \ No newline at end of file