Update docs, some cleanup
This commit is contained in:
parent
f20a29beff
commit
5c913192e9
4 changed files with 307 additions and 241 deletions
3
eo/NEWS
3
eo/NEWS
|
|
@ -1,7 +1,8 @@
|
||||||
* release 0.9.4
|
* release 0.9.4 (not yet released)
|
||||||
- Update introductory pages of documentation and webpage.
|
- Update introductory pages of documentation and webpage.
|
||||||
- Remove support for pre-standard C++ compiler (i.e. gcc-2.x), which allows to
|
- Remove support for pre-standard C++ compiler (i.e. gcc-2.x), which allows to
|
||||||
clean up the code considerably. Assume availability of sstream and limits.
|
clean up the code considerably. Assume availability of sstream and limits.
|
||||||
|
- Implement CMA-ES.
|
||||||
|
|
||||||
|
|
||||||
* release 0.9.3z.1 (1. Oct. 2005)
|
* release 0.9.3z.1 (1. Oct. 2005)
|
||||||
|
|
|
||||||
|
|
@ -102,7 +102,7 @@
|
||||||
</tr>
|
</tr>
|
||||||
<tr>
|
<tr>
|
||||||
<td class="TITLE" align="right" valign="top" width="100">
|
<td class="TITLE" align="right" valign="top" width="100">
|
||||||
<h2>Tutorial</h2>
|
<h2>Documentation</h2>
|
||||||
</td>
|
</td>
|
||||||
<td bgcolor="#ffcc99">
|
<td bgcolor="#ffcc99">
|
||||||
<p>
|
<p>
|
||||||
|
|
@ -120,9 +120,8 @@
|
||||||
The tutorial is also included in
|
The tutorial is also included in
|
||||||
the <a
|
the <a
|
||||||
href="http://sourceforge.net/project/showfiles.php?group_id=9775">released
|
href="http://sourceforge.net/project/showfiles.php?group_id=9775">released
|
||||||
sources</a>.
|
sources</a>.
|
||||||
</p>
|
</p>
|
||||||
|
|
||||||
<p>
|
<p>
|
||||||
The latest
|
The latest
|
||||||
<a
|
<a
|
||||||
|
|
@ -132,6 +131,12 @@
|
||||||
href="http://eodev.sourceforge.net/eo/tutorial/pdf/paradiseoJet7.pdf">introduction
|
href="http://eodev.sourceforge.net/eo/tutorial/pdf/paradiseoJet7.pdf">introduction
|
||||||
to ParadisEO</a>, the parallel version of EO.
|
to ParadisEO</a>, the parallel version of EO.
|
||||||
</p>
|
</p>
|
||||||
|
<p>
|
||||||
|
The complete code is also well documented and you can look at the
|
||||||
|
generated <a
|
||||||
|
href="http://eodev.sourceforge.net/eo/doc/html/index.html">interface
|
||||||
|
documentation</a>.
|
||||||
|
</p>
|
||||||
</td>
|
</td>
|
||||||
</tr>
|
</tr>
|
||||||
<tr>
|
<tr>
|
||||||
|
|
|
||||||
|
|
@ -1,14 +1,14 @@
|
||||||
#ifdef _MSC_VER
|
#ifdef _MSC_VER
|
||||||
// to avoid long name warnings
|
// to avoid long name warnings
|
||||||
#pragma warning(disable:4786)
|
#pragma warning(disable:4786)
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
#include <ctime>
|
#include <ctime>
|
||||||
#include "eoRNG.h"
|
#include "eoRNG.h"
|
||||||
|
|
||||||
namespace eo
|
namespace eo
|
||||||
{
|
{
|
||||||
/// The Global random number generator.
|
/// The Global random number generator.
|
||||||
eoRng rng(time(0));
|
eoRng rng(time(0));
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -1,24 +1,24 @@
|
||||||
/** Random number generator adapted from (see comments below)
|
/** Random number generator adapted from (see comments below)
|
||||||
|
|
||||||
The random number generator is modified into a class
|
The random number generator is modified into a class
|
||||||
by Maarten Keijzer (mak@dhi.dk). Also added the Box-Muller
|
by Maarten Keijzer (mak@dhi.dk). Also added the Box-Muller
|
||||||
transformation to generate normal deviates.
|
transformation to generate normal deviates.
|
||||||
|
|
||||||
This library is free software; you can redistribute it and/or
|
This library is free software; you can redistribute it and/or
|
||||||
modify it under the terms of the GNU Lesser General Public
|
modify it under the terms of the GNU Lesser General Public
|
||||||
License as published by the Free Software Foundation; either
|
License as published by the Free Software Foundation; either
|
||||||
version 2 of the License, or (at your option) any later version.
|
version 2 of the License, or (at your option) any later version.
|
||||||
|
|
||||||
This library is distributed in the hope that it will be useful,
|
This library is distributed in the hope that it will be useful,
|
||||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||||
Lesser General Public License for more details.
|
Lesser General Public License for more details.
|
||||||
|
|
||||||
You should have received a copy of the GNU Lesser General Public
|
You should have received a copy of the GNU Lesser General Public
|
||||||
License along with this library; if not, write to the Free Software
|
License along with this library; if not, write to the Free Software
|
||||||
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
|
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
|
||||||
|
|
||||||
Contact: todos@geneura.ugr.es, http://geneura.ugr.es
|
Contact: todos@geneura.ugr.es, http://geneura.ugr.es
|
||||||
*/
|
*/
|
||||||
|
|
||||||
#ifndef EO_RANDOM_NUMBER_GENERATOR
|
#ifndef EO_RANDOM_NUMBER_GENERATOR
|
||||||
|
|
@ -51,6 +51,15 @@ number generator MT19937 for generating random numbers. The various
|
||||||
member functions implement useful functions for evolutionary
|
member functions implement useful functions for evolutionary
|
||||||
algorithms. Included are: rand(), random(), flip() and normal().
|
algorithms. Included are: rand(), random(), flip() and normal().
|
||||||
|
|
||||||
|
EO provides a global random number generator <tt>rng</tt> that is seeded by the
|
||||||
|
current UNIX time at program start. Moreover some global convenience functions
|
||||||
|
are provided that use the global random number generator: <tt>random</tt>,
|
||||||
|
<tt>normal</tt>.
|
||||||
|
|
||||||
|
@warning If you want to repeatedly generated the same sequence of pseudo-random
|
||||||
|
numbers, you should always reseed the generator at the beginning of your code.
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
<h1>Documentation in original file</h1>
|
<h1>Documentation in original file</h1>
|
||||||
|
|
||||||
|
|
@ -77,19 +86,19 @@ depending on data type sizes, and the code is quite short as well). It generates
|
||||||
random numbers in batches of 624 at a time, so the caching and pipelining of
|
random numbers in batches of 624 at a time, so the caching and pipelining of
|
||||||
modern systems is exploited. It is also divide- and mod-free.
|
modern systems is exploited. It is also divide- and mod-free.
|
||||||
|
|
||||||
The code as Shawn received it included the following notice:
|
The code as Shawn received it included the following notice: <tt>Copyright (C)
|
||||||
- Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura. When you use this,
|
1997 Makoto Matsumoto and Takuji Nishimura. When you use this, send an e-mail to
|
||||||
send an e-mail to <matumoto@math.keio.ac.jp> with an appropriate reference to
|
<matumoto@math.keio.ac.jp> with an appropriate reference to your work.</tt> It
|
||||||
your work.
|
would be nice to Cc: <Cokus@math.washington.edu> and
|
||||||
- It would be nice to CC: <Cokus@math.washington.edu> when you write.
|
<eodev-main@lists.sourceforge.net> when you write.
|
||||||
|
|
||||||
|
|
||||||
<h1>Portability</h1>
|
<h1>Portability</h1>
|
||||||
|
|
||||||
Note for people porting EO to other platforms: please make sure that the type
|
Note for people porting EO to other platforms: please make sure that the type
|
||||||
uint32_t in the file eoRng.h is exactly 32 bits long. It may be longer, but not
|
uint32_t in the file eoRng.h is exactly 32 bits long. It may in principle be
|
||||||
shorter. If it is longer, file compatibility between EO on different platforms
|
longer, but not shorter. If it is longer, file compatibility between EO on
|
||||||
may be broken.
|
different platforms may be broken.
|
||||||
*/
|
*/
|
||||||
class eoRng : public eoObject, public eoPersistent
|
class eoRng : public eoObject, public eoPersistent
|
||||||
{
|
{
|
||||||
|
|
@ -130,85 +139,87 @@ public :
|
||||||
|
|
||||||
/** Re-initializes the Random Number Generator
|
/** Re-initializes the Random Number Generator
|
||||||
|
|
||||||
This is the traditional seeding procedure.
|
This is the traditional seeding procedure. This version is deprecated and
|
||||||
|
only provided for compatibility with old code. In new projects you should
|
||||||
|
use reseed.
|
||||||
|
|
||||||
@see reseed for details on usage of the seeding value.
|
@see reseed for details on usage of the seeding value.
|
||||||
|
|
||||||
@version old version
|
@version old version (deprecated)
|
||||||
*/
|
*/
|
||||||
void oldReseed(uint32_t s)
|
void oldReseed(uint32_t s)
|
||||||
{
|
{
|
||||||
initialize(s);
|
initialize(s);
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
uniform(m = 1.0) returns a random double in the range [0, m)
|
uniform(m = 1.0) returns a random double in the range [0, m)
|
||||||
*/
|
*/
|
||||||
double uniform(double m = 1.0)
|
double uniform(double m = 1.0)
|
||||||
{ // random number between [0, m]
|
{ // random number between [0, m]
|
||||||
return m * double(rand()) / double(1.0 + rand_max());
|
return m * double(rand()) / double(1.0 + rand_max());
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
|
||||||
random() returns a random integer in the range [0, m)
|
|
||||||
*/
|
|
||||||
uint32_t random(uint32_t m)
|
|
||||||
{
|
|
||||||
return uint32_t(uniform() * double(m));
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
flip() tosses a biased coin such that flip(x/100.0) will
|
|
||||||
returns true x% of the time
|
|
||||||
*/
|
|
||||||
bool flip(float bias=0.5)
|
|
||||||
{
|
|
||||||
return uniform() < bias;
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
normal() zero mean gaussian deviate with standard deviation of 1
|
|
||||||
*/
|
|
||||||
double normal(void); // gaussian mutation, stdev 1
|
|
||||||
|
|
||||||
/**
|
|
||||||
normal(stdev) zero mean gaussian deviate with user defined standard deviation
|
|
||||||
*/
|
|
||||||
double normal(double stdev)
|
|
||||||
{
|
|
||||||
return stdev * normal();
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
normal(mean, stdev) user defined mean gaussian deviate with user defined standard deviation
|
|
||||||
*/
|
|
||||||
double normal(double mean, double stdev)
|
|
||||||
{
|
|
||||||
return mean + normal(stdev);
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
Generates random numbers using a negative exponential distribution
|
|
||||||
*/
|
|
||||||
double negexp(double mean)
|
|
||||||
{
|
|
||||||
return ( -mean*log((double)rand() / rand_max()));
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
rand() returns a random number in the range [0, rand_max)
|
|
||||||
*/
|
|
||||||
uint32_t rand();
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
rand_max() the maximum returned by rand()
|
random() returns a random integer in the range [0, m)
|
||||||
|
*/
|
||||||
|
uint32_t random(uint32_t m)
|
||||||
|
{
|
||||||
|
return uint32_t(uniform() * double(m));
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
flip() tosses a biased coin such that flip(x/100.0) will
|
||||||
|
returns true x% of the time
|
||||||
|
*/
|
||||||
|
bool flip(float bias=0.5)
|
||||||
|
{
|
||||||
|
return uniform() < bias;
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
normal() zero mean gaussian deviate with standard deviation of 1
|
||||||
|
*/
|
||||||
|
double normal(void); // gaussian mutation, stdev 1
|
||||||
|
|
||||||
|
/**
|
||||||
|
normal(stdev) zero mean gaussian deviate with user defined standard deviation
|
||||||
|
*/
|
||||||
|
double normal(double stdev)
|
||||||
|
{
|
||||||
|
return stdev * normal();
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
normal(mean, stdev) user defined mean gaussian deviate with user defined standard deviation
|
||||||
|
*/
|
||||||
|
double normal(double mean, double stdev)
|
||||||
|
{
|
||||||
|
return mean + normal(stdev);
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
Generates random numbers using a negative exponential distribution
|
||||||
|
*/
|
||||||
|
double negexp(double mean)
|
||||||
|
{
|
||||||
|
return ( -mean*log((double)rand() / rand_max()));
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
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(void) const { return uint32_t(0xffffffff); }
|
uint32_t rand_max(void) const { return uint32_t(0xffffffff); }
|
||||||
|
|
||||||
/**
|
/**
|
||||||
roulette_wheel(vec, total = 0) does a roulette wheel selection
|
roulette_wheel(vec, total = 0) does a roulette wheel selection
|
||||||
on the input std::vector vec. If the total is not supplied, it is
|
on the input std::vector vec. If the total is not supplied, it is
|
||||||
calculated. It returns an integer denoting the selected argument.
|
calculated. It returns an integer denoting the selected argument.
|
||||||
*/
|
*/
|
||||||
template <typename TYPE>
|
template <typename TYPE>
|
||||||
int roulette_wheel(const std::vector<TYPE>& vec, TYPE total = 0)
|
int roulette_wheel(const std::vector<TYPE>& vec, TYPE total = 0)
|
||||||
|
|
@ -241,79 +252,149 @@ public :
|
||||||
|
|
||||||
@overload
|
@overload
|
||||||
|
|
||||||
|
Provide a version returning a non-const element reference.
|
||||||
|
|
||||||
@return Uniformly chosen element from the vector.
|
@return Uniformly chosen element from the vector.
|
||||||
|
|
||||||
|
@warning Changing the return value does alter the vector.
|
||||||
*/
|
*/
|
||||||
template <typename TYPE>
|
template <typename TYPE>
|
||||||
TYPE& choice(std::vector<TYPE>& vec)
|
TYPE& choice(std::vector<TYPE>& vec)
|
||||||
{ return vec[random(vec.size())]; }
|
{ return vec[random(vec.size())]; }
|
||||||
|
|
||||||
///
|
///
|
||||||
void printOn(std::ostream& _os) const
|
void printOn(std::ostream& _os) const
|
||||||
{
|
{
|
||||||
for (int i = 0; i < N; ++i)
|
for (int i = 0; i < N; ++i)
|
||||||
{
|
{
|
||||||
_os << state[i] << ' ';
|
_os << state[i] << ' ';
|
||||||
}
|
}
|
||||||
_os << int(next - state) << ' ';
|
_os << int(next - state) << ' ';
|
||||||
_os << left << ' ' << cached << ' ' << cacheValue;
|
_os << left << ' ' << cached << ' ' << cacheValue;
|
||||||
}
|
}
|
||||||
|
|
||||||
///
|
///
|
||||||
void readFrom(std::istream& _is)
|
void readFrom(std::istream& _is)
|
||||||
{
|
{
|
||||||
for (int i = 0; i < N; ++i)
|
for (int i = 0; i < N; ++i)
|
||||||
{
|
{
|
||||||
_is >> state[i];
|
_is >> state[i];
|
||||||
}
|
}
|
||||||
|
|
||||||
int n;
|
int n;
|
||||||
_is >> n;
|
_is >> n;
|
||||||
next = state + n;
|
next = state + n;
|
||||||
|
|
||||||
_is >> left;
|
_is >> left;
|
||||||
_is >> cached;
|
_is >> cached;
|
||||||
_is >> cacheValue;
|
_is >> cacheValue;
|
||||||
}
|
}
|
||||||
|
|
||||||
std::string className(void) const { return "Mersenne-Twister"; }
|
std::string className(void) const { return "Mersenne-Twister"; }
|
||||||
|
|
||||||
private :
|
private:
|
||||||
uint32_t restart(void);
|
|
||||||
void initialize(uint32_t seed);
|
|
||||||
|
|
||||||
uint32_t* state; // the array for the state
|
uint32_t restart(void);
|
||||||
uint32_t* next;
|
|
||||||
int left;
|
|
||||||
|
|
||||||
// for normal distribution
|
|
||||||
bool cached;
|
|
||||||
float cacheValue;
|
|
||||||
|
|
||||||
const int N;
|
|
||||||
const int M;
|
|
||||||
const uint32_t K; // a magic constant
|
|
||||||
|
|
||||||
|
|
||||||
/**
|
/* @brief Initialize state
|
||||||
Private copy ctor and assignment operator to make sure that
|
|
||||||
nobody accidentally copies the random number generator.
|
We initialize state[0..(N-1)] via the generator
|
||||||
If you want similar RNG's, make two RNG's and initialize
|
|
||||||
them with the same seed.
|
<tt>x_new = (69069 * x_old) mod 2^32</tt>
|
||||||
*/
|
|
||||||
eoRng (const eoRng&); // no implementation
|
from Line 15 of Table 1, p. 106, Sec. 3.3.4 of Knuth's _The Art of Computer
|
||||||
eoRng& operator=(const eoRng&); // dito
|
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;
|
||||||
|
|
||||||
|
uint32_t* next;
|
||||||
|
|
||||||
|
int left;
|
||||||
|
|
||||||
|
// for normal distribution
|
||||||
|
bool cached;
|
||||||
|
|
||||||
|
float cacheValue;
|
||||||
|
|
||||||
|
const int N;
|
||||||
|
|
||||||
|
const int M;
|
||||||
|
|
||||||
|
/** @brief Magic constant */
|
||||||
|
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 Assignmant operator
|
||||||
|
|
||||||
|
@see Copy constructor eoRng(const eoRng&).
|
||||||
|
*/
|
||||||
|
eoRng& operator=(const eoRng&);
|
||||||
};
|
};
|
||||||
|
|
||||||
/**
|
|
||||||
The one and only global eoRng object
|
|
||||||
*/
|
|
||||||
namespace eo
|
namespace eo
|
||||||
{
|
{
|
||||||
extern eoRng rng;
|
/** The one and only global eoRng object */
|
||||||
|
extern eoRng rng;
|
||||||
}
|
}
|
||||||
|
|
||||||
using eo::rng;
|
using eo::rng;
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
// Implementation of some eoRng members.... Don't mind the mess, it does work.
|
// Implementation of some eoRng members.... Don't mind the mess, it does work.
|
||||||
|
|
||||||
|
|
||||||
|
|
@ -323,139 +404,118 @@ using eo::rng;
|
||||||
#define mixBits(u, v) (hiBit(u)|loBits(v)) // move hi bit of u to hi bit of v
|
#define mixBits(u, v) (hiBit(u)|loBits(v)) // move hi bit of u to hi bit of v
|
||||||
|
|
||||||
inline void eoRng::initialize(uint32_t seed)
|
inline void eoRng::initialize(uint32_t seed)
|
||||||
{
|
{
|
||||||
//
|
left = -1;
|
||||||
// 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.
|
|
||||||
//
|
|
||||||
|
|
||||||
left = -1;
|
register uint32_t x = (seed | 1U) & 0xFFFFFFFFU, *s = state;
|
||||||
|
register int j;
|
||||||
|
|
||||||
register uint32_t x = (seed | 1U) & 0xFFFFFFFFU, *s = state;
|
for(left=0, *s++=x, j=N; --j;
|
||||||
register int j;
|
*s++ = (x*=69069U) & 0xFFFFFFFFU);
|
||||||
|
}
|
||||||
|
|
||||||
for(left=0, *s++=x, j=N; --j;
|
|
||||||
*s++ = (x*=69069U) & 0xFFFFFFFFU);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
inline uint32_t eoRng::restart(void)
|
inline uint32_t eoRng::restart(void)
|
||||||
{
|
{
|
||||||
register uint32_t *p0=state, *p2=state+2, *pM=state+M, s0, s1;
|
register uint32_t *p0=state, *p2=state+2, *pM=state+M, s0, s1;
|
||||||
register int j;
|
register int j;
|
||||||
|
|
||||||
left=N-1, next=state+1;
|
left=N-1, next=state+1;
|
||||||
|
|
||||||
for(s0=state[0], s1=state[1], j=N-M+1; --j; s0=s1, s1=*p2++)
|
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);
|
*p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);
|
||||||
|
|
||||||
for(pM=state, j=M; --j; s0=s1, s1=*p2++)
|
for(pM=state, j=M; --j; s0=s1, s1=*p2++)
|
||||||
*p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);
|
*p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);
|
||||||
|
|
||||||
s1=state[0], *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 >> 11);
|
||||||
s1 ^= (s1 << 7) & 0x9D2C5680U;
|
s1 ^= (s1 << 7) & 0x9D2C5680U;
|
||||||
s1 ^= (s1 << 15) & 0xEFC60000U;
|
s1 ^= (s1 << 15) & 0xEFC60000U;
|
||||||
return(s1 ^ (s1 >> 18));
|
return(s1 ^ (s1 >> 18));
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
inline uint32_t eoRng::rand(void)
|
inline uint32_t eoRng::rand(void)
|
||||||
{
|
{
|
||||||
|
|
||||||
uint32_t y;
|
uint32_t y;
|
||||||
|
|
||||||
|
if(--left < 0)
|
||||||
|
return(restart());
|
||||||
|
|
||||||
|
y = *next++;
|
||||||
|
y ^= (y >> 11);
|
||||||
|
y ^= (y << 7) & 0x9D2C5680U;
|
||||||
|
y ^= (y << 15) & 0xEFC60000U;
|
||||||
|
return(y ^ (y >> 18));
|
||||||
|
}
|
||||||
|
|
||||||
if(--left < 0)
|
|
||||||
return(restart());
|
|
||||||
|
|
||||||
y = *next++;
|
|
||||||
y ^= (y >> 11);
|
|
||||||
y ^= (y << 7) & 0x9D2C5680U;
|
|
||||||
y ^= (y << 15) & 0xEFC60000U;
|
|
||||||
return(y ^ (y >> 18));
|
|
||||||
}
|
|
||||||
|
|
||||||
inline double eoRng::normal(void)
|
inline double eoRng::normal(void)
|
||||||
{
|
{
|
||||||
if (cached)
|
if (cached)
|
||||||
{
|
{
|
||||||
cached = false;
|
cached = false;
|
||||||
return cacheValue;
|
return cacheValue;
|
||||||
}
|
}
|
||||||
|
|
||||||
float rSquare, factor, var1, var2;
|
float rSquare, factor, var1, var2;
|
||||||
|
|
||||||
do
|
do
|
||||||
{
|
{
|
||||||
var1 = 2.0 * uniform() - 1.0;
|
var1 = 2.0 * uniform() - 1.0;
|
||||||
var2 = 2.0 * uniform() - 1.0;
|
var2 = 2.0 * uniform() - 1.0;
|
||||||
|
|
||||||
rSquare = var1 * var1 + var2 * var2;
|
rSquare = var1 * var1 + var2 * var2;
|
||||||
}
|
}
|
||||||
while (rSquare >= 1.0 || rSquare == 0.0);
|
while (rSquare >= 1.0 || rSquare == 0.0);
|
||||||
|
|
||||||
factor = sqrt(-2.0 * log(rSquare) / rSquare);
|
factor = sqrt(-2.0 * log(rSquare) / rSquare);
|
||||||
|
|
||||||
cacheValue = var1 * factor;
|
cacheValue = var1 * factor;
|
||||||
cached = true;
|
cached = true;
|
||||||
|
|
||||||
return (var2 * factor);
|
return (var2 * factor);
|
||||||
}
|
}
|
||||||
|
|
||||||
namespace eo {
|
|
||||||
// a few convenience functions for generating numbers
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Templatized random function, works with most basic types such as:
|
namespace eo {
|
||||||
* char
|
// a few convenience functions for generating numbers
|
||||||
* int
|
|
||||||
* unsigned
|
/** @brief Random function
|
||||||
* float
|
|
||||||
* double
|
Templatized random function, returns a random double in the range [0, max).
|
||||||
*/
|
|
||||||
|
@param max Maximum for distribution
|
||||||
|
|
||||||
|
It works with most basic types such as:
|
||||||
|
- char
|
||||||
|
- int
|
||||||
|
- unsigned
|
||||||
|
- float
|
||||||
|
- double
|
||||||
|
*/
|
||||||
template <typename T>
|
template <typename T>
|
||||||
inline
|
inline T random(const T& max) {
|
||||||
T random(const T& mx) { return static_cast<T>(rng.uniform() * mx); }
|
return static_cast<T>(rng.uniform() * max); }
|
||||||
|
|
||||||
|
/** @brief Random function
|
||||||
|
|
||||||
|
Templatized random function, returns a random double in the range [min, max).
|
||||||
|
|
||||||
|
@param min Minimum for distribution
|
||||||
|
@param max Maximum for distribution
|
||||||
|
|
||||||
|
@see random(const T& max)
|
||||||
|
*/
|
||||||
|
template <typename T>
|
||||||
|
inline T random(const T& min, const T& max) {
|
||||||
|
return static_cast<T>(rng.uniform() * (max-min)) + min; }
|
||||||
|
|
||||||
/** Normal distribution */
|
/** Normal distribution */
|
||||||
inline double normal() { return rng.normal(); }
|
inline double normal() { return rng.normal(); }
|
||||||
|
|
|
||||||
Reference in a new issue