source: branches/0.9-stable/yat/random/random.cc @ 2833

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

fixes #720

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