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

Last change on this file since 2588 was 2588, checked in by Peter, 11 years ago

prefer using standard functions

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 6.6 KB
Line 
1// $Id: random.cc 2588 2011-10-25 21:38:55Z 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  {
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_)
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    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");
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.