yat/random/random.cc

Code
Comments
Other
Rev Date Author Line
362 04 Aug 05 jari 1 // $Id$
362 04 Aug 05 jari 2
561 15 Mar 06 jari 3 /*
2119 12 Dec 09 peter 4   Copyright (C) 2005, 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 5   Copyright (C) 2009, 2011, 2012, 2013, 2015, 2016, 2017, 2018, 2020, 2022, 2023 Peter Johansson
561 15 Mar 06 jari 6
1437 25 Aug 08 peter 7   This file is part of the yat library, http://dev.thep.lu.se/yat
561 15 Mar 06 jari 8
675 10 Oct 06 jari 9   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 10   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 11   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 12   License, or (at your option) any later version.
561 15 Mar 06 jari 13
675 10 Oct 06 jari 14   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 15   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 16   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
561 15 Mar 06 jari 17   General Public License for more details.
561 15 Mar 06 jari 18
561 15 Mar 06 jari 19   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 20   along with yat. If not, see <http://www.gnu.org/licenses/>.
561 15 Mar 06 jari 21 */
561 15 Mar 06 jari 22
2881 18 Nov 12 peter 23 #include <config.h>
2881 18 Nov 12 peter 24
680 11 Oct 06 jari 25 #include "random.h"
675 10 Oct 06 jari 26 #include "yat/statistics/Histogram.h"
738 07 Jan 07 jari 27 #include "yat/utility/Exception.h"
3902 05 May 20 peter 28 #include "yat/utility/utility.h"
371 05 Aug 05 peter 29
1614 05 Nov 08 peter 30 #include <cassert>
3591 20 Jan 17 peter 31 #include <cstddef>
2588 25 Oct 11 peter 32 #include <cstring>
374 05 Aug 05 jari 33 #include <fstream>
3959 30 Jul 20 peter 34 #include <mutex>
738 07 Jan 07 jari 35 #include <sstream>
371 05 Aug 05 peter 36
362 04 Aug 05 jari 37 namespace theplu {
680 11 Oct 06 jari 38 namespace yat {
365 05 Aug 05 peter 39 namespace random {
362 04 Aug 05 jari 40
3591 20 Jan 17 peter 41   std::atomic<RNG*> RNG::instance_ { nullptr };
3959 30 Jul 20 peter 42   std::mutex RNG::init_mutex_;
3959 30 Jul 20 peter 43   thread_local std::unique_ptr<gsl_rng, void(*)(gsl_rng*)> RNG::rng_
3959 30 Jul 20 peter 44   { nullptr, gsl_rng_free};
362 04 Aug 05 jari 45
425 07 Dec 05 jari 46   RNG::RNG(void)
362 04 Aug 05 jari 47   {
738 07 Jan 07 jari 48       // support rng/seed changes through environment vars
738 07 Jan 07 jari 49     if (!gsl_rng_env_setup())
750 17 Feb 07 jari 50       throw utility::GSL_error("RNG::RNG unknown generator");
2800 28 Jul 12 peter 51     seed_ = gsl_rng_default_seed;
2802 29 Jul 12 peter 52     // let's allocate already here, just to behave as yat 0.8
2802 29 Jul 12 peter 53     rng_alloc();
362 04 Aug 05 jari 54   }
362 04 Aug 05 jari 55
362 04 Aug 05 jari 56
362 04 Aug 05 jari 57   RNG::~RNG(void)
362 04 Aug 05 jari 58   {
362 04 Aug 05 jari 59   }
362 04 Aug 05 jari 60
2768 11 Jul 12 peter 61
752 17 Feb 07 jari 62   RNG* RNG::instance(void)
752 17 Feb 07 jari 63   {
3579 16 Jan 17 peter 64     if (instance_==NULL) {
3959 30 Jul 20 peter 65       std::unique_lock<std::mutex> lock(init_mutex_);
3579 16 Jan 17 peter 66       if (instance_==NULL)
3579 16 Jan 17 peter 67         instance_ = new RNG;
3579 16 Jan 17 peter 68     }
2770 11 Jul 12 peter 69     return instance_;
752 17 Feb 07 jari 70   }
752 17 Feb 07 jari 71
752 17 Feb 07 jari 72
1271 09 Apr 08 peter 73   unsigned long RNG::max(void) const
718 26 Dec 06 jari 74   {
2770 11 Jul 12 peter 75     return gsl_rng_max(rng());
718 26 Dec 06 jari 76   }
718 26 Dec 06 jari 77
718 26 Dec 06 jari 78
1271 09 Apr 08 peter 79   unsigned long RNG::min(void) const
718 26 Dec 06 jari 80   {
2770 11 Jul 12 peter 81     return gsl_rng_min(rng());
718 26 Dec 06 jari 82   }
718 26 Dec 06 jari 83
718 26 Dec 06 jari 84
718 26 Dec 06 jari 85   std::string RNG::name(void) const
718 26 Dec 06 jari 86   {
2770 11 Jul 12 peter 87     return gsl_rng_name(rng());
718 26 Dec 06 jari 88   }
718 26 Dec 06 jari 89
718 26 Dec 06 jari 90
718 26 Dec 06 jari 91   const gsl_rng* RNG::rng(void) const
718 26 Dec 06 jari 92   {
2770 11 Jul 12 peter 93     if (rng_.get()==NULL)
2770 11 Jul 12 peter 94       rng_alloc();
2770 11 Jul 12 peter 95     return rng_.get();
718 26 Dec 06 jari 96   }
718 26 Dec 06 jari 97
718 26 Dec 06 jari 98
2770 11 Jul 12 peter 99   void RNG::rng_alloc(void) const
2770 11 Jul 12 peter 100   {
2770 11 Jul 12 peter 101     assert(rng_.get()==NULL);
2770 11 Jul 12 peter 102     gsl_rng* rng = gsl_rng_alloc(gsl_rng_default);
2770 11 Jul 12 peter 103     if (!rng)
2770 11 Jul 12 peter 104       throw utility::GSL_error("RNG failed to allocate memory");
3959 30 Jul 20 peter 105     std::unique_lock<std::mutex> lock(mutex_);
2802 29 Jul 12 peter 106     gsl_rng_set(rng, seed_);
2802 29 Jul 12 peter 107     // bump seed to avoid subsequent gsl_rng to be identical
2802 29 Jul 12 peter 108     ++seed_;
2789 25 Jul 12 peter 109     // rng_ owns rng and takes care of deallocation
2770 11 Jul 12 peter 110     rng_.reset(rng);
2838 16 Sep 12 peter 111   } // lock is released here
2770 11 Jul 12 peter 112
2770 11 Jul 12 peter 113
1271 09 Apr 08 peter 114   void RNG::seed(unsigned long s) const
718 26 Dec 06 jari 115   {
3959 30 Jul 20 peter 116     std::unique_lock<std::mutex> lock(mutex_);
2770 11 Jul 12 peter 117     gsl_rng_set(rng(),s);
2802 29 Jul 12 peter 118     seed_ = s+1;
2838 16 Sep 12 peter 119   } // lock is released here
718 26 Dec 06 jari 120
718 26 Dec 06 jari 121
1271 09 Apr 08 peter 122   unsigned long RNG::seed_from_devurandom(void)
374 05 Aug 05 jari 123   {
2590 25 Oct 11 peter 124     std::ifstream is("/dev/urandom", std::ios::binary);
3902 05 May 20 peter 125     unsigned long s=0;
3902 05 May 20 peter 126     utility::binary_read(is, s);
374 05 Aug 05 jari 127     is.close();
374 05 Aug 05 jari 128     seed(s);
374 05 Aug 05 jari 129     return s;
374 05 Aug 05 jari 130   }
374 05 Aug 05 jari 131
374 05 Aug 05 jari 132
674 10 Oct 06 peter 133   int RNG::set_state(const RNG_state& state)
674 10 Oct 06 peter 134   {
2770 11 Jul 12 peter 135     if (rng_.get()==NULL)
2770 11 Jul 12 peter 136       rng_alloc();
2804 31 Jul 12 peter 137     if (gsl_rng_memcpy(rng_.get(), state.rng()))
2804 31 Jul 12 peter 138       throw utility::GSL_error("yat::random::RNG::set_state failed");
2804 31 Jul 12 peter 139     return 0;
674 10 Oct 06 peter 140   }
674 10 Oct 06 peter 141
718 26 Dec 06 jari 142   // --------------------- RNG_state ----------------------------------
674 10 Oct 06 peter 143
674 10 Oct 06 peter 144   RNG_state::RNG_state(const RNG* rng)
674 10 Oct 06 peter 145   {
1614 05 Nov 08 peter 146     clone(*rng->rng());
674 10 Oct 06 peter 147   }
674 10 Oct 06 peter 148
2770 11 Jul 12 peter 149
1614 05 Nov 08 peter 150   RNG_state::RNG_state(const RNG_state& state)
1614 05 Nov 08 peter 151   {
1614 05 Nov 08 peter 152     clone(*state.rng());
1614 05 Nov 08 peter 153   }
1614 05 Nov 08 peter 154
2770 11 Jul 12 peter 155
674 10 Oct 06 peter 156   RNG_state::~RNG_state(void)
674 10 Oct 06 peter 157   {
674 10 Oct 06 peter 158     gsl_rng_free(rng_);
674 10 Oct 06 peter 159     rng_=NULL;
674 10 Oct 06 peter 160   }
674 10 Oct 06 peter 161
718 26 Dec 06 jari 162   const gsl_rng* RNG_state::rng(void) const
718 26 Dec 06 jari 163   {
718 26 Dec 06 jari 164     return rng_;
718 26 Dec 06 jari 165   }
718 26 Dec 06 jari 166
1614 05 Nov 08 peter 167
1614 05 Nov 08 peter 168   void RNG_state::clone(const gsl_rng& rng)
1614 05 Nov 08 peter 169   {
1614 05 Nov 08 peter 170     assert(rng_!=&rng);
1614 05 Nov 08 peter 171     if (!(rng_ = gsl_rng_clone(&rng)))
2804 31 Jul 12 peter 172       throw utility::GSL_error("RNG_state::clone failed to allocate memory");
1614 05 Nov 08 peter 173   }
1614 05 Nov 08 peter 174
1614 05 Nov 08 peter 175   RNG_state& RNG_state::operator=(const RNG_state& rhs)
1614 05 Nov 08 peter 176   {
1614 05 Nov 08 peter 177     if (this != &rhs) {
1614 05 Nov 08 peter 178       gsl_rng_free(rng_);
1614 05 Nov 08 peter 179       clone(*rhs.rng());
1614 05 Nov 08 peter 180     }
1614 05 Nov 08 peter 181     return *this;
1614 05 Nov 08 peter 182   }
1614 05 Nov 08 peter 183
706 19 Dec 06 jari 184   // --------------------- Discrete distribtuions ---------------------
674 10 Oct 06 peter 185
706 19 Dec 06 jari 186   Discrete::Discrete(void)
706 19 Dec 06 jari 187     : rng_(RNG::instance())
706 19 Dec 06 jari 188   {
706 19 Dec 06 jari 189   }
706 19 Dec 06 jari 190
2804 31 Jul 12 peter 191
706 19 Dec 06 jari 192   Discrete::~Discrete(void)
706 19 Dec 06 jari 193   {
706 19 Dec 06 jari 194   }
706 19 Dec 06 jari 195
706 19 Dec 06 jari 196
1271 09 Apr 08 peter 197   void Discrete::seed(unsigned long s) const
718 26 Dec 06 jari 198   {
718 26 Dec 06 jari 199     rng_->seed(s);
718 26 Dec 06 jari 200   }
718 26 Dec 06 jari 201
718 26 Dec 06 jari 202
2804 31 Jul 12 peter 203   unsigned long Discrete::seed_from_devurandom(void)
2804 31 Jul 12 peter 204   {
2804 31 Jul 12 peter 205     return rng_->seed_from_devurandom();
1271 09 Apr 08 peter 206   }
1271 09 Apr 08 peter 207
1271 09 Apr 08 peter 208
2889 07 Dec 12 peter 209   Binomial::Binomial(double p, unsigned int n)
2890 07 Dec 12 peter 210     : Discrete(), p_(p), n_(n)
2889 07 Dec 12 peter 211   {
2889 07 Dec 12 peter 212   }
2889 07 Dec 12 peter 213
2889 07 Dec 12 peter 214
2889 07 Dec 12 peter 215   unsigned long Binomial::operator()(void) const
2889 07 Dec 12 peter 216   {
2889 07 Dec 12 peter 217     return gsl_ran_binomial(rng_->rng(), p_, n_);
2889 07 Dec 12 peter 218   }
2889 07 Dec 12 peter 219
2889 07 Dec 12 peter 220
4009 19 Oct 20 peter 221   DiscreteGeneral::DiscreteGeneral(void)
4009 19 Oct 20 peter 222     : gen_(nullptr)
4009 19 Oct 20 peter 223   {}
4009 19 Oct 20 peter 224
4009 19 Oct 20 peter 225
371 05 Aug 05 peter 226   DiscreteGeneral::DiscreteGeneral(const statistics::Histogram& hist)
1616 05 Nov 08 peter 227     : gen_(NULL)
371 05 Aug 05 peter 228   {
1616 05 Nov 08 peter 229     p_.reserve(hist.nof_bins());
2804 31 Jul 12 peter 230     for (size_t i=0; i<hist.nof_bins(); i++)
1616 05 Nov 08 peter 231       p_.push_back(hist[i]);
1616 05 Nov 08 peter 232     preproc();
371 05 Aug 05 peter 233   }
371 05 Aug 05 peter 234
374 05 Aug 05 jari 235
1616 05 Nov 08 peter 236   DiscreteGeneral::DiscreteGeneral(const DiscreteGeneral& other)
2090 21 Oct 09 peter 237     : Discrete(other), gen_(NULL), p_(other.p_)
1616 05 Nov 08 peter 238   {
1616 05 Nov 08 peter 239     preproc();
1616 05 Nov 08 peter 240   }
1616 05 Nov 08 peter 241
1616 05 Nov 08 peter 242
4010 19 Oct 20 peter 243   DiscreteGeneral::DiscreteGeneral(DiscreteGeneral&& other)
4382 13 Oct 23 peter 244     : Discrete(std::move(other)), gen_(other.gen_), p_(std::move(other.p_))
4010 19 Oct 20 peter 245   {
4382 13 Oct 23 peter 246     other.gen_ = nullptr;
4010 19 Oct 20 peter 247   }
4010 19 Oct 20 peter 248
4010 19 Oct 20 peter 249
371 05 Aug 05 peter 250   DiscreteGeneral::~DiscreteGeneral(void)
371 05 Aug 05 peter 251   {
1616 05 Nov 08 peter 252     free();
1616 05 Nov 08 peter 253   }
1616 05 Nov 08 peter 254
1616 05 Nov 08 peter 255
1616 05 Nov 08 peter 256   void DiscreteGeneral::free(void)
1616 05 Nov 08 peter 257   {
374 05 Aug 05 jari 258     if (gen_)
371 05 Aug 05 peter 259       gsl_ran_discrete_free( gen_ );
371 05 Aug 05 peter 260     gen_ = NULL;
371 05 Aug 05 peter 261   }
371 05 Aug 05 peter 262
706 19 Dec 06 jari 263
1616 05 Nov 08 peter 264   void DiscreteGeneral::preproc(void)
1616 05 Nov 08 peter 265   {
1616 05 Nov 08 peter 266     assert(!gen_);
1616 05 Nov 08 peter 267     assert(p_.size());
1616 05 Nov 08 peter 268     gen_ = gsl_ran_discrete_preproc( p_.size(), &p_.front() );
1616 05 Nov 08 peter 269     if (!gen_)
1616 05 Nov 08 peter 270       throw utility::GSL_error("DiscreteGeneral failed to setup generator.");
1616 05 Nov 08 peter 271   }
1616 05 Nov 08 peter 272
1616 05 Nov 08 peter 273
1616 05 Nov 08 peter 274   DiscreteGeneral& DiscreteGeneral::operator=(const DiscreteGeneral& rhs)
1616 05 Nov 08 peter 275   {
4010 19 Oct 20 peter 276     DiscreteGeneral tmp(rhs);
4010 19 Oct 20 peter 277     swap(*this, tmp);
1616 05 Nov 08 peter 278     return *this;
1616 05 Nov 08 peter 279   }
1616 05 Nov 08 peter 280
1616 05 Nov 08 peter 281
4010 19 Oct 20 peter 282   DiscreteGeneral& DiscreteGeneral::operator=(DiscreteGeneral&& rhs)
4010 19 Oct 20 peter 283   {
4010 19 Oct 20 peter 284     DiscreteGeneral tmp(std::move(rhs));
4010 19 Oct 20 peter 285     swap(*this, tmp);
4010 19 Oct 20 peter 286     return *this;
4010 19 Oct 20 peter 287   }
4010 19 Oct 20 peter 288
4010 19 Oct 20 peter 289
1271 09 Apr 08 peter 290   unsigned long DiscreteGeneral::operator()(void) const
718 26 Dec 06 jari 291   {
718 26 Dec 06 jari 292     return gsl_ran_discrete(rng_->rng(), gen_);
718 26 Dec 06 jari 293   }
718 26 Dec 06 jari 294
718 26 Dec 06 jari 295
4010 19 Oct 20 peter 296   void swap(DiscreteGeneral& lhs, DiscreteGeneral& rhs)
4010 19 Oct 20 peter 297   {
4010 19 Oct 20 peter 298     std::swap(lhs.gen_, rhs.gen_);
4010 19 Oct 20 peter 299     std::swap(lhs.p_, rhs.p_);
4010 19 Oct 20 peter 300   }
4010 19 Oct 20 peter 301
4010 19 Oct 20 peter 302
1271 09 Apr 08 peter 303   DiscreteUniform::DiscreteUniform(unsigned long n)
718 26 Dec 06 jari 304     : range_(n)
718 26 Dec 06 jari 305   {
738 07 Jan 07 jari 306     if (range_>rng_->max()) {
3074 03 Sep 13 peter 307       std::ostringstream ss;
3074 03 Sep 13 peter 308       ss << "DiscreteUniform::DiscreteUniform: ";
738 07 Jan 07 jari 309       ss << n << " is too large for RNG " << rng_->name();
3074 03 Sep 13 peter 310       ss << "; maximal argument is " << rng_->max();
738 07 Jan 07 jari 311       throw utility::GSL_error(ss.str());
738 07 Jan 07 jari 312     }
718 26 Dec 06 jari 313   }
718 26 Dec 06 jari 314
718 26 Dec 06 jari 315
1271 09 Apr 08 peter 316   unsigned long DiscreteUniform::operator()(void) const
718 26 Dec 06 jari 317   {
738 07 Jan 07 jari 318     return (range_ ?
2804 31 Jul 12 peter 319             gsl_rng_uniform_int(rng_->rng(),range_) : gsl_rng_get(rng_->rng()));
718 26 Dec 06 jari 320   }
718 26 Dec 06 jari 321
718 26 Dec 06 jari 322
1271 09 Apr 08 peter 323   unsigned long DiscreteUniform::operator()(unsigned long n) const
718 26 Dec 06 jari 324   {
738 07 Jan 07 jari 325     // making sure that n is not larger than the range of the
738 07 Jan 07 jari 326     // underlying RNG
738 07 Jan 07 jari 327     if (n>rng_->max()) {
3074 03 Sep 13 peter 328       std::ostringstream ss;
3074 03 Sep 13 peter 329       ss << "DiscreteUniform::operator(unsigned long): ";
738 07 Jan 07 jari 330       ss << n << " is too large for RNG " << rng_->name();
3074 03 Sep 13 peter 331       ss << "; maximal argument is " << rng_->max();
738 07 Jan 07 jari 332       throw utility::GSL_error(ss.str());
738 07 Jan 07 jari 333     }
738 07 Jan 07 jari 334     return gsl_rng_uniform_int(rng_->rng(),n);
718 26 Dec 06 jari 335   }
718 26 Dec 06 jari 336
718 26 Dec 06 jari 337
3894 27 Mar 20 peter 338   unsigned long int DiscreteUniform::min(void) const
3894 27 Mar 20 peter 339   {
3894 27 Mar 20 peter 340     return range_ ? 0 : rng_->min();
3894 27 Mar 20 peter 341   }
3894 27 Mar 20 peter 342
3894 27 Mar 20 peter 343
3894 27 Mar 20 peter 344   unsigned long int DiscreteUniform::max(void) const
3894 27 Mar 20 peter 345   {
3894 27 Mar 20 peter 346     return range_ ? (range_-1) : rng_->max();
3894 27 Mar 20 peter 347   }
3894 27 Mar 20 peter 348
3894 27 Mar 20 peter 349
3446 02 Dec 15 peter 350   Geometric::Geometric(double p)
3446 02 Dec 15 peter 351     : p_(p)
3446 02 Dec 15 peter 352   {}
3446 02 Dec 15 peter 353
3446 02 Dec 15 peter 354
3446 02 Dec 15 peter 355   unsigned long int Geometric::operator()(void) const
3446 02 Dec 15 peter 356   {
4200 19 Aug 22 peter 357     return gsl_ran_geometric (rng_->rng(), p_);
3446 02 Dec 15 peter 358   }
3446 02 Dec 15 peter 359
3446 02 Dec 15 peter 360
3446 02 Dec 15 peter 361   unsigned long int Geometric::operator()(double p) const
3446 02 Dec 15 peter 362   {
3446 02 Dec 15 peter 363     return gsl_ran_geometric (rng_->rng(), p);
3446 02 Dec 15 peter 364   }
3446 02 Dec 15 peter 365
3446 02 Dec 15 peter 366
3465 10 Feb 16 peter 367   HyperGeometric::HyperGeometric(void)
3465 10 Feb 16 peter 368   {}
3465 10 Feb 16 peter 369
3465 10 Feb 16 peter 370
3465 10 Feb 16 peter 371   HyperGeometric::HyperGeometric(unsigned int n1, unsigned int n2,
3465 10 Feb 16 peter 372                                  unsigned int t)
3465 10 Feb 16 peter 373     : n1_(n1), n2_(n2), t_(t)
3465 10 Feb 16 peter 374   {}
3465 10 Feb 16 peter 375
3465 10 Feb 16 peter 376
3465 10 Feb 16 peter 377   unsigned long int HyperGeometric::operator()(void) const
3465 10 Feb 16 peter 378   {
3465 10 Feb 16 peter 379     return (*this)(n1_, n2_, t_);
3465 10 Feb 16 peter 380   }
3465 10 Feb 16 peter 381
3465 10 Feb 16 peter 382
3465 10 Feb 16 peter 383   unsigned long int HyperGeometric::operator()(unsigned int n1,
3465 10 Feb 16 peter 384                                                unsigned int n2,
3465 10 Feb 16 peter 385                                                unsigned int t) const
3465 10 Feb 16 peter 386   {
3465 10 Feb 16 peter 387     return gsl_ran_hypergeometric(rng_->rng(), n1, n2, t);
3465 10 Feb 16 peter 388   }
3465 10 Feb 16 peter 389
3465 10 Feb 16 peter 390
3469 29 Feb 16 peter 391   NegativeHyperGeometric::NegativeHyperGeometric(void)
3469 29 Feb 16 peter 392   {}
3469 29 Feb 16 peter 393
3469 29 Feb 16 peter 394
3469 29 Feb 16 peter 395   NegativeHyperGeometric::NegativeHyperGeometric(unsigned int n1,
3469 29 Feb 16 peter 396                                                  unsigned int n2, unsigned int t)
3469 29 Feb 16 peter 397     : n1_(n1), n2_(n2), t_(t)
3469 29 Feb 16 peter 398   {}
3469 29 Feb 16 peter 399
3469 29 Feb 16 peter 400
3469 29 Feb 16 peter 401   unsigned long int NegativeHyperGeometric::operator()(void) const
3469 29 Feb 16 peter 402   {
3469 29 Feb 16 peter 403     return (*this)(n1_, n2_, t_);
3469 29 Feb 16 peter 404   }
3469 29 Feb 16 peter 405
3469 29 Feb 16 peter 406
3469 29 Feb 16 peter 407   unsigned long int NegativeHyperGeometric::operator()(unsigned int n1,
3469 29 Feb 16 peter 408                                                        unsigned int n2,
3469 29 Feb 16 peter 409                                                        unsigned int t) const
3469 29 Feb 16 peter 410   {
3469 29 Feb 16 peter 411     assert(t <= n2);
3469 29 Feb 16 peter 412
3469 29 Feb 16 peter 413     // NHG can be described as an array with n1 true and n2 false, and
3469 29 Feb 16 peter 414     // NHG(n1, n2, t) is the number of true left of the t:th false. By
3469 29 Feb 16 peter 415     // symmetry number of true right of the t:th false is NHG(n1, n2,
3469 29 Feb 16 peter 416     // n2-t+1) since t:th false counting from left is the (n2-t+1):th
3469 29 Feb 16 peter 417     // false counting from right.
3469 29 Feb 16 peter 418
3469 29 Feb 16 peter 419     // When t is larger than midpoint (2*t > n2+1) we use this
3469 29 Feb 16 peter 420     // symmetry to speed things up.
3469 29 Feb 16 peter 421     if (t > (n2+1)/2)
3469 29 Feb 16 peter 422       return n1 - (*this)(n1, n2, n2-t+1);
3469 29 Feb 16 peter 423
3469 29 Feb 16 peter 424     ContinuousUniform uniform;
3469 29 Feb 16 peter 425     unsigned long int k = 0;
3469 29 Feb 16 peter 426     while (t) {
3469 29 Feb 16 peter 427       assert(n1 + n2);
3469 29 Feb 16 peter 428       double x = uniform();
3469 29 Feb 16 peter 429       if (x * (n1+n2) < n1) {
3469 29 Feb 16 peter 430         --n1;
3469 29 Feb 16 peter 431         ++k;
3469 29 Feb 16 peter 432         if (!n1)
3469 29 Feb 16 peter 433           return k;
3469 29 Feb 16 peter 434       }
3469 29 Feb 16 peter 435       else {
3469 29 Feb 16 peter 436         --t;
3469 29 Feb 16 peter 437         --n2;
3469 29 Feb 16 peter 438       }
3469 29 Feb 16 peter 439     }
3469 29 Feb 16 peter 440     return k;
3469 29 Feb 16 peter 441   }
3469 29 Feb 16 peter 442
3469 29 Feb 16 peter 443
706 19 Dec 06 jari 444   Poisson::Poisson(const double m)
706 19 Dec 06 jari 445     : m_(m)
706 19 Dec 06 jari 446   {
706 19 Dec 06 jari 447   }
706 19 Dec 06 jari 448
1271 09 Apr 08 peter 449   unsigned long Poisson::operator()(void) const
718 26 Dec 06 jari 450   {
718 26 Dec 06 jari 451     return gsl_ran_poisson(rng_->rng(), m_);
718 26 Dec 06 jari 452   }
718 26 Dec 06 jari 453
718 26 Dec 06 jari 454
1271 09 Apr 08 peter 455   unsigned long Poisson::operator()(const double m) const
718 26 Dec 06 jari 456   {
718 26 Dec 06 jari 457     return gsl_ran_poisson(rng_->rng(), m);
718 26 Dec 06 jari 458   }
718 26 Dec 06 jari 459
4147 28 Feb 22 peter 460
4147 28 Feb 22 peter 461   NegativeBinomial::NegativeBinomial(double p, double r)
4147 28 Feb 22 peter 462     : p_(p), r_(r)
4147 28 Feb 22 peter 463   {
4147 28 Feb 22 peter 464   }
4147 28 Feb 22 peter 465
4147 28 Feb 22 peter 466   unsigned long NegativeBinomial::operator()(void) const
4147 28 Feb 22 peter 467   {
4147 28 Feb 22 peter 468     return gsl_ran_negative_binomial(rng_->rng(), p_, r_);
4147 28 Feb 22 peter 469   }
4147 28 Feb 22 peter 470
4147 28 Feb 22 peter 471
4147 28 Feb 22 peter 472   unsigned long NegativeBinomial::operator()(double p, double r) const
4147 28 Feb 22 peter 473   {
4147 28 Feb 22 peter 474     return gsl_ran_negative_binomial(rng_->rng(), p, r);
4147 28 Feb 22 peter 475   }
4147 28 Feb 22 peter 476
706 19 Dec 06 jari 477   // --------------------- Continuous distribtuions ---------------------
706 19 Dec 06 jari 478
706 19 Dec 06 jari 479   Continuous::Continuous(void)
706 19 Dec 06 jari 480     : rng_(RNG::instance())
706 19 Dec 06 jari 481   {
706 19 Dec 06 jari 482   }
706 19 Dec 06 jari 483
706 19 Dec 06 jari 484
706 19 Dec 06 jari 485   Continuous::~Continuous(void)
706 19 Dec 06 jari 486   {
706 19 Dec 06 jari 487   }
706 19 Dec 06 jari 488
706 19 Dec 06 jari 489
1271 09 Apr 08 peter 490   void Continuous::seed(unsigned long s) const
718 26 Dec 06 jari 491   {
718 26 Dec 06 jari 492     rng_->seed(s);
718 26 Dec 06 jari 493   }
718 26 Dec 06 jari 494
718 26 Dec 06 jari 495
2578 04 Oct 11 peter 496   unsigned long Continuous::seed_from_devurandom(void)
2578 04 Oct 11 peter 497   {
2578 04 Oct 11 peter 498     return rng_->seed_from_devurandom();
2578 04 Oct 11 peter 499   }
2578 04 Oct 11 peter 500
2578 04 Oct 11 peter 501
706 19 Dec 06 jari 502   ContinuousGeneral::ContinuousGeneral(const statistics::Histogram& hist)
706 19 Dec 06 jari 503     : discrete_(DiscreteGeneral(hist)), hist_(hist)
706 19 Dec 06 jari 504   {
706 19 Dec 06 jari 505   }
706 19 Dec 06 jari 506
706 19 Dec 06 jari 507
718 26 Dec 06 jari 508   double ContinuousGeneral::operator()(void) const
718 26 Dec 06 jari 509   {
718 26 Dec 06 jari 510     return hist_.observation_value(discrete_())+(u_()-0.5)*hist_.spacing();
718 26 Dec 06 jari 511   }
718 26 Dec 06 jari 512
718 26 Dec 06 jari 513   double ContinuousUniform::operator()(void) const
718 26 Dec 06 jari 514   {
718 26 Dec 06 jari 515     return gsl_rng_uniform(rng_->rng());
718 26 Dec 06 jari 516   }
718 26 Dec 06 jari 517
718 26 Dec 06 jari 518
706 19 Dec 06 jari 519   Exponential::Exponential(const double m)
706 19 Dec 06 jari 520     : m_(m)
706 19 Dec 06 jari 521   {
706 19 Dec 06 jari 522   }
706 19 Dec 06 jari 523
706 19 Dec 06 jari 524
718 26 Dec 06 jari 525   double Exponential::operator()(void) const
718 26 Dec 06 jari 526   {
718 26 Dec 06 jari 527     return gsl_ran_exponential(rng_->rng(), m_);
718 26 Dec 06 jari 528   }
718 26 Dec 06 jari 529
718 26 Dec 06 jari 530
718 26 Dec 06 jari 531   double Exponential::operator()(const double m) const
718 26 Dec 06 jari 532   {
718 26 Dec 06 jari 533     return gsl_ran_exponential(rng_->rng(), m);
718 26 Dec 06 jari 534   }
718 26 Dec 06 jari 535
718 26 Dec 06 jari 536
706 19 Dec 06 jari 537   Gaussian::Gaussian(const double s, const double m)
706 19 Dec 06 jari 538     : m_(m), s_(s)
706 19 Dec 06 jari 539   {
706 19 Dec 06 jari 540   }
706 19 Dec 06 jari 541
718 26 Dec 06 jari 542
718 26 Dec 06 jari 543   double Gaussian::operator()(void) const
718 26 Dec 06 jari 544   {
718 26 Dec 06 jari 545     return gsl_ran_gaussian(rng_->rng(), s_)+m_;
718 26 Dec 06 jari 546   }
718 26 Dec 06 jari 547
718 26 Dec 06 jari 548
718 26 Dec 06 jari 549   double Gaussian::operator()(const double s) const
718 26 Dec 06 jari 550   {
718 26 Dec 06 jari 551     return gsl_ran_gaussian(rng_->rng(), s);
718 26 Dec 06 jari 552   }
718 26 Dec 06 jari 553
718 26 Dec 06 jari 554
718 26 Dec 06 jari 555   double Gaussian::operator()(const double s, const double m) const
718 26 Dec 06 jari 556   {
718 26 Dec 06 jari 557     return gsl_ran_gaussian(rng_->rng(), s)+m;
718 26 Dec 06 jari 558   }
718 26 Dec 06 jari 559
4281 29 Jan 23 peter 560
4281 29 Jan 23 peter 561   MultivariateGaussian::MultivariateGaussian(const utility::VectorBase& m,
4281 29 Jan 23 peter 562                                              const utility::CholeskyDecomposer& C)
4281 29 Jan 23 peter 563     : mean_(m), cholesky_decomposer_(C)
4281 29 Jan 23 peter 564   {}
4281 29 Jan 23 peter 565
4281 29 Jan 23 peter 566
4281 29 Jan 23 peter 567   MultivariateGaussian::MultivariateGaussian(const utility::VectorBase& m,
4281 29 Jan 23 peter 568                                              utility::CholeskyDecomposer&& C)
4281 29 Jan 23 peter 569     : mean_(m), cholesky_decomposer_(std::move(C))
4281 29 Jan 23 peter 570   {}
4281 29 Jan 23 peter 571
4281 29 Jan 23 peter 572
4281 29 Jan 23 peter 573   void MultivariateGaussian::operator()(utility::VectorMutable& x)
4281 29 Jan 23 peter 574   {
4281 29 Jan 23 peter 575     (*this)(mean_, cholesky_decomposer_, x);
4281 29 Jan 23 peter 576   }
4281 29 Jan 23 peter 577
4281 29 Jan 23 peter 578
4281 29 Jan 23 peter 579   void MultivariateGaussian::operator()(const utility::VectorBase& m,
4281 29 Jan 23 peter 580                                         const utility::CholeskyDecomposer& C,
4281 29 Jan 23 peter 581                                         utility::VectorMutable& x)
4281 29 Jan 23 peter 582   {
4281 29 Jan 23 peter 583     // should perhaps be stored as a member variable (in a base class)
4281 29 Jan 23 peter 584     RNG* rng = RNG::instance();
4281 29 Jan 23 peter 585     if (int status = gsl_ran_multivariate_gaussian(rng->rng(),
4281 29 Jan 23 peter 586                                                    m.gsl_vector_p(),
4281 29 Jan 23 peter 587                                                    C.L().gsl_matrix_p(),
4281 29 Jan 23 peter 588                                                    x.gsl_vector_p()))
4281 29 Jan 23 peter 589       throw utility::GSL_error("MultivariateGaussian failed", status);
4281 29 Jan 23 peter 590   }
4281 29 Jan 23 peter 591
680 11 Oct 06 jari 592 }}} // of namespace random, yat, and theplu

Generated on Sat Oct 14 02:23:02 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)