Main MRPT website > C++ reference
MRPT logo

RandomGenerators.h

Go to the documentation of this file.
00001 /* +---------------------------------------------------------------------------+
00002    |          The Mobile Robot Programming Toolkit (MRPT) C++ library          |
00003    |                                                                           |
00004    |                   http://mrpt.sourceforge.net/                            |
00005    |                                                                           |
00006    |   Copyright (C) 2005-2011  University of Malaga                           |
00007    |                                                                           |
00008    |    This software was written by the Machine Perception and Intelligent    |
00009    |      Robotics Lab, University of Malaga (Spain).                          |
00010    |    Contact: Jose-Luis Blanco  <jlblanco@ctima.uma.es>                     |
00011    |                                                                           |
00012    |  This file is part of the MRPT project.                                   |
00013    |                                                                           |
00014    |     MRPT is free software: you can redistribute it and/or modify          |
00015    |     it under the terms of the GNU General Public License as published by  |
00016    |     the Free Software Foundation, either version 3 of the License, or     |
00017    |     (at your option) any later version.                                   |
00018    |                                                                           |
00019    |   MRPT is distributed in the hope that it will be useful,                 |
00020    |     but WITHOUT ANY WARRANTY; without even the implied warranty of        |
00021    |     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         |
00022    |     GNU General Public License for more details.                          |
00023    |                                                                           |
00024    |     You should have received a copy of the GNU General Public License     |
00025    |     along with MRPT.  If not, see <http://www.gnu.org/licenses/>.         |
00026    |                                                                           |
00027    +---------------------------------------------------------------------------+ */
00028 #ifndef RandomGenerator_H
00029 #define RandomGenerator_H
00030 
00031 #include <mrpt/utils/utils_defs.h>
00032 #include <mrpt/math/CMatrixTemplateNumeric.h>
00033 #include <mrpt/math/ops_matrices.h>
00034 
00035 namespace mrpt
00036 {
00037         /** A namespace of pseudo-random numbers genrators of diferent distributions. The central class in this namespace is mrpt::random::CRandomGenerator
00038          */
00039         namespace random
00040         {
00041                 using namespace mrpt::utils;
00042                 using namespace mrpt::math;
00043 
00044                 /** A thred-safe pseudo random number generator, based on an internal MT19937 randomness generator.
00045                   * The base algorithm for randomness is platform-independent. See http://en.wikipedia.org/wiki/Mersenne_twister
00046                   *
00047                   * For real thread-safety, each thread must create and use its own instance of this class.
00048                   *
00049                   * Single-thread programs can use the static object mrpt::random::randomGenerator
00050                   */
00051                 class BASE_IMPEXP CRandomGenerator
00052                 {
00053                 protected:
00054                         /** Data used internally by the MT19937 PRNG algorithm. */
00055                         struct  TMT19937_data
00056                         {
00057                                 TMT19937_data() : index(0), seed_initialized(false)
00058                                 {}
00059                                 uint32_t        MT[624];
00060                                 uint32_t        index;
00061                                 bool            seed_initialized;
00062                         } m_MT19937_data;
00063 
00064                         void MT19937_generateNumbers();
00065                         void MT19937_initializeGenerator(const uint32_t &seed);
00066 
00067                 public:
00068 
00069                         /** @name Initialization
00070                          @{ */
00071 
00072                                 /** Default constructor: initialize random seed based on current time */
00073                                 CRandomGenerator() : m_MT19937_data() { randomize(); }
00074 
00075                                 /** Constructor for providing a custom random seed to initialize the PRNG */
00076                                 CRandomGenerator(const uint32_t seed) : m_MT19937_data() { randomize(seed); }
00077 
00078                                 void randomize(const uint32_t seed);  //!< Initialize the PRNG from the given random seed
00079                                 void randomize();       //!< Randomize the generators, based on current time
00080 
00081                         /** @} */
00082 
00083                         /** @name Uniform pdf
00084                          @{ */
00085 
00086                                 /** Generate a uniformly distributed pseudo-random number using the MT19937 algorithm, in the whole range of 32-bit integers.
00087                                   *  See: http://en.wikipedia.org/wiki/Mersenne_twister */
00088                                 uint32_t drawUniform32bit();
00089 
00090                                 /** Generate a uniformly distributed pseudo-random number using the MT19937 algorithm, scaled to the selected range. */
00091                                 double drawUniform( const double Min, const double Max) {
00092                                         return Min + (Max-Min)* drawUniform32bit() * 2.3283064370807973754314699618685e-10; // 0xFFFFFFFF ^ -1
00093                                 }
00094 
00095                                 /** Fills the given matrix with independent, uniformly distributed samples.
00096                                   * Matrix classes can be CMatrixTemplateNumeric or CMatrixFixedNumeric
00097                                   * \sa drawUniform
00098                                   */
00099                                 template <class MAT>
00100                                 void drawUniformMatrix(
00101                                         MAT &matrix,
00102                                         const  double unif_min = 0,
00103                                         const  double unif_max = 1 )
00104                                 {
00105                                         for (size_t r=0;r<matrix.getRowCount();r++)
00106                                                 for (size_t c=0;c<matrix.getColCount();c++)
00107                                                         matrix.get_unsafe(r,c) = static_cast<typename MAT::value_type>( drawUniform(unif_min,unif_max) );
00108                                 }
00109 
00110                                 /** Fills the given vector with independent, uniformly distributed samples.
00111                                   * \sa drawUniform
00112                                   */
00113                                 template <class VEC>
00114                                 void drawUniformVector(
00115                                         VEC & v,
00116                                         const  double unif_min = 0,
00117                                         const  double unif_max = 1 )
00118                                 {
00119                                         const size_t N = v.size();
00120                                         for (size_t c=0;c<N;c++)
00121                                                 v[c] = static_cast<typename VEC::value_type>( drawUniform(unif_min,unif_max) );
00122                                 }
00123 
00124                         /** @} */
00125 
00126                         /** @name Normal/Gaussian pdf
00127                          @{ */
00128 
00129                                 /** Generate a normalized (mean=0, std=1) normally distributed sample.
00130                                  *  \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.
00131                                  */
00132                                 double drawGaussian1D_normalized( double *likelihood = NULL);
00133 
00134                                 /** Generate a normally distributed pseudo-random number.
00135                                  * \param mean The mean value of desired normal distribution
00136                                  * \param std  The standard deviation value of desired normal distribution
00137                                  */
00138                                 double drawGaussian1D( const double mean, const double std ) {
00139                                         return mean+std*drawGaussian1D_normalized();
00140                                 }
00141 
00142                                 /** Fills the given matrix with independent, 1D-normally distributed samples.
00143                                   * Matrix classes can be CMatrixTemplateNumeric or CMatrixFixedNumeric
00144                                   * \sa drawGaussian1D
00145                                   */
00146                                 template <class MAT>
00147                                 void drawGaussian1DMatrix(
00148                                         MAT &matrix,
00149                                         const double mean = 0,
00150                                         const double std = 1 )
00151                                 {
00152                                         for (size_t r=0;r<matrix.getRowCount();r++)
00153                                                 for (size_t c=0;c<matrix.getColCount();c++)
00154                                                         matrix.get_unsafe(r,c) = static_cast<typename MAT::value_type>( drawGaussian1D(mean,std) );
00155                                 }
00156 
00157                                 /** 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.
00158                                   */
00159                                 CMatrixDouble drawDefinitePositiveMatrix(const size_t dim, const double std_scale = 1.0, const double diagonal_epsilon = 1e-8);
00160 
00161                                 /** Fills the given vector with independent, 1D-normally distributed samples.
00162                                   * \sa drawGaussian1D
00163                                   */
00164                                 template <class VEC>
00165                                 void drawGaussian1DVector(
00166                                         VEC & v,
00167                                         const double mean = 0,
00168                                         const double std = 1 )
00169                                 {
00170                                         const size_t N = v.size();
00171                                         for (size_t c=0;c<N;c++)
00172                                                 v[c] = static_cast<typename VEC::value_type>( drawGaussian1D(mean,std) );
00173                                 }
00174 
00175                                 /** Generate multidimensional random samples according to a given covariance matrix.
00176                                  *  Mean is assumed to be zero if mean==NULL.
00177                                  * \exception std::exception On invalid covariance matrix
00178                                  * \sa drawGaussianMultivariateMany
00179                                  */
00180                                  template <typename T>
00181                                  void  drawGaussianMultivariate(
00182                                         std::vector<T>          &out_result,
00183                                         const CMatrixTemplateNumeric<T> &cov,
00184                                         const std::vector<T>*  mean = NULL
00185                                         );
00186 
00187 
00188                                 /** Generate multidimensional random samples according to a given covariance matrix.
00189                                  *  Mean is assumed to be zero if mean==NULL.
00190                                  * \exception std::exception On invalid covariance matrix
00191                                  * \sa drawGaussianMultivariateMany
00192                                  */
00193                                  template <class VECTORLIKE,class COVMATRIX>
00194                                  void  drawGaussianMultivariate(
00195                                         VECTORLIKE      &out_result,
00196                                         const COVMATRIX &cov,
00197                                         const VECTORLIKE* mean = NULL
00198                                         )
00199                                 {
00200                                         const size_t N = cov.rows();
00201                                         ASSERT_(cov.rows()==cov.cols())
00202                                         if (mean) ASSERT_EQUAL_(size_t(mean->size()),N)
00203 
00204                                         // Compute eigenvalues/eigenvectors of cov:
00205                                         Eigen::SelfAdjointEigenSolver<typename COVMATRIX::PlainObject> eigensolver(cov);
00206 
00207                                         typename Eigen::SelfAdjointEigenSolver<typename COVMATRIX::PlainObject>::MatrixType eigVecs = eigensolver.eigenvectors();
00208                                         typename Eigen::SelfAdjointEigenSolver<typename COVMATRIX::PlainObject>::RealVectorType eigVals = eigensolver.eigenvalues();
00209 
00210                                         // Scale eigenvectors with eigenvalues:
00211                                         // D.Sqrt(); Z = Z * D; (for each column)
00212                                         eigVals = eigVals.array().sqrt();
00213                                         for (typename COVMATRIX::Index i=0;i<eigVecs.cols();i++)
00214                                                 eigVecs.col(i) *= eigVals[i];
00215 
00216                                         // Set size of output vector:
00217                                         out_result.assign(N,0);
00218 
00219                                         for (size_t i=0;i<N;i++)
00220                                         {
00221                                                 typename COVMATRIX::Scalar rnd = drawGaussian1D_normalized();
00222                                                 for (size_t d=0;d<N;d++)
00223                                                         out_result[d]+= eigVecs.coeff(d,i) * rnd;
00224                                         }
00225                                         if (mean)
00226                                                 for (size_t d=0;d<N;d++)
00227                                                         out_result[d]+= (*mean)[d];
00228                                 }
00229 
00230                                 /** Generate a given number of multidimensional random samples according to a given covariance matrix.
00231                                  * \param cov The covariance matrix where to draw the samples from.
00232                                  * \param desiredSamples The number of samples to generate.
00233                                  * \param ret The output list of samples
00234                                  * \param mean The mean, or zeros if mean==NULL.
00235                                  */
00236                                  template <typename VECTOR_OF_VECTORS,typename COVMATRIX>
00237                                  void  drawGaussianMultivariateMany(
00238                                         VECTOR_OF_VECTORS       &ret,
00239                                         size_t               desiredSamples,
00240                                         const COVMATRIX     &cov,
00241                                         const typename VECTOR_OF_VECTORS::value_type *mean = NULL )
00242                                 {
00243                                         ASSERT_EQUAL_(cov.cols(),cov.rows())
00244                                         if (mean) ASSERT_EQUAL_(size_t(mean->size()),size_t(cov.cols()))
00245 
00246                                         // Compute eigenvalues/eigenvectors of cov:
00247                                         Eigen::SelfAdjointEigenSolver<typename COVMATRIX::PlainObject> eigensolver(cov);
00248 
00249                                         typename Eigen::SelfAdjointEigenSolver<typename COVMATRIX::PlainObject>::MatrixType eigVecs = eigensolver.eigenvectors();
00250                                         typename Eigen::SelfAdjointEigenSolver<typename COVMATRIX::PlainObject>::RealVectorType eigVals = eigensolver.eigenvalues();
00251 
00252                                         // Scale eigenvectors with eigenvalues:
00253                                         // D.Sqrt(); Z = Z * D; (for each column)
00254                                         eigVals = eigVals.array().sqrt();
00255                                         for (typename COVMATRIX::Index i=0;i<eigVecs.cols();i++)
00256                                                 eigVecs.col(i) *= eigVals[i];
00257 
00258                                         // Set size of output vector:
00259                                         ret.resize(desiredSamples);
00260                                         const size_t N = cov.cols();
00261                                         for (size_t k=0;k<desiredSamples;k++)
00262                                         {
00263                                                 ret[k].assign(N,0);
00264                                                 for (size_t i=0;i<N;i++)
00265                                                 {
00266                                                         typename COVMATRIX::Scalar rnd = drawGaussian1D_normalized();
00267                                                         for (size_t d=0;d<N;d++)
00268                                                                 ret[k][d]+= eigVecs.coeff(d,i) * rnd;
00269                                                 }
00270                                                 if (mean)
00271                                                         for (size_t d=0;d<N;d++)
00272                                                                 ret[k][d]+= (*mean)[d];
00273                                         }
00274                                 }
00275 
00276 
00277                         /** @} */
00278 
00279 
00280                         /** @name Miscellaneous
00281                          @{ */
00282 
00283                                 /** Returns a random permutation of a vector: all the elements of the input vector are in the output but at random positions.
00284                                   */
00285                                 template <class VEC>
00286                                 void  permuteVector(const VEC &in_vector, VEC &out_result)
00287                                 {
00288                                         out_result = in_vector;
00289                                         const size_t N = out_result.size();
00290                                         if (N>1)
00291                                                 std::random_shuffle( &out_result[0],&out_result[N-1] );
00292                                 }
00293 
00294                         /** @} */
00295 
00296                 }; // end of CRandomGenerator --------------------------------------------------------------
00297 
00298 
00299                 /** A static instance of a CRandomGenerator class, for use in single-thread applications */
00300                 extern BASE_IMPEXP CRandomGenerator randomGenerator;
00301 
00302 
00303                 /** A random number generator for usage in STL algorithms expecting a function like this (eg, random_shuffle):
00304                   */
00305                 inline ptrdiff_t random_generator_for_STL(ptrdiff_t i)
00306                 {
00307                         return randomGenerator.drawUniform32bit() % i;
00308                 }
00309 
00310                 /** Generate a normalized normally distributed pseudo-random number.
00311                  *  \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.
00312                  */
00313                 MRPT_DECLARE_DEPRECATED_FUNCTION("** deprecated **: Use mrpt::random::randomGenerator instead",
00314                 double normalizedGaussian( double *likelihood = NULL) );
00315 
00316                 /** Generate a normally distributed pseudo-random number.
00317                  * \param mean The mean value of desired normal distribution
00318                  * \param std  The standard deviation value of desired normal distribution
00319                  */
00320                 MRPT_DECLARE_DEPRECATED_FUNCTION("** deprecated **: Use mrpt::random::randomGenerator instead",
00321                 double RandomNormal( double mean = 0, double std = 1) );
00322 
00323                 /** Generate a uniformly distributed pseudo-random number using the MT19937 algorithm, in the whole range of 32-bit integers.
00324                   *  See: http://en.wikipedia.org/wiki/Mersenne_twister
00325                   * \sa RandomUni, Randomize
00326                   */
00327                 MRPT_DECLARE_DEPRECATED_FUNCTION("** deprecated **: Use mrpt::random::randomGenerator instead",
00328                 uint32_t RandomUniInt() );
00329 
00330                 /** Generate a uniformly distributed pseudo-random number using the MT19937 algorithm, scaled to the selected range.
00331                   *  This function uses internally RandomUniInt to generate the number, then scales it to the desired range.
00332                   *  Since MRPT 0.6.0 the MT19937 algorithm is used instead of C runtime library "rand" version.
00333                   *  See: http://en.wikipedia.org/wiki/Mersenne_twister
00334                   * \sa RandomUniInt, Randomize
00335                   */
00336                 MRPT_DECLARE_DEPRECATED_FUNCTION("** deprecated **: Use mrpt::random::randomGenerator instead",
00337                 double RandomUni( const double min, const double max) );
00338 
00339                 /** Fills the given matrix with independent, uniformly distributed samples.
00340                   * Matrix classes can be CMatrixTemplateNumeric or CMatrixFixedNumeric
00341                   * \sa matrixRandomNormal
00342                   */
00343                 template <class MAT>
00344                 void matrixRandomUni(
00345                         MAT &matrix,
00346                         const  double unif_min = 0,
00347                         const  double unif_max = 1 )
00348                 {
00349                         for (size_t r=0;r<matrix.getRowCount();r++)
00350                                 for (size_t c=0;c<matrix.getColCount();c++)
00351                                         matrix.get_unsafe(r,c) = static_cast<typename MAT::value_type>( randomGenerator.drawUniform(unif_min,unif_max) );
00352                 }
00353 
00354                 /** Fills the given matrix with independent, uniformly distributed samples.
00355                   * \sa vectorRandomNormal
00356                   */
00357                 template <class T>
00358                 void vectorRandomUni(
00359                         std::vector<T> &v_out,
00360                         const  T& unif_min = 0,
00361                         const  T& unif_max = 1 )
00362                 {
00363                         size_t n = v_out.size();
00364                         for (size_t r=0;r<n;r++)
00365                                 v_out[r] = randomGenerator.drawUniform(unif_min,unif_max);
00366                 }
00367 
00368                 /** Fills the given matrix with independent, normally distributed samples.
00369                   * Matrix classes can be CMatrixTemplateNumeric or CMatrixFixedNumeric
00370                   * \sa matrixRandomUni
00371                   */
00372                 template <class MAT>
00373                 void matrixRandomNormal(
00374                         MAT &matrix,
00375                         const double mean = 0,
00376                         const double std = 1 )
00377                 {
00378                         for (size_t r=0;r<matrix.getRowCount();r++)
00379                                 for (size_t c=0;c<matrix.getColCount();c++)
00380                                         matrix.get_unsafe(r,c) = static_cast<typename MAT::value_type>( mean + std*randomGenerator.drawGaussian1D_normalized() );
00381                 }
00382 
00383                 /** Generates a random vector with independent, normally distributed samples.
00384                   * \sa matrixRandomUni
00385                   */
00386                 template <class T>
00387                 void vectorRandomNormal(
00388                         std::vector<T> &v_out,
00389                         const  T& mean = 0,
00390                         const  T& std = 1 )
00391                 {
00392                         size_t n = v_out.size();
00393                         for (size_t r=0;r<n;r++)
00394                                 v_out[r] = mean + std*randomGenerator.drawGaussian1D_normalized();
00395                 }
00396 
00397                 /** Randomize the generators.
00398                  *   A seed can be providen, or a current-time based seed can be used (default)
00399                  */
00400                 inline void Randomize(const uint32_t seed)  {
00401                         randomGenerator.randomize(seed);
00402                 }
00403                 inline void Randomize()  {
00404                         randomGenerator.randomize();
00405                 }
00406 
00407                 /** Returns a random permutation of a vector: all the elements of the input vector are in the output but at random positions.
00408                   */
00409                 template <class T>
00410                 void  randomPermutation(
00411                         const std::vector<T> &in_vector,
00412                         std::vector<T>       &out_result)
00413                 {
00414                         randomGenerator.permuteVector(in_vector,out_result);
00415                 }
00416 
00417 
00418                 /** Generate multidimensional random samples according to a given covariance matrix.
00419                  * \exception std::exception On invalid covariance matrix
00420                  * \sa randomNormalMultiDimensionalMany
00421                  */
00422                 template <typename T>
00423                 void  randomNormalMultiDimensional(
00424                         const CMatrixTemplateNumeric<T> &cov,
00425                         std::vector<T>          &out_result)
00426                  {
00427                         randomGenerator.drawGaussianMultivariate(out_result,cov);
00428                  }
00429 
00430                  /** Generate a given number of multidimensional random samples according to a given covariance matrix.
00431                  * \param cov The covariance matrix where to draw the samples from.
00432                  * \param desiredSamples The number of samples to generate.
00433                  * \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.
00434                  * \param ret The output list of samples
00435                  *
00436                  * \exception std::exception On invalid covariance matrix
00437                  *
00438                  * \sa randomNormalMultiDimensional
00439                  */
00440                  template <typename T>
00441                  void  randomNormalMultiDimensionalMany(
00442                         const CMatrixTemplateNumeric<T> &cov,
00443                         size_t                                                  desiredSamples,
00444                         std::vector< std::vector<T> >   &ret,
00445                         std::vector<T>                                  *samplesLikelihoods = NULL)
00446                 {
00447                         randomGenerator.drawGaussianMultivariateMany(ret,desiredSamples,cov,static_cast<const std::vector<T>*>(NULL),samplesLikelihoods);
00448                 }
00449 
00450                 /** Generate multidimensional random samples according to a given covariance matrix.
00451                  * \exception std::exception On invalid covariance matrix
00452                  * \sa randomNormalMultiDimensional
00453                  */
00454                  template <typename T,size_t N>
00455                  void  randomNormalMultiDimensionalMany(
00456                         const CMatrixFixedNumeric<T,N,N> &cov,
00457                         size_t                                                  desiredSamples,
00458                         std::vector< std::vector<T> >   &ret )
00459                  {
00460                          randomGenerator.drawGaussianMultivariateMany(ret,desiredSamples,cov);
00461                  }
00462 
00463                 /** Generate multidimensional random samples according to a given covariance matrix.
00464                  * \exception std::exception On invalid covariance matrix
00465                  * \sa randomNormalMultiDimensionalMany
00466                  */
00467                  template <typename T,size_t N>
00468                  void  randomNormalMultiDimensional(
00469                         const CMatrixFixedNumeric<T,N,N> &cov,
00470                         std::vector<T>          &out_result)
00471                 {
00472                         randomGenerator.drawGaussianMultivariate(out_result,cov);
00473                 }
00474 
00475 
00476         }// End of namespace
00477 
00478 } // End of namespace
00479 
00480 #endif



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN: at Sat Mar 26 06:16:28 UTC 2011