nojhan/paradiseo
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Using library <random> in order to re-implement the methods of the class eoRng

Browse source 
(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.
This commit is contained in:
Adèle Harrissart 2014-09-28 23:08:53 +02:00
commit 1353157cb6
9 changed files with 716 additions and 226 deletions
Download patch file Download diff file Expand all files Collapse all files

3
CMakeLists.txt
View file

@ -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()

40
build_gcc_linux_install
View file

@ -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 ..

29
cmake/Config.cmake
View file

@ -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
######################################################################################

6
src/eo/config.h.cmake
View file

@ -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 <inttypes.h> header file. */
#cmakedefine HAVE_INTTYPES_H

10
src/eo/pyeo/random_numbers.cpp
View file

@ -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
{

10
src/eo/utils/eoRNG.cpp
View file

@ -6,10 +6,12 @@
#include <ctime>
#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
{

600
src/eo/utils/eoRNG.h
View file

@ -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 <stdint.h>
#else
#include <inttypes.h>
#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 <stdint.h>
#else
#include <inttypes.h>
#endif
#endif
#include <cmath>
#include <vector>
#include <numeric> // std::accumulate (see roulette_wheel())
#include "../eoPersistent.h"
#include "../eoObject.h"
#ifdef HAVE_RANDOM
#include <random> // 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<double> 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<double> 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<double> 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<TYPE>& 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
<tt>x_new = (69069 * x_old) mod 2^32</tt>
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
<tt>x_new = (69069 * x_old) mod 2^32</tt>
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:

8
test/eo/CMakeLists.txt
View file

@ -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)

292
test/eo/t-eoRNGcxx11.cpp Normal file
View file

@ -0,0 +1,292 @@
#include <paradiseo/eo.h>
#include <paradiseo/eo/utils/eoRNG.h>
#include <sstream> // std::stringstream
#include <assert.h> // assertions are used by serialization tests
#include <vector>
#include <array>
#include <time.h> // in order to reseed generator
#include <Python.h> // statistical tests using Python
#include <numpy/arrayobject.h> // statistical tests using Python
#include <limits>
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<uint_t>(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<int> 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<uint_t>(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<uint_t>(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<uint_t>(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<uint_t>(time(0)));
constexpr double min = std::numeric_limits<double>::min();
constexpr double max = std::numeric_limits<double>::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;
}