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

Last change on this file since 3591 was 3591, checked in by Peter, 6 years ago

fixes #877

Generalize YAT_CXX_RVALUE into YAT_CXX_TRY_CXX11 and use that to
implement new macro YAT_CXX_ATOMIC. Use this macro in configure.ac and
use std::atomic<> conditoionally in RNG singleton.

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