Random.h
Go to the documentation of this file.
00001 
00005 // Random.h: Definition and Implementation of Random Number Distribution Class
00006 //           Ref: Richard Saucier, "Computer Generation of Statistical 
00007 //                Distributions," ARL-TR-2168, US Army Research Laboratory,
00008 //                Aberdeen Proving Ground, MD, 21005-5068, March 2000.
00009 
00010 #ifndef RANDOM_H
00011 #define RANDOM_H
00012 
00013 #include <iostream>
00014 #include <fstream>
00015 #include <vector>
00016 #include <list>
00017 #include <algorithm>
00018 #include <functional>
00019 #include <cassert>
00020 #include <climits>
00021 #include <cmath>
00022 #include <cfloat>     // for FLT_MIN and FLT_MAX
00023 #include <unistd.h>   // for getpid
00024 using namespace std;
00025 
00026 // Constants for Tausworthe random bit generator
00027 // Ref: Tausworthe, Robert C., "Random Numbers Generated by Linear Recurrence
00028 //      Modulo Two," Mathematics of Computation, vol. 19, pp. 201-209, 1965.
00029 
00030 static const unsigned DEGREE_MAX = 32;   // max degree (bits per word)
00031 
00032 static const unsigned BIT[ 1 + DEGREE_MAX ] = {
00033 
00034 // Hexidecimal       Value      Degree   
00035 // -----------       -----      ------
00036    0x00000000,    // 0          0
00037    0x00000001,    // 2^0        1
00038    0x00000002,    // 2^1        2
00039    0x00000004,    // 2^2        3
00040    0x00000008,    // 2^3        4
00041    0x00000010,    // 2^4        5
00042    0x00000020,    // 2^5        6
00043    0x00000040,    // 2^6        7
00044    0x00000080,    // 2^7        8
00045    0x00000100,    // 2^8        9
00046    0x00000200,    // 2^9        10
00047    0x00000400,    // 2^10       11
00048    0x00000800,    // 2^11       12
00049    0x00001000,    // 2^12       13
00050    0x00002000,    // 2^13       14
00051    0x00004000,    // 2^14       15
00052    0x00008000,    // 2^15       16
00053    0x00010000,    // 2^16       17
00054    0x00020000,    // 2^17       18
00055    0x00040000,    // 2^18       19
00056    0x00080000,    // 2^19       20
00057    0x00100000,    // 2^20       21
00058    0x00200000,    // 2^21       22
00059    0x00400000,    // 2^22       23
00060    0x00800000,    // 2^23       24
00061    0x01000000,    // 2^24       25
00062    0x02000000,    // 2^25       26
00063    0x04000000,    // 2^26       27
00064    0x08000000,    // 2^27       28
00065    0x10000000,    // 2^28       29
00066    0x20000000,    // 2^29       30
00067    0x40000000,    // 2^30       31 
00068    0x80000000     // 2^31       32
00069 };
00070 
00071 // Coefficients that define a primitive polynomial (mod 2)
00072 // Ref: Watson, E. J., "Primitive Polynomials (Mod 2),"
00073 //      Mathematics of Computation, vol. 16, pp. 368-369, 1962.
00074 
00075 static const unsigned MASK[ 1 + DEGREE_MAX ] = {
00076 
00077    BIT[0],                                      // 0
00078    BIT[0],                                      // 1
00079    BIT[1],                                      // 2
00080    BIT[1],                                      // 3
00081    BIT[1],                                      // 4
00082    BIT[2],                                      // 5
00083    BIT[1],                                      // 6
00084    BIT[1],                                      // 7
00085    BIT[4] + BIT[3] + BIT[2],                    // 8
00086    BIT[4],                                      // 9
00087    BIT[3],                                      // 10
00088    BIT[2],                                      // 11
00089    BIT[6] + BIT[4] + BIT[1],                    // 12
00090    BIT[4] + BIT[3] + BIT[1],                    // 13
00091    BIT[5] + BIT[3] + BIT[1],                    // 14
00092    BIT[1],                                      // 15
00093    BIT[5] + BIT[3] + BIT[2],                    // 16
00094    BIT[3],                                      // 17
00095    BIT[5] + BIT[2] + BIT[1],                    // 18
00096    BIT[5] + BIT[2] + BIT[1],                    // 19
00097    BIT[3],                                      // 20
00098    BIT[2],                                      // 21
00099    BIT[1],                                      // 22
00100    BIT[5],                                      // 23
00101    BIT[4] + BIT[3] + BIT[1],                    // 24
00102    BIT[3],                                      // 25
00103    BIT[6] + BIT[2] + BIT[1],                    // 26
00104    BIT[5] + BIT[2] + BIT[1],                    // 27
00105    BIT[3],                                      // 28
00106    BIT[2],                                      // 29
00107    BIT[6] + BIT[4] + BIT[1],                    // 30
00108    BIT[3],                                      // 31
00109    BIT[7] + BIT[5] + BIT[3] + BIT[2] + BIT[1]   // 32
00110 };
00111 
00112 // for convenience, define some data structures for bivariate distributions
00113 
00114 struct cartesianCoord {   // cartesian coordinates in 2-D
00115 
00116    float x, y;
00117    cartesianCoord& operator+=( const cartesianCoord& p ) {
00118       x += p.x;
00119       y += p.y;
00120       return *this;
00121    }
00122    cartesianCoord& operator-=( const cartesianCoord& p ) {
00123       x -= p.x;
00124       y -= p.y;
00125       return *this;
00126    }
00127    cartesianCoord& operator*=( float scalar ) {
00128       x *= scalar;
00129       y *= scalar;
00130       return *this;
00131    }
00132    cartesianCoord& operator/=( float scalar ) {
00133       x /= scalar;
00134       y /= scalar;
00135       return *this;
00136    }
00137 };
00138 
00139 struct sphericalCoord {  // spherical coordinates on unit sphere
00140 
00141    float theta, phi;
00142    float x( void ) { return sin( theta ) * cos( phi ); }   // x-coordinate
00143    float y( void ) { return sin( theta ) * sin( phi ); }   // y-coordinate
00144    float z( void ) { return cos( theta ); }                // z-coordinate
00145 };
00146 
00147 class Random {
00148 
00149 // friends list
00150 // overloaded relational operators
00151    
00152    friend bool operator==( const Random& p, const Random& q )
00153    {
00154       bool equal = ( p._seed == q._seed ) && ( p._next == q._next );
00155       for ( int i = 0; i < p._NTAB; i++ )
00156          equal = equal && ( p._table[ i ] == q._table[ i ] );
00157       return equal;
00158    }
00159    
00160    friend bool operator!=( const Random& p, const Random& q )
00161    {
00162       return !( p == q );
00163    }
00164    
00165 // overloaded stream operator
00166    
00167    friend istream& operator>>( istream& is, Random& rv )
00168    {
00169       cout << "Enter a random number seed "
00170            << "(between 1 and " << LONG_MAX - 1 << ", inclusive): " << endl;
00171         
00172       is >> rv._seed;
00173 
00174       assert( rv._seed != 0 && rv._seed != LONG_MAX );
00175    
00176       rv._seedTable();
00177    
00178       return is;
00179    }
00180   
00181 public:
00182 
00183    Random( long seed )   // constructor to set the seed
00184    {
00185      pthread_mutex_init(&random_mutex, NULL);
00186 
00187       assert( seed != 0 && seed != LONG_MAX );
00188       _seed = seed;
00189       _seedTable();
00190       
00191       _seed2 = _seed | 1;  // for tausworthe random bit generation
00192    }
00193    
00194    Random( void )   // default constructor uses process id to set the seed
00195    {
00196      pthread_mutex_init(&random_mutex, NULL);
00197 
00198       _seed = long( getpid() );
00199       _seedTable();
00200       _seed2 = _seed | 1;  // for tausworthe random bit generation
00201    }
00202 
00203   ~Random( void )   // default destructor
00204    {
00205 
00206    }
00207    
00208    Random( const Random& r )   // copy constructor (copies current state)
00209    {
00210      pthread_mutex_init(&random_mutex, NULL);
00211 
00212       _seed  = r._seed;
00213       _seed2 = r._seed2;
00214 
00215       // copy the current state of the shuffle table
00216 
00217       _next = r._next;
00218       for ( int i = 0; i < _NTAB; i++ ) _table[ i ] = r._table[ i ];
00219    }
00220    
00221    Random& operator=( const Random& r )   // overloaded assignment
00222    {
00223       if ( *this == r ) return *this;
00224 
00225       _seed  = r._seed;
00226       _seed2 = r._seed2;
00227    
00228       // copy the current state of the shuffle table
00229 
00230       _next = r._next;
00231       for ( int i = 0; i < _NTAB; i++ ) _table[ i ] = r._table[ i ];
00232 
00233       return *this;
00234    }
00235    
00236 // utility functions
00237 
00238    void reset( long seed )   // reset the seed explicitly
00239    {
00240       assert( seed != 0 && seed != LONG_MAX );
00241       _seed = seed;
00242       _seedTable();
00243       _seed2 = _seed | 1;   // so that all bits cannot be zero
00244    }
00245    
00246    void reset( void )   // reset seed from current process id
00247    {
00248       _seed = long( getpid() );
00249       _seedTable();
00250       _seed2 = _seed | 1;   // so that all bits cannot be zero
00251    }
00252 
00253 // Continuous Distributions
00254 
00255    float arcsine( float xMin = 0., float xMax = 1. )   // Arc Sine
00256    {
00257       float q = sin( M_PI_2 * _u() );
00258       return xMin + ( xMax - xMin ) * q * q;
00259    }
00260    
00261    float beta( float v, float w,                    // Beta
00262                 float xMin = 0., float xMax = 1. )   // (v > 0. and w > 0.)
00263    {
00264       if ( v < w ) return xMax - ( xMax - xMin ) * beta( w, v );
00265       float y1 = gamma( 0., 1., v );
00266       float y2 = gamma( 0., 1., w );
00267       return xMin + ( xMax - xMin ) * y1 / ( y1 + y2 );
00268    }
00269    
00270    float cauchy( float a = 0., float b = 1. )   // Cauchy (or Lorentz)
00271    {
00272       // a is the location parameter and b is the scale parameter
00273       // b is the half width at half maximum (HWHM) and variance doesn't exist
00274    
00275       assert( b > 0. );
00276    
00277       return a + b * tan( M_PI * uniform( -0.5, 0.5 ) );
00278    }
00279 
00280    float chiSquare( int df )   // Chi-Square
00281    {
00282       assert( df >= 1 );
00283    
00284       return gamma( 0., 2., 0.5 * float( df ) );
00285    }
00286 
00287    float cosine( float xMin = 0., float xMax = 1. )   // Cosine
00288    {
00289       assert( xMin < xMax );
00290    
00291       float a = 0.5 * ( xMin + xMax );    // location parameter
00292       float b = ( xMax - xMin ) / M_PI;   // scale parameter
00293    
00294       return a + b * asin( uniform( -1., 1. ) );
00295    }
00296    
00297    float floatLog( float xMin = -1., float xMax = 1. )   // Double Log
00298    {
00299       assert( xMin < xMax );
00300    
00301       float a = 0.5 * ( xMin + xMax );    // location parameter
00302       float b = 0.5 * ( xMax - xMin );    // scale parameter
00303    
00304       if ( bernoulli( 0.5 ) ) return a + b * _u() * _u();
00305       else                    return a - b * _u() * _u();
00306    }
00307 
00308    float erlang( float b, int c )   // Erlang (b > 0. and c >= 1)
00309    {
00310       assert( b > 0. && c >= 1 );
00311   
00312       float prod = 1.;
00313       for ( int i = 0; i < c; i++ ) prod *= _u();
00314       
00315       return -b * log( prod );
00316    }
00317    
00318    float exponential( float a = 0., float c = 1. )   // Exponential
00319    {                                                    // location a, shape c
00320       assert( c > 0.0 );
00321    
00322       return a - c * log( _u() );
00323    }
00324    
00325    float extremeValue( float a = 0., float c = 1. )   // Extreme Value
00326    {                                                     // location a, shape c
00327       assert( c > 0. );
00328    
00329       return a + c * log( -log( _u() ) );
00330    }
00331    
00332    float fRatio( int v, int w )   // F Ratio (v and w >= 1)
00333    {
00334       assert( v >= 1 && w >= 1 );
00335    
00336       return ( chiSquare( v ) / v ) / ( chiSquare( w ) / w );
00337    }
00338    
00339    float gamma( float a, float b, float c )  // Gamma
00340    {                                             // location a, scale b, shape c
00341       assert( b > 0. && c > 0. );
00342    
00343       const float A = 1. / sqrt( 2. * c - 1. );
00344       const float B = c - log( 4. );
00345       const float Q = c + 1. / A;
00346       const float T = 4.5;
00347       const float D = 1. + log( T );
00348       const float C = 1. + c / M_E;
00349    
00350       if ( c < 1. ) {  
00351          while ( true ) {
00352             float p = C * _u();      
00353             if ( p > 1. ) {
00354                float y = -log( ( C - p ) / c );
00355                if ( _u() <= powf( y, c - 1. ) ) return a + b * y;
00356             }
00357             else {
00358                float y = powf( p, 1. / c );
00359                if ( _u() <= exp( -y ) ) return a + b * y;
00360             }
00361          }
00362       }
00363       else if ( c == 1.0 ) return exponential( a, b );
00364       else {
00365          while ( true ) {
00366             float p1 = _u();
00367             float p2 = _u();
00368             float v = A * log( p1 / ( 1. - p1 ) );
00369             float y = c * exp( v );
00370             float z = p1 * p1 * p2;
00371             float w = B + Q * v - y;
00372             if ( w + D - T * z >= 0. || w >= log( z ) ) return a + b * y;
00373          }
00374       }
00375    }
00376    
00377    float laplace( float a = 0., float b = 1. )   // Laplace
00378    {                                                // (or float exponential)
00379       assert( b > 0. );
00380 
00381       // composition method
00382   
00383       if ( bernoulli( 0.5 ) ) return a + b * log( _u() );
00384       else                    return a - b * log( _u() );
00385    }
00386    
00387    float logarithmic( float xMin = 0., float xMax = 1. )   // Logarithmic
00388    {
00389       assert( xMin < xMax );
00390    
00391       float a = xMin;          // location parameter
00392       float b = xMax - xMin;   // scale parameter
00393    
00394       // use convolution formula for product of two IID uniform variates
00395 
00396       return a + b * _u() * _u();
00397    }
00398    
00399    float logistic( float a = 0., float c = 1. )   // Logistic
00400    {
00401       assert( c > 0. );
00402 
00403       return a - c * log( 1. / _u() - 1. );
00404    }
00405    
00406    float lognormal( float a, float mu, float sigma )   // Lognormal
00407    {
00408       return a + exp( normal( mu, sigma ) );
00409    }
00410    
00411    float normal( float mu = 0., float sigma = 1. )   // Normal
00412    {
00413       assert( sigma > 0. );
00414    
00415       static bool f = true;
00416       static float p, p1, p2;
00417       float q;
00418    
00419       if ( f ) {
00420          do {
00421             p1 = uniform( -1., 1. );
00422             p2 = uniform( -1., 1. );
00423             p = p1 * p1 + p2 * p2;
00424          } while ( p >= 1. );
00425          q = p1;
00426       }
00427       else
00428          q = p2;
00429       f = !f;
00430 
00431       return mu + sigma * q * sqrt( -2. * log( p ) / p );
00432    }
00433    
00434    float parabolic( float xMin = 0., float xMax = 1. )   // Parabolic
00435    {  
00436       assert( xMin < xMax );
00437    
00438       float a    = 0.5 * ( xMin + xMax );        // location parameter
00439       float yMax = _parabola( a, xMin, xMax );   // maximum function range
00440    
00441       return userSpecified( _parabola, xMin, xMax, 0., yMax );
00442    }
00443    
00444    float pareto( float c )   // Pareto
00445    {                           // shape c
00446       assert( c > 0. );
00447    
00448       return powf( _u(), -1. / c );
00449    }
00450    
00451    float pearson5( float b, float c )   // Pearson Type 5
00452    {                                       // scale b, shape c
00453       assert( b > 0. && c > 0. );
00454    
00455       return 1. / gamma( 0., 1. / b, c );
00456    }
00457    
00458    float pearson6( float b, float v, float w )   // Pearson Type 6
00459    {                                                 // scale b, shape v & w
00460       assert( b > 0. && v > 0. && w > 0. );
00461    
00462       return gamma( 0., b, v ) / gamma( 0., b, w );
00463    }
00464    
00465    float power( float c )   // Power
00466    {                          // shape c
00467       assert( c > 0. );
00468    
00469       return powf( _u(), 1. / c );
00470    }
00471    
00472    float rayleigh( float a, float b )   // Rayleigh
00473    {                                       // location a, scale b
00474       assert( b > 0. );
00475    
00476       return a + b * sqrt( -log( _u() ) );
00477    }
00478    
00479    float studentT( int df )   // Student's T
00480    {                           // degres of freedom df
00481       assert( df >= 1 );
00482    
00483       return normal() / sqrt( chiSquare( df ) / df );
00484    }
00485    
00486    float triangular( float xMin = 0.,     // Triangular
00487                       float xMax = 1.,     // with default interval [0,1)
00488                       float c    = 0.5 )   // and default mode 0.5
00489    {
00490       assert( xMin < xMax && xMin <= c && c <= xMax );
00491       
00492       float p = _u(), q = 1. - p;
00493       
00494       if ( p <= ( c - xMin ) / ( xMax - xMin ) )
00495          return xMin + sqrt( ( xMax - xMin ) * ( c - xMin ) * p );
00496       else
00497          return xMax - sqrt( ( xMax - xMin ) * ( xMax - c ) * q );
00498    }
00499    
00500    float uniform( float xMin = 0., float xMax = 1. )   // Uniform
00501    {                                                      // on [xMin,xMax)
00502       assert( xMin < xMax );
00503    
00504       return xMin + ( xMax - xMin ) * _u();
00505    }
00506    
00507    float userSpecified(                // User-Specified Distribution
00508         float( *usf )(                 // pointer to user-specified function
00509              float,                    // x
00510              float,                    // xMin
00511              float ),                  // xMax
00512         float xMin, float xMax,       // function domain
00513         float yMin, float yMax )      // function range
00514    {
00515       assert( xMin < xMax && yMin < yMax );
00516    
00517       float x, y, areaMax = ( xMax - xMin ) * ( yMax - yMin ); 
00518 
00519       // acceptance-rejection method
00520    
00521       do {   
00522          x = uniform( 0.0, areaMax ) / ( yMax - yMin ) + xMin;  
00523          y = uniform( yMin, yMax );
00524    
00525       } while ( y > usf( x, xMin, xMax ) );
00526    
00527       return x;
00528    }
00529    
00530    float weibull( float a, float b, float c )   // Weibull
00531    {                                                // location a, scale b,
00532       assert( b > 0. && c > 0. );                   // shape c
00533    
00534       return a + b * powf( -log( _u() ), 1. / c );
00535    }
00536                    
00537 // Discrete Distributions
00538 
00539    bool bernoulli( float p = 0.5 )   // Bernoulli Trial
00540    {
00541       assert( 0. <= p && p <= 1. );
00542    
00543       return _u() < p;
00544    }
00545    
00546    int binomial( int n, float p )   // Binomial
00547    {
00548       assert( n >= 1 && 0. <= p && p <= 1. );
00549    
00550       int sum = 0;
00551       for ( int i = 0; i < n; i++ ) sum += bernoulli( p );
00552       return sum;
00553    }
00554    
00555    int geometric( float p )   // Geometric
00556    {
00557       assert( 0. < p && p < 1. );
00558 
00559       return int( log( _u() ) / log( 1. - p ) );
00560    }
00561    
00562    int hypergeometric( int n, int N, int K )            // Hypergeometric
00563    {                                                    // trials n, size N,
00564       assert( 0 <= n && n <= N && N >= 1 && K >= 0 );   // successes K
00565       
00566       int count = 0;
00567       for ( int i = 0; i < n; i++, N-- ) {
00568    
00569          float p = float( K ) / float( N );
00570          if ( bernoulli( p ) ) { count++; K--; }
00571       }
00572       return count;
00573    }
00574    
00575    void multinomial( int    n,            // Multinomial
00576                      float p[],          // trials n, probability vector p,
00577                      int    count[],      // success vector count,
00578                      int    m )           // number of disjoint events m
00579    {
00580       assert( m >= 2 );   // at least 2 events
00581       float sum = 0.;
00582       for ( int bin = 0; bin < m; bin++ ) sum += p[ bin ];    // probabilities
00583       assert( sum == 1. );                                    // must sum to 1
00584       
00585       for ( int bin = 0; bin < m; bin++ ) count[ bin ] = 0;   // initialize
00586    
00587       // generate n uniform variates in the interval [0,1) and bin the results
00588    
00589       for ( int i = 0; i < n; i++ ) {
00590 
00591          float lower = 0., upper = 0., u = _u();
00592 
00593          for ( int bin = 0; bin < m; bin++ ) {
00594 
00595          // locate subinterval, which is of length p[ bin ],
00596          // that contains the variate and increment the corresponding counter
00597       
00598             lower = upper;
00599             upper += p[ bin ];
00600             if ( lower <= u && u < upper ) { count[ bin ]++; break; }
00601          }
00602       }
00603    }
00604    
00605    int negativeBinomial( int s, float p )   // Negative Binomial
00606    {                                         // successes s, probability p
00607       assert( s >= 1 && 0. < p && p < 1. );
00608    
00609       int sum = 0;
00610       for ( int i = 0; i < s; i++ ) sum += geometric( p );
00611       return sum;
00612    }
00613    
00614    int pascal( int s, float p )              // Pascal
00615    {                                          // successes s, probability p
00616       return negativeBinomial( s, p ) + s;
00617    }
00618    
00619    int poisson( float mu )   // Poisson
00620    {
00621       assert ( mu > 0. );
00622    
00623       float a = exp( -mu );
00624       float b = 1.;
00625    
00626       int i;
00627       for ( i = 0; b >= a; i++ ) b *= _u();   
00628       return i - 1;
00629    }
00630    
00631    int uniformDiscrete( int i, int j )   // Uniform Discrete
00632    {                                     // inclusive i to j
00633       assert( i < j );
00634 
00635       return i + int( ( j - i + 1 ) * _u() );
00636    }
00637 
00638 // Empirical and Data-Driven Distributions
00639 
00640    float empirical( void )   // Empirical Continuous
00641    {
00642       static vector< float > x, cdf;
00643       static int              n;
00644       static bool             init = false;
00645    
00646       if ( !init ) {
00647          ifstream in( "empiricalDistribution" );
00648          if ( !in ) {
00649             cerr << "Cannot open \"empiricalDistribution\" input file" << endl;
00650             exit( 1 );
00651          }
00652          float value, prob;
00653          while ( in >> value >> prob ) {   // read in empirical distribution
00654             x.push_back( value );
00655             cdf.push_back( prob );
00656          }
00657          n = x.size();
00658          init = true;
00659       
00660          // check that this is indeed a cumulative distribution
00661 
00662          assert( 0. == cdf[ 0 ] && cdf[ n - 1 ] == 1. );
00663          for ( int i = 1; i < n; i++ ) assert( cdf[ i - 1 ] < cdf[ i ] );
00664       }
00665 
00666       float p = _u();
00667       for ( int i = 0; i < n - 1; i++ )
00668          if ( cdf[ i ] <= p && p < cdf[ i + 1 ] )
00669             return x[ i ] + ( x[ i + 1 ] - x[ i ] ) *
00670                             ( p - cdf[ i ] ) / ( cdf[ i + 1 ] - cdf[ i ] );
00671       return x[ n - 1 ];
00672    }
00673    
00674    int empiricalDiscrete( void )   // Empirical Discrete
00675    {
00676       static vector< int >    k;
00677       static vector< float > f[ 2 ];   // f[ 0 ] is pdf and f[ 1 ] is cdf
00678       static float           max;
00679       static int              n;
00680       static bool             init = false;
00681    
00682       if ( !init ) {
00683          ifstream in ( "empiricalDiscrete" );
00684          if ( !in ) {
00685             cerr << "Cannot open \"empiricalDiscrete\" input file" << endl;
00686             exit( 1 );
00687          }
00688          int value;
00689          float freq;
00690          while ( in >> value >> freq ) {   // read in empirical data
00691             k.push_back( value );
00692             f[ 0 ].push_back( freq );
00693          }
00694          n = k.size();
00695          init = true;
00696    
00697          // form the cumulative distribution
00698 
00699          f[ 1 ].push_back( f[ 0 ][ 0 ] );
00700          for ( int i = 1; i < n; i++ )
00701             f[ 1 ].push_back( f[ 1 ][ i - 1 ] + f[ 0 ][ i ] );
00702 
00703          // check that the integer points are in ascending order
00704 
00705          for ( int i = 1; i < n; i++ ) assert( k[ i - 1 ] < k[ i ] );
00706       
00707          max = f[ 1 ][ n - 1 ];
00708       }
00709    
00710       // select a uniform variate between 0 and the max value of the cdf
00711 
00712       float p = uniform( 0., max );
00713 
00714       // locate and return the corresponding index
00715 
00716       for ( int i = 0; i < n; i++ ) if ( p <= f[ 1 ][ i ] ) return k[ i ];
00717       return k[ n - 1 ];
00718    }
00719    
00720    float sample( bool replace = true )  // Sample w or w/o replacement from a
00721    {                                     // distribution of 1-D data in a file
00722       static vector< float > v;         // vector for sampling with replacement
00723       static bool init = false;          // flag that file has been read in
00724       static int n;                      // number of data elements in the file
00725       static int index = 0;              // subscript in the sequential order
00726    
00727       if ( !init ) {
00728          ifstream in( "sampleData" );
00729          if ( !in ) {
00730             cerr << "Cannot open \"sampleData\" file" << endl;
00731             exit( 1 );
00732          }
00733          float d;
00734          while ( in >> d ) v.push_back( d );
00735          in.close();
00736          n = v.size();
00737          init = true;
00738          if ( replace == false ) {   // sample without replacement
00739          
00740             // shuffle contents of v once and for all
00741             // Ref: Knuth, D. E., The Art of Computer Programming, Vol. 2:
00742             //      Seminumerical Algorithms. London: Addison-Wesley, 1969.
00743 
00744             for ( int i = n - 1; i > 0; i-- ) {
00745                int j = int( ( i + 1 ) * _u() );
00746                swap( v[ i ], v[ j ] );
00747             }
00748          }
00749       }
00750 
00751       // return a random sample
00752 
00753       if ( replace )                                // sample w/ replacement
00754          return v[ uniformDiscrete( 0, n - 1 ) ];
00755       else {                                        // sample w/o replacement
00756          assert( index < n );                       // retrieve elements
00757          return v[ index++ ];                       // in sequential order
00758       }
00759    }
00760    
00761    void sample( float x[], int ndim )   // Sample from a given distribution
00762    {                                     // of multi-dimensional data
00763       static const int N_DIM = 6;
00764       assert( ndim <= N_DIM );
00765    
00766       static vector< float > v[ N_DIM ];
00767       static bool init = false;
00768       static int n;
00769    
00770       if ( !init )  {
00771          ifstream in( "sampleData" );
00772          if ( !in ) {
00773             cerr << "Cannot open \"sampleData\" file" << endl;
00774             exit( 1 );
00775          }
00776          float d;
00777          while ( !in.eof() ) {
00778             for ( int i = 0; i < ndim; i++ ) {
00779                in >> d;
00780                v[ i ].push_back( d );
00781             }
00782          }
00783          in.close();
00784          n = v[ 0 ].size();
00785          init = true;
00786       }
00787       int index = uniformDiscrete( 0, n - 1 );
00788       for ( int i = 0; i < ndim; i++ ) x[ i ] = v[ i ][ index ];
00789    }
00790 
00791    // comparison functor for determining the neighborhood of a data point
00792 
00793    struct dSquared :
00794    public binary_function< cartesianCoord, cartesianCoord, bool > {
00795       bool operator()( cartesianCoord p, cartesianCoord q ) {
00796          return p.x * p.x + p.y * p.y < q.x * q.x + q.y * q.y;
00797       }
00798    };
00799 
00800    cartesianCoord stochasticInterpolation( void )   // Stochastic Interpolation
00801 
00802    // Refs: Taylor, M. S. and J. R. Thompson, Computational Statistics & Data 
00803    //       Analysis, Vol. 4, pp. 93-101, 1986; Thompson, J. R., Empirical Model
00804    //       Building, pp. 108-114, Wiley, 1989; Bodt, B. A. and M. S. Taylor,
00805    //       A Data Based Random Number Generator for A Multivariate Distribution
00806    //       - A User's Manual, ARBRL-TR-02439, BRL, APG, MD, Nov. 1982.
00807    {
00808       static vector< cartesianCoord > data;
00809       static cartesianCoord           min, max;
00810       static int                      m;
00811       static float                   lower, upper;
00812       static bool                     init = false;
00813 
00814       if ( !init ) {
00815          ifstream in( "stochasticData" );
00816          if ( !in ) {
00817             cerr << "Cannot open \"stochasticData\" input file" << endl;
00818             exit( 1 );
00819          }
00820    
00821          // read in the data and set min and max values
00822    
00823          min.x = min.y = FLT_MAX;
00824          max.x = max.y = FLT_MIN;
00825          cartesianCoord p;
00826          while ( in >> p.x >> p.y ) {
00827 
00828             min.x = ( p.x < min.x ? p.x : min.x );
00829             min.y = ( p.y < min.y ? p.y : min.y );
00830             max.x = ( p.x > max.x ? p.x : max.x );
00831             max.y = ( p.y > max.y ? p.y : max.y );
00832       
00833             data.push_back( p );
00834          }
00835          in.close();
00836          init = true;
00837    
00838          // scale the data so that each dimension will have equal weight
00839 
00840          for ( unsigned i = 0; i < data.size(); i++ ) {
00841        
00842             data[ i ].x = ( data[ i ].x - min.x ) / ( max.x - min.x );
00843             data[ i ].y = ( data[ i ].y - min.y ) / ( max.y - min.y );
00844          }
00845    
00846          // set m, the number of points in a neighborhood of a given point
00847    
00848          m = data.size() / 20;       // 5% of all the data points
00849          if ( m < 5  ) m = 5;        // but no less than 5
00850          if ( m > 20 ) m = 20;       // and no more than 20
00851    
00852          lower = ( 1. - sqrt( 3. * ( float( m ) - 1. ) ) ) / float( m );
00853          upper = ( 1. + sqrt( 3. * ( float( m ) - 1. ) ) ) / float( m );
00854       }
00855    
00856       // uniform random selection of a data point (with replacement)
00857       
00858       cartesianCoord origin = data[ uniformDiscrete( 0, data.size() - 1 ) ];
00859 
00860       // make this point the origin of the coordinate system
00861 
00862       for ( unsigned n = 0; n < data.size(); n++ ) data[ n ] -= origin;
00863       
00864       // sort the data with respect to its distance (squared) from this origin
00865       
00866       sort( data.begin(), data.end(), dSquared() );
00867       
00868       // find the mean value of the data in the neighborhood about this point
00869       
00870       cartesianCoord mean;
00871       mean.x = mean.y = 0.;
00872       for ( int n = 0; n < m; n++ ) mean += data[ n ];
00873       mean /= float( m );
00874    
00875       // select a random linear combination of the points in this neighborhood
00876 
00877       cartesianCoord p;
00878       p.x = p.y = 0.;
00879       for ( int n = 0; n < m; n++ ) {
00880          
00881          float rn;
00882          if ( m == 1 ) rn = 1.;
00883          else          rn = uniform( lower, upper );
00884          p.x += rn * ( data[ n ].x - mean.x );
00885          p.y += rn * ( data[ n ].y - mean.y );
00886       }
00887       
00888       // restore the data to its original form
00889 
00890       for ( unsigned n = 0; n < data.size(); n++ ) data[ n ] += origin;
00891       
00892       // use mean and original point to translate the randomly-chosen point
00893 
00894       p += mean;
00895       p += origin;
00896 
00897       // scale randomly-chosen point to the dimensions of the original data
00898       
00899       p.x = p.x * ( max.x - min.x ) + min.x;
00900       p.y = p.y * ( max.y - min.y ) + min.y;
00901 
00902       return p;
00903    }
00904 
00905 // Multivariate Distributions
00906 
00907    cartesianCoord bivariateNormal( float muX    = 0.,   // Bivariate Gaussian
00908                                    float sigmaX = 1.,
00909                                    float muY    = 0., 
00910                                    float sigmaY = 1. )
00911    {
00912       assert( sigmaX > 0. && sigmaY > 0. );
00913    
00914       cartesianCoord p;
00915       p.x = normal( muX, sigmaX );
00916       p.y = normal( muY, sigmaY );
00917       return p;
00918    }
00919    
00920    cartesianCoord bivariateUniform( float xMin = -1.,    // Bivariate Uniform
00921                                     float xMax =  1.,
00922                                     float yMin = -1.,
00923                                     float yMax =  1. )
00924    {
00925       assert( xMin < xMax && yMin < yMax );
00926    
00927       float x0 = 0.5 * ( xMin + xMax );
00928       float y0 = 0.5 * ( yMin + yMax );
00929       float a  = 0.5 * ( xMax - xMin );
00930       float b  = 0.5 * ( yMax - yMin );
00931       float x, y;
00932    
00933       do {
00934          x = uniform( -1., 1. );
00935          y = uniform( -1., 1. );
00936       
00937       } while( x * x + y * y > 1. );
00938       
00939       cartesianCoord p;
00940       p.x = x0 + a * x;
00941       p.y = y0 + b * y;
00942       return p;
00943    }
00944    
00945    cartesianCoord corrNormal( float r,              // Correlated Normal
00946                               float muX    = 0.,
00947                               float sigmaX = 1.,
00948                               float muY    = 0.,
00949                               float sigmaY = 1. )
00950    {
00951       assert( -1. <= r && r <= 1. &&          // bounds on correlation coeff
00952               sigmaX > 0. && sigmaY > 0. );   // positive std dev
00953    
00954       float x = normal();
00955       float y = normal();
00956    
00957       y = r * x + sqrt( 1. - r * r ) * y;     // correlate the variables
00958    
00959       cartesianCoord p;
00960       p.x = muX + sigmaX * x;                 // translate and scale
00961       p.y = muY + sigmaY * y;
00962       return p;
00963    }
00964    
00965    cartesianCoord corrUniform( float r,        // Correlated Uniform
00966                                float xMin = 0.,
00967                                float xMax = 1.,
00968                                float yMin = 0.,
00969                                float yMax = 1. )
00970    {
00971       assert( -1. <= r && r <= 1. &&          // bounds on correlation coeff
00972               xMin < xMax && yMin < yMax );   // bounds on domain
00973 
00974       float x0 = 0.5 * ( xMin + xMax );
00975       float y0 = 0.5 * ( yMin + yMax );
00976       float a  = 0.5 * ( xMax - xMin );
00977       float b  = 0.5 * ( yMax - yMin );
00978       float x, y;
00979    
00980       do {
00981          x = uniform( -1., 1. );
00982          y = uniform( -1., 1. );
00983       
00984       } while ( x * x + y * y > 1. );
00985    
00986       y = r * x + sqrt( 1. - r * r ) * y;   // correlate the variables
00987    
00988       cartesianCoord p;
00989       p.x = x0 + a * x;                     // translate and scale
00990       p.y = y0 + b * y;
00991       return p;
00992    }
00993    
00994    sphericalCoord spherical( float thMin = 0.,     // Uniform Spherical
00995                              float thMax = M_PI,
00996                              float phMin = 0.,
00997                              float phMax = 2. * M_PI )
00998    {
00999       assert( 0. <= thMin && thMin < thMax && thMax <= M_PI &&
01000               0. <= phMin && phMin < phMax && phMax <= 2. * M_PI );
01001    
01002       sphericalCoord p;                       // azimuth
01003       p.theta = acos( uniform( cos( thMax ), cos( thMin ) ) );   // polar angle
01004       p.phi   = uniform( phMin, phMax );                         // azimuth
01005       return p;
01006    }
01007    
01008    void sphericalND( float x[], int n )   // Uniform over the surface of
01009                                            // an n-dimensional unit sphere
01010    {
01011       // generate a point inside the unit n-sphere by normal polar method
01012 
01013       float r2 = 0.;
01014       for ( int i = 0; i < n; i++ ) {
01015          x[ i ] = normal();
01016          r2 += x[ i ] * x[ i ];
01017       }
01018    
01019       // project the point onto the surface of the unit n-sphere by scaling
01020 
01021       const float A = 1. / sqrt( r2 );
01022       for ( int i = 0; i < n; i++ ) x[ i ] *= A;
01023    }
01024 
01025 // Number Theoretic Distributions
01026 
01027    float avoidance( void )   // Maximal Avoidance (1-D)
01028    {                          // overloaded for convenience
01029       float x[ 1 ];
01030       avoidance( x, 1 );
01031       return x[ 0 ];
01032    }
01033    
01034    void avoidance( float x[], unsigned ndim )   // Maximal Avoidance (N-D)
01035    {
01036       static const unsigned MAXBIT = 30;
01037       static const unsigned MAXDIM = 6;
01038    
01039       assert( ndim <= MAXDIM );
01040 
01041       static unsigned long ix[ MAXDIM + 1 ] = { 0 };
01042       static unsigned long *u[ MAXBIT + 1 ];
01043       static unsigned long mdeg[ MAXDIM + 1 ] = { // degree of
01044          0,                                       // primitive polynomial
01045          1, 2, 3, 3, 4, 4
01046       };
01047       static unsigned long p[ MAXDIM + 1 ] = {   // decimal encoded
01048          0,                                      // interior bits
01049          0, 1, 1, 2, 1, 4
01050       };
01051       static unsigned long v[ MAXDIM * MAXBIT + 1 ] = {
01052           0,
01053           1,  1, 1,  1,  1,  1,
01054           3,  1, 3,  3,  1,  1,
01055           5,  7, 7,  3,  3,  5,
01056          15, 11, 5, 15, 13,  9
01057       };
01058 
01059       static float fac;
01060       static int in = -1;
01061       unsigned j, k;
01062       unsigned long i, m, pp;
01063       
01064       if ( in == -1 ) {
01065          in = 0;
01066          fac = 1. / ( 1L << MAXBIT );
01067          for ( j = 1, k = 0; j <= MAXBIT; j++, k += MAXDIM ) u[ j ] = &v[ k ];
01068          for ( k = 1; k <= MAXDIM; k++ ) {
01069             for ( j = 1; j <= mdeg[ k ]; j++ ) u[ j ][ k ] <<= ( MAXBIT - j );
01070             for ( j = mdeg[ k ] + 1; j <= MAXBIT; j++ ) {
01071                pp = p[ k ];
01072                i = u[ j - mdeg[ k ] ][ k ];
01073                i ^= ( i >> mdeg[ k ] );
01074                for ( int n = mdeg[ k ] - 1; n >= 1; n-- ) {
01075                   if ( pp & 1 ) i ^= u[ j - n ][ k ];
01076                   pp >>= 1;
01077                }
01078                u[ j ][ k ] = i;
01079             }
01080          }
01081       }
01082       m = in++;
01083       for ( j = 0; j < MAXBIT; j++, m >>= 1 ) if ( !( m & 1 ) ) break;
01084       if ( j >= MAXBIT ) exit( 1 );
01085       m = j * MAXDIM;
01086       for ( k = 0; k < ndim; k++ ) {
01087          ix[ k + 1 ] ^= v[ m + k + 1 ];
01088          x[ k ] = ix[ k + 1 ] * fac;
01089       }
01090    }
01091    
01092    bool tausworthe( unsigned n )   // Tausworthe random bit generator
01093    {                               // returns a single random bit
01094       assert( 1 <= n && n <= 32 );
01095 
01096       if ( _seed2 & BIT[ n ] ) {
01097          _seed2 = ( ( _seed2 ^ MASK[ n ] ) << 1 ) | BIT[ 1 ];
01098          return true;
01099       }
01100       else {
01101          _seed2 <<= 1;
01102          return false;
01103       }
01104    }
01105    
01106    void tausworthe( bool*    bitvec,   // Tausworthe random bit generator
01107                     unsigned n )       // returns a bit vector of length n
01108                     
01109    // It is guaranteed to cycle through all possible combinations of n bits
01110    // (except all zeros) before repeating, i.e., cycle has maximal length 2^n-1.
01111    // Ref: Press, W. H., B. P. Flannery, S. A. Teukolsky and W. T. Vetterling,
01112    //      Numerical Recipes in C, Cambridge Univ. Press, Cambridge, 1988.
01113    {
01114       assert( 1 <= n && n <= 32 );   // length of bit vector
01115 
01116       if ( _seed2 & BIT[ n ] )
01117          _seed2 = ( ( _seed2 ^ MASK[ n ] ) << 1 ) | BIT[ 1 ];
01118       else
01119          _seed2 <<= 1;
01120       for ( unsigned i = 0; i < n; i++ ) bitvec[ i ] = _seed2 & ( BIT[ n ] >> i );
01121    }
01122                                  
01123 private:
01124 
01125    static const long   _M    = 0x7fffffff; // 2147483647 (Mersenne prime 2^31-1)
01126    static const long   A_    = 0x10ff5;    // 69621
01127    static const long   Q_    = 0x787d;     // 30845
01128    static const long   R_    = 0x5d5e;     // 23902
01129    static const float _F    = 1. / _M;
01130    static const short  _NTAB = 32;         // arbitrary length of shuffle table
01131    static const long   _DIV  = 1+(_M-1)/_NTAB;
01132    long         _table[ _NTAB ];          // shuffle table of seeds
01133    long         _next;                    // seed to be used as index into table
01134    long         _seed;                    // current random number seed
01135    unsigned     _seed2;                   // seed for tausworthe random bit
01136    pthread_mutex_t random_mutex;
01137 
01138    void _seedTable( void )                          // seeds the shuffle table
01139    {
01140      pthread_mutex_lock(&random_mutex);
01141      
01142      for ( int i = _NTAB + 7; i >= 0; i-- ) {      // first perform 8 warm-ups
01143        
01144        long k = _seed / Q_;                       // seed = ( A * seed ) % M
01145        _seed = A_ * ( _seed - k * Q_ ) - k * R_;  // without overflow by
01146        if ( _seed < 0 ) _seed += _M;              // Schrage's method
01147        
01148        if ( i < _NTAB ) _table[ i ] = _seed;      // store seeds into table
01149      }
01150      _next = _table[ 0 ];                          // used as index next time
01151      
01152      pthread_mutex_unlock(&random_mutex);
01153    }
01154    
01155    float _u( void )                                // uniform rng
01156    { 
01157      pthread_mutex_lock(&random_mutex);
01158      
01159      long k = _seed / Q_;                          // seed = ( A*seed ) % M
01160      _seed = A_ * ( _seed - k * Q_ ) - k * R_;     // without overflow by
01161      if ( _seed < 0 ) _seed += _M;                 // Schrage's method
01162      
01163      int index = _next / _DIV;                     // Bays-Durham shuffle
01164      _next = _table[ index ];                      // seed used for next time
01165      _table[ index ] = _seed;                      // replace with new seed
01166      float retval = _next * _F;                   // scale value within [0,1)
01167      
01168      pthread_mutex_unlock(&random_mutex);
01169      
01170      return retval;                           
01171    }
01172    
01173    static float _parabola( float x, float xMin, float xMax )  // parabola
01174    {
01175       if ( x < xMin || x > xMax ) return 0.0;
01176    
01177       float a    = 0.5 * ( xMin + xMax );   // location parameter
01178       float b    = 0.5 * ( xMax - xMin );   // scale parameter
01179       float yMax = 0.75 / b;
01180       
01181       return yMax * ( 1. - ( x - a ) * ( x - a ) / ( b * b ) );
01182    }
01183 };
01184 
01185 #endif


rl_common
Author(s):
autogenerated on Thu Jun 6 2019 22:00:08