dlib C++ Library - rand_kernel_1.h

// Copyright (C) 2007 Davis E. King (davis@dlib.net)
// License: Boost Software License See LICENSE.txt for the full license.
#ifndef DLIB_RAND_KERNEl_1_
#define DLIB_RAND_KERNEl_1_
#include <string>
#include <complex>
#include "../algs.h"
#include "rand_kernel_abstract.h"
#include "mersenne_twister.h"
#include "../is_kind.h"
#include <iostream>
#include "../serialize.h"
#include "../string.h"
namespace dlib
{
 class rand
 {
 /*! 
 INITIAL VALUE
 - seed == ""
 CONVENTION
 - the random numbers come from the boost mersenne_twister code
 - get_seed() == seed
 !*/
 
 public:
 // These typedefs are here for backwards compatibility with older versions of dlib.
 typedef rand kernel_1a;
 typedef rand float_1a;
 rand(
 ) 
 {
 init();
 }
 rand (
 time_t seed_value
 )
 {
 init();
 set_seed(cast_to_string(seed_value));
 }
 rand (
 const std::string& seed_value
 )
 {
 init();
 set_seed(seed_value);
 }
 virtual ~rand(
 )
 {}
 void clear(
 )
 {
 mt.seed();
 seed.clear();
 has_gaussian = false;
 next_gaussian = 0;
 // prime the generator a bit
 for (int i = 0; i < 10000; ++i)
 mt();
 }
 
 const std::string& get_seed (
 )
 {
 return seed;
 }
 void set_seed (
 const std::string& value
 )
 {
 seed = value;
 // make sure we do the seeding so that using a seed of "" gives the same
 // state as calling this->clear()
 if (value.size() != 0)
 {
 uint32 s = 0;
 for (std::string::size_type i = 0; i < seed.size(); ++i)
 {
 s = (s*37) + static_cast<uint32>(seed[i]);
 }
 mt.seed(s);
 }
 else
 {
 mt.seed();
 }
 // prime the generator a bit
 for (int i = 0; i < 10000; ++i)
 mt();
 has_gaussian = false;
 next_gaussian = 0;
 }
 unsigned char get_random_8bit_number (
 )
 {
 return static_cast<unsigned char>(mt());
 }
 uint16 get_random_16bit_number (
 )
 {
 return static_cast<uint16>(mt());
 }
 inline uint32 get_random_32bit_number (
 )
 {
 return mt();
 }
 inline uint64 get_random_64bit_number (
 )
 {
 const uint64 a = get_random_32bit_number();
 const uint64 b = get_random_32bit_number();
 return (a<<32)|b;
 }
 double get_double_in_range (
 double begin,
 double end
 )
 {
 DLIB_ASSERT(begin <= end);
 return begin + get_random_double()*(end-begin);
 }
 long long get_integer_in_range(
 long long begin,
 long long end
 )
 {
 DLIB_ASSERT(begin <= end);
 if (begin == end)
 return begin;
 auto r = get_random_64bit_number();
 const auto limit = std::numeric_limits<decltype(r)>::max();
 const auto range = end-begin;
 // Use rejection sampling to remove the biased sampling you would get with
 // the naive get_random_64bit_number()%range sampling. 
 while(r >= (limit/range)*range)
 r = get_random_64bit_number();
 return begin + static_cast<long long>(r%range);
 }
 long long get_integer(
 long long end
 )
 {
 DLIB_ASSERT(end >= 0);
 return get_integer_in_range(0,end);
 }
 double get_random_double (
 )
 {
 uint32 temp;
 temp = rand::get_random_32bit_number();
 temp &= 0xFFFFFF;
 double val = static_cast<double>(temp);
 val *= 0x1000000;
 temp = rand::get_random_32bit_number();
 temp &= 0xFFFFFF;
 val += temp;
 val /= max_val;
 if (val < 1.0)
 {
 return val;
 }
 else
 {
 // return a value slightly less than 1.0
 return 1.0 - std::numeric_limits<double>::epsilon();
 }
 }
 float get_random_float (
 )
 {
 uint32 temp;
 temp = rand::get_random_32bit_number();
 temp &= 0xFFFFFF;
 const float scale = 1.0/0x1000000;
 const float val = static_cast<float>(temp)*scale;
 if (val < 1.0f)
 {
 return val;
 }
 else
 {
 // return a value slightly less than 1.0
 return 1.0f - std::numeric_limits<float>::epsilon();
 }
 }
 
 std::complex<double> get_random_complex_gaussian (
 )
 {
 double x1, x2, w;
 const double rndmax = std::numeric_limits<dlib::uint32>::max();
 // Generate a pair of Gaussian random numbers using the Box-Muller transformation.
 do 
 {
 const double rnd1 = get_random_32bit_number()/rndmax;
 const double rnd2 = get_random_32bit_number()/rndmax;
 x1 = 2.0 * rnd1 - 1.0;
 x2 = 2.0 * rnd2 - 1.0;
 w = x1 * x1 + x2 * x2;
 } while ( w >= 1.0 );
 w = std::sqrt( (-2.0 * std::log( w ) ) / w );
 return std::complex<double>(x1 * w, x2 * w);
 }
 double get_random_gaussian (
 )
 {
 if (has_gaussian)
 {
 has_gaussian = false;
 return next_gaussian;
 }
 
 std::complex<double> r = get_random_complex_gaussian();
 next_gaussian = r.imag();
 has_gaussian = true;
 return r.real();
 }
 double get_random_exponential (
 double lambda
 )
 {
 DLIB_ASSERT(lambda > 0, "lambda must be greater than zero");
 double u = 0.0;
 while (u == 0.0)
 u = get_random_double();
 return -std::log( u ) / lambda;
 }
 double get_random_weibull (
 double lambda,
 double k,
 double gamma
 )
 {
 DLIB_ASSERT(k > 0, "k must be greater than zero");
 DLIB_ASSERT(lambda > 0, "lambda must be greater than zero");
 double u = 0.0;
 while (u == 0.0)
 u = get_random_double();
 return gamma + lambda*std::pow(-std::log(u), 1.0 / k);
 }
 double get_random_beta (
 double alpha,
 double beta
 )
 {
 DLIB_CASSERT(alpha > 0, "alpha must be greater than zero")
 DLIB_CASSERT(beta > 0, "beta must be greater than zero");
 auto u = std::pow(get_random_double(), 1 / alpha);
 auto v = std::pow(get_random_double(), 1 / beta);
 while ((u + v) > 1 || (u == 0 && v == 0))
 {
 u = std::pow(get_random_double(), 1 / alpha);
 v = std::pow(get_random_double(), 1 / beta);
 }
 return u / (u + v);
 }
 void swap (
 rand& item
 )
 {
 exchange(mt,item.mt);
 exchange(seed, item.seed);
 exchange(has_gaussian, item.has_gaussian);
 exchange(next_gaussian, item.next_gaussian);
 }
 
 friend void serialize(
 const rand& item, 
 std::ostream& out
 );
 friend void deserialize(
 rand& item, 
 std::istream& in 
 );
 private:
 void init()
 {
 // prime the generator a bit
 for (int i = 0; i < 10000; ++i)
 mt();
 max_val = 0xFFFFFF;
 max_val *= 0x1000000;
 max_val += 0xFFFFFF;
 max_val += 0.05;
 has_gaussian = false;
 next_gaussian = 0;
 }
 mt19937 mt;
 std::string seed;
 double max_val;
 bool has_gaussian;
 double next_gaussian;
 };
 inline void swap (
 rand& a, 
 rand& b 
 ) { a.swap(b); } 
 template <>
 struct is_rand<rand>
 {
 static const bool value = true; 
 };
 inline void serialize(
 const rand& item, 
 std::ostream& out
 )
 {
 int version = 1;
 serialize(version, out);
 serialize(item.mt, out);
 serialize(item.seed, out);
 serialize(item.has_gaussian, out);
 serialize(item.next_gaussian, out);
 }
 inline void deserialize(
 rand& item, 
 std::istream& in 
 )
 {
 int version;
 deserialize(version, in);
 if (version != 1)
 throw serialization_error("Error deserializing object of type rand: unexpected version."); 
 deserialize(item.mt, in);
 deserialize(item.seed, in);
 deserialize(item.has_gaussian, in);
 deserialize(item.next_gaussian, in);
 }
}
#endif // DLIB_RAND_KERNEl_1_