mixture.h
Go to the documentation of this file.
00001 // $Id: mixture.h 2009-01-21 tdelaet $
00002 // Copyright (C)  2009 Tinne De Laet <first dot last at mech dot kuleuven dot be>
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 MIXTURE_H
00031 #define MIXTURE_H
00032 
00033 #include "pdf.h"
00034 #include "discretepdf.h"
00035 #include "../wrappers/matrix/vector_wrapper.h"
00036 #include "../wrappers/matrix/matrix_wrapper.h"
00037 #include "../wrappers/rng/rng.h"
00038 #include <vector>
00039 
00040 namespace BFL
00041 {
00043   //kind of pdfs but they should all share the same state space i.e. they all
00044   //have to be of the same template type Pdf<T> 
00048   template <typename T> class Mixture : public Pdf<T> // inherit abstract_template_class
00049     {
00050     protected:
00052       unsigned int _numComponents;
00053 
00055       vector<Probability> *_componentWeights;
00056 
00058       vector<Pdf<T>* > *_componentPdfs;
00059 
00061       bool NormalizeWeights();
00062 
00064       // BEWARE: this vector has length _numComponents +1 (first element is
00065       // always 0, last element is always 1)
00066       vector<double> _cumWeights;
00067 
00069       bool CumWeightsUpdate();
00070 
00072       //requires number of components>0
00073       void TestNotInit() const;
00074 
00075     public:
00077 
00079       Mixture(const unsigned int dimension=0);
00080 
00082 
00084       template <typename U> Mixture(const U &componentVector);
00085 
00087 
00088       Mixture(const Mixture & my_mixture);
00089 
00091       virtual ~Mixture();
00092 
00094       virtual Mixture* Clone() const;
00095 
00097       unsigned int NumComponentsGet()const;
00098 
00100       Probability ProbabilityGet(const T& state) const;
00101 
00102       bool SampleFrom (vector<Sample<T> >& list_samples,
00103                        const unsigned int num_samples,
00104                        int method = DEFAULT, 
00105                        void * args = NULL) const;
00106 
00107       bool SampleFrom (Sample<T>& one_sample, int method = DEFAULT, void * args = NULL) const;
00108 
00109       T ExpectedValueGet() const;
00110 
00111       MatrixWrapper::SymmetricMatrix CovarianceGet (  ) const;
00112 
00114       vector<Probability> WeightsGet() const;
00115 
00117 
00120       Probability WeightGet(unsigned int componentNumber) const;
00121 
00123 
00128       bool WeightsSet(vector<Probability> & weights);
00129 
00131 
00140       bool WeightSet(unsigned int componentNumber, Probability w);
00141 
00143       //equally probable, the index of the first most probable one is returned
00144       int MostProbableComponentGet() const;
00145 
00147        // the weight of the new component is set to 0 (except when the number of
00148        // components is zero, then the weight is set to 1)!!
00151       bool AddComponent(Pdf<T>& pdf );
00152 
00154 
00157       bool AddComponent(Pdf<T>& pdf, Probability w);
00158 
00160 
00163       bool DeleteComponent(unsigned int componentNumber );
00164 
00166       vector<Pdf<T>* > ComponentsGet() const;
00167 
00169 
00172       Pdf<T>*  ComponentGet(unsigned int componentNumber) const;
00173     };
00174 
00176 // Template Code here
00178 
00179 // Constructor
00180 //TODO: is this usefull because pointers to components point to nothing!
00181 template<typename T>
00182 Mixture<T>::Mixture(const unsigned int dimension):
00183   Pdf<T>(dimension)
00184   , _numComponents(0)
00185   , _cumWeights(_numComponents+1)
00186   {
00187     //create pointer to vector of component weights
00188     _componentWeights = new vector<Probability>(this->NumComponentsGet());
00189 
00190     //create pointer to vector of pointers to pdfs
00191     _componentPdfs = new vector< Pdf<T>* >(NumComponentsGet());
00192 #ifdef __CONSTRUCTOR__
00193     cout << "Mixture constructor\n";
00194 #endif // __CONSTRUCTOR__
00195   }
00196 
00197 // Constructor
00198 template<typename T> template <typename U>
00199 Mixture<T>::Mixture(const U &componentVector): Pdf<T>(componentVector[0]->DimensionGet() )
00200     , _numComponents(componentVector.size())
00201 {
00202     //create pointer to vector of component weights
00203     _componentWeights = new vector<Probability>(NumComponentsGet());
00204     for (int i=0; i<NumComponentsGet();i++)
00205     {
00206         (*_componentWeights)[i] = (Probability)(1.0/NumComponentsGet());
00207     }
00208     _cumWeights.insert(_cumWeights.begin(),NumComponentsGet()+1,0.0);
00209     CumWeightsUpdate();
00210 
00211     //create pointer to vector of pointers to weights
00212     _componentPdfs = new vector< Pdf<T>* >(NumComponentsGet());
00213     //TODO: take copy or point to same???
00214     for (int i=0; i<NumComponentsGet();i++)
00215     {
00216         //TODO: will this call the constructor of e.g. Gaussian or always the
00217         //general one?
00218         (*_componentPdfs)[i] = (componentVector[i])->Clone();
00219     }
00220 #ifdef __CONSTRUCTOR__
00221     cout << "Mixture constructor\n";
00222 #endif // __CONSTRUCTOR__
00223 }
00224 
00225 template<typename T > 
00226 Mixture<T>::Mixture(const Mixture & my_mixture):Pdf<T>(my_mixture.DimensionGet() )
00227         ,_numComponents(my_mixture.NumComponentsGet())
00228   {
00229     //create pointer to vector of component weights
00230     _componentWeights = new vector<Probability>(this->NumComponentsGet());
00231     (*_componentWeights) = my_mixture.WeightsGet();
00232     _cumWeights.insert(_cumWeights.begin(),NumComponentsGet()+1,0.0);
00233     CumWeightsUpdate();
00234 
00235     //create pointer to vector of pointers to weights
00236     _componentPdfs = new vector< Pdf<T>* >(NumComponentsGet());
00237     for (int i=0; i<NumComponentsGet();i++)
00238     {
00239         (*_componentPdfs)[i] = (my_mixture.ComponentGet(i))->Clone();
00240     }
00241 #ifdef __CONSTRUCTOR__
00242     cout << "Mixture copy constructor\n";
00243 #endif // __CONSTRUCTOR__
00244   }
00245 
00246 template<typename T>
00247   Mixture<T>::~Mixture()
00248   {
00249 #ifdef __CONSTRUCTOR__
00250     cout << "Mixture destructor\n";
00251 #endif
00252     // Release memory!
00253     delete _componentWeights;
00254     for (int i=0; i<NumComponentsGet();i++)
00255     {
00256         delete (*_componentPdfs)[i];
00257     }
00258     delete _componentPdfs;
00259   }
00260 
00261 template<typename T>
00262   Mixture<T>* Mixture<T>::Clone() const
00263   {
00264       return new Mixture(*this);
00265   }
00266 
00267 template<typename T>
00268   unsigned int Mixture<T>::NumComponentsGet() const
00269   {
00270     return _numComponents;
00271   }
00272 
00273 template<typename T>
00274   Probability Mixture<T>::ProbabilityGet(const T& state) const
00275   {
00276     TestNotInit(); 
00277     Probability prob(0.0);
00278     for (int i=0; i<NumComponentsGet();i++)
00279     {
00280         prob= prob + (*_componentWeights)[i] * (*_componentPdfs)[i]->ProbabilityGet(state);
00281     }
00282     return prob;
00283   }
00284 
00285 template<typename T>
00286   bool Mixture<T>::SampleFrom (vector<Sample<T> >& list_samples,
00287                            const unsigned int num_samples,
00288                            int method,
00289                            void * args) const
00290   {
00291     TestNotInit(); 
00292     switch(method)
00293     {
00294       case DEFAULT: // O(N log(N) efficiency)
00295           return Pdf<T>::SampleFrom(list_samples, num_samples,method,args);
00296 
00297       case RIPLEY: // See mcpdf.cpp for more explanation
00298           {
00299             list_samples.resize(num_samples);
00300             // GENERATE N ORDERED IID UNIFORM SAMPLES
00301             std::vector<double> unif_samples(num_samples);
00302             for ( unsigned int i = 0; i < num_samples ; i++)
00303               unif_samples[i] = runif();
00304 
00305             /* take n-th racine of u_N */
00306             unif_samples[num_samples-1] = pow(unif_samples[num_samples-1],
00307                                            double (1.0/num_samples));
00308             /* rescale samples */
00309             for ( int i = num_samples-2; i >= 0 ; i--)
00310                 unif_samples[i] = pow(unif_samples[i], double (1.0/(i+1))) * unif_samples[i+1];
00311 
00312             // CHECK WHERE THESE SAMPLES ARE IN _CUMWEIGHTS
00313             unsigned int index = 0;
00314             unsigned int num_states = NumComponentsGet();
00315             vector<double>::const_iterator CumPDFit = _cumWeights.begin();
00316             typename vector<Sample<T> >::iterator sit = list_samples.begin();
00317 
00318             for ( unsigned int i = 0; i < num_samples ; i++)
00319               {
00320                 while ( unif_samples[i] > *CumPDFit )
00321                 {
00322                 // check for internal error
00323                 assert(index <= num_states);
00324                 index++; CumPDFit++;
00325                 }
00326           // index-1 is a sample of the discrete pdf of the mixture weights
00327           // get a sample from the index-1'th mixture component
00328               (*_componentPdfs)[index-1]->SampleFrom(*sit,method,args);
00329               sit++;
00330               }
00331             return true;
00332           }
00333       default:
00334             cerr << "Mixture::Samplefrom(T, void *): No such sampling method" << endl;
00335             return false;
00336     }
00337   }
00338 template<typename T>
00339   bool Mixture<T>::SampleFrom (Sample<T>& one_sample, int method, void * args) const
00340   {
00341     TestNotInit(); 
00342     switch(method)
00343       {
00344       case DEFAULT:
00345         {
00346           // Sample from univariate uniform rng between 0 and 1;
00347           double unif_sample; unif_sample = runif();
00348           // Compare where we should be
00349           unsigned int index = 0;
00350           while ( unif_sample > _cumWeights[index] )
00351             {
00352               assert(index <= NumComponentsGet());
00353               index++;
00354             }
00355          // index-1 is a sample of the discrete pdf of the mixture weights
00356       // get a sample from the index-1'th mixture component
00357           (*_componentPdfs)[index-1]->SampleFrom(one_sample,method,args);
00358           return true;
00359         }
00360       default:
00361         cerr << "Mixture::Samplefrom(T, void *): No such sampling method"
00362              << endl;
00363         return false;
00364       }
00365   }
00366 
00367 template<typename T>
00368   T Mixture<T>::ExpectedValueGet() const
00369   {
00370     cerr << "Mixture ExpectedValueGet: not implemented for the template parameters you use." 
00371          << endl << "Use template specialization as shown in mixture.cpp " << endl; 
00372     assert(0);
00373     T expectedValue;
00374     return expectedValue;
00375   }
00376 
00377 template <typename T>
00378   MatrixWrapper::SymmetricMatrix Mixture<T>::CovarianceGet (  ) const
00379   {
00380     TestNotInit(); 
00381     cerr << "Mixture CovarianceGet: not implemented since so far I don't believe its usefull"
00382      << endl << "If you decide to implement is: Use template specialization as shown in mcpdf.cpp " << endl; 
00383 
00384     assert(0);
00385     MatrixWrapper::SymmetricMatrix result;
00386     return result;
00387   }
00388 
00389 template<typename T>
00390   vector<Probability> Mixture<T>::WeightsGet() const
00391   {
00392     TestNotInit(); 
00393     return *_componentWeights;
00394    }
00395 
00396 template<typename T>
00397   Probability Mixture<T>::WeightGet(unsigned int componentNumber) const
00398   {
00399     TestNotInit(); 
00400     assert((int)componentNumber >= 0 && componentNumber < NumComponentsGet());
00401     return (*_componentWeights)[componentNumber];
00402    }
00403 
00404 template<typename T>
00405   bool Mixture<T>::WeightsSet(vector<Probability> & weights)
00406   {
00407     TestNotInit(); 
00408     assert(weights.size() == NumComponentsGet());
00409     *_componentWeights = weights;
00410     //normalize the probabilities and update the cumulative sum
00411     return (NormalizeWeights() && CumWeightsUpdate());
00412   }
00413 
00414 template<typename T>
00415   bool Mixture<T>::WeightSet(unsigned int componentNumber, Probability weight)
00416   {
00417     TestNotInit(); 
00418     assert((int)componentNumber >= 0 && componentNumber < NumComponentsGet());
00419     assert((double)weight<=1.0);
00420     
00421     if (NumComponentsGet() == 1)
00422     {
00423         (*_componentWeights)[0] = weight;
00424     }
00425     else
00426     {
00427         // renormalize other weights such that sum of probabilities will be
00428         // one after the weight of the component is set to weight
00429         // This should keep the probabilities normalized
00430         Probability old_weight = WeightGet(componentNumber);
00431         if ((double)old_weight!=1.0) {
00432             double normalization_factor = (1-weight)/(1-old_weight);
00433             for (int i=0; i<NumComponentsGet();i++)
00434             {
00435                 (*_componentWeights)[i] = (Probability)( (double)( (*_componentWeights)[i] )* normalization_factor);
00436             }
00437         }
00438         else{
00439             for (int i=0; i<NumComponentsGet();i++)
00440             {
00441                 (*_componentWeights)[i] = (Probability)( (1-weight)/(NumComponentsGet()-1));
00442             }
00443         }
00444         (*_componentWeights)[componentNumber] = weight;
00445     }
00446     return CumWeightsUpdate();
00447   }
00448 
00449 template<typename T>
00450   int Mixture<T>::MostProbableComponentGet() const
00451   {
00452     TestNotInit(); 
00453     int index_mostProbable= -1;
00454     Probability prob_mostProbable= 0.0;
00455     for (int component = 0 ; component < NumComponentsGet() ; component++)
00456     {
00457        if ( (*_componentWeights)[component] > prob_mostProbable)
00458        {
00459             index_mostProbable= component;
00460             prob_mostProbable= (*_componentWeights)[component];
00461        }
00462     }
00463     return index_mostProbable;
00464   }
00465 
00466 template<typename T>
00467   bool Mixture<T>::AddComponent(Pdf<T>& pdf) 
00468   {
00469     if (NumComponentsGet()==0)
00470         return AddComponent(pdf, Probability(1.0));
00471     else
00472     {
00473         _numComponents++;
00474         (*_componentPdfs).push_back(pdf.Clone() );
00475 
00476         (*_componentWeights).push_back(Probability(0.0));
00477             _cumWeights.push_back(0.0);
00478         //assert length of vectors
00479         assert(NumComponentsGet()==(*_componentPdfs).size());
00480         assert(NumComponentsGet()==(*_componentWeights).size());
00481         assert(NumComponentsGet()+1==_cumWeights.size());
00482         return (NormalizeWeights() && CumWeightsUpdate());
00483     }
00484   }
00485 
00486 template<typename T>
00487   bool Mixture<T>::AddComponent(Pdf<T>& pdf, Probability w) 
00488   {
00489     if (NumComponentsGet()==0 && w!=1.0)
00490         return AddComponent(pdf, Probability(1.0));
00491     else
00492     {
00493         _numComponents++;
00494         (*_componentPdfs).push_back(pdf.Clone() );
00495         (*_componentWeights).push_back(Probability(0.0));
00496         _cumWeights.resize(NumComponentsGet()+1);
00497         //assert length of vectors
00498         assert(NumComponentsGet()==(*_componentPdfs).size());
00499         assert(NumComponentsGet()==(*_componentWeights).size());
00500         assert(NumComponentsGet()+1==_cumWeights.size());
00501         WeightSet(_numComponents-1,w);
00502         return (NormalizeWeights() && CumWeightsUpdate());
00503      }
00504   }
00505 
00506 template<typename T>
00507   bool Mixture<T>::DeleteComponent(unsigned int componentNumber) 
00508   {
00509     //assert length of vectors
00510     assert(NumComponentsGet()==(*_componentPdfs).size());
00511     assert(NumComponentsGet()==(*_componentWeights).size());
00512     assert(NumComponentsGet()+1==_cumWeights.size());
00513 
00514     TestNotInit(); 
00515     assert((int)componentNumber >= 0 && componentNumber < NumComponentsGet());
00516     _numComponents--;
00517     Pdf<T>* pointer = (*_componentPdfs)[componentNumber];
00518     delete pointer;
00519     (*_componentPdfs).erase((*_componentPdfs).begin()+componentNumber);
00520     (*_componentWeights).erase((*_componentWeights).begin()+componentNumber);
00521     _cumWeights.resize(NumComponentsGet()+1);
00522     //assert length of vectors
00523     assert(NumComponentsGet()==(*_componentPdfs).size());
00524     assert(NumComponentsGet()==(*_componentWeights).size());
00525     assert(NumComponentsGet()+1==_cumWeights.size());
00526     if(_numComponents==0) //don't do normalization if numComponents == 0
00527         return true;
00528     else
00529         return (NormalizeWeights() && CumWeightsUpdate());
00530   }
00531 
00532 template<typename T>
00533   vector<Pdf<T>*> Mixture<T>::ComponentsGet() const
00534   {
00535     TestNotInit();
00536     return _componentPdfs;
00537   }
00538 
00539 template<typename T>
00540   Pdf<T>* Mixture<T>::ComponentGet(unsigned int componentNumber) const
00541   {
00542     TestNotInit();
00543     return (*_componentPdfs)[componentNumber];
00544   }
00545 
00546 template<typename T>
00547   void Mixture<T>::TestNotInit() const
00548   {
00549     if (NumComponentsGet() == 0)
00550     {
00551     cerr << "Mixture method called which requires that the number of components is not zero"
00552          << endl << "Current number of components: " << NumComponentsGet() << endl;
00553     assert(0);
00554     }
00555   }
00556 
00557 template<typename T>
00558   bool Mixture<T>::NormalizeWeights()
00559   {
00560     double SumOfWeights = 0.0;
00561     for ( unsigned int i = 0; i < NumComponentsGet() ; i++){
00562         SumOfWeights += (*_componentWeights)[i];
00563     }
00564     if (SumOfWeights > 0){
00565       for ( unsigned int i = 0; i < NumComponentsGet() ; i++){
00566           (*_componentWeights)[i] = (Probability)( (double) ( (*_componentWeights)[i]) /SumOfWeights);
00567       }
00568       return true;
00569     }
00570     else{
00571       cerr << "Mixture::NormalizeProbs(): SumOfWeights = " << SumOfWeights << endl;
00572       return false;
00573     }
00574   }
00575 
00576 template<typename T>
00577   bool Mixture<T>::CumWeightsUpdate()
00578   {
00579     // precondition: sum of probabilities should be 1
00580     double CumSum=0.0;
00581     static vector<double>::iterator CumWeightsit;
00582     CumWeightsit = _cumWeights.begin();
00583     *CumWeightsit = 0.0;
00584 
00585     // Calculate the Cumulative PDF
00586     for ( unsigned int i = 0; i < NumComponentsGet(); i++)
00587     {
00588             CumWeightsit++;
00589             // Calculate the __normalised__ Cumulative sum!!!
00590             CumSum += ( (*_componentWeights)[i] );
00591             *CumWeightsit = CumSum;
00592     }
00593     // Check if last element of valuelist is +- 1
00594     assert( (_cumWeights[NumComponentsGet()] >= 1.0 - NUMERIC_PRECISION) &&
00595             (_cumWeights[NumComponentsGet()] <= 1.0 + NUMERIC_PRECISION) );
00596 
00597     _cumWeights[NumComponentsGet()]=1;
00598     return true;
00599   }
00600 
00601 } // End namespace
00602 
00603 #include "mixture.cpp"
00604 
00605 #endif // MIXTURE_H


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