LAL  7.1.7.1-15d842a
Random.c
Go to the documentation of this file.
1 /*
2 * Copyright (C) 2007 Alexander Dietz, Jolien Creighton, Kipp Cannon
3 *
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 2 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with with program; see the file COPYING. If not, write to the
16 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17 * MA 02110-1301 USA
18 */
19 
20 #include <time.h>
21 #include <math.h>
22 #include <lal/LALStdlib.h>
23 #include <lal/Random.h>
24 #include <lal/Sequence.h>
25 #include <lal/XLALError.h>
26 
27 /**
28  * \defgroup Random_c Module Random.c
29  * \ingroup Random_h
30  *
31  * \brief Functions for generating random numbers.
32  *
33  * ### Description ###
34  *
35  * The routines <tt>LALCreateRandomParams()</tt> and <tt>LALDestroyRandomParams()</tt>
36  * create and destroy a parameter structure for the generation of random
37  * variables. The creation routine requires a random number seed \c seed.
38  * If the seed is zero then a seed is generated using the current time.
39  *
40  * The routine <tt>LALUniformDeviate()</tt> returns a single random deviate
41  * distributed uniformly between zero and unity.
42  *
43  * The routine <tt>LALNormalDeviates()</tt> fills a vector with normal (Gaussian)
44  * deviates with zero mean and unit variance, whereas the function\c XLALNormalDeviate just returns one normal distributed random number.
45  *
46  * ### Operating Instructions ###
47  *
48  * \code
49  * static LALStatus status;
50  * static RandomParams *params;
51  * static REAL4Vector *vector;
52  * UINT4 i;
53  * INT4 seed = 0;
54  *
55  * LALCreateVector( &status, &vector, 9999 );
56  * LALCreateRandomParams( &status, &params, seed );
57  *
58  * /\* fill vector with uniform deviates *\/
59  * for ( i = 0; i < vector->length; ++i )
60  * {
61  * LALUniformDeviate( &status, vector->data + i, params );
62  * }
63  *
64  * /\* fill vector with normal deviates *\/
65  * LALNormalDeviates( &status, vector, params );
66  *
67  * LALDestroyRandomParams( &status, &params );
68  * LALDestroyVector( &status, &vector );
69  * \endcode
70  *
71  * ### Algorithm ###
72  *
73  * This is an implementation of the random number generators \c ran1 and
74  * \c gasdev described in Numerical Recipes \cite ptvf1992 .
75  *
76  */
77 /** @{ */
78 
79 static const INT4 a = 16807;
80 static const INT4 m = 2147483647;
81 static const INT4 q = 127773;
82 static const INT4 r = 2836;
83 
84 static const REAL4 eps = 1.2e-7;
85 
86 /*
87  *
88  * XLAL Routines.
89  *
90  */
91 
93 {
94  INT4 k;
95  k = i/q;
96  i = a*(i - k*q) - r*k;
97  if (i < 0)
98  i += m;
99  return i;
100 }
101 
103 {
104  RandomParams *params;
105 
106  params = XLALMalloc( sizeof( *params) );
107  if ( ! params )
109 
110  while ( seed == 0 ) /* use system clock to get seed */
111  seed = time( NULL );
112 
113  if ( seed < 0 )
114  seed = -seed;
115 
116  XLALResetRandomParams( params, seed );
117 
118  return params;
119 }
120 
121 
123 {
124  UINT4 n;
125 
126  params->i = seed;
127  for ( n = 0; n < 8; ++n ) /* shuffle 8 times */
128  params->i = XLALBasicRandom( params->i );
129 
130  /* populate vector of random numbers */
131  for ( n = 0; n < sizeof( params->v )/sizeof( *params->v ); ++n )
132  params->v[n] = params->i = XLALBasicRandom( params->i );
133 
134  /* first random number is the 0th element of v */
135  params->y = params->v[0];
136 
137 }
138 
139 
141 {
142  XLALFree( params );
143 }
144 
145 
147 {
148  REAL4 ans;
149  INT4 ndiv;
150  INT4 n;
151 
152  if ( ! params )
154 
155  /* randomly choose which element of the vector of random numbers to use */
156  ndiv = 1 + (m - 1)/(XLAL_NUM_ELEM(params->v));
157  n = params->y/ndiv;
158  params->y = params->v[n];
159 
160  /* repopulate this element */
161  params->v[n] = params->i = XLALBasicRandom( params->i );
162 
163  ans = params->y/(REAL4)m;
164  if ( ans > 1 - eps ) /* make sure it is not exactly 1 */
165  ans = 1 - eps;
166 
167  return ans;
168 }
169 
170 
171 int XLALNormalDeviates( REAL4Vector *deviates, RandomParams *params )
172 {
173  REAL4 *data;
174  INT4 half;
175 
176  if ( ! deviates || ! deviates->data || ! params )
178  if ( ! deviates->length )
180 
181  data = deviates->data;
182  half = deviates->length/2;
183 
184  while (half-- > 0)
185  {
186  REAL4 u;
187  REAL4 v;
188  REAL4 x;
189  REAL4 y;
190  REAL4 rsq;
191  REAL4 fac;
192 
193  do {
194  u = XLALUniformDeviate( params );
195  v = XLALUniformDeviate( params );
196  x = 2*u - 1;
197  y = 2*v - 1;
198  rsq = x*x + y*y;
199  }
200  while (rsq >= 1 || rsq == 0);
201 
202  fac = sqrt(-2*log(rsq)/rsq);
203  *data++ = fac*x;
204  *data++ = fac*y;
205  }
206 
207  /* do it again if there is an odd amount of data */
208  if (deviates->length % 2)
209  {
210  REAL4 u;
211  REAL4 v;
212  REAL4 x;
213  REAL4 y;
214  REAL4 rsq;
215  REAL4 fac;
216 
217  do {
218  u = XLALUniformDeviate( params );
219  v = XLALUniformDeviate( params );
220  x = 2*u - 1;
221  y = 2*v - 1;
222  rsq = x*x + y*y;
223  }
224  while (rsq >= 1 || rsq == 0);
225 
226  fac = sqrt(-2*log(rsq)/rsq);
227  *data = fac*x;
228  /* throw away y */
229  }
230 
231  return XLAL_SUCCESS;
232 }
233 
235 {
236  REAL4Sequence *deviates;
237  REAL4 deviate;
238 
239  if ( ! params )
241 
242  /* create a vector */
243  deviates = XLALCreateREAL4Sequence(1);
244  if(!deviates)
246 
247  /* call the actual function */
248  XLALNormalDeviates( deviates, params );
249  deviate = deviates->data[0];
250 
251  /* destroy the vector */
252  XLALDestroyREAL4Sequence(deviates);
253 
254  return deviate;
255 }
256 
257 /*
258  *
259  * LAL Routines.
260  *
261  */
262 
263 
264 
265 void
267  LALStatus *status,
268  RandomParams **params,
269  INT4 seed
270  )
271 {
272  INITSTATUS(status);
273 
274  ASSERT (params, status, RANDOMH_ENULL, RANDOMH_MSGENULL);
275  ASSERT (!*params, status, RANDOMH_ENNUL, RANDOMH_MSGENNUL);
276 
277  *params = XLALCreateRandomParams( seed );
278  if ( ! params )
279  {
280  XLALClearErrno();
281  ABORT( status, RANDOMH_ENULL, RANDOMH_MSGENULL );
282  }
283 
284  RETURN (status);
285 }
286 
287 
288 
289 void
291  LALStatus *status,
292  RandomParams **params
293  )
294 {
295  INITSTATUS(status);
296 
297  ASSERT (params, status, RANDOMH_ENULL, RANDOMH_MSGENULL);
298  ASSERT (*params, status, RANDOMH_ENULL, RANDOMH_MSGENULL);
299 
300  XLALDestroyRandomParams( *params );
301  *params = NULL;
302 
303  RETURN (status);
304 }
305 
306 
307 
308 void
310  LALStatus *status,
311  REAL4 *deviate,
312  RandomParams *params
313  )
314 {
315  INITSTATUS(status);
316 
317  ASSERT (deviate, status, RANDOMH_ENULL, RANDOMH_MSGENULL);
318  ASSERT (params, status, RANDOMH_ENULL, RANDOMH_MSGENULL);
319 
320  *deviate = XLALUniformDeviate( params );
321  if ( XLAL_IS_REAL4_FAIL_NAN( *deviate ) )
322  {
323  XLALClearErrno();
324  ABORT( status, RANDOMH_ENULL, RANDOMH_MSGENULL );
325  }
326 
327  RETURN (status);
328 }
329 
330 
331 
332 void
334  LALStatus *status,
335  REAL4Vector *deviates,
336  RandomParams *params
337  )
338 {
339  INITSTATUS(status);
340 
341  ASSERT (params, status, RANDOMH_ENULL, RANDOMH_MSGENULL);
342  ASSERT (deviates, status, RANDOMH_ENULL, RANDOMH_MSGENULL);
343  ASSERT (deviates->data, status, RANDOMH_ENULL, RANDOMH_MSGENULL);
344  ASSERT (deviates->length > 0, status, RANDOMH_ESIZE, RANDOMH_MSGESIZE);
345 
346  if ( XLALNormalDeviates( deviates, params ) != XLAL_SUCCESS )
347  {
348  int errnum = xlalErrno;
349  XLALClearErrno();
350  switch ( errnum )
351  {
352  case XLAL_EFAULT:
353  ABORT( status, RANDOMH_ENULL, RANDOMH_MSGENULL );
354  case XLAL_EBADLEN:
355  ABORT( status, RANDOMH_ESIZE, RANDOMH_MSGESIZE );
356  default:
357  ABORTXLAL( status );
358  }
359  }
360 
361  RETURN (status);
362 }
363 /** @} */
static const REAL4 eps
Definition: Random.c:84
RandomParams * XLALCreateRandomParams(INT4 seed)
Definition: Random.c:102
#define XLALMalloc(n)
Definition: LALMalloc.h:44
#define xlalErrno
Modifiable lvalue containing the XLAL error number.
Definition: XLALError.h:572
#define XLALFree(p)
Definition: LALMalloc.h:47
void LALNormalDeviates(LALStatus *status, REAL4Vector *deviates, RandomParams *params)
Definition: Random.c:333
#define XLAL_ERROR(...)
Macro to invoke a failure from a XLAL routine returning an integer.
Definition: XLALError.h:701
#define RANDOMH_ENNUL
Non-null pointer.
Definition: Random.h:52
#define RANDOMH_ENULL
Null pointer.
Definition: Random.h:51
static const INT4 a
Definition: Random.c:79
float REAL4
Single precision real floating-point number (4 bytes).
static const INT4 m
Definition: Random.c:80
Vector of type REAL4, see DATATYPE-Vector types for more details.
Definition: LALDatatypes.h:145
void LALDestroyRandomParams(LALStatus *status, RandomParams **params)
Definition: Random.c:290
REAL4Sequence * XLALCreateREAL4Sequence(size_t length)
void LALCreateRandomParams(LALStatus *status, RandomParams **params, INT4 seed)
Definition: Random.c:266
#define RETURN(statusptr)
#define INITSTATUS(statusptr)
int XLALClearErrno(void)
Clears the XLAL error number, returns the old value.
Definition: XLALError.c:363
Invalid pointer.
Definition: XLALError.h:409
int XLALNormalDeviates(REAL4Vector *deviates, RandomParams *params)
Definition: Random.c:171
#define RANDOMH_ESIZE
Invalid size.
Definition: Random.h:53
void XLALDestroyREAL4Sequence(REAL4Sequence *sequence)
#define XLAL_ERROR_NULL(...)
Macro to invoke a failure from a XLAL routine returning a pointer.
Definition: XLALError.h:714
void XLALResetRandomParams(RandomParams *params, INT4 seed)
Definition: Random.c:122
#define ABORT(statusptr, code, mesg)
REAL4 * data
Pointer to the data array.
Definition: LALDatatypes.h:150
INT4 i
Definition: Random.h:87
LAL status structure, see The LALStatus structure for more details.
Definition: LALDatatypes.h:947
Success return value (not an error number)
Definition: XLALError.h:402
REAL4 XLALUniformDeviate(RandomParams *params)
Definition: Random.c:146
static const INT4 q
Definition: Random.c:81
uint32_t UINT4
Four-byte unsigned integer.
INT4 y
Definition: Random.h:88
#define XLAL_IS_REAL4_FAIL_NAN(val)
Tests if val is a XLAL REAL4 failure NaN.
Definition: XLALError.h:393
#define ABORTXLAL(sp)
Memory allocation error.
Definition: XLALError.h:408
static const INT4 r
Definition: Random.c:82
int32_t INT4
Four-byte signed integer.
INT4 v[32]
Definition: Random.h:89
Internal function call failed bit: "or" this with existing error number.
Definition: XLALError.h:463
REAL4 XLALNormalDeviate(RandomParams *params)
Definition: Random.c:234
#define ASSERT(assertion, statusptr, code, mesg)
#define XLAL_ERROR_REAL4(...)
Macro to invoke a failure from a XLAL routine returning a REAL4.
Definition: XLALError.h:740
UINT4 length
Number of elements in array.
Definition: LALDatatypes.h:149
void LALUniformDeviate(LALStatus *status, REAL4 *deviate, RandomParams *params)
Definition: Random.c:309
Inconsistent or invalid length.
Definition: XLALError.h:420
void XLALDestroyRandomParams(RandomParams *params)
Definition: Random.c:140
#define XLAL_NUM_ELEM(x)
MACRO which gives the number of elements in a fixed-size array.
This structure contains the parameters necessary for generating the current sequence of random number...
Definition: Random.h:84
INT4 XLALBasicRandom(INT4 i)
Definition: Random.c:92