mixture.h
Go to the documentation of this file.
1 // $Id: mixture.h 2009-01-21 tdelaet $
2 // Copyright (C) 2009 Tinne De Laet <first dot last at mech dot kuleuven dot be>
3  /***************************************************************************
4  * This library is free software; you can redistribute it and/or *
5  * modify it under the terms of the GNU General Public *
6  * License as published by the Free Software Foundation; *
7  * version 2 of the License. *
8  * *
9  * As a special exception, you may use this file as part of a free *
10  * software library without restriction. Specifically, if other files *
11  * instantiate templates or use macros or inline functions from this *
12  * file, or you compile this file and link it with other files to *
13  * produce an executable, this file does not by itself cause the *
14  * resulting executable to be covered by the GNU General Public *
15  * License. This exception does not however invalidate any other *
16  * reasons why the executable file might be covered by the GNU General *
17  * Public License. *
18  * *
19  * This library is distributed in the hope that it will be useful, *
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
22  * Lesser General Public License for more details. *
23  * *
24  * You should have received a copy of the GNU General Public *
25  * License along with this library; if not, write to the Free Software *
26  * Foundation, Inc., 59 Temple Place, *
27  * Suite 330, Boston, MA 02111-1307 USA *
28  * *
29  ***************************************************************************/
30 #ifndef MIXTURE_H
31 #define MIXTURE_H
32 
33 #include "pdf.h"
34 #include "discretepdf.h"
35 #include "../wrappers/matrix/vector_wrapper.h"
36 #include "../wrappers/matrix/matrix_wrapper.h"
37 #include "../wrappers/rng/rng.h"
38 #include <vector>
39 
40 namespace BFL
41 {
43  //kind of pdfs but they should all share the same state space i.e. they all
44  //have to be of the same template type Pdf<T>
48  template <typename T> class Mixture : public Pdf<T> // inherit abstract_template_class
49  {
50  protected:
52  unsigned int _numComponents;
53 
55  vector<Probability> *_componentWeights;
56 
58  vector<Pdf<T>* > *_componentPdfs;
59 
61  bool NormalizeWeights();
62 
64  // BEWARE: this vector has length _numComponents +1 (first element is
65  // always 0, last element is always 1)
66  vector<double> _cumWeights;
67 
69  bool CumWeightsUpdate();
70 
72  //requires number of components>0
73  void TestNotInit() const;
74 
75  public:
77 
79  Mixture(const unsigned int dimension=0);
80 
82 
84  template <typename U> Mixture(const U &componentVector);
85 
87 
88  Mixture(const Mixture & my_mixture);
89 
91  virtual ~Mixture();
92 
94  virtual Mixture* Clone() const;
95 
97  unsigned int NumComponentsGet()const;
98 
100  Probability ProbabilityGet(const T& state) const;
101 
102  bool SampleFrom (vector<Sample<T> >& list_samples,
103  const unsigned int num_samples,
104  int method = DEFAULT,
105  void * args = NULL) const;
106 
107  bool SampleFrom (Sample<T>& one_sample, int method = DEFAULT, void * args = NULL) const;
108 
109  T ExpectedValueGet() const;
110 
111  MatrixWrapper::SymmetricMatrix CovarianceGet ( ) const;
112 
114  vector<Probability> WeightsGet() const;
115 
117 
120  Probability WeightGet(unsigned int componentNumber) const;
121 
123 
128  bool WeightsSet(vector<Probability> & weights);
129 
131 
140  bool WeightSet(unsigned int componentNumber, Probability w);
141 
143  //equally probable, the index of the first most probable one is returned
144  int MostProbableComponentGet() const;
145 
147  // the weight of the new component is set to 0 (except when the number of
148  // components is zero, then the weight is set to 1)!!
151  bool AddComponent(Pdf<T>& pdf );
152 
154 
157  bool AddComponent(Pdf<T>& pdf, Probability w);
158 
160 
163  bool DeleteComponent(unsigned int componentNumber );
164 
166  vector<Pdf<T>* > ComponentsGet() const;
167 
169 
172  Pdf<T>* ComponentGet(unsigned int componentNumber) const;
173  };
174 
176 // Template Code here
178 
179 // Constructor
180 //TODO: is this usefull because pointers to components point to nothing!
181 template<typename T>
182 Mixture<T>::Mixture(const unsigned int dimension):
183  Pdf<T>(dimension)
184  , _numComponents(0)
186  {
187  //create pointer to vector of component weights
188  _componentWeights = new vector<Probability>(this->NumComponentsGet());
189 
190  //create pointer to vector of pointers to pdfs
191  _componentPdfs = new vector< Pdf<T>* >(NumComponentsGet());
192 #ifdef __CONSTRUCTOR__
193  cout << "Mixture constructor\n";
194 #endif // __CONSTRUCTOR__
195  }
196 
197 // Constructor
198 template<typename T> template <typename U>
199 Mixture<T>::Mixture(const U &componentVector): Pdf<T>(componentVector[0]->DimensionGet() )
200  , _numComponents(componentVector.size())
201 {
202  //create pointer to vector of component weights
203  _componentWeights = new vector<Probability>(NumComponentsGet());
204  for (int i=0; i<NumComponentsGet();i++)
205  {
206  (*_componentWeights)[i] = (Probability)(1.0/NumComponentsGet());
207  }
208  _cumWeights.insert(_cumWeights.begin(),NumComponentsGet()+1,0.0);
210 
211  //create pointer to vector of pointers to weights
212  _componentPdfs = new vector< Pdf<T>* >(NumComponentsGet());
213  //TODO: take copy or point to same???
214  for (int i=0; i<NumComponentsGet();i++)
215  {
216  //TODO: will this call the constructor of e.g. Gaussian or always the
217  //general one?
218  (*_componentPdfs)[i] = (componentVector[i])->Clone();
219  }
220 #ifdef __CONSTRUCTOR__
221  cout << "Mixture constructor\n";
222 #endif // __CONSTRUCTOR__
223 }
224 
225 template<typename T >
226 Mixture<T>::Mixture(const Mixture & my_mixture):Pdf<T>(my_mixture.DimensionGet() )
227  ,_numComponents(my_mixture.NumComponentsGet())
228  {
229  //create pointer to vector of component weights
230  _componentWeights = new vector<Probability>(this->NumComponentsGet());
231  (*_componentWeights) = my_mixture.WeightsGet();
232  _cumWeights.insert(_cumWeights.begin(),NumComponentsGet()+1,0.0);
234 
235  //create pointer to vector of pointers to weights
236  _componentPdfs = new vector< Pdf<T>* >(NumComponentsGet());
237  for (int i=0; i<NumComponentsGet();i++)
238  {
239  (*_componentPdfs)[i] = (my_mixture.ComponentGet(i))->Clone();
240  }
241 #ifdef __CONSTRUCTOR__
242  cout << "Mixture copy constructor\n";
243 #endif // __CONSTRUCTOR__
244  }
245 
246 template<typename T>
248  {
249 #ifdef __CONSTRUCTOR__
250  cout << "Mixture destructor\n";
251 #endif
252  // Release memory!
253  delete _componentWeights;
254  for (int i=0; i<NumComponentsGet();i++)
255  {
256  delete (*_componentPdfs)[i];
257  }
258  delete _componentPdfs;
259  }
260 
261 template<typename T>
263  {
264  return new Mixture(*this);
265  }
266 
267 template<typename T>
268  unsigned int Mixture<T>::NumComponentsGet() const
269  {
270  return _numComponents;
271  }
272 
273 template<typename T>
274  Probability Mixture<T>::ProbabilityGet(const T& state) const
275  {
276  TestNotInit();
277  Probability prob(0.0);
278  for (int i=0; i<NumComponentsGet();i++)
279  {
280  prob= prob + (*_componentWeights)[i] * (*_componentPdfs)[i]->ProbabilityGet(state);
281  }
282  return prob;
283  }
284 
285 template<typename T>
286  bool Mixture<T>::SampleFrom (vector<Sample<T> >& list_samples,
287  const unsigned int num_samples,
288  int method,
289  void * args) const
290  {
291  TestNotInit();
292  switch(method)
293  {
294  case DEFAULT: // O(N log(N) efficiency)
295  return Pdf<T>::SampleFrom(list_samples, num_samples,method,args);
296 
297  case RIPLEY: // See mcpdf.cpp for more explanation
298  {
299  list_samples.resize(num_samples);
300  // GENERATE N ORDERED IID UNIFORM SAMPLES
301  std::vector<double> unif_samples(num_samples);
302  for ( unsigned int i = 0; i < num_samples ; i++)
303  unif_samples[i] = runif();
304 
305  /* take n-th racine of u_N */
306  unif_samples[num_samples-1] = pow(unif_samples[num_samples-1],
307  double (1.0/num_samples));
308  /* rescale samples */
309  for ( int i = num_samples-2; i >= 0 ; i--)
310  unif_samples[i] = pow(unif_samples[i], double (1.0/(i+1))) * unif_samples[i+1];
311 
312  // CHECK WHERE THESE SAMPLES ARE IN _CUMWEIGHTS
313  unsigned int index = 0;
314  unsigned int num_states = NumComponentsGet();
315  vector<double>::const_iterator CumPDFit = _cumWeights.begin();
316  typename vector<Sample<T> >::iterator sit = list_samples.begin();
317 
318  for ( unsigned int i = 0; i < num_samples ; i++)
319  {
320  while ( unif_samples[i] > *CumPDFit )
321  {
322  // check for internal error
323  assert(index <= num_states);
324  index++; CumPDFit++;
325  }
326  // index-1 is a sample of the discrete pdf of the mixture weights
327  // get a sample from the index-1'th mixture component
328  (*_componentPdfs)[index-1]->SampleFrom(*sit,method,args);
329  sit++;
330  }
331  return true;
332  }
333  default:
334  cerr << "Mixture::Samplefrom(T, void *): No such sampling method" << endl;
335  return false;
336  }
337  }
338 template<typename T>
339  bool Mixture<T>::SampleFrom (Sample<T>& one_sample, int method, void * args) const
340  {
341  TestNotInit();
342  switch(method)
343  {
344  case DEFAULT:
345  {
346  // Sample from univariate uniform rng between 0 and 1;
347  double unif_sample; unif_sample = runif();
348  // Compare where we should be
349  unsigned int index = 0;
350  while ( unif_sample > _cumWeights[index] )
351  {
352  assert(index <= NumComponentsGet());
353  index++;
354  }
355  // index-1 is a sample of the discrete pdf of the mixture weights
356  // get a sample from the index-1'th mixture component
357  (*_componentPdfs)[index-1]->SampleFrom(one_sample,method,args);
358  return true;
359  }
360  default:
361  cerr << "Mixture::Samplefrom(T, void *): No such sampling method"
362  << endl;
363  return false;
364  }
365  }
366 
367 template<typename T>
369  {
370  cerr << "Mixture ExpectedValueGet: not implemented for the template parameters you use."
371  << endl << "Use template specialization as shown in mixture.cpp " << endl;
372  assert(0);
373  T expectedValue;
374  return expectedValue;
375  }
376 
377 template <typename T>
378  MatrixWrapper::SymmetricMatrix Mixture<T>::CovarianceGet ( ) const
379  {
380  TestNotInit();
381  cerr << "Mixture CovarianceGet: not implemented since so far I don't believe its usefull"
382  << endl << "If you decide to implement is: Use template specialization as shown in mcpdf.cpp " << endl;
383 
384  assert(0);
385  MatrixWrapper::SymmetricMatrix result;
386  return result;
387  }
388 
389 template<typename T>
390  vector<Probability> Mixture<T>::WeightsGet() const
391  {
392  TestNotInit();
393  return *_componentWeights;
394  }
395 
396 template<typename T>
397  Probability Mixture<T>::WeightGet(unsigned int componentNumber) const
398  {
399  TestNotInit();
400  assert((int)componentNumber >= 0 && componentNumber < NumComponentsGet());
401  return (*_componentWeights)[componentNumber];
402  }
403 
404 template<typename T>
405  bool Mixture<T>::WeightsSet(vector<Probability> & weights)
406  {
407  TestNotInit();
408  assert(weights.size() == NumComponentsGet());
409  *_componentWeights = weights;
410  //normalize the probabilities and update the cumulative sum
411  return (NormalizeWeights() && CumWeightsUpdate());
412  }
413 
414 template<typename T>
415  bool Mixture<T>::WeightSet(unsigned int componentNumber, Probability weight)
416  {
417  TestNotInit();
418  assert((int)componentNumber >= 0 && componentNumber < NumComponentsGet());
419  assert((double)weight<=1.0);
420 
421  if (NumComponentsGet() == 1)
422  {
423  (*_componentWeights)[0] = weight;
424  }
425  else
426  {
427  // renormalize other weights such that sum of probabilities will be
428  // one after the weight of the component is set to weight
429  // This should keep the probabilities normalized
430  Probability old_weight = WeightGet(componentNumber);
431  if ((double)old_weight!=1.0) {
432  double normalization_factor = (1-weight)/(1-old_weight);
433  for (int i=0; i<NumComponentsGet();i++)
434  {
435  (*_componentWeights)[i] = (Probability)( (double)( (*_componentWeights)[i] )* normalization_factor);
436  }
437  }
438  else{
439  for (int i=0; i<NumComponentsGet();i++)
440  {
441  (*_componentWeights)[i] = (Probability)( (1-weight)/(NumComponentsGet()-1));
442  }
443  }
444  (*_componentWeights)[componentNumber] = weight;
445  }
446  return CumWeightsUpdate();
447  }
448 
449 template<typename T>
451  {
452  TestNotInit();
453  int index_mostProbable= -1;
454  Probability prob_mostProbable= 0.0;
455  for (int component = 0 ; component < NumComponentsGet() ; component++)
456  {
457  if ( (*_componentWeights)[component] > prob_mostProbable)
458  {
459  index_mostProbable= component;
460  prob_mostProbable= (*_componentWeights)[component];
461  }
462  }
463  return index_mostProbable;
464  }
465 
466 template<typename T>
467  bool Mixture<T>::AddComponent(Pdf<T>& pdf)
468  {
469  if (NumComponentsGet()==0)
470  return AddComponent(pdf, Probability(1.0));
471  else
472  {
473  _numComponents++;
474  (*_componentPdfs).push_back(pdf.Clone() );
475 
476  (*_componentWeights).push_back(Probability(0.0));
477  _cumWeights.push_back(0.0);
478  //assert length of vectors
479  assert(NumComponentsGet()==(*_componentPdfs).size());
480  assert(NumComponentsGet()==(*_componentWeights).size());
481  assert(NumComponentsGet()+1==_cumWeights.size());
482  return (NormalizeWeights() && CumWeightsUpdate());
483  }
484  }
485 
486 template<typename T>
488  {
489  if (NumComponentsGet()==0 && w!=1.0)
490  return AddComponent(pdf, Probability(1.0));
491  else
492  {
493  _numComponents++;
494  (*_componentPdfs).push_back(pdf.Clone() );
495  (*_componentWeights).push_back(Probability(0.0));
496  _cumWeights.resize(NumComponentsGet()+1);
497  //assert length of vectors
498  assert(NumComponentsGet()==(*_componentPdfs).size());
499  assert(NumComponentsGet()==(*_componentWeights).size());
500  assert(NumComponentsGet()+1==_cumWeights.size());
502  return (NormalizeWeights() && CumWeightsUpdate());
503  }
504  }
505 
506 template<typename T>
507  bool Mixture<T>::DeleteComponent(unsigned int componentNumber)
508  {
509  //assert length of vectors
510  assert(NumComponentsGet()==(*_componentPdfs).size());
511  assert(NumComponentsGet()==(*_componentWeights).size());
512  assert(NumComponentsGet()+1==_cumWeights.size());
513 
514  TestNotInit();
515  assert((int)componentNumber >= 0 && componentNumber < NumComponentsGet());
516  _numComponents--;
517  Pdf<T>* pointer = (*_componentPdfs)[componentNumber];
518  delete pointer;
519  (*_componentPdfs).erase((*_componentPdfs).begin()+componentNumber);
520  (*_componentWeights).erase((*_componentWeights).begin()+componentNumber);
521  _cumWeights.resize(NumComponentsGet()+1);
522  //assert length of vectors
523  assert(NumComponentsGet()==(*_componentPdfs).size());
524  assert(NumComponentsGet()==(*_componentWeights).size());
525  assert(NumComponentsGet()+1==_cumWeights.size());
526  if(_numComponents==0) //don't do normalization if numComponents == 0
527  return true;
528  else
529  return (NormalizeWeights() && CumWeightsUpdate());
530  }
531 
532 template<typename T>
533  vector<Pdf<T>*> Mixture<T>::ComponentsGet() const
534  {
535  TestNotInit();
536  return _componentPdfs;
537  }
538 
539 template<typename T>
540  Pdf<T>* Mixture<T>::ComponentGet(unsigned int componentNumber) const
541  {
542  TestNotInit();
543  return (*_componentPdfs)[componentNumber];
544  }
545 
546 template<typename T>
547  void Mixture<T>::TestNotInit() const
548  {
549  if (NumComponentsGet() == 0)
550  {
551  cerr << "Mixture method called which requires that the number of components is not zero"
552  << endl << "Current number of components: " << NumComponentsGet() << endl;
553  assert(0);
554  }
555  }
556 
557 template<typename T>
559  {
560  double SumOfWeights = 0.0;
561  for ( unsigned int i = 0; i < NumComponentsGet() ; i++){
562  SumOfWeights += (*_componentWeights)[i];
563  }
564  if (SumOfWeights > 0){
565  for ( unsigned int i = 0; i < NumComponentsGet() ; i++){
566  (*_componentWeights)[i] = (Probability)( (double) ( (*_componentWeights)[i]) /SumOfWeights);
567  }
568  return true;
569  }
570  else{
571  cerr << "Mixture::NormalizeProbs(): SumOfWeights = " << SumOfWeights << endl;
572  return false;
573  }
574  }
575 
576 template<typename T>
578  {
579  // precondition: sum of probabilities should be 1
580  double CumSum=0.0;
581  static vector<double>::iterator CumWeightsit;
582  CumWeightsit = _cumWeights.begin();
583  *CumWeightsit = 0.0;
584 
585  // Calculate the Cumulative PDF
586  for ( unsigned int i = 0; i < NumComponentsGet(); i++)
587  {
588  CumWeightsit++;
589  // Calculate the __normalised__ Cumulative sum!!!
590  CumSum += ( (*_componentWeights)[i] );
591  *CumWeightsit = CumSum;
592  }
593  // Check if last element of valuelist is +- 1
594  assert( (_cumWeights[NumComponentsGet()] >= 1.0 - NUMERIC_PRECISION) &&
596 
598  return true;
599  }
600 
601 } // End namespace
602 
603 #include "mixture.cpp"
604 
605 #endif // MIXTURE_H
bool WeightSet(unsigned int componentNumber, Probability w)
Function to change/set the weigth of a single component.
int MostProbableComponentGet() const
Get the index of the most probable component, if a few component are.
Class PDF: Virtual Base class representing Probability Density Functions.
Definition: pdf.h:53
#define DEFAULT
Probability WeightGet(unsigned int componentNumber) const
Get the component weight of component "componentNumber".
Mixture(const unsigned int dimension=0)
Constructor: An equal weight is set for all components.
vector< Probability > WeightsGet() const
Get all component weights.
MatrixWrapper::SymmetricMatrix CovarianceGet() const
Get the Covariance Matrix E[(x - E[x])^2] of the Analytic pdf.
vector< Pdf< T > *> * _componentPdfs
Pointer to the vector of component pdfs.
Definition: mixture.h:58
bool WeightsSet(vector< Probability > &weights)
Set all component weights.
bool AddComponent(Pdf< T > &pdf)
Add a component pdf: THIS IS A NON-REALTIME OPERATION.
void TestNotInit() const
Called when a the number of components=0 and if method is called which.
unsigned int _numComponents
The number of components.
Definition: mixture.h:52
vector< double > _cumWeights
Vector containing the cumulative component weights (for efficient sampling)
Definition: mixture.h:66
Class representing a mixture of PDFs, the mixture can contain different.
Definition: mixture.h:48
bool CumWeightsUpdate()
Updates the cumWeights.
bool SampleFrom(vector< Sample< T > > &list_samples, const unsigned int num_samples, int method=DEFAULT, void *args=NULL) const
unsigned int NumComponentsGet() const
Get the number of components.
virtual Pdf< T > * Clone() const =0
Pure virtual clone function.
Probability ProbabilityGet(const T &state) const
Implementation of virtual base class method.
unsigned int DimensionGet() const
Get the dimension of the argument.
#define RIPLEY
virtual bool SampleFrom(vector< Sample< T > > &list_samples, const unsigned int num_samples, int method=DEFAULT, void *args=NULL) const
Draw multiple samples from the Pdf (overloaded)
vector< Probability > * _componentWeights
Pointer to the vector of mixture weights, the sum of the elements = 1.
Definition: mixture.h:55
virtual Mixture * Clone() const
Clone function.
double runif()
Class representing a probability (a double between 0 and 1)
Definition: bfl_constants.h:39
Pdf< T > * ComponentGet(unsigned int componentNumber) const
Get the pointer to the component pdf of component "componentNumber".
#define NUMERIC_PRECISION
bool NormalizeWeights()
Normalize the component weigths (eg. after setting a component weight)
bool DeleteComponent(unsigned int componentNumber)
Delete a component pdf: THIS IS A NON_REALTIME OPERATION.
virtual ~Mixture()
Destructor.
T ExpectedValueGet() const
Get the expected value E[x] of the pdf.
vector< Pdf< T > *> ComponentsGet() const
Get the vector of pointers to the component pdfs.


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 Mon Feb 28 2022 21:56:33