Commit 87c3ec6f authored by Ogier Maitre's avatar Ogier Maitre

Main modification : using MersenTwister random number generator.

Add some destructors to suppress warnings.
parent faabf922
......@@ -105,6 +105,9 @@ CEvolutionaryAlgorithm::CEvolutionaryAlgorithm(Parameters* params){
#endif
}
CEvolutionaryAlgorithm::~CEvolutionaryAlgorithm(){
delete population;
}
void CEvolutionaryAlgorithm::addStoppingCriterion(CStoppingCriterion* sc){
this->stoppingCriteria.push_back(sc);
}
......
......@@ -110,8 +110,13 @@ void parseArguments(const char* parametersFileName, int ac, char** av,
po::variables_map& vm, po::variables_map& vm_file){
char** argv;
int argc = loadParametersFile(parametersFileName,&argv);
int argc;
if( parametersFileName )
argc = loadParametersFile(parametersFileName,&argv);
else{
argc = 0;
argv = NULL;
}
po::options_description desc("Allowed options ");
desc.add_options()
("help", "produce help message")
......@@ -149,10 +154,12 @@ void parseArguments(const char* parametersFileName, int ac, char** av,
("u2",po::value<string>(),"User defined parameter 2")
("u3",po::value<int>(),"User defined parameter 3")
("u4",po::value<int>(),"User defined parameter 4")
("u5",po::value<int>(),"User defined parameter 5")
;
try{
po::store(po::parse_command_line(ac, av, desc,0), vm);
if( parametersFileName )
po::store(po::parse_command_line(argc, argv, desc,0), vm_file);
}
catch(po::unknown_option& e){
......@@ -171,6 +178,7 @@ void parseArguments(const char* parametersFileName, int ac, char** av,
for( int i = 0 ; i<argc ; i++ )
free(argv[i]);
if( argv )
free(argv);
}
......
......@@ -418,6 +418,9 @@ void CPopulation::weakElitism(size_t elitismSize, CIndividual** parentsPopulatio
void CPopulation::addIndividualParentPopulation(CIndividual* indiv, size_t id){
parents[id] = indiv;
}
void CPopulation::addIndividualParentPopulation(CIndividual* indiv){
parents[actualParentPopulationSize++] = indiv;
}
std::ostream& operator << (std::ostream& O, const CPopulation& B)
{
......
......@@ -7,62 +7,92 @@
#include "include/CRandomGenerator.h"
#include "include/global.h"
#include <stdlib.h>
//#include <stdlib.h>
#include <stdio.h>
#include <math.h>
CRandomGenerator::CRandomGenerator(unsigned int seed){
srand(seed);
this->seed = seed;
this->mt_rnd = new MTRand(seed);
//srand(seed);
}
int CRandomGenerator::randInt(){
return rand();
//return rand();
return mt_rnd->randInt();
}
bool CRandomGenerator::tossCoin(){
int rVal = rand();
if( rVal >=(RNDMAX/2))
return true;
else return false;
int rVal = mt_rnd->randInt(1);
return (rVal==1?true:false);
}
bool CRandomGenerator::tossCoin(float bias){
int rVal = rand();
if( rVal <=(RNDMAX*bias) )
double rVal = mt_rnd->rand(1.);
if( rVal <=bias )
return true;
else return false;
}
int CRandomGenerator::randInt(int min, int max){
int rValue = (((float)rand()/RNDMAX))*(max-min);
//DEBUG_PRT("Int Random Value : %d",min+rValue);
return rValue+min;
max--; // exclude upper bound
return min+mt_rnd->randInt(max-min);
}
int CRandomGenerator::random(int min, int max){
return randInt(min,max);
return this->randInt(min,max);
}
float CRandomGenerator::randFloat(float min, float max){
float rValue = (((float)rand()/RNDMAX))*(max-min);
//DEBUG_PRT("Float Random Value : %f",min+rValue);
return rValue+min;
return min+mt_rnd->randExc(max-min);
}
float CRandomGenerator::random(float min, float max){
return randFloat(min,max);
return this->randFloat(min,max);
}
double CRandomGenerator::random(double min, double max){
return randFloat(min,max);
return this->randFloat(min,max);
}
int CRandomGenerator::getRandomIntMax(int max){
double r = rand();
r = r / RNDMAX;
r = r * max;
return r;
return this->randInt(0,max);
}
/**
Box-Muller method for gaussian random distribution
Not sure, this function is really working.
*/
void CRandomGenerator::random_gauss(float mean, float std_dev, float* z_0, float* z_1){
float x1, x2, w;
do {
x1 = 2.0 * this->random(0.,1.) - 1.0;
x2 = 2.0 * this->random(0.,1.) - 1.0;
w = x1 * x1 + x2 * x2;
} while ( w >= 1.0 );
w = sqrt( (-2.0 * log( w ) ) / w );
*z_0 = (x1 * w)*std_dev+mean;
*z_1 = (x2 * w)*std_dev+mean;
}
float CRandomGenerator::random_gauss(float mean, float std_dev){
float z_0,z_1;
this->random_gauss(mean, std_dev, &z_0,&z_1);
return (this->tossCoin(0.5)?z_0:z_1);
}
std::ostream & operator << (std::ostream & os, const CRandomGenerator& rg) {
os<< "s : " << rg.seed << std::endl;
return os;
}
......@@ -62,6 +62,7 @@ public:
float currentAverageFitness;
float currentSTDEV;
virtual ~CEvolutionaryAlgorithm();
std::vector<CStoppingCriterion*> stoppingCriteria;
......
......@@ -64,6 +64,7 @@ public:
//virtual void initializeParentPopulation() = 0;
void addIndividualParentPopulation(CIndividual* indiv, size_t id);
void addIndividualParentPopulation(CIndividual* indiv);
void evaluatePopulation(CIndividual** population, size_t populationSize);
virtual void optimisePopulation(CIndividual** population, size_t populationSize);
virtual void evaluateParentPopulation();
......@@ -100,7 +101,7 @@ public:
void sortParentPopulation(){ CPopulation::sortPopulation(parents,actualParentPopulationSize);}
void produceOffspringPopulation();
virtual void produceOffspringPopulation();
friend std::ostream& operator << (std::ostream& O, const CPopulation& B);
......
......@@ -8,7 +8,12 @@
#ifndef CRANDOMGENERATOR_H_
#define CRANDOMGENERATOR_H_
#include <iostream>
#include "MersenneTwister.h"
class CRandomGenerator {
private:
unsigned seed;
MTRand* mt_rnd;
public:
CRandomGenerator(unsigned int seed);
int randInt();
......@@ -20,6 +25,13 @@ public:
int random(int min, int max);
float random(float min, float max);
double random(double min, double max);
void random_gauss(float min, float max, float* z_0, float* z_1);
float random_gauss(float mean, float std_dev);
unsigned get_seed()const {return this->seed;}
friend std::ostream & operator << (std::ostream & os, const CRandomGenerator& rg);
};
#endif /* CRANDOMGENERATOR_H_ */
......@@ -17,6 +17,7 @@ public:
virtual void initialize(CIndividual** population, float selectionPressure, size_t populationSize);
virtual size_t selectNext(size_t populationSize);
virtual float getExtremum() = 0 ;
virtual ~CSelectionOperator(){;}
protected:
CIndividual** population;
float currentSelectionPressure;
......
......@@ -27,7 +27,7 @@ class CStoppingCriterion {
public:
virtual bool reached() = 0;
virtual ~CStoppingCriterion(){;}
};
......
// MersenneTwister.h
// Mersenne Twister random number generator -- a C++ class MTRand
// Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
// Richard J. Wagner v1.1 28 September 2009 wagnerr@umich.edu
// The Mersenne Twister is an algorithm for generating random numbers. It
// was designed with consideration of the flaws in various other generators.
// The period, 2^19937-1, and the order of equidistribution, 623 dimensions,
// are far greater. The generator is also fast; it avoids multiplication and
// division, and it benefits from caches and pipelines. For more information
// see the inventors' web page at
// http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
// Reference
// M. Matsumoto and T. Nishimura, "Mersenne Twister: A 623-Dimensionally
// Equidistributed Uniform Pseudo-Random Number Generator", ACM Transactions on
// Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp 3-30.
// Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
// Copyright (C) 2000 - 2009, Richard J. Wagner
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions
// are met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. The names of its contributors may not be used to endorse or promote
// products derived from this software without specific prior written
// permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
// The original code included the following notice:
//
// When you use this, send an email to: m-mat@math.sci.hiroshima-u.ac.jp
// with an appropriate reference to your work.
//
// It would be nice to CC: wagnerr@umich.edu and Cokus@math.washington.edu
// when you write.
#ifndef MERSENNETWISTER_H
#define MERSENNETWISTER_H
// Not thread safe (unless auto-initialization is avoided and each thread has
// its own MTRand object)
#include <iostream>
#include <climits>
#include <cstdio>
#include <ctime>
#include <cmath>
class MTRand {
// Data
public:
typedef unsigned long uint32; // unsigned integer type, at least 32 bits
enum { N = 624 }; // length of state vector
enum { SAVE = N + 1 }; // length of array for save()
protected:
enum { M = 397 }; // period parameter
uint32 state[N]; // internal state
uint32 *pNext; // next value to get from state
int left; // number of values left before reload needed
// Methods
public:
MTRand( const uint32 oneSeed ); // initialize with a simple uint32
MTRand( uint32 *const bigSeed, uint32 const seedLength = N ); // or array
MTRand(); // auto-initialize with /dev/urandom or time() and clock()
MTRand( const MTRand& o ); // copy
// Do NOT use for CRYPTOGRAPHY without securely hashing several returned
// values together, otherwise the generator state can be learned after
// reading 624 consecutive values.
// Access to 32-bit random numbers
uint32 randInt(); // integer in [0,2^32-1]
uint32 randInt( const uint32 n ); // integer in [0,n] for n < 2^32
double rand(); // real number in [0,1]
double rand( const double n ); // real number in [0,n]
double randExc(); // real number in [0,1)
double randExc( const double n ); // real number in [0,n)
double randDblExc(); // real number in (0,1)
double randDblExc( const double n ); // real number in (0,n)
double operator()(); // same as rand()
// Access to 53-bit random numbers (capacity of IEEE double precision)
double rand53(); // real number in [0,1)
// Access to nonuniform random number distributions
double randNorm( const double mean = 0.0, const double stddev = 1.0 );
// Re-seeding functions with same behavior as initializers
void seed( const uint32 oneSeed );
void seed( uint32 *const bigSeed, const uint32 seedLength = N );
void seed();
// Saving and loading generator state
void save( uint32* saveArray ) const; // to array of size SAVE
void load( uint32 *const loadArray ); // from such array
friend std::ostream& operator<<( std::ostream& os, const MTRand& mtrand );
friend std::istream& operator>>( std::istream& is, MTRand& mtrand );
MTRand& operator=( const MTRand& o );
protected:
void initialize( const uint32 oneSeed );
void reload();
uint32 hiBit( const uint32 u ) const { return u & 0x80000000UL; }
uint32 loBit( const uint32 u ) const { return u & 0x00000001UL; }
uint32 loBits( const uint32 u ) const { return u & 0x7fffffffUL; }
uint32 mixBits( const uint32 u, const uint32 v ) const
{ return hiBit(u) | loBits(v); }
uint32 magic( const uint32 u ) const
{ return loBit(u) ? 0x9908b0dfUL : 0x0UL; }
uint32 twist( const uint32 m, const uint32 s0, const uint32 s1 ) const
{ return m ^ (mixBits(s0,s1)>>1) ^ magic(s1); }
static uint32 hash( time_t t, clock_t c );
};
// Functions are defined in order of usage to assist inlining
inline MTRand::uint32 MTRand::hash( time_t t, clock_t c )
{
// Get a uint32 from t and c
// Better than uint32(x) in case x is floating point in [0,1]
// Based on code by Lawrence Kirby (fred@genesis.demon.co.uk)
static uint32 differ = 0; // guarantee time-based seeds will change
uint32 h1 = 0;
unsigned char *p = (unsigned char *) &t;
for( size_t i = 0; i < sizeof(t); ++i )
{
h1 *= UCHAR_MAX + 2U;
h1 += p[i];
}
uint32 h2 = 0;
p = (unsigned char *) &c;
for( size_t j = 0; j < sizeof(c); ++j )
{
h2 *= UCHAR_MAX + 2U;
h2 += p[j];
}
return ( h1 + differ++ ) ^ h2;
}
inline void MTRand::initialize( const uint32 seed )
{
// Initialize generator state with seed
// See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
// In previous versions, most significant bits (MSBs) of the seed affect
// only MSBs of the state array. Modified 9 Jan 2002 by Makoto Matsumoto.
register uint32 *s = state;
register uint32 *r = state;
register int i = 1;
*s++ = seed & 0xffffffffUL;
for( ; i < N; ++i )
{
*s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL;
r++;
}
}
inline void MTRand::reload()
{
// Generate N new values in state
// Made clearer and faster by Matthew Bellew (matthew.bellew@home.com)
static const int MmN = int(M) - int(N); // in case enums are unsigned
register uint32 *p = state;
register int i;
for( i = N - M; i--; ++p )
*p = twist( p[M], p[0], p[1] );
for( i = M; --i; ++p )
*p = twist( p[MmN], p[0], p[1] );
*p = twist( p[MmN], p[0], state[0] );
left = N, pNext = state;
}
inline void MTRand::seed( const uint32 oneSeed )
{
// Seed the generator with a simple uint32
initialize(oneSeed);
reload();
}
inline void MTRand::seed( uint32 *const bigSeed, const uint32 seedLength )
{
// Seed the generator with an array of uint32's
// There are 2^19937-1 possible initial states. This function allows
// all of those to be accessed by providing at least 19937 bits (with a
// default seed length of N = 624 uint32's). Any bits above the lower 32
// in each element are discarded.
// Just call seed() if you want to get array from /dev/urandom
initialize(19650218UL);
register int i = 1;
register uint32 j = 0;
register int k = ( N > seedLength ? N : seedLength );
for( ; k; --k )
{
state[i] =
state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL );
state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
state[i] &= 0xffffffffUL;
++i; ++j;
if( i >= N ) { state[0] = state[N-1]; i = 1; }
if( j >= seedLength ) j = 0;
}
for( k = N - 1; k; --k )
{
state[i] =
state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL );
state[i] -= i;
state[i] &= 0xffffffffUL;
++i;
if( i >= N ) { state[0] = state[N-1]; i = 1; }
}
state[0] = 0x80000000UL; // MSB is 1, assuring non-zero initial array
reload();
}
inline void MTRand::seed()
{
// Seed the generator with an array from /dev/urandom if available
// Otherwise use a hash of time() and clock() values
// First try getting an array from /dev/urandom
FILE* urandom = fopen( "/dev/urandom", "rb" );
if( urandom )
{
uint32 bigSeed[N];
register uint32 *s = bigSeed;
register int i = N;
register bool success = true;
while( success && i-- )
success = fread( s++, sizeof(uint32), 1, urandom );
fclose(urandom);
if( success ) { seed( bigSeed, N ); return; }
}
// Was not successful, so use time() and clock() instead
seed( hash( time(NULL), clock() ) );
}
inline MTRand::MTRand( const uint32 oneSeed )
{ seed(oneSeed); }
inline MTRand::MTRand( uint32 *const bigSeed, const uint32 seedLength )
{ seed(bigSeed,seedLength); }
inline MTRand::MTRand()
{ seed(); }
inline MTRand::MTRand( const MTRand& o )
{
register const uint32 *t = o.state;
register uint32 *s = state;
register int i = N;
for( ; i--; *s++ = *t++ ) {}
left = o.left;
pNext = &state[N-left];
}
inline MTRand::uint32 MTRand::randInt()
{
// Pull a 32-bit integer from the generator state
// Every other access function simply transforms the numbers extracted here
if( left == 0 ) reload();
--left;
register uint32 s1;
s1 = *pNext++;
s1 ^= (s1 >> 11);
s1 ^= (s1 << 7) & 0x9d2c5680UL;
s1 ^= (s1 << 15) & 0xefc60000UL;
return ( s1 ^ (s1 >> 18) );
}
inline MTRand::uint32 MTRand::randInt( const uint32 n )
{
// Find which bits are used in n
// Optimized by Magnus Jonsson (magnus@smartelectronix.com)
uint32 used = n;
used |= used >> 1;
used |= used >> 2;
used |= used >> 4;
used |= used >> 8;
used |= used >> 16;
// Draw numbers until one is found in [0,n]
uint32 i;
do
i = randInt() & used; // toss unused bits to shorten search
while( i > n );
return i;
}
inline double MTRand::rand()
{ return double(randInt()) * (1.0/4294967295.0); }
inline double MTRand::rand( const double n )
{ return rand() * n; }
inline double MTRand::randExc()
{ return double(randInt()) * (1.0/4294967296.0); }
inline double MTRand::randExc( const double n )
{ return randExc() * n; }
inline double MTRand::randDblExc()
{ return ( double(randInt()) + 0.5 ) * (1.0/4294967296.0); }
inline double MTRand::randDblExc( const double n )
{ return randDblExc() * n; }
inline double MTRand::rand53()
{
uint32 a = randInt() >> 5, b = randInt() >> 6;
return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0); // by Isaku Wada
}
inline double MTRand::randNorm( const double mean, const double stddev )
{
// Return a real number from a normal (Gaussian) distribution with given
// mean and standard deviation by polar form of Box-Muller transformation
double x, y, r;
do
{
x = 2.0 * rand() - 1.0;
y = 2.0 * rand() - 1.0;
r = x * x + y * y;
}
while ( r >= 1.0 || r == 0.0 );
double s = sqrt( -2.0 * log(r) / r );
return mean + x * s * stddev;
}
inline double MTRand::operator()()