Bayesian Filtering Library Generated from SVN r
mixture.h
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., 51 Franklin Street, Fifth Floor, *
27 * Boston, MA 02110-1301 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
40namespace 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 const SampleMthd method = SampleMthd::DEFAULT,
105 void * args = NULL) const;
106
107 bool SampleFrom (Sample<T>& one_sample, const SampleMthd method = SampleMthd::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!
181template<typename T>
182Mixture<T>::Mixture(const unsigned int dimension):
183 Pdf<T>(dimension)
184 , _numComponents(0)
185 , _cumWeights(_numComponents+1)
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
198template<typename T> template <typename U>
199Mixture<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);
209 CumWeightsUpdate();
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
225template<typename T >
226Mixture<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);
233 CumWeightsUpdate();
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
246template<typename T>
247 Mixture<T>::~Mixture()
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
261template<typename T>
262 Mixture<T>* Mixture<T>::Clone() const
263 {
264 return new Mixture(*this);
265 }
266
267template<typename T>
268 unsigned int Mixture<T>::NumComponentsGet() const
269 {
270 return _numComponents;
271 }
272
273template<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
285template<typename T>
286 bool Mixture<T>::SampleFrom (vector<Sample<T> >& list_samples,
287 const unsigned int num_samples,
288 const SampleMthd method,
289 void * args) const
290 {
291 TestNotInit();
292 switch(method)
293 {
294 case SampleMthd::DEFAULT: // O(N log(N) efficiency)
295 return Pdf<T>::SampleFrom(list_samples, num_samples,method,args);
296
297 case SampleMthd::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 }
338template<typename T>
339 bool Mixture<T>::SampleFrom (Sample<T>& one_sample, const SampleMthd method, void * args) const
340 {
341 TestNotInit();
342 switch(method)
343 {
344 case SampleMthd::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
367template<typename T>
368 T Mixture<T>::ExpectedValueGet() const
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
377template <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
389template<typename T>
390 vector<Probability> Mixture<T>::WeightsGet() const
391 {
392 TestNotInit();
393 return *_componentWeights;
394 }
395
396template<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
404template<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
414template<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
449template<typename T>
450 int Mixture<T>::MostProbableComponentGet() const
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
466template<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
486template<typename T>
487 bool Mixture<T>::AddComponent(Pdf<T>& pdf, Probability w)
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());
501 WeightSet(_numComponents-1,w);
502 return (NormalizeWeights() && CumWeightsUpdate());
503 }
504 }
505
506template<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
532template<typename T>
533 vector<Pdf<T>*> Mixture<T>::ComponentsGet() const
534 {
535 TestNotInit();
536 return _componentPdfs;
537 }
538
539template<typename T>
540 Pdf<T>* Mixture<T>::ComponentGet(unsigned int componentNumber) const
541 {
542 TestNotInit();
543 return (*_componentPdfs)[componentNumber];
544 }
545
546template<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
557template<typename T>
558 bool Mixture<T>::NormalizeWeights()
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
576template<typename T>
577 bool Mixture<T>::CumWeightsUpdate()
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) &&
595 (_cumWeights[NumComponentsGet()] <= 1.0 + NUMERIC_PRECISION) );
596
597 _cumWeights[NumComponentsGet()]=1;
598 return true;
599 }
600
601} // End namespace
602
603#include "mixture.cpp"
604
605#endif // MIXTURE_H