00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #ifndef MIXTURE_H
00031 #define MIXTURE_H
00032
00033 #include "pdf.h"
00034 #include "../wrappers/matrix/vector_wrapper.h"
00035 #include "../wrappers/matrix/matrix_wrapper.h"
00036 #include "../wrappers/rng/rng.h"
00037 #include <vector>
00038
00039 namespace BFL
00040 {
00042
00043
00047 template <typename T> class Mixture : public Pdf<T>
00048 {
00049 protected:
00051 unsigned int _numComponents;
00052
00054 vector<Probability> *_componentWeights;
00055
00057 vector<Pdf<T>* > *_componentPdfs;
00058
00060 bool NormalizeWeights();
00061
00063
00064
00065 vector<double> _cumWeights;
00066
00068 bool CumWeightsUpdate();
00069
00071
00072 void TestNotInit() const;
00073
00074 public:
00076
00078 Mixture(unsigned int dimension=0);
00079
00081
00083 template <typename U> Mixture(U &componentVector);
00084
00086
00087 Mixture(const Mixture & my_mixture);
00088
00090 virtual ~Mixture();
00091
00093 virtual Mixture* Clone() const;
00094
00096 unsigned int NumComponentsGet()const;
00097
00099 Probability ProbabilityGet(const T& state) const;
00100
00101 bool SampleFrom (vector<Sample<T> >& list_samples,
00102 const unsigned int num_samples,
00103 int method = DEFAULT,
00104 void * args = NULL) const;
00105
00106 bool SampleFrom (Sample<T>& one_sample, int method = DEFAULT, void * args = NULL) const;
00107
00108 T ExpectedValueGet() const;
00109
00110 MatrixWrapper::SymmetricMatrix CovarianceGet ( ) const;
00111
00113 vector<Probability> WeightsGet() const;
00114
00116
00119 Probability WeightGet(unsigned int componentNumber) const;
00120
00122
00127 bool WeightsSet(vector<Probability> & weights);
00128
00130
00139 bool WeightSet(unsigned int componentNumber, Probability w);
00140
00142
00143 int MostProbableComponentGet() const;
00144
00146
00147
00150 bool AddComponent(Pdf<T>& pdf );
00151
00153
00156 bool AddComponent(Pdf<T>& pdf, Probability w);
00157
00159
00162 bool DeleteComponent(unsigned int componentNumber );
00163
00165 vector<Pdf<T>* > ComponentsGet() const;
00166
00168
00171 Pdf<T>* ComponentGet(unsigned int componentNumber) const;
00172 };
00173
00175
00177
00178
00179
00180 template<typename T>
00181 Mixture<T>::Mixture( unsigned int dimension):
00182 Pdf<T>(dimension)
00183 , _numComponents(0)
00184 , _cumWeights(_numComponents+1)
00185 {
00186
00187 _componentWeights = new vector<Probability>(this->NumComponentsGet());
00188
00189
00190 _componentPdfs = new vector< Pdf<T>* >(NumComponentsGet());
00191 #ifdef __CONSTRUCTOR__
00192 cout << "Mixture constructor\n";
00193 #endif // __CONSTRUCTOR__
00194 }
00195
00196
00197 template<typename T> template <typename U>
00198 Mixture<T>::Mixture(U &componentVector): Pdf<T>(componentVector[0]->DimensionGet() )
00199 , _numComponents(componentVector.size())
00200 {
00201
00202 _componentWeights = new vector<Probability>(NumComponentsGet());
00203 for (int i=0; i<NumComponentsGet();i++)
00204 {
00205 (*_componentWeights)[i] = (Probability)(1.0/NumComponentsGet());
00206 }
00207 _cumWeights.insert(_cumWeights.begin(),NumComponentsGet()+1,0.0);
00208 CumWeightsUpdate();
00209
00210
00211 _componentPdfs = new vector< Pdf<T>* >(NumComponentsGet());
00212
00213 for (int i=0; i<NumComponentsGet();i++)
00214 {
00215
00216
00217 (*_componentPdfs)[i] = (componentVector[i])->Clone();
00218 }
00219 #ifdef __CONSTRUCTOR__
00220 cout << "Mixture constructor\n";
00221 #endif // __CONSTRUCTOR__
00222 }
00223
00224 template<typename T>
00225 Mixture<T>::Mixture(const Mixture & my_mixture):Pdf<T>(my_mixture.DimensionGet() )
00226 ,_numComponents(my_mixture.NumComponentsGet())
00227 {
00228
00229 _componentWeights = new vector<Probability>(this->NumComponentsGet());
00230 (*_componentWeights) = my_mixture.WeightsGet();
00231 _cumWeights.insert(_cumWeights.begin(),NumComponentsGet()+1,0.0);
00232 CumWeightsUpdate();
00233
00234
00235 _componentPdfs = new vector< Pdf<T>* >(NumComponentsGet());
00236
00237
00238 for (int i=0; i<NumComponentsGet();i++)
00239 {
00240
00241
00242
00243
00244 (*_componentPdfs)[i] = (my_mixture.ComponentGet(i))->Clone();
00245 }
00246 #ifdef __CONSTRUCTOR__
00247 cout << "Mixture copy constructor\n";
00248 #endif // __CONSTRUCTOR__
00249 }
00250
00251 template<typename T>
00252 Mixture<T>::~Mixture()
00253 {
00254 #ifdef __CONSTRUCTOR__
00255 cout << "Mixture destructor\n";
00256 #endif
00257
00258 delete _componentWeights;
00259 for (int i=0; i<NumComponentsGet();i++)
00260 {
00261 delete (*_componentPdfs)[i];
00262 }
00263 delete _componentPdfs;
00264 }
00265
00266 template<typename T>
00267 Mixture<T>* Mixture<T>::Clone() const
00268 {
00269 return new Mixture(*this);
00270 }
00271
00272 template<typename T>
00273 unsigned int Mixture<T>::NumComponentsGet() const
00274 {
00275 return _numComponents;
00276 }
00277
00278 template<typename T>
00279 Probability Mixture<T>::ProbabilityGet(const T& state) const
00280 {
00281 TestNotInit();
00282 Probability prob(0.0);
00283 for (int i=0; i<NumComponentsGet();i++)
00284 {
00285 prob= prob + (*_componentWeights)[i] * (*_componentPdfs)[i]->ProbabilityGet(state);
00286 }
00287 return prob;
00288 }
00289
00290 template<typename T>
00291 bool Mixture<T>::SampleFrom (vector<Sample<T> >& list_samples,
00292 const unsigned int num_samples,
00293 int method,
00294 void * args) const
00295 {
00296 TestNotInit();
00297 switch(method)
00298 {
00299 case DEFAULT:
00300 return Pdf<T>::SampleFrom(list_samples, num_samples,method,args);
00301
00302 case RIPLEY:
00303 {
00304 list_samples.resize(num_samples);
00305
00306 std::vector<double> unif_samples(num_samples);
00307 for ( unsigned int i = 0; i < num_samples ; i++)
00308 unif_samples[i] = runif();
00309
00310
00311 unif_samples[num_samples-1] = pow(unif_samples[num_samples-1],
00312 double (1.0/num_samples));
00313
00314 for ( int i = num_samples-2; i >= 0 ; i--)
00315 unif_samples[i] = pow(unif_samples[i], double (1.0/(i+1))) * unif_samples[i+1];
00316
00317
00318 unsigned int index = 0;
00319 unsigned int num_states = NumComponentsGet();
00320 vector<double>::const_iterator CumPDFit = _cumWeights.begin();
00321 typename vector<Sample<T> >::iterator sit = list_samples.begin();
00322
00323 for ( unsigned int i = 0; i < num_samples ; i++)
00324 {
00325 while ( unif_samples[i] > *CumPDFit )
00326 {
00327
00328 assert(index <= num_states);
00329 index++; CumPDFit++;
00330 }
00331
00332
00333 (*_componentPdfs)[index-1]->SampleFrom(*sit,method,args);
00334 sit++;
00335 }
00336 return true;
00337 }
00338 default:
00339 cerr << "Mixture::Samplefrom(T, void *): No such sampling method" << endl;
00340 return false;
00341 }
00342 }
00343 template<typename T>
00344 bool Mixture<T>::SampleFrom (Sample<T>& one_sample, int method, void * args) const
00345 {
00346 TestNotInit();
00347 switch(method)
00348 {
00349 case DEFAULT:
00350 {
00351
00352 double unif_sample; unif_sample = runif();
00353
00354 unsigned int index = 0;
00355 while ( unif_sample > _cumWeights[index] )
00356 {
00357 assert(index <= NumComponentsGet());
00358 index++;
00359 }
00360
00361
00362 (*_componentPdfs)[index-1]->SampleFrom(one_sample,method,args);
00363 return true;
00364 }
00365 default:
00366 cerr << "Mixture::Samplefrom(T, void *): No such sampling method"
00367 << endl;
00368 return false;
00369 }
00370 }
00371
00372 template<typename T>
00373 T Mixture<T>::ExpectedValueGet() const
00374 {
00375 cerr << "Mixture ExpectedValueGet: not implemented for the template parameters you use."
00376 << endl << "Use template specialization as shown in mixture.cpp " << endl;
00377 assert(0);
00378 T expectedValue;
00379 return expectedValue;
00380 }
00381
00382 template <typename T>
00383 MatrixWrapper::SymmetricMatrix Mixture<T>::CovarianceGet ( ) const
00384 {
00385 TestNotInit();
00386 cerr << "Mixture CovarianceGet: not implemented since so far I don't believe its usefull"
00387 << endl << "If you decide to implement is: Use template specialization as shown in mcpdf.cpp " << endl;
00388
00389 assert(0);
00390 MatrixWrapper::SymmetricMatrix result;
00391 return result;
00392 }
00393
00394 template<typename T>
00395 vector<Probability> Mixture<T>::WeightsGet() const
00396 {
00397 TestNotInit();
00398 return *_componentWeights;
00399 }
00400
00401 template<typename T>
00402 Probability Mixture<T>::WeightGet(unsigned int componentNumber) const
00403 {
00404 TestNotInit();
00405 assert((int)componentNumber >= 0 && componentNumber < NumComponentsGet());
00406 return (*_componentWeights)[componentNumber];
00407 }
00408
00409 template<typename T>
00410 bool Mixture<T>::WeightsSet(vector<Probability> & weights)
00411 {
00412 TestNotInit();
00413 assert(weights.size() == NumComponentsGet());
00414 *_componentWeights = weights;
00415
00416 return (NormalizeWeights() && CumWeightsUpdate());
00417 }
00418
00419 template<typename T>
00420 bool Mixture<T>::WeightSet(unsigned int componentNumber, Probability weight)
00421 {
00422 TestNotInit();
00423 assert((int)componentNumber >= 0 && componentNumber < NumComponentsGet());
00424 assert((double)weight<=1.0);
00425
00426 if (NumComponentsGet() == 1)
00427 {
00428 (*_componentWeights)[0] = weight;
00429 }
00430 else
00431 {
00432
00433
00434
00435 Probability old_weight = WeightGet(componentNumber);
00436 if ((double)old_weight!=1.0) {
00437 double normalization_factor = (1-weight)/(1-old_weight);
00438 for (int i=0; i<NumComponentsGet();i++)
00439 {
00440 (*_componentWeights)[i] = (Probability)( (double)( (*_componentWeights)[i] )* normalization_factor);
00441 }
00442 }
00443 else{
00444 for (int i=0; i<NumComponentsGet();i++)
00445 {
00446 (*_componentWeights)[i] = (Probability)( (1-weight)/(NumComponentsGet()-1));
00447 }
00448 }
00449 (*_componentWeights)[componentNumber] = weight;
00450 }
00451 return CumWeightsUpdate();
00452 }
00453
00454 template<typename T>
00455 int Mixture<T>::MostProbableComponentGet() const
00456 {
00457 TestNotInit();
00458 int index_mostProbable= -1;
00459 Probability prob_mostProbable= 0.0;
00460 for (int component = 0 ; component < NumComponentsGet() ; component++)
00461 {
00462 if ( (*_componentWeights)[component] > prob_mostProbable)
00463 {
00464 index_mostProbable= component;
00465 prob_mostProbable= (*_componentWeights)[component];
00466 }
00467 }
00468 return index_mostProbable;
00469 }
00470
00471 template<typename T>
00472 bool Mixture<T>::AddComponent(Pdf<T>& pdf)
00473 {
00474 if (NumComponentsGet()==0)
00475 return AddComponent(pdf, Probability(1.0));
00476 else
00477 {
00478 _numComponents++;
00479 (*_componentPdfs).push_back(pdf.Clone() );
00480
00481 (*_componentWeights).push_back(Probability(0.0));
00482 _cumWeights.push_back(0.0);
00483
00484 assert(NumComponentsGet()==(*_componentPdfs).size());
00485 assert(NumComponentsGet()==(*_componentWeights).size());
00486 assert(NumComponentsGet()+1==_cumWeights.size());
00487 return (NormalizeWeights() && CumWeightsUpdate());
00488 }
00489 }
00490
00491 template<typename T>
00492 bool Mixture<T>::AddComponent(Pdf<T>& pdf, Probability w)
00493 {
00494 if (NumComponentsGet()==0 && w!=1.0)
00495 return AddComponent(pdf, Probability(1.0));
00496 else
00497 {
00498 _numComponents++;
00499 (*_componentPdfs).push_back(pdf.Clone() );
00500 (*_componentWeights).push_back(Probability(0.0));
00501 _cumWeights.resize(NumComponentsGet()+1);
00502
00503 assert(NumComponentsGet()==(*_componentPdfs).size());
00504 assert(NumComponentsGet()==(*_componentWeights).size());
00505 assert(NumComponentsGet()+1==_cumWeights.size());
00506 WeightSet(_numComponents-1,w);
00507 return (NormalizeWeights() && CumWeightsUpdate());
00508 }
00509 }
00510
00511 template<typename T>
00512 bool Mixture<T>::DeleteComponent(unsigned int componentNumber)
00513 {
00514
00515 assert(NumComponentsGet()==(*_componentPdfs).size());
00516 assert(NumComponentsGet()==(*_componentWeights).size());
00517 assert(NumComponentsGet()+1==_cumWeights.size());
00518
00519 TestNotInit();
00520 assert((int)componentNumber >= 0 && componentNumber < NumComponentsGet());
00521 _numComponents--;
00522 Pdf<T>* pointer = (*_componentPdfs)[componentNumber];
00523 delete pointer;
00524 (*_componentPdfs).erase((*_componentPdfs).begin()+componentNumber);
00525 (*_componentWeights).erase((*_componentWeights).begin()+componentNumber);
00526 _cumWeights.resize(NumComponentsGet()+1);
00527
00528 assert(NumComponentsGet()==(*_componentPdfs).size());
00529 assert(NumComponentsGet()==(*_componentWeights).size());
00530 assert(NumComponentsGet()+1==_cumWeights.size());
00531 if(_numComponents==0)
00532 return true;
00533 else
00534 return (NormalizeWeights() && CumWeightsUpdate());
00535 }
00536
00537 template<typename T>
00538 vector<Pdf<T>*> Mixture<T>::ComponentsGet() const
00539 {
00540 TestNotInit();
00541 return _componentPdfs;
00542 }
00543
00544 template<typename T>
00545 Pdf<T>* Mixture<T>::ComponentGet(unsigned int componentNumber) const
00546 {
00547 TestNotInit();
00548 return (*_componentPdfs)[componentNumber];
00549 }
00550
00551 template<typename T>
00552 void Mixture<T>::TestNotInit() const
00553 {
00554 if (NumComponentsGet() == 0)
00555 {
00556 cerr << "Mixture method called which requires that the number of components is not zero"
00557 << endl << "Current number of components: " << NumComponentsGet() << endl;
00558 assert(0);
00559 }
00560 }
00561
00562 template<typename T>
00563 bool Mixture<T>::NormalizeWeights()
00564 {
00565 double SumOfWeights = 0.0;
00566 for ( unsigned int i = 0; i < NumComponentsGet() ; i++){
00567 SumOfWeights += (*_componentWeights)[i];
00568 }
00569 if (SumOfWeights > 0){
00570 for ( unsigned int i = 0; i < NumComponentsGet() ; i++){
00571 (*_componentWeights)[i] = (Probability)( (double) ( (*_componentWeights)[i]) /SumOfWeights);
00572 }
00573 return true;
00574 }
00575 else{
00576 cerr << "Mixture::NormalizeProbs(): SumOfWeights = " << SumOfWeights << endl;
00577 return false;
00578 }
00579 }
00580
00581 template<typename T>
00582 bool Mixture<T>::CumWeightsUpdate()
00583 {
00584
00585 double CumSum=0.0;
00586 static vector<double>::iterator CumWeightsit;
00587 CumWeightsit = _cumWeights.begin();
00588 *CumWeightsit = 0.0;
00589
00590
00591 for ( unsigned int i = 0; i < NumComponentsGet(); i++)
00592 {
00593 CumWeightsit++;
00594
00595 CumSum += ( (*_componentWeights)[i] );
00596 *CumWeightsit = CumSum;
00597 }
00598
00599 assert( (_cumWeights[NumComponentsGet()] >= 1.0 - NUMERIC_PRECISION) &&
00600 (_cumWeights[NumComponentsGet()] <= 1.0 + NUMERIC_PRECISION) );
00601
00602 _cumWeights[NumComponentsGet()]=1;
00603 return true;
00604 }
00605
00606 }
00607
00608 #include "mixture.cpp"
00609
00610 #endif // MIXTURE_H