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