ARTS  2.4.0(git:4fb77825)
rng.cc
Go to the documentation of this file.
1 /* Copyright (C) 2003-2012 Cory Davis <cory@met.ed.ac.uk>
2 
3  This program is free software; you can redistribute it and/or modify it
4  under the terms of the GNU General Public License as published by the
5  Free Software Foundation; either version 2, or (at your option) any
6  later version.
7 
8  This program is distributed in the hope that it will be useful,
9  but WITHOUT ANY WARRANTY; without even the implied warranty of
10  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11  GNU General Public License for more details.
12 
13  You should have received a copy of the GNU General Public License
14  along with this program; if not, write to the Free Software
15  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16  USA. */
17 
31 #include "rng.h"
32 #include <algorithm>
33 #include <climits>
34 #include <cstddef>
35 #include <cstdio>
36 #include <cstdlib>
37 #include <cstring>
38 #include <iostream>
39 #include <vector>
40 
41 #include "arts.h"
42 #include "messages.h"
43 
48 
53 
54 void Rng::seed(unsigned long int n, const Verbosity &verbosity) {
55  // Static pool of previously used seeds.
56  static vector<unsigned long int> seeds;
57 
58 #pragma omp critical(Rng_seed)
59  {
60  unsigned long int n_orig = n;
61  //cout << "Requested seed: " << n;
62  while (find(seeds.begin(), seeds.end(), n) != seeds.end()) {
63  if (n >= ULONG_MAX - 1)
64  n = 0;
65  else
66  n++;
67 
68  // If all possible seeds were already used, we empty the pool and
69  // start over.
70  if (n == n_orig) {
72  out0
73  << "Rng Warning: Couldn't find an unused seed. Clearing seed pool.\n";
74  seeds.empty();
75  break;
76  }
77  }
78  seeds.push_back(n);
79  seed_no = n;
80  //cout << " Got seed: " << seed_no << endl;
81  }
82 
84 }
85 
89 void Rng::force_seed(unsigned long int n) {
90  seed_no = n;
92 }
93 
97 double Rng::draw() { return gsl_rng_uniform(r); }
98 
102 unsigned long int Rng::showseed() const { return seed_no; }
103 
104 /*
105  * rng/mt.c
106  *
107  * This program is free software; you can redistribute it and/or
108  * modify it under the terms of the GNU General Public License as
109  * published by the Free Software Foundation; either version 2 of the
110  * License, or (at your option) any later version.
111  *
112  * This program is distributed in the hope that it will be useful, but
113  * WITHOUT ANY WARRANTY; without even the implied warranty of
114  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
115  * General Public License for more details. You should have received
116  * a copy of the GNU General Public License along with this program;
117  * if not, write to the Free Foundation, Inc., 59 Temple Place, Suite
118  * 330, Boston, MA 02111-1307 USA.
119  *
120  * Original implementation was copyright (C) 1997 Makoto Matsumoto and
121  * Takuji Nishimura. Coded by Takuji Nishimura, considering the
122  * suggestions by Topher Cooper and Marc Rieffel in July-Aug. 1997, "A
123  * C-program for MT19937: Integer version (1998/4/6)".
124  *
125  * This implementation copyright (C) 1998 Brian Gough. I reorganized
126  * the code to use the module framework of GSL. The license on this
127  * implementation was changed from LGPL to GPL, following paragraph 3
128  * of the LGPL, version 2.
129  *
130  * Update:
131  *
132  * The seeding procedure has been updated to match the 10/99 release
133  * of MT19937.
134  *
135  * Update:
136  *
137  * The seeding procedure has been updated again to match the 2002
138  * release of MT19937.
139  *
140  * The original code included the comment: "When you use this, send an
141  * email to: matumoto@math.keio.ac.jp with an appropriate reference to
142  * your work".
143  *
144  * Makoto Matsumoto has a web page with more information about the
145  * generator, http://www.math.keio.ac.jp/~matumoto/emt.html.
146  *
147  * The paper below has details of the algorithm.
148  *
149  * From: Makoto Matsumoto and Takuji Nishimura, "Mersenne Twister: A
150  * 623-dimensionally equidistributerd uniform pseudorandom number
151  *generator". ACM Transactions on Modeling and Computer Simulation,
152  * Vol. 8, No. 1 (Jan. 1998), Pages 3-30.
153  *
154  * You can obtain the paper directly from Makoto Matsumoto's web page.
155  *
156  * The period of this generator is 2^{19937} - 1.
157  *
158  */
159 
160 static inline unsigned long int mt_get(void *vstate);
161 static double mt_get_double(void *vstate);
162 static void mt_set(void *state, unsigned long int s);
163 
164 #define N 624 /* Period parameters */
165 #define M 397
166 
167 /* most significant w-r bits */
168 static const unsigned long UPPER_MASK = 0x80000000UL;
169 
170 /* least significant r bits */
171 static const unsigned long LOWER_MASK = 0x7fffffffUL;
172 
173 typedef struct {
174  unsigned long mt[N];
175  int mti;
176 } mt_state_t;
177 
178 static inline unsigned long mt_get(void *vstate) {
179  mt_state_t *state = (mt_state_t *)vstate;
180 
181  unsigned long k;
182  unsigned long int *const mt = state->mt;
183 
184 #define MAGIC(y) (((y)&0x1) ? 0x9908b0dfUL : 0)
185 
186  if (state->mti >= N) { /* generate N words at one time */
187  int kk;
188 
189  for (kk = 0; kk < N - M; kk++) {
190  unsigned long y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
191  mt[kk] = mt[kk + M] ^ (y >> 1) ^ MAGIC(y);
192  }
193  for (; kk < N - 1; kk++) {
194  unsigned long y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
195  mt[kk] = mt[kk + (M - N)] ^ (y >> 1) ^ MAGIC(y);
196  }
197 
198  {
199  unsigned long y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
200  mt[N - 1] = mt[M - 1] ^ (y >> 1) ^ MAGIC(y);
201  }
202 
203  state->mti = 0;
204  }
205 
206  /* Tempering */
207 
208  k = mt[state->mti];
209  k ^= (k >> 11);
210  k ^= (k << 7) & 0x9d2c5680UL;
211  k ^= (k << 15) & 0xefc60000UL;
212  k ^= (k >> 18);
213 
214  state->mti++;
215 
216  return k;
217 }
218 
219 static double mt_get_double(void *vstate) {
220  return (double)mt_get(vstate) / 4294967296.0;
221 }
222 
223 static void mt_set(void *vstate, unsigned long int s) {
224  mt_state_t *state = (mt_state_t *)vstate;
225  int i;
226 
227  if (s == 0) s = 4357; /* the default seed is 4357 */
228 
229  state->mt[0] = s & 0xffffffffUL;
230 
231  for (i = 1; i < N; i++) {
232  /* See Knuth's "Art of Computer Programming" Vol. 2, 3rd
233  Ed. p.106 for multiplier. */
234 
235  state->mt[i] =
236  (1812433253UL * (state->mt[i - 1] ^ (state->mt[i - 1] >> 30)) + i);
237 
238  state->mt[i] &= 0xffffffffUL;
239  }
240 
241  state->mti = i;
242 }
243 
244 static const gsl_rng_type mt_type = {"mt19937", /* name */
245  0xffffffffUL, /* RAND_MAX */
246  0, /* RAND_MIN */
247  sizeof(mt_state_t),
248  &mt_set,
249  &mt_get,
250  &mt_get_double};
251 
252 const gsl_rng_type *gsl_rng_mt19937 = &mt_type;
253 unsigned long int gsl_rng_default_seed = 0;
254 /* rng/types.c
255  *
256  * Copyright (C) 2001 Brian Gough
257  *
258  * This program is free software; you can redistribute it and/or modify
259  * it under the terms of the GNU General Public License as published by
260  * the Free Software Foundation; either version 2 of the License, or (at
261  * your option) any later version.
262  *
263  * This program is distributed in the hope that it will be useful, but
264  * WITHOUT ANY WARRANTY; without even the implied warranty of
265  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
266  * General Public License for more details.
267  *
268  * You should have received a copy of the GNU General Public License
269  * along with this program; if not, write to the Free Software
270  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
271  */
272 
273 #define N1 100
274 
276 
277 #define ADD(t) \
278  { \
279  if (i == N1) abort(); \
280  gsl_rng_generator_types[i] = (t); \
281  i++; \
282  };
283 
285  int i = 0;
286 
288  ADD(0);
289 
291 }
292 
293 /* err/stream.c
294  *
295  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
296  *
297  * This program is free software; you can redistribute it and/or modify
298  * it under the terms of the GNU General Public License as published by
299  * the Free Software Foundation; either version 2 of the License, or (at
300  * your option) any later version.
301  *
302  * This program is distributed in the hope that it will be useful, but
303  * WITHOUT ANY WARRANTY; without even the implied warranty of
304  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
305  * General Public License for more details.
306  *
307  * You should have received a copy of the GNU General Public License
308  * along with this program; if not, write to the Free Software
309  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
310  */
311 
312 FILE *gsl_stream = NULL;
314 
315 void gsl_stream_printf(const char *label,
316  const char *file,
317  int line,
318  const char *reason) {
319  if (gsl_stream == NULL) {
320  gsl_stream = stderr;
321  }
322  if (gsl_stream_handler) {
323  (*gsl_stream_handler)(label, file, line, reason);
324  return;
325  }
326  fprintf(gsl_stream, "gsl: %s:%d: %s: %s\n", file, line, label, reason);
327 }
328 
330  gsl_stream_handler_t *new_handler) {
331  gsl_stream_handler_t *previous_handler = gsl_stream_handler;
332  gsl_stream_handler = new_handler;
333  return previous_handler;
334 }
335 
336 FILE *gsl_set_stream(FILE *new_stream) {
337  FILE *previous_stream;
338  if (gsl_stream == NULL) {
339  gsl_stream = stderr;
340  }
341  previous_stream = gsl_stream;
342  gsl_stream = new_stream;
343  return previous_stream;
344 }
345 
346 /* err/error.c
347  *
348  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
349  *
350  * This program is free software; you can redistribute it and/or modify
351  * it under the terms of the GNU General Public License as published by
352  * the Free Software Foundation; either version 2 of the License, or (at
353  * your option) any later version.
354  *
355  * This program is distributed in the hope that it will be useful, but
356  * WITHOUT ANY WARRANTY; without even the implied warranty of
357  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
358  * General Public License for more details.
359  *
360  * You should have received a copy of the GNU General Public License
361  * along with this program; if not, write to the Free Software
362  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
363  */
364 
366 
367 static void no_error_handler(const char *reason,
368  const char *file,
369  int line,
370  int gsl_errno);
371 
372 void gsl_error(const char *reason, const char *file, int line, int gsl_errno) {
373  if (gsl_error_handler) {
374  (*gsl_error_handler)(reason, file, line, gsl_errno);
375  return;
376  }
377 
378  gsl_stream_printf("ERROR", file, line, reason);
379 
380  fprintf(stderr, "Default GSL error handler invoked.\n");
381  abort();
382 }
383 
385  gsl_error_handler_t *previous_handler = gsl_error_handler;
386  gsl_error_handler = new_handler;
387  return previous_handler;
388 }
389 
391  gsl_error_handler_t *previous_handler = gsl_error_handler;
392  gsl_error_handler = no_error_handler;
393  return previous_handler;
394 }
395 
396 static void no_error_handler(const char *reason _U_,
397  const char *file _U_,
398  int line _U_,
399  int gsl_errno _U_) {
400  /* do nothing */
401  return;
402 }
403 
404 /* rng/rng.c
405  *
406  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
407  *
408  * This program is free software; you can redistribute it and/or modify
409  * it under the terms of the GNU General Public License as published by
410  * the Free Software Foundation; either version 2 of the License, or (at
411  * your option) any later version.
412  *
413  * This program is distributed in the hope that it will be useful, but
414  * WITHOUT ANY WARRANTY; without even the implied warranty of
415  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
416  * General Public License for more details.
417  *
418  * You should have received a copy of the GNU General Public License
419  * along with this program; if not, write to the Free Software
420  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
421  */
422 
424  gsl_rng *r = (gsl_rng *)malloc(sizeof(gsl_rng));
425 
426  if (r == 0) {
427  GSL_ERROR_VAL("failed to allocate space for rng struct", GSL_ENOMEM, 0);
428  };
429 
430  r->state = malloc(T->size);
431 
432  if (r->state == 0) {
433  free(r); /* exception in constructor, avoid memory leak */
434 
435  GSL_ERROR_VAL("failed to allocate space for rng state", GSL_ENOMEM, 0);
436  };
437 
438  r->type = T;
439 
440  gsl_rng_set(r, gsl_rng_default_seed); /* seed the generator */
441 
442  return r;
443 }
444 
445 int gsl_rng_memcpy(gsl_rng *dest, const gsl_rng *src) {
446  if (dest->type != src->type) {
447  GSL_ERROR("generators must be of the same type", GSL_EINVAL);
448  }
449 
450  memcpy(dest->state, src->state, src->type->size);
451 
452  return GSL_SUCCESS;
453 }
454 
456  gsl_rng *r = (gsl_rng *)malloc(sizeof(gsl_rng));
457 
458  if (r == 0) {
459  GSL_ERROR_VAL("failed to allocate space for rng struct", GSL_ENOMEM, 0);
460  };
461 
462  r->state = malloc(q->type->size);
463 
464  if (r->state == 0) {
465  free(r); /* exception in constructor, avoid memory leak */
466 
467  GSL_ERROR_VAL("failed to allocate space for rng state", GSL_ENOMEM, 0);
468  };
469 
470  r->type = q->type;
471 
472  memcpy(r->state, q->state, q->type->size);
473 
474  return r;
475 }
476 
477 void gsl_rng_set(const gsl_rng *r, unsigned long int seed) {
478  (r->type->set)(r->state, seed);
479 }
480 
481 #ifndef HIDE_INLINE_STATIC
482 unsigned long int gsl_rng_get(const gsl_rng *r) {
483  return (r->type->get)(r->state);
484 }
485 
486 double gsl_rng_uniform(const gsl_rng *r) {
487  return (r->type->get_double)(r->state);
488 }
489 
490 double gsl_rng_uniform_pos(const gsl_rng *r) {
491  double x;
492  do {
493  x = (r->type->get_double)(r->state);
494  } while (x == 0);
495 
496  return x;
497 }
498 
499 unsigned long int gsl_rng_uniform_int(const gsl_rng *r, unsigned long int n) {
500  unsigned long int offset = r->type->min;
501  unsigned long int range = r->type->max - offset;
502  unsigned long int scale = range / n;
503  unsigned long int k;
504 
505  if (n > range) {
506  GSL_ERROR_VAL("n exceeds maximum value of generator", GSL_EINVAL, 0);
507  }
508 
509  do {
510  k = (((r->type->get)(r->state)) - offset) / scale;
511  } while (k >= n);
512 
513  return k;
514 }
515 #endif
516 
517 unsigned long int gsl_rng_max(const gsl_rng *r) { return r->type->max; }
518 
519 unsigned long int gsl_rng_min(const gsl_rng *r) { return r->type->min; }
520 
521 const char *gsl_rng_name(const gsl_rng *r) { return r->type->name; }
522 
523 size_t gsl_rng_size(const gsl_rng *r) { return r->type->size; }
524 
525 void *gsl_rng_state(const gsl_rng *r) { return r->state; }
526 
527 void gsl_rng_print_state(const gsl_rng *r) {
528  size_t i;
529  unsigned char *p = (unsigned char *)(r->state);
530  const size_t n = r->type->size;
531 
532  for (i = 0; i < n; i++) {
533  /* FIXME: we're assuming that a char is 8 bits */
534  printf("%.2x", *(p + i));
535  }
536 }
537 
539  free(r->state);
540  free(r);
541 }
gsl_rng_print_state
void gsl_rng_print_state(const gsl_rng *r)
Definition: rng.cc:527
mt_state_t::mt
unsigned long mt[N]
Definition: rng.cc:174
gsl_rng_max
unsigned long int gsl_rng_max(const gsl_rng *r)
Definition: rng.cc:517
gsl_rng_min
unsigned long int gsl_rng_min(const gsl_rng *r)
Definition: rng.cc:519
Rng::showseed
unsigned long int showseed() const
Returns the seed number.
Definition: rng.cc:102
gsl_set_stream_handler
gsl_stream_handler_t * gsl_set_stream_handler(gsl_stream_handler_t *new_handler)
Definition: rng.cc:329
Rng::~Rng
~Rng()
Destructor frees memory allocated to gsl_rng.
Definition: rng.cc:52
Rng::seed_no
unsigned long int seed_no
Definition: rng.h:557
gsl_set_stream
FILE * gsl_set_stream(FILE *new_stream)
Definition: rng.cc:336
gsl_rng::type
const gsl_rng_type * type
Definition: rng.h:129
gsl_rng_types_setup
const gsl_rng_type ** gsl_rng_types_setup(void)
Definition: rng.cc:284
N1
#define N1
Definition: rng.cc:273
ARTS::Var::verbosity
Verbosity verbosity(Workspace &ws) noexcept
Definition: autoarts.h:7112
gsl_rng_mt19937
const gsl_rng_type * gsl_rng_mt19937
Definition: rng.cc:252
GSL_SUCCESS
@ GSL_SUCCESS
Definition: rng.h:313
Rng::seed
void seed(unsigned long int n, const Verbosity &verbosity)
Seeds the Rng with the integer argument.
Definition: rng.cc:54
gsl_rng_type::get
unsigned long int(* get)(void *state)
Definition: rng.h:124
ARTS::Var::y
Vector y(Workspace &ws) noexcept
Definition: autoarts.h:7401
gsl_set_error_handler
gsl_error_handler_t * gsl_set_error_handler(gsl_error_handler_t *new_handler)
Definition: rng.cc:384
gsl_rng_uniform_pos
double gsl_rng_uniform_pos(const gsl_rng *r)
Definition: rng.cc:490
gsl_set_error_handler_off
gsl_error_handler_t * gsl_set_error_handler_off(void)
Definition: rng.cc:390
mt_state_t
Definition: rng.cc:173
gsl_rng_get
unsigned long int gsl_rng_get(const gsl_rng *r)
Definition: rng.cc:482
gsl_rng_default_seed
unsigned long int gsl_rng_default_seed
Definition: rng.cc:253
gsl_stream
FILE * gsl_stream
Definition: rng.cc:312
gsl_rng_type::get_double
double(* get_double)(void *state)
Definition: rng.h:125
gsl_rng_clone
gsl_rng * gsl_rng_clone(const gsl_rng *q)
Definition: rng.cc:455
rng.h
member functions of the Rng class and gsl_rng code
GSL_EINVAL
@ GSL_EINVAL
Definition: rng.h:319
Rng::draw
double draw()
Draws a double from the uniform distribution [0,1)
Definition: rng.cc:97
gsl_rng::state
void * state
Definition: rng.h:130
Rng::force_seed
void force_seed(unsigned long int n)
Seeds the Rng with the integer argument.
Definition: rng.cc:89
gsl_stream_handler
gsl_stream_handler_t * gsl_stream_handler
Definition: rng.cc:313
gsl_stream_handler_t
void gsl_stream_handler_t(const char *label, const char *file, int line, const char *reason)
Definition: rng.h:366
CREATE_OUT0
#define CREATE_OUT0
Definition: messages.h:204
gsl_rng_type::min
unsigned long int min
Definition: rng.h:121
messages.h
Declarations having to do with the four output streams.
gsl_rng_type::max
unsigned long int max
Definition: rng.h:120
_U_
#define _U_
Definition: config.h:183
gsl_error_handler
gsl_error_handler_t * gsl_error_handler
Definition: rng.cc:365
gsl_error_handler_t
void gsl_error_handler_t(const char *reason, const char *file, int line, int gsl_errno)
Definition: rng.h:361
GSL_ERROR_VAL
#define GSL_ERROR_VAL(reason, gsl_errno, value)
Definition: rng.h:389
gsl_rng_state
void * gsl_rng_state(const gsl_rng *r)
Definition: rng.cc:525
GSL_ERROR
#define GSL_ERROR(reason, gsl_errno)
Definition: rng.h:381
Verbosity
Definition: messages.h:49
gsl_rng_type::name
const char * name
Definition: rng.h:119
gsl_rng_type::set
void(* set)(void *state, unsigned long int seed)
Definition: rng.h:123
gsl_rng_alloc
gsl_rng * gsl_rng_alloc(const gsl_rng_type *T)
Definition: rng.cc:423
gsl_rng_uniform_int
unsigned long int gsl_rng_uniform_int(const gsl_rng *r, unsigned long int n)
Definition: rng.cc:499
mt_state_t::mti
int mti
Definition: rng.cc:175
N
#define N
Definition: rng.cc:164
Rng::Rng
Rng()
Constructor creates instance of gsl_rng of type gsl_rng_mt19937.
Definition: rng.cc:47
gsl_rng_size
size_t gsl_rng_size(const gsl_rng *r)
Definition: rng.cc:523
gsl_error
void gsl_error(const char *reason, const char *file, int line, int gsl_errno)
Definition: rng.cc:372
q
#define q
Definition: legacy_continua.cc:21712
gsl_rng_generator_types
const gsl_rng_type * gsl_rng_generator_types[N1]
Definition: rng.cc:275
gsl_rng_name
const char * gsl_rng_name(const gsl_rng *r)
Definition: rng.cc:521
GSL_ENOMEM
@ GSL_ENOMEM
Definition: rng.h:323
gsl_rng_free
void gsl_rng_free(gsl_rng *r)
Definition: rng.cc:538
gsl_stream_printf
void gsl_stream_printf(const char *label, const char *file, int line, const char *reason)
Definition: rng.cc:315
Rng::r
gsl_rng * r
Definition: rng.h:555
MAGIC
#define MAGIC(y)
gsl_rng_type::size
size_t size
Definition: rng.h:122
ARTS::Var::x
Vector x(Workspace &ws) noexcept
Definition: autoarts.h:7346
gsl_rng_type
Definition: rng.h:118
ADD
#define ADD(t)
Definition: rng.cc:277
M
#define M
Definition: rng.cc:165
gsl_rng
Definition: rng.h:128
gsl_rng_uniform
double gsl_rng_uniform(const gsl_rng *r)
Definition: rng.cc:486
arts.h
The global header file for ARTS.
gsl_rng_set
void gsl_rng_set(const gsl_rng *r, unsigned long int seed)
Definition: rng.cc:477
gsl_rng_memcpy
int gsl_rng_memcpy(gsl_rng *dest, const gsl_rng *src)
Definition: rng.cc:445