LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
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 <limits>
30#include <string>
31#include <exception>
32
33#include "boost/format.hpp"
34
35#include "gsl/gsl_randist.h"
36
37#include "lsst/pex/exceptions.h"
38
40
41namespace ex = lsst::pex::exceptions;
42
43namespace lsst {
44namespace afw {
45namespace 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
54char const *const Random::_algorithmNames[Random::NUM_ALGORITHMS] = {
55 "MT19937", "RANLXS0", "RANLXS1", "RANLXS2", "RANLXD1", "RANLXD2", "RANLUX",
56 "RANLUX389", "CMRG", "MRG", "TAUS", "TAUS2", "GFSR4"};
57
58char const *const Random::_algorithmEnvVarName = "LSST_RNG_ALGORITHM";
59char const *const Random::_seedEnvVarName = "LSST_RNG_SEED";
60
61// -- Private helper functions --------
62
63void Random::initialize() {
64 ::gsl_rng *rng = ::gsl_rng_alloc(_gslRngTypes[_algorithm]);
65 if (rng == nullptr) {
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
75void 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
89Random::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
96Random::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
113void 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
124std::size_t Random::getStateSize() const { return ::gsl_rng_size(_rng.get()); }
125
126// -- Accessors --------
127
128Random::Algorithm Random::getAlgorithm() const { return _algorithm; }
129
130std::string Random::getAlgorithmName() const { return std::string(_algorithmNames[_algorithm]); }
131
133 static std::vector<std::string> names;
134 if (names.size() == 0) {
135 for (auto _algorithmName : _algorithmNames) {
136 names.emplace_back(_algorithmName);
137 }
138 }
139 return names;
140}
141
142unsigned long Random::getSeed() const { return _seed; }
143
144// -- Mutators: generating random numbers --------
145
146double Random::uniform() { return ::gsl_rng_uniform(_rng.get()); }
147
148double Random::uniformPos() { return ::gsl_rng_uniform_pos(_rng.get()); }
149
150unsigned 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
159double Random::flat(double const a, double const b) { return ::gsl_ran_flat(_rng.get(), a, b); }
160
161double Random::gaussian() { return ::gsl_ran_gaussian_ziggurat(_rng.get(), 1.0); }
162
163double Random::chisq(double nu) { return ::gsl_ran_chisq(_rng.get(), nu); }
164
165double 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
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
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::string State
Accessors for the opaque state of the random number generator.
Definition Random.h:154
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)
T reset(T... args)
T size(T... args)