Bayesian Filtering Library Generated from SVN r
mcpdf.h
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., 51 Franklin Street, Fifth Floor, *
27 * Boston, MA 02110-1301 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
41namespace BFL
42{
43
45
49 template <typename T> class MCPdf: public Pdf<T>
50 {
51 protected:
52
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, 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;
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;
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>*
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 const SampleMthd method,
231 void * args) const
232 {
233 list_samples.resize(numsamples);
234 switch(method)
235 {
236 case SampleMthd::DEFAULT: // O(N log(N) efficiency)
237 {
238 return Pdf<T>::SampleFrom(list_samples, numsamples,method,args);
239 }
240 case SampleMthd::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, const SampleMthd method, void * args) const
296 {
297 switch(method)
298 {
299 case SampleMthd::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
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
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
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>
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
Monte Carlo Pdf: Sample based implementation of Pdf.
Definition mcpdf.h:50
void NumSamplesSet(unsigned int num_samples)
Set number of samples.
Definition mcpdf.h:344
const WeightedSample< T > & SampleGet(unsigned int i) const
Get one sample.
Definition mcpdf.h:336
vector< double > _CumPDF
STL-iterator.
Definition mcpdf.h:60
double _SumWeights
Sum of all weights: used for normalising purposes.
Definition mcpdf.h:54
bool SumWeightsUpdate()
STL-iterator for cumulative PDF list.
Definition mcpdf.h:466
virtual ~MCPdf()
destructor
Definition mcpdf.h:195
MCPdf(const MCPdf< T > &)
copy constructor
bool SampleFrom(Sample< T > &one_sample, const SampleMthd method=SampleMthd::DEFAULT, void *args=NULL) const
Draw 1 sample from the Pdf:
Definition mcpdf.h:295
void CumPDFUpdate()
After updating weights, we have to update the cumPDF.
Definition mcpdf.h:510
bool ListOfSamplesUpdate(const vector< WeightedSample< T > > &list_of_samples)
Update the list of samples (overloaded)
Definition mcpdf.h:429
const vector< WeightedSample< T > > & ListOfSamplesGet() const
Get the list of weighted samples.
Definition mcpdf.h:422
bool ListOfSamplesSet(const vector< WeightedSample< T > > &list_of_samples)
Set the list of weighted samples.
Definition mcpdf.h:382
T ExpectedValueGet() const
Get the expected value E[x] of the pdf.
Definition mcpdf.h:529
MatrixWrapper::SymmetricMatrix CovarianceGet() const
Get the Covariance Matrix E[(x - E[x])^2] of the Analytic pdf.
Definition mcpdf.h:541
bool NormalizeWeights()
Normalizing the weights.
Definition mcpdf.h:492
virtual MCPdf< T > * Clone() const
Clone function.
Definition mcpdf.h:222
MCPdf(unsigned int num_samples=0, unsigned int dimension=0)
Constructor.
Definition mcpdf.h:171
vector< WeightedSample< T > > _listOfSamples
STL-list containing the list of samples.
Definition mcpdf.h:56
unsigned int NumSamplesGet() const
Get number of samples.
Definition mcpdf.h:330
vector< double > & CumulativePDFGet()
Add a sample to the list.
Definition mcpdf.h:554
Class PDF: Virtual Base class representing Probability Density Functions.
Definition pdf.h:51
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)
Definition pdf.h:179