mcpdf.h
Go to the documentation of this file.
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


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 Sun Oct 5 2014 22:29:53