LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
Random.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010 LSST Corporation.
6  *
7  * This product includes software developed by the
8  * LSST Project (http://www.lsst.org/).
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the LSST License Statement and
21  * the GNU General Public License along with this program. If not,
22  * see <http://www.lsstcorp.org/LegalNotices/>.
23  */
24 
25 
32 #include <cstdlib>
33 
34 #include "boost/format.hpp"
35 #include "boost/lexical_cast.hpp"
36 
37 #include "gsl/gsl_errno.h"
38 #include "gsl/gsl_randist.h"
39 
40 #include "lsst/pex/exceptions.h"
41 
42 #include "lsst/afw/math/Random.h"
43 
45 
46 namespace ex = lsst::pex::exceptions;
47 namespace math = lsst::afw::math;
48 
49 
50 // -- Static data --------
51 
52 ::gsl_rng_type const * const math::Random::_gslRngTypes[math::Random::NUM_ALGORITHMS] = {
53  ::gsl_rng_mt19937,
54  ::gsl_rng_ranlxs0,
55  ::gsl_rng_ranlxs1,
56  ::gsl_rng_ranlxs2,
57  ::gsl_rng_ranlxd1,
58  ::gsl_rng_ranlxd2,
59  ::gsl_rng_ranlux,
60  ::gsl_rng_ranlux389,
61  ::gsl_rng_cmrg,
62  ::gsl_rng_mrg,
63  ::gsl_rng_taus,
64  ::gsl_rng_taus2,
65  ::gsl_rng_gfsr4
66 };
67 
69  "MT19937",
70  "RANLXS0",
71  "RANLXS1",
72  "RANLXS2",
73  "RANLXD1",
74  "RANLXD2",
75  "RANLUX",
76  "RANLUX389",
77  "CMRG",
78  "MRG",
79  "TAUS",
80  "TAUS2",
81  "GFSR4"
82 };
83 
84 char const * const math::Random::_algorithmEnvVarName = "LSST_RNG_ALGORITHM";
85 char const * const math::Random::_seedEnvVarName = "LSST_RNG_SEED";
86 
87 
88 // -- Private helper functions --------
89 
97  if (_seed == 0) {
98  throw LSST_EXCEPT(ex::InvalidParameterError,
99  (boost::format("Invalid RNG seed: %lu") % _seed).str());
100  }
101  ::gsl_rng * rng = ::gsl_rng_alloc(_gslRngTypes[_algorithm]);
102  if (rng == 0) {
103  throw LSST_EXCEPT(ex::MemoryError, "gsl_rng_alloc() failed");
104  }
105  ::gsl_rng_set(rng, _seed);
106  _rng.reset(rng, ::gsl_rng_free);
107 }
108 
119 void math::Random::initialize(std::string const & algorithm) {
120  // linear search (the number of algorithms is small)
121  for (int i = 0; i < NUM_ALGORITHMS; ++i) {
122  if (_algorithmNames[i] == algorithm) {
123  _algorithm = static_cast<Algorithm>(i);
124  initialize();
125  return;
126  }
127  }
128  throw LSST_EXCEPT(ex::InvalidParameterError, "RNG algorithm " +
129  algorithm + " is not supported");
130 }
131 
132 
133 // -- Constructor --------
134 
151 math::Random::Random(Algorithm const algorithm, unsigned long seed)
152  : _rng(), _seed(seed), _algorithm(algorithm)
153 {
154  if (_algorithm < 0 || _algorithm >= NUM_ALGORITHMS) {
155  throw LSST_EXCEPT(ex::InvalidParameterError, "Invalid RNG algorithm");
156  }
157  initialize();
158 }
159 
160 
175 math::Random::Random(std::string const & algorithm, unsigned long seed)
176  : _rng(), _seed(seed)
177 {
178  initialize(algorithm);
179 }
180 
181 
201  : _rng(), _seed()
202 {
203  std::string const seed(policy->getString("rngSeed"));
204  try {
205  _seed = boost::lexical_cast<unsigned long>(seed);
206  } catch(boost::bad_lexical_cast &) {
207  throw LSST_EXCEPT(ex::RuntimeError,
208  (boost::format("Invalid \"rngSeed\" policy value: \"%1%\"") % seed).str());
209  }
210  initialize(policy->getString("rngAlgorithm"));
211 }
212 
213 
224  Random rng = *this;
225  rng._rng.reset(::gsl_rng_clone(_rng.get()), ::gsl_rng_free);
226  if (!rng._rng) {
227  throw LSST_EXCEPT(ex::MemoryError, "gsl_rng_clone() failed");
228  }
229  return rng;
230 }
231 
233  return State(static_cast<char*>(::gsl_rng_state(_rng.get())), getStateSize());
234 }
235 
236 void math::Random::setState(State const & state) {
237  if (state.size() != getStateSize()) {
238  throw LSST_EXCEPT(
239  pex::exceptions::LengthError,
240  (boost::format("Size of given state vector (%d) does not match expected size (%d)")
241  % state.size() % getStateSize()).str()
242  );
243  }
244  std::copy(state.begin(), state.end(), static_cast<char*>(::gsl_rng_state(_rng.get())));
245 }
246 
247 std::size_t math::Random::getStateSize() const {
248  return ::gsl_rng_size(_rng.get());
249 }
250 
251 // -- Accessors --------
252 
257  return _algorithm;
258 }
259 
263 std::string math::Random::getAlgorithmName() const {
264  return std::string(_algorithmNames[_algorithm]);
265 }
266 
270 std::vector<std::string> const & math::Random::getAlgorithmNames() {
271  static std::vector<std::string> names;
272  if (names.size() == 0) {
273  for (int i = 0; i < NUM_ALGORITHMS; ++i) {
274  names.push_back(_algorithmNames[i]);
275  }
276  }
277  return names;
278 }
279 
285 unsigned long math::Random::getSeed() const {
286  return _seed;
287 }
288 
289 
290 // -- Mutators: generating random numbers --------
291 
305  return ::gsl_rng_uniform(_rng.get());
306 }
307 
320  return ::gsl_rng_uniform_pos(_rng.get());
321 }
322 
341 unsigned long math::Random::uniformInt(unsigned long n) {
342  if (n > ::gsl_rng_max(_rng.get()) - ::gsl_rng_min(_rng.get())) {
343  throw LSST_EXCEPT(ex::RangeError,
344  "Desired random number range exceeds generator range");
345  }
346  return ::gsl_rng_uniform_int(_rng.get(), n);
347 }
348 
349 // -- Mutators: computing random variates for various distributions --------
350 
358 double math::Random::flat(double const a, double const b) {
359  return ::gsl_ran_flat(_rng.get(), a, b);
360 }
361 
371  return ::gsl_ran_gaussian_ziggurat(_rng.get(), 1.0);
372 }
373 
380 double math::Random::chisq(double nu) {
381  return ::gsl_ran_chisq(_rng.get(), nu);
382 }
383 
384 
390 double math::Random::poisson(double mu
391  ) {
392  return ::gsl_ran_poisson(_rng.get(), mu);
393 }
394 
std::string getAlgorithmName() const
Definition: Random.cc:263
static char const *const _seedEnvVarName
Definition: Random.h:147
Algorithm getAlgorithm() const
Definition: Random.cc:256
double chisq(double const nu)
Definition: Random.cc:380
unsigned long getSeed() const
Definition: Random.cc:285
double flat(double const a, double const b)
Definition: Random.cc:358
a container for holding hierarchical configuration data in memory.
Definition: Policy.h:169
unsigned long _seed
Definition: Random.h:141
SelectEigenView< T >::Type copy(Eigen::EigenBase< T > const &other)
Copy an arbitrary Eigen expression into a new EigenView.
Definition: eigen.h:390
std::size_t getStateSize() const
Definition: Random.cc:247
boost::shared_ptr< Policy > Ptr
Definition: Policy.h:172
std::string State
Definition: Random.h:116
::gsl_rng_type const *const _gslRngTypes[NUM_ALGORITHMS]
Definition: Random.h:144
double poisson(double const nu)
Definition: Random.cc:390
Random(Algorithm algorithm=MT19937, unsigned long seed=1)
Definition: Random.cc:151
unsigned long uniformInt(unsigned long n)
Definition: Random.cc:341
static char const *const _algorithmEnvVarName
Definition: Random.h:146
double uniformPos()
Definition: Random.cc:319
void setState(State const &state)
Definition: Random.cc:236
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
Random deepCopy() const
Definition: Random.cc:223
Random number generator class.
static std::vector< std::string > const & getAlgorithmNames()
Definition: Random.cc:270
static char const *const _algorithmNames[NUM_ALGORITHMS]
Definition: Random.h:145
afw::table::Key< double > b
State getState() const
Definition: Random.cc:232
Include files required for standard LSST Exception handling.
boost::shared_ptr< ::gsl_rng > _rng
Definition: Random.h:140
Algorithm _algorithm
Definition: Random.h:142