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
00089 MCPdf(unsigned int num_samples = 0, unsigned int dimension=0);
00091 virtual ~MCPdf();
00093 MCPdf(const MCPdf<T> &);
00094
00096 virtual MCPdf<T>* Clone() const;
00097
00098
00099 bool SampleFrom (Sample<T>& one_sample, int method = DEFAULT, void * args = NULL) const;
00100 bool SampleFrom (vector<Sample<T> >& list_samples, const unsigned int num_samples, int method = DEFAULT,
00101 void * args = NULL) const;
00102 T ExpectedValueGet() const;
00103 MatrixWrapper::SymmetricMatrix CovarianceGet() const;
00104
00105
00107
00110 void NumSamplesSet(unsigned int num_samples);
00111
00112
00114
00116 unsigned int NumSamplesGet() const;
00117
00119
00122 const WeightedSample<T>& SampleGet(unsigned int i) const;
00123
00125
00128 bool ListOfSamplesSet(const vector<WeightedSample<T> > & list_of_samples);
00130
00133 bool ListOfSamplesSet(const vector<Sample<T> > & list_of_samples);
00134
00136
00138 const vector<WeightedSample<T> > & ListOfSamplesGet() const;
00139
00141
00145 bool ListOfSamplesUpdate(const vector<WeightedSample<T> > & list_of_samples);
00146
00148
00152 bool ListOfSamplesUpdate(const vector<Sample<T> > & list_of_samples);
00153
00155
00158
00159
00161
00163 vector<double> & CumulativePDFGet();
00164
00165 };
00166
00168
00170
00171
00172 template <typename T> MCPdf<T>::MCPdf(unsigned int num_samples, unsigned int dimension) :
00173 Pdf<T>(dimension)
00174 , _CumSum(dimension)
00175 , _mean(dimension)
00176 , _diff(dimension)
00177 , _covariance(dimension)
00178 , _diffsum(dimension,dimension)
00179 {
00180 _SumWeights = 0;
00181 WeightedSample<T> my_sample(dimension);
00182 _listOfSamples.insert(_listOfSamples.begin(),num_samples,my_sample);
00183 _CumPDF.insert(_CumPDF.begin(),num_samples+1,0.0);
00184
00185 _los.assign(num_samples,WeightedSample<T>(dimension));
00186 _it_los = _los.begin();
00187 #ifdef __CONSTRUCTOR__
00188
00189 cout << "MCPDF Constructor: NumSamples = " << _listOfSamples.size()
00190 << ", CumPDF Samples = " << _CumPDF.size()
00191 << ", _SumWeights = " << _SumWeights << endl;
00192 #endif // __CONSTRUCTOR__
00193 }
00194
00195
00196 template <typename T>
00197 MCPdf<T>::~MCPdf()
00198 {
00199 #ifdef __DESTRUCTOR__
00200 cout << "MCPDF::Destructor" << endl;
00201 #endif // __DESTRUCTOR__
00202 }
00203
00204
00205 template <typename T>
00206 MCPdf<T>::MCPdf(const MCPdf & pdf) : Pdf<T>(pdf)
00207 , _CumSum(pdf.DimensionGet())
00208 , _mean(pdf.DimensionGet())
00209 , _diff(pdf.DimensionGet())
00210 , _covariance(pdf.DimensionGet())
00211 , _diffsum(pdf.DimensionGet(),pdf.DimensionGet())
00212 {
00213 this->_listOfSamples = pdf._listOfSamples;
00214 this->_CumPDF = pdf._CumPDF;
00215 _SumWeights = pdf._SumWeights;
00216 this->_los = pdf._listOfSamples;
00217 _it_los = _los.begin();
00218 #ifdef __CONSTRUCTOR__
00219 cout << "MCPDF Copy Constructor: NumSamples = " << _listOfSamples.size()
00220 << ", CumPDF Samples = " << _CumPDF.size()
00221 << ", SumWeights = " << _SumWeights << endl;
00222 #endif // __CONSTRUCTOR__
00223 }
00224
00225
00226 template <typename T> MCPdf<T>*
00227 MCPdf<T>::Clone() const
00228 {
00229 return new MCPdf<T>(*this);
00230 }
00231
00232 template <typename T> bool
00233 MCPdf<T>::SampleFrom (vector<Sample<T> >& list_samples,
00234 const unsigned int numsamples,
00235 int method,
00236 void * args) const
00237 {
00238 list_samples.resize(numsamples);
00239 switch(method)
00240 {
00241 case DEFAULT:
00242 {
00243 return Pdf<T>::SampleFrom(list_samples, numsamples,method,args);
00244 }
00245 case RIPLEY:
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256 {
00257 std::vector<double> unif_samples(numsamples);
00258 for ( unsigned int i = 0; i < numsamples ; i++)
00259 unif_samples[i] = runif();
00260
00261
00262 unif_samples[numsamples-1] = pow(unif_samples[numsamples-1], double (1.0/numsamples));
00263
00264
00265 if (numsamples > 1)
00266 for ( int i = numsamples-2; i >= 0 ; i--)
00267 unif_samples[i] = pow(unif_samples[i],double (1.0/(i+1))) * unif_samples[i+1];
00268
00269
00270 unsigned int index = 0;
00271 unsigned int size;
00272 size = _listOfSamples.size();
00273 typename vector<WeightedSample<T> >::const_iterator it = _listOfSamples.begin();
00274 typename vector<double>::const_iterator CumPDFit = _CumPDF.begin();
00275 typename vector<Sample<T> >::iterator sit = list_samples.begin();
00276
00277 for ( unsigned int i = 0; i < numsamples ; i++)
00278 {
00279 while ( unif_samples[i] > *CumPDFit )
00280 {
00281 assert(index <= size);
00282 index++; it++; CumPDFit++;
00283 }
00284 it--;
00285 *sit = *it;
00286 it++;
00287 sit++;
00288 }
00289 return true;
00290 }
00291 default:
00292 {
00293 cerr << "MCPdf::Samplefrom(int, void *): No such sampling method" << endl;
00294 return false;
00295 }
00296 }
00297 }
00298
00299 template <typename T> bool
00300 MCPdf<T>::SampleFrom(Sample<T>& one_sample, int method, void * args) const
00301 {
00302 switch(method)
00303 {
00304 case DEFAULT:
00305 {
00306
00307 double unif_sample; unif_sample = runif();
00308
00309 unsigned int index = 0;
00310 unsigned int size; size = _listOfSamples.size();
00311 typename vector<WeightedSample<T> >::const_iterator it;
00312 it = _listOfSamples.begin();
00313 typename vector<double>::const_iterator CumPDFit;
00314 CumPDFit = _CumPDF.begin();
00315
00316 while ( unif_sample > *CumPDFit )
00317 {
00318
00319 assert(index <= size);
00320 index++; it++; CumPDFit++;
00321 }
00322 it--;
00323 one_sample = *it;
00324 return true;
00325 }
00326 default:
00327 {
00328 cerr << "MCPdf::Samplefrom(int, void *): No such sampling method" << endl;
00329 return false;
00330 }
00331 }
00332 }
00333
00334
00335 template <typename T> unsigned int MCPdf<T>::NumSamplesGet() const
00336 {
00337 return _listOfSamples.size();
00338 }
00339
00340 template <typename T> const WeightedSample<T>&
00341 MCPdf<T>::SampleGet(unsigned int i) const
00342 {
00343 assert(i < NumSamplesGet());
00344 return _listOfSamples[i];
00345 }
00346
00347
00348 template <typename T> void
00349 MCPdf<T>::NumSamplesSet(unsigned int num_samples)
00350 {
00351 #ifdef __MCPDF_DEBUG__
00352 cout << "MCPDF::NumSamplesSet BEFORE: NumSamples " << _listOfSamples.size() << endl;
00353 cout << "MCPDF::NumSamplesSet BEFORE: CumPDF Samples " << _CumPDF.size() << endl;
00354 #endif // __MCPDF_DEBUG__
00355 unsigned int ns = num_samples;
00356 unsigned int size = _listOfSamples.size();
00357 static typename vector<double>::iterator CumPDFit;
00358 static typename vector<WeightedSample<T> >::iterator it;
00359 if (size < ns)
00360 {
00361 WeightedSample<T> ws;
00362 _listOfSamples.insert(_listOfSamples.end(),(ns - size),ws);
00363 _CumPDF.insert(_CumPDF.end(),(ns - size),0.0);
00364 }
00365 else if (size > ns)
00366 {
00367 it = _listOfSamples.begin(); CumPDFit = _CumPDF.begin();
00368 for ( unsigned int index = 0; index < (size-ns); index++ )
00369 {
00370 it = _listOfSamples.erase(it);
00371 CumPDFit = _CumPDF.erase(CumPDFit);
00372 }
00373 #ifdef __MCPDF_DEBUG__
00374 cout << "MCPDF::NumSamplesSet: WARNING DELETING SAMPLES!!" << endl;
00375 #endif // __MCPDF_DEBUG__
00376 }
00377 else {;}
00378 #ifdef __MCPDF_DEBUG__
00379 cout << "MCPDF::NumSamplesSet: Setting NumSamples to " << _listOfSamples.size() << endl;
00380 cout << "MCPDF::NumSamplesSet: Setting CumPDF Samples to " << _CumPDF.size() << endl;
00381 #endif // __MCPDF_DEBUG__
00382 }
00383
00384
00385
00386 template <typename T> bool
00387 MCPdf<T>::ListOfSamplesSet(const vector<WeightedSample<T> > & los)
00388 {
00389
00390 this->NumSamplesSet(los.size());
00391 _listOfSamples = los;
00392 #ifdef __MCPDF_DEBUG__
00393 cout << "MCPDF::ListOfSamplesSet: NumSamples = " << ListOfSamples.size() << endl;
00394 #endif // __MCPDF_DEBUG__
00395 return this->NormalizeWeights();
00396 }
00397
00398
00399 template <typename T> bool
00400 MCPdf<T>::ListOfSamplesSet(const vector<Sample<T> > & los)
00401 {
00402 unsigned int numsamples = los.size();
00403 typename vector<Sample<T> >::const_iterator lit; lit=los.begin();
00404 static typename vector<WeightedSample<T> >::iterator it;
00405
00406 this->NumSamplesSet(numsamples);
00407
00408 for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ )
00409 {
00410 *it = *lit; ;
00411 it->WeightSet(1.0/numsamples);
00412 lit++;
00413 }
00414 _SumWeights = 1.0;
00415
00416 this->CumPDFUpdate();
00417
00418 #ifdef __MCPDF_DEBUG__
00419 cout << "MCPDF ListOfSamplesSet: NumSamples = " << _listOfSamples.size()
00420 << " SumWeights = " << _SumWeights << endl;
00421 #endif // __MCPDF_DEBUG__
00422
00423 return true;
00424 }
00425
00426 template <typename T> const vector<WeightedSample<T> > &
00427 MCPdf<T>::ListOfSamplesGet() const
00428 {
00429 return _listOfSamples;
00430 }
00431
00432
00433 template <typename T> bool
00434 MCPdf<T>::ListOfSamplesUpdate(const vector<WeightedSample<T> > & los)
00435 {
00436 assert (los.size() == _listOfSamples.size());
00437 if (los.size() != 0)
00438 {
00439 _listOfSamples = los;
00440 return this->NormalizeWeights();
00441 }
00442 return true;
00443 }
00444
00445 template <typename T> bool
00446 MCPdf<T>::ListOfSamplesUpdate(const vector<Sample<T> > & los)
00447 {
00448 unsigned int numsamples = _listOfSamples.size();
00449 if ((numsamples = los.size()) == _listOfSamples.size())
00450 {
00451 assert (numsamples != 0);
00452 typename vector<Sample<T> >::const_iterator lit; lit=los.begin();
00453 static typename vector<WeightedSample<T> >::iterator it;
00454
00455 this->NumSamplesSet(numsamples);
00456
00457 for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ )
00458 {
00459 *it = *lit; ;
00460 it->WeightSet(1.0/numsamples);
00461 lit++;
00462 }
00463 _SumWeights = 1.0;
00464 this->CumPDFUpdate();
00465 }
00466 return true;
00467 }
00468
00469
00470 template <typename T> bool
00471 MCPdf<T>::SumWeightsUpdate()
00472 {
00473 double SumOfWeights = 0.0;
00474 double current_weight;
00475 static typename vector<WeightedSample<T> >::iterator it;
00476 for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ )
00477 {
00478 current_weight = it->WeightGet();
00479 SumOfWeights += current_weight;
00480 }
00481
00482 #ifdef __MCPDF_DEBUG__
00483 cout << "MCPDF::SumWeightsUpdate: SumWeights = " << SumOfWeights << endl;
00484 #endif // __MCPDF_DEBUG__
00485
00486 if (SumOfWeights > 0){
00487 this->_SumWeights = SumOfWeights;
00488 return true;
00489 }
00490 else{
00491 cerr << "MCPDF::SumWeightsUpdate: SumWeights = " << SumOfWeights << endl;
00492 return false;
00493 }
00494 }
00495
00496 template <typename T> bool
00497 MCPdf<T>::NormalizeWeights()
00498 {
00499 static typename vector<WeightedSample<T> >::iterator it;
00500
00501
00502 if (!this->SumWeightsUpdate()) return false;
00503
00504 for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ )
00505 {
00506 it->WeightSet(it->WeightGet() / _SumWeights);
00507 }
00508 this->_SumWeights = 1.0;
00509 this->CumPDFUpdate();
00510 return true;
00511 }
00512
00513
00514 template <typename T> void
00515 MCPdf<T>::CumPDFUpdate()
00516 {
00517 double CumSum=0.0;
00518 static typename vector<double>::iterator CumPDFit;
00519 static typename vector<WeightedSample<T> >::iterator it;
00520 CumPDFit = _CumPDF.begin(); *CumPDFit = 0.0;
00521
00522
00523 for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ )
00524 {
00525 CumPDFit++;
00526
00527 CumSum += ( it->WeightGet() / _SumWeights);
00528 *CumPDFit = CumSum;
00529 }
00530 }
00531
00532
00533 template <typename T>
00534 T MCPdf<T>::ExpectedValueGet ( ) const
00535 {
00536 cerr << "MCPDF ExpectedValueGet: not implemented for the template parameters you use."
00537 << endl << "Use template specialization as shown in mcpdf.cpp " << endl;
00538
00539 assert(0);
00540 T result;
00541 return result;
00542 }
00543
00544
00545 template <typename T>
00546 MatrixWrapper::SymmetricMatrix MCPdf<T>::CovarianceGet ( ) const
00547 {
00548 cerr << "MCPDF CovarianceGet: not implemented for the template parameters you use."
00549 << endl << "Use template specialization as shown in mcpdf.cpp " << endl;
00550
00551 assert(0);
00552 MatrixWrapper::SymmetricMatrix result;
00553 return result;
00554 }
00555
00556
00557
00558 template <typename T>
00559 vector<double> & MCPdf<T>::CumulativePDFGet()
00560 {
00561 return _CumPDF;
00562 }
00563
00564
00565
00566 }
00567
00568 #include "mcpdf.cpp"
00569
00570 #endif