random_number_generators.hpp

Go to the documentation of this file.

#ifndef RANDOM_NUMBER_GENERATORS_HPP_INCLUDED
#define RANDOM_NUMBER_GENERATORS_HPP_INCLUDED

// Standard Library
#include <cmath>    // pow
#include <cstdlib>  // abort
#include <limits>   // 
#include <random>   // mt19937, uniform_real_distribution, normal_distribution
#include <vector>


namespace rng
{

  // ==============================================================================================
  double powerlaw_rand(std::vector<double> &param, std::mt19937& _engine)
  {
    std::uniform_real_distribution<double> uniform_01;
    return param[1] * std::pow(1-( uniform_01(_engine) ),-1./(param[0]-1));
  }

  // ==============================================================================================
  double powerlaw_m1(std::vector<double> &param)
  {
    if(param[0] > 2)
    {
      return ( ( param[0]-1 ) * param[1] ) / ( param[0]-2 );
    }
    else
    {
      return std::numeric_limits<double>::infinity();
      std::cout << "First moment of power-law distribution infinite (gamma < 2)." << std::endl;
      // abort();
    }
  }

  // ==============================================================================================
  double powerlaw_hardcutoff_rand(std::vector<double> &param, std::mt19937& _engine)
  {
    std::uniform_real_distribution<double> uniform_01;
    return param[1] * std::pow(1-( uniform_01(_engine) * (1-std::pow(param[2]/param[1],(1.-param[0])))),-1./(param[0]-1));
  }

  // ==============================================================================================
  double powerlaw_hardcutoff_m1(std::vector<double> &param)
  {
    return ( ( param[0]-1 ) * param[1] * ( 1 - std::pow(param[2]/param[1],2-param[0]) ) ) / ( ( param[0]-2 ) * ( 1 - std::pow(param[2]/param[1],1-param[0]) ) );
  }

  // ==============================================================================================
  double exponential_rand(std::vector<double> &param, std::mt19937& _engine)
  {
    std::uniform_real_distribution<double> uniform_01;
    return -param[0] * std::log(1-uniform_01(_engine));
  }

  // ==============================================================================================
  double exponential_m1(double param)
  {
    return param;
  }

  // ==============================================================================================
  double gamma_rand(std::vector<double> &param, std::mt19937& _engine)
  {
    std::uniform_real_distribution<double> uniform_01;
    std::normal_distribution<double> normal_01(0.0,1.0);
    double d = param[0] - 1./3.;
    double lhs = 5.;
    double rhs = 1.;
    double x,v;
    while( lhs > rhs )
    {
      lhs = std::log( uniform_01(_engine) );
      x = normal_01(_engine);
      v = std::pow(1+x/(std::sqrt(9.*d)),3);
      rhs = std::pow(x,2)/2. + d - d*v + d*std::log(v);
    }
    return param[1]*d*v;
  }

  // ==============================================================================================
  double gamma_m1(std::vector<double> param)
  {
    return param[0]*param[1];
  }
 // End group RandomNumberGenerators

} // End of namespace rng.



#endif // RANDOM_NUMBER_GENERATORS_HPP_INCLUDED