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

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

change back to RNG being a singleton, which now stores gsl_rng* in a thread_specific_ptr<gsl_rng>. refs #568

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 6.8 KB
Line 
1// $Id: random.cc 2770 2012-07-11 14:44:45Z 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  RNG* RNG::instance_=NULL;
37
38  RNG::RNG(void)
39    : rng_(gsl_rng_free)
40  {
41      // support rng/seed changes through environment vars
42    if (!gsl_rng_env_setup())
43      throw utility::GSL_error("RNG::RNG unknown generator");
44  }
45
46
47  RNG::~RNG(void)
48  {
49  }
50
51
52  RNG* RNG::instance(void)
53  {
54    if (instance_==NULL)
55      instance_ = new RNG;
56    return instance_;
57  }
58
59
60  unsigned long RNG::max(void) const
61  {
62    return gsl_rng_max(rng());
63  }
64
65
66  unsigned long RNG::min(void) const
67  {
68    return gsl_rng_min(rng());
69  }
70
71
72  std::string RNG::name(void) const
73  {
74    return gsl_rng_name(rng());
75  }
76
77
78  const gsl_rng* RNG::rng(void) const
79  {
80    if (rng_.get()==NULL)
81      rng_alloc();
82    return rng_.get();
83  }
84
85
86  void RNG::rng_alloc(void) const
87  {
88    assert(rng_.get()==NULL);
89    gsl_rng* rng = gsl_rng_alloc(gsl_rng_default);
90    if (!rng)
91      throw utility::GSL_error("RNG failed to allocate memory");
92    rng_.reset(rng);
93  }
94
95
96  void RNG::seed(unsigned long s) const
97  {
98    gsl_rng_set(rng(),s);
99  }
100
101
102  unsigned long RNG::seed_from_devurandom(void)
103  {
104    unsigned char ulongsize=sizeof(unsigned long);
105    char* buffer=new char[ulongsize];
106    std::ifstream is("/dev/urandom", std::ios::binary);
107    is.read(buffer,ulongsize);
108    is.close();
109    unsigned long s=0;
110    memcpy(&s, buffer, ulongsize);
111    delete[] buffer;
112    seed(s);
113    return s;
114  }
115
116
117  int RNG::set_state(const RNG_state& state)
118  {
119    if (rng_.get()==NULL)
120      rng_alloc();
121    return gsl_rng_memcpy(rng_.get(), state.rng());
122  }
123
124  // --------------------- RNG_state ----------------------------------
125
126  RNG_state::RNG_state(const RNG* rng)
127  {
128    clone(*rng->rng());
129  }
130
131
132  RNG_state::RNG_state(const RNG_state& state)
133  {
134    clone(*state.rng());
135  }
136
137
138  RNG_state::~RNG_state(void)
139  {
140    gsl_rng_free(rng_);
141    rng_=NULL;
142  }
143
144  const gsl_rng* RNG_state::rng(void) const
145  {
146    return rng_;
147  }
148
149
150  void RNG_state::clone(const gsl_rng& rng)
151  {
152    assert(rng_!=&rng);
153    if (!(rng_ = gsl_rng_clone(&rng)))
154      throw utility::GSL_error("RNG_state::RNG_state failed to allocate memory");
155  }
156
157  RNG_state& RNG_state::operator=(const RNG_state& rhs)
158  {
159    if (this != &rhs) {
160      gsl_rng_free(rng_);
161      clone(*rhs.rng());
162    }
163    return *this;
164  }
165
166  // --------------------- Discrete distribtuions ---------------------
167
168  Discrete::Discrete(void)
169    : rng_(RNG::instance())
170  {
171  }
172
173 
174  Discrete::~Discrete(void)
175  {
176  }
177
178
179  void Discrete::seed(unsigned long s) const
180  {
181    rng_->seed(s);
182  }
183
184
185  unsigned long Discrete::seed_from_devurandom(void) 
186  { 
187    return rng_->seed_from_devurandom(); 
188  }
189
190
191  DiscreteGeneral::DiscreteGeneral(const statistics::Histogram& hist)
192    : gen_(NULL)
193  {
194    p_.reserve(hist.nof_bins());
195    for (size_t i=0; i<hist.nof_bins(); i++) 
196      p_.push_back(hist[i]);
197    preproc();
198  }
199
200
201  DiscreteGeneral::DiscreteGeneral(const DiscreteGeneral& other)
202    : Discrete(other), gen_(NULL), p_(other.p_)
203  {
204    preproc();
205  }
206
207
208  DiscreteGeneral::~DiscreteGeneral(void)
209  {
210    free();
211  }
212
213
214  void DiscreteGeneral::free(void)
215  {
216    if (gen_)
217      gsl_ran_discrete_free( gen_ );
218    gen_ = NULL;
219  }
220
221
222  void DiscreteGeneral::preproc(void)
223  {
224    assert(!gen_);
225    assert(p_.size());
226    gen_ = gsl_ran_discrete_preproc( p_.size(), &p_.front() );
227    if (!gen_)
228      throw utility::GSL_error("DiscreteGeneral failed to setup generator.");
229  }
230
231
232  DiscreteGeneral& DiscreteGeneral::operator=(const DiscreteGeneral& rhs)
233  {
234    free();
235    p_ = rhs.p_;
236    preproc();
237    return *this;
238  }
239
240
241  unsigned long DiscreteGeneral::operator()(void) const
242  {
243    return gsl_ran_discrete(rng_->rng(), gen_);
244  }
245
246
247  DiscreteUniform::DiscreteUniform(unsigned long n)
248    : range_(n)
249  {
250    if (range_>rng_->max()) {
251      std::stringstream ss("DiscreteUniform::DiscreteUniform: ");
252      ss << n << " is too large for RNG " << rng_->name();
253      throw utility::GSL_error(ss.str());
254    }
255  }
256
257
258  unsigned long DiscreteUniform::operator()(void) const
259  {
260    return (range_ ?
261            gsl_rng_uniform_int(rng_->rng(), range_) : gsl_rng_get(rng_->rng()));
262  }
263
264
265  unsigned long DiscreteUniform::operator()(unsigned long n) const
266  {
267    // making sure that n is not larger than the range of the
268    // underlying RNG
269    if (n>rng_->max()) {
270      std::stringstream ss("DiscreteUniform::operator(unsigned long): ");
271      ss << n << " is too large for RNG " << rng_->name();
272      throw utility::GSL_error(ss.str());
273    }
274    return gsl_rng_uniform_int(rng_->rng(),n);
275  }
276
277
278  Poisson::Poisson(const double m)
279    : m_(m)
280  {
281  }
282
283  unsigned long Poisson::operator()(void) const
284  {
285    return gsl_ran_poisson(rng_->rng(), m_);
286  }
287
288
289  unsigned long Poisson::operator()(const double m) const
290  {
291    return gsl_ran_poisson(rng_->rng(), m);
292  }
293
294  // --------------------- Continuous distribtuions ---------------------
295
296  Continuous::Continuous(void)
297    : rng_(RNG::instance())
298  {
299  }
300
301
302  Continuous::~Continuous(void)
303  {
304  }
305
306
307  void Continuous::seed(unsigned long s) const
308  {
309    rng_->seed(s);
310  }
311
312
313  unsigned long Continuous::seed_from_devurandom(void)
314  {
315    return rng_->seed_from_devurandom();
316  }
317
318
319  ContinuousGeneral::ContinuousGeneral(const statistics::Histogram& hist)
320    : discrete_(DiscreteGeneral(hist)), hist_(hist)
321  {
322  }
323
324
325  double ContinuousGeneral::operator()(void) const
326  {
327    return hist_.observation_value(discrete_())+(u_()-0.5)*hist_.spacing();
328  }
329
330  double ContinuousUniform::operator()(void) const
331  {
332    return gsl_rng_uniform(rng_->rng());
333  }
334
335
336  Exponential::Exponential(const double m)
337    : m_(m)
338  {
339  }
340
341
342  double Exponential::operator()(void) const
343  {
344    return gsl_ran_exponential(rng_->rng(), m_);
345  }
346
347
348  double Exponential::operator()(const double m) const
349  {
350    return gsl_ran_exponential(rng_->rng(), m);
351  }
352
353
354  Gaussian::Gaussian(const double s, const double m)
355    : m_(m), s_(s)
356  {
357  }
358
359
360  double Gaussian::operator()(void) const
361  {
362    return gsl_ran_gaussian(rng_->rng(), s_)+m_;
363  }
364
365
366  double Gaussian::operator()(const double s) const
367  {
368    return gsl_ran_gaussian(rng_->rng(), s);
369  }
370
371
372  double Gaussian::operator()(const double s, const double m) const
373  {
374    return gsl_ran_gaussian(rng_->rng(), s)+m;
375  }
376
377}}} // of namespace random, yat, and theplu
Note: See TracBrowser for help on using the repository browser.