Main MRPT website > C++ reference
MRPT logo
RandomGenerators.h
Go to the documentation of this file.
1 /* +---------------------------------------------------------------------------+
2  | Mobile Robot Programming Toolkit (MRPT) |
3  | http://www.mrpt.org/ |
4  | |
5  | Copyright (c) 2005-2014, Individual contributors, see AUTHORS file |
6  | See: http://www.mrpt.org/Authors - All rights reserved. |
7  | Released under BSD License. See details in http://www.mrpt.org/License |
8  +---------------------------------------------------------------------------+ */
9 #ifndef RandomGenerator_H
10 #define RandomGenerator_H
11 
12 #include <mrpt/utils/utils_defs.h>
14 
15 namespace mrpt
16 {
17  namespace math { template <class T> class CMatrixTemplateNumeric; } // Frwd. decl.
18 
19  /** A namespace of pseudo-random numbers genrators of diferent distributions. The central class in this namespace is mrpt::random::CRandomGenerator
20  * \ingroup mrpt_base_grp
21  */
22  namespace random
23  {
24  using namespace mrpt::math;
25 
26  /** A thred-safe pseudo random number generator, based on an internal MT19937 randomness generator.
27  * The base algorithm for randomness is platform-independent. See http://en.wikipedia.org/wiki/Mersenne_twister
28  *
29  * For real thread-safety, each thread must create and use its own instance of this class.
30  *
31  * Single-thread programs can use the static object mrpt::random::randomGenerator
32  * \ingroup mrpt_base_grp
33  */
35  {
36  protected:
37  /** Data used internally by the MT19937 PRNG algorithm. */
39  {
40  TMT19937_data() : index(0), seed_initialized(false)
41  {}
42  uint32_t MT[624];
43  uint32_t index;
45  } m_MT19937_data;
46 
49 
50  void MT19937_generateNumbers();
51  void MT19937_initializeGenerator(const uint32_t &seed);
52 
53  public:
54 
55  /** @name Initialization
56  @{ */
57 
58  /** Default constructor: initialize random seed based on current time */
59  CRandomGenerator() : m_MT19937_data(),m_std_gauss_set(false) { randomize(); }
60 
61  /** Constructor for providing a custom random seed to initialize the PRNG */
62  CRandomGenerator(const uint32_t seed) : m_MT19937_data() { randomize(seed); }
63 
64  void randomize(const uint32_t seed); //!< Initialize the PRNG from the given random seed
65  void randomize(); //!< Randomize the generators, based on current time
66 
67  /** @} */
68 
69  /** @name Uniform pdf
70  @{ */
71 
72  /** Generate a uniformly distributed pseudo-random number using the MT19937 algorithm, in the whole range of 32-bit integers.
73  * See: http://en.wikipedia.org/wiki/Mersenne_twister */
74  uint32_t drawUniform32bit();
75 
76  /** Returns a uniformly distributed pseudo-random number by joining two 32bit numbers from \a drawUniform32bit() */
77  uint64_t drawUniform64bit();
78 
79  /** You can call this overloaded method with either 32 or 64bit unsigned ints for the sake of general coding. */
80  void drawUniformUnsignedInt(uint32_t &ret_number) { ret_number=drawUniform32bit(); }
81  void drawUniformUnsignedInt(uint64_t &ret_number) { ret_number=drawUniform64bit(); }
82 
83  /** Return a uniform unsigned integer in the range [min_val,max_val] (both inclusive) */
84  template <typename T, typename U,typename V>
85  void drawUniformUnsignedIntRange(T &ret_number,const U min_val,const V max_val)
86  {
87  const T range = max_val-min_val+1;
88  T rnd;
89  drawUniformUnsignedInt(rnd);
90  ret_number=min_val+ (rnd%range);
91  }
92 
93  /** Generate a uniformly distributed pseudo-random number using the MT19937 algorithm, scaled to the selected range. */
94  double drawUniform( const double Min, const double Max) {
95  return Min + (Max-Min)* drawUniform32bit() * 2.3283064370807973754314699618685e-10; // 0xFFFFFFFF ^ -1
96  }
97 
98  /** Fills the given matrix with independent, uniformly distributed samples.
99  * Matrix classes can be CMatrixTemplateNumeric or CMatrixFixedNumeric
100  * \sa drawUniform
101  */
102  template <class MAT>
104  MAT &matrix,
105  const double unif_min = 0,
106  const double unif_max = 1 )
107  {
108  for (size_t r=0;r<matrix.getRowCount();r++)
109  for (size_t c=0;c<matrix.getColCount();c++)
110  matrix.get_unsafe(r,c) = static_cast<typename MAT::Scalar>( drawUniform(unif_min,unif_max) );
111  }
112 
113  /** Fills the given vector with independent, uniformly distributed samples.
114  * \sa drawUniform
115  */
116  template <class VEC>
118  VEC & v,
119  const double unif_min = 0,
120  const double unif_max = 1 )
121  {
122  const size_t N = v.size();
123  for (size_t c=0;c<N;c++)
124  v[c] = static_cast<typename mrpt::math::ContainerType<VEC>::element_t>( drawUniform(unif_min,unif_max) );
125  }
126 
127  /** @} */
128 
129  /** @name Normal/Gaussian pdf
130  @{ */
131 
132  /** Generate a normalized (mean=0, std=1) normally distributed sample.
133  * \param likelihood If desired, pass a pointer to a double which will receive the likelihood of the given sample to have been obtained, that is, the value of the normal pdf at the sample value.
134  */
135  double drawGaussian1D_normalized( double *likelihood = NULL);
136 
137  /** Generate a normally distributed pseudo-random number.
138  * \param mean The mean value of desired normal distribution
139  * \param std The standard deviation value of desired normal distribution
140  */
141  double drawGaussian1D( const double mean, const double std ) {
142  return mean+std*drawGaussian1D_normalized();
143  }
144 
145  /** Fills the given matrix with independent, 1D-normally distributed samples.
146  * Matrix classes can be CMatrixTemplateNumeric or CMatrixFixedNumeric
147  * \sa drawGaussian1D
148  */
149  template <class MAT>
151  MAT &matrix,
152  const double mean = 0,
153  const double std = 1 )
154  {
155  for (size_t r=0;r<matrix.getRowCount();r++)
156  for (size_t c=0;c<matrix.getColCount();c++)
157  matrix.get_unsafe(r,c) = static_cast<typename MAT::Scalar>( drawGaussian1D(mean,std) );
158  }
159 
160  /** Generates a random definite-positive matrix of the given size, using the formula C = v*v^t + epsilon*I, with "v" being a vector of gaussian random samples.
161  */
162  mrpt::math::CMatrixDouble drawDefinitePositiveMatrix(const size_t dim, const double std_scale = 1.0, const double diagonal_epsilon = 1e-8);
163 
164  /** Fills the given vector with independent, 1D-normally distributed samples.
165  * \sa drawGaussian1D
166  */
167  template <class VEC>
169  VEC & v,
170  const double mean = 0,
171  const double std = 1 )
172  {
173  const size_t N = v.size();
174  for (size_t c=0;c<N;c++)
175  v[c] = static_cast<typename mrpt::math::ContainerType<VEC>::element_t>( drawGaussian1D(mean,std) );
176  }
177 
178  /** Generate multidimensional random samples according to a given covariance matrix.
179  * Mean is assumed to be zero if mean==NULL.
180  * \exception std::exception On invalid covariance matrix
181  * \sa drawGaussianMultivariateMany
182  */
183  template <typename T>
184  void drawGaussianMultivariate(
185  std::vector<T> &out_result,
187  const std::vector<T>* mean = NULL
188  );
189 
190 
191  /** Generate multidimensional random samples according to a given covariance matrix.
192  * Mean is assumed to be zero if mean==NULL.
193  * \exception std::exception On invalid covariance matrix
194  * \sa drawGaussianMultivariateMany
195  */
196  template <class VECTORLIKE,class COVMATRIX>
198  VECTORLIKE &out_result,
199  const COVMATRIX &cov,
200  const VECTORLIKE* mean = NULL
201  )
202  {
203  const size_t N = cov.rows();
204  ASSERT_(cov.rows()==cov.cols())
205  if (mean) ASSERT_EQUAL_(size_t(mean->size()),N)
206 
207  // Compute eigenvalues/eigenvectors of cov:
208  Eigen::SelfAdjointEigenSolver<typename COVMATRIX::PlainObject> eigensolver(cov);
209 
210  typename Eigen::SelfAdjointEigenSolver<typename COVMATRIX::PlainObject>::MatrixType eigVecs = eigensolver.eigenvectors();
211  typename Eigen::SelfAdjointEigenSolver<typename COVMATRIX::PlainObject>::RealVectorType eigVals = eigensolver.eigenvalues();
212 
213  // Scale eigenvectors with eigenvalues:
214  // D.Sqrt(); Z = Z * D; (for each column)
215  eigVals = eigVals.array().sqrt();
216  for (typename COVMATRIX::Index i=0;i<eigVecs.cols();i++)
217  eigVecs.col(i) *= eigVals[i];
218 
219  // Set size of output vector:
220  out_result.assign(N,0);
221 
222  for (size_t i=0;i<N;i++)
223  {
224  typename COVMATRIX::Scalar rnd = drawGaussian1D_normalized();
225  for (size_t d=0;d<N;d++)
226  out_result[d]+= eigVecs.coeff(d,i) * rnd;
227  }
228  if (mean)
229  for (size_t d=0;d<N;d++)
230  out_result[d]+= (*mean)[d];
231  }
232 
233  /** Generate a given number of multidimensional random samples according to a given covariance matrix.
234  * \param cov The covariance matrix where to draw the samples from.
235  * \param desiredSamples The number of samples to generate.
236  * \param ret The output list of samples
237  * \param mean The mean, or zeros if mean==NULL.
238  */
239  template <typename VECTOR_OF_VECTORS,typename COVMATRIX>
241  VECTOR_OF_VECTORS &ret,
242  size_t desiredSamples,
243  const COVMATRIX &cov,
244  const typename VECTOR_OF_VECTORS::value_type *mean = NULL )
245  {
246  ASSERT_EQUAL_(cov.cols(),cov.rows())
247  if (mean) ASSERT_EQUAL_(size_t(mean->size()),size_t(cov.cols()))
248 
249  // Compute eigenvalues/eigenvectors of cov:
250  Eigen::SelfAdjointEigenSolver<typename COVMATRIX::PlainObject> eigensolver(cov);
251 
252  typename Eigen::SelfAdjointEigenSolver<typename COVMATRIX::PlainObject>::MatrixType eigVecs = eigensolver.eigenvectors();
253  typename Eigen::SelfAdjointEigenSolver<typename COVMATRIX::PlainObject>::RealVectorType eigVals = eigensolver.eigenvalues();
254 
255  // Scale eigenvectors with eigenvalues:
256  // D.Sqrt(); Z = Z * D; (for each column)
257  eigVals = eigVals.array().sqrt();
258  for (typename COVMATRIX::Index i=0;i<eigVecs.cols();i++)
259  eigVecs.col(i) *= eigVals[i];
260 
261  // Set size of output vector:
262  ret.resize(desiredSamples);
263  const size_t N = cov.cols();
264  for (size_t k=0;k<desiredSamples;k++)
265  {
266  ret[k].assign(N,0);
267  for (size_t i=0;i<N;i++)
268  {
269  typename COVMATRIX::Scalar rnd = drawGaussian1D_normalized();
270  for (size_t d=0;d<N;d++)
271  ret[k][d]+= eigVecs.coeff(d,i) * rnd;
272  }
273  if (mean)
274  for (size_t d=0;d<N;d++)
275  ret[k][d]+= (*mean)[d];
276  }
277  }
278 
279 
280  /** @} */
281 
282 
283  /** @name Miscellaneous
284  @{ */
285 
286  /** Returns a random permutation of a vector: all the elements of the input vector are in the output but at random positions.
287  */
288  template <class VEC>
289  void permuteVector(const VEC &in_vector, VEC &out_result)
290  {
291  out_result = in_vector;
292  const size_t N = out_result.size();
293  if (N>1)
294  std::random_shuffle( &out_result[0],&out_result[N-1] );
295  }
296 
297  /** @} */
298 
299  }; // end of CRandomGenerator --------------------------------------------------------------
300 
301 
302  /** A static instance of a CRandomGenerator class, for use in single-thread applications */
303  extern BASE_IMPEXP CRandomGenerator randomGenerator;
304 
305 
306  /** A random number generator for usage in STL algorithms expecting a function like this (eg, random_shuffle):
307  */
308  inline ptrdiff_t random_generator_for_STL(ptrdiff_t i)
309  {
310  return randomGenerator.drawUniform32bit() % i;
311  }
312 
313  /** Fills the given matrix with independent, uniformly distributed samples.
314  * Matrix classes can be CMatrixTemplateNumeric or CMatrixFixedNumeric
315  * \sa matrixRandomNormal
316  */
317  template <class MAT>
319  MAT &matrix,
320  const double unif_min = 0,
321  const double unif_max = 1 )
322  {
323  for (size_t r=0;r<matrix.getRowCount();r++)
324  for (size_t c=0;c<matrix.getColCount();c++)
325  matrix.get_unsafe(r,c) = static_cast<typename MAT::Scalar>( randomGenerator.drawUniform(unif_min,unif_max) );
326  }
327 
328  /** Fills the given matrix with independent, uniformly distributed samples.
329  * \sa vectorRandomNormal
330  */
331  template <class T>
333  std::vector<T> &v_out,
334  const T& unif_min = 0,
335  const T& unif_max = 1 )
336  {
337  size_t n = v_out.size();
338  for (size_t r=0;r<n;r++)
339  v_out[r] = randomGenerator.drawUniform(unif_min,unif_max);
340  }
341 
342  /** Fills the given matrix with independent, normally distributed samples.
343  * Matrix classes can be CMatrixTemplateNumeric or CMatrixFixedNumeric
344  * \sa matrixRandomUni
345  */
346  template <class MAT>
348  MAT &matrix,
349  const double mean = 0,
350  const double std = 1 )
351  {
352  for (size_t r=0;r<matrix.getRowCount();r++)
353  for (size_t c=0;c<matrix.getColCount();c++)
354  matrix.get_unsafe(r,c) = static_cast<typename MAT::Scalar>( mean + std*randomGenerator.drawGaussian1D_normalized() );
355  }
356 
357  /** Generates a random vector with independent, normally distributed samples.
358  * \sa matrixRandomUni
359  */
360  template <class T>
362  std::vector<T> &v_out,
363  const T& mean = 0,
364  const T& std = 1 )
365  {
366  size_t n = v_out.size();
367  for (size_t r=0;r<n;r++)
368  v_out[r] = mean + std*randomGenerator.drawGaussian1D_normalized();
369  }
370 
371  /** Randomize the generators.
372  * A seed can be providen, or a current-time based seed can be used (default)
373  */
374  inline void Randomize(const uint32_t seed) {
375  randomGenerator.randomize(seed);
376  }
377  inline void Randomize() {
378  randomGenerator.randomize();
379  }
380 
381  /** Returns a random permutation of a vector: all the elements of the input vector are in the output but at random positions.
382  */
383  template <class T>
385  const std::vector<T> &in_vector,
386  std::vector<T> &out_result)
387  {
388  randomGenerator.permuteVector(in_vector,out_result);
389  }
390 
391 
392  /** Generate multidimensional random samples according to a given covariance matrix.
393  * \exception std::exception On invalid covariance matrix
394  * \sa randomNormalMultiDimensionalMany
395  */
396  template <typename T>
399  std::vector<T> &out_result)
400  {
401  randomGenerator.drawGaussianMultivariate(out_result,cov);
402  }
403 
404  /** Generate a given number of multidimensional random samples according to a given covariance matrix.
405  * \param cov The covariance matrix where to draw the samples from.
406  * \param desiredSamples The number of samples to generate.
407  * \param samplesLikelihoods If desired, set to a valid pointer to a vector, where it will be stored the likelihoods of having obtained each sample: the product of the gaussian-pdf for each independent variable.
408  * \param ret The output list of samples
409  *
410  * \exception std::exception On invalid covariance matrix
411  *
412  * \sa randomNormalMultiDimensional
413  */
414  template <typename T>
417  size_t desiredSamples,
418  std::vector< std::vector<T> > &ret,
419  std::vector<T> *samplesLikelihoods = NULL)
420  {
421  randomGenerator.drawGaussianMultivariateMany(ret,desiredSamples,cov,static_cast<const std::vector<T>*>(NULL),samplesLikelihoods);
422  }
423 
424  /** Generate multidimensional random samples according to a given covariance matrix.
425  * \exception std::exception On invalid covariance matrix
426  * \sa randomNormalMultiDimensional
427  */
428  template <typename T,typename MATRIXLIKE>
430  const MATRIXLIKE &cov,
431  size_t desiredSamples,
432  std::vector< std::vector<T> > &ret )
433  {
434  randomGenerator.drawGaussianMultivariateMany(ret,desiredSamples,cov);
435  }
436 
437  /** Generate multidimensional random samples according to a given covariance matrix.
438  * \exception std::exception On invalid covariance matrix
439  * \sa randomNormalMultiDimensionalMany
440  */
441  template <typename T,typename MATRIXLIKE>
443  const MATRIXLIKE &cov,
444  std::vector<T> &out_result)
445  {
446  randomGenerator.drawGaussianMultivariate(out_result,cov);
447  }
448 
449 
450  }// End of namespace
451 
452 } // End of namespace
453 
454 #endif
#define ASSERT_EQUAL_(__A, __B)
double drawUniform(const double Min, const double Max)
Generate a uniformly distributed pseudo-random number using the MT19937 algorithm, scaled to the selected range.
void randomNormalMultiDimensional(const CMatrixTemplateNumeric< T > &cov, std::vector< T > &out_result)
Generate multidimensional random samples according to a given covariance matrix.
void permuteVector(const VEC &in_vector, VEC &out_result)
Returns a random permutation of a vector: all the elements of the input vector are in the output but ...
void drawUniformUnsignedInt(uint32_t &ret_number)
You can call this overloaded method with either 32 or 64bit unsigned ints for the sake of general cod...
BASE_IMPEXP CRandomGenerator randomGenerator
A static instance of a CRandomGenerator class, for use in single-thread applications.
void drawGaussian1DVector(VEC &v, const double mean=0, const double std=1)
Fills the given vector with independent, 1D-normally distributed samples.
void randomNormalMultiDimensionalMany(const CMatrixTemplateNumeric< T > &cov, size_t desiredSamples, std::vector< std::vector< T > > &ret, std::vector< T > *samplesLikelihoods=NULL)
Generate a given number of multidimensional random samples according to a given covariance matrix...
void drawGaussian1DMatrix(MAT &matrix, const double mean=0, const double std=1)
Fills the given matrix with independent, 1D-normally distributed samples.
STL namespace.
A thred-safe pseudo random number generator, based on an internal MT19937 randomness generator...
void randomPermutation(const std::vector< T > &in_vector, std::vector< T > &out_result)
Returns a random permutation of a vector: all the elements of the input vector are in the output but ...
void drawUniformUnsignedIntRange(T &ret_number, const U min_val, const V max_val)
Return a uniform unsigned integer in the range [min_val,max_val] (both inclusive) ...
void drawGaussianMultivariateMany(VECTOR_OF_VECTORS &ret, size_t desiredSamples, const COVMATRIX &cov, const typename VECTOR_OF_VECTORS::value_type *mean=NULL)
Generate a given number of multidimensional random samples according to a given covariance matrix...
double drawGaussian1D(const double mean, const double std)
Generate a normally distributed pseudo-random number.
void drawUniformMatrix(MAT &matrix, const double unif_min=0, const double unif_max=1)
Fills the given matrix with independent, uniformly distributed samples.
This base provides a set of functions for maths stuff.
Definition: CArray.h:18
Eigen::Matrix< typename MATRIX::Scalar, MATRIX::ColsAtCompileTime, MATRIX::ColsAtCompileTime > cov(const MATRIX &v)
Computes the covariance matrix from a list of samples in an NxM matrix, where each row is a sample...
Definition: ops_matrices.h:135
Data used internally by the MT19937 PRNG algorithm.
void drawGaussianMultivariate(VECTORLIKE &out_result, const COVMATRIX &cov, const VECTORLIKE *mean=NULL)
Generate multidimensional random samples according to a given covariance matrix.
void Randomize(const uint32_t seed)
Randomize the generators.
void vectorRandomNormal(std::vector< T > &v_out, const T &mean=0, const T &std=1)
Generates a random vector with independent, normally distributed samples.
CRandomGenerator(const uint32_t seed)
Constructor for providing a custom random seed to initialize the PRNG.
ptrdiff_t random_generator_for_STL(ptrdiff_t i)
A random number generator for usage in STL algorithms expecting a function like this (eg...
This is the global namespace for all Mobile Robot Programming Toolkit (MRPT) libraries.
void matrixRandomNormal(MAT &matrix, const double mean=0, const double std=1)
Fills the given matrix with independent, normally distributed samples.
void matrixRandomUni(MAT &matrix, const double unif_min=0, const double unif_max=1)
Fills the given matrix with independent, uniformly distributed samples.
#define ASSERT_(f)
CRandomGenerator()
Default constructor: initialize random seed based on current time.
CONTAINER::value_type element_t
Definition: math_frwds.h:85
void vectorRandomUni(std::vector< T > &v_out, const T &unif_min=0, const T &unif_max=1)
Fills the given matrix with independent, uniformly distributed samples.
EIGEN_STRONG_INLINE double mean() const
Computes the mean of the entire matrix.
void drawUniformVector(VEC &v, const double unif_min=0, const double unif_max=1)
Fills the given vector with independent, uniformly distributed samples.
void drawUniformUnsignedInt(uint64_t &ret_number)



Page generated by Doxygen 1.8.8 for MRPT 1.2.2 SVN:Unversioned directory at Tue Oct 14 02:14:08 UTC 2014