/*****************************************************************************
 This code is available for academic use
 under the LESSER GENERAL PUBLIC LICENSE 

 Weight Perturbation: enhancing local search strategies
 by perturbing the weights of training instances.
 Copyright (C) 2002  Gal Elidan, Matan Ninio, Nir Friedman and Dale Schuurmans

 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.1 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
 
 Please cite using this refrence:
 
@incollection{Elidan+al:2002,
   author = "Gal Elidan and Matan Ninio and Nir Friedman and Dale Schuurmans",
   booktitle = "Proc. National Conference on Artificial Intelligence (AAAI-02)",
   pages = "132-139",
   year = "2002",
   title = "Data Perturbation for Escaping Local Maxima in Learning",
 }
 
 You can contact the authors at annealing@cs.huji.ac.il
 
*****************************************************************************/
#include "RandomProb.h"
#include <iostream>
#include <stdlib.h>
#include <assert.h>
#include <math.h>

// you may get a warning like 
// "warning: decimal constant is so large that it is  unsigned"
// on the followoing line. this is OK, and is due to a bag in Gcc
static unsigned int MAX_INT=4294967291;

tRandomGenerator _RandomProbGenerator;

unsigned long
tMarsagliaGenerator::Next()
{
  register unsigned long hold;
  if (--i==0) i=43;
  if (--j==0) j=43;
  hold = word1[i]+carry;
  if (word1[j] < hold)
    {
      word1[i] = MAX_INT - (hold-word1[j]);
      carry = 1;
    }
  else
    {
      word1[i] = word1[j]-hold;
      carry=0;
    }
  weyl-=362436069;
  return word1[i] - weyl;
}

tMarsagliaGenerator::tMarsagliaGenerator(unsigned query)
{
  unsigned start=1;  // default

  if (query)
    {
      cout << "Enter random number seed, 0<seed<65000." << endl;
      cin >> start;
    }

  srand(start);

  for (j=1;j<=43;j++)        		// make initial numbers.
    {
      word1[j]=rand();
      word1[j]=((word1[j] << 15 ) + rand());
      word1[j]=((word1[j] << 2 ) + rand()%4);
      if (word1[j]>MAX_INT) word1[j] = MAX_INT;
    }

  // initialize markers
  i=44;  j=23;  carry=1; weyl=rand();

  for (unsigned long a=1,garbage; a<100000; a++)       	// Warm-up
    garbage = Next();
}

void
tMarsagliaGenerator::Initialize(unsigned start)
{
  srand(start);

  for (j=1;j<=43;j++)        		// make initial numbers.
    {
      word1[j]=rand();
      word1[j]=((word1[j] << 15 ) + rand());
      word1[j]=((word1[j] << 2 ) + rand()%4);
      if (word1[j]>MAX_INT) word1[j] = MAX_INT;
    }

  i=44;  j=23;  carry=1; weyl=rand();				// initialize markers

  for (unsigned long a=1,garbage; a<100000; a++)       	// Warm-up
    garbage = Next();
}

double tMarsagliaGenerator::RandomDouble(double range)
{
  //  return ((((double)Next())/4294967295)*range);
  return ((((double)Next())/MAX_INT)*range);
}


double
tRandomGenerator::RandomProb(void)
{
  return RandomDouble(1.0);
}

// Code adopted from David Heckerman
//-----------------------------------------------------------
//	DblGammaGreaterThanOne(dblAlpha)
//
//	routine to generate a gamma random variable with unit scale and
//      alpha > 1
//	reference: Ripley, Stochastic Simulation, p.90 
//	Chang and Feast, Appl.Stat. (28) p.290
//-----------------------------------------------------------
double tRandomGenerator::DblGammaGreaterThanOne(double dblAlpha)
{
  double rgdbl[6];

  rgdbl[1] = dblAlpha - 1.0;
  rgdbl[2] = (dblAlpha - (1.0 / (6.0 * dblAlpha))) / rgdbl[1];
  rgdbl[3] = 2.0 / rgdbl[1];
  rgdbl[4] = rgdbl[3] + 2.0;
  rgdbl[5] = 1.0 / sqrt(dblAlpha);

  for (;;)
  {
    double  dblRand1;
    double  dblRand2;
    do
    {
      dblRand1 = RandomProb();
      dblRand2 = RandomProb();

      if (dblAlpha > 2.5)
	dblRand1 = dblRand2 + rgdbl[5] * (1.0 - 1.86 * dblRand1);

    } while (!(0.0 < dblRand1 && dblRand1 < 1.0));

    double dblTemp = rgdbl[2] * dblRand2 / dblRand1;

    if (rgdbl[3] * dblRand1 + dblTemp + 1.0 / dblTemp <= rgdbl[4] ||
	rgdbl[3] * log(dblRand1) + dblTemp - log(dblTemp) < 1.0)
    {
      return dblTemp * rgdbl[1];
    }
  }
  assert(false);
  return 0.0;
} 

/* routine to generate a gamma random variable with unit scale and alpha
< 1

   reference: Ripley, Stochastic Simulation, p.88 */

double
tRandomGenerator::DblGammaLessThanOne(double dblAlpha)
{
  double dblTemp;

  const double	dblexp = exp(1.0);

  for (;;)
  {
    double dblRand0 = RandomProb();
    double dblRand1 = RandomProb();
    if (dblRand0 <= (dblexp / (dblAlpha + dblexp))) 
    {
      dblTemp = pow(((dblAlpha + dblexp) * dblRand0) /
		    dblexp, 1.0 / dblAlpha);
      if (dblRand1 <= exp(-1.0 * dblTemp)) 
	return dblTemp;
    }
    else 
    {
      dblTemp = -1.0 * log((dblAlpha + dblexp) * (1.0 - dblRand0) /
			   (dblAlpha * dblexp)); 
      if (dblRand1 <= pow(dblTemp,dblAlpha - 1.0)) 
	return dblTemp;
    }				
  }
  assert(false);
  return 0.0;
}  /* DblGammaLessThanOne */

// Routine to generate a gamma random variable with unit scale (beta = 1)
double
tRandomGenerator::DblRanGamma(double dblAlpha)
{
  assert(dblAlpha > 0.0);
  if( dblAlpha < 1.0 )
    return DblGammaLessThanOne(dblAlpha);
  else
    if( dblAlpha > 1.0 )
      return DblGammaGreaterThanOne(dblAlpha);

  return -log(RandomProb());
}  /* gamma */