1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2017 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This 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 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24
25 #include <queso/Defines.h>
26 #include <queso/RngGsl.h>
27 #include <gsl/gsl_randist.h>
28
29 namespace QUESO {
30
31 //! Constructor with seed ---------------------------
RngGsl(int seed,int worldRank)32 RngGsl::RngGsl(int seed, int worldRank)
33 :
34 RngBase(seed,worldRank),
35 m_rng (NULL)
36 {
37 gsl_rng_default_seed = (unsigned long int) m_seed;
38 m_rng = gsl_rng_alloc(gsl_rng_ranlxd2);
39 queso_require_msg(m_rng, "null m_rng");
40
41 //gsl_rng_set(m_rng, gsl_rng_default_seed);
42 #if 0
43 if (m_worldRank == 0) {
44 std::cout << "In RngGsl::constructor():"
45 << "\n m_seed = " << m_seed
46 << "\n internal seed = " << gsl_rng_default_seed
47 //<< "\n first generated sample from uniform distribution = " << gsl_rng_uniform(m_rng)
48 //<< "\n first generated sample from std normal distribution = " << gsl_ran_gaussian(m_rng,1.)
49 << std::endl;
50 }
51 #endif
52 }
53
54 // Destructor ---------------------------------------
~RngGsl()55 RngGsl::~RngGsl()
56 {
57 if (m_rng) gsl_rng_free(m_rng);
58 }
59
60 // Sampling methods ---------------------------------
61 void
resetSeed(int newSeed)62 RngGsl::resetSeed(int newSeed)
63 {
64 RngBase::resetSeed(newSeed);
65 gsl_rng_free(m_rng);
66
67 gsl_rng_default_seed = (unsigned long int) m_seed;
68 m_rng = gsl_rng_alloc(gsl_rng_ranlxd2);
69 queso_require_msg(m_rng, "null m_rng");
70 return;
71 }
72
73 // --------------------------------------------------
74 double
uniformSample()75 RngGsl::uniformSample() const
76 {
77 return gsl_rng_uniform(m_rng);
78 }
79
80 // --------------------------------------------------
81 double
gaussianSample(double stdDev)82 RngGsl::gaussianSample(double stdDev) const
83 {
84 return gsl_ran_gaussian(m_rng,stdDev);
85 }
86
87 // --------------------------------------------------
88 double
betaSample(double alpha,double beta)89 RngGsl::betaSample(double alpha, double beta) const
90 {
91 return gsl_ran_beta(m_rng,alpha,beta);
92 }
93
94 // --------------------------------------------------
95 double
gammaSample(double a,double b)96 RngGsl::gammaSample(double a, double b) const
97 {
98 return gsl_ran_gamma(m_rng,a,b);
99 }
100
101 const gsl_rng*
rng()102 RngGsl::rng() const
103 {
104 return m_rng;
105 }
106
107 } // End namespace QUESO
108