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 "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
00044
00048 template <typename T> class Mixture : public Pdf<T>
00049 {
00050 protected:
00052 unsigned int _numComponents;
00053
00055 vector<Probability> *_componentWeights;
00056
00058 vector<Pdf<T>* > *_componentPdfs;
00059
00061 bool NormalizeWeights();
00062
00064
00065
00066 vector<double> _cumWeights;
00067
00069 bool CumWeightsUpdate();
00070
00072
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
00144 int MostProbableComponentGet() const;
00145
00147
00148
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
00178
00179
00180
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
00188 _componentWeights = new vector<Probability>(this->NumComponentsGet());
00189
00190
00191 _componentPdfs = new vector< Pdf<T>* >(NumComponentsGet());
00192 #ifdef __CONSTRUCTOR__
00193 cout << "Mixture constructor\n";
00194 #endif // __CONSTRUCTOR__
00195 }
00196
00197
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
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
00212 _componentPdfs = new vector< Pdf<T>* >(NumComponentsGet());
00213
00214 for (int i=0; i<NumComponentsGet();i++)
00215 {
00216
00217
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
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
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
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:
00295 return Pdf<T>::SampleFrom(list_samples, num_samples,method,args);
00296
00297 case RIPLEY:
00298 {
00299 list_samples.resize(num_samples);
00300
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
00306 unif_samples[num_samples-1] = pow(unif_samples[num_samples-1],
00307 double (1.0/num_samples));
00308
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
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
00323 assert(index <= num_states);
00324 index++; CumPDFit++;
00325 }
00326
00327
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
00347 double unif_sample; unif_sample = runif();
00348
00349 unsigned int index = 0;
00350 while ( unif_sample > _cumWeights[index] )
00351 {
00352 assert(index <= NumComponentsGet());
00353 index++;
00354 }
00355
00356
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
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
00428
00429
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
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
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
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
00523 assert(NumComponentsGet()==(*_componentPdfs).size());
00524 assert(NumComponentsGet()==(*_componentWeights).size());
00525 assert(NumComponentsGet()+1==_cumWeights.size());
00526 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
00580 double CumSum=0.0;
00581 static vector<double>::iterator CumWeightsit;
00582 CumWeightsit = _cumWeights.begin();
00583 *CumWeightsit = 0.0;
00584
00585
00586 for ( unsigned int i = 0; i < NumComponentsGet(); i++)
00587 {
00588 CumWeightsit++;
00589
00590 CumSum += ( (*_componentWeights)[i] );
00591 *CumWeightsit = CumSum;
00592 }
00593
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 }
00602
00603 #include "mixture.cpp"
00604
00605 #endif // MIXTURE_H