Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/RNG/RNGSSE.hpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Interface to RNGSSE random number generators
9 : : \details Interface to RNGSSE random number generators
10 : : */
11 : : // *****************************************************************************
12 : : #ifndef RNGSSE_h
13 : : #define RNGSSE_h
14 : :
15 : : #include <cstring>
16 : : #include <random>
17 : : #include <memory>
18 : :
19 : : #include "NoWarning/beta_distribution.hpp"
20 : : #include <boost/random/gamma_distribution.hpp>
21 : :
22 : : #include "Exception.hpp"
23 : : #include "Macro.hpp"
24 : : #include "Options/RNGSSESeqLen.hpp"
25 : :
26 : : namespace tk {
27 : :
28 : : //! RNGSSE-based random number generator used polymorphically with tk::RNG
29 : : template< class State, typename SeqNumType, unsigned int (*Generate)(State*) >
30 : : class RNGSSE {
31 : :
32 : : private:
33 : : using InitFn = void (*)( State*, SeqNumType );
34 : : using ncomp_t = kw::ncomp::info::expect::type;
35 : :
36 : : //! Adaptor to use a std distribution with the RNGSSE generator
37 : : //! \see C++ concepts: UniformRandomNumberGenerator
38 : : struct Adaptor {
39 : : using result_type = unsigned int;
40 : 392000792 : Adaptor( const std::unique_ptr< State[] >& s, int t ) : str(s), tid(t) {}
41 : 2103515772 : static constexpr result_type min() { return 0u; }
42 : 33845532 : static constexpr result_type max() { return 4294967295u; }
43 : 2069670240 : result_type operator()()
44 : 2069670240 : { return Generate( &str[ static_cast<std::size_t>(tid) ] ); }
45 : : const std::unique_ptr< State[] >& str;
46 : : int tid;
47 : : };
48 : :
49 : : public:
50 : : //! Constructor
51 : : //! \param[in] n Initialize RNG using this many independent streams
52 : : //! \param[in] fnShort RNG initializer function for short streams
53 : : //! \param[in] seqlen Sequence length enum: short, medium or long
54 : : //! \param[in] fnLong RNG initializer function for long streams
55 : : //! \param[in] fnMed RNG initializer function for medium streams
56 : 1060 : explicit RNGSSE( SeqNumType n,
57 : : InitFn fnShort,
58 : : ctr::RNGSSESeqLenType seqlen = ctr::RNGSSESeqLenType::SHORT,
59 : : InitFn fnLong = nullptr,
60 : : InitFn fnMed = nullptr) :
61 : : m_nthreads( n ),
62 : 1060 : m_init( seqlen == ctr::RNGSSESeqLenType::LONG ? fnLong :
63 [ - + ]: 1060 : seqlen == ctr::RNGSSESeqLenType::MEDIUM ? fnMed : fnShort ),
64 [ + - ]: 2131 : m_stream()
65 : : {
66 [ - + ][ - - ]: 1060 : Assert( m_init != nullptr, "nullptr passed to RNGSSE constructor" );
[ - - ][ - - ]
67 [ + + ][ + - ]: 1060 : Assert( n > 0, "Need at least one thread" );
[ + - ][ + - ]
68 : : // Allocate array of stream-pointers for threads
69 [ + - ]: 1049 : m_stream = std::make_unique< State[] >( n );
70 : : // Initialize thread-streams
71 [ + + ][ + - ]: 5416 : for (SeqNumType i=0; i<n; ++i) m_init( &m_stream[i], i );
72 : 1049 : }
73 : :
74 : : //! Uniform RNG: Generate uniform random numbers
75 : : //! \param[in] tid Thread (or more precisely) stream ID
76 : : //! \param[in] num Number of RNGs to generate
77 : : //! \param[in,out] r Pointer to memory to write the random numbers to
78 : 907324433 : void uniform( int tid, ncomp_t num, double* r ) const {
79 [ + + ]: 1816848778 : for (ncomp_t i=0; i<num; ++i)
80 : 909524345 : r[i] = static_cast<double>(
81 : 909524345 : Generate( &m_stream[ static_cast<std::size_t>(tid) ] ) )
82 : 909524345 : / 4294967296.0;
83 : 907324433 : }
84 : :
85 : : //! Gaussian RNG: Generate Gaussian random numbers
86 : : //! \param[in] tid Thread (or more precisely stream) ID
87 : : //! \param[in] num Number of RNGs to generate
88 : : //! \param[in,out] r Pointer to memory to write the random numbers to
89 : : //! \details Generating Gaussian random numbers is implemented via an
90 : : //! adaptor, modeling std::UniformRandomNumberGenerator, outsourcing the
91 : : //! transformation of uniform random numbers to Gaussian ones, to the
92 : : //! standard library. The adaptor is instantiated here because a standard
93 : : //! distribution, such as e.g., std::normal_distribution, generates
94 : : //! numbers using operator() with no arguments, thus the RNG state and the
95 : : //! thread ID (this latter only known here) must be stored in the adaptor
96 : : //! functor's state. Even though creating the adaptor seems like a
97 : : //! potentially costly operation for every call, using the standard
98 : : //! library implementation is still faster than a hand-coded
99 : : //! implementation of the Box-Muller algorithm. Note that libc++ uses a
100 : : //! cache, as Box-Muller, implemented using the polar algorithm generates
101 : : //! 2 Gaussian numbers for each pair of uniform ones, caching every 2nd.
102 : 392000616 : void gaussian( int tid, ncomp_t num, double* r ) const {
103 : 392000616 : Adaptor generator( m_stream, tid );
104 [ + - ]: 392000616 : std::normal_distribution<> gauss_dist( 0.0, 1.0 );
105 [ + + ][ + - ]: 1191400616 : for (ncomp_t i=0; i<num; ++i) r[i] = gauss_dist( generator );
106 : 392000616 : }
107 : :
108 : : //! \brief Multi-variate Gaussian RNG: Generate multi-variate Gaussian
109 : : //! random numbers
110 : : //! \param[in] tid Thread (or more precisely stream) ID
111 : : //! \param[in] num Number of RNGs to generate
112 : : //! \param[in] d Dimension d ( d ≥ 1) of output random vectors
113 : : //! \param[in] mean Mean vector of dimension d
114 : : //! \param[in] cov Lower triangle of covariance matrix, stored as a vector
115 : : //! of length d(d+1)/2
116 : : //! \param[in,out] r Pointer to memory to write the random numbers to
117 : : //! \warning Not implemented!
118 : 0 : void gaussianmv( [[maybe_unused]] int tid,
119 : : [[maybe_unused]] ncomp_t num,
120 : : [[maybe_unused]] ncomp_t d,
121 : : [[maybe_unused]] const double* const mean,
122 : : [[maybe_unused]] const double* const cov,
123 : : [[maybe_unused]] double* r ) const
124 : : {
125 : 0 : }
126 : :
127 : : //! Beta RNG: Generate beta random numbers
128 : : //! \param[in] tid Thread (or more precisely stream) ID
129 : : //! \param[in] num Number of RNGs to generate
130 : : //! \param[in] p First beta shape parameter
131 : : //! \param[in] q Second beta shape parameter
132 : : //! \param[in] a Beta displacement parameter
133 : : //! \param[in] b Beta scale factor
134 : : //! \param[in,out] r Pointer to memory to write the random numbers to
135 : : //! \details Generating beta-distributed random numbers is implemented via
136 : : //! an adaptor, modeling boost::UniformRandomNumberGenerator, outsourcing
137 : : //! the transformation of uniform random numbers to beta-distributed ones,
138 : : //! to boost::random. The adaptor is instantiated here because a boost
139 : : //! random number distribution, such as e.g.,
140 : : //! boost::random::beta_distribution, generates numbers using operator()
141 : : //! with no arguments, thus the RNG state and the thread ID (this latter
142 : : //! only known here) must be stored in the adaptor functor's state.
143 : 176 : void beta( int tid, ncomp_t num, double p, double q, double a, double b,
144 : : double* r ) const {
145 : 176 : Adaptor generator( m_stream, tid );
146 : 176 : boost::random::beta_distribution<> beta_dist( p, q );
147 [ + + ][ + - ]: 4400176 : for (ncomp_t i=0; i<num; ++i) r[i] = beta_dist( generator ) * b + a;
148 : 176 : }
149 : :
150 : : //! Gamma RNG: Generate gamma random numbers
151 : : //! \param[in] tid Thread (or more precisely stream) ID
152 : : //! \param[in] num Number of RNGs to generate
153 : : //! \param[in] a Gamma shape parameter
154 : : //! \param[in] b Gamma scale factor
155 : : //! \param[in,out] r Pointer to memory to write the random numbers to
156 : : //! \details Generating gamma-distributed random numbers is implemented via
157 : : //! an adaptor, modeling boost::UniformRandomNumberGenerator, outsourcing
158 : : //! the transformation of uniform random numbers to gamma-distributed
159 : : //! ones, to boost::random. The adaptor is instantiated here because a
160 : : //! boost random number distribution, such as e.g.,
161 : : //! boost::random::gamma_distribution, generates numbers using operator()
162 : : //! with no arguments, thus the RNG state and the thread ID (this latter
163 : : //! only known here) must be stored in the adaptor functor's state.
164 : 0 : void gamma( int tid, ncomp_t num, double a, double b, double* r ) const {
165 : 0 : Adaptor generator( m_stream, tid );
166 [ - - ]: 0 : boost::random::gamma_distribution<> gamma_dist( a, b );
167 [ - - ][ - - ]: 0 : for (ncomp_t i=0; i<num; ++i) r[i] = gamma_dist( generator );
168 : 0 : }
169 : :
170 : : //! Copy assignment
171 : 110 : RNGSSE& operator=( const RNGSSE& x ) {
172 : 110 : m_nthreads = x.m_nthreads;
173 : 110 : m_init = x.m_init;
174 : 110 : m_stream = std::make_unique< State[] >( x.m_nthreads );
175 [ + + ]: 550 : for (SeqNumType i=0; i<x.m_nthreads; ++i) m_init( &m_stream[i], i );
176 : 110 : return *this;
177 : : }
178 : :
179 : : //! Copy constructor: in terms of copy assignment
180 [ + - ]: 99 : RNGSSE( const RNGSSE& x ) { operator=(x); }
181 : :
182 : : //! Move assignment
183 : 2846 : RNGSSE& operator=( RNGSSE&& x ) {
184 : 2846 : m_nthreads = x.m_nthreads;
185 : 2846 : m_init = x.m_init;
186 : 2846 : m_stream = std::make_unique< State[] >( x.m_nthreads );
187 [ + + ]: 14572 : for (SeqNumType i=0; i<x.m_nthreads; ++i) {
188 : 11726 : m_stream[i] = x.m_stream[i];
189 : 11726 : std::memset( &x.m_stream[i], 0, sizeof(x.m_stream[i]) );
190 : : }
191 : 2846 : x.m_nthreads = 0;
192 : 2846 : x.m_init = nullptr;
193 : 2846 : x.m_stream.reset( nullptr );
194 : 2846 : return *this;
195 : : }
196 : :
197 : : //! Move constructor: in terms of move assignment
198 : 2846 : RNGSSE( RNGSSE&& x ) :
199 : : m_nthreads( 0 ),
200 : : m_init( nullptr ),
201 : 2846 : m_stream( nullptr )
202 [ + - ]: 2846 : { *this = std::move( x ); }
203 : :
204 : : //! Accessor to the number of threads we operate on
205 : 220 : SeqNumType nthreads() const noexcept { return m_nthreads; }
206 : :
207 : : private:
208 : : SeqNumType m_nthreads; //!< Number of threads
209 : : InitFn m_init; //!< Sequence length initializer
210 : : std::unique_ptr< State[] > m_stream; //!< Random number stream for threads
211 : : };
212 : :
213 : : } // tk::
214 : :
215 : : #endif // RNGSSE_h
|
