[73] | 1 | // This file is part of the PACPUS framework distributed under the
|
---|
| 2 | // CECILL-C License, Version 1.0.
|
---|
| 3 | //
|
---|
| 4 | /// @file
|
---|
| 5 | /// @author Firstname Surname <firstname.surname@utc.fr>
|
---|
| 6 | /// @date Month, Year
|
---|
| 7 | /// @version $Id: particle_filtering.hpp 73 2013-01-10 16:56:42Z kurdejma $
|
---|
| 8 | /// @copyright Copyright (c) UTC/CNRS Heudiasyc 2006 - 2013. All rights reserved.
|
---|
| 9 | /// @brief Brief description.
|
---|
| 10 | ///
|
---|
| 11 | /// Detailed description.
|
---|
| 12 |
|
---|
| 13 | #ifndef __PARTICLE_FILTERING_BASE__
|
---|
| 14 | #define __PARTICLE_FILTERING_BASE__
|
---|
| 15 |
|
---|
| 16 | #include "bayes_filtering.hpp"
|
---|
| 17 | #include "../math/ublas.hpp"
|
---|
| 18 | #include "../math/rng.hpp"
|
---|
| 19 | #include "../math/distributions.hpp"
|
---|
| 20 |
|
---|
| 21 | namespace filter{
|
---|
| 22 | namespace particle {
|
---|
| 23 |
|
---|
| 24 | using namespace math;
|
---|
| 25 | using namespace ublas;
|
---|
| 26 | using namespace rng;
|
---|
| 27 | using namespace distributions;
|
---|
| 28 |
|
---|
| 29 |
|
---|
| 30 | /*!
|
---|
| 31 | * \class Particle
|
---|
| 32 | * \brief This class describe a particle \n
|
---|
| 33 | * A particle is constituted by : 1) a state vector and 2) importance weight \n
|
---|
| 34 | */
|
---|
| 35 | class Particle {
|
---|
| 36 | public :
|
---|
| 37 | /*!
|
---|
| 38 | * \brief Constructor
|
---|
| 39 | * \param state_size : the size of the state vector
|
---|
| 40 | * \param weight_ : the initial weight of the particle
|
---|
| 41 | */
|
---|
| 42 | Particle(const size_t state_size, const double & weight_){
|
---|
| 43 | X=ZeroVector(state_size);
|
---|
| 44 | weight=weight_;
|
---|
| 45 | }
|
---|
| 46 |
|
---|
| 47 | Vector X; /*!< state vector */
|
---|
| 48 |
|
---|
| 49 | double weight; /*!< particle weight */
|
---|
| 50 |
|
---|
| 51 | };
|
---|
| 52 |
|
---|
| 53 | /*! \class ParticleSet
|
---|
| 54 | * \brief This class describe a set of particles \n
|
---|
| 55 | * A set of particle is reprented by a vector of particles \n
|
---|
| 56 | * somes methods can be applied to the set of particles like : \n
|
---|
| 57 | * estimate computation, resampling scheme or normalization method \n
|
---|
| 58 | */
|
---|
| 59 |
|
---|
| 60 | template<class P=Particle> class ParticleSet {
|
---|
| 61 | public :
|
---|
| 62 |
|
---|
| 63 | std::vector<P> particles; /*!< vector of particles */
|
---|
| 64 |
|
---|
| 65 |
|
---|
| 66 | std::vector<P> particles_; /*!< temporary vector of particles used in resample scheme */
|
---|
| 67 |
|
---|
| 68 |
|
---|
| 69 | /*!
|
---|
| 70 | * \brief Allocate the set of paticles
|
---|
| 71 | * \param nparticle : number of particles
|
---|
| 72 | */
|
---|
| 73 | void Allocator(const size_t &nparticle);
|
---|
| 74 |
|
---|
| 75 | /*!
|
---|
| 76 | * \brief Resample the set of particles \n
|
---|
| 77 | * \param threshold : threshold parameter \n
|
---|
| 78 | * The resampling scheme is the reasmpling scheme proposed by Kitawaga \n
|
---|
| 79 | * This can be overloaded to use an another resampling scheme \n
|
---|
| 80 | */
|
---|
| 81 | virtual void Resample(const double & threshold=1.0);
|
---|
| 82 |
|
---|
| 83 | /*!
|
---|
| 84 | * \brief Compute the number of effective samples in the particle set
|
---|
| 85 | * \return number of effective samples
|
---|
| 86 | */
|
---|
| 87 | double Ness();
|
---|
| 88 |
|
---|
| 89 | /*!
|
---|
| 90 | * \brief Get the number of particles
|
---|
| 91 | * \return number of particles
|
---|
| 92 | */
|
---|
| 93 | size_t ParticlesNo(){ return particles.size(); }
|
---|
| 94 |
|
---|
| 95 | /*!
|
---|
| 96 | * \brief Compute the estimate
|
---|
| 97 | * \return a vector double describing the estimate
|
---|
| 98 | */
|
---|
| 99 | Vector Estimate ();
|
---|
| 100 |
|
---|
| 101 | /*!
|
---|
| 102 | * \brief Normalize the weights of each particles
|
---|
| 103 | */
|
---|
| 104 | void NormalizeWeights();
|
---|
| 105 |
|
---|
| 106 | /*!
|
---|
| 107 | * \brief Destructor
|
---|
| 108 | */
|
---|
| 109 | virtual ~ParticleSet(){};
|
---|
| 110 |
|
---|
| 111 | protected :
|
---|
| 112 |
|
---|
| 113 | /*!
|
---|
| 114 | * \brief Compute the cumulative sum of the weight of each particles
|
---|
| 115 | * \return a vector of double describing this cumulative sum
|
---|
| 116 | */
|
---|
| 117 | std::vector<double> WeightCumSum();
|
---|
| 118 |
|
---|
| 119 | /*!
|
---|
| 120 | * \brief Create the a random cumulative sum
|
---|
| 121 | * \return a vector of double describing this cumulative sum
|
---|
| 122 | */
|
---|
| 123 | std::vector<double> RandomCumSum();
|
---|
| 124 |
|
---|
| 125 | };
|
---|
| 126 |
|
---|
| 127 |
|
---|
| 128 | //Particle set memeber functions
|
---|
| 129 | template<class P> void ParticleSet<P>::Allocator(const size_t &nparticles){
|
---|
| 130 | particles.resize(nparticles);
|
---|
| 131 | particles_.resize(nparticles);
|
---|
| 132 | }
|
---|
| 133 |
|
---|
| 134 | template<class P> double ParticleSet<P>::Ness(){
|
---|
| 135 | double ness =0;
|
---|
| 136 | for(typename std::vector<P>::iterator I=particles.begin();I!=particles.end();I++)
|
---|
| 137 | ness+=(*I).weight*(*I).weight;
|
---|
| 138 |
|
---|
| 139 | if(ness==0)throw filter_error("In particle filter :: number effective sample computation -> ness is equal to infinity");
|
---|
| 140 |
|
---|
| 141 | return 1/ness;
|
---|
| 142 | }
|
---|
| 143 |
|
---|
| 144 |
|
---|
| 145 | template<class P> std::vector<double> ParticleSet<P>::WeightCumSum(){
|
---|
| 146 | double sum=0;
|
---|
| 147 | std::vector<double> cumsum(particles.size());
|
---|
| 148 | for(size_t i=0;i<particles.size();i++){
|
---|
| 149 | sum+=particles[i].weight;
|
---|
| 150 | cumsum[i]=sum;
|
---|
| 151 | }
|
---|
| 152 | return cumsum;
|
---|
| 153 | }
|
---|
| 154 |
|
---|
| 155 | template<class P> void ParticleSet<P>::NormalizeWeights(){
|
---|
| 156 |
|
---|
| 157 | double sum=0;
|
---|
| 158 | for(typename std::vector<P>::iterator I=particles.begin();I!=particles.end();I++)
|
---|
| 159 | sum+=(*I).weight;
|
---|
| 160 |
|
---|
| 161 |
|
---|
| 162 | if(sum==0) throw filter_error("In particle filter :: normalization of weight -> weight sum is equal to 0");
|
---|
| 163 |
|
---|
| 164 | for(typename std::vector<P>::iterator I=particles.begin();I!=particles.end();I++)
|
---|
| 165 | (*I).weight=(*I).weight/sum;
|
---|
| 166 |
|
---|
| 167 | }
|
---|
| 168 |
|
---|
| 169 | template<class P> std::vector<double> ParticleSet<P>::RandomCumSum(){
|
---|
| 170 | double sum=0;
|
---|
| 171 | std::vector<double> cumsum(particles.size());
|
---|
| 172 | uniform_rng UD;
|
---|
| 173 | for(size_t i=0;i<particles.size();i++){
|
---|
| 174 | sum+=UD();
|
---|
| 175 | cumsum[i]=sum;
|
---|
| 176 | }
|
---|
| 177 |
|
---|
| 178 | //normalization
|
---|
| 179 | for(size_t i=0;i<particles.size();i++)cumsum[i]=cumsum[i]/sum;
|
---|
| 180 |
|
---|
| 181 | return cumsum;
|
---|
| 182 | }
|
---|
| 183 |
|
---|
| 184 |
|
---|
| 185 | template<class P> void ParticleSet<P>::Resample(const double & threshold){
|
---|
| 186 | size_t N=ParticlesNo();
|
---|
| 187 | if(Ness()<=N*threshold){
|
---|
| 188 | // resample method by Carpenter
|
---|
| 189 | std::vector<double> Q=WeightCumSum();
|
---|
| 190 | std::vector<double> T=RandomCumSum();
|
---|
| 191 |
|
---|
| 192 | // select the best particle
|
---|
| 193 | std::vector <size_t> ind(N);
|
---|
| 194 | unsigned int j=0;
|
---|
| 195 | for(unsigned int i=0; i <N; i++){
|
---|
| 196 | do{
|
---|
| 197 | if(T[i] < Q[j]){ind[i] = j;break;}
|
---|
| 198 | j++;
|
---|
| 199 | }while(j<N);
|
---|
| 200 | }
|
---|
| 201 |
|
---|
| 202 | // build the new particle set
|
---|
| 203 | for(unsigned int i=0; i <N; i++){
|
---|
| 204 | particles_[i]=particles[ind[i]];
|
---|
| 205 | particles_[i].weight=1.0/static_cast<double>(N);
|
---|
| 206 | }
|
---|
| 207 |
|
---|
| 208 | // swap the particle set
|
---|
| 209 | particles.swap(particles_);
|
---|
| 210 | }
|
---|
| 211 | }
|
---|
| 212 |
|
---|
| 213 |
|
---|
| 214 |
|
---|
| 215 | template<class P> Vector ParticleSet<P>::Estimate(){
|
---|
| 216 | Vector estimate = ZeroVector(particles[0]->X.size());
|
---|
| 217 | for(typename std::vector<P>::iterator I=particles.begin();I!=particles.end();I++)
|
---|
| 218 | estimate+=(*I).X*(*I).weight;
|
---|
| 219 | return estimate;
|
---|
| 220 | }
|
---|
| 221 |
|
---|
| 222 |
|
---|
| 223 | /*!
|
---|
| 224 | * \class DynamicEquation
|
---|
| 225 | * \brief This class describe a basic dynamic equation for particle filtering \n
|
---|
| 226 | * where : class S is a set particle \n
|
---|
| 227 | * : class I is a particle \n
|
---|
| 228 | * : class D is a random generator used to draw input data \n
|
---|
| 229 | */
|
---|
| 230 | template <class D, template <class> class S=ParticleSet, class P=Particle> class DynamicEquation:public BayesDynamicEquation< S<P>, P >{
|
---|
| 231 | public :
|
---|
| 232 |
|
---|
| 233 | /*!
|
---|
| 234 | * \brief virtual method where matrices of state system must be allocated
|
---|
| 235 | * \param state_size : the size of the state vector of each particle
|
---|
| 236 | * \param data_size : the size of the input vector
|
---|
| 237 | */
|
---|
| 238 | virtual void Allocator(const size_t &state_size,const size_t &data_size)=0;
|
---|
| 239 |
|
---|
| 240 | /*!
|
---|
| 241 | * \brief virtual method where parameters of the dynamic equation must be evaluated
|
---|
| 242 | * \param s : set of particles at time k-1
|
---|
| 243 | *
|
---|
| 244 | */
|
---|
| 245 | virtual void EvaluateParameters( P *s)=0;
|
---|
| 246 |
|
---|
| 247 | /*!
|
---|
| 248 | * \brief virtual method where the a priori state vector must be computed
|
---|
| 249 | * \param in : the set of particles at time k-1
|
---|
| 250 | * \param out : the set of particles at time k
|
---|
| 251 | */
|
---|
| 252 | virtual void Predict(S<P> *in, S<P> *out)=0;
|
---|
| 253 |
|
---|
| 254 |
|
---|
| 255 | /*!
|
---|
| 256 | * \brief Destructor
|
---|
| 257 | */
|
---|
| 258 | virtual ~DynamicEquation(){}
|
---|
| 259 |
|
---|
| 260 | /*!
|
---|
| 261 | * \brief Get/Set random generator used to draw input random U
|
---|
| 262 | */
|
---|
| 263 | D & URNG(){return _URNG;}
|
---|
| 264 |
|
---|
| 265 | protected :
|
---|
| 266 |
|
---|
| 267 | D _URNG; /*!< Random generator used to drawing input data*/
|
---|
| 268 |
|
---|
| 269 | Vector _Urand; /*!< Vector where random input data are stored */
|
---|
| 270 |
|
---|
| 271 | };
|
---|
| 272 |
|
---|
| 273 |
|
---|
| 274 |
|
---|
| 275 |
|
---|
| 276 | /*!
|
---|
| 277 | * \class LinearDynamicEquation
|
---|
| 278 | * \brief This class describe a linear dynamic equation \n
|
---|
| 279 | * X(k+1)(i) =A*X(k)(i)+B*U(k)(i) \n
|
---|
| 280 | * X(.)(i) = the state of the particle i \n
|
---|
| 281 | * A = state matrix \n
|
---|
| 282 | * B = entrie matrix \n
|
---|
| 283 | * U = input random vector \n
|
---|
| 284 | */
|
---|
| 285 | template <class D,template <class> class S=ParticleSet, class P=Particle> class LinearDynamicEquation : public DynamicEquation<D,S,P>{
|
---|
| 286 | public :
|
---|
| 287 | /*!
|
---|
| 288 | * \brief virtual method where matrices of state system must be allocated
|
---|
| 289 | * \param state_size : the size of the state vector of each particle
|
---|
| 290 | * \param data_size : the size of the input vector
|
---|
| 291 | */
|
---|
| 292 | void Allocator(const size_t &state_size,const size_t &data_size);
|
---|
| 293 |
|
---|
| 294 | /*!
|
---|
| 295 | * \brief virtual method where parameters of the dynamic equation must be evaluated
|
---|
| 296 | * \param s : set of particles at time k-1
|
---|
| 297 | */
|
---|
| 298 | virtual void EvaluateParameters(P *s)=0;
|
---|
| 299 |
|
---|
| 300 | /*!
|
---|
| 301 | * \brief virtual method where the a priori state vector must be computed
|
---|
| 302 | * \param in : the set of particles at time k-1
|
---|
| 303 | * \param out : the set of particles at time k
|
---|
| 304 | */
|
---|
| 305 | virtual void Predict(S<P> *in,S<P> *out);
|
---|
| 306 |
|
---|
| 307 | /*!
|
---|
| 308 | * \brief Destructor
|
---|
| 309 | */
|
---|
| 310 | virtual ~LinearDynamicEquation(){}
|
---|
| 311 |
|
---|
| 312 | /*!
|
---|
| 313 | * \brief Get/Set a constant data in A matrix
|
---|
| 314 | */
|
---|
| 315 | double & A(const int &row, const int &column){return _A(row,column);}
|
---|
| 316 |
|
---|
| 317 | /*!
|
---|
| 318 | * \brief Get/Set a constant data in B matrix
|
---|
| 319 | */
|
---|
| 320 | double & B(const int &row, const int &column){return _B(row,column);}
|
---|
| 321 |
|
---|
| 322 | protected :
|
---|
| 323 |
|
---|
| 324 | Matrix _A; /*!< A matrix */
|
---|
| 325 |
|
---|
| 326 | Matrix _B; /*!< A matrix */
|
---|
| 327 |
|
---|
| 328 | };
|
---|
| 329 |
|
---|
| 330 |
|
---|
| 331 | // Particle linear dynamic equation member functions
|
---|
| 332 | template <class D, template <class> class S,class P> void LinearDynamicEquation<D,S,P>::Allocator(const size_t &state_size,const size_t &data_size){
|
---|
| 333 | DynamicEquation<D,S,P>::_Urand=ZeroVector(data_size);
|
---|
| 334 | _A=ZeroMatrix(state_size,state_size);
|
---|
| 335 | _B=ZeroMatrix(state_size,data_size);
|
---|
| 336 | }
|
---|
| 337 |
|
---|
| 338 |
|
---|
| 339 | template <class D, template <class> class S,class P> void LinearDynamicEquation<D,S,P>::Predict(S<P> *in,S<P> *out){
|
---|
| 340 |
|
---|
| 341 | for(size_t i=0;i<in->particles.size();i++){
|
---|
| 342 | EvaluateParameters(&in->particles[i]);
|
---|
| 343 |
|
---|
| 344 | out->particles[i]->X=_A*in->particles[i]->X+_B*DynamicEquation<D,S,P>::_Urand;
|
---|
| 345 | out->particles[i]->weight=in->particles[i]->weight;
|
---|
| 346 | }
|
---|
| 347 | }
|
---|
| 348 |
|
---|
| 349 |
|
---|
| 350 | /*!
|
---|
| 351 | * \class NonLinearDynamicEquation
|
---|
| 352 | * \brief This class describe a linear dynamic equation
|
---|
| 353 | * \n
|
---|
| 354 | * X(k+1)(i) =f(X(k)(i)+U(k)(i)) \n
|
---|
| 355 | * X(.)(i) = the state of the particle i \n
|
---|
| 356 | * A = state matrix \n
|
---|
| 357 | * B = entrie matrix \n
|
---|
| 358 | * U = input random vector \n
|
---|
| 359 | */
|
---|
| 360 | template <class D, template <class> class S,class P> class NonLinearDynamicEquation : public DynamicEquation<D,S,P>{
|
---|
| 361 | public :
|
---|
| 362 | /*!
|
---|
| 363 | * \brief virtual method where matrices of state system must be allocated
|
---|
| 364 | * \param state_size : the size of the state vector of each particle
|
---|
| 365 | * \param data_size : the size of the input vector
|
---|
| 366 | */
|
---|
| 367 | void Allocator(const size_t &state_size,const size_t &data_size);
|
---|
| 368 |
|
---|
| 369 | /*!
|
---|
| 370 | * \brief virtual method where parameters of the dynamic equation must be evaluated
|
---|
| 371 | * \param s : set of particles at time k-1
|
---|
| 372 | * f=
|
---|
| 373 | */
|
---|
| 374 | virtual void EvaluateParameters(P *s)=0;
|
---|
| 375 |
|
---|
| 376 | /*!
|
---|
| 377 | * \brief virtual method where the a priori state vector must be computed
|
---|
| 378 | * \param in : the set of particles at time k-1
|
---|
| 379 | * \param out : the set of particles at time k
|
---|
| 380 | */
|
---|
| 381 | virtual void Predict(S<P> *in,S<P> *out);
|
---|
| 382 |
|
---|
| 383 | /*!
|
---|
| 384 | * \brief Destructor
|
---|
| 385 | */
|
---|
| 386 | virtual ~NonLinearDynamicEquation(){}
|
---|
| 387 |
|
---|
| 388 |
|
---|
| 389 | protected:
|
---|
| 390 |
|
---|
| 391 | Vector _f; /*!< Matrix f where result of (f(X,U)) is stored */
|
---|
| 392 |
|
---|
| 393 | };
|
---|
| 394 |
|
---|
| 395 | // Particle non linear dynamic equation member functions
|
---|
| 396 | template <class D,template <class> class S,class P> void NonLinearDynamicEquation<D,S,P>::Allocator(const size_t &state_size,const size_t &data_size){
|
---|
| 397 | DynamicEquation<D,S,P>::_Urand=ZeroVector(data_size);
|
---|
| 398 | _f=ZeroVector(state_size);
|
---|
| 399 | }
|
---|
| 400 |
|
---|
| 401 |
|
---|
| 402 |
|
---|
| 403 | template <class D,template <class> class S,class P> void NonLinearDynamicEquation<D,S,P>::Predict(S<P> *in,S<P> *out){
|
---|
| 404 | for(size_t i=0;i<in->particles.size();i++){
|
---|
| 405 | EvaluateParameters(&in->particles[i]);
|
---|
| 406 |
|
---|
| 407 | out->particles[i].X=_f;
|
---|
| 408 | out->particles[i].weight=in->particles[i].weight;
|
---|
| 409 |
|
---|
| 410 | }
|
---|
| 411 | }
|
---|
| 412 |
|
---|
| 413 | /*!
|
---|
| 414 | * \class MeasureEquation
|
---|
| 415 | * \brief This clas describe a basic measure equation particle filtering
|
---|
| 416 | * \n
|
---|
| 417 | * where : class S is a set particle \n
|
---|
| 418 | * : class I is a particle \n
|
---|
| 419 | * : class D is the distribution of probability of observation data \n
|
---|
| 420 | */
|
---|
| 421 | template <class D , template <class> class S=ParticleSet, class P=Particle> class MeasureEquation: public BayesMeasureEquation< S<P>, P >{
|
---|
| 422 | public :
|
---|
| 423 |
|
---|
| 424 | /*!
|
---|
| 425 | * \brief virtual method where matrices of state system must be allocated
|
---|
| 426 | * \param state_size : the size of the state vector of each particle
|
---|
| 427 | * \param data_size : the size of the input vector
|
---|
| 428 | */
|
---|
| 429 | virtual void Allocator(const size_t &state_size,const size_t &data_size)=0;
|
---|
| 430 |
|
---|
| 431 |
|
---|
| 432 | /*!
|
---|
| 433 | * \brief virtual method where parameters of the measure equation must be evaluated
|
---|
| 434 | * \param s : set of particles at time k
|
---|
| 435 | */
|
---|
| 436 | virtual void EvaluateParameters(P *s)=0;
|
---|
| 437 |
|
---|
| 438 | /*!
|
---|
| 439 | * \brief virtual method where the a posteriori state vector must be computed
|
---|
| 440 | * \param in : the a priori set of particles at time k
|
---|
| 441 | * \param out : the a posteriori set of particles at time k
|
---|
| 442 | */
|
---|
| 443 | virtual void Update(S<P> *in,S<P> *out)=0;
|
---|
| 444 |
|
---|
| 445 | /*!
|
---|
| 446 | * \brief Destructor
|
---|
| 447 | */
|
---|
| 448 | virtual ~MeasureEquation(){};
|
---|
| 449 |
|
---|
| 450 | /*!
|
---|
| 451 | * \brief Get/Set the distribution of probability of data observation Z
|
---|
| 452 | */
|
---|
| 453 | D & ZDistribution(){return _ZDistribution;}
|
---|
| 454 |
|
---|
| 455 | protected :
|
---|
| 456 |
|
---|
| 457 | D _ZDistribution; /*!< The distribution of probability of observation data Z*/
|
---|
| 458 | };
|
---|
| 459 |
|
---|
| 460 |
|
---|
| 461 |
|
---|
| 462 | /*!
|
---|
| 463 | * \class LinearMeasureEquation
|
---|
| 464 | * \brief This class describe a linear measure equation
|
---|
| 465 | * \n
|
---|
| 466 | * X(k)(i)=H*Z(k) \n
|
---|
| 467 | * Z(k) = observation data \n
|
---|
| 468 | * => weigth update \n
|
---|
| 469 | * w(k+1)(i)=p(Z|HX(i))w(k)(i) \n
|
---|
| 470 | * w(.)(i)= weight of particle i \n
|
---|
| 471 | */
|
---|
| 472 | template <class D, template <class> class S=ParticleSet, class P=Particle> class LinearMeasureEquation : public MeasureEquation< D,S,P >{
|
---|
| 473 | public :
|
---|
| 474 | /*!
|
---|
| 475 | * \brief virtual method where matrices of state system must be allocated
|
---|
| 476 | * \param state_size : the size of the state vector of each particle
|
---|
| 477 | * \param data_size : the size of the input vector
|
---|
| 478 | */
|
---|
| 479 | void Allocator(const size_t &state_size,const size_t &data_size);
|
---|
| 480 |
|
---|
| 481 |
|
---|
| 482 | /*!
|
---|
| 483 | * \brief virtual method where parameters of the measure equation must be evaluated
|
---|
| 484 | * \param s : set of particles at time k
|
---|
| 485 | */
|
---|
| 486 | virtual void EvaluateParameters(P *s)=0;
|
---|
| 487 |
|
---|
| 488 | /*!
|
---|
| 489 | * \brief virtual method where the a posteriori state vector must be computed
|
---|
| 490 | * \param in : the a priori set of particles at time k
|
---|
| 491 | * \param out : the a posteriori set of particles at time k
|
---|
| 492 | */
|
---|
| 493 | virtual void Update(S<P> *in,S<P> *out);
|
---|
| 494 |
|
---|
| 495 | /*!
|
---|
| 496 | * \brief Destructor
|
---|
| 497 | */
|
---|
| 498 | virtual ~LinearMeasureEquation(){}
|
---|
| 499 |
|
---|
| 500 | /*!
|
---|
| 501 | * \brief Get/Set a constant data in observation matrix H
|
---|
| 502 | */
|
---|
| 503 | double & H(int row, int column){return _H(row,column);}
|
---|
| 504 |
|
---|
| 505 | protected:
|
---|
| 506 |
|
---|
| 507 | Matrix _H; /*!< Observation matrix */
|
---|
| 508 | };
|
---|
| 509 |
|
---|
| 510 |
|
---|
| 511 | // Particle linear measure equation member functions
|
---|
| 512 | template <class D,template <class> class S, class P> void LinearMeasureEquation<D,S,P>::Allocator(const size_t &state_size,const size_t &data_size){
|
---|
| 513 | _H=ZeroMatrix(data_size,state_size);
|
---|
| 514 | }
|
---|
| 515 |
|
---|
| 516 | template <class D,template <class> class S, class P> void LinearMeasureEquation<D,S,P>::Update(S<P> *in,S<P> *out){
|
---|
| 517 | for(size_t i=0;i<in->particles.size();i++){
|
---|
| 518 | EvaluateParameters(&in->particles[i]);
|
---|
| 519 |
|
---|
| 520 | out->particles[i].weight=pdf(MeasureEquation<D,S,P>::_ZDistribution, _H*in->particles[i].X)*in->particles[i].weight;
|
---|
| 521 |
|
---|
| 522 | }
|
---|
| 523 | out->NormalizeWeights();
|
---|
| 524 | }
|
---|
| 525 |
|
---|
| 526 |
|
---|
| 527 |
|
---|
| 528 | /*!
|
---|
| 529 | * \class NonLinearMeasureEquation
|
---|
| 530 | * \brief This class describe a non linear measure equation
|
---|
| 531 | * \n
|
---|
| 532 | * Z(k)=h(X)(i) \n
|
---|
| 533 | * Z = observation data \n
|
---|
| 534 | * => weigth update \n
|
---|
| 535 | * w(k+1)(i)=p(Z|h(X(k)(i)))w(k)(i) \n
|
---|
| 536 | * w(.)(i)= weight of particle i \n
|
---|
| 537 | */
|
---|
| 538 | template <class D,template <class> class S=ParticleSet, class P=Particle> class NonLinearMeasureEquation : public MeasureEquation<D,S,P>{
|
---|
| 539 | public :
|
---|
| 540 | /*!
|
---|
| 541 | * \brief virtual method where matrices of state system must be allocated
|
---|
| 542 | * \param state_size : the size of the state vector of each particle
|
---|
| 543 | * \param data_size : the size of the input vector
|
---|
| 544 | */
|
---|
| 545 | void Allocator(const size_t &state_size,const size_t &data_size);
|
---|
| 546 |
|
---|
| 547 | /*!
|
---|
| 548 | * \brief virtual method where parameters of the measure equation must be evaluated
|
---|
| 549 | * \param s : set of particles at time k
|
---|
| 550 | * h=
|
---|
| 551 | * H=
|
---|
| 552 | */
|
---|
| 553 | virtual void EvaluateParameters(P *s )=0;
|
---|
| 554 |
|
---|
| 555 | /*!
|
---|
| 556 | * \brief virtual method where the a posteriori state vector must be computed
|
---|
| 557 | * \param in : the a priori set of particles at time k
|
---|
| 558 | * \param out : the a posteriori set of particles at time k
|
---|
| 559 | */
|
---|
| 560 | virtual void Update(S<P> *in,S<P> *out);
|
---|
| 561 |
|
---|
| 562 | /*!
|
---|
| 563 | * \brief Destructor
|
---|
| 564 | */
|
---|
| 565 | virtual ~NonLinearMeasureEquation(){}
|
---|
| 566 |
|
---|
| 567 | protected :
|
---|
| 568 |
|
---|
| 569 | Vector _h; /*!< vector where h(X) is stored */
|
---|
| 570 |
|
---|
| 571 | };
|
---|
| 572 |
|
---|
| 573 | // Particle non linear measure equation member functions
|
---|
| 574 | template <class D,template <class> class S, class P> void NonLinearMeasureEquation<D,S,P>::Allocator(const size_t &state_size,const size_t &data_size){
|
---|
| 575 | _h=ZeroVector(data_size);
|
---|
| 576 | }
|
---|
| 577 |
|
---|
| 578 |
|
---|
| 579 | template <class D,template <class> class S, class P> void NonLinearMeasureEquation<D,S,P>::Update(S<P> *in,S<P> *out){
|
---|
| 580 |
|
---|
| 581 | for(size_t i=0;i<in->particles.size();i++){
|
---|
| 582 | EvaluateParameters(&in->particles[i]);
|
---|
| 583 |
|
---|
| 584 | out->particles[i].weight=pdf(MeasureEquation<D,S,P>::_ZDistribution,_h)*in->particles[i].weight;
|
---|
| 585 |
|
---|
| 586 | }
|
---|
| 587 |
|
---|
| 588 | out->NormalizeWeights();
|
---|
| 589 | }
|
---|
| 590 |
|
---|
| 591 | } // namespace particle
|
---|
| 592 | } // namespace filter
|
---|
| 593 |
|
---|
| 594 |
|
---|
| 595 | #endif // __PARTICLE_FILTERING_BASE__
|
---|