LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
Random.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008-2016 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 /*
26  * Random number generator implementaion.
27  */
28 
29 #include <cstdlib>
30 #include <limits>
31 
32 #include "boost/format.hpp"
33 
34 #include "gsl/gsl_errno.h"
35 #include "gsl/gsl_randist.h"
36 
37 #include "lsst/pex/exceptions.h"
38 
39 #include "lsst/afw/math/Random.h"
40 
41 namespace ex = lsst::pex::exceptions;
42 
43 namespace lsst {
44 namespace afw {
45 namespace math {
46 
47 // -- Static data --------
48 
49 ::gsl_rng_type const *const Random::_gslRngTypes[Random::NUM_ALGORITHMS] = {
50  ::gsl_rng_mt19937, ::gsl_rng_ranlxs0, ::gsl_rng_ranlxs1, ::gsl_rng_ranlxs2, ::gsl_rng_ranlxd1,
51  ::gsl_rng_ranlxd2, ::gsl_rng_ranlux, ::gsl_rng_ranlux389, ::gsl_rng_cmrg, ::gsl_rng_mrg,
52  ::gsl_rng_taus, ::gsl_rng_taus2, ::gsl_rng_gfsr4};
53 
54 char const *const Random::_algorithmNames[Random::NUM_ALGORITHMS] = {
55  "MT19937", "RANLXS0", "RANLXS1", "RANLXS2", "RANLXD1", "RANLXD2", "RANLUX",
56  "RANLUX389", "CMRG", "MRG", "TAUS", "TAUS2", "GFSR4"};
57 
58 char const *const Random::_algorithmEnvVarName = "LSST_RNG_ALGORITHM";
59 char const *const Random::_seedEnvVarName = "LSST_RNG_SEED";
60 
61 // -- Private helper functions --------
62 
63 void Random::initialize() {
64  ::gsl_rng *rng = ::gsl_rng_alloc(_gslRngTypes[_algorithm]);
65  if (rng == 0) {
66  throw std::bad_alloc();
67  }
68  // This seed is guaranteed to be non-zero.
69  // We want to give a non-zero seed to GSL to avoid it choosing its own.
70  unsigned long int useSeed = _seed == 0 ? std::numeric_limits<unsigned long int>::max() : _seed;
71  ::gsl_rng_set(rng, useSeed);
72  _rng.reset(rng, ::gsl_rng_free);
73 }
74 
75 void Random::initialize(std::string const &algorithm) {
76  // linear search (the number of algorithms is small)
77  for (int i = 0; i < NUM_ALGORITHMS; ++i) {
78  if (_algorithmNames[i] == algorithm) {
79  _algorithm = static_cast<Algorithm>(i);
80  initialize();
81  return;
82  }
83  }
84  throw LSST_EXCEPT(ex::InvalidParameterError, "RNG algorithm " + algorithm + " is not supported");
85 }
86 
87 // -- Constructor --------
88 
89 Random::Random(Algorithm const algorithm, unsigned long seed) : _rng(), _seed(seed), _algorithm(algorithm) {
90  if (_algorithm < 0 || _algorithm >= NUM_ALGORITHMS) {
91  throw LSST_EXCEPT(ex::InvalidParameterError, "Invalid RNG algorithm");
92  }
93  initialize();
94 }
95 
96 Random::Random(std::string const &algorithm, unsigned long seed) : _rng(), _seed(seed) {
97  initialize(algorithm);
98 }
99 
101  Random rng = *this;
102  rng._rng.reset(::gsl_rng_clone(_rng.get()), ::gsl_rng_free);
103  if (!rng._rng) {
104  throw std::bad_alloc();
105  }
106  return rng;
107 }
108 
110  return State(static_cast<char *>(::gsl_rng_state(_rng.get())), getStateSize());
111 }
112 
113 void Random::setState(State const &state) {
114  if (state.size() != getStateSize()) {
115  throw LSST_EXCEPT(
117  (boost::format("Size of given state vector (%d) does not match expected size (%d)") %
118  state.size() % getStateSize())
119  .str());
120  }
121  std::copy(state.begin(), state.end(), static_cast<char *>(::gsl_rng_state(_rng.get())));
122 }
123 
124 std::size_t Random::getStateSize() const { return ::gsl_rng_size(_rng.get()); }
125 
126 // -- Accessors --------
127 
128 Random::Algorithm Random::getAlgorithm() const { return _algorithm; }
129 
130 std::string Random::getAlgorithmName() const { return std::string(_algorithmNames[_algorithm]); }
131 
133  static std::vector<std::string> names;
134  if (names.size() == 0) {
135  for (int i = 0; i < NUM_ALGORITHMS; ++i) {
136  names.push_back(_algorithmNames[i]);
137  }
138  }
139  return names;
140 }
141 
142 unsigned long Random::getSeed() const { return _seed; }
143 
144 // -- Mutators: generating random numbers --------
145 
146 double Random::uniform() { return ::gsl_rng_uniform(_rng.get()); }
147 
148 double Random::uniformPos() { return ::gsl_rng_uniform_pos(_rng.get()); }
149 
150 unsigned long Random::uniformInt(unsigned long n) {
151  if (n > ::gsl_rng_max(_rng.get()) - ::gsl_rng_min(_rng.get())) {
152  throw LSST_EXCEPT(ex::RangeError, "Desired random number range exceeds generator range");
153  }
154  return ::gsl_rng_uniform_int(_rng.get(), n);
155 }
156 
157 // -- Mutators: computing random variates for various distributions --------
158 
159 double Random::flat(double const a, double const b) { return ::gsl_ran_flat(_rng.get(), a, b); }
160 
161 double Random::gaussian() { return ::gsl_ran_gaussian_ziggurat(_rng.get(), 1.0); }
162 
163 double Random::chisq(double nu) { return ::gsl_ran_chisq(_rng.get(), nu); }
164 
165 double Random::poisson(double mu) { return ::gsl_ran_poisson(_rng.get(), mu); }
166 } // namespace math
167 } // namespace afw
168 } // namespace lsst
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
table::Key< int > b
table::Key< int > a
T begin(T... args)
A class that can be used to generate sequences of random numbers according to a number of different a...
Definition: Random.h:57
Random(Algorithm algorithm=MT19937, unsigned long seed=1)
Creates a random number generator that uses the given algorithm to produce random numbers,...
Definition: Random.cc:89
Algorithm
Identifiers for the list of supported algorithms.
Definition: Random.h:60
@ NUM_ALGORITHMS
Number of supported algorithms.
Definition: Random.h:92
double uniformPos()
Returns a uniformly distributed random double precision floating point number from the generator.
Definition: Random.cc:148
double chisq(double const nu)
Returns a random variate from the chi-squared distribution with nu degrees of freedom.
Definition: Random.cc:163
Random deepCopy() const
Creates a deep copy of this random number generator.
Definition: Random.cc:100
std::string State
Accessors for the opaque state of the random number generator.
Definition: Random.h:154
static std::vector< std::string > const & getAlgorithmNames()
Definition: Random.cc:132
double poisson(double const mu)
Returns a random variate from the poisson distribution with mean mu.
Definition: Random.cc:165
Algorithm getAlgorithm() const
Definition: Random.cc:128
std::string getAlgorithmName() const
Definition: Random.cc:130
unsigned long getSeed() const
Definition: Random.cc:142
double gaussian()
Returns a gaussian random variate with mean 0 and standard deviation 1
Definition: Random.cc:161
void setState(State const &state)
Definition: Random.cc:113
std::size_t getStateSize() const
Definition: Random.cc:124
State getState() const
Definition: Random.cc:109
double uniform()
Returns a uniformly distributed random double precision floating point number from the generator.
Definition: Random.cc:146
unsigned long uniformInt(unsigned long n)
Returns a uniformly distributed random integer from 0 to n-1.
Definition: Random.cc:150
double flat(double const a, double const b)
Returns a random variate from the flat (uniform) distribution on [a, b).
Definition: Random.cc:159
Reports invalid arguments.
Definition: Runtime.h:66
Reports attempts to exceed implementation-defined length limits for some classes.
Definition: Runtime.h:76
Reports when the result of an operation cannot be represented by the destination type.
Definition: Runtime.h:115
T copy(T... args)
T end(T... args)
T get(T... args)
T max(T... args)
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
T push_back(T... args)
T reset(T... args)
T size(T... args)