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 _MRPT_MONTE_CARLO_H_
00029 #define _MRPT_MONTE_CARLO_H_
00030
00031 #include <map>
00032 #include <cmath>
00033 #include <vector>
00034 #include <numeric>
00035 #include <algorithm>
00036 #include <stdexcept>
00037 #include <functional>
00038
00039 #include <mrpt/random.h>
00040 #include <mrpt/utils/CTicTac.h>
00041
00042 namespace mrpt { namespace math {
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057 template<typename T,typename NUM,typename OTHER> class CMonteCarlo {
00058 private:
00059 mrpt::random::CRandomGenerator gen;
00060 class CStatisticalAnalyzer {
00061 private:
00062 mrpt_base_vector<NUM> data;
00063 public:
00064 template<typename VEC> inline CStatisticalAnalyzer(const VEC &v1):data(v1.begin(),v1.end()) {}
00065 template<typename VEC> inline void setData(const VEC &v1) {
00066 data.assign(v1.begin(),v1.end());
00067 }
00068 template<typename VEC> inline void getData(VEC &v1) const {
00069 v1.assign(data.begin(),data.end());
00070 }
00071 template<typename VEC1,typename VEC2> inline void getDistribution(VEC1 &vx,VEC2 &vy,const NUM width=1.0) const {
00072 mrpt::vector_double vvx,vvy;
00073 getDistribution(vvx,vvy,width);
00074 vx.assign(vvx.begin(),vvx.end());
00075 vy.assign(vvy.begin(),vvy.end());
00076 }
00077
00078 inline void getDistribution(mrpt::vector_double &vx,mrpt::vector_double &vy,const NUM width=1.0) const {
00079 CHistogram hist(CHistogram::createWithFixedWidth(0,*max_element(data.begin(),data.end()),width));
00080 hist.add(data);
00081 hist.getHistogram(vx,vy);
00082 }
00083
00084 };
00085 public:
00086
00087
00088 T (*valueGenerator)(mrpt::random::CRandomGenerator &);
00089
00090 NUM (*errorFun1)(const T &);
00091
00092 OTHER (*intermediateFun)(const T &);
00093 NUM (*errorFun2)(const T &,const OTHER &);
00094 inline CMonteCarlo():gen(),valueGenerator(NULL),errorFun1(NULL),intermediateFun(NULL),errorFun2(NULL) {}
00095 NUM doExperiment(size_t N,double &time,bool showInWindow=false) {
00096 if (!valueGenerator) throw std::logic_error("Value generator function is not set.");
00097 std::vector<T> baseData(N);
00098 std::vector<NUM> errorData(N);
00099 mrpt::utils::CTicTac meter;
00100 for (size_t i=0;i<N;++i) baseData[i]=valueGenerator(gen);
00101 if (errorFun1) {
00102 meter.Tic();
00103 std::transform(baseData.begin(),baseData.end(),errorData.begin(),errorFun1);
00104 time=meter.Tac();
00105 } else {
00106 if (!intermediateFun||!errorFun2) throw std::logic_error("Experiment-related functions are not set.");
00107 std::vector<OTHER> intermediate(N);
00108 transform(baseData.begin(),baseData.end(),intermediate.begin(),intermediateFun);
00109 meter.Tic();
00110 for (size_t i=0;i<N;++i) errorData[i]=errorFun2(baseData[i],intermediate[i]);
00111 time=meter.Tac();
00112 }
00113 NUM res=accumulate(errorData.begin(),errorData.end(),NUM(0))/errorData.size();
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125 return res;
00126 }
00127 };
00128
00129 }}
00130 #endif