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

Last change on this file since 3579 was 3579, checked in by Peter, 5 years ago

merge release 0.14 into trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 9.7 KB
Line 
1// $Id: random.cc 3579 2017-01-16 03:54:43Z peter $
2
3/*
4  Copyright (C) 2005, 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
5  Copyright (C) 2009, 2011, 2012, 2013, 2015, 2016, 2017 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 <config.h>
24
25#include "random.h"
26#include "yat/statistics/Histogram.h"
27#include "yat/utility/Exception.h"
28
29#include <boost/thread/locks.hpp>
30
31#include <cassert>
32#include <cstring>
33#include <fstream>
34#include <sstream>
35
36namespace theplu {
37namespace yat {
38namespace random {
39
40  RNG* RNG::instance_=NULL;
41  boost::mutex RNG::init_mutex_;
42
43  RNG::RNG(void)
44    : rng_(gsl_rng_free)
45  {
46      // support rng/seed changes through environment vars
47    if (!gsl_rng_env_setup())
48      throw utility::GSL_error("RNG::RNG unknown generator");
49    seed_ = gsl_rng_default_seed;
50    // let's allocate already here, just to behave as yat 0.8
51    rng_alloc();
52  }
53
54
55  RNG::~RNG(void)
56  {
57  }
58
59
60  RNG* RNG::instance(void)
61  {
62    if (instance_==NULL) {
63      boost::unique_lock<boost::mutex> lock(init_mutex_);
64      if (instance_==NULL)
65        instance_ = new RNG;
66    }
67    return instance_;
68  }
69
70
71  unsigned long RNG::max(void) const
72  {
73    return gsl_rng_max(rng());
74  }
75
76
77  unsigned long RNG::min(void) const
78  {
79    return gsl_rng_min(rng());
80  }
81
82
83  std::string RNG::name(void) const
84  {
85    return gsl_rng_name(rng());
86  }
87
88
89  const gsl_rng* RNG::rng(void) const
90  {
91    if (rng_.get()==NULL)
92      rng_alloc();
93    return rng_.get();
94  }
95
96
97  void RNG::rng_alloc(void) const
98  {
99    assert(rng_.get()==NULL);
100    gsl_rng* rng = gsl_rng_alloc(gsl_rng_default);
101    if (!rng)
102      throw utility::GSL_error("RNG failed to allocate memory");
103    boost::unique_lock<boost::mutex> lock(mutex_);
104    gsl_rng_set(rng, seed_);
105    // bump seed to avoid subsequent gsl_rng to be identical
106    ++seed_;
107    // rng_ owns rng and takes care of deallocation
108    rng_.reset(rng);
109  } // lock is released here
110
111
112  void RNG::seed(unsigned long s) const
113  {
114    boost::unique_lock<boost::mutex> lock(mutex_);
115    gsl_rng_set(rng(),s);
116    seed_ = s+1;
117  } // lock is released here
118
119
120  unsigned long RNG::seed_from_devurandom(void)
121  {
122    unsigned char ulongsize=sizeof(unsigned long);
123    char* buffer=new char[ulongsize];
124    std::ifstream is("/dev/urandom", std::ios::binary);
125    is.read(buffer,ulongsize);
126    is.close();
127    unsigned long s=0;
128    memcpy(&s, buffer, ulongsize);
129    delete[] buffer;
130    seed(s);
131    return s;
132  }
133
134
135  int RNG::set_state(const RNG_state& state)
136  {
137    if (rng_.get()==NULL)
138      rng_alloc();
139    if (gsl_rng_memcpy(rng_.get(), state.rng()))
140      throw utility::GSL_error("yat::random::RNG::set_state failed");
141    return 0;
142  }
143
144  // --------------------- RNG_state ----------------------------------
145
146  RNG_state::RNG_state(const RNG* rng)
147  {
148    clone(*rng->rng());
149  }
150
151
152  RNG_state::RNG_state(const RNG_state& state)
153  {
154    clone(*state.rng());
155  }
156
157
158  RNG_state::~RNG_state(void)
159  {
160    gsl_rng_free(rng_);
161    rng_=NULL;
162  }
163
164  const gsl_rng* RNG_state::rng(void) const
165  {
166    return rng_;
167  }
168
169
170  void RNG_state::clone(const gsl_rng& rng)
171  {
172    assert(rng_!=&rng);
173    if (!(rng_ = gsl_rng_clone(&rng)))
174      throw utility::GSL_error("RNG_state::clone failed to allocate memory");
175  }
176
177  RNG_state& RNG_state::operator=(const RNG_state& rhs)
178  {
179    if (this != &rhs) {
180      gsl_rng_free(rng_);
181      clone(*rhs.rng());
182    }
183    return *this;
184  }
185
186  // --------------------- Discrete distribtuions ---------------------
187
188  Discrete::Discrete(void)
189    : rng_(RNG::instance())
190  {
191  }
192
193
194  Discrete::~Discrete(void)
195  {
196  }
197
198
199  void Discrete::seed(unsigned long s) const
200  {
201    rng_->seed(s);
202  }
203
204
205  unsigned long Discrete::seed_from_devurandom(void)
206  {
207    return rng_->seed_from_devurandom();
208  }
209
210
211  Binomial::Binomial(double p, unsigned int n)
212    : Discrete(), p_(p), n_(n)
213  {
214  }
215
216
217  unsigned long Binomial::operator()(void) const
218  {
219    return gsl_ran_binomial(rng_->rng(), p_, n_);
220  }
221
222
223  DiscreteGeneral::DiscreteGeneral(const statistics::Histogram& hist)
224    : gen_(NULL)
225  {
226    p_.reserve(hist.nof_bins());
227    for (size_t i=0; i<hist.nof_bins(); i++)
228      p_.push_back(hist[i]);
229    preproc();
230  }
231
232
233  DiscreteGeneral::DiscreteGeneral(const DiscreteGeneral& other)
234    : Discrete(other), gen_(NULL), p_(other.p_)
235  {
236    preproc();
237  }
238
239
240  DiscreteGeneral::~DiscreteGeneral(void)
241  {
242    free();
243  }
244
245
246  void DiscreteGeneral::free(void)
247  {
248    if (gen_)
249      gsl_ran_discrete_free( gen_ );
250    gen_ = NULL;
251  }
252
253
254  void DiscreteGeneral::preproc(void)
255  {
256    assert(!gen_);
257    assert(p_.size());
258    gen_ = gsl_ran_discrete_preproc( p_.size(), &p_.front() );
259    if (!gen_)
260      throw utility::GSL_error("DiscreteGeneral failed to setup generator.");
261  }
262
263
264  DiscreteGeneral& DiscreteGeneral::operator=(const DiscreteGeneral& rhs)
265  {
266    free();
267    p_ = rhs.p_;
268    preproc();
269    return *this;
270  }
271
272
273  unsigned long DiscreteGeneral::operator()(void) const
274  {
275    return gsl_ran_discrete(rng_->rng(), gen_);
276  }
277
278
279  DiscreteUniform::DiscreteUniform(unsigned long n)
280    : range_(n)
281  {
282    if (range_>rng_->max()) {
283      std::ostringstream ss;
284      ss << "DiscreteUniform::DiscreteUniform: ";
285      ss << n << " is too large for RNG " << rng_->name();
286      ss << "; maximal argument is " << rng_->max();
287      throw utility::GSL_error(ss.str());
288    }
289  }
290
291
292  unsigned long DiscreteUniform::operator()(void) const
293  {
294    return (range_ ?
295            gsl_rng_uniform_int(rng_->rng(),range_) : gsl_rng_get(rng_->rng()));
296  }
297
298
299  unsigned long DiscreteUniform::operator()(unsigned long n) const
300  {
301    // making sure that n is not larger than the range of the
302    // underlying RNG
303    if (n>rng_->max()) {
304      std::ostringstream ss;
305      ss << "DiscreteUniform::operator(unsigned long): ";
306      ss << n << " is too large for RNG " << rng_->name();
307      ss << "; maximal argument is " << rng_->max();
308      throw utility::GSL_error(ss.str());
309    }
310    return gsl_rng_uniform_int(rng_->rng(),n);
311  }
312
313
314  Geometric::Geometric(double p)
315    : p_(p)
316  {}
317
318
319  unsigned long int Geometric::operator()(void) const
320  {
321    return gsl_ran_geometric (rng_->rng(), p_); 
322  }
323
324
325  unsigned long int Geometric::operator()(double p) const
326  {
327    return gsl_ran_geometric (rng_->rng(), p);
328  }
329
330
331  HyperGeometric::HyperGeometric(void)
332  {}
333
334
335  HyperGeometric::HyperGeometric(unsigned int n1, unsigned int n2,
336                                 unsigned int t)
337    : n1_(n1), n2_(n2), t_(t)
338  {}
339
340
341  unsigned long int HyperGeometric::operator()(void) const
342  {
343    return (*this)(n1_, n2_, t_);
344  }
345
346
347  unsigned long int HyperGeometric::operator()(unsigned int n1,
348                                               unsigned int n2,
349                                               unsigned int t) const
350  {
351    return gsl_ran_hypergeometric(rng_->rng(), n1, n2, t);
352  }
353
354
355  NegativeHyperGeometric::NegativeHyperGeometric(void)
356  {}
357
358
359  NegativeHyperGeometric::NegativeHyperGeometric(unsigned int n1,
360                                                 unsigned int n2, unsigned int t)
361    : n1_(n1), n2_(n2), t_(t)
362  {}
363
364
365  unsigned long int NegativeHyperGeometric::operator()(void) const
366  {
367    return (*this)(n1_, n2_, t_);
368  }
369
370
371  unsigned long int NegativeHyperGeometric::operator()(unsigned int n1,
372                                                       unsigned int n2,
373                                                       unsigned int t) const
374  {
375    assert(t <= n2);
376
377    // NHG can be described as an array with n1 true and n2 false, and
378    // NHG(n1, n2, t) is the number of true left of the t:th false. By
379    // symmetry number of true right of the t:th false is NHG(n1, n2,
380    // n2-t+1) since t:th false counting from left is the (n2-t+1):th
381    // false counting from right.
382
383    // When t is larger than midpoint (2*t > n2+1) we use this
384    // symmetry to speed things up.
385    if (t > (n2+1)/2)
386      return n1 - (*this)(n1, n2, n2-t+1);
387
388    ContinuousUniform uniform;
389    unsigned long int k = 0;
390    while (t) {
391      assert(n1 + n2);
392      double x = uniform();
393      if (x * (n1+n2) < n1) {
394        --n1;
395        ++k;
396        if (!n1)
397          return k;
398      }
399      else {
400        --t;
401        --n2;
402      }
403    }
404    return k;
405  }
406
407
408  Poisson::Poisson(const double m)
409    : m_(m)
410  {
411  }
412
413  unsigned long Poisson::operator()(void) const
414  {
415    return gsl_ran_poisson(rng_->rng(), m_);
416  }
417
418
419  unsigned long Poisson::operator()(const double m) const
420  {
421    return gsl_ran_poisson(rng_->rng(), m);
422  }
423
424  // --------------------- Continuous distribtuions ---------------------
425
426  Continuous::Continuous(void)
427    : rng_(RNG::instance())
428  {
429  }
430
431
432  Continuous::~Continuous(void)
433  {
434  }
435
436
437  void Continuous::seed(unsigned long s) const
438  {
439    rng_->seed(s);
440  }
441
442
443  unsigned long Continuous::seed_from_devurandom(void)
444  {
445    return rng_->seed_from_devurandom();
446  }
447
448
449  ContinuousGeneral::ContinuousGeneral(const statistics::Histogram& hist)
450    : discrete_(DiscreteGeneral(hist)), hist_(hist)
451  {
452  }
453
454
455  double ContinuousGeneral::operator()(void) const
456  {
457    return hist_.observation_value(discrete_())+(u_()-0.5)*hist_.spacing();
458  }
459
460  double ContinuousUniform::operator()(void) const
461  {
462    return gsl_rng_uniform(rng_->rng());
463  }
464
465
466  Exponential::Exponential(const double m)
467    : m_(m)
468  {
469  }
470
471
472  double Exponential::operator()(void) const
473  {
474    return gsl_ran_exponential(rng_->rng(), m_);
475  }
476
477
478  double Exponential::operator()(const double m) const
479  {
480    return gsl_ran_exponential(rng_->rng(), m);
481  }
482
483
484  Gaussian::Gaussian(const double s, const double m)
485    : m_(m), s_(s)
486  {
487  }
488
489
490  double Gaussian::operator()(void) const
491  {
492    return gsl_ran_gaussian(rng_->rng(), s_)+m_;
493  }
494
495
496  double Gaussian::operator()(const double s) const
497  {
498    return gsl_ran_gaussian(rng_->rng(), s);
499  }
500
501
502  double Gaussian::operator()(const double s, const double m) const
503  {
504    return gsl_ran_gaussian(rng_->rng(), s)+m;
505  }
506
507}}} // of namespace random, yat, and theplu
Note: See TracBrowser for help on using the repository browser.