mcpdf.h
Go to the documentation of this file.
1 // $Id$
2 // Copyright (C) 2002 Klaas Gadeyne <first dot last at gmail dot com>
3  /***************************************************************************
4  * This library is free software; you can redistribute it and/or *
5  * modify it under the terms of the GNU General Public *
6  * License as published by the Free Software Foundation; *
7  * version 2 of the License. *
8  * *
9  * As a special exception, you may use this file as part of a free *
10  * software library without restriction. Specifically, if other files *
11  * instantiate templates or use macros or inline functions from this *
12  * file, or you compile this file and link it with other files to *
13  * produce an executable, this file does not by itself cause the *
14  * resulting executable to be covered by the GNU General Public *
15  * License. This exception does not however invalidate any other *
16  * reasons why the executable file might be covered by the GNU General *
17  * Public License. *
18  * *
19  * This library is distributed in the hope that it will be useful, *
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
22  * Lesser General Public License for more details. *
23  * *
24  * You should have received a copy of the GNU General Public *
25  * License along with this library; if not, write to the Free Software *
26  * Foundation, Inc., 59 Temple Place, *
27  * Suite 330, Boston, MA 02111-1307 USA *
28  * *
29  ***************************************************************************/
30 #ifndef MCPDF_H
31 #define MCPDF_H
32 
33 #include "pdf.h"
34 #include "../sample/weightedsample.h"
35 #include "../wrappers/rng/rng.h"
36 #include "../bfl_err.h"
37 #include <list>
38 #include <vector>
39 #include <cassert>
40 
41 namespace BFL
42 {
43 
45 
49  template <typename T> class MCPdf: public Pdf<T>
50  {
51  protected:
52 
54  double _SumWeights;
56  vector<WeightedSample<T> > _listOfSamples;
58  // vector<WeightedSample<T> >::iterator _it;
60  vector<double> _CumPDF;
62  // typename vector<double>::iterator CumPDFit;
63 
64 
66  bool SumWeightsUpdate();
68  bool NormalizeWeights();
70  void CumPDFUpdate();
71 
72  private:
73  // Variables to avoid allocation during call of
74  //expectedvalueget/covarianceget
75  mutable T _CumSum;
76  mutable vector<WeightedSample<T> > _los;
77  mutable T _mean;
78  mutable T _diff;
79  mutable SymmetricMatrix _covariance;
80  mutable Matrix _diffsum;
81  mutable typename vector<WeightedSample<T> >::iterator _it_los;
82 
83  public:
85 
88  MCPdf(unsigned int num_samples = 0, unsigned int dimension=0);
90  virtual ~MCPdf();
92  MCPdf(const MCPdf<T> &);
93 
95  virtual MCPdf<T>* Clone() const;
96 
97  // implemented virtual functions
98  bool SampleFrom (Sample<T>& one_sample, int method = DEFAULT, void * args = NULL) const;
99  bool SampleFrom (vector<Sample<T> >& list_samples, const unsigned int num_samples, int method = DEFAULT,
100  void * args = NULL) const;
101  T ExpectedValueGet() const;
102  MatrixWrapper::SymmetricMatrix CovarianceGet() const;
103 
104 
106 
109  void NumSamplesSet(unsigned int num_samples);
110 
111 
113 
115  unsigned int NumSamplesGet() const;
116 
118 
121  const WeightedSample<T>& SampleGet(unsigned int i) const;
122 
124 
127  bool ListOfSamplesSet(const vector<WeightedSample<T> > & list_of_samples);
129 
132  bool ListOfSamplesSet(const vector<Sample<T> > & list_of_samples);
133 
135 
137  const vector<WeightedSample<T> > & ListOfSamplesGet() const;
138 
140 
144  bool ListOfSamplesUpdate(const vector<WeightedSample<T> > & list_of_samples);
145 
147 
151  bool ListOfSamplesUpdate(const vector<Sample<T> > & list_of_samples);
152 
154 
157  // void SampleAdd(WeightedSample<T> sample);
158 
160 
162  vector<double> & CumulativePDFGet();
163 
164  };
165 
167  // Template Code here
169 
170  // Constructor
171  template <typename T> MCPdf<T>::MCPdf(unsigned int num_samples, unsigned int dimension) :
172  Pdf<T>(dimension)
173  , _covariance(dimension)
174  , _diffsum(dimension,dimension)
175  {
176  _SumWeights = 0;
177  WeightedSample<T> my_sample(dimension);
178  _listOfSamples.insert(_listOfSamples.begin(),num_samples,my_sample);
179  _CumPDF.insert(_CumPDF.begin(),num_samples+1,0.0);
180 
181  _los.assign(num_samples,WeightedSample<T>(dimension));
182  _it_los = _los.begin();
183 #ifdef __CONSTRUCTOR__
184  // if (num_samples > 0)
185  cout << "MCPDF Constructor: NumSamples = " << _listOfSamples.size()
186  << ", CumPDF Samples = " << _CumPDF.size()
187  << ", _SumWeights = " << _SumWeights << endl;
188 #endif // __CONSTRUCTOR__
189  }
190 
191 
192 
193  // Destructor
194  template <typename T>
196  {
197 #ifdef __DESTRUCTOR__
198  cout << "MCPDF::Destructor" << endl;
199 #endif // __DESTRUCTOR__
200  }
201 
202  // Copy constructor
203  template <typename T>
204  MCPdf<T>::MCPdf(const MCPdf & pdf) : Pdf<T>(pdf)
205  , _covariance(pdf.DimensionGet())
206  , _diffsum(pdf.DimensionGet(),pdf.DimensionGet())
207  {
208  this->_listOfSamples = pdf._listOfSamples;
209  this->_CumPDF = pdf._CumPDF;
210  _SumWeights = pdf._SumWeights;
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()
216  << ", SumWeights = " << _SumWeights << endl;
217 #endif // __CONSTRUCTOR__
218  }
219 
220  //Clone function
221  template <typename T> MCPdf<T>*
222  MCPdf<T>::Clone() const
223  {
224  return new MCPdf<T>(*this);
225  }
226 
227  template <typename T> bool
228  MCPdf<T>::SampleFrom (vector<Sample<T> >& list_samples,
229  const unsigned int numsamples,
230  int method,
231  void * args) const
232  {
233  list_samples.resize(numsamples);
234  switch(method)
235  {
236  case DEFAULT: // O(N log(N) efficiency)
237  {
238  return Pdf<T>::SampleFrom(list_samples, numsamples,method,args);
239  }
240  case RIPLEY: // Only possible here ( O(N) efficiency )
241  /* See
242  @Book{ ripley87,
243  author = {Ripley, Brian D.},
244  title = {Stochastic Simulation},
245  publisher = {John Wiley and Sons},
246  year = {1987},
247  annote = {ISBN 0271-6356, WBIB 1 519.245}
248  }
249  */
250  // GENERATE N ORDERED IID UNIFORM SAMPLES
251  {
252  std::vector<double> unif_samples(numsamples);
253  for ( unsigned int i = 0; i < numsamples ; i++)
254  unif_samples[i] = runif();
255 
256  /* take n-th racine of u_N */
257  unif_samples[numsamples-1] = pow(unif_samples[numsamples-1], double (1.0/numsamples));
258  /* rescale other samples */
259  // only resample if more than one sample
260  if (numsamples > 1)
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];
263 
264  // CHECK WHERE THESE SAMPLES ARE IN _CUMPDF
265  unsigned int index = 0;
266  unsigned int size;
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();
271 
272  for ( unsigned int i = 0; i < numsamples ; i++)
273  {
274  while ( unif_samples[i] > *CumPDFit )
275  {
276  assert(index <= size);
277  index++; it++; CumPDFit++;
278  }
279  it--;
280  *sit = *it;
281  it++;
282  sit++;
283  }
284  return true;
285  }
286  default:
287  {
288  cerr << "MCPdf::Samplefrom(int, void *): No such sampling method" << endl;
289  return false;
290  }
291  }
292  }
293 
294  template <typename T> bool
295  MCPdf<T>::SampleFrom(Sample<T>& one_sample, int method, void * args) const
296  {
297  switch(method)
298  {
299  case DEFAULT:
300  {
301  // Sample from univariate uniform rng between 0 and 1;
302  double unif_sample; unif_sample = runif();
303  // Compare where we should be:
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();
310 
311  while ( unif_sample > *CumPDFit )
312  {
313  // check for internal error
314  assert(index <= size);
315  index++; it++; CumPDFit++;
316  }
317  it--;
318  one_sample = *it;
319  return true;
320  }
321  default:
322  {
323  cerr << "MCPdf::Samplefrom(int, void *): No such sampling method" << endl;
324  return false;
325  }
326  }
327  }
328 
329 
330  template <typename T> unsigned int MCPdf<T>::NumSamplesGet() const
331  {
332  return _listOfSamples.size();
333  }
334 
335  template <typename T> const WeightedSample<T>&
336  MCPdf<T>::SampleGet(unsigned int i) const
337  {
338  assert(i < NumSamplesGet());
339  return _listOfSamples[i];
340  }
341 
342  // Get and set number of samples
343  template <typename T> void
344  MCPdf<T>::NumSamplesSet(unsigned int num_samples)
345  {
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;
354  if (size < ns) // Add samples
355  {
357  _listOfSamples.insert(_listOfSamples.end(),(ns - size),ws);
358  _CumPDF.insert(_CumPDF.end(),(ns - size),0.0);
359  }
360  else if (size > ns) // Delete some samples
361  {
362  it = _listOfSamples.begin(); CumPDFit = _CumPDF.begin();
363  for ( unsigned int index = 0; index < (size-ns); index++ )
364  {
365  it = _listOfSamples.erase(it);
366  CumPDFit = _CumPDF.erase(CumPDFit);
367  }
368 #ifdef __MCPDF_DEBUG__
369  cout << "MCPDF::NumSamplesSet: WARNING DELETING SAMPLES!!" << endl;
370 #endif // __MCPDF_DEBUG__
371  }
372  else {;} // Do nothing (number of samples are equal)
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__
377  }
378 
379 
380  // Get and set the list of samples
381  template <typename T> bool
382  MCPdf<T>::ListOfSamplesSet(const vector<WeightedSample<T> > & los)
383  {
384  // Allocate necessary memory
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();
391  }
392 
393 
394  template <typename T> bool
395  MCPdf<T>::ListOfSamplesSet(const vector<Sample<T> > & los)
396  {
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;
400  // Allocate necessary memory
401  this->NumSamplesSet(numsamples);
402  // Update the list of samples
403  for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ )
404  {
405  *it = *lit; ;
406  it->WeightSet(1.0/numsamples);
407  lit++;
408  }
409  _SumWeights = 1.0;
410  // Update Cum PDF
411  this->CumPDFUpdate();
412 
413 #ifdef __MCPDF_DEBUG__
414  cout << "MCPDF ListOfSamplesSet: NumSamples = " << _listOfSamples.size()
415  << " SumWeights = " << _SumWeights << endl;
416 #endif // __MCPDF_DEBUG__
417 
418  return true;
419  }
420 
421  template <typename T> const vector<WeightedSample<T> > &
423  {
424  return _listOfSamples;
425  }
426 
427 
428  template <typename T> bool
430  {
431  assert (los.size() == _listOfSamples.size());
432  if (los.size() != 0)
433  {
434  _listOfSamples = los;
435  return this->NormalizeWeights();
436  }
437  return true;
438  }
439 
440  template <typename T> bool
441  MCPdf<T>::ListOfSamplesUpdate(const vector<Sample<T> > & los)
442  {
443  unsigned int numsamples = _listOfSamples.size();
444  if ((numsamples = los.size()) == _listOfSamples.size())
445  {
446  assert (numsamples != 0);
447  typename vector<Sample<T> >::const_iterator lit; lit=los.begin();
448  static typename vector<WeightedSample<T> >::iterator it;
449  // Allocate necessary memory
450  this->NumSamplesSet(numsamples);
451  // Update the sumweights
452  for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ )
453  {
454  *it = *lit; ;
455  it->WeightSet(1.0/numsamples);
456  lit++;
457  }
458  _SumWeights = 1.0;
459  this->CumPDFUpdate();
460  }
461  return true;
462  }
463 
464 
465  template <typename T> bool
467  {
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++ )
472  {
473  current_weight = it->WeightGet();
474  SumOfWeights += current_weight;
475  }
476 
477 #ifdef __MCPDF_DEBUG__
478  cout << "MCPDF::SumWeightsUpdate: SumWeights = " << SumOfWeights << endl;
479 #endif // __MCPDF_DEBUG__
480 
481  if (SumOfWeights > 0){
482  this->_SumWeights = SumOfWeights;
483  return true;
484  }
485  else{
486  cerr << "MCPDF::SumWeightsUpdate: SumWeights = " << SumOfWeights << endl;
487  return false;
488  }
489  }
490 
491  template <typename T> bool
493  {
494  static typename vector<WeightedSample<T> >::iterator it;
495 
496  // if sumweights = 0, something is wrong
497  if (!this->SumWeightsUpdate()) return false;
498 
499  for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ )
500  {
501  it->WeightSet(it->WeightGet() / _SumWeights);
502  }
503  this->_SumWeights = 1.0;
504  this->CumPDFUpdate();
505  return true;
506  }
507 
508 
509  template <typename T> void
511  {
512  double CumSum=0.0;
513  static typename vector<double>::iterator CumPDFit;
514  static typename vector<WeightedSample<T> >::iterator it;
515  CumPDFit = _CumPDF.begin(); *CumPDFit = 0.0;
516 
517  // Calculate the Cumulative PDF
518  for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ )
519  {
520  CumPDFit++;
521  // Calculate the __normalised__ Cumulative sum!!!
522  CumSum += ( it->WeightGet() / _SumWeights);
523  *CumPDFit = CumSum;
524  }
525  }
526 
527 
528  template <typename T>
529  T MCPdf<T>::ExpectedValueGet ( ) const
530  {
531  cerr << "MCPDF ExpectedValueGet: not implemented for the template parameters you use."
532  << endl << "Use template specialization as shown in mcpdf.cpp " << endl;
533 
534  assert(0);
535  T result;
536  return result;
537  }
538 
539 
540  template <typename T>
541  MatrixWrapper::SymmetricMatrix MCPdf<T>::CovarianceGet ( ) const
542  {
543  cerr << "MCPDF CovarianceGet: not implemented for the template parameters you use."
544  << endl << "Use template specialization as shown in mcpdf.cpp " << endl;
545 
546  assert(0);
547  MatrixWrapper::SymmetricMatrix result;
548  return result;
549  }
550 
551 
552 
553  template <typename T>
554  vector<double> & MCPdf<T>::CumulativePDFGet()
555  {
556  return _CumPDF;
557  }
558 
559 
560 
561 } // End namespace BFL
562 
563 #include "mcpdf.cpp"
564 
565 #endif
virtual ~MCPdf()
destructor
Class PDF: Virtual Base class representing Probability Density Functions.
Definition: pdf.h:53
#define DEFAULT
virtual MCPdf< T > * Clone() const
Clone function.
unsigned int NumSamplesGet() const
Get number of samples.
vector< WeightedSample< T > >::iterator _it_los
Definition: mcpdf.h:81
bool ListOfSamplesUpdate(const vector< WeightedSample< T > > &list_of_samples)
Update the list of samples (overloaded)
vector< WeightedSample< T > > _los
Definition: mcpdf.h:76
bool SumWeightsUpdate()
STL-iterator for cumulative PDF list.
vector< WeightedSample< T > > _listOfSamples
STL-list containing the list of samples.
Definition: mcpdf.h:56
void NumSamplesSet(unsigned int num_samples)
Set number of samples.
T _diff
Definition: mcpdf.h:78
vector< double > & CumulativePDFGet()
Add a sample to the list.
const vector< WeightedSample< T > > & ListOfSamplesGet() const
Get the list of weighted samples.
Monte Carlo Pdf: Sample based implementation of Pdf.
Definition: mcpdf.h:49
const WeightedSample< T > & SampleGet(unsigned int i) const
Get one sample.
MatrixWrapper::SymmetricMatrix CovarianceGet() const
Get the Covariance Matrix E[(x - E[x])^2] of the Analytic pdf.
unsigned int DimensionGet() const
Get the dimension of the argument.
#define RIPLEY
Matrix _diffsum
Definition: mcpdf.h:80
virtual bool SampleFrom(vector< Sample< T > > &list_samples, const unsigned int num_samples, int method=DEFAULT, void *args=NULL) const
Draw multiple samples from the Pdf (overloaded)
vector< double > _CumPDF
STL-iterator.
Definition: mcpdf.h:60
T _mean
Definition: mcpdf.h:77
void CumPDFUpdate()
After updating weights, we have to update the cumPDF.
T ExpectedValueGet() const
Get the expected value E[x] of the pdf.
double runif()
bool NormalizeWeights()
Normalizing the weights.
bool SampleFrom(Sample< T > &one_sample, int method=DEFAULT, void *args=NULL) const
double _SumWeights
Sum of all weights: used for normalising purposes.
Definition: mcpdf.h:54
T _CumSum
Definition: mcpdf.h:75
bool ListOfSamplesSet(const vector< WeightedSample< T > > &list_of_samples)
Set the list of weighted samples.
MCPdf(unsigned int num_samples=0, unsigned int dimension=0)
Constructor.
SymmetricMatrix _covariance
Definition: mcpdf.h:79


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 Mon Feb 28 2022 21:56:33