34 #include "../sample/weightedsample.h" 35 #include "../wrappers/rng/rng.h" 36 #include "../bfl_err.h" 49 template <
typename T>
class MCPdf:
public Pdf<T>
76 mutable vector<WeightedSample<T> > _los;
79 mutable SymmetricMatrix _covariance;
80 mutable Matrix _diffsum;
81 mutable typename vector<WeightedSample<T> >::iterator _it_los;
88 MCPdf(
unsigned int num_samples = 0,
unsigned int dimension=0);
98 bool SampleFrom (
Sample<T>& one_sample,
const SampleMthd method = SampleMthd::DEFAULT,
void * args = NULL)
const;
99 bool SampleFrom (vector<
Sample<T> >& list_samples,
const unsigned int num_samples,
const SampleMthd method = SampleMthd::DEFAULT,
100 void * args = NULL)
const;
171 template <
typename T>
MCPdf<T>::MCPdf(
unsigned int num_samples,
unsigned int dimension) :
173 , _covariance(dimension)
174 , _diffsum(dimension,dimension)
182 _it_los = _los.begin();
183 #ifdef __CONSTRUCTOR__ 185 cout <<
"MCPDF Constructor: NumSamples = " <<
_listOfSamples.size()
186 <<
", CumPDF Samples = " <<
_CumPDF.size()
188 #endif // __CONSTRUCTOR__ 194 template <
typename T>
197 #ifdef __DESTRUCTOR__ 198 cout <<
"MCPDF::Destructor" << endl;
199 #endif // __DESTRUCTOR__ 203 template <
typename T>
205 , _covariance(pdf.DimensionGet())
206 , _diffsum(pdf.DimensionGet(),pdf.DimensionGet())
211 this->_los = pdf._listOfSamples;
212 _it_los = _los.begin();
213 #ifdef __CONSTRUCTOR__ 214 cout <<
"MCPDF Copy Constructor: NumSamples = " <<
_listOfSamples.size()
215 <<
", CumPDF Samples = " <<
_CumPDF.size()
217 #endif // __CONSTRUCTOR__ 221 template <
typename T> MCPdf<T>*
224 return new MCPdf<T>(*this);
227 template <
typename T>
bool 228 MCPdf<T>::SampleFrom (vector<Sample<T> >& list_samples,
229 const unsigned int numsamples,
230 const SampleMthd method,
233 list_samples.resize(numsamples);
236 case SampleMthd::DEFAULT:
240 case SampleMthd::RIPLEY:
252 std::vector<double> unif_samples(numsamples);
253 for (
unsigned int i = 0; i < numsamples ; i++)
254 unif_samples[i] = runif();
257 unif_samples[numsamples-1] = pow(unif_samples[numsamples-1],
double (1.0/numsamples));
261 for (
int i = numsamples-2; i >= 0 ; i--)
262 unif_samples[i] = pow(unif_samples[i],
double (1.0/(i+1))) * unif_samples[i+1];
265 unsigned int index = 0;
267 size = _listOfSamples.size();
268 typename vector<WeightedSample<T> >::const_iterator it = _listOfSamples.begin();
269 typename vector<double>::const_iterator CumPDFit = _CumPDF.begin();
270 typename vector<Sample<T> >::iterator sit = list_samples.begin();
272 for (
unsigned int i = 0; i < numsamples ; i++)
274 while ( unif_samples[i] > *CumPDFit )
276 assert(index <= size);
277 index++; it++; CumPDFit++;
288 cerr <<
"MCPdf::Samplefrom(int, void *): No such sampling method" << endl;
294 template <
typename T>
bool 295 MCPdf<T>::SampleFrom(Sample<T>& one_sample,
const SampleMthd method,
void * args)
const 299 case SampleMthd::DEFAULT:
302 double unif_sample; unif_sample = runif();
304 unsigned int index = 0;
305 unsigned int size; size = _listOfSamples.size();
306 typename vector<WeightedSample<T> >::const_iterator it;
307 it = _listOfSamples.begin();
308 typename vector<double>::const_iterator CumPDFit;
309 CumPDFit = _CumPDF.begin();
311 while ( unif_sample > *CumPDFit )
314 assert(index <= size);
315 index++; it++; CumPDFit++;
323 cerr <<
"MCPdf::Samplefrom(int, void *): No such sampling method" << endl;
332 return _listOfSamples.size();
335 template <
typename T>
const WeightedSample<T>&
338 assert(i < NumSamplesGet());
339 return _listOfSamples[i];
343 template <
typename T>
void 346 #ifdef __MCPDF_DEBUG__ 347 cout <<
"MCPDF::NumSamplesSet BEFORE: NumSamples " << _listOfSamples.size() << endl;
348 cout <<
"MCPDF::NumSamplesSet BEFORE: CumPDF Samples " << _CumPDF.size() << endl;
349 #endif // __MCPDF_DEBUG__ 350 unsigned int ns = num_samples;
351 unsigned int size = _listOfSamples.size();
352 static typename vector<double>::iterator CumPDFit;
353 static typename vector<WeightedSample<T> >::iterator it;
356 WeightedSample<T> ws;
357 _listOfSamples.insert(_listOfSamples.end(),(ns - size),ws);
358 _CumPDF.insert(_CumPDF.end(),(ns - size),0.0);
362 it = _listOfSamples.begin(); CumPDFit = _CumPDF.begin();
363 for (
unsigned int index = 0; index < (size-ns); index++ )
365 it = _listOfSamples.erase(it);
366 CumPDFit = _CumPDF.erase(CumPDFit);
368 #ifdef __MCPDF_DEBUG__ 369 cout <<
"MCPDF::NumSamplesSet: WARNING DELETING SAMPLES!!" << endl;
370 #endif // __MCPDF_DEBUG__ 373 #ifdef __MCPDF_DEBUG__ 374 cout <<
"MCPDF::NumSamplesSet: Setting NumSamples to " << _listOfSamples.size() << endl;
375 cout <<
"MCPDF::NumSamplesSet: Setting CumPDF Samples to " << _CumPDF.size() << endl;
376 #endif // __MCPDF_DEBUG__ 381 template <
typename T>
bool 385 this->NumSamplesSet(los.size());
386 _listOfSamples = los;
387 #ifdef __MCPDF_DEBUG__ 388 cout <<
"MCPDF::ListOfSamplesSet: NumSamples = " << ListOfSamples.size() << endl;
389 #endif // __MCPDF_DEBUG__ 390 return this->NormalizeWeights();
394 template <
typename T>
bool 397 unsigned int numsamples = los.size();
398 typename vector<Sample<T> >::const_iterator lit; lit=los.begin();
399 static typename vector<WeightedSample<T> >::iterator it;
401 this->NumSamplesSet(numsamples);
403 for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ )
406 it->WeightSet(1.0/numsamples);
411 this->CumPDFUpdate();
413 #ifdef __MCPDF_DEBUG__ 414 cout <<
"MCPDF ListOfSamplesSet: NumSamples = " << _listOfSamples.size()
415 <<
" SumWeights = " << _SumWeights << endl;
416 #endif // __MCPDF_DEBUG__ 421 template <
typename T>
const vector<WeightedSample<T> > &
424 return _listOfSamples;
428 template <
typename T>
bool 431 assert (los.size() == _listOfSamples.size());
434 _listOfSamples = los;
435 return this->NormalizeWeights();
440 template <
typename T>
bool 443 unsigned int numsamples = _listOfSamples.size();
444 if ((numsamples = los.size()) == _listOfSamples.size())
446 assert (numsamples != 0);
447 typename vector<Sample<T> >::const_iterator lit; lit=los.begin();
448 static typename vector<WeightedSample<T> >::iterator it;
450 this->NumSamplesSet(numsamples);
452 for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ )
455 it->WeightSet(1.0/numsamples);
459 this->CumPDFUpdate();
465 template <
typename T>
bool 468 double SumOfWeights = 0.0;
469 double current_weight;
470 static typename vector<WeightedSample<T> >::iterator it;
471 for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ )
473 current_weight = it->WeightGet();
474 SumOfWeights += current_weight;
477 #ifdef __MCPDF_DEBUG__ 478 cout <<
"MCPDF::SumWeightsUpdate: SumWeights = " << SumOfWeights << endl;
479 #endif // __MCPDF_DEBUG__ 481 if (SumOfWeights > 0){
482 this->_SumWeights = SumOfWeights;
486 cerr <<
"MCPDF::SumWeightsUpdate: SumWeights = " << SumOfWeights << endl;
491 template <
typename T>
bool 494 static typename vector<WeightedSample<T> >::iterator it;
497 if (!this->SumWeightsUpdate())
return false;
499 for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ )
501 it->WeightSet(it->WeightGet() / _SumWeights);
503 this->_SumWeights = 1.0;
504 this->CumPDFUpdate();
509 template <
typename T>
void 513 static typename vector<double>::iterator CumPDFit;
514 static typename vector<WeightedSample<T> >::iterator it;
515 CumPDFit = _CumPDF.begin(); *CumPDFit = 0.0;
518 for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ )
522 CumSum += ( it->WeightGet() / _SumWeights);
528 template <
typename T>
531 cerr <<
"MCPDF ExpectedValueGet: not implemented for the template parameters you use." 532 << endl <<
"Use template specialization as shown in mcpdf.cpp " << endl;
540 template <
typename T>
543 cerr <<
"MCPDF CovarianceGet: not implemented for the template parameters you use." 544 << endl <<
"Use template specialization as shown in mcpdf.cpp " << endl;
553 template <
typename T>
vector< double > & CumulativePDFGet()
Add a sample to the list.
Class PDF: Virtual Base class representing Probability Density Functions.
virtual bool SampleFrom(vector< Sample< T > > &list_samples, const unsigned int num_samples, const SampleMthd method=SampleMthd::DEFAULT, void *args=NULL) const
Draw multiple samples from the Pdf (overloaded)
vector< double > _CumPDF
STL-iterator.
unsigned int NumSamplesGet() const
Get number of samples.
bool ListOfSamplesUpdate(const vector< WeightedSample< T > > &list_of_samples)
Update the list of samples (overloaded)
virtual MCPdf< T > * Clone() const
Clone function.
MCPdf(unsigned int num_samples=0, unsigned int dimension=0)
Constructor.
void NumSamplesSet(unsigned int num_samples)
Set number of samples.
const WeightedSample< T > & SampleGet(unsigned int i) const
Get one sample.
Monte Carlo Pdf: Sample based implementation of Pdf.
MatrixWrapper::SymmetricMatrix CovarianceGet() const
Get the Covariance Matrix E[(x - E[x])^2] of the Analytic pdf.
void CumPDFUpdate()
After updating weights, we have to update the cumPDF.
const vector< WeightedSample< T > > & ListOfSamplesGet() const
Get the list of weighted samples.
double _SumWeights
Sum of all weights: used for normalising purposes.
T ExpectedValueGet() const
Get the expected value E[x] of the pdf.
bool SumWeightsUpdate()
STL-iterator for cumulative PDF list.
bool NormalizeWeights()
Normalizing the weights.
virtual ~MCPdf()
destructor
bool ListOfSamplesSet(const vector< WeightedSample< T > > &list_of_samples)
Set the list of weighted samples.
vector< WeightedSample< T > > _listOfSamples
STL-list containing the list of samples.