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


Page generated by Doxygen 1.8.3 for MRPT 0.9.6 SVN: at Fri Feb 15 22:05:02 EST 2013