Nemo  2.3.46
Uniform.h
Go to the documentation of this file.
1 
28 #ifndef __UNIFORM_H
29 #define __UNIFORM_H
30 #include <cmath>
31 #include <limits.h>
32 #include <assert.h>
33 #include "output.h"
34 #include <iostream>
35 
36 #ifdef HAS_SPRNG
37  #define SIMPLE_SPRNG
38  #include <sprng.h>
39 #endif
40 
41 #ifdef HAS_GSL
42  #include <gsl/gsl_rng.h>
43  #include <gsl/gsl_randist.h>
44  #include <gsl/gsl_blas.h>
45  #include <gsl/gsl_math.h>
46  #include <gsl/gsl_vector.h>
47  #include <gsl/gsl_matrix.h>
48  #include <gsl/gsl_permutation.h>
49 #endif
50 
52 class RAND {
53 private:
54 
55  RAND();
56 
57 public:
58 
59 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
60 
61  static gsl_rng * r;
62 
63  static void init_gsl(const gsl_rng_type* T, unsigned long seed)
64  {
65  r = gsl_rng_alloc (T);
66  gsl_rng_set(r,seed);
67  }
68 
69  static void free()
70  {
71  gsl_rng_free(r);
72  }
73 
74 #elif !defined(HAS_GSL) && !defined(HAS_SPRNG)
75 
76  static long Seed1,Seed2;
77 
78 #endif
79 
80  static void init (unsigned long seed) {
81 
82 #if defined(HAS_SPRNG)
83 
84  init_sprng( SPRNG_LFG, seed, SPRNG_DEFAULT );
85 
86 #elif defined(HAS_GSL)
87 
88  init_gsl( gsl_rng_mt19937, seed );
89 
90 #else
91 
92  Seed1 = seed;
93 
94 #endif
95  }
101  static inline double Uniform () {
102 
103 #ifdef HAS_SPRNG
104  return sprng();
105 #elif defined(HAS_GSL)
106  return gsl_rng_uniform(r);
107 #else
108  register long z, w;
109 
110  do{
111  w = Seed1 / 53668;
112 
113  Seed1 = 40014 * (Seed1 - w * 53668) - w * 12211;
114 
115  if (Seed1 < 0) Seed1 += 2147483563;
116 
117  w = (Seed2 / 52774);
118 
119  Seed2 = 40692 * (Seed2 - w * 52774) - w * 3791;
120 
121  if (Seed2 < 0) Seed2 += 2147483399;
122 
123  z = Seed1 - Seed2;
124 
125  if (z < 1) z += 2147483562;
126 
127  }while (!((z * 4.656613e-10) < 1.0));
128 
129  return (z * 4.656613e-10);
130 #endif
131  }
133  static inline unsigned int Uniform (unsigned int max)
134  {
135  return (unsigned int)(Uniform() * max);
136  }
137 
139  static inline bool RandBool() {
140  //generate a random int (or long)
141 #ifdef HAS_SPRNG
142  static int intrand = isprng();
143 #elif defined(HAS_GSL)
144  static unsigned long int intrand = gsl_rng_get(r);
145 #else
146  static int intrand = (int)(Uniform()*(double)std::numeric_limits<int>::max());
147 #endif
148  //read up to the first 16 bits
149  static unsigned int num = 0;
150  //redraw an number after that limit
151  if ( ++num > 16 ) {
152  num = 0;
153 
154 #ifdef HAS_SPRNG
155  intrand = isprng();
156 #elif defined(HAS_GSL)
157  intrand = gsl_rng_get(r);
158 #else
159  intrand = (int)(Uniform()*(double)std::numeric_limits<int>::max());
160 #endif
161 
162  }
163  return (intrand & (1 << num) );
164 
165  }
166 
168  static inline unsigned long RandULong() {
169 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
170  return gsl_rng_get(r);
171 #else
172  unsigned long rnd, limit = 0x10000000; //=2^28 this gives ~7% redraws
173  do{
174  rnd = (unsigned long)(Uniform()*ULONG_MAX);
175  }while(rnd < limit);
176 
177  return rnd;
178 #endif
179  }
181  static inline double gammln (double xx) {
182  double x,y,tmp,ser=1.000000000190015;
183  static double cof[6]={76.18009172947146,-86.50532032941677,
184  24.01409824083091,-1.231739572450155,
185  0.1208650973866179e-2,-0.5395239384953e-5};
186  int j;
187  y=x=xx;
188  tmp=x+5.5;
189  tmp -= (x+0.5)*log(tmp);
190  for (j = 0; j < 6; ++j) ser += cof[j]/++y;
191 
192  return -tmp+log(2.5066282746310005*ser/x);
193 
194  }
196  static inline double Poisson (double mean) {
197 
198 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
199  return gsl_ran_poisson(r, mean);
200 #else
201  static double sq,alxm,g,oldm=(-1.0);
202  double em,t,y;
203 
204  if (mean < 12.0)
205  {
206  if (mean != oldm){
207  oldm=mean;
208  g=exp(-mean);
209  }
210  em = -1;
211  t=1.0;
212  do {
213  ++em;
214  t *= Uniform();
215  } while (t > g);
216  }else
217  {
218  if (mean != oldm)
219  {
220  oldm=mean;
221  sq=sqrt(2.0*mean);
222  alxm=log(mean);
223  g=mean*alxm-gammln(mean+1.0);
224  }
225  do {
226  do {
227  y=tan(M_PI*Uniform());
228  em=sq*y+mean;
229  } while (em < 0.0);
230  em=floor(em);
231  t=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);
232  } while (Uniform() > t);
233  }
234  return em;
235 #endif
236  }
237 
238  static inline double Gaussian(double sigma)
239  {
240 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
241 
242  return gsl_ran_gaussian_ratio_method (r, sigma);
243 
244 #else
245 
246 // double x, y, r2;
247 //
248 // do
249 // {
250 // /* choose x,y in uniform square (-1,-1) to (+1,+1) */
251 //
252 // x = -1 + 2 * Uniform ( );
253 // y = -1 + 2 * Uniform ( );
254 //
255 // /* see if it is in the unit circle */
256 // r2 = x * x + y * y;
257 // }
258 // while (r2 > 1.0 || r2 == 0);
259 //
260 // /* Box-Muller transform */
261 // return sigma * y * sqrt (-2.0 * log (r2) / r2);
262 
263  /*trying the ratio method from GSL*/
264  /* see code in gsl/randist/gauss.c */
265  double u, v, x, y, Q;
266  const double s = 0.449871; /* Constants from Leva */
267  const double t = -0.386595;
268  const double a = 0.19600;
269  const double b = 0.25472;
270  const double r1 = 0.27597;
271  const double r2 = 0.27846;
272 
273  do /* This loop is executed 1.369 times on average */
274  {
275  /* Generate a point P = (u, v) uniform in a rectangle enclosing
276  the K+M region v^2 <= - 4 u^2 log(u). */
277 
278  /* u in (0, 1] to avoid singularity at u = 0 */
279  u = 1 - RAND::Uniform();
280 
281  /* v is in the asymmetric interval [-0.5, 0.5). However v = -0.5
282  is rejected in the last part of the while clause. The
283  resulting normal deviate is strictly symmetric about 0
284  (provided that v is symmetric once v = -0.5 is excluded). */
285  v = RAND::Uniform() - 0.5;
286 
287  /* Constant 1.7156 > sqrt(8/e) (for accuracy); but not by too
288  much (for efficiency). */
289  v *= 1.7156;
290 
291  /* Compute Leva's quadratic form Q */
292  x = u - s;
293  y = fabs (v) - t;
294  Q = x * x + y * (a * y - b * x);
295 
296  /* Accept P if Q < r1 (Leva) */
297  /* Reject P if Q > r2 (Leva) */
298  /* Accept if v^2 <= -4 u^2 log(u) (K+M) */
299  /* This final test is executed 0.012 times on average. */
300  }
301  while (Q >= r1 && (Q > r2 || v * v > -4 * u * u * log (u)));
302 
303  return sigma * (v / u);
304 #endif
305  }
306 
307  static inline void BivariateGaussian(double sigma1,double sigma2, double rho, double *out1,double *out2) {
308 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
309  gsl_ran_bivariate_gaussian(r,sigma1,sigma2,rho,out1,out2);
310 #else
311  //gsl code:
312  double u, v, r2, scale;
313  //double *x = out, *y = (++out);
314  do
315  {
316  /* choose x,y in uniform square (-1,-1) to (+1,+1) */
317 
318  u = -1 + 2 * Uniform ();
319  v = -1 + 2 * Uniform ();
320 
321  /* see if it is in the unit circle */
322  r2 = u * u + v * v;
323  }
324  while (r2 > 1.0 || r2 == 0);
325 
326  scale = sqrt (-2.0 * log (r2) / r2);
327 
328  *out1 = sigma1 * u * scale;
329  *out2 = sigma2 * (rho * u + sqrt(1 - rho*rho) * v) * scale;
330 
331 #endif
332  }
333 
334  static inline double LogNormal (double zeta, double sigma) {
335 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
336  return gsl_ran_lognormal(r,zeta,sigma);
337 #else
338  //this is the GSL code:
339  double u, v, r2, normal, z;
340 
341  do
342  {
343  /* choose x,y in uniform square (-1,-1) to (+1,+1) */
344 
345  u = -1 + 2 * Uniform();
346  v = -1 + 2 * Uniform();
347 
348  /* see if it is in the unit circle */
349  r2 = u * u + v * v;
350  }
351  while (r2 > 1.0 || r2 == 0);
352 
353  normal = u * sqrt (-2.0 * log (r2) / r2);
354 
355  z = exp (sigma * normal + zeta);
356 
357  return z;
358 
359 #endif
360  }
361 
362  static inline double Gamma (double a, double b) {
363 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
364  return gsl_ran_gamma(r, a, b);
365 #else
366  /*pasting GSL code in here*/
367 
368  /* assume a > 0 */
369 
370  if (a < 1)
371  {
372  double u = Uniform();//might give singularity at 0....
373  return Gamma(1.0 + a, b) * pow (u, 1.0 / a);
374  }
375 
376  {
377  double x, v, u;
378  double d = a - 1.0 / 3.0;
379  double c = (1.0 / 3.0) / sqrt (d);
380 
381  while (1)
382  {
383  do
384  {
385  x = Gaussian(1.0); //should use the Ziggurat method?
386  v = 1.0 + c * x;
387  }
388  while (v <= 0);
389 
390  v = v * v * v;
391  u = Uniform();//might give singularity at 0....
392 
393  if (u < 1 - 0.0331 * x * x * x * x)
394  break;
395 
396  if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
397  break;
398  }
399 
400  return b * d * v;
401  }
402 
403 #endif
404  }
405 
406  static inline double Bernoulli (double p)
407  {
408 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
409  return gsl_ran_bernoulli(r, p);
410 #else
411  double u = RAND::Uniform() ;
412 
413  if (u < p)
414  {
415  return 1 ;
416  }
417  else
418  {
419  return 0 ;
420  }
421 #endif
422  }
423 
424  static inline double Exponential (double mu) {
425  return -mu * log(RAND::Uniform());
426  }
427 
428 #ifdef HAS_GSL
429  static inline void MultivariateGaussian(gsl_vector *eval, gsl_matrix *evec,
430  gsl_vector *workspace, gsl_vector *px)
431  {
432  size_t i;
433 
434  for (i=0; i<eval->size; i++)
435  gsl_vector_set (workspace, i,
436  Gaussian(gsl_vector_get (eval, i)));
437 
438  gsl_blas_dgemv (CblasNoTrans, 1.0,
439  evec, workspace, 0.0, px); /* px = evec * px */
440 
441  }
442 
443  static inline void MultivariateGaussianCholesky(gsl_vector *sigma, gsl_matrix *M, gsl_vector *px)
444  {//generate n i.i.d. random normal variates
445  for (size_t i = 0; i < sigma->size; i++)
446  gsl_vector_set (px, i, Gaussian(gsl_vector_get(sigma, i)));
447 
448  gsl_blas_dtrmv (CblasLower, CblasNoTrans, CblasNonUnit, M, px);
449  }
450 
451  static inline long double factoriel(unsigned int x) {
452  long double f = 1;
453  for(unsigned int i=1;i<=x;++i)
454  f *= i;
455  return f;
456  }
457 
458 // static inline double Binomial (double p,unsigned int k,long double N,
459 // const unsigned int n) {
460 // register long double bincoeff = N/(factoriel(k)*factoriel(n - k));
461 // return (bincoeff*pow(p,(int)k)*pow(1-p,(int)(n - k)));
462 // }
463 
464  static inline double Binomial(double p, unsigned int n)
465  {
466  /*implements Knuth method*/
467  unsigned int i, a, b, k = 0;
468 
469  while (n > 10) /* This parameter is tunable */
470  {
471  double X;
472  a = 1 + (n / 2);
473  b = 1 + n - a;
474 
475  X = RAND::Beta((double) a, (double) b);
476 
477  if (X >= p)
478  {
479  n = a - 1;
480  p /= X;
481  }
482  else
483  {
484  k += a;
485  n = b - 1;
486  p = (p - X) / (1 - X);
487  }
488  }
489 
490  for (i = 0; i < n; i++)
491  {
492  double u = RAND::Uniform();
493  if (u < p)
494  k++;
495  }
496 
497  return k;
498 
499  }
500 
501  static inline double Beta (const double a, const double b)
502  {
503  /*from Knuth*/
504  double x1 = RAND::Gamma(a, 1.0);
505  double x2 = RAND::Gamma(b, 1.0);
506 
507  return x1 / (x1 + x2);
508  }
509 #endif
510 
511 
512  static inline void Sample (int from, int to, unsigned int num, int* result, bool replace)
513  {
514  assert(from < to);
515 
516  unsigned int seq_length = to - from; //don't include last (to)
517 
518  assert(num <= seq_length);
519 
520  int *seq = new int [seq_length];
521 
522  seq[0] = from;
523 
524  for(unsigned int i = 1; seq[i-1] < to && i < seq_length; ++i) seq[i] = seq[i-1] + 1;
525 
526  if(!replace) { //without replacement
527 
528  unsigned int size = seq_length, pos, last = seq_length - 1;
529 
530  for (unsigned int i = 0; i < num; i++) {
531  pos = RAND::Uniform(size);
532  result[i] = seq[pos];
533  seq[pos] = seq[last];
534 // seq[last] = result[i]; useless operation
535  size--; last--;
536  }
537 
538  } else { //with replacement
539  for (unsigned int i = 0; i < num; i++) result[i] = seq[ RAND::Uniform(seq_length) ];
540  }
541 
542  delete [] seq;
543  }
544 
545  static inline void SampleSeq(int from, int to, int by, unsigned int num, int* result, bool replace = false)
546  {
547  assert(from < to && by < to && by > 0);
548 
549  unsigned int seq_length = (int)((to - from) / by); //don't include last (to)
550 
551  assert(num <= seq_length);
552 
553  int *seq = new int [seq_length];
554 
555  seq[0] = from;
556 
557  for(unsigned int i = 1; seq[i-1] + by < to && i < seq_length; ++i) seq[i] = seq[i-1] + by;
558 
559  if(!replace) { //without replacement
560 
561  unsigned int size = seq_length, pos, last = seq_length - 1;
562 
563  for (unsigned int i = 0; i < num; i++) {
564  pos = RAND::Uniform(size);
565  result[i] = seq[pos];
566  seq[pos] = seq[last];
567 // seq[last] = result[i]; useless operation
568  size--; last--;
569  }
570 
571  } else { //with replacement
572  for (unsigned int i = 0; i < num; i++) result[i] = seq[ RAND::Uniform(seq_length) ];
573  }
574 
575  delete [] seq;
576  }
577 
578  static inline void SampleSeqWithReciprocal(int from, int to, int by, unsigned int num1, int* result1, unsigned int num2, int* result2)
579  {
580  assert(from < to && by < to && by > 0);
581 
582  unsigned int seq_length = (int)((to - from) / by); //don't include last (to)
583 
584  assert(num1 + num2 == seq_length);
585 
586  int *seq = new int [seq_length];
587 
588  seq[0] = from;
589 
590  for(unsigned int i = 1; seq[i-1] + by < to && i < seq_length; ++i) seq[i] = seq[i-1] + by;
591 
592  unsigned int size = seq_length, pos, last = seq_length - 1;
593 
594  for (unsigned int i = 0; i < num1; i++) {
595  pos = RAND::Uniform(size);
596  result1[i] = seq[pos];
597  seq[pos] = seq[last];
598  seq[last] = result1[i];
599  size--; last--;
600  }
601 
602  for (unsigned int i = 0; i < num2 && i < size; i++)
603  result2[i] = seq[i];
604 
605  delete [] seq;
606  }
607 
608  };
609 
610 #endif
static double Exponential(double mu)
Definition: Uniform.h:424
static double LogNormal(double zeta, double sigma)
Definition: Uniform.h:334
static void BivariateGaussian(double sigma1, double sigma2, double rho, double *out1, double *out2)
Definition: Uniform.h:307
static double Poisson(double mean)
From the Numerical Recieps.
Definition: Uniform.h:196
static unsigned int Uniform(unsigned int max)
Returns a uniformly distributed random number from [0.0, max[.
Definition: Uniform.h:133
static void init(unsigned long seed)
Initialize the random generator's seed.
Definition: Uniform.h:80
static void SampleSeqWithReciprocal(int from, int to, int by, unsigned int num1, int *result1, unsigned int num2, int *result2)
Definition: Uniform.h:578
static unsigned long RandULong()
Return a random unsigned long, from uniform distribution.
Definition: Uniform.h:168
static double Gaussian(double sigma)
Definition: Uniform.h:238
static double Uniform()
Generates a random number from [0.0, 1.0[ uniformly distributed.
Definition: Uniform.h:101
static long Seed2
Definition: Uniform.h:76
static long Seed1
Definition: Uniform.h:76
static double Bernoulli(double p)
Definition: Uniform.h:406
static void Sample(int from, int to, unsigned int num, int *result, bool replace)
Definition: Uniform.h:512
static void SampleSeq(int from, int to, int by, unsigned int num, int *result, bool replace=false)
Definition: Uniform.h:545
static double Gamma(double a, double b)
Definition: Uniform.h:362
static bool RandBool()
Returns a random boolean.
Definition: Uniform.h:139
static double gammln(double xx)
From the Numerical Recieps.
Definition: Uniform.h:181
Random number generation class, uses various types of random generators depending on the implementati...
Definition: Uniform.h:52

Generated for Nemo v2.3.0 by  doxygen 1.8.8 --
Catalogued on GSR