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