source: trunk/yat/random/random.cc @ 2768

Last change on this file since 2768 was 2768, checked in by Peter, 10 years ago

First version of redesigned RNG that provides one instance per thread rather than one global instance. refs #568

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 6.6 KB
Line 
1// $Id: random.cc 2768 2012-07-11 10:59:22Z peter $
2
3/*
4  Copyright (C) 2005, 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
5  Copyright (C) 2009, 2011 Peter Johansson
6
7  This file is part of the yat library, http://dev.thep.lu.se/yat
8
9  The yat library is free software; you can redistribute it and/or
10  modify it under the terms of the GNU General Public License as
11  published by the Free Software Foundation; either version 3 of the
12  License, or (at your option) any later version.
13
14  The yat library is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17  General Public License for more details.
18
19  You should have received a copy of the GNU General Public License
20  along with yat. If not, see <http://www.gnu.org/licenses/>.
21*/
22
23#include "random.h"
24#include "yat/statistics/Histogram.h"
25#include "yat/utility/Exception.h"
26
27#include <cassert>
28#include <cstring>
29#include <fstream>
30#include <sstream>
31
32namespace theplu {
33namespace yat {
34namespace random {
35
36  boost::thread_specific_ptr<RNG> RNG::instance_;
37
38  RNG::RNG(void)
39  {
40      // support rng/seed changes through environment vars
41    if (!gsl_rng_env_setup())
42      throw utility::GSL_error("RNG::RNG unknown generator");
43    if (!(rng_=gsl_rng_alloc(gsl_rng_default)))
44      throw utility::GSL_error("RNG::RNG failed to allocate memory");
45  }
46
47
48  RNG::~RNG(void)
49  {
50    gsl_rng_free(rng_);
51    rng_=NULL;
52  }
53
54
55  RNG* RNG::instance(void)
56  {
57    if (instance_.get()==NULL)
58      instance_.reset(new RNG);
59    return instance_.get();
60  }
61
62
63  unsigned long RNG::max(void) const
64  {
65    return gsl_rng_max(rng_);
66  }
67
68
69  unsigned long RNG::min(void) const
70  {
71    return gsl_rng_min(rng_);
72  }
73
74
75  std::string RNG::name(void) const
76  {
77    return gsl_rng_name(rng_);
78  }
79
80
81  const gsl_rng* RNG::rng(void) const
82  {
83    return rng_;
84  }
85
86
87  void RNG::seed(unsigned long s) const
88  {
89    gsl_rng_set(rng_,s);
90  }
91
92
93  unsigned long RNG::seed_from_devurandom(void)
94  {
95    unsigned char ulongsize=sizeof(unsigned long);
96    char* buffer=new char[ulongsize];
97    std::ifstream is("/dev/urandom", std::ios::binary);
98    is.read(buffer,ulongsize);
99    is.close();
100    unsigned long s=0;
101    memcpy(&s, buffer, ulongsize);
102    delete[] buffer;
103    seed(s);
104    return s;
105  }
106
107
108  int RNG::set_state(const RNG_state& state)
109  {
110    return gsl_rng_memcpy(rng_, state.rng());
111  }
112
113  // --------------------- RNG_state ----------------------------------
114
115  RNG_state::RNG_state(const RNG* rng)
116  {
117    clone(*rng->rng());
118  }
119 
120
121  RNG_state::RNG_state(const RNG_state& state)
122  {
123    clone(*state.rng());
124  }
125 
126
127  RNG_state::~RNG_state(void)
128  {
129    gsl_rng_free(rng_);
130    rng_=NULL;
131  }
132
133  const gsl_rng* RNG_state::rng(void) const
134  {
135    return rng_;
136  }
137
138
139  void RNG_state::clone(const gsl_rng& rng)
140  {
141    assert(rng_!=&rng);
142    if (!(rng_ = gsl_rng_clone(&rng)))
143      throw utility::GSL_error("RNG_state::RNG_state failed to allocate memory");
144  }
145
146  RNG_state& RNG_state::operator=(const RNG_state& rhs)
147  {
148    if (this != &rhs) {
149      gsl_rng_free(rng_);
150      clone(*rhs.rng());
151    }
152    return *this;
153  }
154
155  // --------------------- Discrete distribtuions ---------------------
156
157  Discrete::Discrete(void)
158    : rng_(RNG::instance())
159  {
160  }
161
162 
163  Discrete::~Discrete(void)
164  {
165  }
166
167
168  void Discrete::seed(unsigned long s) const
169  {
170    rng_->seed(s);
171  }
172
173
174  unsigned long Discrete::seed_from_devurandom(void) 
175  { 
176    return rng_->seed_from_devurandom(); 
177  }
178
179
180  DiscreteGeneral::DiscreteGeneral(const statistics::Histogram& hist)
181    : gen_(NULL)
182  {
183    p_.reserve(hist.nof_bins());
184    for (size_t i=0; i<hist.nof_bins(); i++) 
185      p_.push_back(hist[i]);
186    preproc();
187  }
188
189
190  DiscreteGeneral::DiscreteGeneral(const DiscreteGeneral& other)
191    : Discrete(other), gen_(NULL), p_(other.p_)
192  {
193    preproc();
194  }
195
196
197  DiscreteGeneral::~DiscreteGeneral(void)
198  {
199    free();
200  }
201
202
203  void DiscreteGeneral::free(void)
204  {
205    if (gen_)
206      gsl_ran_discrete_free( gen_ );
207    gen_ = NULL;
208  }
209
210
211  void DiscreteGeneral::preproc(void)
212  {
213    assert(!gen_);
214    assert(p_.size());
215    gen_ = gsl_ran_discrete_preproc( p_.size(), &p_.front() );
216    if (!gen_)
217      throw utility::GSL_error("DiscreteGeneral failed to setup generator.");
218  }
219
220
221  DiscreteGeneral& DiscreteGeneral::operator=(const DiscreteGeneral& rhs)
222  {
223    free();
224    p_ = rhs.p_;
225    preproc();
226    return *this;
227  }
228
229
230  unsigned long DiscreteGeneral::operator()(void) const
231  {
232    return gsl_ran_discrete(rng_->rng(), gen_);
233  }
234
235
236  DiscreteUniform::DiscreteUniform(unsigned long n)
237    : range_(n)
238  {
239    if (range_>rng_->max()) {
240      std::stringstream ss("DiscreteUniform::DiscreteUniform: ");
241      ss << n << " is too large for RNG " << rng_->name();
242      throw utility::GSL_error(ss.str());
243    }
244  }
245
246
247  unsigned long DiscreteUniform::operator()(void) const
248  {
249    return (range_ ?
250            gsl_rng_uniform_int(rng_->rng(), range_) : gsl_rng_get(rng_->rng()));
251  }
252
253
254  unsigned long DiscreteUniform::operator()(unsigned long n) const
255  {
256    // making sure that n is not larger than the range of the
257    // underlying RNG
258    if (n>rng_->max()) {
259      std::stringstream ss("DiscreteUniform::operator(unsigned long): ");
260      ss << n << " is too large for RNG " << rng_->name();
261      throw utility::GSL_error(ss.str());
262    }
263    return gsl_rng_uniform_int(rng_->rng(),n);
264  }
265
266
267  Poisson::Poisson(const double m)
268    : m_(m)
269  {
270  }
271
272  unsigned long Poisson::operator()(void) const
273  {
274    return gsl_ran_poisson(rng_->rng(), m_);
275  }
276
277
278  unsigned long Poisson::operator()(const double m) const
279  {
280    return gsl_ran_poisson(rng_->rng(), m);
281  }
282
283  // --------------------- Continuous distribtuions ---------------------
284
285  Continuous::Continuous(void)
286    : rng_(RNG::instance())
287  {
288  }
289
290
291  Continuous::~Continuous(void)
292  {
293  }
294
295
296  void Continuous::seed(unsigned long s) const
297  {
298    rng_->seed(s);
299  }
300
301
302  unsigned long Continuous::seed_from_devurandom(void)
303  {
304    return rng_->seed_from_devurandom();
305  }
306
307
308  ContinuousGeneral::ContinuousGeneral(const statistics::Histogram& hist)
309    : discrete_(DiscreteGeneral(hist)), hist_(hist)
310  {
311  }
312
313
314  double ContinuousGeneral::operator()(void) const
315  {
316    return hist_.observation_value(discrete_())+(u_()-0.5)*hist_.spacing();
317  }
318
319  double ContinuousUniform::operator()(void) const
320  {
321    return gsl_rng_uniform(rng_->rng());
322  }
323
324
325  Exponential::Exponential(const double m)
326    : m_(m)
327  {
328  }
329
330
331  double Exponential::operator()(void) const
332  {
333    return gsl_ran_exponential(rng_->rng(), m_);
334  }
335
336
337  double Exponential::operator()(const double m) const
338  {
339    return gsl_ran_exponential(rng_->rng(), m);
340  }
341
342
343  Gaussian::Gaussian(const double s, const double m)
344    : m_(m), s_(s)
345  {
346  }
347
348
349  double Gaussian::operator()(void) const
350  {
351    return gsl_ran_gaussian(rng_->rng(), s_)+m_;
352  }
353
354
355  double Gaussian::operator()(const double s) const
356  {
357    return gsl_ran_gaussian(rng_->rng(), s);
358  }
359
360
361  double Gaussian::operator()(const double s, const double m) const
362  {
363    return gsl_ran_gaussian(rng_->rng(), s)+m;
364  }
365
366}}} // of namespace random, yat, and theplu
Note: See TracBrowser for help on using the repository browser.