STK_RandBase.h

Go to the documentation of this file.

/*--------------------------------------------------------------------*/
/*     Copyright (C) 2004-2016  Serge Iovleff, Université Lille 1, Inria
 
    This program 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 program 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 program; if not, write to the
    Free Software Foundation, Inc.,
    59 Temple Place,
    Suite 330,
    Boston, MA 02111-1307
    USA
 
    Contact : S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
*/
 
/* Project:  stkpp::STatistiK
 * Purpose:  Main pseudo-random uniform, Gaussian and exponential generators.
 * Author:   Serge Iovleff, S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
 **/
 
#ifndef STK_RANDBASE_H
#define STK_RANDBASE_H
 
#ifdef IS_RTKPP_LIB /* is rtkpp lib */
 
#include <R_ext/Random.h>
#include <Rmath.h>
#include <Rcpp/sugar/undoRmath.h>
 
#include <STKernel.h>
 
namespace STK
{
class RandBase
{
  public:
    inline RandBase(){}
    inline ~RandBase(){}
    inline Real randUnif()
    { GetRNGstate();
      Real s = Rf_runif(0, 1);
      PutRNGstate();
      return s;
    }
    inline Real operator()() { return randUnif(); }
    inline Real randGauss( Real const& mu = 0, Real const& sigma = 1)
    { GetRNGstate();
      Real s = Rf_rnorm(mu, sigma);
      PutRNGstate();
      return s;
    }
    inline Real rand( Real const& n )
    { GetRNGstate();
      Real s = Rf_runif(0, n);
      PutRNGstate();
      return s;
    }
    inline Real randDblExc( Real const& n ) { return rand(n);}
};
} // namespace STK
 
 
 
#else
 
// MersenneTwister header.
#include "MersenneTwister.h"
 
namespace STK
{
 
class RandBase: protected MTRand
{
  public:
      RandBase( Real const& glimit = 3.442619855899
              , Real const& gvol   = 9.91256303526217e-3
              , int gsize  = 128
              );
    RandBase( int oneSeed
            , Real const& glimit = 3.442619855899
            , Real const& gvol   = 9.91256303526217e-3
            , int gsize  = 128
            );
 
    template< class TContainer1D>
    RandBase( TContainer1D const& bigSeed
            , Real const& glimit = 3.442619855899
            , Real const& gvol   = 9.91256303526217e-3
            , int gsize  = 128
            )
            : gsize_(gsize), glimit_(glimit), gvol_(gvol)
    {
      // dimension
      int first = bigSeed.begin(), size = bigSeed.size();
      uint32* arraySeed = new uint32[size];
      // cast int in uint32
      for (int i=first; i<bigSeed.end(); i++)
      { arraySeed[i-first] = uint32(bigSeed[i]);}
      // Re-seeding functions with same behavior as initializers
      seed(arraySeed, size);
      delete[] arraySeed;
      // initialize zigourat method for gaussian random generator
      gaussInit();
    }
    ~RandBase();
    inline int randDiscreteUnif() { return int(Real(randInt()));}
    inline Real randUnif() { return Real(randDblExc());}
 
    inline Real operator()() { return (Real)MTRand::operator()(); }
    inline Real rand() { return (Real)MTRand::rand(); }
    inline Real rand( Real const& n ) { return (Real)MTRand::rand((double)n); }
    inline Real randExc() { return (Real)MTRand::randExc(); }
    inline Real randExc( Real const& n ) { return (Real)MTRand::randExc((double)n); }
    inline Real randDblExc() { return (Real)MTRand::randDblExc(); }
    inline Real randDblExc( Real const& n ) { return (Real)MTRand::randDblExc((double)n); }
 
    Real randGauss(Real const& mu = 0, Real const& sigma = 1);
    Real randExp();
 
  private:
    /*  Gauss parameters. */
    const int gsize_;
    Real const   glimit_;
    Real const   gvol_;
    Real *kn, *wn, *fn;
    void gaussInit();
};
 
 
/* Initialize with a simple int seed. */
inline RandBase::RandBase( int oneSeed, Real const& glimit, Real const& gvol, int gsize)
                         : MTRand( (uint32)oneSeed ), gsize_(gsize), glimit_(glimit), gvol_(gvol)
{ gaussInit();}
 
/* auto-initialize with /dev/urandom or time() and clock() */
inline RandBase::RandBase( Real const& glimit, Real const& gvol, int gsize)
                         : MTRand(), gsize_(gsize), glimit_(glimit), gvol_(gvol)
{ gaussInit();}
 
/* destructor. */
inline RandBase::~RandBase()
{
  delete [] kn;
  delete [] wn;
  delete [] fn;
}
 
/* Pseudo-random gaussian generator.
 * \f[     f(x) = \frac{1}{\sqrt{2\pi}}
 *         \exp\left\{-\frac{x^2}{2}\right\}
 * \f]
 *  Return a real number from a normal (Gaussian) normalized
 *  distribution.
**/
inline Real RandBase::randGauss(Real const& mu, Real const& sigma)
{
  while(1)
  {
    // uniforms number
    Real u = 2.0 * randUnif() - 1.0;
    // random box
    uint32 i = randInt() & 0x7F; //
    // squeeze step : try the rectangular boxes
    if (std::abs(u) < wn[i]) return mu + sigma * u * kn[i];
    // bottom box: we have to sample from the tail
    if (i == 0)
    {
      Real x;
      // loop
      do { x = randExp() / glimit_;}
      while ( 2.0 * randExp() < x * x);
      // result sign(u) * (x+dr)
      return mu + sigma * sign(u, x+glimit_);
    }
    else // other box
    { // x : result, f1 : intermediary result
      Real x = u * kn[i], f1 = fn[i+1];
      // reject step : is this a sample from the wedges ?
      if ( f1 + randUnif()*( fn[i]*exp(x*x/2.)-f1) < 1.0)
      { return mu + sigma * x;}
    }
  }
  // avoid warning at compilation
  return 0.;
}
 
/* Pseudo-random exponential generator.
 * \f[
 *      f(x) = \exp\left\{ -x \right\}
 * \f]
 *  Return a real number from an exponential normalized
 *  distribution.
**/
inline Real RandBase::randExp()
{ return -Real(std::log(randUnif()));}
 
/* Initialization (Zigourrat method).  */
inline void RandBase::gaussInit()
{
  kn = new Real[gsize_+1];
  wn = new Real[gsize_];
  fn = new Real[gsize_+1];
 
  Real f = exp(-0.5 * glimit_ * glimit_);
  kn[0]  = gvol_ / f; // [0] is bottom block: gvol/f(gbound)
 
  kn[gsize_] = 0;
  fn[gsize_] = 1;
 
  kn[1]  = glimit_;
  fn[1]  = f;
 
  for (int i = 2; i < gsize_; ++i)
  {
    kn[i] = sqrt(-2 * log(gvol_ / kn[i-1] + f));
    fn[i] = (f = exp(-0.5 * kn[i] * kn[i]));
  }
  for (int i = 0; i < gsize_; ++i)
    wn[i] = kn[i + 1] / kn[i];
}
 
} // namespace STK
 
#endif
 
#endif //STK_RANDBASE_H