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

Last change on this file since 718 was 718, checked in by Jari Häkkinen, 16 years ago

Addresses #170.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 4.9 KB
Line 
1// $Id: random.cc 718 2006-12-26 09:56:26Z jari $
2
3/*
4  Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson
5
6  This file is part of the yat library, http://lev.thep.lu.se/trac/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 2 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 this program; if not, write to the Free Software
20  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
21  02111-1307, USA.
22*/
23
24#include "random.h"
25#include "yat/statistics/Histogram.h"
26
27#include <fstream>
28#include <iostream>
29
30namespace theplu {
31namespace yat {
32namespace random {
33
34  RNG* RNG::instance_ = NULL;
35
36  RNG::RNG(void)
37  {
38    gsl_rng_env_setup();  // support rng/seed changes through environment vars
39    rng_ = gsl_rng_alloc(gsl_rng_default);  // Memory is allocated here!
40  }
41
42
43  RNG::~RNG(void)
44  {
45    gsl_rng_free(rng_);
46    rng_=NULL;
47    delete instance_;
48  }
49
50 
51  u_long RNG::max(void) const
52  {
53    return gsl_rng_max(rng_);
54  }
55
56
57  u_long RNG::min(void) const
58  {
59    return gsl_rng_min(rng_);
60  }
61
62
63  std::string RNG::name(void) const
64  {
65    return gsl_rng_name(rng_);
66  }
67
68
69  const gsl_rng* RNG::rng(void) const
70  {
71    return rng_;
72  }
73
74
75  void RNG::seed(u_long s) const
76  {
77    gsl_rng_set(rng_,s);
78  }
79
80
81  u_long RNG::seed_from_devurandom(void)
82  {
83    u_char ulongsize=sizeof(u_long);
84    char* buffer=new char[ulongsize];
85    std::ifstream is("/dev/urandom");
86    is.read(buffer,ulongsize);
87    is.close();
88    u_long s=0;
89    for (u_int i=0; i<ulongsize; i++) {
90      u_char ch=buffer[i];
91      s+=ch<<((ulongsize-1-i)*8);
92    }
93    seed(s);
94    return s;
95  }
96
97
98  int RNG::set_state(const RNG_state& state)
99  {
100    return gsl_rng_memcpy(rng_, state.rng());
101  }
102
103  // --------------------- RNG_state ----------------------------------
104
105  RNG_state::RNG_state(const RNG* rng)
106  {
107    rng_ = gsl_rng_clone(rng->rng()); 
108  }
109 
110
111  RNG_state::~RNG_state(void)
112  {
113    gsl_rng_free(rng_);
114    rng_=NULL;
115  }
116
117  const gsl_rng* RNG_state::rng(void) const
118  {
119    return rng_;
120  }
121
122  // --------------------- Discrete distribtuions ---------------------
123
124  Discrete::Discrete(void)
125    : rng_(RNG::instance())
126  {
127  }
128
129 
130  Discrete::~Discrete(void)
131  {
132  }
133
134
135  void Discrete::seed(u_long s) const
136  {
137    rng_->seed(s);
138  }
139
140
141  DiscreteGeneral::DiscreteGeneral(const statistics::Histogram& hist)
142  {
143    double* p = new double[ hist.nof_bins() ];
144    for (size_t i=0; i<hist.nof_bins(); i++) 
145      p[ i ] = hist[i];
146    gen_ = gsl_ran_discrete_preproc( hist.nof_bins(), p );
147    delete p;
148  }
149
150
151  DiscreteGeneral::~DiscreteGeneral(void)
152  {
153    if (gen_)
154      gsl_ran_discrete_free( gen_ );
155    gen_ = NULL;
156  }
157
158
159  u_long DiscreteGeneral::operator()(void) const
160  {
161    return gsl_ran_discrete(rng_->rng(), gen_);
162  }
163
164
165  DiscreteUniform::DiscreteUniform(void)
166    : range_(rng_->max())
167  {
168  }
169
170  DiscreteUniform::DiscreteUniform(const u_long n)
171    : range_(n)
172  {
173    if (range_>rng_->max())
174      range_=rng_->max();
175  }
176
177
178  u_long DiscreteUniform::operator()(void) const
179  {
180    return gsl_rng_uniform_int(rng_->rng(), range_);
181  }
182
183
184  u_long DiscreteUniform::operator()(const u_long n) const
185  {
186    return gsl_rng_uniform_int(rng_->rng(), n);
187  }
188
189
190  Poisson::Poisson(const double m)
191    : m_(m)
192  {
193  }
194
195  u_long Poisson::operator()(void) const
196  {
197    return gsl_ran_poisson(rng_->rng(), m_);
198  }
199
200
201  u_long Poisson::operator()(const double m) const
202  {
203    return gsl_ran_poisson(rng_->rng(), m);
204  }
205
206  // --------------------- Continuous distribtuions ---------------------
207
208  Continuous::Continuous(void)
209    : rng_(RNG::instance())
210  {
211  }
212
213
214  Continuous::~Continuous(void)
215  {
216  }
217
218
219  void Continuous::seed(u_long s) const
220  {
221    rng_->seed(s);
222  }
223
224
225  ContinuousGeneral::ContinuousGeneral(const statistics::Histogram& hist)
226    : discrete_(DiscreteGeneral(hist)), hist_(hist)
227  {
228  }
229
230
231  double ContinuousGeneral::operator()(void) const
232  {
233    return hist_.observation_value(discrete_())+(u_()-0.5)*hist_.spacing();
234  }
235
236  double ContinuousUniform::operator()(void) const
237  {
238    return gsl_rng_uniform(rng_->rng());
239  }
240
241
242  Exponential::Exponential(const double m)
243    : m_(m)
244  {
245  }
246
247
248  double Exponential::operator()(void) const
249  {
250    return gsl_ran_exponential(rng_->rng(), m_);
251  }
252
253
254  double Exponential::operator()(const double m) const
255  {
256    return gsl_ran_exponential(rng_->rng(), m);
257  }
258
259
260  Gaussian::Gaussian(const double s, const double m)
261    : m_(m), s_(s)
262  {
263  }
264
265
266  double Gaussian::operator()(void) const
267  {
268    return gsl_ran_gaussian(rng_->rng(), s_)+m_;
269  }
270
271
272  double Gaussian::operator()(const double s) const
273  {
274    return gsl_ran_gaussian(rng_->rng(), s);
275  }
276
277
278  double Gaussian::operator()(const double s, const double m) const
279  {
280    return gsl_ran_gaussian(rng_->rng(), s)+m;
281  }
282
283}}} // of namespace random, yat, and theplu
Note: See TracBrowser for help on using the repository browser.