Bayesian Filtering Library Generated from SVN r
|
00001 // $Id: mcpdf.h 30606 2009-10-02 10:01:02Z tdelaet $ 00002 // Copyright (C) 2002 Klaas Gadeyne <first dot last at gmail dot com> 00003 /*************************************************************************** 00004 * This library is free software; you can redistribute it and/or * 00005 * modify it under the terms of the GNU General Public * 00006 * License as published by the Free Software Foundation; * 00007 * version 2 of the License. * 00008 * * 00009 * As a special exception, you may use this file as part of a free * 00010 * software library without restriction. Specifically, if other files * 00011 * instantiate templates or use macros or inline functions from this * 00012 * file, or you compile this file and link it with other files to * 00013 * produce an executable, this file does not by itself cause the * 00014 * resulting executable to be covered by the GNU General Public * 00015 * License. This exception does not however invalidate any other * 00016 * reasons why the executable file might be covered by the GNU General * 00017 * Public License. * 00018 * * 00019 * This library is distributed in the hope that it will be useful, * 00020 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 00021 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * 00022 * Lesser General Public License for more details. * 00023 * * 00024 * You should have received a copy of the GNU General Public * 00025 * License along with this library; if not, write to the Free Software * 00026 * Foundation, Inc., 59 Temple Place, * 00027 * Suite 330, Boston, MA 02111-1307 USA * 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 // vector<WeightedSample<T> >::iterator _it; 00060 vector<double> _CumPDF; 00062 // typename vector<double>::iterator CumPDFit; 00063 00064 00066 bool SumWeightsUpdate(); 00068 bool NormalizeWeights(); 00070 void CumPDFUpdate(); 00071 00072 private: 00073 // Variables to avoid allocation during call of 00074 //expectedvalueget/covarianceget 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 // implemented virtual functions 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 // void SampleAdd(WeightedSample<T> sample); 00158 00160 00162 vector<double> & CumulativePDFGet(); 00163 00164 }; 00165 00167 // Template Code here 00169 00170 // Constructor 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 // if (num_samples > 0) 00185 cout << "MCPDF Constructor: NumSamples = " << _listOfSamples.size() 00186 << ", CumPDF Samples = " << _CumPDF.size() 00187 << ", _SumWeights = " << _SumWeights << endl; 00188 #endif // __CONSTRUCTOR__ 00189 } 00190 00191 00192 00193 // Destructor 00194 template <typename T> 00195 MCPdf<T>::~MCPdf() 00196 { 00197 #ifdef __DESTRUCTOR__ 00198 cout << "MCPDF::Destructor" << endl; 00199 #endif // __DESTRUCTOR__ 00200 } 00201 00202 // Copy constructor 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 //Clone function 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: // O(N log(N) efficiency) 00237 { 00238 return Pdf<T>::SampleFrom(list_samples, numsamples,method,args); 00239 } 00240 case RIPLEY: // Only possible here ( O(N) efficiency ) 00241 /* See 00242 @Book{ ripley87, 00243 author = {Ripley, Brian D.}, 00244 title = {Stochastic Simulation}, 00245 publisher = {John Wiley and Sons}, 00246 year = {1987}, 00247 annote = {ISBN 0271-6356, WBIB 1 519.245} 00248 } 00249 */ 00250 // GENERATE N ORDERED IID UNIFORM SAMPLES 00251 { 00252 std::vector<double> unif_samples(numsamples); 00253 for ( unsigned int i = 0; i < numsamples ; i++) 00254 unif_samples[i] = runif(); 00255 00256 /* take n-th racine of u_N */ 00257 unif_samples[numsamples-1] = pow(unif_samples[numsamples-1], double (1.0/numsamples)); 00258 /* rescale other samples */ 00259 // only resample if more than one sample 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 // CHECK WHERE THESE SAMPLES ARE IN _CUMPDF 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 // Sample from univariate uniform rng between 0 and 1; 00302 double unif_sample; unif_sample = runif(); 00303 // Compare where we should be: 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 // check for internal error 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 // Get and set number of samples 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) // Add samples 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) // Delete some samples 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 {;} // Do nothing (number of samples are equal) 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 // Get and set the list of samples 00381 template <typename T> bool 00382 MCPdf<T>::ListOfSamplesSet(const vector<WeightedSample<T> > & los) 00383 { 00384 // Allocate necessary memory 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 // Allocate necessary memory 00401 this->NumSamplesSet(numsamples); 00402 // Update the list of samples 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 // Update Cum PDF 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 // Allocate necessary memory 00450 this->NumSamplesSet(numsamples); 00451 // Update the sumweights 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 // if sumweights = 0, something is wrong 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 // Calculate the Cumulative PDF 00518 for ( it = _listOfSamples.begin() ; it != _listOfSamples.end() ; it++ ) 00519 { 00520 CumPDFit++; 00521 // Calculate the __normalised__ Cumulative sum!!! 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 } // End namespace BFL 00562 00563 #include "mcpdf.cpp" 00564 00565 #endif