A.I.M.S algorithms


random.h
Go to the documentation of this file.
1 /* This software and supporting documentation are distributed by
2  * Institut Federatif de Recherche 49
3  * CEA/NeuroSpin, Batiment 145,
4  * 91191 Gif-sur-Yvette cedex
5  * France
6  *
7  * This software is governed by the CeCILL-B license under
8  * French law and abiding by the rules of distribution of free software.
9  * You can use, modify and/or redistribute the software under the
10  * terms of the CeCILL-B license as circulated by CEA, CNRS
11  * and INRIA at the following URL "http://www.cecill.info".
12  *
13  * As a counterpart to the access to the source code and rights to copy,
14  * modify and redistribute granted by the license, users are provided only
15  * with a limited warranty and the software's author, the holder of the
16  * economic rights, and the successive licensors have only limited
17  * liability.
18  *
19  * In this respect, the user's attention is drawn to the risks associated
20  * with loading, using, modifying and/or developing or reproducing the
21  * software by the user in light of its specific status of free software,
22  * that may mean that it is complicated to manipulate, and that also
23  * therefore means that it is reserved for developers and experienced
24  * professionals having in-depth computer knowledge. Users are therefore
25  * encouraged to load and test the software's suitability as regards their
26  * requirements in conditions enabling the security of their systems and/or
27  * data to be ensured and, more generally, to use and operate it in the
28  * same conditions as regards security.
29  *
30  * The fact that you are presently reading this means that you have had
31  * knowledge of the CeCILL-B license and that you accept its terms.
32  */
33 
34 
35 #ifndef AIMS_MATH_RANDOM_H
36 #define AIMS_MATH_RANDOM_H
37 
38 #include <cstdlib>
40 #include <stdlib.h>
41 #include <time.h>
42 #include <math.h>
43 #include <complex>
44 
52 
55 template <class T>
57 T UniformRandom(const T &min,const T &max);
60 std::complex<float> UniformRandom(const std::complex<float> &min,
61  const std::complex<float> &max);
64 std::complex<double> UniformRandom(const std::complex<double> &min,
65  const std::complex<double> &max);
67 
68 
71 template <class T>
73 T GaussianRandom(const T &mean,const T &sigma);
76 std::complex<float> GaussianRandom(const std::complex<float> &mean,
77  const std::complex<float> &sigma);
80 std::complex<double> GaussianRandom(const std::complex<double> &mean,
81  const std::complex<double> &sigma);
83 
84 
85 template <class T> class AimsData;
86 
92 
93 template <class T> inline
94 T UniformRandom(const T &min,const T &max)
95 { return((T) ((double)min+(double)(max-min)*(double)UniformRandom()));
96 }
97 
98 
99 inline
100 short UniformRandom( const short& min, const short& max )
101 {
102  return min + short( UniformRandom() * (max-min+0.99999) );
103 }
104 
105 
106 inline
107 unsigned short
108 UniformRandom( const unsigned short& min, const unsigned short& max )
109 {
110  return min + (unsigned short)( UniformRandom() * (max-min+0.99999) );
111 }
112 
113 
114 inline
115 int UniformRandom( const int& min, const int& max )
116 {
117  return min + int( UniformRandom() * (max-min+0.99999) );
118 }
119 
120 
121 inline
122 unsigned int
123 UniformRandom( const unsigned int& min, const unsigned int& max )
124 {
125  return min + (unsigned int)( UniformRandom() * (max-min+0.99999) );
126 }
127 
128 
129 inline
130 long UniformRandom( const long& min, const long& max )
131 {
132  return min + long( UniformRandom() * (max-min+0.99999) );
133 }
134 
135 
136 inline
137 unsigned long
138 UniformRandom( const unsigned long& min, const unsigned long& max )
139 {
140  return min + (unsigned long)( UniformRandom() * (max-min+0.99999) );
141 }
142 
143 
144 template <class T> inline
145 T GaussianRandom(const T &mean,const T &sigma)
146 { return((T)((double)mean + (double)sigma * NormalRandom()));
147 }
148 
149 
150 
151 
152 
154 
155 
157 template< class T >
159 {
160 public:
161 
164  UniformDeviates() { }
167  virtual ~UniformDeviates() { }
169 
177  T ran0( int *idum=NULL );
178 
187  T ran1( int *idum=NULL );
188 
197  T ran2( int *idum=NULL );
198 
204  T ran3( int *idum=NULL );
205 };
206 
207 
208 template< class T > inline
210 {
211  int k;
212  T ans;
213  int idumm;
214 
215  if ( idum == NULL )
216  {
217  srand( time( NULL ) );
218  idumm = - rand();
219  }
220  else idumm = *idum;
221 
222  idumm ^= 123459876;
223  k = idumm / 127773;
224  idumm = 16807 * ( idumm - k * 127773 ) - 2836 * k;
225 
226  if ( idumm < 0 ) idumm += 2147483647;
227 
228  ans = (T)( ( 1.0 / 2147483647 ) * idumm );
229  idumm ^= 123459876;
230 
231  if ( idum ) *idum = idumm;
232 
233  return ans;
234 }
235 
236 
237 template< class T > inline
239 {
240  int j;
241  int k, idumm;
242  static int iy = 0;
243  static int iv[ 32 ];
244  T temp;
245 
246  if ( idum == NULL )
247  {
248  srand( time( NULL ) );
249  idumm = - rand();
250  }
251  else idumm = *idum;
252 
253  if ( idumm <= 0 || !iy )
254  {
255  if ( -idumm < 1 ) idumm = 1;
256  else idumm = -idumm;
257 
258  for ( j=39; j>=0; j-- )
259  {
260  k = idumm / 127773;
261  idumm = 16807 * ( idumm - k * 127773 ) - 2836 * k;
262 
263  if ( idumm < 0 ) idumm += 2147483647;
264  if ( j < 32 ) iv[ j ] = idumm;
265  }
266 
267  iy = iv[ 0 ];
268  }
269 
270  k = idumm / 127773;
271  idumm = 16807 * ( idumm - k * 127773 ) - 2836 * k;
272 
273  if ( idumm < 0 ) idumm += 2147483647;
274 
275  j = iy / ( 1 + 2147483646 / 32 );
276  iy = iv[ j ];
277  iv[ j ] = idumm;
278 
279  if ( idum ) *idum = idumm;
280 
281  if ( ( temp = (T)( 1.0 * iy / 2147483647 ) ) > (T)( 1.0 - 1.2e-7 ) )
282  return (T)( 1.0 - 1.2e-7 );
283  else return (T)temp;
284 }
285 
286 
287 template< class T > inline
289 {
290  int j;
291  int k, idumm;
292  static int idum2 = 123456789;
293  static int iy = 0;
294  static int iv[ 32 ];
295  T temp;
296 
297  if ( idum == NULL )
298  {
299  srand( time( NULL ) );
300  idumm = - rand();
301  }
302  else idumm = *idum;
303 
304  if ( idumm <= 0 )
305  {
306  if ( -idumm < 1 ) idumm = 1;
307  else idumm = -idumm;
308 
309  idum2 = idumm;
310 
311  for ( j = 39; j>=0; j-- )
312  {
313  k = idumm / 53668;
314  idumm = 40014 * ( idumm - k * 53668 ) - k * 12211;
315 
316  if ( idumm < 0 ) idumm += 2147483563;
317  if ( j < 32 ) iv[ j ] = idumm;
318  }
319 
320  iy = iv[ 0 ];
321  }
322 
323  k = idumm / 53668;
324  idumm = 40014 * ( idumm - k * 53668 ) - k * 12211;
325 
326  if ( idumm < 0 ) idumm += 2147483563;
327 
328  k = idum2 / 52774;
329  idum2 = 40692 * ( idum2 - k * 52774 ) - k * 3791;
330 
331  if ( idum2 < 0 ) idum2 += 2147483399;
332 
333  j = iy / ( 1 + 2147483562 / 32 );
334  iy = iv[ j ] - idum2;
335  iv[ j ] = idumm;
336  if ( iy < 1 ) iy += 2147483562;
337 
338  if ( idum ) *idum = idumm;
339 
340  if ( ( temp = (T)( 1.0 * iy / 2147483563 ) ) > (T)( 1.0 - 1.2e-7 ) )
341  return (T)( 1.0 - 1.2e-7 );
342  else return (T)temp;
343 }
344 
345 
346 template< class T > inline
348 {
349  static int inext, inextp;
350  static int ma[ 56 ];
351  static int iff = 0;
352  int mj, mk, idumm;
353  int i, ii, k;
354 
355  if ( idum == NULL )
356  {
357  srand( time( NULL ) );
358  idumm = - rand();
359  }
360  else idumm = *idum;
361 
362  if ( idumm < 0 || iff == 0 )
363  {
364  iff = 1;
365  mj = 161803398 - ( idumm < 0 ? -idumm : idumm );
366  mj %= 1000000000;
367  ma[ 55 ] = mj;
368  mk = 1;
369 
370  for ( i=1; i<=54; i++ )
371  {
372  ii = ( 21 * i ) % 55;
373  ma[ ii ] = mk;
374  mk = mj - mk;
375  if ( mk < 0 ) mk += 1000000000;
376  mj = ma[ ii ];
377  }
378 
379  for ( k=1; k<=4; k++ )
380  for ( i=1; i<=55; i++ )
381  {
382  ma[ i ] -= ma[ 1 + ( i + 30 ) % 55 ];
383 
384  if ( ma[ i ] < 0 ) ma[ i ] += 1000000000;
385  }
386 
387  inext = 0;
388  inextp = 31;
389  idumm = 1;
390  }
391 
392  if ( ++inext == 56 ) inext = 1;
393  if ( ++inextp == 56 ) inextp = 1;
394 
395  mj = ma[ inext ] - ma[ inextp ];
396 
397  if ( mj < 0 ) mj += 1000000000;
398 
399  ma[ inext ] = mj;
400 
401  if ( idum ) *idum = idumm;
402 
403  return (T)( 1.0 * mj / 1000000000 );
404 }
405 
406 
408 template< class T >
410 {
411 public:
412 
415  ExponentialDeviates() { }
418  virtual ~ExponentialDeviates() { }
420 
427  T doit( int *idum=NULL );
428 };
429 
430 
431 template< class T > inline
433 {
434  T dum;
436 
437  do dum = udev.ran1( idum );
438  while ( dum == (T)0 );
439 
440  return (T)( -log( (double)dum ) );
441 }
442 
443 
445 template< class T >
447 {
448 public:
449 
452  NormalDeviates() { }
455  virtual ~NormalDeviates() { }
457 
464  T doit( int *idum=NULL );
465 };
466 
467 
468 template< class T > inline
470 {
471  static int iset = 0;
472  static T gset;
473  T fac, rsq, v1, v2;
475 
476  if ( iset == 0 )
477  {
478  do
479  {
480  v1 = (T)( 2.0 * udev.ran1( idum ) - 1.0 );
481  v2 = (T)( 2.0 * udev.ran1( idum ) - 1.0 );
482  rsq = v1 * v1 + v2 * v2;
483  }
484  while ( rsq >= (T)1 || rsq == (T)0 );
485 
486  fac = (T)sqrt( -2.0 * log( (double)rsq ) / rsq );
487  gset = v1 * fac;
488  iset = 1;
489  return v2 * fac;
490  }
491  else
492  {
493  iset = 0;
494  return gset;
495  }
496 }
497 
498 #endif
float min(float x, float y)
Definition: thickness.h:105
T ran0(int *idum=NULL)
'Minimal' random number generator of Park and Miller.
Definition: random.h:209
NormalDeviates()
Constructors and destructors.
Definition: random.h:453
Exponential random numbers.
Definition: random.h:409
float max(float x, float y)
Definition: thickness.h:97
AIMSALGOPUB_API AimsData< int > AimsRandomList(int size)
generate a random list of int between 0 and (size-1)
AIMSALGOPUB_API double UniformRandom()
Uniform distribution between 0.0 and 1.0.
T GaussianRandom(const T &mean, const T &sigma)
generate gaussian random number with mean,sigma
Definition: random.h:145
T doit(int *idum=NULL)
Normal random numbers generator.
Definition: random.h:469
UniformDeviates()
Constructor and destructors.
Definition: random.h:165
ExponentialDeviates()
Constructor and destructor.
Definition: random.h:416
T ran1(int *idum=NULL)
'Minimal' random number generator of Park and Miller with Bays-Durham shuffle and added safeguards...
Definition: random.h:238
virtual ~ExponentialDeviates()
destructor
Definition: random.h:418
#define AIMSALGOPUB_API
T ran3(int *idum=NULL)
Long period random number generator of Knuth.
Definition: random.h:347
T doit(int *idum=NULL)
Exponential random number generator.
Definition: random.h:432
Normal random numbers.
Definition: random.h:446
T ran2(int *idum=NULL)
Long period random number generator of L'Ecuyer with Bays-Derham shuffle and added safeguards...
Definition: random.h:288
AIMSALGOPUB_API double NormalRandom()
Normal gaussian distribution with mean = 0.0 and sigma = 1.0.
virtual ~NormalDeviates()
destructor
Definition: random.h:455
virtual ~UniformDeviates()
destructor
Definition: random.h:167