mxnet
random.h
Go to the documentation of this file.
1 
8 #ifndef MSHADOW_RANDOM_H_
9 #define MSHADOW_RANDOM_H_
10 
11 #include <cstdlib>
12 #include <algorithm>
13 #include <random>
14 #include "./base.h"
15 #include "./tensor.h"
16 #include "./tensor_container.h"
17 
18 #if MSHADOW_IN_CXX11
19 #include <random> // use cxx11 random by default
20 #endif
21 
22 #if _MSC_VER
23 #define rand_r(x) rand()
24 #endif
25 
26 
27 namespace mshadow {
33 template<typename Device, typename DType MSHADOW_DEFAULT_DTYPE>
34 class Random {};
35 
37 template<typename DType>
38 class Random<cpu, DType> {
39  public:
44  explicit Random(int seed) {
45  this->Seed(seed);
46  buffer_.Resize(Shape1(kRandBufferSize));
47  }
48  ~Random(void) {
49  }
54  inline void Seed(int seed) {
55 #if MSHADOW_IN_CXX11
56  rnd_engine_.seed(seed);
57 #endif
58  this->rseed_ = static_cast<unsigned>(seed);
59  }
64  inline unsigned GetSeed() const {
65  return rseed_;
66  }
71  inline void set_stream(Stream<cpu> *stream) {
72  }
73 
74 // These samplers are only avail in C++11.
75 #if MSHADOW_IN_CXX11
76 
81  inline unsigned GetRandInt() {
82  return rnd_engine_();
83  }
84 
88  inline void GetRandInt(const Tensor<cpu, 1, unsigned>& dst) {
89  std::generate_n(dst.dptr_, dst.size(0), [&](){ return rnd_engine_(); });
90  }
91 
98  template<int dim, class Sampler>
99  inline void SampleDistribution(Tensor<cpu, dim, DType> *dst, Sampler sampler) {
100  if (dst->CheckContiguous()) {
101  std::generate_n(dst->dptr_, dst->shape_.Size(), sampler);
102  } else {
103  Tensor<cpu, 2, DType> mat = dst->FlatTo2D();
104  for (index_t i = 0; i < mat.size(0); ++i) {
105  std::generate_n(mat[i].dptr_, mat.size(1), sampler);
106  }
107  }
108  }
109 
117  template<int dim, typename PType>
118  inline void SampleUniform(Tensor<cpu, dim, DType> *dst,
119  PType a = 0.0f , PType b = 1.0f ) {
120  // Ensure that half_t is handled correctly.
121  typedef typename std::conditional<std::is_floating_point<DType>::value,
122  DType, double>::type FType;
123  typedef typename std::conditional<std::is_integral<DType>::value,
124  std::uniform_int_distribution<DType>,
125  std::uniform_real_distribution<FType>>::type GType;
126  GType dist_uniform(a, b);
127  SampleDistribution(dst, [&](){ return dist_uniform(rnd_engine_);});
128  }
129 
137  template<int dim, typename PType>
138  inline void SampleGaussian(Tensor<cpu, dim, DType> *dst,
139  PType mu = 0.0f, PType sigma = 1.0f ) {
140  if (sigma <= 0) {
141  *dst = mu; return;
142  }
143  typedef typename std::conditional<std::is_floating_point<DType>::value,
144  DType, double>::type GType;
145  std::normal_distribution<GType> dist_normal(mu, sigma);
146  SampleDistribution(dst, [&](){ return dist_normal(rnd_engine_);});
147  }
148 
156  template<int dim, typename PType>
157  inline void SampleGamma(Tensor<cpu, dim, DType> *dst,
158  PType alpha, PType beta) {
159  typedef typename std::conditional<std::is_floating_point<DType>::value,
160  DType, double>::type GType;
161  std::gamma_distribution<GType> dist_gamma(alpha, beta);
162  SampleDistribution(dst, [&](){ return dist_gamma(rnd_engine_);});
163  }
164 
171  template<int dim, typename PType>
172  inline void SampleExponential(Tensor<cpu, dim, DType> *dst, PType lambda ) {
173  typedef typename std::conditional<std::is_floating_point<DType>::value,
174  DType, double>::type GType;
175  std::exponential_distribution<GType> dist_exp(lambda);
176  SampleDistribution(dst, [&](){ return dist_exp(rnd_engine_);});
177  }
178 
185  template<int dim, typename PType>
186  inline void SamplePoisson(Tensor<cpu, dim, DType> *dst, PType lambda) {
187  typedef typename std::conditional<std::is_integral<DType>::value, DType, int>::type GType;
188  std::poisson_distribution<GType> dist_poisson(lambda);
189  SampleDistribution(dst, [&](){ return static_cast<DType>(dist_poisson(rnd_engine_));});
190  }
191 
199  template<int dim, typename PType1, typename PType2>
200  inline void SampleNegativeBinomial(Tensor<cpu, dim, DType> *dst, PType1 k, PType2 p) {
201  typedef typename std::conditional<std::is_integral<DType>::value, DType, int>::type GType;
202  std::negative_binomial_distribution<GType> dist_negbinomial(k, p);
203  SampleDistribution(dst, [&](){ return static_cast<DType>(dist_negbinomial(rnd_engine_));});
204  }
205 
214  template<int dim, typename PType>
215  inline void SampleGeneralizedNegativeBinomial(Tensor<cpu, dim, DType> *dst,
216  PType mu, PType alpha) {
217  if (alpha == PType(0)) {
218  SamplePoisson(dst, mu); // limit of Poisson
219  } else {
220  PType r(PType(1) / alpha);
221  PType beta = mu * alpha;
222  std::gamma_distribution<> dist_gamma(r, beta);
223  typedef typename std::conditional<std::is_integral<DType>::value, DType, int>::type GType;
224  SampleDistribution(dst,
225  [&](){ std::poisson_distribution<GType> dist_poisson(dist_gamma(rnd_engine_));
226  return static_cast<DType>(dist_poisson(rnd_engine_));});
227  }
228  }
229 #endif
230 
242  template<int dim>
243  inline expr::ReshapeExp<Tensor<cpu, 1, DType>, DType, dim, 1>
245  buffer_.Resize(Shape1(shape.Size()));
246  this->SampleGaussian(&buffer_, 0.0f, 1.0f);
247  return expr::reshape(buffer_, shape);
248  }
260  template<int dim>
261  inline expr::ReshapeExp<Tensor<cpu, 1, DType>, DType, dim, 1>
263  buffer_.Resize(Shape1(shape.Size()));
264  this->SampleUniform(&buffer_, 0.0f, 1.0f);
265  return expr::reshape(buffer_, shape);
266  }
267 
268  std::mt19937 &GetRndEngine() {
269  return rnd_engine_;
270  }
271 
272  private:
273 #if MSHADOW_IN_CXX11
274 
275  std::mt19937 rnd_engine_;
277  unsigned rseed_;
278 
279 #else
280 
282  unsigned rseed_;
283  // functions
284  template<int dim>
285  inline void SampleUniform(Tensor<cpu, dim, DType> *dst,
286  DType a = 0.0f, DType b = 1.0f) {
287  if (dst->CheckContiguous()) {
288  this->GenUniform(dst->dptr_, dst->shape_.Size(), a, b);
289  } else {
290  Tensor<cpu, 2, DType> mat = dst->FlatTo2D();
291  for (index_t i = 0; i < mat.size(0); ++i) {
292  this->GenUniform(mat[i].dptr_, mat.size(1), a, b);
293  }
294  }
295  }
296  template<int dim>
297  inline void SampleGaussian(Tensor<cpu, dim, DType> *dst,
298  DType mu = 0.0f, DType sigma = 1.0f) {
299  if (sigma <= 0.0f) {
300  *dst = mu; return;
301  }
302  if (dst->CheckContiguous()) {
303  this->GenGaussian(dst->dptr_, dst->shape_.Size(), mu, sigma);
304  } else {
305  Tensor<cpu, 2, DType> mat = dst->FlatTo2D();
306  for (index_t i = 0; i < mat.size(0); ++i) {
307  this->GenGaussian(mat[i].dptr_, mat.size(1), mu, sigma);
308  }
309  }
310  }
311  inline void GenUniform(float *dptr, index_t size, float a, float b) {
312  for (index_t j = 0; j < size; ++j) {
313  dptr[j] = static_cast<float>(RandNext()) * (b - a) + a;
314  }
315  }
316  inline void GenUniform(double *dptr, index_t size, double a, double b) {
317  for (index_t j = 0; j < size; ++j) {
318  dptr[j] = static_cast<double>(RandNext()) * (b - a) + a;
319  }
320  }
321  inline void GenGaussian(float *dptr, index_t size, float mu, float sigma) {
322  this->GenGaussianX(dptr, size, mu, sigma);
323  }
324  inline void GenGaussian(double *dptr, index_t size, double mu, double sigma) {
325  this->GenGaussianX(dptr, size, mu, sigma);
326  }
327  inline void GenGaussianX(DType *dptr, index_t size, DType mu, DType sigma) {
328  DType g1 = 0.0f, g2 = 0.0f;
329  for (index_t j = 0; j < size; ++j) {
330  if ((j & 1) == 0) {
331  this->SampleNormal2D(&g1, &g2);
332  dptr[j] = mu + g1 * sigma;
333  } else {
334  dptr[j] = mu + g2 * sigma;
335  }
336  }
337  }
339  inline DType RandNext(void) {
340  return static_cast<DType>(rand_r(&rseed_)) /
341  (static_cast<DType>(RAND_MAX) + 1.0f);
342  }
344  inline DType RandNext2(void) {
345  return (static_cast<DType>(rand_r(&rseed_)) + 1.0f) /
346  (static_cast<DType>(RAND_MAX) + 2.0f);
347  }
353  inline void SampleNormal2D(DType *xx_, DType *yy_) {
354  DType &xx = *xx_, &yy = *yy_;
355  DType x, y, s;
356  do {
357  x = 2.0f * RandNext2() - 1.0f;
358  y = 2.0f * RandNext2() - 1.0f;
359  s = x * x + y * y;
360  } while (s >= 1.0f || s == 0.0f);
361  DType t = std::sqrt(-2.0f * std::log(s) / s);
362  xx = x * t; yy = y * t;
363  }
364 #endif
365 
367 }; // class Random<cpu, DType>
368 
369 // only allow GPU PRNG when cuda is enabled
370 #if MSHADOW_USE_CUDA
371 
372 template<typename DType>
373 class Random<gpu, DType> {
374  public:
379  explicit Random(int seed) : gen_(NULL) {
380  this->Seed(seed);
381  buffer_.Resize(Shape1(kRandBufferSize));
382  }
384  DeleteGenerator();
385  }
390  inline void set_stream(Stream<gpu> *stream) {
391  curandStatus_t status;
392  status = curandSetStream(gen_, Stream<gpu>::GetStream(stream));
393 
394  CHECK_EQ(status, CURAND_STATUS_SUCCESS) << "set_stream CURAND failed";
395  }
400  inline void Seed(int seed) {
401  // Create a new rng, either initially or if the RNG type can't reset its offset.
402  if (gen_ == NULL || (curandSetGeneratorOffset(gen_, 0ULL) != CURAND_STATUS_SUCCESS))
403  CreateGenerator();
404  // Now set the seed.
405  curandStatus_t status;
406  status = curandSetPseudoRandomGeneratorSeed(gen_, static_cast<uint64_t>(seed));
407  CHECK_EQ(status, CURAND_STATUS_SUCCESS) << "Set CURAND seed failed.";
408  }
412  inline void GetRandInt(const Tensor<gpu, 1, unsigned>& dst) {
413  curandStatus_t status = curandGenerate(gen_, dst.dptr_, dst.size(0));
414  CHECK_EQ(status, CURAND_STATUS_SUCCESS) << "CURAND Gen rand ints failed.";
415  }
423  template<int dim>
424  inline void SampleUniform(Tensor<gpu, dim, DType> *dst,
425  DType a = 0.0f, DType b = 1.0f);
426 
434  template<int dim>
435  inline void SampleGaussian(Tensor<gpu, dim, DType> *dst,
436  DType mu = 0.0f, DType sigma = 1.0f);
450  template<int dim>
451  inline expr::ReshapeExp<Tensor<gpu, 1, DType>, DType, dim, 1>
452  gaussian(Shape<dim> shape, DType mu = 0.0f, DType sigma = 1.0f);
464  template<int dim>
465  inline expr::ReshapeExp<Tensor<gpu, 1, DType>, DType, dim, 1>
466  uniform(Shape<dim> shape);
467 
468  private:
469  inline void GenGaussian(float *dptr, size_t size, float mu, float sigma) {
470  curandStatus_t status;
471  status = curandGenerateNormal(gen_, dptr, size, mu, sigma);
472  CHECK_EQ(status, CURAND_STATUS_SUCCESS) << "CURAND Gen Normal float failed."
473  << " size = " << size
474  << ",mu = " << mu
475  << ",sigma = " << sigma;
476  }
477  inline void GenGaussian(double *dptr, size_t size, double mu, double sigma) {
478  curandStatus_t status;
479  status = curandGenerateNormalDouble(gen_, dptr, size, mu, sigma);
480  CHECK_EQ(status, CURAND_STATUS_SUCCESS) << "CURAND Gen Normal double failed."
481  << " size = " << size
482  << ",mu = " << mu
483  << ",sigma = " << sigma;
484  }
485  inline void GenUniform(float *dptr, size_t size) {
486  curandStatus_t status;
487  status = curandGenerateUniform(gen_, dptr, size);
488  CHECK_EQ(status, CURAND_STATUS_SUCCESS) << "CURAND Gen Uniform float failed."
489  << " size = " << size;
490  }
491  inline void GenUniform(double *dptr, size_t size) {
492  curandStatus_t status;
493  status = curandGenerateUniformDouble(gen_, dptr, size);
494  CHECK_EQ(status, CURAND_STATUS_SUCCESS) << "CURAND Gen Uniform double failed."
495  << " size = " << size;
496  }
497  inline void CreateGenerator() {
498  if (gen_ != NULL)
499  DeleteGenerator();
500  curandStatus_t status;
501  status = curandCreateGenerator(&gen_, CURAND_RNG_PSEUDO_DEFAULT);
502  CHECK_EQ(status, CURAND_STATUS_SUCCESS) << "Cannot create CURAND Generator";
503  }
504  inline void DeleteGenerator() {
505  if (gen_ != NULL) {
506  curandStatus_t status;
507  status = curandDestroyGenerator(gen_);
508  CHECK_EQ(status, CURAND_STATUS_SUCCESS) << "Destory CURAND Gen failed";
509  gen_ = NULL;
510  }
511  }
513  curandGenerator_t gen_;
516 }; // class Random<gpu, DType>
517 #endif // MSHADOW_USE_CUDA
518 
519 #ifdef __CUDACC__
520 // implementations that depends on cuda kernels
521 template<typename DType>
522 template<int dim>
524  Tensor<gpu, dim, DType> *dst, DType a, DType b) {
525  if (a == 0.0f && b == 1.0f) {
526  if (dst->CheckContiguous()) {
527  this->GenUniform(dst->dptr_, dst->shape_.Size());
528  } else {
529  *dst = this->uniform(dst->shape_);
530  }
531  } else {
532  *dst = this->uniform(dst->shape_) * (b - a) + a;
533  }
534 }
535 template<typename DType>
536 template<int dim>
538  Tensor<gpu, dim, DType> *dst, DType mu, DType sigma) {
539  // We need to check whether the shape size is even since CuRand supports only normal distribution
540  // generation of even number of elements.
541  if (dst->CheckContiguous() && (dst->shape_.Size() % 2 == 0)) {
542  this->GenGaussian(dst->dptr_, dst->shape_.Size(), mu, sigma);
543  } else {
544  *dst = this->gaussian(dst->shape_, mu, sigma);
545  }
546 }
547 
548 template<typename DType>
549 template<int dim>
550 inline expr::ReshapeExp<Tensor<gpu, 1, DType>, DType, dim, 1>
551 Random<gpu, DType>::gaussian(Shape<dim> shape, DType mu, DType sigma) {
552  size_t aligned_sz = ((shape.Size() + 1UL) >> 1) << 1;
553  // allocate alligned size
554  buffer_.Resize(Shape1(aligned_sz));
555  buffer_.Resize(Shape1(shape.Size()));
556  this->GenGaussian(buffer_.dptr_, aligned_sz, mu, sigma);
557  return expr::reshape(buffer_, shape);
558 }
559 
560 template<typename DType>
561 template<int dim>
562 inline expr::ReshapeExp<Tensor<gpu, 1, DType>, DType, dim, 1>
564  buffer_.Resize(Shape1(shape.Size()));
565  this->GenUniform(buffer_.dptr_, buffer_.size(0));
566  return expr::reshape(buffer_, shape);
567 }
568 #endif // __CUDACC__
569 } // namespace mshadow
570 #endif // MSHADOW_RANDOM_H_
expr::ReshapeExp< Tensor< cpu, 1, DType >, DType, dim, 1 > gaussian(Shape< dim > shape)
return a temporal expression storing standard gaussian random variables the temporal tensor is only v...
Definition: random.h:244
random number generator
Definition: random.h:34
void Seed(int seed)
seed random number generator using this seed
Definition: random.h:400
MSHADOW_XINLINE index_t Size(void) const
Definition: tensor.h:126
DType * dptr_
pointer to the data
Definition: tensor.h:416
unsigned GetSeed() const
get random seed used in random generator
Definition: random.h:64
void Seed(int seed)
seed random number generator using this seed
Definition: random.h:54
Definition: stream_gpu-inl.h:19
~Random(void) MSHADOW_THROW_EXCEPTION
Definition: random.h:383
Shape< dimension > shape_
shape of the tensor
Definition: tensor.h:418
~Random(void)
Definition: random.h:48
expr::ReshapeExp< Tensor< cpu, 1, DType >, DType, dim, 1 > uniform(Shape< dim > shape)
return a temporal expression storing standard uniform [0,1) the temporal tensor is only valid before ...
Definition: random.h:262
device name CPU
Definition: tensor.h:21
device name GPU
Definition: tensor.h:28
MSHADOW_XINLINE Tensor< Device, 2, DType > FlatTo2D(void) const
flatten the tensor to 2 dimension, collapse the higher dimensions together
Definition: tensor.h:501
const unsigned kRandBufferSize
buffer size for each random number generator
Definition: base.h:284
void SampleExponential(real_t lambda, NDArray *out)
Sample exponential distribution for each elements of out.
int32_t index_t
type that will be used for index
Definition: base.h:291
ReshapeExp< SrcExp, DType, dimdst, ExpInfo< SrcExp >::kDim > reshape(const Exp< SrcExp, DType, etype > &src, Shape< dimdst > oshape)
a expression that reshapes a tensor to another shape
Definition: reshape.h:48
void SampleGaussian(real_t mu, real_t sigma, NDArray *out)
Sample gaussian distribution for each elements of out.
tensor container that does memory allocation and resize like STL, use it to save the lines of FreeSpa...
Definition: tensor_container.h:22
MSHADOW_XINLINE Shape< 1 > Shape1(index_t s0)
construct a one dimension shape, stride will equal s0
Definition: tensor.h:188
void GetRandInt(const Tensor< gpu, 1, unsigned > &dst)
get a set of random integers
Definition: random.h:412
void SampleUniform(real_t begin, real_t end, NDArray *out)
Sample uniform distribution for each elements of out.
MSHADOW_XINLINE bool CheckContiguous(void) const
Definition: tensor.h:473
reshape the content to another shape input: Tensor<Device,dimsrc>: ishape output: Tensor<Device...
Definition: reshape.h:21
std::mt19937 & GetRndEngine()
Definition: random.h:268
Random(int seed)
constructor of random engine
Definition: random.h:379
tensor container that does memory allocation and resize like STL
void set_stream(Stream< cpu > *stream)
set the stream of computation
Definition: random.h:71
Random(int seed)
constructor of random engine
Definition: random.h:44
namespace for mshadow
Definition: base.h:282
#define MSHADOW_THROW_EXCEPTION
Definition: base.h:234
MSHADOW_XINLINE index_t size(int idx) const
return size of i-th dimension, start counting from highest dimension
Definition: tensor.h:487
general tensor
Definition: tensor.h:402
void SamplePoisson(real_t lambda, NDArray *out)
Sample Poisson distribution for each elements of out.
void set_stream(Stream< gpu > *stream)
set the stream of computation
Definition: random.h:390
void SampleGamma(real_t alpha, real_t beta, NDArray *out)
Sample gamma distribution for each elements of out.
computaion stream structure, used for asynchronous computations
Definition: tensor.h:365