Bayesian Filtering Library
Generated from SVN r
|
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