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__
227 template <
typename T>
bool
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
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;
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>