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

Last change on this file since 2804 was 2804, checked in by Peter, 9 years ago

RNG::set_state: prefer throw on error.

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