/** Random number generator adapted from (see comments below)

    The random number generator is modified into a class
    by Maarten Keijzer (mak@dhi.dk). Also added the Box-Muller
    transformation to generate normal deviates.

    This library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
    License as published by the Free Software Foundation; either
    version 2 of the License, or (at your option) any later version.

    This library is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with this library; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

    Contact: todos@geneura.ugr.es, http://geneura.ugr.es
*/

#ifndef EO_RANDOM_NUMBER_GENERATOR
#define EO_RANDOM_NUMBER_GENERATOR

// uint32_t is an unsigned integer type capable of holding 32 bits.
//
// In the applicatione here exactly 32 but should typically be fastest, but 64
// might be better on an Alpha with GCC at -O3 optimization so try your options
// and see what's best for you.
//
// The C99-standard defines uint32_t to be declared in stdint.h, but some
// systmes don't have that and implement it in inttypes.h.
#if (! defined __sun)
#include <stdint.h>
#else
#include <inttypes.h>
#endif

#include <vector>
#include "eoPersistent.h"
#include "eoObject.h"

/** Random Number Generator

@class eoRng eoRNG.h utils/eoRNG.h

eoRng is a persistent class that uses the ``Mersenne Twister'' random
number generator MT19937 for generating random numbers. The various
member functions implement useful functions for evolutionary
algorithms. Included are: rand(), random(), flip() and normal().


<h1>Documentation in original file</h1>

This is the ``Mersenne Twister'' random number generator MT19937, which
generates pseudorandom integers uniformly distributed in 0..(2^32 - 1) starting
from any odd seed in 0..(2^32 - 1). This version is a recode by Shawn Cokus
(Cokus@math.washington.edu) on March 8, 1998 of a version by Takuji Nishimura
(who had suggestions from Topher Cooper and Marc Rieffel in July-August 1997).

Effectiveness of the recoding (on Goedel2.math.washington.edu, a DEC Alpha
running OSF/1) using GCC -O3 as a compiler: before recoding: 51.6 sec. to
generate 300 million random numbers; after recoding: 24.0 sec. for the same
(i.e., 46.5% of original time), so speed is now about 12.5 million random number
generations per second on this machine.

According to the URL <http://www.math.keio.ac.jp/~matumoto/emt.html> (and
paraphrasing a bit in places), the Mersenne Twister is ``designed with
consideration of the flaws of various existing generators,'' has a period of
2^19937 - 1, gives a sequence that is 623-dimensionally equidistributed, and
``has passed many std::stringent tests, including the die-hard test of G.
Marsaglia and the load test of P. Hellekalek and S. Wegenkittl.'' It is
efficient in memory usage (typically using 2506 to 5012 bytes of static data,
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
modern systems is exploited. It is also divide- and mod-free.

The code as Shawn received it included the following notice:
- Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura. When you use this,
  send an e-mail to <matumoto@math.keio.ac.jp> with an appropriate reference to
  your work.
- It would be nice to CC: <Cokus@math.washington.edu> when you write.


<h1>Portability</h1>

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
shorter. If it is longer, file compatibility between EO on different platforms
may be broken.
*/
class eoRng : public eoObject, public eoPersistent
{
public :

    /** Constructor

    @param s Random seed; if you want another seed, use reseed.

    @see reseed for details on usage of the seeding value.
    */
    eoRng(uint32_t s)
        : state(0), next(0), left(-1), cached(false), N(624), M(397), K(0x9908B0DFU)
        {
            state = new uint32_t[N+1];
            initialize(2*s);
        }

    ~eoRng(void)
        {
            delete [] state;
        }

    /** Re-initializes the Random Number Generator.

    WARNING: Jeroen Eggermont <jeggermo@liacs.nl> noticed that initialize does
    not differentiate between odd and even numbers, therefore the argument to
    reseed is now doubled before being passed on.

    Manually divide the seed by 2 if you want to re-run old runs

    @version MS. 5 Oct. 2001
    */
    void reseed(uint32_t s)
        {
            initialize(2*s);
        }

    /** Re-initializes the Random Number Generator

    This is the traditional seeding procedure.

    @see reseed for details on usage of the seeding value.

    @version old version
    */
    void oldReseed(uint32_t s)
        {
            initialize(s);
        }

 /**
    uniform(m = 1.0) returns a random double in the range [0, m)
 */
 double uniform(double m = 1.0)
   { // random number between [0, m]
     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()
    */
    uint32_t rand_max(void) const { return uint32_t(0xffffffff); }

    /**
       roulette_wheel(vec, total = 0) does a roulette wheel selection
       on the input std::vector vec. If the total is not supplied, it is
       calculated. It returns an integer denoting the selected argument.
    */
    template <typename TYPE>
    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];
            }
            double fortune = uniform() * total;
            int i = 0;
            while (fortune > 0)
            {
                fortune -= vec[i++];
            }
            return --i;
        }


    /** Randomly select element from vector.

    @return Uniformly chosen element from the vector.
    */
    template <typename TYPE>
    const TYPE& choice(const std::vector<TYPE>& vec) const
        { return vec[random(vec.size())]; }


    /** Randomly select element from vector.

    @overload

    @return Uniformly chosen element from the vector.
    */
    template <typename TYPE>
    TYPE& choice(std::vector<TYPE>& vec)
        { return vec[random(vec.size())]; }

 ///
 void printOn(std::ostream& _os) const
   {
     for (int i = 0; i < N; ++i)
       {
	 _os << state[i] << ' ';
       }
     _os << int(next - state) << ' ';
     _os << left << ' ' << cached << ' ' << cacheValue;
   }

 ///
 void readFrom(std::istream& _is)
   {
     for (int i = 0; i < N; ++i)
       {
	 _is >> state[i];
       }

     int n;
     _is >> n;
     next = state + n;

     _is >> left;
     _is >> cached;
     _is >> cacheValue;
   }

 std::string className(void) const { return "Mersenne-Twister"; }

private :
  uint32_t restart(void);
 void initialize(uint32_t seed);

 uint32_t* state; // the array for the state
 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


 /**
	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.
 */
  eoRng (const eoRng&);				// no implementation
  eoRng& operator=(const eoRng&);	// dito
};

/**
	The one and only global eoRng object
*/
namespace eo
{
extern eoRng rng;
}

using eo::rng;

// 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

inline void eoRng::initialize(uint32_t seed)
 {
    //
    // 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;

   for(left=0, *s++=x, j=N; --j;
       *s++ = (x*=69069U) & 0xFFFFFFFFU);
 }


inline uint32_t eoRng::restart(void)
{
  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(void)
 {

   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));
 }

inline double eoRng::normal(void)
{
  if (cached)
    {
      cached = false;
      return cacheValue;
    }

  float rSquare, factor, 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);

  factor = sqrt(-2.0 * log(rSquare) / rSquare);

  cacheValue = var1 * factor;
  cached = true;

  return (var2 * factor);
}

namespace eo {
// a few convenience functions for generating numbers

    /**
     * Templatized random function, works with most basic types such as:
     *	char
     *	int
     *	unsigned
     *	float
     *	double
     */
    template <typename T>
	inline
	T random(const T& mx) { return static_cast<T>(rng.uniform() * mx); }

    /** Normal distribution */
    inline double normal() { return rng.normal(); }
}


#endif


// Local Variables:
// coding: iso-8859-1
// mode: C++
// c-file-style: "Stroustrup"
// fill-column: 80
// End: