// Copyright (C) 2002 Yoshua Bengio (bengioy@iro.umontreal.ca)
//                and Samy Bengio (bengio@idiap.ch)
//                and Ronan Collobert (collober@iro.umontreal.ca)
//                
//
// This file is part of Torch. Release II.
// [The Ultimate Machine Learning Library]
//
// Torch is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// Torch 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 General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Torch; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

#include "random.h"

namespace Torch {

void getShuffledIndices(int *indices, int n_indices)
{
  for(int i = 0; i < n_indices; i++)
    indices[i] = i;

  shuffle(indices, sizeof(int), n_indices);
}

void shuffle(void *tabular, int size_elem, int n_elems)
{
  void *save = xalloc(size_elem);
  char *tab = (char *)tabular;

  for(int i = 0; i < n_elems; i++)
  {
    int z = (int)(uniform() * ((real)(n_elems-i-1)) );
    memcpy(save, tab+i*size_elem, size_elem);
    memcpy(tab+i*size_elem, tab+(z+i)*size_elem, size_elem);
    memcpy(tab+(z+i)*size_elem, save, size_elem);
  }
  free(save);
}

/*  
  The static data to store the seed used by the random number generators.
  */

static long  the_seed;
static int   iset=0;
static real gset;

/*  
  Special functions.
  =================
  */

/*  
  gammln(): returns the natural logarithm of the gamma function from the 'numerical recipes'.
  */

real  gammln(real xx)
{
  double x,y,tmp,ser;
  static double cof[6]={ 76.18009172947146     ,
                        -86.50532032941677     ,
	                 24.01409824083091     ,
                         -1.231739572450155    ,
		          0.1208650973866179e-2,
                         -0.5395239384953e-5   };
  int j;

  y=x=xx;
  tmp=x+5.5;
  tmp -= (x+0.5)*log(tmp);
  ser=1.000000000190015;
  for (j=0;j<=5;j++) ser += cof[j]/++y;
  return -tmp+log(2.5066282746310005*ser/x);
}

/*   
  Utilities for random numbers generation. 
  =======================================
  */

/*  
  manual_seed(): gives a seed for random number generators.

  Rem: - The stored value is negative.
  */

void  manual_seed(long x)
{
  the_seed = - labs(x);
  iset     = 0;
}

/*  
  seed(): generates a seed for random number generators, using the cpu time.
  */

void  seed()
{
  time_t  ltime;
  struct  tm *today;
  time(&ltime);
  today = localtime(&ltime);
  manual_seed((long)today->tm_sec);
}

/*  
  get_seed(): returns the current value of the 'seed'.
  */

long  get_seed()
{
  long seed = the_seed;
  return seed;
}

/*  
  Constants used by the 'numerical recipes' random number generators.
  */

#define IA 16807                /*   needs for ran1          */
#define IM 2147483647           /*   needs for ran1          */
#define AM (1.0/IM)             /*   needs for ran1          */
#define IQ 127773               /*   needs for ran1          */
#define IR 2836                 /*   needs for ran1          */
#define NTAB 32                 /*   needs for ran1 & ran2   */
#define NDIV (1+(IM-1)/NTAB)    /*   needs for ran1          */
#define EPS 1.2e-7              /*   needs for ran1 & ran2   */
#define RNMX (1.0-EPS)          /*   needs for ran1 & ran2   */
#define IM1 2147483563          /*   needs for ran2          */
#define IM2 2147483399          /*   needs for ran2          */
#define AM1 (1.0/IM1)           /*   needs for ran2          */
#define IMM1 (IM1-1)            /*   needs for ran2          */
#define IA1 40014               /*   needs for ran2          */
#define IA2 40692               /*   needs for ran2          */
#define IQ1 53668               /*   needs for ran2          */
#define IQ2 52774               /*   needs for ran2          */
#define IR1 12211               /*   needs for ran2          */
#define IR2 3791                /*   needs for ran2          */
#define NDIV1 (1+IMM1/NTAB)     /*   needs for ran2          */

/*  
  ran1(): minimal ramdom number generator from the 'numerical recipes'.

  Rem: - It is the 'Minimal' random number generator of Park and Miller
         with Bays-Durham shuffle and added safeguards.

       - Returns a uniform random deviate between 0.0 and 1.0
         (exclusive of the endpoint values).

       - Initilized with a negative seed.
  */

real  ran1()
{
  int j;
  long k;
  static long iy=0;
  static long iv[NTAB];
  real temp;

  if (the_seed <= 0 || !iy) {
    if (-the_seed < 1) the_seed=1;
    else the_seed = -the_seed;
    for (j=NTAB+7;j>=0;j--) {
      k=the_seed/IQ;
      the_seed=IA*(the_seed-k*IQ)-IR*k;
      if (the_seed < 0) the_seed += IM;
      if (j < NTAB) iv[j] = the_seed;
    }
    iy=iv[0];
  }
  k=the_seed/IQ;
  the_seed=IA*(the_seed-k*IQ)-IR*k;
  if (the_seed < 0) the_seed += IM;
  j=iy/NDIV;
  iy=iv[j];
  iv[j] = the_seed;
  if ((temp=AM*iy) > RNMX) return RNMX;
  else return temp;
}

/*  
  ran2(): long period ramdom number generator from the 'numerical recipes'.

  Rem: - It is a long period (> 2 x 10^18) random number generator of L'Ecuyer
         with Bays-Durham shuffle and added safeguards.

       - Returns a uniform random deviate between 0.0 and 1.0
         (exclusive of the endpoint values).

       - Initilized with a negative seed.
  */

real uniform()  /*   real ran2()   */
{
  int j;
  long k;
  static long idum2=123456789;
  static long iy=0;
  static long iv[NTAB];
  real temp;

  if (the_seed <= 0) {
    if (-the_seed < 1) the_seed=1;
    else the_seed = -the_seed;
    idum2=the_seed;
    for (j=NTAB+7;j>=0;j--) {
      k=the_seed/IQ1;
      the_seed=IA1*(the_seed-k*IQ1)-k*IR1;
      if (the_seed < 0) the_seed += IM1;
      if (j < NTAB) iv[j] = the_seed;
    }
    iy=iv[0];
  }
  k=the_seed/IQ1;
  the_seed=IA1*(the_seed-k*IQ1)-k*IR1;
  if (the_seed < 0) the_seed += IM1;
  k=idum2/IQ2;
  idum2=IA2*(idum2-k*IQ2)-k*IR2;
  if (idum2 < 0) idum2 += IM2;
  j=iy/NDIV1;
  iy=iv[j]-idum2;
  iv[j] = the_seed;
  if (iy < 1) iy += IMM1;
  if ((temp=AM1*iy) > RNMX) return RNMX;
  else return temp;
}

/*  
  bounded_uniform(): return an uniform random generator in [a,b].
  */

real  bounded_uniform(real a,real b)
{
  real res = uniform()*(b-a) + a;
  if (res >= b) return b*RNMX;
  else return res;
}

#undef IA
#undef IM
#undef AM
#undef IQ
#undef IR
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX
#undef IM1
#undef IM2
#undef AM1
#undef IMM1
#undef IA1
#undef IA2
#undef IQ1
#undef IQ2
#undef IR1
#undef IR2
#undef NDIV1

/*  
  TRANSFORMATION METHOD:
  ---------------------
  */

/*  
  expdev(): exponential deviate random number from the 'numerical recipes'.
  */

real  expdev()
{
  real dum;

  do
    dum=uniform();
  while (dum == 0.0);
  return -log(dum);
}

/*  
  gasdev(): gaussian (normal) distributed deviate from the 'numerical recioes'.

  Rem: - The current gaussian deviate has zero as mean and unit as variance.
         It uses ran1() as source of uniform deviates.

       - i.e. N(0,1).
  */

real gaussian_01()  /*   real gasdev()   */
{
  real fac,rsq,v1,v2;

  if(the_seed < 0) iset=0;
  if (iset == 0) {
    do {
      v1=2.0*uniform()-1.0;
      v2=2.0*uniform()-1.0;
      rsq=v1*v1+v2*v2;
    } while (rsq >= 1.0 || rsq == 0.0);
    fac=sqrt(-2.0*log(rsq)/rsq);
    gset=v1*fac;
    iset=1;
    return v2*fac;
  } else {
    iset=0;
    return gset;
  }
}

/*  
  gaussian_mu_sigma(): returns a gaussian distributed random number
                       with mean 'mu' and standard deviation 'sigma'.

  Rem: - i.e. N(mu,sigma).
  */

real  gaussian_mu_sigma(real mu, real sigma)
{
  return gaussian_01() * sigma + mu;
}

/*  
  gaussian_misture_mu_sigma(): returns a random number with mixture of gaussians,
                               'w' is the mixture (must be positive summing to 1).

  Rem: - i.e. SUM w[i] * N(mu[i],sigma[i])
  */


/*  
  REJECTION METHOD:
  ----------------
  */

/*  
  gamdev(): returns a deviate distributed as a gamma distribution from the 'numerical recipes'.
  */

real  gamdev(int ia)
{
  int j;
  real am,e,s,v1,v2,x,y;

  if (ia < 1) error("random: error in routine gamdev");
  if (ia < 6) {
    x=1.0;
    for (j=1;j<=ia;j++) x *= uniform();
    x = -log(x);
  } else {
    do {
      do {
        do {
	  v1=uniform();
	  v2=2.0*uniform()-1.0;
	} while (v1*v1+v2*v2 > 1.0);
	y=v2/v1;
	am=ia-1;
	s=sqrt(2.0*am+1.0);
	x=s*y+am;
      } while (x <= 0.0);
      e=(1.0+y*y)*exp(am*log(x/am)-s*y);
    } while (uniform() > e);
  }
  return x;
}

/*  
  poidev(): returns a deviate distributed as a poisson distribution of mean (lambda) 'xm'
            from the 'numerical recipes'.
  */

real  poidev(real xm)
{
  static real sq,alxm,g,oldm=(-1.0);
  real em,t,y;

  if (xm < 12.0) {
    if (xm != oldm) {
      oldm=xm;
      g=exp(-xm);
    }
    em = -1;
    t=1.0;
    do {
      ++em;
      t *= uniform();
    } while (t > g);
  } else {
    if (xm != oldm) {
      oldm=xm;
      sq=sqrt(2.0*xm);
      alxm=log(xm);
      g=xm*alxm-gammln(xm+1.0);
    }
    do {
      do {
        y=tan(M_PI*uniform());
	em=sq*y+xm;
      } while (em < 0.0);
      em=floor(em);
      t=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);
    } while (uniform() > t);
  }
  return em;
}

/*  
  bnldev(): return a random deviate drawn from a binomial distribution of 'n' trials
            each of probability 'pp', from 'numerical recipes'.

  Rem: - The returned type is an real although a binomial random variable is an integer.
  */

real  bnldev(real pp, int n)
{
  int j;
  static int nold=(-1);
  real am,em,g,angle,p,bnl,sq,t,y;
  static real pold=(-1.0),pc,plog,pclog,en,oldg;

  p=(pp <= 0.5 ? pp : 1.0-pp);
  am=n*p;
  if (n < 25) {
    bnl=0.0;
    for (j=1;j<=n;j++)
      if (uniform() < p) ++bnl;
  } else if (am < 1.0) {
    g=exp(-am);
    t=1.0;
    for (j=0;j<=n;j++) {
      t *= uniform();
      if (t < g) break;
    }
    bnl=(j <= n ? j : n);
  } else {
    if (n != nold) {
      en=n;
      oldg=gammln(en+1.0);
      nold=n;
    } if (p != pold) {
      pc=1.0-p;
      plog=log(p);
      pclog=log(pc);
      pold=p;
    }
    sq=sqrt(2.0*am*pc);
    do {
      do {
        angle=M_PI*uniform();
	y=tan(angle);
	em=sq*y+am;
      } while (em < 0.0 || em >= (en+1.0));
      em=floor(em);
      t=1.2*sq*(1.0+y*y)*exp(oldg-gammln(em+1.0)
        -gammln(en-em+1.0)+em*plog+(en-em)*pclog);
    } while (uniform() > t);
    bnl=em;
  }
  if (p != pp) bnl=n-bnl;
  return bnl;
}



}