Random.hpp

Go to the documentation of this file.

/****************************************************************
 * RANDOM.HPP
 * 
 * This file contains random number generators related functions.
 * 
 * Authors: J. Barthelemy
 * Date   : 17 july 2012
 ****************************************************************/

#ifndef RANDOM_HPP_
#define RANDOM_HPP_

#include <iostream>
#include <vector>
#include <cmath>
#include <stdexcept>


typedef struct draw_2d draw_2d;
struct draw_2d {

  float x1;  
  float x2;  

  draw_2d(float aX1, float aX2) : x1(aX1), x2(aX2) {};

};



typedef struct dist_param dist_param;
struct dist_param {
  float mu;    
  float sigma; 
  float max;   
};



typedef struct dist_param_mixture dist_param_mixture;
struct dist_param_mixture {
  std::vector<float> mu;    
  std::vector<float> sigma; 
  std::vector<float> p;     
  float              max;   
  float              min;   
};



typedef struct dist_param_2d dist_param_2d;
struct dist_param_2d {
  float mu[2];    
  float sigma[3]; 
};



typedef struct dist_param_mixture_2d dist_param_mixture_2d;
struct dist_param_mixture_2d {
  std::vector<dist_param_2d > components;  
  std::vector<float> p;                    
  float min[2];                            
  float max[2];                            
};


struct Ranq1 {

  unsigned long long v;  


  Ranq1( unsigned long long j) : v(4101842887655102017LL) {
    v ^= j;
    v = int64();
  }


  inline unsigned long long int64() {
    v ^= v >> 21;
    v ^= v << 35;
    v ^= 4;
    return v * 2685821657736338717LL;
  }


  inline double doub() {
    return 5.42101086242752217E-20 * int64();
  }


  inline float fl() {
    return (float)doub();
  }


  inline unsigned int int32() {
    return (unsigned int)int64();
  }

};



struct Ranfib {

  double dtab[55];      
  double dd;            
  int inext;            
  int inextp;           


  Ranfib( unsigned long j ) : inext(0), inextp(31) {

    dd = 0.0;

    // Using Ranq1 for initialization.
    Ranq1 init(j);

    for( int k = 0; k < 55; k++) {
      dtab[k] = init.doub();
    }

  }


  inline double doub() {

    if( ++inext  == 55 ) inext  = 0;
    if( ++inextp == 55 ) inextp = 0;
    dd = dtab[inext] - dtab[inextp];
    if( dd < 0 ) dd += 1.0;
    return ( dtab[inext] = dd );

  }


  inline float fl() {
    return (float)doub();
  }

};

struct Normaldev : Ranfib {


  Normaldev(unsigned long long i) : Ranfib(i) {};


  inline float dev(double mu, double sigma) {

    float u, v, x, y, q;

    do {
      u = fl();
      v = 1.7156 * ( fl() - 0.5 );
      x = u - 0.449871;
      y = abs(v) + 0.386595;
      q = x * x + y * ( 0.19600 * y - 0.25472 * x );
    } while ( q > 0.27597 && ( q > 0.27846 || v * v > -4.0 * log(u) * u * u ) );

    return mu + sigma * v / u;

  }


  inline float dev(double mu, double sigma, float max) {

      float result;

      do {
        result = dev(mu, sigma);
      } while ( result > max );
      return result;

  }


  inline float dev(double mu, double sigma, float min, float max) {

    float result;

    do {
      result = dev(mu, sigma);
    } while ( result > max || result < min );
    return result;

  }


};

struct LogNormaldev : Ranfib {


  LogNormaldev(unsigned long long i) : Ranfib(i) {};


  inline float dev(double mu, double sigma) {

    float u, v, x, y, q;

    // normal draw
    do {
      u = fl();
      v = 1.7156 * ( fl() - 0.5 );
      x = u - 0.449871;
      y = abs(v) + 0.386595;
      q = x * x + y * ( 0.19600 * y - 0.25472 * x );
    } while ( q > 0.27597 && ( q > 0.27846 || v * v > -4.0 * log(u) * u * u ) );

    return exp(mu + sigma * v / u);

  }


  inline float dev(double mu, double sigma, float max) {

    float result;

    do {
      result = dev(mu, sigma);
    } while ( result > max );
    return result;

  }


  inline float dev(double mu, double sigma, float min, float max) {

    float result;

    do {
      result = dev(mu, sigma);
    } while ( result > max || result < min );
    return result;

  }


};

struct MixtureNormal : Normaldev {

  MixtureNormal(unsigned long long i) : Normaldev(i) {};


  inline float dev(std::vector<float> mu, std::vector<float> sigma, std::vector<float> p) {

    float  a_p, prop_cum;
    int    comp;
    float u, v, x, y, q;

    // looking for the right component of the mixture
    a_p      = fl();
    comp     = 0;
    prop_cum = p[comp];
    while( prop_cum < a_p ) {
      comp++;                               // moving to the next component
      prop_cum = prop_cum + p[comp];        // computation of the cumulative proportion
    }

    // normal draw
    do {
      u = fl();
      v = 1.7156 * ( fl() - 0.5 );
      x = u - 0.449871;
      y = abs(v) + 0.386595;
      q = x * x + y * ( 0.19600 * y - 0.25472 * x );
    } while ( q > 0.27597 && ( q > 0.27846 || v * v > -4.0 * log(u) * u * u ) );    

    return (mu[comp] + sigma[comp] * v / u);

  }


  inline float dev( std::vector<float> mu, std::vector<float> sigma, std::vector<float> p, float max) {

    float result;

    do {
      result = dev(mu, sigma, p);
      std::cout << "TEST MIXTURE NORM " << result << std::endl;
    } while ( result > max );
    return result;

  }


  inline float dev( std::vector<float> mu, std::vector<float> sigma, std::vector<float> p, float min, float max) {

    float result;

    do {
      result = dev(mu, sigma, p);
      //std::cout << "TEST MIXTURE NORM " << result << " MIN "<< min << " MAX " << max << std::endl;
    } while ( result > max || result < min );
    return result;

  }

};

struct MixtureLogNormal : LogNormaldev {

  MixtureLogNormal(unsigned long long i) : LogNormaldev(i) {};


  inline float dev(std::vector<float> mu, std::vector<float> sigma, std::vector<float> p) {

    float  a_p, prop_cum;
    int    comp;
    float u, v, x, y, q;

    // looking for the right component of the mixture
    a_p      = fl();
    comp     = 0;
    prop_cum = p[comp];
    while( prop_cum < a_p ) {
      comp++;                               // moving to the next component
      prop_cum = prop_cum + p[comp];        // computation of the cumulative proportion
    }

    // log-normal draw
    do {
      u = doub();
      v = 1.7156 * ( fl() - 0.5 );
      x = u - 0.449871;
      y = abs(v) + 0.386595;
      q = x * x + y * ( 0.19600 * y - 0.25472 * x );
    } while ( q > 0.27597 && ( q > 0.27846 || v * v > -4.0 * log(u) * u * u ) );

    return exp(mu[comp] + sigma[comp] * v / u);

  }


  inline float dev( std::vector<float> mu, std::vector<float> sigma, std::vector<float> p, float max) {

    float result;

    do {
      result = dev(mu, sigma, p);
    } while ( result > max );
    return result;

  }


  inline float dev( std::vector<float> mu, std::vector<float> sigma, std::vector<float> p, float min, float max) {

    float result;

    do {
      result = dev(mu, sigma, p);
    } while ( result > max || result < min );
    return result;

  }

};


struct MixtureLogNormal2D : Normaldev {

  MixtureLogNormal2D(unsigned long long i) : Normaldev(i) {};

  inline draw_2d dev(dist_param_mixture_2d distrib) {

    float   a_p, prop_cum;
    int     comp;
    draw_2d result(0.0,0.0);

    // looking for the right component of the mixture
    a_p      = fl();
    comp     = 0;
    prop_cum = distrib.p[comp];
    while( prop_cum < a_p ) {
      comp++;                                       // moving to the next component
      prop_cum = prop_cum + distrib.p[comp];        // computation of the cumulative proportion
    }

    do {

      result.x1 = Normaldev::dev(0,1);
      result.x2 = Normaldev::dev(0,1);

      // generating the draws
      result.x1 = exp(distrib.components[comp].mu[0] + distrib.components[comp].sigma[0] * result.x1 + distrib.components[comp].sigma[1] * result.x2);
      result.x2 = exp(distrib.components[comp].mu[1] + distrib.components[comp].sigma[2] * result.x2);

    } while (result.x1 > distrib.max[0] || result.x2 > distrib.max[1] );

    return result;

  }

};

template <typename T>
class SingletonRnd {

protected:

  SingletonRnd () {}

  ~SingletonRnd () {}

public:


  static void makeInstance ( unsigned long long i ) {
    if (nullptr == _singleton)  {
      _singleton = new T(i);
    }
  }


  static T *getInstance () {
    if (nullptr == _singleton)    {
      std::cerr << "No instance already created!" << std::endl;
    }
    return (static_cast<T*> (_singleton));
  }

  static void kill () {
    if (nullptr != _singleton) {
      delete _singleton;
      _singleton = nullptr;
    }
  }

private:

  static T *_singleton;    

};

template <typename T>
T *SingletonRnd<T>::_singleton = nullptr;


class RandomGenerators : public SingletonRnd<RandomGenerators> {

  friend class SingletonRnd<RandomGenerators>;

public:

  Ranq1              unif;                
  Ranfib             fast_unif;           
  Normaldev          norm_dev;            
  LogNormaldev       lognorm_dev;         
  MixtureNormal      mixt_norm_dev;       
  MixtureLogNormal   mixt_lognorm_dev;    
  MixtureLogNormal2D mixt_lognorm_dev_2d; 


  RandomGenerators(unsigned long long i) : unif(i + 10000), fast_unif(i+10), norm_dev(i + 100), lognorm_dev(i + 1000),
                                           mixt_norm_dev(i + 10000), mixt_lognorm_dev(i + 100000), mixt_lognorm_dev_2d(i + 1000000) {};

  virtual ~RandomGenerators() {};

};


template <typename T> int draw_discrete(std::vector<T> freq) {

  if( freq.size() < 1 ) throw std::logic_error( "draw_discrete: empty vector!" );

  int    N = freq.size();                                                        // size of freq
  double cumfreq [N];                                                            // cumulative distribution
  double draw;                                                                   // random draw
  double icum;                                                                   // cumulative total
  double low;                                                                    // lower bound
  int    i;

  // Build the cumulative distribution.

  icum  = 0.0;
  for( i = 0; i < N; i++ ) {
    icum=icum+freq[i];
    cumfreq[i]=icum;
  }

  // ...  distribution is non-trivial
  if(icum>0.0) {
    for(i=0;i<N;i++) {
      cumfreq[i]=cumfreq[i]/icum;
    }
  }
  // ... distribution is trivial
  else {
    for(i=0;i<N;i++) {
      cumfreq[i]=(i+1.0)/(N+1.0);
    }
  }

  //  Draw a random number between 0 and 1.

  draw = RandomGenerators::getInstance()->fast_unif.doub();

  //  Find the interval in which this number falls
  //  and return the associated class identifier.

  low  = 0.0;
  for ( i = 0; i < N; i++ ) {
    if ( (low <= draw) && (draw <= cumfreq[i]) ) {
      return i;
    }
    low = cumfreq[i];
  }

  std::cerr << "Could not find a class identifier!" << std::endl;
  return -1;

}



template <typename T> int draw_discrete( std::vector<T> freq, int lo_bound, int up_bound ) {

  int   N = ( up_bound - lo_bound ) + 1;                                        // compute the size of the interval
  double cumfreq [N];                                                            // cumulative distribution
  double draw;                                                                   // random draw
  double icum;
  double low;
  int   i, j;

  // Build the cumulative distribution.

  icum  = 0;
  j     = 0;
  for( i = lo_bound; i <= up_bound; i++ ) {
    icum=icum+freq[i];
    cumfreq[j]=icum;
    j++;
  }

  // ... distribution is non-trivial
  if(icum>0) {
    for(i = 0; i < N; i++ ) {
      cumfreq[i]=cumfreq[i]/icum;
    }
  }
  // ... distribution is trivial
  else {
    std::cout << "Error: trivial distribution" << std::endl;
    for( i = 0; i < N; i++ ) {
      cumfreq[i]=(i+1)/(N+1);
    }
  }

  //  Draw a random number between 0 and 1.

  draw = RandomGenerators::getInstance()->fast_unif.doub();

  //  Find the interval in which this number falls
  //  and return the associated class identifier.

  low  = 0.0;
  for ( i = 0; i < N; i++ ) {
    if ( (low <= draw) && (draw <= cumfreq[i]) ) {
      return i + lo_bound;
    }
    low = cumfreq[i];
  }

  std::cerr << "Could not find a class identifier!" << std::endl;
  return -1;

}

// Functions below are depreciated!

/*
 \param mu mean of the distribution
 \param sigma standard deviation of the distribution

 \return a float randomly generated
 */
float norm_rand( float mu, float sigma );

/*
 \param mu mean of the distribution
 \param sigma standard deviation of the distribution
 \param max the upper bound of the distribution

 \return a float randomly generated
 */

float norm_rand( float mu, float sigma, float max );

/*
 \param mu mean of the distribution
 \param sigma standard deviation of the distribution

 \return a float randomly generated
 */
float log_rand( float mu, float sigma );

/*
 \param mu mean of the distribution
 \param sigma standard deviation of the distribution
 \param max the upper bound of the distribution

 \return a float randomly generated
 */
float log_rand( float mu, float sigma, float max );


#endif /* RANDOM_HPP_ */