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

Last change on this file since 4010 was 4010, checked in by Peter, 12 months ago

Implement move constructo and move-assignment operator for
DiscreteGeneral? (closes #964). Using swap idiom also for
copy-assignment as it gives stronger exception safety, i.e., if
preproc throws, lhs is not modified (whereas before ::p_ was).

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