00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include <cmath>
00037 #include <ctime>
00038 #include <cstdio>
00039 #ifndef _WIN32
00040 #include <sys/time.h>
00041 #else
00042 #include <Windows.h>
00043 #if defined(_MSC_VER) || defined(_MSC_EXTENSIONS)
00044 #define DELTA_EPOCH_IN_MICROSECS 11644473600000000Ui64
00045 #else
00046 #define DELTA_EPOCH_IN_MICROSECS 11644473600000000ULL
00047 #endif
00048 #endif
00049 #include "randnum.h"
00050
00051 using namespace EMAN;
00052
00053 namespace {
00054 #ifndef _WIN32
00055
00058 unsigned long long random_seed()
00059 {
00060 unsigned int seed;
00061 struct timeval tv;
00062 FILE *devrandom;
00063
00064 if ((devrandom = fopen("/dev/urandom","r")) == NULL) {
00065 gettimeofday(&tv,0);
00066 seed = tv.tv_sec + tv.tv_usec;
00067
00068 }
00069 else {
00070 fread(&seed,sizeof(seed),1,devrandom);
00071
00072 fclose(devrandom);
00073 }
00074
00075 return seed;
00076 }
00077 #else
00078
00079 unsigned long long random_seed()
00080 {
00081 FILETIME ft;
00082 unsigned __int64 tmpres = 0;
00083 static int tzflag;
00084 struct timeval tv;
00085
00086 GetSystemTimeAsFileTime(&ft);
00087 tmpres |= ft.dwHighDateTime;
00088 tmpres <<= 32;
00089 tmpres |= ft.dwLowDateTime;
00090
00091
00092 tmpres /= 10;
00093 tmpres -= DELTA_EPOCH_IN_MICROSECS;
00094 tv.tv_sec = (long long)(tmpres / 1000000UL);
00095 tv.tv_usec = (long long)(tmpres % 1000000UL);
00096
00097 unsigned long long seed = tv.tv_sec + tv.tv_usec;
00098 return seed;
00099 }
00100 #endif
00101 }
00102
00103 Randnum * Randnum::_instance = 0;
00104 const gsl_rng_type * Randnum::T = gsl_rng_default;
00105 gsl_rng * Randnum::r = 0;
00106 unsigned long long Randnum::_seed = random_seed();
00107
00108 Randnum * Randnum::Instance() {
00109 if(_instance == 0) {
00110 _instance = new Randnum();
00111 }
00112
00113 return _instance;
00114 }
00115
00116 Randnum * Randnum::Instance(const gsl_rng_type * _t) {
00117 if(_instance == 0) {
00118 _instance = new Randnum(_t);
00119 }
00120 else if(_t != _instance->T) {
00121 gsl_rng_free (_instance->r);
00122 _instance->r = gsl_rng_alloc (_t);
00123 gsl_rng_set(_instance->r, _seed );
00124 }
00125
00126 return _instance;
00127 }
00128
00129 Randnum::Randnum()
00130 {
00131 r = gsl_rng_alloc (T);
00132 gsl_rng_set(r, _seed );
00133 }
00134
00135 Randnum::Randnum(const gsl_rng_type * _t)
00136 {
00137 r = gsl_rng_alloc (_t);
00138 gsl_rng_set(r, _seed );
00139 }
00140
00141 Randnum::~Randnum()
00142 {
00143 gsl_rng_free (r);
00144 }
00145
00146 void Randnum::set_seed(unsigned long long seed)
00147 {
00148 _seed = seed;
00149 gsl_rng_set(r, _seed);
00150 }
00151
00152 unsigned long long Randnum::get_seed()
00153 {
00154 return _seed;
00155 }
00156
00157 long long Randnum::get_irand(long long lo, long long hi) const
00158 {
00159 return gsl_rng_uniform_int(r, hi-lo+1)+lo;
00160 }
00161
00162 float Randnum::get_frand(double lo, double hi) const
00163 {
00164 return static_cast<float>(gsl_rng_uniform(r) * (hi -lo) + lo);
00165 }
00166
00167 float Randnum::get_frand_pos(double lo, double hi) const
00168 {
00169 return static_cast<float>(gsl_rng_uniform_pos(r) * (hi -lo) + lo);
00170 }
00171
00172 float Randnum::get_gauss_rand(float mean, float sigma) const
00173 {
00174 float x = 0.0f;
00175 float y = 0.0f;
00176 float r = 0.0f;
00177 bool valid = true;
00178
00179 while (valid) {
00180 x = get_frand_pos(-1.0, 1.0);
00181 y = get_frand_pos(-1.0, 1.0);
00182 r = x * x + y * y;
00183
00184 if (r <= 1.0 && r > 0) {
00185 valid = false;
00186 }
00187 }
00188
00189 float f = std::sqrt(-2.0f * std::log(r) / r);
00190 float result = x * f * sigma + mean;
00191 return result;
00192 }
00193
00194 void Randnum::print_generator_type() const
00195 {
00196 const gsl_rng_type **t, **t0;
00197
00198 t0 = gsl_rng_types_setup ();
00199
00200 printf ("Available generators:\n");
00201
00202 for (t = t0; *t != 0; t++) {
00203 printf ("%s\n", (*t)->name);
00204 }
00205 }