EMAN2
randnum.cpp
Go to the documentation of this file.
00001 
00005 /*
00006 * Author: Grant Tang, 03/10/2008 (gtang@bcm.edu)
00007 * Copyright (c) 2000-2006 Baylor College of Medicine
00008 * 
00009 * This software is issued under a joint BSD/GNU license. You may use the
00010 * source code in this file under either license. However, note that the
00011 * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012 * so you are responsible for compliance with the licenses of these packages
00013 * if you opt to use BSD licensing. The warranty disclaimer below holds
00014 * in either instance.
00015 * 
00016 * This complete copyright notice must be included in any revised version of the
00017 * source code. Additional authorship citations may be added, but existing
00018 * author citations must be preserved.
00019 * 
00020 * This program is free software; you can redistribute it and/or modify
00021 * it under the terms of the GNU General Public License as published by
00022 * the Free Software Foundation; either version 2 of the License, or
00023 * (at your option) any later version.
00024 * 
00025 * This program is distributed in the hope that it will be useful,
00026 * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028 * GNU General Public License for more details.
00029 * 
00030 * You should have received a copy of the GNU General Public License
00031 * along with this program; if not, write to the Free Software
00032 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
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                         //printf("Got seed %u from gettimeofday()\n",seed);
00068                 } 
00069                 else {
00070                         fread(&seed,sizeof(seed),1,devrandom);
00071                         //printf("Got seed %u from /dev/random\n",seed);
00072                         fclose(devrandom);
00073                 }
00074         
00075                 return seed;
00076         }
00077 #else
00078         /* this function works on Windows */
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         /*converting file time to unix epoch*/
00092         tmpres /= 10;  /*convert into microseconds*/
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 }