00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
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
00038
00039 namespace random
00040 {
00041 using namespace mrpt::utils;
00042 using namespace mrpt::math;
00043
00044
00045
00046
00047
00048
00049
00050
00051 class BASE_IMPEXP CRandomGenerator
00052 {
00053 protected:
00054
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
00070
00071
00072
00073 CRandomGenerator() : m_MT19937_data() { randomize(); }
00074
00075
00076 CRandomGenerator(const uint32_t seed) : m_MT19937_data() { randomize(seed); }
00077
00078 void randomize(const uint32_t seed);
00079 void randomize();
00080
00081
00082
00083
00084
00085
00086
00087
00088 uint32_t drawUniform32bit();
00089
00090
00091 double drawUniform( const double Min, const double Max) {
00092 return Min + (Max-Min)* drawUniform32bit() * 2.3283064370807973754314699618685e-10;
00093 }
00094
00095
00096
00097
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
00111
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
00127
00128
00129
00130
00131
00132 double drawGaussian1D_normalized( double *likelihood = NULL);
00133
00134
00135
00136
00137
00138 double drawGaussian1D( const double mean, const double std ) {
00139 return mean+std*drawGaussian1D_normalized();
00140 }
00141
00142
00143
00144
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
00158
00159 CMatrixDouble drawDefinitePositiveMatrix(const size_t dim, const double std_scale = 1.0, const double diagonal_epsilon = 1e-8);
00160
00161
00162
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
00176
00177
00178
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
00189
00190
00191
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
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
00211
00212 eigVals = eigVals.array().sqrt();
00213 for (typename COVMATRIX::Index i=0;i<eigVecs.cols();i++)
00214 eigVecs.col(i) *= eigVals[i];
00215
00216
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
00231
00232
00233
00234
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
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
00253
00254 eigVals = eigVals.array().sqrt();
00255 for (typename COVMATRIX::Index i=0;i<eigVecs.cols();i++)
00256 eigVecs.col(i) *= eigVals[i];
00257
00258
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
00281
00282
00283
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 };
00297
00298
00299
00300 extern BASE_IMPEXP CRandomGenerator randomGenerator;
00301
00302
00303
00304
00305 inline ptrdiff_t random_generator_for_STL(ptrdiff_t i)
00306 {
00307 return randomGenerator.drawUniform32bit() % i;
00308 }
00309
00310
00311
00312
00313 MRPT_DECLARE_DEPRECATED_FUNCTION("** deprecated **: Use mrpt::random::randomGenerator instead",
00314 double normalizedGaussian( double *likelihood = NULL) );
00315
00316
00317
00318
00319
00320 MRPT_DECLARE_DEPRECATED_FUNCTION("** deprecated **: Use mrpt::random::randomGenerator instead",
00321 double RandomNormal( double mean = 0, double std = 1) );
00322
00323
00324
00325
00326
00327 MRPT_DECLARE_DEPRECATED_FUNCTION("** deprecated **: Use mrpt::random::randomGenerator instead",
00328 uint32_t RandomUniInt() );
00329
00330
00331
00332
00333
00334
00335
00336 MRPT_DECLARE_DEPRECATED_FUNCTION("** deprecated **: Use mrpt::random::randomGenerator instead",
00337 double RandomUni( const double min, const double max) );
00338
00339
00340
00341
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
00355
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
00369
00370
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
00384
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
00398
00399
00400 inline void Randomize(const uint32_t seed) {
00401 randomGenerator.randomize(seed);
00402 }
00403 inline void Randomize() {
00404 randomGenerator.randomize();
00405 }
00406
00407
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
00419
00420
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
00431
00432
00433
00434
00435
00436
00437
00438
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
00451
00452
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
00464
00465
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 }
00477
00478 }
00479
00480 #endif