TensorRandom.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2016 Benoit Steiner <benoit.steiner.goog@gmail.com>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_CXX11_TENSOR_TENSOR_RANDOM_H
11 #define EIGEN_CXX11_TENSOR_TENSOR_RANDOM_H
12 
13 namespace Eigen {
14 namespace internal {
15 
16 namespace {
17 
18 EIGEN_DEVICE_FUNC uint64_t get_random_seed() {
19 #ifdef __CUDA_ARCH__
20  // We don't support 3d kernels since we currently only use 1 and
21  // 2d kernels.
22  assert(threadIdx.z == 0);
23  return clock64() +
24  blockIdx.x * blockDim.x + threadIdx.x +
25  gridDim.x * blockDim.x * (blockIdx.y * blockDim.y + threadIdx.y);
26 
27 #elif defined _WIN32
28  // Use the current time as a baseline.
29  SYSTEMTIME st;
30  GetSystemTime(&st);
31  int time = st.wSecond + 1000 * st.wMilliseconds;
32  // Mix in a random number to make sure that we get different seeds if
33  // we try to generate seeds faster than the clock resolution.
34  // We need 2 random values since the generator only generate 16 bits at
35  // a time (https://msdn.microsoft.com/en-us/library/398ax69y.aspx)
36  int rnd1 = ::rand();
37  int rnd2 = ::rand();
38  uint64_t rnd = (rnd1 | rnd2 << 16) ^ time;
39  return rnd;
40 
41 #elif defined __APPLE__
42  // Same approach as for win32, except that the random number generator
43  // is better (// https://developer.apple.com/legacy/library/documentation/Darwin/Reference/ManPages/man3/random.3.html#//apple_ref/doc/man/3/random).
44  uint64_t rnd = ::random() ^ mach_absolute_time();
45  return rnd;
46 
47 #else
48  // Augment the current time with pseudo random number generation
49  // to ensure that we get different seeds if we try to generate seeds
50  // faster than the clock resolution.
51  timespec ts;
52  clock_gettime(CLOCK_REALTIME, &ts);
53  uint64_t rnd = ::random() ^ ts.tv_nsec;
54  return rnd;
55 #endif
56 }
57 
58 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE unsigned PCG_XSH_RS_generator(uint64_t* state) {
59  // TODO: Unify with the implementation in the non blocking thread pool.
60  uint64_t current = *state;
61  // Update the internal state
62  *state = current * 6364136223846793005ULL + 0xda3e39cb94b95bdbULL;
63  // Generate the random output (using the PCG-XSH-RS scheme)
64  return static_cast<unsigned>((current ^ (current >> 22)) >> (22 + (current >> 61)));
65 }
66 
67 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE uint64_t PCG_XSH_RS_state(uint64_t seed) {
68  seed = seed ? seed : get_random_seed();
69  return seed * 6364136223846793005ULL + 0xda3e39cb94b95bdbULL;
70 }
71 
72 } // namespace
73 
74 
75 template <typename T> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
77  unsigned rnd = PCG_XSH_RS_generator(state);
78  return static_cast<T>(rnd);
79 }
80 
81 
82 template <> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
83 Eigen::half RandomToTypeUniform<Eigen::half>(uint64_t* state) {
85  // Generate 10 random bits for the mantissa
86  unsigned rnd = PCG_XSH_RS_generator(state);
87  result.x = static_cast<uint16_t>(rnd & 0x3ffu);
88  // Set the exponent
89  result.x |= (static_cast<uint16_t>(15) << 10);
90  // Return the final result
91  return result - Eigen::half(1.0f);
92 }
93 
94 
95 template <> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
97  typedef union {
98  uint32_t raw;
99  float fp;
100  } internal;
101  internal result;
102  // Generate 23 random bits for the mantissa mantissa
103  const unsigned rnd = PCG_XSH_RS_generator(state);
104  result.raw = rnd & 0x7fffffu;
105  // Set the exponent
106  result.raw |= (static_cast<uint32_t>(127) << 23);
107  // Return the final result
108  return result.fp - 1.0f;
109 }
110 
111 template <> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
113  typedef union {
114  uint64_t raw;
115  double dp;
116  } internal;
117  internal result;
118  result.raw = 0;
119  // Generate 52 random bits for the mantissa
120  // First generate the upper 20 bits
121  unsigned rnd1 = PCG_XSH_RS_generator(state) & 0xfffffu;
122  // The generate the lower 32 bits
123  unsigned rnd2 = PCG_XSH_RS_generator(state);
124  result.raw = (static_cast<uint64_t>(rnd1) << 32) | rnd2;
125  // Set the exponent
126  result.raw |= (static_cast<uint64_t>(1023) << 52);
127  // Return the final result
128  return result.dp - 1.0;
129 }
130 
131 template <> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
132 std::complex<float> RandomToTypeUniform<std::complex<float> >(uint64_t* state) {
133  return std::complex<float>(RandomToTypeUniform<float>(state),
135 }
136 template <> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
137 std::complex<double> RandomToTypeUniform<std::complex<double> >(uint64_t* state) {
138  return std::complex<double>(RandomToTypeUniform<double>(state),
140 }
141 
142 template <typename T> class UniformRandomGenerator {
143  public:
144  static const bool PacketAccess = true;
145 
146  // Uses the given "seed" if non-zero, otherwise uses a random seed.
148  uint64_t seed = 0) {
149  m_state = PCG_XSH_RS_state(seed);
150  }
152  const UniformRandomGenerator& other) {
153  m_state = other.m_state;
154  }
155 
156  template<typename Index> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
157  T operator()(Index i) const {
158  uint64_t local_state = m_state + i;
159  T result = RandomToTypeUniform<T>(&local_state);
160  m_state = local_state;
161  return result;
162  }
163 
164  template<typename Packet, typename Index> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
166  const int packetSize = internal::unpacket_traits<Packet>::size;
167  EIGEN_ALIGN_MAX T values[packetSize];
168  uint64_t local_state = m_state + i;
169  for (int j = 0; j < packetSize; ++j) {
170  values[j] = RandomToTypeUniform<T>(&local_state);
171  }
172  m_state = local_state;
173  return internal::pload<Packet>(values);
174  }
175 
176  private:
177  mutable uint64_t m_state;
178 };
179 
180 template <typename Scalar>
182  enum {
183  // Rough estimate for floating point, multiplied by ceil(sizeof(T) / sizeof(float)).
185  ((sizeof(Scalar) + sizeof(float) - 1) / sizeof(float)),
187  };
188 };
189 
190 
191 
192 template <typename T> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
194  // Use the ratio of uniform method to generate numbers following a normal
195  // distribution. See for example Numerical Recipes chapter 7.3.9 for the
196  // details.
197  T u, v, q;
198  do {
199  u = RandomToTypeUniform<T>(state);
200  v = T(1.7156) * (RandomToTypeUniform<T>(state) - T(0.5));
201  const T x = u - T(0.449871);
202  const T y = numext::abs(v) + T(0.386595);
203  q = x*x + y * (T(0.196)*y - T(0.25472)*x);
204  } while (q > T(0.27597) &&
205  (q > T(0.27846) || v*v > T(-4) * numext::log(u) * u*u));
206 
207  return v/u;
208 }
209 
210 template <> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
211 std::complex<float> RandomToTypeNormal<std::complex<float> >(uint64_t* state) {
212  return std::complex<float>(RandomToTypeNormal<float>(state),
213  RandomToTypeNormal<float>(state));
214 }
215 template <> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
216 std::complex<double> RandomToTypeNormal<std::complex<double> >(uint64_t* state) {
217  return std::complex<double>(RandomToTypeNormal<double>(state),
218  RandomToTypeNormal<double>(state));
219 }
220 
221 
222 template <typename T> class NormalRandomGenerator {
223  public:
224  static const bool PacketAccess = true;
225 
226  // Uses the given "seed" if non-zero, otherwise uses a random seed.
227  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NormalRandomGenerator(uint64_t seed = 0) {
228  m_state = PCG_XSH_RS_state(seed);
229  }
231  const NormalRandomGenerator& other) {
232  m_state = other.m_state;
233  }
234 
235  template<typename Index> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
236  T operator()(Index i) const {
237  uint64_t local_state = m_state + i;
238  T result = RandomToTypeNormal<T>(&local_state);
239  m_state = local_state;
240  return result;
241  }
242 
243  template<typename Packet, typename Index> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
245  const int packetSize = internal::unpacket_traits<Packet>::size;
246  EIGEN_ALIGN_MAX T values[packetSize];
247  uint64_t local_state = m_state + i;
248  for (int j = 0; j < packetSize; ++j) {
249  values[j] = RandomToTypeNormal<T>(&local_state);
250  }
251  m_state = local_state;
252  return internal::pload<Packet>(values);
253  }
254 
255  private:
256  mutable uint64_t m_state;
257 };
258 
259 
260 template <typename Scalar>
262  enum {
263  // On average, we need to generate about 3 random numbers
264  // 15 mul, 8 add, 1.5 logs
267  3 * functor_traits<scalar_log_op<Scalar> >::Cost / 2,
269  };
270 };
271 
272 
273 } // end namespace internal
274 } // end namespace Eigen
275 
276 #endif // EIGEN_CXX11_TENSOR_TENSOR_RANDOM_H
SCALAR Scalar
Definition: bench_gemm.cpp:33
#define EIGEN_STRONG_INLINE
Definition: Macros.h:494
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NormalRandomGenerator(const NormalRandomGenerator &other)
Definition: TensorRandom.h:230
ArrayXcf v
Definition: Cwise_arg.cpp:1
dim3 threadIdx
Definition: cuda_common.h:11
leaf::MyValues values
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
unsigned short uint16_t
Definition: ms_stdint.h:84
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(Index i) const
Definition: TensorRandom.h:165
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE internal::enable_if< NumTraits< T >::IsSigned||NumTraits< T >::IsComplex, typename NumTraits< T >::Real >::type abs(const T &x)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T RandomToTypeUniform(uint64_t *state)
Definition: TensorRandom.h:76
Values result
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE UniformRandomGenerator(uint64_t seed=0)
Definition: TensorRandom.h:147
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NormalRandomGenerator(uint64_t seed=0)
Definition: TensorRandom.h:227
unsigned int uint32_t
Definition: ms_stdint.h:85
unsigned __int64 uint64_t
Definition: ms_stdint.h:95
Eigen::Triplet< double > T
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
#define time
EIGEN_DEVICE_FUNC const Scalar & q
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE UniformRandomGenerator(const UniformRandomGenerator &other)
Definition: TensorRandom.h:151
A small structure to hold a non zero as a triplet (i,j,value).
Definition: SparseUtil.h:154
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T operator()(Index i) const
Definition: TensorRandom.h:236
dim3 blockIdx
Definition: cuda_common.h:11
#define EIGEN_ALIGN_MAX
Definition: Macros.h:757
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float RandomToTypeUniform< float >(uint64_t *state)
Definition: TensorRandom.h:96
const mpreal random(unsigned int seed=0)
Definition: mpreal.h:2614
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(Index i) const
Definition: TensorRandom.h:244
dim3 blockDim
Definition: cuda_common.h:11
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double RandomToTypeUniform< double >(uint64_t *state)
Definition: TensorRandom.h:112
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T log(const T &x)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T RandomToTypeNormal(uint64_t *state)
Definition: TensorRandom.h:193
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T operator()(Index i) const
Definition: TensorRandom.h:157
std::ptrdiff_t j
unsigned short x
Definition: Half.h:57


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:45:40