Go to the documentation of this file.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
00029
00030 #ifndef MCPDF_H
00031 #define MCPDF_H
00032
00033 #include "pdf.h"
00034 #include "../sample/weightedsample.h"
00035 #include "../wrappers/rng/rng.h"
00036 #include "../bfl_err.h"
00037 #include <list>
00038 #include <vector>
00039 #include <cassert>
00040
00041 namespace BFL
00042 {
00043
00045
00049 template <typename T> class MCPdf: public Pdf<T>
00050 {
00051 protected:
00052
00054 double _SumWeights;
00056 vector<WeightedSample<T> > _listOfSamples;
00058
00060 vector<double> _CumPDF;
00062
00063
00064
00066 bool SumWeightsUpdate();
00068 bool NormalizeWeights();
00070 void CumPDFUpdate();
00071
00072 private:
00073
00074
00075 mutable T _CumSum;
00076 mutable vector<WeightedSample<T> > _los;
00077 mutable T _mean;
00078 mutable T _diff;
00079 mutable SymmetricMatrix _covariance;
00080 mutable Matrix _diffsum;
00081 mutable typename vector<WeightedSample<T> >::iterator _it_los;
00082
00083 public:
00085
00088 MCPdf(unsigned int num_samples = 0, unsigned int dimension=0);
00090 virtual ~MCPdf();
00092 MCPdf(const MCPdf<T> &);
00093
00095 virtual MCPdf<T>* Clone() const;
00096
00097
00098 bool SampleFrom (Sample<T>& one_sample, int method = DEFAULT, void * args = NULL) const;
00099 bool SampleFrom (vector<Sample<T> >& list_samples, const unsigned int num_samples, int method = DEFAULT,
00100 void * args = NULL) const;
00101 T ExpectedValueGet() const;
00102 MatrixWrapper::SymmetricMatrix CovarianceGet() const;
00103
00104
00106
00109 void NumSamplesSet(unsigned int num_samples);
00110
00111
00113
00115 unsigned int NumSamplesGet() const;
00116
00118
00121 const WeightedSample<T>& SampleGet(unsigned int i) const;
00122
00124
00127 bool ListOfSamplesSet(const vector<WeightedSample<T> > & list_of_samples);
00129
00132 bool ListOfSamplesSet(const vector<Sample<T> > & list_of_samples);
00133
00135
00137 const vector<WeightedSample<T> > & ListOfSamplesGet() const;
00138
00140
00144 bool ListOfSamplesUpdate(const vector<WeightedSample<T> > & list_of_samples);
00145
00147
00151 bool ListOfSamplesUpdate(const vector<Sample<T> > & list_of_samples);
00152
00154
00157
00158
00160
00162 vector<double> & CumulativePDFGet();
00163
00164 };
00165
00167
00169
00170
00171 template <typename T> MCPdf<T>::MCPdf(unsigned int num_samples, unsigned int dimension) :
00172 Pdf<T>(dimension)
00173 , _covariance(dimension)
00174 , _diffsum(dimension,dimension)
00175 {
00176 _SumWeights = 0;
00177 WeightedSample<T> my_sample(dimension);
00178 _listOfSamples.insert(_listOfSamples.begin(),num_samples,my_sample);
00179 _CumPDF.insert(_CumPDF.begin(),num_samples+1,0.0);
00180
00181 _los.assign(num_samples,WeightedSample<T>(dimension));
00182 _it_los = _los.begin();
00183 #ifdef __CONSTRUCTOR__
00184
00185 cout << "MCPDF Constructor: NumSamples = " << _listOfSamples.size()
00186 << ", CumPDF Samples = " << _CumPDF.size()
00187 << ", _SumWeights = " << _SumWeights << endl;
00188 #endif // __CONSTRUCTOR__
00189 }
00190
00191
00192
00193
00194 template <typename T>
00195 MCPdf<T>::~MCPdf()
00196 {
00197 #ifdef __DESTRUCTOR__
00198 cout << "MCPDF::Destructor" << endl;
00199 #endif // __DESTRUCTOR__
00200 }
00201
00202
00203 template <typename T>
00204 MCPdf<T>::MCPdf(const MCPdf & pdf) : Pdf<T>(pdf)
00205 , _covariance(pdf.DimensionGet())
00206 , _diffsum(pdf.DimensionGet(),pdf.DimensionGet())
00207 {
00208 this->_listOfSamples = pdf._listOfSamples;
00209 this->_CumPDF = pdf._CumPDF;
00210 _SumWeights = pdf._SumWeights;
00211 this->_los = pdf._listOfSamples;
00212 _it_los = _los.begin();
00213 #ifdef __CONSTRUCTOR__
00214 cout << "MCPDF Copy Constructor: NumSamples = " << _listOfSamples.size()
00215 << ", CumPDF Samples = " << _CumPDF.size()
00216 << ", SumWeights = " << _SumWeights << endl;
00217 #endif // __CONSTRUCTOR__
00218 }
00219
00220
00221 template <typename T> MCPdf<T>*
00222 MCPdf<T>::Clone() const
00223 {
00224 return new MCPdf<T>(*this);
00225 }
00226
00227 template <typename T> bool
00228 MCPdf<T>::SampleFrom (vector<Sample<T> >& list_samples,
00229 const unsigned int numsamples,
00230 int method,
00231 void * args) const
00232 {
00233 list_samples.resize(numsamples);
00234 switch(method)
00235 {
00236 case DEFAULT:
00237 {
00238 return Pdf<T>::SampleFrom(list_samples, numsamples,method,args);
00239 }
00240 case RIPLEY:
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251 {
00252 std::vector<double> unif_samples(numsamples);
00253 for ( unsigned int i = 0; i < numsamples ; i++)
00254 unif_samples[i] = runif();
00255
00256
00257 unif_samples[numsamples-1] = pow(unif_samples[numsamples-1], double (1.0/numsamples));
00258
00259
00260 if (numsamples > 1)
00261 for ( int i = numsamples-2; i >= 0 ; i--)
00262 unif_samples[i] = pow(unif_samples[i],double (1.0/(i+1))) * unif_samples[i+1];
00263
00264
00265 unsigned int index = 0;
00266 unsigned int size;
00267 size = _listOfSamples.size();
00268 typename vector<WeightedSample<T> >::const_iterator it = _listOfSamples.begin();
00269 typename vector<double>::const_iterator CumPDFit = _CumPDF.begin();
00270 typename vector<Sample<T> >::iterator sit = list_samples.begin();
00271
00272 for ( unsigned int i = 0; i < numsamples ; i++)
00273 {
00274 while ( unif_samples[i] > *CumPDFit )
00275 {
00276 assert(index <= size);
00277 index++; it++; CumPDFit++;
00278 }
00279 it--;
00280 *sit = *it;
00281 it++;
00282 sit++;
00283 }
00284 return true;
00285 }
00286 default:
00287 {
00288 cerr << "MCPdf::Samplefrom(int, void *): No such sampling method" << endl;
00289 return false;
00290 }
00291 }
00292 }
00293
00294 template <typename T> bool
00295 MCPdf<T>::SampleFrom(Sample<T>& one_sample, int method, void * args) const
00296 {
00297 switch(method)
00298 {
00299 case DEFAULT:
00300 {
00301
00302 double unif_sample; unif_sample = runif();
00303
00304 unsigned int index = 0;
00305 unsigned int size; size = _listOfSamples.size();
00306 typename vector<WeightedSample<T> >::const_iterator it;
00307 it = _listOfSamples.begin();
00308 typename vector<double>::const_iterator CumPDFit;
00309 CumPDFit = _CumPDF.begin();
00310
00311 while ( unif_sample > *CumPDFit )
00312 {
00313
00314 assert(index <= size);
00315 index++; it++; CumPDFit++;
00316 }
00317 it--;
00318 one_sample = *it;
00319 return true;
00320 }
00321 default:
00322 {
00323 cerr << "MCPdf::Samplefrom(int, void *): No such sampling method" << endl;
00324 return false;
00325 }
00326 }
00327 }
00328
00329
00330 template <typename T> unsigned int MCPdf<T>::NumSamplesGet() const
00331 {
00332 return _listOfSamples.size();
00333 }
00334
00335 template <typename T> const WeightedSample<T>&
00336 MCPdf<T>::SampleGet(unsigned int i) const
00337 {
00338 assert(i < NumSamplesGet());
00339 return _listOfSamples[i];
00340 }
00341
00342
00343 template <typename T> void
00344 MCPdf<T>::NumSamplesSet(unsigned int num_samples)
00345 {
00346 #ifdef __MCPDF_DEBUG__
00347 cout << "MCPDF::NumSamplesSet BEFORE: NumSamples " << _listOfSamples.size() << endl;
00348 cout << "MCPDF::NumSamplesSet BEFORE: CumPDF Samples " << _CumPDF.size() << endl;
00349 #endif // __MCPDF_DEBUG__
00350 unsigned int ns = num_samples;
00351 unsigned int size = _listOfSamples.size();
00352 static typename vector<double>::iterator CumPDFit;
00353 static typename vector<WeightedSample<T> >::iterator it;
00354 if (size < ns)
00355 {
00356 WeightedSample<T> ws;
00357 _listOfSamples.insert(_listOfSamples.end(),(ns - size),ws);
00358 _CumPDF.insert(_CumPDF.end(),(ns - size),0.0);
00359 }
00360 else if (size > ns)
00361 {
00362 it = _listOfSamples.begin(); CumPDFit = _CumPDF.begin();
00363 for ( unsigned int index = 0; index < (size-ns); index++ )
00364 {
00365 it = _listOfSamples.erase(it);
00366 CumPDFit = _CumPDF.erase(CumPDFit);
00367 }
00368 #ifdef __MCPDF_DEBUG__
00369 cout << "MCPDF::NumSamplesSet: WARNING DELETING SAMPLES!!" << endl;
00370 #endif // __MCPDF_DEBUG__
00371 }
00372 else {;}
00373 #ifdef __MCPDF_DEBUG__
00374 cout << "MCPDF::NumSamplesSet: Setting NumSamples to " << _listOfSamples.size() << endl;
00375 cout << "MCPDF::NumSamplesSet: Setting CumPDF Samples to " << _CumPDF.size() << endl;
00376 #endif // __MCPDF_DEBUG__
00377 }
00378
00379
00380
00381 template <typename T> bool
00382 MCPdf<T>::ListOfSamplesSet(const vector<WeightedSample<T> > & los)
00383 {
00384
00385 this->NumSamplesSet(los.size());
00386 _listOfSamples = los;
00387 #ifdef __MCPDF_DEBUG__
00388 cout << "MCPDF::ListOfSamplesSet: NumSamples = " << ListOfSamples.size() << endl;
00389 #endif // __MCPDF_DEBUG__
00390 return this->NormalizeWeights();
00391 }
00392
00393
00394 template <typename T> bool
00395 MCPdf<T>::ListOfSamplesSet(const vector<Sample<T> > & los)
00396 {
00397 unsigned int numsamples = los.size();
00398 typename vector<Sample<T> >::const_iterator lit; lit=los.begin();
00399 static typename vector<WeightedSample<T> >::iterator it;
00400
00401 this->NumSamplesSet(numsamples);
00402
00403 for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ )
00404 {
00405 *it = *lit; ;
00406 it->WeightSet(1.0/numsamples);
00407 lit++;
00408 }
00409 _SumWeights = 1.0;
00410
00411 this->CumPDFUpdate();
00412
00413 #ifdef __MCPDF_DEBUG__
00414 cout << "MCPDF ListOfSamplesSet: NumSamples = " << _listOfSamples.size()
00415 << " SumWeights = " << _SumWeights << endl;
00416 #endif // __MCPDF_DEBUG__
00417
00418 return true;
00419 }
00420
00421 template <typename T> const vector<WeightedSample<T> > &
00422 MCPdf<T>::ListOfSamplesGet() const
00423 {
00424 return _listOfSamples;
00425 }
00426
00427
00428 template <typename T> bool
00429 MCPdf<T>::ListOfSamplesUpdate(const vector<WeightedSample<T> > & los)
00430 {
00431 assert (los.size() == _listOfSamples.size());
00432 if (los.size() != 0)
00433 {
00434 _listOfSamples = los;
00435 return this->NormalizeWeights();
00436 }
00437 return true;
00438 }
00439
00440 template <typename T> bool
00441 MCPdf<T>::ListOfSamplesUpdate(const vector<Sample<T> > & los)
00442 {
00443 unsigned int numsamples = _listOfSamples.size();
00444 if ((numsamples = los.size()) == _listOfSamples.size())
00445 {
00446 assert (numsamples != 0);
00447 typename vector<Sample<T> >::const_iterator lit; lit=los.begin();
00448 static typename vector<WeightedSample<T> >::iterator it;
00449
00450 this->NumSamplesSet(numsamples);
00451
00452 for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ )
00453 {
00454 *it = *lit; ;
00455 it->WeightSet(1.0/numsamples);
00456 lit++;
00457 }
00458 _SumWeights = 1.0;
00459 this->CumPDFUpdate();
00460 }
00461 return true;
00462 }
00463
00464
00465 template <typename T> bool
00466 MCPdf<T>::SumWeightsUpdate()
00467 {
00468 double SumOfWeights = 0.0;
00469 double current_weight;
00470 static typename vector<WeightedSample<T> >::iterator it;
00471 for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ )
00472 {
00473 current_weight = it->WeightGet();
00474 SumOfWeights += current_weight;
00475 }
00476
00477 #ifdef __MCPDF_DEBUG__
00478 cout << "MCPDF::SumWeightsUpdate: SumWeights = " << SumOfWeights << endl;
00479 #endif // __MCPDF_DEBUG__
00480
00481 if (SumOfWeights > 0){
00482 this->_SumWeights = SumOfWeights;
00483 return true;
00484 }
00485 else{
00486 cerr << "MCPDF::SumWeightsUpdate: SumWeights = " << SumOfWeights << endl;
00487 return false;
00488 }
00489 }
00490
00491 template <typename T> bool
00492 MCPdf<T>::NormalizeWeights()
00493 {
00494 static typename vector<WeightedSample<T> >::iterator it;
00495
00496
00497 if (!this->SumWeightsUpdate()) return false;
00498
00499 for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ )
00500 {
00501 it->WeightSet(it->WeightGet() / _SumWeights);
00502 }
00503 this->_SumWeights = 1.0;
00504 this->CumPDFUpdate();
00505 return true;
00506 }
00507
00508
00509 template <typename T> void
00510 MCPdf<T>::CumPDFUpdate()
00511 {
00512 double CumSum=0.0;
00513 static typename vector<double>::iterator CumPDFit;
00514 static typename vector<WeightedSample<T> >::iterator it;
00515 CumPDFit = _CumPDF.begin(); *CumPDFit = 0.0;
00516
00517
00518 for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ )
00519 {
00520 CumPDFit++;
00521
00522 CumSum += ( it->WeightGet() / _SumWeights);
00523 *CumPDFit = CumSum;
00524 }
00525 }
00526
00527
00528 template <typename T>
00529 T MCPdf<T>::ExpectedValueGet ( ) const
00530 {
00531 cerr << "MCPDF ExpectedValueGet: not implemented for the template parameters you use."
00532 << endl << "Use template specialization as shown in mcpdf.cpp " << endl;
00533
00534 assert(0);
00535 T result;
00536 return result;
00537 }
00538
00539
00540 template <typename T>
00541 MatrixWrapper::SymmetricMatrix MCPdf<T>::CovarianceGet ( ) const
00542 {
00543 cerr << "MCPDF CovarianceGet: not implemented for the template parameters you use."
00544 << endl << "Use template specialization as shown in mcpdf.cpp " << endl;
00545
00546 assert(0);
00547 MatrixWrapper::SymmetricMatrix result;
00548 return result;
00549 }
00550
00551
00552
00553 template <typename T>
00554 vector<double> & MCPdf<T>::CumulativePDFGet()
00555 {
00556 return _CumPDF;
00557 }
00558
00559
00560
00561 }
00562
00563 #include "mcpdf.cpp"
00564
00565 #endif
bfl
Author(s): Klaas Gadeyne, Wim Meeussen, Tinne Delaet and many others. See web page for a full contributor list. ROS package maintained by Wim Meeussen.
autogenerated on Sun Oct 5 2014 22:29:53