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

Last change on this file since 2090 was 2090, checked in by Peter, 12 years ago

initialize base class explicitely in copy constructor. doesn't really matter here but is just a matter of style.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 6.5 KB
Line 
1// $Id: random.cc 2090 2009-10-21 21:53:33Z peter $
2
3/*
4  Copyright (C) 2005, 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
5
6  This file is part of the yat library, http://dev.thep.lu.se/yat
7
8  The yat library is free software; you can redistribute it and/or
9  modify it under the terms of the GNU General Public License as
10  published by the Free Software Foundation; either version 3 of the
11  License, or (at your option) any later version.
12
13  The yat library is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16  General Public License for more details.
17
18  You should have received a copy of the GNU General Public License
19  along with yat. If not, see <http://www.gnu.org/licenses/>.
20*/
21
22#include "random.h"
23#include "yat/statistics/Histogram.h"
24#include "yat/utility/Exception.h"
25
26#include <cassert>
27#include <fstream>
28#include <sstream>
29
30namespace theplu {
31namespace yat {
32namespace random {
33
34  RNG* RNG::instance_ = NULL;
35
36  RNG::RNG(void)
37  {
38      // support rng/seed changes through environment vars
39    if (!gsl_rng_env_setup())
40      throw utility::GSL_error("RNG::RNG unknown generator");
41    if (!(rng_=gsl_rng_alloc(gsl_rng_default)))
42      throw utility::GSL_error("RNG::RNG failed to allocate memory");
43  }
44
45
46  RNG::~RNG(void)
47  {
48    gsl_rng_free(rng_);
49    rng_=NULL;
50  }
51
52 
53  RNG* RNG::instance(void)
54  {
55    if (!instance_)
56      instance_=new RNG;
57    return instance_;
58  }
59
60
61  unsigned long RNG::max(void) const
62  {
63    return gsl_rng_max(rng_);
64  }
65
66
67  unsigned long RNG::min(void) const
68  {
69    return gsl_rng_min(rng_);
70  }
71
72
73  std::string RNG::name(void) const
74  {
75    return gsl_rng_name(rng_);
76  }
77
78
79  const gsl_rng* RNG::rng(void) const
80  {
81    return rng_;
82  }
83
84
85  void RNG::seed(unsigned long s) const
86  {
87    gsl_rng_set(rng_,s);
88  }
89
90
91  unsigned long RNG::seed_from_devurandom(void)
92  {
93    unsigned char ulongsize=sizeof(unsigned long);
94    char* buffer=new char[ulongsize];
95    std::ifstream is("/dev/urandom");
96    is.read(buffer,ulongsize);
97    is.close();
98    unsigned long s=0;
99    for (unsigned int i=0; i<ulongsize; i++) {
100      unsigned char ch=buffer[i];
101      s+=ch<<((ulongsize-1-i)*8);
102    }
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  ContinuousGeneral::ContinuousGeneral(const statistics::Histogram& hist)
303    : discrete_(DiscreteGeneral(hist)), hist_(hist)
304  {
305  }
306
307
308  double ContinuousGeneral::operator()(void) const
309  {
310    return hist_.observation_value(discrete_())+(u_()-0.5)*hist_.spacing();
311  }
312
313  double ContinuousUniform::operator()(void) const
314  {
315    return gsl_rng_uniform(rng_->rng());
316  }
317
318
319  Exponential::Exponential(const double m)
320    : m_(m)
321  {
322  }
323
324
325  double Exponential::operator()(void) const
326  {
327    return gsl_ran_exponential(rng_->rng(), m_);
328  }
329
330
331  double Exponential::operator()(const double m) const
332  {
333    return gsl_ran_exponential(rng_->rng(), m);
334  }
335
336
337  Gaussian::Gaussian(const double s, const double m)
338    : m_(m), s_(s)
339  {
340  }
341
342
343  double Gaussian::operator()(void) const
344  {
345    return gsl_ran_gaussian(rng_->rng(), s_)+m_;
346  }
347
348
349  double Gaussian::operator()(const double s) const
350  {
351    return gsl_ran_gaussian(rng_->rng(), s);
352  }
353
354
355  double Gaussian::operator()(const double s, const double m) const
356  {
357    return gsl_ran_gaussian(rng_->rng(), s)+m;
358  }
359
360}}} // of namespace random, yat, and theplu
Note: See TracBrowser for help on using the repository browser.